From efb22b7ca4405f6d6378b16edc8dd006e8f9ab7f Mon Sep 17 00:00:00 2001 From: Rob Falck Date: Wed, 11 Dec 2024 14:16:39 -0500 Subject: [PATCH] Trajectory results report plots now plot vector variables, respect the x_name argument, and avoid noisy plots when values are nearly constant. (#1135) * Bokeh timeseries reports now respect the x_name argument that may be passed to run_problem via plot_kwargs * Bokeh traj results report tests now run with reports enabled. Changed handling of indexing with *idxs to be compliant with older versions of python. * more cleanup of test_pycodestyle --- .../test/ex_brachistochrone_vector_states.py | 11 +- .../test_ex_brachistochrone_vector_states.py | 54 ++++++- ...est_ex_two_burn_orbit_raise_bokeh_plots.py | 1 - .../test/test_ex_min_time_climb.py | 100 ++++++++++--- dymos/run_problem.py | 7 +- dymos/test/test_pycodestyle.py | 33 ++++- .../timeseries/bokeh_timeseries_report.py | 137 ++++++++++++------ 7 files changed, 266 insertions(+), 77 deletions(-) diff --git a/dymos/examples/brachistochrone/test/ex_brachistochrone_vector_states.py b/dymos/examples/brachistochrone/test/ex_brachistochrone_vector_states.py index 8840492ba..acbdcf691 100644 --- a/dymos/examples/brachistochrone/test/ex_brachistochrone_vector_states.py +++ b/dymos/examples/brachistochrone/test/ex_brachistochrone_vector_states.py @@ -7,7 +7,8 @@ def brachistochrone_min_time(transcription='gauss-lobatto', num_segments=8, transcription_order=3, grid_type='lgl', compressed=True, optimizer='SLSQP', dynamic_simul_derivs=True, force_alloc_complex=False, - solve_segments=False, run_driver=True): + solve_segments=False, run_driver=True, simulate=False, + make_plots=False): p = om.Problem(model=om.Group()) if optimizer == 'SNOPT': @@ -78,9 +79,11 @@ def brachistochrone_min_time(transcription='gauss-lobatto', num_segments=8, tran phase.set_state_val('v', [0, 9.9]) phase.set_control_val('theta', [5, 100]) phase.set_parameter_val('g', 9.80665) - p.run_model() - if run_driver: - p.run_driver() + + dm.run_problem(p, + run_driver=run_driver, + simulate=simulate, + make_plots=make_plots) return p diff --git a/dymos/examples/brachistochrone/test/test_ex_brachistochrone_vector_states.py b/dymos/examples/brachistochrone/test/test_ex_brachistochrone_vector_states.py index d50fbc136..f2ef1f443 100644 --- a/dymos/examples/brachistochrone/test/test_ex_brachistochrone_vector_states.py +++ b/dymos/examples/brachistochrone/test/test_ex_brachistochrone_vector_states.py @@ -1,11 +1,20 @@ +import importlib +import os +import pathlib import unittest from numpy.testing import assert_almost_equal +import sys + from openmdao.utils.general_utils import set_pyoptsparse_opt, printoptions -from openmdao.utils.testing_utils import use_tempdirs +from openmdao.utils.testing_utils import use_tempdirs, set_env_vars_context +from openmdao.utils.tests.test_hooks import hooks_active +import dymos as dm import dymos.examples.brachistochrone.test.ex_brachistochrone_vector_states as ex_brachistochrone_vs -from dymos.utils.testing_utils import assert_check_partials +from dymos.utils.testing_utils import assert_check_partials, _get_reports_dir + +bokeh_available = importlib.util.find_spec('bokeh') is not None OPT, OPTIMIZER = set_pyoptsparse_opt('SNOPT') @@ -13,6 +22,14 @@ @use_tempdirs class TestBrachistochroneVectorStatesExample(unittest.TestCase): + def setUp(self): + self.testflo_running = os.environ.pop('TESTFLO_RUNNING', None) + + def tearDown(self): + # restore what was there before running the test + if self.testflo_running is not None: + os.environ['TESTFLO_RUNNING'] = self.testflo_running + def assert_results(self, p): t_initial = p.get_val('traj0.phase0.timeseries.time')[0] t_final = p.get_val('traj0.phase0.timeseries.time')[-1] @@ -45,7 +62,7 @@ def assert_results(self, p): def assert_partials(self, p): with printoptions(linewidth=1024, edgeitems=100): - cpd = p.check_partials(method='cs', compact_print=True) + cpd = p.check_partials(method='cs', compact_print=True, out_stream=None) assert_check_partials(cpd) def test_ex_brachistochrone_vs_radau_compressed(self): @@ -97,6 +114,37 @@ def test_ex_brachistochrone_vs_birkhoff(self): self.assert_results(p) self.assert_partials(p) + @unittest.skipIf(not bokeh_available, 'bokeh unavailable') + @hooks_active + def test_bokeh_plots(self): + with set_env_vars_context(OPENMDAO_REPORTS='1'): + with dm.options.temporary(plots='bokeh'): + p = ex_brachistochrone_vs.brachistochrone_min_time(transcription='radau-ps', + compressed=False, + force_alloc_complex=True, + run_driver=True, + simulate=True, + make_plots=True) + + self.assert_results(p) + self.assert_partials(p) + + html_file = pathlib.Path(_get_reports_dir(p)) / 'traj0_results_report.html' + self.assertTrue(html_file.exists(), msg=f'{html_file} does not exist!') + + with open(html_file) as f: + html_data = f.read() + + expected_labels = ['"axis_label":"pos[0] (m)"', + '"axis_label":"pos[1] (m)"', + '"axis_label":"v (m/s)"', + '"axis_label":"theta (deg)"'] + + for label in expected_labels: + self.assertIn(label, html_data) + + self.assertNotIn('"axis_label":"pos (m)"', html_data) + if __name__ == "__main__": unittest.main() diff --git a/dymos/examples/finite_burn_orbit_raise/test/test_ex_two_burn_orbit_raise_bokeh_plots.py b/dymos/examples/finite_burn_orbit_raise/test/test_ex_two_burn_orbit_raise_bokeh_plots.py index 9519e7ca9..86ec7533c 100644 --- a/dymos/examples/finite_burn_orbit_raise/test/test_ex_two_burn_orbit_raise_bokeh_plots.py +++ b/dymos/examples/finite_burn_orbit_raise/test/test_ex_two_burn_orbit_raise_bokeh_plots.py @@ -13,7 +13,6 @@ from openmdao.utils.general_utils import set_pyoptsparse_opt from openmdao.utils.testing_utils import use_tempdirs, set_env_vars_context - import dymos as dm from dymos.examples.finite_burn_orbit_raise.finite_burn_orbit_raise_problem import two_burn_orbit_raise_problem from dymos.utils.testing_utils import _get_reports_dir diff --git a/dymos/examples/min_time_climb/test/test_ex_min_time_climb.py b/dymos/examples/min_time_climb/test/test_ex_min_time_climb.py index 32dfaa1b3..4b6c8a257 100644 --- a/dymos/examples/min_time_climb/test/test_ex_min_time_climb.py +++ b/dymos/examples/min_time_climb/test/test_ex_min_time_climb.py @@ -1,3 +1,5 @@ +import importlib +import os import unittest import numpy as np from numpy.polynomial import Polynomial as P @@ -9,18 +11,21 @@ import openmdao.api as om from openmdao.utils.assert_utils import assert_near_equal -from dymos.utils.testing_utils import assert_timeseries_near_equal +from dymos.utils.testing_utils import assert_timeseries_near_equal, _get_reports_dir from dymos.utils.introspection import get_promoted_vars import dymos as dm from dymos.examples.min_time_climb.min_time_climb_ode import MinTimeClimbODE from dymos.utils.misc import om_version -from openmdao.utils.testing_utils import use_tempdirs, require_pyoptsparse +from openmdao.utils.testing_utils import use_tempdirs, require_pyoptsparse, set_env_vars_context + + +bokeh_available = importlib.util.find_spec('bokeh') is not None def min_time_climb(optimizer='SLSQP', num_seg=3, transcription='gauss-lobatto', transcription_order=3, force_alloc_complex=False, add_rate=False, time_name='time', - simulate=True, path_constraints=True): + simulate=True, path_constraints=True, make_plots=False): p = om.Problem(model=om.Group()) @@ -139,7 +144,7 @@ def min_time_climb(optimizer='SLSQP', num_seg=3, transcription='gauss-lobatto', with np.printoptions(linewidth=1024, edgeitems=1024): p.check_partials(compact_print=False, method='cs', show_only_incorrect=True, out_stream=f) - dm.run_problem(p, simulate=simulate, make_plots=False) + dm.run_problem(p, simulate=simulate, make_plots=make_plots, plot_kwargs={'x_name': time_name}) return p @@ -293,35 +298,90 @@ def test_results_birkhoff(self): self._test_timeseries_units(p) + +@use_tempdirs +class TestMinTimeClimbWithReports(TestMinTimeClimb): + + def setUp(self): + self.testflo_running = os.environ.pop('TESTFLO_RUNNING', None) + + def tearDown(self): + # restore what was there before running the test + if self.testflo_running is not None: + os.environ['TESTFLO_RUNNING'] = self.testflo_running + + def _test_traj_results_report(self, p): + html_file = _get_reports_dir(p) / 'traj_results_report.html' + self.assertTrue(html_file.exists(), msg=f'{html_file} does not exist!') + + with open(html_file) as f: + html_data = f.read() + + expected_labels = ['"axis_label":"alpha (deg)"', + '"axis_label":"t (s)"', + '"axis_label":"CD (None)"', + '"axis_label":"CD0 (None)"', + '"axis_label":"CL (None)"', + '"axis_label":"CLa (None)"', + '"axis_label":"f_drag (N)"', + '"axis_label":"f_lift (lbf)"', + '"axis_label":"gam (rad)"', + '"axis_label":"h (m)"', + '"axis_label":"kappa (None)"', + '"axis_label":"m (kg)"', + '"axis_label":"m_dot (kg/s)"', + '"axis_label":"mach (None)"', + '"axis_label":"mach_rate (None)"', + '"axis_label":"q (N/m**2)"', + '"axis_label":"r (m)"', + '"axis_label":"thrust (lbf)"', + '"axis_label":"v (m/s)"'] + + for label in expected_labels: + self.assertIn(label, html_data) + @require_pyoptsparse(optimizer='IPOPT') + @unittest.skipIf(not bokeh_available, 'bokeh is not available') def test_results_gauss_lobatto_renamed_time(self): - NUM_SEG = 12 - ORDER = 3 - p = min_time_climb(optimizer='IPOPT', num_seg=NUM_SEG, transcription_order=ORDER, - transcription='gauss-lobatto', add_rate=True, time_name='t') + with set_env_vars_context(OPENMDAO_REPORTS='1'): + with dm.options.temporary(plots='bokeh'): + NUM_SEG = 12 + ORDER = 3 + p = min_time_climb(optimizer='IPOPT', num_seg=NUM_SEG, transcription_order=ORDER, + force_alloc_complex=True, + transcription='gauss-lobatto', add_rate=True, time_name='t', + make_plots=True) - self._test_results(p, time_name='t') + self._test_results(p, time_name='t') - self._test_wilcard_outputs(p) + self._test_wilcard_outputs(p) - self._test_timeseries_units(p) + self._test_timeseries_units(p) + + self._test_mach_rate(p, time_name='t') - self._test_mach_rate(p, time_name='t') + self._test_traj_results_report(p) @require_pyoptsparse(optimizer='IPOPT') + @unittest.skipIf(not bokeh_available, 'bokeh is not available') def test_results_radau_renamed_time(self): - NUM_SEG = 15 - ORDER = 3 - p = min_time_climb(optimizer='IPOPT', num_seg=NUM_SEG, transcription_order=ORDER, - transcription='radau-ps', add_rate=True, time_name='t', force_alloc_complex=True) + with set_env_vars_context(OPENMDAO_REPORTS='1'): + with dm.options.temporary(plots='bokeh'): + NUM_SEG = 15 + ORDER = 3 + p = min_time_climb(optimizer='IPOPT', num_seg=NUM_SEG, transcription_order=ORDER, + transcription='radau-ps', add_rate=True, time_name='t', + force_alloc_complex=True, make_plots=True) - self._test_results(p, time_name='t') + self._test_results(p, time_name='t') - self._test_wilcard_outputs(p) + self._test_wilcard_outputs(p) - self._test_timeseries_units(p) + self._test_timeseries_units(p) + + self._test_mach_rate(p, plot=False, time_name='t') - self._test_mach_rate(p, plot=False, time_name='t') + self._test_traj_results_report(p) if __name__ == '__main__': # pragma: no cover diff --git a/dymos/run_problem.py b/dymos/run_problem.py index 8a3fbfd26..63af00dbb 100755 --- a/dymos/run_problem.py +++ b/dymos/run_problem.py @@ -122,7 +122,6 @@ def run_problem(problem, refine_method='hp', refine_iteration_limit=0, run_drive sims[subsys.pathname] = sim_prob if make_plots: - if om_version()[0] > (3, 34, 2): outputs_dir = problem.get_outputs_dir() if os.sep in str(solution_record_file): @@ -142,13 +141,15 @@ def run_problem(problem, refine_method='hp', refine_iteration_limit=0, run_drive _sol_record_file = solution_record_file _sim_record_file = None if not simulate else simulation_record_file + _plot_kwargs = plot_kwargs if plot_kwargs is not None else {} + if dymos_options['plots'] == 'bokeh': from dymos.visualization.timeseries.bokeh_timeseries_report import make_timeseries_report make_timeseries_report(prob=problem, solution_record_file=_sol_record_file, - simulation_record_file=_sim_record_file) + simulation_record_file=_sim_record_file, + **_plot_kwargs) else: - _plot_kwargs = plot_kwargs if plot_kwargs is not None else {} plots_dir = problem.get_reports_dir() / 'plots' timeseries_plots(_sol_record_file, simulation_record_file=_sim_record_file, diff --git a/dymos/test/test_pycodestyle.py b/dymos/test/test_pycodestyle.py index 2325a31b4..05d4f9caf 100644 --- a/dymos/test/test_pycodestyle.py +++ b/dymos/test/test_pycodestyle.py @@ -1,4 +1,6 @@ import os +import pathlib +import subprocess import sys import unittest @@ -33,6 +35,33 @@ def _discover_python_files(path): return python_files +def _get_tracked_python_files(git_root): + """ + Given a git root directory + + Parameters + ---------- + git_root : str or Path + The root directory of the git repository. + + Returns + ------- + set + Python files in the given git directory that a tracked by git. + """ + _git_root = str(git_root) + file_list = _discover_python_files(_git_root) + untracked_set = set( + subprocess.run( + ['git', 'ls-files', _git_root, '--exclude-standard', '--others'], + stdout=subprocess.PIPE, + text=True + ).stdout.splitlines() + ) + untracked_set = {str(pathlib.Path(f).absolute()) for f in untracked_set if f.endswith('.py')} + return set(file_list) - untracked_set + + @unittest.skipIf(pycodestyle is None, "This test requires pycodestyle") class TestPyCodeStyle(unittest.TestCase): @@ -47,6 +76,8 @@ def test_pycodestyle(self): dymos_path = os.path.split(dymos.__file__)[0] pyfiles = _discover_python_files(dymos_path) + files_to_check = _get_tracked_python_files(dymos_path) + style = pycodestyle.StyleGuide(ignore=['E226', # missing whitespace around arithmetic operator 'E241', # multiple spaces after ',' 'W504', # line break after binary operator @@ -65,7 +96,7 @@ def test_pycodestyle(self): try: sys.stdout = buff_out = StringIO() sys.stdout = buff_err = StringIO() - report = style.check_files(pyfiles) + report = style.check_files(files_to_check) finally: sys.stdout = save_out sys.stderr = save_err diff --git a/dymos/visualization/timeseries/bokeh_timeseries_report.py b/dymos/visualization/timeseries/bokeh_timeseries_report.py index ca7ef5bb4..9c7999412 100644 --- a/dymos/visualization/timeseries/bokeh_timeseries_report.py +++ b/dymos/visualization/timeseries/bokeh_timeseries_report.py @@ -1,16 +1,20 @@ from collections import ChainMap import datetime +import itertools from pathlib import Path import os.path +import numpy as np + from dymos.trajectory.trajectory import Trajectory from dymos.phase.phase import Phase try: from bokeh.io import save from bokeh.layouts import column, grid, row - from bokeh.models import Legend, DataTable, Div, ColumnDataSource, TableColumn, \ - TabPanel, Tabs, CheckboxButtonGroup, CustomJS, MultiChoice + from bokeh.models import Legend, DataRange1d, DataTable, Div, \ + ColumnDataSource, TableColumn, TabPanel, Tabs, \ + CheckboxButtonGroup, CustomJS, MultiChoice from bokeh.plotting import figure, curdoc import bokeh.palettes as bp import bokeh.resources as bokeh_resources @@ -306,8 +310,6 @@ def _load_data_sources(traj_and_phase_meta=None, solution_record_file=None, simu phase_sol_data = data_dict[traj_path]['sol_data_by_phase'][phase_name] = {} phase_sim_data = data_dict[traj_path]['sim_data_by_phase'][phase_name] = {} - # param_outputs = {op: meta for op, meta in outputs.items() - # if op.startswith(f'{phase_path}.param_comp.parameter_vals:')} ts_outputs = {op: meta for op, meta in outputs.items() if op.startswith(f'{phase_path}.timeseries.')} @@ -342,9 +344,11 @@ def _load_data_sources(traj_and_phase_meta=None, solution_record_file=None, simu var_name = prom_name.split('.')[-1] if sol_case: - phase_sol_data[var_name] = sol_case.get_val(prom_name, units=ts_units_dict[var_name]) + data = sol_case.get_val(prom_name, units=ts_units_dict[var_name]) + phase_sol_data[var_name] = data if sim_case: - phase_sim_data[var_name] = sim_case.get_val(prom_name, units=ts_units_dict[var_name]) + data = sim_case.get_val(prom_name, units=ts_units_dict[var_name]) + phase_sim_data[var_name] = data return data_dict @@ -383,6 +387,26 @@ def _gather_system_options(model, options, sys_cls=None, rank=0,): return system_options +def _new_figure(x_name, y_name, x_units, y_units, margin, x_range=None): + fig_kwargs = {'x_range': x_range} if x_range is not None else {} + + tool_tips = [(f'{x_name}', f'@{x_name}'), (f'{y_name}', f'@{y_name}')] + + fig = figure(tools='pan,box_zoom,xwheel_zoom,hover,undo,reset,save', + tooltips=tool_tips, + x_axis_label=f'{x_name} ({x_units})', + y_axis_label=f'{y_name} ({y_units})', + toolbar_location='above', + sizing_mode='stretch_both', + min_height=250, max_height=300, + margin=margin, + **fig_kwargs) + fig.xaxis.axis_label_text_font_size = '10pt' + fig.yaxis.axis_label_text_font_size = '10pt' + fig.toolbar.autohide = True + return fig + + def make_timeseries_report(prob, solution_record_file=None, simulation_record_file=None, x_name='time', ncols=2, margin=10, theme='light_minimal'): """ @@ -454,55 +478,79 @@ def make_timeseries_report(prob, solution_record_file=None, simulation_record_fi # Plot the timeseries ts_units_dict = source_data[traj_path]['timeseries_units'] - figures = [] + figures = {} + legend_data_per_figure = {} x_range = None - for var_name in sorted(ts_units_dict.keys(), key=str.casefold): - fig_kwargs = {'x_range': x_range} if x_range is not None else {} + # var_name is the actual dymos variable name, without any index information. + # var_name_with_idxs is the variable name with index information. - tool_tips = [(f'{x_name}', f'@{x_name}'), (f'{var_name}', f'@{var_name}')] - - fig = figure(tools='pan,box_zoom,xwheel_zoom,hover,undo,reset,save', - tooltips=tool_tips, - x_axis_label=f'{x_name} ({ts_units_dict[x_name]})', - y_axis_label=f'{var_name} ({ts_units_dict[var_name]})', - toolbar_location='above', - sizing_mode='stretch_both', - min_height=250, max_height=300, - margin=margin, - **fig_kwargs) - fig.xaxis.axis_label_text_font_size = '10pt' - fig.yaxis.axis_label_text_font_size = '10pt' - fig.toolbar.autohide = True - legend_data = [] - if x_range is None: - x_range = fig.x_range + for var_name in sorted(ts_units_dict.keys(), key=str.casefold): for i, phase_name in enumerate(phase_names): color = colors[i % 20] sol_data = source_data[traj_path]['sol_data_by_phase'][phase_name] sim_data = source_data[traj_path]['sim_data_by_phase'][phase_name] - sol_source = ColumnDataSource(sol_data) - sim_source = ColumnDataSource(sim_data) + if x_name in sol_data and var_name in sol_data: - legend_items = [] - if sol_data: - sol_plot = fig.scatter(x='time', y=var_name, source=sol_source, - color=color, size=5) - sol_plot.tags.extend(['sol', f'phase:{phase_name}']) - legend_items.append(sol_plot) - if sim_data: - sim_plot = fig.line(x='time', y=var_name, source=sim_source, color=color) - sim_plot.tags.extend(['sim', f'phase:{phase_name}']) - legend_items.append(sim_plot) - legend_data.append((phase_name, legend_items)) - - legend = Legend(items=legend_data, location='center', label_text_font_size='8pt') + shape = sol_data[var_name].shape[1:] + indices = list(itertools.product(*(range(dim) for dim in shape))) + if np.prod(shape) > 1: + sources = {} + for idxs in indices: + str_idxs = ','.join([str(i) for i in idxs]) + # Bokeh ColumnDataSource doesn't allow special characters in keys, + # but we want the y_axis label to show the indices of the columns + # being plotted as 'varname[i,j,k]'. + sources[f'{var_name}[{str_idxs}]'] = s = f'{var_name}_{str_idxs.replace(",", "_")}' + sol_data_column = sol_data[var_name][(slice(None), *idxs)] + sol_data[s] = sol_data_column + sim_data_column = sim_data[var_name][(slice(None), *idxs)] + sim_data[s] = sim_data_column + else: + sources = {var_name: var_name} + sol_source = ColumnDataSource(sol_data) + sim_source = ColumnDataSource(sim_data) + + for var_name_with_idxs, _source in sources.items(): + legend_items = [] + + if var_name_with_idxs in figures: + fig = figures[var_name_with_idxs] + fig_legend_data = legend_data_per_figure[var_name_with_idxs] + else: + fig = _new_figure(x_name=x_name, y_name=var_name_with_idxs, + x_units=ts_units_dict[x_name], + y_units=ts_units_dict[var_name], + margin=margin, + x_range=x_range) + figures[var_name_with_idxs] = fig + fig_legend_data = legend_data_per_figure[var_name_with_idxs] = [] + if sol_data: + sol_plot = fig.scatter(x=x_name, y=_source, source=sol_source, + color=color, size=5) + sol_plot.tags.extend(['sol', f'phase:{phase_name}']) + legend_items.append(sol_plot) + if sim_data: + sim_plot = fig.line(x=x_name, y=_source, source=sim_source, color=color) + sim_plot.tags.extend(['sim', f'phase:{phase_name}']) + legend_items.append(sim_plot) + fig_legend_data.append((phase_name, legend_items)) + if x_range is None: + x_range = fig.x_range + fig.y_range = DataRange1d(min_interval=1.0E-12) + + figures[var_name_with_idxs] = fig + + # Add the legend data to each figure + for var_name_with_idxs, fig in figures.items(): + legend_data = legend_data_per_figure[var_name_with_idxs] + legend = Legend(items=legend_data, location='center', label_text_font_size='8pt', + orientation='vertical') fig.add_layout(legend, 'right') - figures.append(fig) # Since we're putting figures in two columns, make sure we have an even number of things to put in the layout. if len(figures) % 2 == 1: - figures.append(None) + figures['__NONE__'] = None param_panels = [TabPanel(child=table, title=f'{phase_names[i]} parameters') for i, table in enumerate(param_tables)] @@ -522,7 +570,7 @@ def make_timeseries_report(prob, solution_record_file=None, simulation_record_fi phase_select_row = row(children=[Div(text='Plot phases:'), phase_select], sizing_mode='stretch_width') - figures_grid = grid(children=figures, ncols=ncols, sizing_mode='stretch_both') + figures_grid = grid(children=list(figures.values()), ncols=ncols, sizing_mode='stretch_both') ts_layout = column(children=[sol_sim_row, phase_select_row, figures_grid], sizing_mode='stretch_both') @@ -550,6 +598,5 @@ def make_timeseries_report(prob, solution_record_file=None, simulation_record_fi phase_select=phase_select))) # Save - save(report_layout, filename=report_path, title=f'trajectory results for {traj_path}', resources=bokeh_resources.INLINE)