updating

updating

.gitignore

@@ -1,23 +1,20 @@
 # Excluding datasets of exercises from commit
-/examples/Cylinder_DATASET_CFD/
-/examples/Ex_2_TR_PIV_Jet
-/examples/Ex_4_TR_PIV_Jet
-/examples/Ex_4_TR_PIV.zip
-/modulo
+#/examples/Cylinder_DATASET_CFD/
+#/examples/Ex_2_TR_PIV_Jet
+#/examples/Ex_4_TR_PIV_Jet
+#/examples/Ex_4_TR_PIV.zip
 
 # Excluding draft of Spods
-/new_functions/
+#/new_functions/
 
 # Exclude MODULO tmp
 /MODULO_tmp/
 
 # Exclude all _py_chache
 modulo/__pycache__/
-modulo/__pycache__/
-
 
 # Exclude development tests
-/tests/
+#/tests/
 
 # Extra files here and there 
 #dev_tests

README.md

@@ -25,7 +25,7 @@ For a more comprehensive overview of data driven decompositions, we refer to the
 This version expands considerably the version v1 in "Old_Python_Implementation", for which a first tutorial was provided by L. Schena in https://www.youtube.com/watch?v=y2uSvdxAwHk. 
 The major updates are the following :
 
-1. Faster EIG/SVD algorithms, using powerful randomized svd solvers from scikit_learn (see [this]() and [this]() ). It is now possible to select various options as "eig_solver" and "svd_solver", offering different trade offs in terms of accuracy vs computational time.
+1. Faster EIG/SVD algorithms, using powerful randomized svd solvers from scikit_learn (see [this](https://scikit-learn.org/stable/modules/generated/sklearn.decomposition.TruncatedSVD.html) and [this](https://scikit-learn.org/stable/modules/generated/sklearn.utils.extmath.randomized_svd.html) ). It is now possible to select various options as "eig_solver" and "svd_solver", offering different trade offs in terms of accuracy vs computational time.
 2. Faster subscale estimators for the mPOD: the previous version used the rank of the correlation matrix in each scale to define the number of modes to be computed in each portion of the splitting vector before assemblying the full basis. This is computationally very demanding. This estimation has been replaced by a frequency-based threshold (i.e based on the frequency bins within each portion), since 
 
 3. Implementation of Dynamic Mode Decomposition (DMD) from [Schmid, P.J 2010](https://www.cambridge.org/core/journals/journal-of-fluid-mechanics/article/dynamic-mode-decomposition-of-numerical-and-experimental-data/AA4C763B525515AD4521A6CC5E10DBD4)

modulo/__init__.py

@@ -0,0 +1,16 @@
+
+#from ._version import get_versions
+#__version__ = get_versions()['version']
+#del get_versions
+
+
+from .utils.read_db import *
+from .utils._data_matrix import *
+from .core._k_matrix import *
+from .utils._utils import *
+from .core._mpod_time import *
+from .core._mpod_space import *
+from .core._pod_time import *
+from .core._pod_space import *
+from .core._dft import *
+

modulo/core/_dft.py

@@ -0,0 +1,46 @@
+import os
+
+import numpy as np
+from tqdm import tqdm
+
+
+def dft_fit(N_T, F_S, D, FOLDER_OUT, SAVE_DFT=False):
+    """
+    :param N_T:
+    :param F_S:
+    :param D:
+    :param FOLDER_OUT:
+    :param SAVE_DFT:
+    :return:
+    """
+    n_t = int(N_T)
+    Freqs = np.fft.fftfreq(n_t) * F_S  # Compute the frequency bins
+    # PSI_F = np.conj(np.fft.fft(np.eye(n_t)) / np.sqrt(n_t))  # Prepare the Fourier Matrix.
+
+    # Method 1 (didactic!)
+    # PHI_SIGMA = np.dot(D, np.conj(PSI_F))  # This is PHI * SIGMA
+
+    # Method 2
+    PHI_SIGMA = (np.fft.fft(D, n_t, 1)) / (n_t ** 0.5)
+
+    PHI_F = np.zeros((D.shape[0], n_t), dtype=complex)  # Initialize the PHI_F MATRIX
+    SIGMA_F = np.zeros(n_t)  # Initialize the SIGMA_F MATRIX
+
+    # Now we proceed with the normalization. This is also intense so we time it
+    for r in tqdm(range(0, n_t)):  # Loop over the PHI_SIGMA to normalize
+        # MEX = 'Proj ' + str(r + 1) + ' /' + str(n_t)
+        # print(MEX)
+        SIGMA_F[r] = abs(np.vdot(PHI_SIGMA[:, r], PHI_SIGMA[:, r])) ** 0.5
+        PHI_F[:, r] = PHI_SIGMA[:, r] / SIGMA_F[r]
+
+    Indices = np.flipud(np.argsort(SIGMA_F))  # find indices for sorting in decreasing order
+    Sorted_Sigmas = SIGMA_F[Indices]  # Sort all the sigmas
+    Sorted_Freqs = Freqs[Indices]  # Sort all the frequencies accordingly.
+    Phi_F = PHI_F[:, Indices]  # Sorted Spatial Structures Matrix
+    SIGMA_F = Sorted_Sigmas  # Sorted Amplitude Matrix (vector)
+
+    if SAVE_DFT:
+        os.makedirs(FOLDER_OUT + 'DFT', exist_ok=True)
+        np.savez(FOLDER_OUT + 'DFT/dft_fitted', Freqs=Sorted_Freqs, Phis=Phi_F, Sigmas=SIGMA_F)
+
+    return Sorted_Freqs, Phi_F, SIGMA_F

modulo/core/_dmd_s.py

@@ -0,0 +1,84 @@
+import os
+import numpy as np
+
+# import jax.numpy as jnp
+from sklearn.decomposition import TruncatedSVD
+# For efficient linear algebra
+from numpy import linalg as LA
+# For Timing
+import time
+from ..utils._utils import switch_svds
+
+
+def dmd_s(D_1, D_2, n_Modes, F_S,
+          SAVE_T_DMD=False,
+          FOLDER_OUT='./',
+          svd_solver: str = 'svd_sklearn_truncated'):
+    """
+    This method computes the Dynamic Mode Decomposition (DMD). 
+
+    --------------------------------------------------------------------------------------------------------------------
+    Parameters:
+    ----------
+    :param D_2: np.array
+           Second portion of the data, i.e. D[:,1:n_t]
+    :param Phi_P, Psi_P, Sigma_P: np.arrays
+            POD decomposition of D1
+    :param F_S: float
+           Sampling frequency in Hz
+    :param FOLDER_OUT: str
+            Folder in which the results will be saved (if SAVE_T_DMD=True)
+    :param K: np.array
+            Temporal correlation matrix
+    :param SAVE_T_POD: bool
+            A flag deciding whether the results are saved on disk or not. If the MEMORY_SAVING feature is active, it is
+            switched True by default.
+    :param n_Modes: int
+           number of modes that will be computed
+    :param svd_solver: str,
+            svd solver to be used
+    --------------------------------------------------------------------------------------------------------------------
+    Returns:
+    --------
+    :return Phi_D: np.array
+            DMD Psis
+
+    :return Lambda_D: np.array
+            DMD Eigenvalues (of the reduced propagator)
+
+    :return freqs: np.array
+            Frequencies (in Hz, associated to the DMD modes)
+
+     :return a0s: np.array
+            Initial Coefficients of the Modes
+    """   
+
+    Phi_P, Psi_P, Sigma_P = switch_svds(D_1, n_Modes, svd_solver)
+    print('SVD of D1 rdy')
+    Sigma_inv = np.diag(1 / Sigma_P)
+    dt = 1 / F_S
+    # %% Step 3: Compute approximated propagator
+    P_A = LA.multi_dot([np.transpose(Phi_P), D_2, Psi_P, Sigma_inv])
+    print('reduced propagator rdy')
+
+    # %% Step 4: Compute eigenvalues of the system
+    Lambda, Q = LA.eig(P_A) # not necessarily symmetric def pos! Avoid eigsh, eigh
+    freqs = np.imag(np.log(Lambda)) / (2 * np.pi * dt)
+    print(' lambdas and freqs rdy')
+
+    # %% Step 5: Spatial structures of the DMD in the PIP style
+    Phi_D = LA.multi_dot([D_2, Psi_P, Sigma_inv, Q])
+    print('Phi_D rdy')
+
+    # %% Step 6: Compute the initial coefficients
+    # a0s=LA.lstsq(Phi_D, D_1[:,0],rcond=None)
+    a0s = LA.pinv(Phi_D).dot(D_1[:, 0])
+    print('Sigma_D rdy')
+
+    if SAVE_T_DMD:
+        os.makedirs(FOLDER_OUT + "/DMD/", exist_ok=True)
+        print("Saving DMD results")
+        np.savez(FOLDER_OUT + '/DMD/dmd_decomposition',
+                 Phi_D=Phi_D, Lambda=Lambda, freqs=freqs, a0s=a0s)
+
+    return Phi_D, Lambda, freqs, a0s

modulo/core/_k_matrix.py

@@ -0,0 +1,98 @@
+import os
+from tqdm import tqdm
+import numpy as np
+import math
+
+
+def CorrelationMatrix(N_T, N_PARTITIONS=1, MEMORY_SAVING=False, FOLDER_OUT='./', SAVE_K=True, D=None,weights = np.array([])):
+    """
+    This method computes the temporal correlation matrix, given a data matrix as input. It's possible to use memory saving
+    then splitting the computing in different tranches if computationally heavy. If D has been computed using MODULO
+    then the dimension dim_col and N_PARTITIONS is automatically loaded
+    --------------------------------------------------------------------------------------------------------------------
+    Parameters
+    ----------
+
+    :param N_T: int
+                Number of temporal snapshots
+    :param D: np.array
+                Data matrix
+    :param SAVE_K: bool
+                If SAVE_K=True, the matrix K is saved on disk. If the MEMORY_SAVING feature is active, this is done
+                by default.
+    :param MEMORY_SAVING: bool
+                If MEMORY_SAVING = True, the computation of the correlation matrix is done by steps. It requires the
+                data matrix to be partitioned, following algorithm in MODULO._data_processing.
+    :param FOLDER_OUT: str
+                Folder in which the temporal correlation matrix will be stored
+    :param N_PARTITIONS: int
+                Number of partitions to be read in computing the correlation matrix. If _data_processing is used to
+                partition the data matrix, this is inherited from the main class
+    :param weights: weight vector [w_i,....,w_{N_s}] where w_i = area_cell_i/area_grid
+                    Only needed if grid is non-uniform & MEMORY_SAVING== True
+    --------------------------------------------------------------------------------------------------------------------
+    Returns
+    -------
+
+    :return: K (: np.array) if the memory saving is not active. None type otherwise.
+    """
+
+    if not MEMORY_SAVING:
+        print("\n Computing Temporal correlation matrix K ...")
+        K = np.dot(D.T, D)
+        print("\n Done.")
+
+    else:
+        SAVE_K = True
+        print("\n Using Memory Saving feature...")
+        K = np.zeros((N_T, N_T))
+        dim_col = math.floor(N_T / N_PARTITIONS)
+
+        if N_T % N_PARTITIONS != 0:
+            tot_blocks_col = N_PARTITIONS + 1
+        else:
+            tot_blocks_col = N_PARTITIONS
+
+        for k in tqdm(range(tot_blocks_col)):
+
+            di = np.load(FOLDER_OUT + f"/data_partitions/di_{k + 1}.npz")['di']
+            if weights.size != 0:
+                di = np.transpose(np.transpose(di) * np.sqrt(weights))
+
+            ind_start = k * dim_col
+            ind_end = ind_start + dim_col
+
+            if (k == tot_blocks_col - 1) and (N_T - dim_col * N_PARTITIONS > 0):
+                dim_col = N_T - dim_col * N_PARTITIONS
+                ind_end = ind_start + dim_col
+
+            K[ind_start:ind_end, ind_start:ind_end] = np.dot(di.transpose(), di)
+
+            block = k + 2
+
+            while block <= tot_blocks_col:
+                dj = np.load(FOLDER_OUT + f"/data_partitions/di_{block}.npz")['di']
+                if weights.size != 0:
+                    dj = np.transpose(np.transpose(dj) * np.sqrt(weights))
+
+                ind_start_out = (block - 1) * dim_col
+                ind_end_out = ind_start_out + dim_col
+
+                if (block == tot_blocks_col) and (N_T - dim_col * N_PARTITIONS > 0):
+                    dim_col = N_T - dim_col * N_PARTITIONS
+                    ind_end_out = ind_start_out + dim_col
+                    dj = dj[:, :dim_col]
+
+                K[ind_start:ind_end, ind_start_out:ind_end_out] = np.dot(di.T, dj)
+
+                K[ind_start_out:ind_end_out, ind_start:ind_end] = K[ind_start:ind_end, ind_start_out:ind_end_out].T
+
+                block += 1
+
+                dim_col = math.floor(N_T / N_PARTITIONS)
+
+    if SAVE_K:
+        os.makedirs(FOLDER_OUT + '/correlation_matrix', exist_ok=True)
+        np.savez(FOLDER_OUT + "/correlation_matrix/k_matrix", K=K)
+
+    return K if not MEMORY_SAVING else None