working reconstructions with lag parameter

Author: nmank

pydmd/dmdc.py

@@ -191,9 +191,12 @@ 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, lag = 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
@@ -203,14 +206,14 @@ def __init__(self, svd_rank=0, tlsq_rank=0, opt=False, svd_rank_omega=-1, lag =
             "tlsq_rank": tlsq_rank,
         }
 
-        self._lag = lag
         self._opt = opt
         self._exact = False
 
         self._B = None
         self._snapshots_holder = None
         self._controlin = None
         self._basis = None
+        self._lag = lag
 
         self._modes_activation_bitmask_proxy = None
 
@@ -256,7 +259,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"
@@ -269,7 +272,10 @@ def reconstructed_data(self, control_input=None):
             [self.modes, np.diag(eigs), np.linalg.pinv(self.modes)]
         )
 
-        data = [self.snapshots[:, 0]]
+        data = []
+        for i in range(self._lag):
+            data.append(self.snapshots[:, i])
+
         expected_shape = data[0].shape
 
         for i, u in enumerate(controlin.T):
@@ -299,13 +305,18 @@ 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(
+            f"Time lag must be positive."
+                )
+        
         X = self.snapshots[:, :-self._lag]
         Y = self.snapshots[:, self._lag:]
-        n_samples = X.shape[1]+1
 
         self._set_initial_time_dictionary(
             {"t0": 0, "tend": n_samples - 1, "dt": 1}

tests/test_dmdc.py

@@ -12,16 +12,30 @@ def create_system_with_B():
     return {"snapshots": snapshots, "u": u, "B": B}
 
 
-def create_system_without_B(lag = 1):
+def create_system_without_B():
+    n = 5  # dimension snapshots
+    m = 15  # number snapshots
+    A = scipy.linalg.helmert(n, True)
+    B = np.random.rand(n, n) - 0.5
+    x0 = np.array([0.25] * n)
+    u = np.random.rand(n, m - 1) - 0.5
+    snapshots = [x0]
+    for i in range(m - 1):
+        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 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 i in range(lag):
+    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):
+    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}
@@ -42,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)
@@ -74,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)