GLAPPO

pydmd/cdmd.py

@@ -7,12 +7,12 @@
 
 import numpy as np
 import scipy.sparse
-from scipy.linalg import sqrtm
 
 from .dmdbase import DMDBase
 from .dmdoperator import DMDOperator
-from .snapshots import Snapshots
+from pydmd.linalg import build_linalg_module
 from .utils import compute_svd, compute_tlsq
+from .snapshots import Snapshots
 
 
 class CDMDOperator(DMDOperator):
@@ -78,6 +78,7 @@ def compute_operator(self, compressedX, compressedY, nonCompressedY):
         :rtype: numpy.ndarray, numpy.ndarray, numpy.ndarray
         """
 
+        linalg_module = build_linalg_module(compressedX)
         U, s, V = compute_svd(compressedX, svd_rank=self._svd_rank)
 
         atilde = self._least_square_operator(U, s, V, compressedY)
@@ -86,7 +87,9 @@ def compute_operator(self, compressedX, compressedY, nonCompressedY):
             # b stands for "backward"
             bU, bs, bV = compute_svd(compressedY, svd_rank=self._svd_rank)
             atilde_back = self._least_square_operator(bU, bs, bV, compressedX)
-            atilde = sqrtm(atilde.dot(np.linalg.inv(atilde_back)))
+            atilde_back_inv = linalg_module.inv(atilde_back)
+            atilde_dotted = linalg_module.dot(atilde, atilde_back_inv)
+            atilde = linalg_module.matrix_sqrt(atilde_dotted)
 
         self._Atilde = atilde
         self._compute_eigenquantities()
@@ -196,7 +199,7 @@ def _compress_snapshots(self):
         :rtype: numpy.ndarray
         """
 
-        C_shape = (self.snapshots.shape[1], self.snapshots.shape[0])
+        C_shape = (self.snapshots.shape[-1], self.snapshots.shape[-2])
         if isinstance(self.compression_matrix, np.ndarray):
             C = self.compression_matrix
         elif self.compression_matrix == "uniform":
@@ -212,29 +215,34 @@ def _compress_snapshots(self):
                 np.random.choice(*self.snapshots.shape, replace=False),
             ] = 1.0
 
+        linalg_module = build_linalg_module(self.snapshots)
+        C = linalg_module.to(self.snapshots, C)
+
         # compress the matrix
-        Y = C.dot(self.snapshots)
+        Y = linalg_module.dot(C, self.snapshots)
 
         return Y
 
-    def fit(self, X):
+    def fit(self, X, batch=False):
         """
         Compute the Dynamic Modes Decomposition to the input data.
 
         :param X: the input snapshots.
         :type X: numpy.ndarray or iterable
+        :param batch: If `True`, the first dimension is dedicated to batching.
+        :type batch: bool
         """
         self._reset()
 
-        self._snapshots_holder = Snapshots(X)
+        self._snapshots_holder = Snapshots(X, batch=batch)
         compressed_snapshots = self._compress_snapshots()
 
-        n_samples = compressed_snapshots.shape[1]
-        X = compressed_snapshots[:, :-1]
-        Y = compressed_snapshots[:, 1:]
+        n_samples = compressed_snapshots.shape[-1]
+        X = compressed_snapshots[..., :-1]
+        Y = compressed_snapshots[..., 1:]
 
         X, Y = compute_tlsq(X, Y, self._tlsq_rank)
