Merge branch 'master' of https://github.com/OpenMDAO/dymos into radau…

…_refactor

.github/workflows/dymos_docs_workflow.yml

@@ -84,13 +84,11 @@ jobs:
         with:
           python-version: ${{ matrix.PY }}
           channels: conda-forge
+          conda-remove-defaults: true
 
       - name: Install Numpy/Scipy
         shell: bash -l {0}
         run: |
-          echo "Make sure we are not using anaconda packages"
-          conda config --remove channels defaults
-
           echo "============================================================="
           echo "Install Numpy/Scipy"
           echo "============================================================="

.github/workflows/dymos_tests_workflow.yml

@@ -222,14 +222,12 @@ jobs:
         with:
           python-version: ${{ matrix.PY }}
           channels: conda-forge
+          conda-remove-defaults: true
 
       - name: Install Numpy/Scipy
         if: ${{ ! matrix.EXCLUDE }}
         shell: bash -l {0}
         run: |
-          echo "Make sure we are not using anaconda packages"
-          conda config --remove channels defaults
-
           echo "============================================================="
           echo "Install Numpy/Scipy"
           echo "============================================================="

docs/dymos_book/examples/multi_phase_cannonball/multi_phase_cannonball.ipynb

@@ -167,10 +167,10 @@
    "outputs": [],
    "source": [
     "import numpy as np\n",
-    "from scipy.interpolate import interp1d\n",
     "import matplotlib.pyplot as plt\n",
     "\n",
     "import openmdao.api as om\n",
+    "from openmdao.components.interp_util.interp import InterpND\n",
     "\n",
     "import dymos as dm\n",
     "from dymos.models.atmosphere.atmos_1976 import USatm1976Data\n",
@@ -271,9 +271,9 @@
     "\n",
     "        alt_data = USatm1976Data.alt * om.unit_conversion('ft', 'm')[0]\n",
     "        rho_data = USatm1976Data.rho * om.unit_conversion('slug/ft**3', 'kg/m**3')[0]\n",
-    "        self.rho_interp = interp1d(np.array(alt_data, dtype=complex),\n",
-    "                                   np.array(rho_data, dtype=complex),\n",
-    "                                   kind='linear')\n",
+    "        self.rho_interp = InterpND(points=np.array(alt_data),\n",
+    "                                   values=np.array(rho_data),\n",
+    "                                   method='slinear').interpolate\n",
     "\n",
     "    def compute(self, inputs, outputs):\n",
     "\n",
@@ -560,7 +560,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.11.4"
+   "version": "3.12.6"
   }
  },
  "nbformat": 4,

dymos/examples/brachistochrone/test/test_brachistochrone_integrated_control.py

@@ -2,9 +2,9 @@
 import unittest
 
 import numpy as np
-from scipy.interpolate import interp1d
 
 import openmdao.api as om
+from openmdao.components.interp_util.interp import InterpND
 from openmdao.utils.testing_utils import use_tempdirs
 
 import dymos as dm
@@ -150,17 +150,20 @@ def test_brachistochrone_integrated_control_gauss_lobatto(self):
         theta_dot_sim = sim_out.get_val('phase0.timeseries.theta_dot')
         time_sim = sim_out.get_val('phase0.timeseries.time')
 
-        x_interp = interp1d(time_sim[:, 0], x_sim[:, 0])
-        y_interp = interp1d(time_sim[:, 0], y_sim[:, 0])
-        v_interp = interp1d(time_sim[:, 0], v_sim[:, 0])
-        theta_interp = interp1d(time_sim[:, 0], theta_sim[:, 0])
-        theta_dot_interp = interp1d(time_sim[:, 0], theta_dot_sim[:, 0])
+        # need unique (monotonically increasing) times for interpolation
+        times, idxs = np.unique(time_sim[:, 0], return_index=True)
 
