Repository: OpenMDAO/dymos

Issue:
Timeseries report should display each element of vectorized states/controls/parameters separately.
### Description

When a state is a vector/matrix/tensor quantity at each given point in time the resulting timeseries plot is nonsensical.
Dymos should extract each element from the state (or control/parameter/etc) in this case and plot it along with a legend to indicate which plot belongs to which index.

### Example

See the matrix-states cannonball test case.

Ideally the z plot would have 4 series in it, labeled in the legend as `z[0, 0]`, `z[0, 1]`, `z[1, 0]`, and `z[1, 1]`.

### Dymos Version

1.9.2-dev

### Relevant environment information

_No response_

Discussion:
coveralls: 
[![Coverage Status](https://coveralls.io/builds/71216281/badge)](https://coveralls.io/builds/71216281)

coverage: 93.605% (-0.2%) from 93.802%
when pulling **3f6e85b6a84537e5d30924c7c6434d63d8f75829 on robfalck:timeseries_fixes**
into **d324563198e0c455ee2d38b69a282a16d3cacbe7 on OpenMDAO:master**.


Patch:
diff --git a/dymos/run_problem.py b/dymos/run_problem.py
index 8a3fbfd2..63af00db 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/visualization/timeseries/bokeh_timeseries_report.py b/dymos/visualization/timeseries/bokeh_timeseries_report.py
index ca7ef5bb..9c799941 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)


Test patch:
diff --git a/dymos/examples/brachistochrone/test/ex_brachistochrone_vector_states.py b/dymos/examples/brachistochrone/test/ex_brachistochrone_vector_states.py
index 8840492b..acbdcf69 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 @@ from dymos.examples.brachistochrone.brachistochrone_vector_states_ode \
 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 d50fbc13..f2ef1f44 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 @@ OPT, OPTIMIZER = set_pyoptsparse_opt('SNOPT')
 @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 @@ class TestBrachistochroneVectorStatesExample(unittest.TestCase):
 
     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 @@ class TestBrachistochroneVectorStatesExample(unittest.TestCase):
         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 9519e7ca..86ec7533 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.testing_utils import require_pyoptsparse
 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 32dfaa1b..4b6c8a25 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 @@ except ImportError:
 
 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 @@ class TestMinTimeClimb(unittest.TestCase):
 
         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/test/test_pycodestyle.py b/dymos/test/test_pycodestyle.py
index 2325a31b..05d4f9ca 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 @@ class TestPyCodeStyle(unittest.TestCase):
         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 @@ class TestPyCodeStyle(unittest.TestCase):
         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