Merge pull request PyDMD#515 from nmank/lagged_dmdc

Lagged dmdc

Author: mtezzele

pydmd/dmdc.py

@@ -191,9 +191,14 @@ class DMDc(DMDBase):
         equal than `svd_rank`. For the possible values please refer to the
         `svd_rank` parameter description above.
     :type svd_rank_omega: int or float
+    :param lag: the time lag for the snapshots. Used in fit method to generate
+        X and Y.
+    :type lag: int
     """
 
-    def __init__(self, svd_rank=0, tlsq_rank=0, opt=False, svd_rank_omega=-1):
+    def __init__(
+        self, svd_rank=0, tlsq_rank=0, opt=False, svd_rank_omega=-1, lag=1
+    ):
         # we're going to initialize Atilde when we know if B is known
         self._Atilde = None
         # remember the arguments for when we'll need them
@@ -210,6 +215,7 @@ def __init__(self, svd_rank=0, tlsq_rank=0, opt=False, svd_rank_omega=-1):
         self._snapshots_holder = None
         self._controlin = None
         self._basis = None
+        self._lag = lag
 
         self._modes_activation_bitmask_proxy = None
 
@@ -255,7 +261,7 @@ def reconstructed_data(self, control_input=None):
             else self._controlin
         )
 
-        if controlin.shape[-1] != self.dynamics.shape[-1] - 1:
+        if controlin.shape[-1] != self.dynamics.shape[-1] - self._lag:
             raise RuntimeError(
                 "The number of control inputs and the number of snapshots to "
                 "reconstruct has to be the same"
@@ -268,7 +274,8 @@ def reconstructed_data(self, control_input=None):
             [self.modes, np.diag(eigs), np.linalg.pinv(self.modes)]
         )
 
-        data = [self.snapshots[:, 0]]
+        data = [self.snapshots[:, i] for i in range(self._lag)]
+
         expected_shape = data[0].shape
 
         for i, u in enumerate(controlin.T):
@@ -298,13 +305,16 @@ def fit(self, X, I, B=None):
         :type B: numpy.ndarray or iterable
         """
         self._reset()
+        self._controlin = np.atleast_2d(np.asarray(I))
 
         self._snapshots_holder = Snapshots(X)
-        self._controlin = np.atleast_2d(np.asarray(I))
+        n_samples = self.snapshots.shape[-1]
+
+        if self._lag < 1:
+            raise ValueError("Time lag must be positive.")
 
-        n_samples = self.snapshots.shape[1]
-        X = self.snapshots[:, :-1]
-        Y = self.snapshots[:, 1:]
+        X = self.snapshots[:, : -self._lag]
+        Y = self.snapshots[:, self._lag :]
 
         self._set_initial_time_dictionary(
             {"t0": 0, "tend": n_samples - 1, "dt": 1}

tests/test_dmdc.py

@@ -26,6 +26,21 @@ def create_system_without_B():
     return {"snapshots": snapshots, "u": u, "B": B, "A": A}
 
 
+def create_system_without_B_lag(lag=1):
+    n = 5  # dimension snapshots
+    m = 15  # number snapshots
+    A = scipy.linalg.helmert(n, True)
+    B = np.random.rand(n, n) - 0.5
+    snapshots = []
+    for _ in range(lag):
+        snapshots.append(np.array([0.25] * n))
+    u = np.random.rand(n, m - lag) - 0.5
+    for i in range(m - lag):
+        snapshots.append(A.dot(snapshots[i]) + B.dot(u[:, i]))
+    snapshots = np.array(snapshots).T
+    return {"snapshots": snapshots, "u": u, "B": B, "A": A}
+
+
 def test_eigs_b_known():
     system = create_system_with_B()
     dmdc = DMDc(svd_rank=-1)
@@ -41,13 +56,29 @@ def test_eigs_b_unknown():
     assert dmdc.eigs.shape[0] == 3
 
 
+def test_eigs_b_unknown_lag():
+    lag = 3
+    system = create_system_without_B_lag(lag=lag)
+    dmdc = DMDc(svd_rank=3, opt=False, svd_rank_omega=4, lag=lag)
+    dmdc.fit(system["snapshots"], system["u"])
+    assert dmdc.eigs.shape[0] == 3
+
+
 def test_modes_b_unknown():
     system = create_system_without_B()
     dmdc = DMDc(svd_rank=3, opt=False, svd_rank_omega=4)
     dmdc.fit(system["snapshots"], system["u"])
     assert dmdc.modes.shape[1] == 3
 
 
+def test_modes_b_unknown_lag():
+    lag = 3
+    system = create_system_without_B_lag(lag=lag)
+    dmdc = DMDc(svd_rank=3, opt=False, svd_rank_omega=4, lag=lag)
+    dmdc.fit(system["snapshots"], system["u"])
+    assert dmdc.modes.shape[1] == 3
+
+
 def test_reconstruct_b_known():
     system = create_system_with_B()
     dmdc = DMDc(svd_rank=-1)
@@ -73,6 +104,16 @@ def test_reconstruct_b_unknown():
     )
 
 
+def test_reconstruct_b_unknown_lag():
+    lag = 3
+    system = create_system_without_B_lag(lag=lag)
+    dmdc = DMDc(svd_rank=-1, opt=True, lag=lag)
+    dmdc.fit(system["snapshots"], system["u"])
+    np.testing.assert_array_almost_equal(
+        dmdc.reconstructed_data(), system["snapshots"], decimal=6
+    )
+
+
 def test_atilde_b_unknown():
     system = create_system_without_B()
     dmdc = DMDc(svd_rank=-1, opt=True)
@@ -232,3 +273,10 @@ def test_correct_amplitudes():
     dmd = DMDc(svd_rank=-1)
     dmd.fit(system["snapshots"], system["u"], system["B"])
     np.testing.assert_array_almost_equal(dmd.amplitudes, dmd._b)
+
+
+def test_lag_param_b_unknown_raises():
+    system = create_system_without_B_lag(lag=3)
+    dmdc = DMDc(svd_rank=-1, opt=True, lag=0)
+    with raises(ValueError):
+        dmdc.fit(system["snapshots"], system["u"])