-        assert_near_equal(x_interp(time_sol), x_sol, tolerance=1.0E-4)
-        assert_near_equal(y_interp(time_sol), y_sol, tolerance=1.0E-4)
-        assert_near_equal(v_interp(time_sol), v_sol, tolerance=1.0E-4)
-        assert_near_equal(theta_interp(time_sol), theta_sol, tolerance=1.0E-4)
-        assert_near_equal(theta_dot_interp(time_sol), theta_dot_sol, tolerance=1.0E-4)
+        x_interp = InterpND('slinear', times, x_sim[:, 0][idxs]).interpolate
+        y_interp = InterpND('slinear', times, y_sim[:, 0][idxs]).interpolate
+        v_interp = InterpND('slinear', times, v_sim[:, 0][idxs]).interpolate
+        theta_interp = InterpND('slinear', times, theta_sim[:, 0][idxs]).interpolate
+        theta_dot_interp = InterpND('slinear', times, theta_dot_sim[:, 0][idxs]).interpolate
+
+        assert_near_equal(x_interp(time_sol), x_sol[:, 0], tolerance=1.0E-4)
+        assert_near_equal(y_interp(time_sol), y_sol[:, 0], tolerance=1.0E-4)
+        assert_near_equal(v_interp(time_sol), v_sol[:, 0], tolerance=1.0E-4)
+        assert_near_equal(theta_interp(time_sol), theta_sol[:, 0], tolerance=1.0E-4)
+        assert_near_equal(theta_dot_interp(time_sol), theta_dot_sol[:, 0], tolerance=1.0E-4)
 
     def test_brachistochrone_integrated_control_radau_ps(self):
         import openmdao.api as om
@@ -233,17 +236,20 @@ def test_brachistochrone_integrated_control_radau_ps(self):
         theta_dot_sim = sim_case.get_val('traj.phase0.timeseries.theta_dot')
         time_sim = sim_case.get_val('traj.phase0.timeseries.time')
 
-        x_interp = interp1d(time_sim[:, 0], x_sim[:, 0])
-        y_interp = interp1d(time_sim[:, 0], y_sim[:, 0])
-        v_interp = interp1d(time_sim[:, 0], v_sim[:, 0])
-        theta_interp = interp1d(time_sim[:, 0], theta_sim[:, 0])
-        theta_dot_interp = interp1d(time_sim[:, 0], theta_dot_sim[:, 0])
-
-        assert_near_equal(x_interp(time_sol), x_sol, tolerance=1.0E-4)
-        assert_near_equal(y_interp(time_sol), y_sol, tolerance=1.0E-4)
-        assert_near_equal(v_interp(time_sol), v_sol, tolerance=1.0E-4)
-        assert_near_equal(theta_interp(time_sol), theta_sol, tolerance=1.0E-4)
-        assert_near_equal(theta_dot_interp(time_sol), theta_dot_sol, tolerance=1.0E-4)
+        # need unique (monotonically increasing) times for interpolation
+        times, idxs = np.unique(time_sim[:, 0], return_index=True)
+
+        x_interp = InterpND('slinear', times, x_sim[:, 0][idxs]).interpolate
+        y_interp = InterpND('slinear', times, y_sim[:, 0][idxs]).interpolate
+        v_interp = InterpND('slinear', times, v_sim[:, 0][idxs]).interpolate
+        theta_interp = InterpND('slinear', times, theta_sim[:, 0][idxs]).interpolate
+        theta_dot_interp = InterpND('slinear', times, theta_dot_sim[:, 0][idxs]).interpolate
+
+        assert_near_equal(x_interp(time_sol), x_sol[:, 0], tolerance=1.0E-4)
+        assert_near_equal(y_interp(time_sol), y_sol[:, 0], tolerance=1.0E-4)
+        assert_near_equal(v_interp(time_sol), v_sol[:, 0], tolerance=1.0E-4)
+        assert_near_equal(theta_interp(time_sol), theta_sol[:, 0], tolerance=1.0E-4)
+        assert_near_equal(theta_dot_interp(time_sol), theta_dot_sol[:, 0], tolerance=1.0E-4)
 
 
 if __name__ == '__main__':  # pragma: no cover