-        self.operator.compute_operator(X, Y, self.snapshots[:, 1:])
+        self.operator.compute_operator(X, Y, self.snapshots[..., 1:])
 
         # Default timesteps
         self._set_initial_time_dictionary(

pydmd/dmd.py

@@ -2,10 +2,8 @@
 Derived module from dmdbase.py for classic dmd.
 """
 
-import numpy as np
-from scipy.linalg import pinv
-
 from .dmdbase import DMDBase
+from pydmd.linalg import assert_same_linalg_type, build_linalg_module
 from .snapshots import Snapshots
 from .utils import compute_tlsq
 
@@ -46,20 +44,22 @@ class DMD(DMDBase):
     :type tikhonov_regularization: int or float
     """
 
-    def fit(self, X):
+    def fit(self, X, batch=False):
         """
         Compute the Dynamic Modes Decomposition to the input data.
 
         :param X: the input snapshots.
         :type X: numpy.ndarray or iterable
+        :param batch: If `True`, the first dimension is dedicated to batching.
+        :type batch: bool
         """
         self._reset()
 
-        self._snapshots_holder = Snapshots(X)
+        self._snapshots_holder = Snapshots(X, batch=batch)
 
-        n_samples = self.snapshots.shape[1]
-        X = self.snapshots[:, :-1]
-        Y = self.snapshots[:, 1:]
+        n_samples = self.snapshots.shape[-1]
+        X = self.snapshots[..., :-1]
+        Y = self.snapshots[..., 1:]
 
         X, Y = compute_tlsq(X, Y, self._tlsq_rank)
         self._svd_modes, _, _ = self.operator.compute_operator(X, Y)
@@ -81,6 +81,10 @@ def predict(self, X):
         :return: one time-step ahead predicted output.
         :rtype: numpy.ndarray
         """
-        return np.linalg.multi_dot(
-            [self.modes, np.diag(self.eigs), pinv(self.modes), X]
+        assert_same_linalg_type(X, self.modes)
+
+        linalg_module = build_linalg_module(X)
+        return linalg_module.multi_dot(
+            (self.modes, linalg_module.diag_matrix(self.eigs),
+                linalg_module.pinv(self.modes), X)
         )

pydmd/dmd_modes_tuner.py

@@ -13,7 +13,6 @@ def select_modes(
     criteria,
     in_place=True,
     return_indexes=False,
-    nullify_amplitudes=False,
 ):
     """
     Select the DMD modes by using the given `criteria`.
@@ -45,10 +44,6 @@ def select_modes(
     :param bool return_indexes: If `True`, this function returns the indexes
         corresponding to DMD modes cut using the given `criteria` (default
         `False`).
-    :param bool nullify_amplitudes: If `True`, the amplitudes associated with
-        DMD modes to be removed are set to 0, therefore the number of DMD
-        modes remains constant. If `False` (default) DMD modes are actually
-        removed, therefore the number of DMD modes in the instance decreases.
     :returns: If `return_indexes` is `True`, the returned value is a tuple
         whose items are:
 

pydmd/dmdbase.py

@@ -5,11 +5,12 @@
 
 import pickle
 from copy import copy, deepcopy
-from os.path import splitext
 
 import numpy as np
+from past.utils import old_div
 
 from .dmdoperator import DMDOperator
+from pydmd.linalg import build_linalg_module
 from .utils import compute_svd
 
 
@@ -39,12 +40,13 @@ class ActivationBitmaskProxy:
     """
 
     def __init__(self, dmd_operator, amplitudes):
+        linalg_module = build_linalg_module(dmd_operator.modes)
         self._original_modes = dmd_operator.modes
-        self._original_eigs = np.atleast_1d(dmd_operator.eigenvalues)
-        self._original_amplitudes = np.atleast_1d(amplitudes)
+        self._original_eigs = linalg_module.atleast_1d(dmd_operator.eigenvalues)
+        self._original_amplitudes = linalg_module.atleast_1d(amplitudes)
 
         self.old_bitmask = None
-        self.change_bitmask(np.full(len(dmd_operator.eigenvalues), True))
+        self.change_bitmask(np.full(dmd_operator.eigenvalues.shape[-1], True))
 
     def change_bitmask(self, value):
         """
@@ -58,16 +60,9 @@ def change_bitmask(self, value):
             `bool` whose size is the same of the number of DMD modes.
         :type value: np.ndarray
         """
-
-        # apply changes made on the proxied values to the original values
-        if self.old_bitmask is not None:
-            self._original_modes[:, self.old_bitmask] = self.modes
-            self._original_eigs[self.old_bitmask] = self.eigs
-            self._original_amplitudes[self.old_bitmask] = self.amplitudes
-
-        self._modes = np.array(self._original_modes)[:, value]
-        self._eigs = np.array(self._original_eigs)[value]
-        self._amplitudes = np.array(self._original_amplitudes)[value]
+        self._modes = self._original_modes[..., value]
+        self._eigs = self._original_eigs[..., value]
+        self._amplitudes = self._original_amplitudes[..., value]
 
         self.old_bitmask = value
 
@@ -205,10 +200,12 @@ def dmd_timesteps(self):
         :return: the time intervals of the original snapshots.
         :rtype: numpy.ndarray
         """
-        return np.arange(
+        linalg_module = build_linalg_module(self.eigs)
+        return linalg_module.arange(
             self.dmd_time["t0"],
             self.dmd_time["tend"] + self.dmd_time["dt"],
             self.dmd_time["dt"],
+            device=linalg_module.device(self.snapshots),
         )
 
     @property
@@ -219,10 +216,12 @@ def original_timesteps(self):
         :return: the time intervals of the original snapshots.
         :rtype: numpy.ndarray
         """
-        return np.arange(
+        linalg_module = build_linalg_module(self.eigs)
+        return linalg_module.arange(
             self.original_time["t0"],
             self.original_time["tend"] + self.original_time["dt"],
             self.original_time["dt"],
+            device=linalg_module.device(self.snapshots),
         )
 
     @property
@@ -285,12 +284,14 @@ def dynamics(self):
             row.
         :rtype: numpy.ndarray
         """
-        temp = np.repeat(
-            self.eigs[:, None], self.dmd_timesteps.shape[0], axis=1
+        linalg_module = build_linalg_module(self.eigs)
+        temp = linalg_module.repeat(
+            self.eigs[..., None], self.dmd_timesteps.shape[0], axis=-1
+        )
+        tpow = old_div(
+            self.dmd_timesteps - self.original_time["t0"],
+            self.original_time["dt"],
         )
-        tpow = (
-            self.dmd_timesteps - self.original_time["t0"]
-        ) / self.original_time["dt"]
 
         # The new formula is x_(k+j) = \Phi \Lambda^k \Phi^(-1) x_j.
         # Since j is fixed, for a given snapshot "u" we have the following
@@ -299,7 +300,7 @@ def dynamics(self):
         # Therefore tpow must be scaled appropriately.
         tpow = self._translate_eigs_exponent(tpow)
 
-        return np.power(temp, tpow) * self.amplitudes[:, None]
+        return linalg_module.pow(temp, tpow) * self.amplitudes[..., None]
 
     def _translate_eigs_exponent(self, tpow):
         """
@@ -332,7 +333,8 @@ def reconstructed_data(self):
         :return: the matrix that contains the reconstructed snapshots.
         :rtype: numpy.ndarray
         """
-        return self.modes.dot(self.dynamics)
+        linalg_module = build_linalg_module(self.modes)
+        return linalg_module.dot(self.modes, self.dynamics)
 
     @property
     def snapshots(self):
@@ -366,7 +368,9 @@ def frequency(self):
         :return: the array that contains the frequencies of the eigenvalues.
         :rtype: numpy.ndarray
         """
-        return np.log(self.eigs).imag / (2 * np.pi * self.original_time["dt"])
+        linalg_module = build_linalg_module(self.eigs)
+        div = 2 * np.pi * self.original_time["dt"]
+        return linalg_module.log(self.eigs).imag / div
 
     @property
     def growth_rate(self):  # To check
@@ -443,7 +447,7 @@ def modes_activation_bitmask(self):
 
         bitmask = self._modes_activation_bitmask_proxy.old_bitmask
         # make sure that the array is immutable
-        bitmask.flags.writeable = False
+        build_linalg_module(bitmask).make_not_writeable(bitmask)
         return bitmask
 
     @modes_activation_bitmask.setter
@@ -452,7 +456,6 @@ def modes_activation_bitmask(self, value):
         if not self.fitted:
             raise RuntimeError("This DMD instance has not been fitted yet.")
 
-        value = np.array(value)
         if value.dtype != bool:
             raise RuntimeError(
                 "Unxpected dtype, expected bool, got {}.".format(value.dtype)
@@ -659,36 +662,31 @@ def load(fname):
 
     def _optimal_dmd_matrices(self):
         # compute the vandermonde matrix
-        vander = np.vander(self.eigs, len(self.dmd_timesteps), True)
+        linalg_module = build_linalg_module(self.eigs)
+        vander = linalg_module.vander(self.eigs, len(self.dmd_timesteps), True)
 
-        P = np.multiply(
-            np.dot(self.modes.conj().T, self.modes),
-            np.conj(np.dot(vander, vander.conj().T)),
-        )
+        a = linalg_module.dot(self.modes.conj().swapaxes(-1, -2), self.modes)
+        b = linalg_module.dot(vander, vander.conj().swapaxes(-1, -2)).conj()
+        P = linalg_module.multiply_elementwise(a, b)
 
         if self._exact:
-            q = np.conj(
-                np.diag(
-                    np.linalg.multi_dot(
-                        [vander, self.snapshots.conj().T, self.modes]
-                    )
-                )
+            vsm = linalg_module.multi_dot(
+                (vander, self.snapshots.conj().swapaxes(-1, -2), self.modes)
             )
+            q = linalg_module.extract_diagonal(vsm).conj()
         else:
-            _, s, V = compute_svd(self.snapshots[:, :-1], self.modes.shape[-1])
-
-            q = np.conj(
-                np.diag(
-                    np.linalg.multi_dot(
-                        [
-                            vander[:, :-1],
-                            V,
-                            np.diag(s).conj(),
-                            self.operator.eigenvectors,
-                        ]
-                    )
-                )
+            _, s, V = compute_svd(
+                self.snapshots[..., :-1], self.modes.shape[-1]
+            )
+
+            s_conj = linalg_module.diag_matrix(s).conj()
+            s_conj, V, vander = linalg_module.to(
+                self.operator.eigenvectors, s_conj, V, vander
+            )
+            vVse = linalg_module.multi_dot(
+                (vander[..., :-1], V, s_conj, self.operator.eigenvectors)
             )
+            q = linalg_module.extract_diagonal(vVse).conj()
 
         return P, q
 
@@ -710,20 +708,22 @@ def _compute_amplitudes(self):
         Jovanovic et al. 2014, Sparsity-promoting dynamic mode decomposition,
         https://hal-polytechnique.archives-ouvertes.fr/hal-00995141/document
         """
+        linalg_module = build_linalg_module(self.modes)
         if isinstance(self._opt, bool) and self._opt:
-            # b optimal
-            a = np.linalg.solve(*self._optimal_dmd_matrices())
+            A, b = self._optimal_dmd_matrices()
+            a = linalg_module.solve(A, b)
         else:
             if isinstance(self._opt, bool):
                 amplitudes_snapshot_index = 0
             else:
                 amplitudes_snapshot_index = self._opt
 
-            a = np.linalg.lstsq(
+            selected_snapshots = self.snapshots[..., amplitudes_snapshot_index]
+            a = linalg_module.lstsq(
                 self.modes,
-                self.snapshots.T[amplitudes_snapshot_index],
+                linalg_module.to(self.modes, selected_snapshots),
                 rcond=None,
-            )[0]
+            )
 
         return a