dymos/examples/shuttle_reentry/test/test_reentry.py

@@ -283,16 +283,20 @@ def test_reentry_mixed_controls(self):
         sol = om.CaseReader(sol_db).get_case('final')
         sim = om.CaseReader(sim_db).get_case('final')
 
-        from scipy.interpolate import interp1d
+        from openmdao.components.interp_util.interp import InterpND
 
         t_sol = sol.get_val('traj.phase0.timeseries.time')
         beta_sol = sol.get_val('traj.phase0.timeseries.beta', units='deg')
 
         t_sim = sim.get_val('traj.phase0.timeseries.time')
         beta_sim = sim.get_val('traj.phase0.timeseries.beta', units='deg')
 
-        sol_interp = interp1d(t_sol.ravel(), beta_sol.ravel())
-        sim_interp = interp1d(t_sim.ravel(), beta_sim.ravel())
+        # need unique (monotonically increasing) times for interpolation
+        t_sol_u, t_sol_idx = np.unique(t_sol, return_index=True)
+        t_sim_u, t_sim_idx = np.unique(t_sim, return_index=True)
+
+        sol_interp = InterpND('slinear', t_sol_u.ravel(), beta_sol[t_sol_idx].ravel()).interpolate
+        sim_interp = InterpND('slinear', t_sim_u.ravel(), beta_sim[t_sim_idx].ravel()).interpolate
 
         t = np.linspace(0, t_sol.ravel()[-1], 1000)
 

dymos/test/test_load_case.py

@@ -123,6 +123,7 @@ def test_load_case_unchanged_grid_polynomial_control(self):
                           case.get_val('phase0.controls:theta'))
 
     def test_load_case_lgl_to_radau(self):
+        import numpy as np
         import openmdao.api as om
         from openmdao.utils.assert_utils import assert_near_equal
         import dymos as dm
@@ -155,11 +156,15 @@ def test_load_case_lgl_to_radau(self):
         time_val = outputs['phase0.timeseries.timeseries_comp.time']['val']
         theta_val = outputs['phase0.timeseries.timeseries_comp.theta']['val']
 
+        time_val_uniq, idx = np.unique(time_val, return_index=True)
+        theta_val_uniq = theta_val[idx]
+
         assert_near_equal(q['phase0.timeseries.timeseries_comp.theta'],
-                          q.model.phase0.interp(xs=time_val, ys=theta_val, nodes='all'),
+                          q.model.phase0.interp(xs=time_val_uniq, ys=theta_val_uniq, nodes='all'),
                           tolerance=1.0E-3)
 
     def test_load_case_radau_to_lgl(self):
+        import numpy as np
         import openmdao.api as om
         from openmdao.utils.assert_utils import assert_near_equal
         import dymos as dm
@@ -192,11 +197,18 @@ def test_load_case_radau_to_lgl(self):
         time_q = q.get_val('phase0.timeseries.time')
         theta_q = q.get_val('phase0.timeseries.theta')
 
-        assert_near_equal(q.model.phase0.interp(xs=time_p, ys=theta_p, nodes='all'),
-                          q.model.phase0.interp(xs=time_q, ys=theta_q, nodes='all'),
+        time_p_unique, p_idx = np.unique(time_p, return_index=True)
+        theta_p_unique = theta_p[p_idx]
+
+        time_q_unique, q_idx = np.unique(time_q, return_index=True)
+        theta_q_unique = theta_q[q_idx]
+
+        assert_near_equal(q.model.phase0.interp(xs=time_p_unique, ys=theta_p_unique, nodes='all'),
+                          q.model.phase0.interp(xs=time_q_unique, ys=theta_q_unique, nodes='all'),
                           tolerance=1.0E-2)
 
     def test_load_case_radau_to_birkhoff(self):
+        import numpy as np
         import openmdao.api as om
         from openmdao.utils.assert_utils import assert_near_equal
         import dymos as dm
@@ -258,8 +270,14 @@ def test_load_case_radau_to_birkhoff(self):
         v0_q = q.get_val('phase0.initial_states:v')
         vf_q = q.get_val('phase0.final_states:v')
 
-        assert_near_equal(q.model.phase0.interp(xs=time_p, ys=theta_p, nodes='all'),
-                          q.model.phase0.interp(xs=time_q, ys=theta_q, nodes='all'),
+        time_p_unique, p_idx = np.unique(time_p, return_index=True)
+        theta_p_unique = theta_p[p_idx]
+
+        time_q_unique, q_idx = np.unique(time_q, return_index=True)
+        theta_q_unique = theta_q[q_idx]
+
+        assert_near_equal(q.model.phase0.interp(xs=time_p_unique, ys=theta_p_unique, nodes='all'),
+                          q.model.phase0.interp(xs=time_q_unique, ys=theta_q_unique, nodes='all'),
                           tolerance=1.0E-2)
 
         assert_near_equal(x_p[0, ...], x0_q, tolerance=1.0E-5)

dymos/test/test_run_problem.py

@@ -562,7 +562,7 @@ def test_run_brachistochrone_problem_with_simulate(self):
     def test_restart_from_file(self):
         from dymos.examples.vanderpol.vanderpol_dymos import vanderpol
         from dymos.run_problem import run_problem
-        from scipy.interpolate import interp1d
+        from openmdao.components.interp_util.interp import InterpND
         from numpy.testing import assert_almost_equal
 
         # Create the Dymos problem instance
@@ -602,9 +602,9 @@ def test_restart_from_file(self):
         us = s.get_val('traj.phase0.timeseries.u')[:, 0][nodup]
 
         # create interpolation functions so that values can be looked up at matching time points
-        fx1s = interp1d(ts, x1s, kind='cubic')
-        fx0s = interp1d(ts, x0s, kind='cubic')
-        fus = interp1d(ts, us, kind='cubic')
+        fx1s = InterpND('cubic', ts, x1s).interpolate
+        fx0s = InterpND('cubic', ts, x0s).interpolate
+        fus = InterpND('cubic', ts, us).interpolate
 
         assert_almost_equal(x1q, fx1s(tq), decimal=2)
         assert_almost_equal(x0q, fx0s(tq), decimal=2)
@@ -613,7 +613,7 @@ def test_restart_from_file(self):
     def test_restart_from_case(self):
         from dymos.examples.vanderpol.vanderpol_dymos import vanderpol
         from dymos.run_problem import run_problem
-        from scipy.interpolate import interp1d
+        from openmdao.components.interp_util.interp import InterpND
         from numpy.testing import assert_almost_equal
 
         # Create the Dymos problem instance
@@ -653,9 +653,9 @@ def test_restart_from_case(self):
         us = s.get_val('traj.phase0.timeseries.u')[:, 0][nodup]
 
         # create interpolation functions so that values can be looked up at matching time points
-        fx1s = interp1d(ts, x1s, kind='cubic')
-        fx0s = interp1d(ts, x0s, kind='cubic')
-        fus = interp1d(ts, us, kind='cubic')
+        fx1s = InterpND('cubic', ts, x1s).interpolate
+        fx0s = InterpND('cubic', ts, x0s).interpolate
+        fus = InterpND('cubic', ts, us).interpolate
 
         assert_almost_equal(x1q, fx1s(tq), decimal=2)
         assert_almost_equal(x0q, fx0s(tq), decimal=2)

dymos/trajectory/test/test_trajectory.py

@@ -1396,9 +1396,9 @@ def test_invalid_linkage_args(self):
 
     def test_linkage_units(self):
         import numpy as np
-        from scipy.interpolate import interp1d
 
         import openmdao.api as om
+        from openmdao.components.interp_util.interp import InterpND
 
         import dymos as dm
         from dymos.models.atmosphere.atmos_1976 import USatm1976Data
@@ -1464,9 +1464,7 @@ def setup(self):
 
                 alt_data = USatm1976Data.alt * om.unit_conversion('ft', 'm')[0]
                 rho_data = USatm1976Data.rho * om.unit_conversion('slug/ft**3', 'kg/m**3')[0]
-                self.rho_interp = interp1d(np.array(alt_data, dtype=complex),
-                                           np.array(rho_data, dtype=complex),
-                                           kind='linear')
+                self.rho_interp = InterpND('slinear', np.array(alt_data), np.array(rho_data)).interpolate
 
             def compute(self, inputs, outputs):
 
@@ -1660,9 +1658,9 @@ def compute(self, inputs, outputs):
 
     def test_linkage_units_connected(self):
         import numpy as np
-        from scipy.interpolate import interp1d
 
         import openmdao.api as om
+        from openmdao.components.interp_util.interp import InterpND
 
         import dymos as dm
         from dymos.models.atmosphere.atmos_1976 import USatm1976Data
@@ -1728,9 +1726,7 @@ def setup(self):
 
                 alt_data = USatm1976Data.alt * om.unit_conversion('ft', 'm')[0]
                 rho_data = USatm1976Data.rho * om.unit_conversion('slug/ft**3', 'kg/m**3')[0]
-                self.rho_interp = interp1d(np.array(alt_data, dtype=complex),
-                                           np.array(rho_data, dtype=complex),
-                                           kind='linear')
+                self.rho_interp = InterpND('slinear', np.array(alt_data), np.array(rho_data)).interpolate
 
             def compute(self, inputs, outputs):
 

dymos/utils/testing_utils.py

@@ -5,7 +5,7 @@
 
 import numpy as np
 
-from scipy.interpolate import interp1d
+from scipy.interpolate import Akima1DInterpolator
 
 import openmdao.api as om
 import openmdao.utils.assert_utils as _om_assert_utils
@@ -273,7 +273,7 @@ def assert_timeseries_near_equal(t_ref, x_ref, t_check, x_check, abs_tolerance=N
     x_to_interp = x_ref_data_flattened[idxs_unique_ref, ...]
     t_check = t_check.ravel()
 
-    interp = interp1d(x=t_ref_unique, y=x_to_interp, kind='slinear', axis=0)
+    interp = Akima1DInterpolator(t_ref_unique, x_to_interp)
 
     # only want t_check in the overlapping range of t_begin and t_end
     t_check_in_range_condition = np.logical_and(t_check >= t_begin, t_check <= t_end)

joss/test/test_cannonball_for_joss.py

@@ -16,21 +16,21 @@ class TestCannonballForJOSS(unittest.TestCase):
     def test_results(self):
         # begin code for paper
         import numpy as np
-        from scipy.interpolate import interp1d
         import matplotlib
         import matplotlib.pyplot as plt
 
         import openmdao.api as om
+        from openmdao.components.interp_util.interp import InterpND
 
         import dymos as dm
         from dymos.models.atmosphere.atmos_1976 import USatm1976Data as USatm1976Data
 
         # CREATE an atmosphere interpolant
         english_to_metric_rho = om.unit_conversion('slug/ft**3', 'kg/m**3')[0]
         english_to_metric_alt = om.unit_conversion('ft', 'm')[0]
-        rho_interp = interp1d(np.array(USatm1976Data.alt * english_to_metric_alt, dtype=complex),
-                              np.array(USatm1976Data.rho * english_to_metric_rho, dtype=complex),
-                              kind='linear')
+        rho_interp = InterpND('slinear',
+                              USatm1976Data.alt * english_to_metric_alt,
+                              USatm1976Data.rho * english_to_metric_rho)
 
         class CannonballSize(om.ExplicitComponent):
             """
@@ -102,9 +102,9 @@ def compute(self, inputs, outputs):
 
                 # handle complex-step gracefully from the interpolant
                 if np.iscomplexobj(alt):
-                    rho = rho_interp(inputs['alt'])
+                    rho = rho_interp.interpolate(inputs['alt'])
                 else:
-                    rho = rho_interp(inputs['alt']).real
+                    rho = rho_interp.interpolate(inputs['alt']).real
 
                 q = 0.5 * rho * inputs['v'] ** 2
                 qS = q * S