diff --git a/.gitignore b/.gitignore index 0be9221..6b22cb7 100644 --- a/.gitignore +++ b/.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 diff --git a/README.md b/README.md index d0d4d98..1989ab4 100755 --- a/README.md +++ b/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) diff --git a/modulo/__init__.py b/modulo/__init__.py new file mode 100644 index 0000000..6e48909 --- /dev/null +++ b/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 * + diff --git a/modulo/core/__pycache__/_dft.cpython-310.pyc b/modulo/core/__pycache__/_dft.cpython-310.pyc new file mode 100644 index 0000000..40f8998 Binary files /dev/null and b/modulo/core/__pycache__/_dft.cpython-310.pyc differ diff --git a/modulo/core/__pycache__/_dft.cpython-39.pyc b/modulo/core/__pycache__/_dft.cpython-39.pyc new file mode 100644 index 0000000..1e0f8b3 Binary files /dev/null and b/modulo/core/__pycache__/_dft.cpython-39.pyc differ diff --git a/modulo/core/__pycache__/_dmd_s.cpython-310.pyc b/modulo/core/__pycache__/_dmd_s.cpython-310.pyc new file mode 100644 index 0000000..2cc7794 Binary files /dev/null and b/modulo/core/__pycache__/_dmd_s.cpython-310.pyc differ diff --git a/modulo/core/__pycache__/_dmd_s.cpython-39.pyc b/modulo/core/__pycache__/_dmd_s.cpython-39.pyc new file mode 100644 index 0000000..b80ab16 Binary files /dev/null and b/modulo/core/__pycache__/_dmd_s.cpython-39.pyc differ diff --git a/modulo/core/__pycache__/_k_matrix.cpython-310.pyc b/modulo/core/__pycache__/_k_matrix.cpython-310.pyc new file mode 100644 index 0000000..d9ae3ce Binary files /dev/null and b/modulo/core/__pycache__/_k_matrix.cpython-310.pyc differ diff --git a/modulo/core/__pycache__/_k_matrix.cpython-39.pyc b/modulo/core/__pycache__/_k_matrix.cpython-39.pyc new file mode 100644 index 0000000..fe85d1d Binary files /dev/null and b/modulo/core/__pycache__/_k_matrix.cpython-39.pyc differ diff --git a/modulo/core/__pycache__/_mpod_space.cpython-310.pyc b/modulo/core/__pycache__/_mpod_space.cpython-310.pyc new file mode 100644 index 0000000..e4264bd Binary files /dev/null and b/modulo/core/__pycache__/_mpod_space.cpython-310.pyc differ diff --git a/modulo/core/__pycache__/_mpod_space.cpython-39.pyc b/modulo/core/__pycache__/_mpod_space.cpython-39.pyc new file mode 100644 index 0000000..57d90e9 Binary files /dev/null and b/modulo/core/__pycache__/_mpod_space.cpython-39.pyc differ diff --git a/modulo/core/__pycache__/_mpod_time.cpython-310.pyc b/modulo/core/__pycache__/_mpod_time.cpython-310.pyc new file mode 100644 index 0000000..fe332a2 Binary files /dev/null and b/modulo/core/__pycache__/_mpod_time.cpython-310.pyc differ diff --git a/modulo/core/__pycache__/_mpod_time.cpython-39.pyc b/modulo/core/__pycache__/_mpod_time.cpython-39.pyc new file mode 100644 index 0000000..419e738 Binary files /dev/null and b/modulo/core/__pycache__/_mpod_time.cpython-39.pyc differ diff --git a/modulo/core/__pycache__/_pod_space.cpython-310.pyc b/modulo/core/__pycache__/_pod_space.cpython-310.pyc new file mode 100644 index 0000000..4c3851a Binary files /dev/null and b/modulo/core/__pycache__/_pod_space.cpython-310.pyc differ diff --git a/modulo/core/__pycache__/_pod_space.cpython-39.pyc b/modulo/core/__pycache__/_pod_space.cpython-39.pyc new file mode 100644 index 0000000..9e69b2f Binary files /dev/null and b/modulo/core/__pycache__/_pod_space.cpython-39.pyc differ diff --git a/modulo/core/__pycache__/_pod_time.cpython-310.pyc b/modulo/core/__pycache__/_pod_time.cpython-310.pyc new file mode 100644 index 0000000..b801acc Binary files /dev/null and b/modulo/core/__pycache__/_pod_time.cpython-310.pyc differ diff --git a/modulo/core/__pycache__/_pod_time.cpython-39.pyc b/modulo/core/__pycache__/_pod_time.cpython-39.pyc new file mode 100644 index 0000000..4058eb1 Binary files /dev/null and b/modulo/core/__pycache__/_pod_time.cpython-39.pyc differ diff --git a/modulo/core/__pycache__/_spod_s.cpython-310.pyc b/modulo/core/__pycache__/_spod_s.cpython-310.pyc new file mode 100644 index 0000000..79f58ef Binary files /dev/null and b/modulo/core/__pycache__/_spod_s.cpython-310.pyc differ diff --git a/modulo/core/__pycache__/_spod_s.cpython-39.pyc b/modulo/core/__pycache__/_spod_s.cpython-39.pyc new file mode 100644 index 0000000..3ac79fe Binary files /dev/null and b/modulo/core/__pycache__/_spod_s.cpython-39.pyc differ diff --git a/modulo/core/__pycache__/_spod_t.cpython-310.pyc b/modulo/core/__pycache__/_spod_t.cpython-310.pyc new file mode 100644 index 0000000..113af9d Binary files /dev/null and b/modulo/core/__pycache__/_spod_t.cpython-310.pyc differ diff --git a/modulo/core/__pycache__/_spod_t.cpython-39.pyc b/modulo/core/__pycache__/_spod_t.cpython-39.pyc new file mode 100644 index 0000000..4016d28 Binary files /dev/null and b/modulo/core/__pycache__/_spod_t.cpython-39.pyc differ diff --git a/modulo/core/_dft.py b/modulo/core/_dft.py new file mode 100644 index 0000000..cca4ff8 --- /dev/null +++ b/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 diff --git a/modulo/core/_dmd_s.py b/modulo/core/_dmd_s.py new file mode 100644 index 0000000..43585c4 --- /dev/null +++ b/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 diff --git a/modulo/core/_k_matrix.py b/modulo/core/_k_matrix.py new file mode 100644 index 0000000..0be199f --- /dev/null +++ b/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 diff --git a/modulo/core/_mpod_space.py b/modulo/core/_mpod_space.py new file mode 100644 index 0000000..e206b32 --- /dev/null +++ b/modulo/core/_mpod_space.py @@ -0,0 +1,164 @@ +from scipy import signal +import numpy as np +import os +import sys +import gc +from tqdm import tqdm +from scipy.sparse import diags +from scipy.sparse.linalg import svds +from scipy.signal import firwin # To create FIR kernels +import math + + +def spatial_basis_mPOD(D, PSI_M, N_T, N_PARTITIONS, N_S, MEMORY_SAVING, FOLDER_OUT, SAVE: bool = False): + """ + Given the temporal basis of the mPOD now the spatial ones are computed + + :param D: + :param PSI_M: + :return: + """ + R1 = 0 + R2 = 0 + if MEMORY_SAVING: + SAVE = True + os.makedirs(FOLDER_OUT + '/mPOD/', exist_ok=True) + dim_col = math.floor(N_T / N_PARTITIONS) + dim_row = math.floor(N_S / N_PARTITIONS) + dr = np.zeros((dim_row, N_T)) + + # 1 --- Converting partitions dC to dR + if N_S % N_PARTITIONS != 0: + tot_blocks_row = N_PARTITIONS + 1 + else: + tot_blocks_row = N_PARTITIONS + + if N_T % N_PARTITIONS != 0: + tot_blocks_col = N_PARTITIONS + 1 + else: + tot_blocks_col = N_PARTITIONS + + fixed = 0 + + for i in range(1, tot_blocks_row + 1): + # --- Check if dim_row has to be fixed: + if i == tot_blocks_row and (N_S - dim_row * N_PARTITIONS > 0): + dim_row_fix = N_S - dim_row * N_PARTITIONS + dr = np.zeros((dim_row_fix, N_T)) + + for cont in range(1, tot_blocks_col + 1): + di = np.load(FOLDER_OUT + f"/data_partitions/di_{cont}.npz")['di'] + + if i == tot_blocks_row and (N_S - dim_row * N_PARTITIONS > 0) and fixed == 0: + R1 = R2 + R2 = R1 + (N_S - dim_row * N_PARTITIONS) + fixed = 1 + elif fixed == 0: + R1 = (i - 1) * dim_row + R2 = i * dim_row + + # Same as before, but we don't need the variable fixed because if + # % the code runs this loop, it will be the last time + + if cont == tot_blocks_col and (N_T - dim_col * N_PARTITIONS > 0): + C1 = C2 + C2 = C1 + (N_T - dim_col * N_PARTITIONS) + else: + C1 = (cont - 1) * dim_col + C2 = cont * dim_col + + dr[:, C1:C2] = di[R1:R2, :] + + # 2 --- Computing partitions R of PHI_SIGMA + PHI_SIGMA_BLOCK = dr @ PSI_M + np.savez(FOLDER_OUT + f'/mPOD/phi_sigma_{i}', PHI_SIGMA_BLOCK) + + # 3 --- Convert partitions R to partitions C and get SIGMA + R = PSI_M.shape[1] + dim_col = math.floor(R / N_PARTITIONS) + dim_row = math.floor(N_S / N_PARTITIONS) + dps = np.zeros((N_S, dim_col)) + SIGMA_M = [] + PHI_M = [] + + if R % N_PARTITIONS != 0: + tot_blocks_col = N_PARTITIONS + 1 + else: + tot_blocks_col = N_PARTITIONS + + fixed = 0 + + # Here we apply the same logic of the loop before + + for j in range(1, tot_blocks_col + 1): + + if j == tot_blocks_col and (R - dim_col * N_PARTITIONS > 0): + dim_col_fix = R - dim_col * N_PARTITIONS + dps = np.zeros((N_S, dim_col_fix)) + + for k in range(1, tot_blocks_row + 1): + PHI_SIGMA_BLOCK = np.load(FOLDER_OUT + f"/mPOD/phi_sigma_{k}.npz")['arr_0'] + + if j == tot_blocks_col and (R - dim_col * N_PARTITIONS > 0) and fixed == 0: + R1 = R2 + R2 = R1 + (R - dim_col * N_PARTITIONS) + fixed = 1 + elif fixed == 0: + R1 = (j - 1) * dim_col + R2 = j * dim_col + + if k == tot_blocks_row and (N_S - dim_row * N_PARTITIONS > 0): + C1 = C2 + C2 = C1 + (N_S - dim_row * N_PARTITIONS) + else: + C1 = (k - 1) * dim_row + C2 = k * dim_row + + dps[C1:C2, :] = PHI_SIGMA_BLOCK[:, R1:R2] + + # Getting sigmas and phis + for z in range(R1, R2): + zz = z - R1 + SIGMA_M.append(np.linalg.norm(dps[:, zz])) + tmp = dps[:, zz] / SIGMA_M[z] + #print(f'Shape tmp = {np.shape(tmp)}') + PHI_M.append(tmp) + np.savez(FOLDER_OUT + f'/mPOD/phi_{z + 1}', tmp) + + Indices = np.argsort(SIGMA_M)[::-1] # find indices for sorting in decreasing order + SIGMA_M = np.asarray(SIGMA_M) + PHI_M = np.asarray(PHI_M).T + PSI_M = np.asarray(PSI_M) + Sorted_Sigmas = SIGMA_M[Indices] # Sort all the sigmas + Phi_M = PHI_M[:, Indices] # Sorted Spatial Structures Matrix + Psi_M = PSI_M[:, Indices] # Sorted Temporal Structures Matrix + Sigma_M = Sorted_Sigmas # Sorted Amplitude Matrix + + else: + R = PSI_M.shape[1] + PHI_M_SIGMA_M = np.dot(D, (PSI_M)) + # Initialize the output + PHI_M = np.zeros((N_S, R)) + SIGMA_M = np.zeros((R)) + + for i in tqdm(range(0, R)): + # print('Completing mPOD Mode ' + str(i)) + # Assign the norm as amplitude + SIGMA_M[i] = np.linalg.norm(PHI_M_SIGMA_M[:, i]) + # Normalize the columns of C to get spatial modes + PHI_M[:, i] = PHI_M_SIGMA_M[:, i] / SIGMA_M[i] + + Indices = np.flipud(np.argsort(SIGMA_M)) # find indices for sorting in decreasing order + Sorted_Sigmas = SIGMA_M[Indices] # Sort all the sigmas + Phi_M = PHI_M[:, Indices] # Sorted Spatial Structures Matrix + Psi_M = PSI_M[:, Indices] # Sorted Temporal Structures Matrix + Sigma_M = Sorted_Sigmas # Sorted Amplitude Matrix + + if SAVE: + '''Saving results in MODULO tmp proper folder''' + os.makedirs(FOLDER_OUT + '/mPOD/', exist_ok=True) + np.savez(FOLDER_OUT + "/mPOD/sorted_phis", Phi_M) + np.savez(FOLDER_OUT + "/mPOD/sorted_psis", Psi_M) + np.savez(FOLDER_OUT + "/mPOD/sorted_sigma", Sorted_Sigmas) + + return Phi_M, Psi_M, Sigma_M diff --git a/modulo/core/_mpod_time.py b/modulo/core/_mpod_time.py new file mode 100644 index 0000000..23cee04 --- /dev/null +++ b/modulo/core/_mpod_time.py @@ -0,0 +1,180 @@ +import os + +import numpy as np +from scipy.signal import firwin # To create FIR kernels +#from scipy.sparse.linalg import svds +from tqdm import tqdm + +from modulo.utils._utils import conv_m, switch_eigs + + +def temporal_basis_mPOD(K, Nf, Ex, F_V, Keep, boundaries, MODE='reduced', dt=1, + FOLDER_OUT: str = "./", + MEMORY_SAVING: bool = False, + SAT: int = 100, + n_Modes=10, + eig_solver: str = 'svd_sklearn_randomized'): + ''' + This function computes the PSIs for the mPOD. In this implementation, a "dft-trick" is proposed, in order to avoid + expansive SVDs. Randomized SVD is used instead. + -------------------------------------------------------------------------------------------------------------------- + Parameters: + ----------- + + :param K: np.array + Temporal correlation matrix + :param dt: float + Time step + :param Nf: np.array + Filters + :param Ex: int + Extension. Must have Ex >= Nf + :param F_V: np.array + Filter array + :param Keep: np.array + vector defining which scale to keep. + :param boundaries: str + In order to avoid 'edge effects' if the time correlation matrix is not periodic, several boundary conditions + can be used. + Options are (from scipy.ndimage.convolve): + ‘reflect’ (d c b a | a b c d | d c b a) The input is extended by reflecting about the edge of the last pixel. + ‘nearest’ (a a a a | a b c d | d d d d) The input is extended by replicating the last pixel. + ‘wrap’ (a b c d | a b c d | a b c d) The input is extended by wrapping around to the opposite edge. + :param MODE: str + As a final step of this algorithm, the orthogonality is imposed via a QR-factorization. This parameter + define how to perform such factorization, according to numpy. + Options: this is a wrapper to np.linalg.qr(_, mode=MODE). + Check numpy's documentation. + if ‘reduced’ The final basis will not necessarely be full + if ‘complete’ The final basis will always be full + :param FOLDER_OUT: str + This is the directory where intermediate results will be stored if the memory saving is active + It will be ignored if MEMORY_SAVING=False. + + + :param MEMORY_SAVING: Bool + If memory saving is active, the results will be saved locally + Nevertheless, since Psi_M is usually not expensive, it will be returned. + :param SAT: int + Maximum number of modes per scale. + The user can decide how many modes to compute; otherwise, modulo set the default SAT=100. + :param n_Modes: int + Total number of modes that will be finally exported + + :param eig_solver: str + This is the eigenvalue solver that will be used. Refer to eigs_swith for the options. + + -------------------------------------------------------------------------------------------------------------------- + Returns: + -------- + :return PSI_M: np.array + mPOD PSIs + ''' + + if Ex < np.max(Nf): + raise RuntimeError("For the mPOD temporal basis computation Ex must be larger than or equal to Nf") + + '''Converting F_V in radiants and initialise number of scales M''' + Fs = 1 / dt + F_Bank_r = F_V * 2 / Fs + M = len(F_Bank_r) + + # Loop over the scales to show the transfer functions + Psi_M = np.array([]) + Lambda_M = np.array([]) + n_t = K.shape[1] + + # if K_S: + # Ks = np.zeros((n_t, n_t, M + 1)) + + '''DFT-trick below: computing frequency bins.''' + Freqs = np.fft.fftfreq(n_t) * Fs + + print("Filtering and Diagonalizing H scale: \n") + + '''Filtering and computing eigenvectors''' + + for m in tqdm(range(0, M)): + # Generate the 1d filter for this + if m < 1: + if Keep[m] == 1: + # Low Pass Filter + h_A = firwin(Nf[m], F_Bank_r[m], window='hamming') + # Filter K_LP + print('\n Filtering Largest Scale') + K_L = conv_m(K=K, h=h_A, Ex=Ex, boundaries=boundaries) + # R_K = np.linalg.matrix_rank(K_L, tol=None, hermitian=True) + '''We replace it with an estimation based on the non-zero freqs the cut off frequency of the scale is ''' + F_CUT = F_Bank_r[m] * Fs / 2 + Indices = np.argwhere(np.abs(Freqs) < F_CUT) + R_K = np.min([len(Indices), SAT]) + print(str(len(Indices)) + ' Modes Estimated') + print('\n Diagonalizing Largest Scale') + Psi_P, Lambda_P = switch_eigs(K_L, R_K, eig_solver) #svds_RND(K_L, R_K) + Psi_M=Psi_P; Lambda_M=Lambda_P + else: + print('\n Scale '+str(m)+' jumped (keep['+str(m)+']=0)') + # if K_S: + # Ks[:, :, m] = K_L # First large scale + + # method = signal.choose_conv_method(K, h2d, mode='same') + elif m > 0 and m < M - 1: + if Keep[m] == 1: + # print(m) + print('\n Working on Scale '+str(m)+'/'+str(M)) + # This is the 1d Kernel for Band pass + h1d_H = firwin(Nf[m], [F_Bank_r[m], F_Bank_r[m + 1]], pass_zero=False) # Band-pass + F_CUT1 = F_Bank_r[m] * Fs / 2 + F_CUT2 = F_Bank_r[m + 1] * Fs / 2 + Indices = np.argwhere((np.abs(Freqs) > F_CUT1) & (np.abs(Freqs) < F_CUT2)) + R_K = np.min([len(Indices), SAT]) # number of frequencies here + print(str(len(Indices)) + ' Modes Estimated') + # print('Filtering H Scale ' + str(m + 1) + '/' + str(M)) + K_H = conv_m(K, h1d_H, Ex, boundaries) + # Ks[:, :, m + 1] = K_H # Intermediate band-pass + print('Diagonalizing H Scale ' + str(m + 1) + '/' + str(M)) + # R_K = np.linalg.matrix_rank(K_H, tol=None, hermitian=True) + Psi_P, Lambda_P = switch_eigs(K_H, R_K, eig_solver) #svds_RND(K_H, R_K) # Diagonalize scale + if np.shape(Psi_M)[0]==0: # if this is the first contribute to the basis + Psi_M=Psi_P; Lambda_M=Lambda_P + else: + Psi_M = np.concatenate((Psi_M, Psi_P), axis=1) # append to the previous + Lambda_M = np.concatenate((Lambda_M, Lambda_P), axis=0) + else: + print('\n Scale '+str(m)+' jumped (keep['+str(m)+']=0)') + + else: # this is the case m=M: this is a high pass + if Keep[m] == 1: + print('Working on Scale '+str(m)+'/'+str(M)) + # This is the 1d Kernel for High Pass (last scale) + h1d_H = firwin(Nf[m], F_Bank_r[m], pass_zero=False) + F_CUT1 = F_Bank_r[m] * Fs / 2 + Indices = np.argwhere((np.abs(Freqs) > F_CUT1)) + R_K = len(Indices) + R_K = np.min([len(Indices), SAT]) # number of frequencies here + print(str(len(Indices)) + ' Modes Estimated') + print('Filtering H Scale ' + str(m + 1) + '/ ' + str(M)) + K_H = conv_m(K, h1d_H, Ex, boundaries) + # Ks[:, :, m + 1] = K_H # Last (high pass) scale + print('Diagonalizing H Scale ' + str(m + 1) + '/ ' + str(M)) + # R_K = np.linalg.matrix_rank(K_H, tol=None, hermitian=True) + Psi_P, Lambda_P = switch_eigs(K_H, R_K, eig_solver) #svds_RND(K_H, R_K) # Diagonalize scale + Psi_M = np.concatenate((Psi_M, Psi_P), axis=1) # append to the previous + Lambda_M = np.concatenate((Lambda_M, Lambda_P), axis=0) + else: + print('\n Scale '+str(m)+' jumped (keep['+str(m)+']=0)') + + # Now Order the Scales + Indices = np.flip(np.argsort(Lambda_M)) # find indices for sorting in decreasing order + Psi_M = Psi_M[:, Indices] # Sort the temporal structures + #print(f"Size psis in mpodtime = {np.shape(Psi_M)}") + # Now we complete the basis via re-orghotonalization + print('\n QR Polishing...') + PSI_M, R = np.linalg.qr(Psi_M, mode=MODE) + print('Done!') + + if MEMORY_SAVING: + os.makedirs(FOLDER_OUT + '/mPOD', exist_ok=True) + np.savez(FOLDER_OUT + '/mPOD/Psis', Psis=PSI_M) + + return PSI_M[:,0:n_Modes] diff --git a/modulo/core/_pod_space.py b/modulo/core/_pod_space.py new file mode 100644 index 0000000..b3009b5 --- /dev/null +++ b/modulo/core/_pod_space.py @@ -0,0 +1,194 @@ +import math +import os + +import numpy as np +from tqdm import tqdm + + +def Spatial_basis_POD(D, PSI_P, Sigma_P, MEMORY_SAVING, + N_T, FOLDER_OUT='./', N_PARTITIONS=1, + SAVE_SPATIAL_POD=False): + """ + Given the temporal basis now the POD spatial ones are computed + -------------------------------------------------------------------------------------------------------------------- + Parameters: + ---------- + :param N_T: int + Number of temporal snapshots + + :param FOLDER_OUT: str + Folder in which the results are saved if SAVE_SPATIAL_POD = True + + :param SAVE_SPATIAL_POD: bool + If True, results are saved on disk and released from memory + + :param N_PARTITIONS: int + Number of partitions to be loaded. If D has been partitioned using MODULO, this parameter is automatically + inherited from the main class. To be specified otherwise. + + :param MEMORY_SAVING: bool + Inherited from main class, if True turns on the MEMORY_SAVING feature, loading the partitions and starting + the proper algorithm + + :param D: np.array + Data matrix on which to project the temporal basis + + :param PSI_P: np.array + POD Psis + :param Sigma_P: np.array + POD Sigmas + -------------------------------------------------------------------------------------------------------------------- + Returns: + -------- + + :return Phi_P: np.array + POD Phis + """ + R = PSI_P.shape[1] + + if not MEMORY_SAVING: + N_S = D.shape[0] + + + #The following is the general normalization approach. + #not needed for POD. I leave it commented for the moment. + + # Phi_P = np.zeros((N_S, R)) + + # # N_S = D.shape[0] + # PHI_P_SIGMA_P = np.dot(D, PSI_P) + + # print("Completing Spatial Structures Modes: \n") + + # for i in tqdm(range(0, R)): + # # Normalize the columns of C to get spatial modes + # Phi_P[:, i] = PHI_P_SIGMA_P[:, i] / Sigma_P[i] + + # if SAVE_SPATIAL_POD: + # os.makedirs(FOLDER_OUT + 'POD', exist_ok=True) + # np.savez(FOLDER_OUT + '/POD/pod_spatial_basis', + # phis=Phi_P, PHI_P_SIGMA_P=PHI_P_SIGMA_P) + + # return Phi_P + + # We take only the first R modes. + Sigma_P_t=Sigma_P[0:R]; Sigma_P_Inv_V=1/Sigma_P_t + # So we have the inverse + Sigma_P_Inv=np.diag(Sigma_P_Inv_V) + # Here is the one shot projection: + Phi_P=np.linalg.multi_dot([D,PSI_P[:,0:R],Sigma_P_Inv]) + + return Phi_P + + + else: + + N_S = np.shape(np.load(FOLDER_OUT + "/data_partitions/di_1.npz")['di'])[0] + + dim_col = math.floor(N_T / N_PARTITIONS) + dim_row = math.floor(N_S / N_PARTITIONS) + dr = np.zeros((dim_row, N_T)) + + # 1 -- Converting partitions dC to dR + if N_S % N_PARTITIONS != 0: + tot_blocks_row = N_PARTITIONS + 1 + else: + tot_blocks_row = N_PARTITIONS + + if N_T % N_PARTITIONS != 0: + tot_blocks_col = N_PARTITIONS + 1 + else: + tot_blocks_col = N_PARTITIONS + + # --- Loading Psi_P + fixed = 0 + R1 = 0 + R2 = 0 + C1 = 0 + C2 = 0 + + for i in range(1, tot_blocks_row + 1): + + if (i == tot_blocks_row) and (N_S - dim_row * N_PARTITIONS > 0): + dim_row_fix = N_S - dim_row * N_PARTITIONS + dr = np.zeros((dim_row_fix, N_T)) + + for b in range(1, tot_blocks_col + 1): + di = np.load(FOLDER_OUT + f"/data_partitions/di_{b}.npz")['di'] + if (i == tot_blocks_row) and (N_S - dim_row * N_PARTITIONS > 0) and fixed == 0: + R1 = R2 + R2 = R1 + (N_S - dim_row * N_PARTITIONS) + fixed = 1 + elif fixed == 0: + R1 = (i - 1) * dim_row + R2 = i * dim_row + + if (b == tot_blocks_col) and (N_T - dim_col * N_PARTITIONS > 0): + C1 = C2 + C2 = C1 + (N_T - dim_col * N_PARTITIONS) + else: + C1 = (b - 1) * dim_col + C2 = b * dim_col + + np.copyto(dr[:, C1:C2], di[R1:R2, :]) + + PHI_SIGMA_BLOCK = np.dot(dr, PSI_P) + np.savez(FOLDER_OUT + f"/POD/PHI_SIGMA_{i}", + phi_sigma=PHI_SIGMA_BLOCK) + + # 3 - Converting partitions R to partitions C and get Sigmas + dim_col = math.floor(R / N_PARTITIONS) + dim_row = math.floor(N_S / N_PARTITIONS) + dps = np.zeros((N_S, dim_col)) + Phi_P = np.zeros((N_S, R)) + + if R % N_PARTITIONS != 0: + tot_blocks_col = N_PARTITIONS + 1 + else: + tot_blocks_col = N_PARTITIONS + + fixed = 0 + + for i in range(1, tot_blocks_col + 1): + + if (i == tot_blocks_col) and (R - dim_col * N_PARTITIONS > 0): + dim_col_fix = R - dim_col * N_PARTITIONS + dps = np.zeros((N_S, dim_col_fix)) + + for b in range(1, tot_blocks_row + 1): + + PHI_SIGMA_BLOCK = np.load(FOLDER_OUT + f"/POD/PHI_SIGMA_{b}.npz")['phi_sigma'] + + if (i == tot_blocks_col) and (R - dim_col * N_PARTITIONS > 0) and fixed == 0: + R1 = R2 + R2 = R1 + (R - dim_col * N_PARTITIONS) + fixed = 1 + elif fixed == 0: + R1 = (i - 1) * dim_col + R2 = i * dim_col + + if (b == tot_blocks_col) and (N_S - dim_row * N_PARTITIONS > 0): + C1 = C2 + C2 = C1 + (N_S - dim_row * N_PARTITIONS) + else: + C1 = (b - 1) * dim_row + C2 = b * dim_row + + dps[C1:C2, :] = PHI_SIGMA_BLOCK[:, R1:R2] + + # Computing Sigmas and Phis: + # TODO: np.linalg.norm(dps[:, jj]) to substitute with Sigma (that is already available from _pod_time) + for j in range(R1, R2): + jj = j - R1 + Phi_P = dps[:, jj] / np.linalg.norm(dps[:, jj]) + np.savez(FOLDER_OUT + f"/POD/phi_{j + 1}", phi_p=Phi_P) + + # Read the temporary files to build Phi_P_Matrix (Lorenzo pls fix this :D) + # TODO + Phi_P_M = np.zeros((N_S, R)) + for j in range(R): + Phi_P_V = np.load(FOLDER_OUT + f"/POD/phi_{j + 1}.npz")['phi_p'] + Phi_P_M[:, j] = Phi_P_V + + + return Phi_P_M diff --git a/modulo/core/_pod_time.py b/modulo/core/_pod_time.py new file mode 100644 index 0000000..3f4ae57 --- /dev/null +++ b/modulo/core/_pod_time.py @@ -0,0 +1,68 @@ +import os +import numpy as np +from sklearn.decomposition import TruncatedSVD +from ..utils._utils import switch_eigs +from scipy.sparse.linalg import svds + + +def Temporal_basis_POD(K, SAVE_T_POD=False, FOLDER_OUT='./', n_Modes=10, + eig_solver: str = 'eigh'): + """ + This method computes the POD basis. For some theoretical insights, you can find + the theoretical background of the proper orthogonal decomposition in a nutshell here: + + https://youtu.be/8fhupzhAR_M + -------------------------------------------------------------------------------------------------------------------- + Parameters: + ---------- + :param FOLDER_OUT: str + Folder in which the results will be saved (if SAVE_T_POD=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 throughout the computation + -------------------------------------------------------------------------------------------------------------------- + Returns: + -------- + + :return: Psi_P: np.array + POD Psis + :return: Sigma_P: np.array + POD Sigmas + """ + # Solver 1: Use the standard SVD + # Psi_P, Lambda_P, _ = np.linalg.svd(K) + # Sigma_P = np.sqrt(Lambda_P) + + # Solver 2: Use randomized SVD ############## WARNING ################# + # if svd_solver.lower() == 'svd_sklearn_truncated': + # svd = TruncatedSVD(n_Modes) + # svd.fit_transform(K) + # Psi_P = svd.components_.T + # Lambda_P = svd.singular_values_ + # Sigma_P = np.sqrt(Lambda_P) + # elif svd_solver.lower() == 'svd_numpy': + # Psi_P, Lambda_P, _ = np.linalg.svd(K) + # Sigma_P = np.sqrt(Lambda_P) + # elif svd_solver.lower() == 'svd_sklearn_randomized': + # Psi_P, Lambda_P, _ = svds_RND(K, n_Modes) + # Sigma_P = np.sqrt(Lambda_P) + # elif svd_solver.lower() == 'svd_scipy_sparse': + # Psi_P, Lambda_P, _ = svds(K, k=n_Modes) + # Sigma_P = np.sqrt(Lambda_P) + + print("diagonalizing K....") + Psi_P, Sigma_P = switch_eigs(K, n_Modes, eig_solver) + + + if SAVE_T_POD: + os.makedirs(FOLDER_OUT + "/POD/", exist_ok=True) + print("Saving POD temporal basis") + np.savez(FOLDER_OUT + '/POD/temporal_basis', Psis=Psi_P, Sigmas=Sigma_P) + + return Psi_P, Sigma_P diff --git a/modulo/core/_spod_s.py b/modulo/core/_spod_s.py new file mode 100644 index 0000000..4555f29 --- /dev/null +++ b/modulo/core/_spod_s.py @@ -0,0 +1,97 @@ +import numpy as np +#from ._k_matrix import CorrelationMatrix +from scipy import signal +from scipy.signal import firwin +from ._pod_time import Temporal_basis_POD +from ._pod_space import Spatial_basis_POD + +def compute_SPOD_s(D, K, F_S, n_s, n_t, + N_o=100, f_c=0.3, + n_Modes=10, SAVE_SPOD=True, + FOLDER_OUT='./', + MEMORY_SAVING=False, N_PARTITIONS=1): + """ + This method computes the Spectral POD of your data. + This is the one by Sieber + et al (https://www.cambridge.org/core/journals/journal-of-fluid-mechanics/article/abs/spectral-proper-orthogonal-decomposition/DCD8A6EDEFD56F5A9715DBAD38BD461A) + + :param F_S: float, + Sampling Frequency [Hz] + :param N_o: float, + Semi-Order of the diagonal filter. + Note that the filter order will be 2 N_o +1 (to make sure it is odd) + :param f_c: float, + cut-off frequency of the diagonal filter + :param n_Modes: float, + number of modes to be computed + :param SAVE_SPOD: bool, + If True, MODULO will save the output in self.FOLDER OUT/MODULO_tmp + :return Psi_P: np.array + SPOD Psis + :return Sigma_P: np.array + SPOD Sigmas. + :return Phi_P: np.array + SPOD Phis + """ + # if self.D is None: + # D = np.load(self.FOLDER_OUT + '/MODULO_tmp/data_matrix/database.npz')['D'] + # SAVE_SPOD = True + # # TODO : Lorenzo check this stuff + # else: + # D = self.D + # + # n_s = self.N_S # Repeat variable for debugging compatibility + # n_t = self.N_T + # + # print('Computing Correlation Matrix \n') + + # The first step is the same as the POD: we compute the correlation matrix + # K = CorrelationMatrix(self.N_T, self.N_PARTITIONS, self.MEMORY_SAVING, + # self.FOLDER_OUT, D=self.D) + + # 1. Initialize the extended + K_e = np.zeros((n_t + 2 * N_o, n_t + 2 * N_o)) + # From which we clearly know that: + K_e[N_o:n_t + N_o, N_o:n_t + N_o] = K + + # 2. We fill the edges ( a bit of repetition but ok.. ) + + # Row-wise, Upper part + for i in range(0, N_o): + K_e[i, i:i + n_t] = K[0, :] + + # Row-wise, bottom part + for i in range(N_o + n_t, n_t + 2 * N_o): + K_e[i, i - n_t + 1:i + 1] = K[-1, :] + + # Column-wise, left part + for j in range(0, N_o): + K_e[j:j + n_t, j] = K[:, 0] + + # Column-wise, right part + for j in range(N_o + n_t, 2 * N_o + n_t): + K_e[j - n_t + 1:j + 1, j] = K[:, -1] + + # Now you create the diagonal kernel in 2D + h_f = firwin(N_o, f_c) # Kernel in 1D + # This is also something that must be put in a separate file: + # To cancel the phase lag we make this non-causal with a symmetric + # shift, hence with zero padding as equal as possible on both sides + n_padd_l = round((n_t - N_o) / 2); + n_padd_r = n_t - N_o - n_padd_l + + h_f_pad = np.pad(h_f, (n_padd_l, n_padd_r)) # symmetrically padded kernel in 1D + h_f_2 = np.diag(h_f_pad) + + # Finally the filtered K is just + K_F = signal.fftconvolve(K_e, h_f_2, mode='same')[N_o:n_t + N_o, N_o:n_t + N_o] + # plt.plot(np.diag(K),'b--'); plt.plot(np.diag(K_F_e),'r') + + # From now on it's just POD: + Psi_P, Sigma_P = Temporal_basis_POD(K_F, SAVE_SPOD, FOLDER_OUT, n_Modes) + + Phi_P = Spatial_basis_POD(D, N_T=n_t, PSI_P=Psi_P, Sigma_P=Sigma_P, + MEMORY_SAVING=MEMORY_SAVING, FOLDER_OUT=FOLDER_OUT, + N_PARTITIONS=N_PARTITIONS) + + return Phi_P, Psi_P, Sigma_P \ No newline at end of file diff --git a/modulo/core/_spod_t.py b/modulo/core/_spod_t.py new file mode 100644 index 0000000..07937c3 --- /dev/null +++ b/modulo/core/_spod_t.py @@ -0,0 +1,103 @@ +import numpy as np +from modulo.utils._utils import overlap +#from sklearn.utils.extmath import randomized_svd +from tqdm import tqdm +import os + +from modulo.utils._utils import switch_svds + + + +def compute_SPOD_t(D, F_S, L_B=500, O_B=250, + n_Modes=10, SAVE_SPOD=True, FOLDER_OUT='/', + possible_svds='svd_sklearn_truncated'): + """ + This method computes the Spectral POD of your data. + :param D: array, + snapshot matrix to decompose, of size N_S,N_T + :param F_S: float, + Sampling Frequency [Hz] + :param L_B: float, + lenght of the chunks + :param O_B: float, + Overlapping between blocks in the chunk + :param n_Modes: float, + number of modes to be computed FOR EACH FREQUENCY + :param SAVE_SPOD: bool, + If True, MODULO will save the output in FOLDER OUT/MODULO_tmp + :param possible_svds + :return Psi_P_hat: np.array + Spectra of the SPOD Modes + :return Sigma_P: np.array + Amplitudes of the SPOD Modes. + :return Phi_P: np.array + SPOD Phis + :return freq: float + frequency bins for the Spectral POD + + + """ + + # if D is None: + # D = np.load(FOLDER_OUT + '/MODULO_tmp/data_matrix/database.npz')['D'] + # SAVE_SPOD = True + # else: + # D = D + # + # n_s = N_S # Repeat variable for debugging compatibility + # n_t = N_T + # + # # First comput the PS in each point (this is very time consuming and should be parallelized) + # # Note: this can be improved a lot...! ok for the moment + print('Computing PSD at all points\n') + N_S,N_T=np.shape(D) + + # Step 1 : Partition the data into blocks ( potentially overlapping) + Ind = np.arange(N_T) + Indices = overlap(Ind, len_chunk=L_B, len_sep=O_B) + + N_B = np.shape(Indices)[1] + N_P = np.shape(Indices)[0] + print('Partitioned into blocks of length n_B=' + str(N_B)) + print('Number of partitions retained is n_P=' + str(N_P)) + + # The frequency bins are thus defined: + Freqs = np.fft.fftfreq(N_B) * F_S # Compute the frequency bins + Keep_IND = np.where(Freqs >= 0) + N_B2 = len(Keep_IND[0]) # indexes for positive frequencies + Freqs_Pos = Freqs[Keep_IND] # positive frequencies + + # Step 2 : Construct the D_hats in each partition + D_P_hat_Tens = np.zeros((N_S, N_B, N_P)) + print('Computing DFTs in each partition') + for k in tqdm(range(0, N_P)): # Loop over the partitions + D_p = D[:, Indices[k]] # Take the portion of data + D_P_hat_Tens[:, :, k] = np.fft.fft(D_p, N_B, 1) + + # This would be the mean over the frequencies + # D_hat_Mean=np.mean(D_P_hat_Tens,axis=1) + + # Initialize the outputs + Sigma_SP = np.zeros((n_Modes, N_B2)) + Phi_SP = np.zeros((N_S, n_Modes, N_B2)) + + # Step 3: Loop over frequencies to build the modes. + # Note: you only care about half of these frequencies. + # This is why you loop over N_B2, not N_B + print('Computing POD for each frequency') + for j in tqdm(range(0, N_B2)): + # Get D_hat of the chunk + D_hat_f = D_P_hat_Tens[:, j, :] + # Go for the SVD + + U,V,Sigma=switch_svds(D_hat_f,n_Modes,svd_solver=possible_svds) + + Phi_SP[:, :, j] = U + Sigma_SP[:, j] = Sigma / (N_S * N_B) + + if SAVE_SPOD: + folder_dir = FOLDER_OUT + 'SPOD_T' + os.makedirs(folder_dir, exist_ok=True) + np.savez(folder_dir + 'spod_t.npz', Phi=Phi_SP, Sigma=Sigma_SP, Freqs=Freqs_Pos) + + return Phi_SP, Sigma_SP, Freqs_Pos diff --git a/modulo/modulo.py b/modulo/modulo.py new file mode 100644 index 0000000..13cbd25 --- /dev/null +++ b/modulo/modulo.py @@ -0,0 +1,881 @@ +# Functional ones: +import os +import numpy as np +from scipy import linalg +from sklearn.metrics.pairwise import pairwise_kernels +# To have fancy loading bar +from tqdm import tqdm + +# All the functions from the modulo package +from modulo.core._dft import dft_fit +from modulo.core._dmd_s import dmd_s +from modulo.core._k_matrix import CorrelationMatrix +from modulo.core._mpod_space import spatial_basis_mPOD +from modulo.core._mpod_time import temporal_basis_mPOD +from modulo.core._pod_space import Spatial_basis_POD +from modulo.core._pod_time import Temporal_basis_POD +from modulo.core._spod_s import compute_SPOD_s +from modulo.core._spod_t import compute_SPOD_t +from modulo.utils._data_matrix import DataMatrix +from modulo.utils._utils import switch_svds + + + + +class MODULO: + """ + MODULO (MODal mULtiscale pOd) is a software developed at the von Karman Institute to perform Multiscale + Modal Analysis of numerical and experimental data using the Multiscale Proper Orthogonal Decomposition (mPOD). + + Theoretical foundation can be found at: + https://arxiv.org/abs/1804.09646 + + Presentation of the MODULO framework available here: + https://arxiv.org/pdf/2004.12123.pdf + + YouTube channel with hands-on tutorials can be found at: + https://youtube.com/playlist?list=PLEJZLD0-4PeKW6Ze984q08bNz28GTntkR + + All the codes so far assume that the dataset is equally spaced both in space (i.e. along a Cartesian grid) + and in time. The extension to non-uniformly sampled data will be included in future releases. + + + """ + + def __init__(self, data: np.array, + N_PARTITIONS: int = 1, + FOLDER_OUT='./', + N_T: int =100, + N_S: int= 200, + n_Modes: int = 10, + dtype: str = 'float32', + eig_solver: str = 'eigh', + svd_solver: str = 'svd_sklearn_truncated', + weights: np.array = np.array([])): + """ + This function initializes the main parameters needed by MODULO. + + Attributes: + + :param data: This is the data matrix to factorize. It is a np.array with + shape ((N_S, N_T)). If the data has not yet been prepared in the form of a np.array, + the method ReadData in MODULO can be used (see ReadData). + If the memory saving is active (N_PARTITIONS >1), the folder with partitions should be prepared. + If the memory saving is active, this entry = None. The data matrix is assumed to big to be saved and the + + + :param N_PARTITIONS: If memory saving feature is active, this parameter sets the number of partitions + that will be used to store the data matrices during the computations. + + :param FOLDER_OUT: Folder in which the output will be stored. + The output includes the matrices Phi, Sigma and Psi (optional) and temporary files + used for some of the calculations (e.g.: for memory saving). + + :param N_T: Number of time steps, must be given when N_PARTITIONS >1 + + :param N_S: Number of grid points, must be given when N_PARTITIONS >1 + + :param n_Modes: Number of Modes to be computed + + :param dtype: Cast "data" with type dtype + + :param eig_solver: Numerical solver to compute the eigen values + + :param svd_solver: Numerical solver to compute the Single Value Decomposition + + :param weights: weight vector [w_i,....,w_{N_s}] where w_i = area_cell_i/area_grid + Only needed if grid is non-uniform. + + + """ + + print("MODULO (MODal mULtiscale pOd) is a software developed at the von Karman Institute to perform " + "data driven modal decomposition of numerical and experimental data. \n" + "\n") + + + if not isinstance(data, np.ndarray) and N_PARTITIONS==1: + raise TypeError("Please check that your database is in an numpy array format. If D=None, then you must have memory saving (N_PARTITIONS>1)") + + + # Load the data matrix + if isinstance(data,np.ndarray): + # Number of points in time and space + self.N_T = data.shape[1] + self.N_S = data.shape[0] + # Check the data type + self.D = data.astype(dtype) + else: + self.D = None # D is never saved when N_partitions >1 + self.N_S = N_S; # so N_S and N_t must be given as parameters of modulo + self.N_T = N_T + + + # Load and applied the weights to the D matrix + if weights.size != 0: + if len(weights) == self.N_S: + print("The weights you have input have the size of the columns of D \n" + "MODULO has considered that you have already duplicated the dimensions of the weights " + "to match the dimensions of the D columns \n") + self.weights = weights + elif 2*len(weights) == self.N_S: # 2D computation only + self.weights = np.concatenate((weights,weights)) + print("Modulo assumes you have a 2D domain and has duplicated the weight " + "array to match the size of the D columns \n") + else: + raise AttributeError("Make sure the size of the weight array is twice smaller than the size of D") + # Dstar is used to compute the K matrix + self.Dstar = np.transpose(np.transpose(self.D) * np.sqrt(self.weights)) + else: + print("Modulo assumes you have a uniform grid. " + "If not, please give the weights as parameters of MODULO!") + self.weights = weights + self.Dstar = self.D + + if N_PARTITIONS>1: + self.MEMORY_SAVING = True + else: + self.MEMORY_SAVING = False + + # Assign the number of modes + self.n_Modes = n_Modes + # If particular needs, override choice for svd and eigen solve + self.svd_solver = svd_solver.lower() + self.eig_solver = eig_solver.lower() + possible_svds = ['svd_numpy', 'svd_scipy_sparse', 'svd_sklearn_randomized', 'svd_sklearn_truncated'] + possible_eigs = ['svd_sklearn_randomized', 'eigsh', 'eigh'] + + if self.svd_solver not in possible_svds: + raise NotImplementedError("The requested SVD solver is not implemented. Please pick one of the following:" + "which belongs to: \n {}".format(possible_svds)) + + if self.eig_solver not in possible_eigs: + raise NotImplementedError("The requested EIG solver is not implemented. Please pick one of the following: " + " \n {}".format(possible_eigs)) + + # if N_PARTITIONS >= self.N_T: + # raise AttributeError("The number of requested partitions is greater of the total columns (N_T). Please," + # "try again.") + + self.N_PARTITIONS = N_PARTITIONS + + self.FOLDER_OUT = FOLDER_OUT + 'MODULO_tmp' + + if self.MEMORY_SAVING: + os.makedirs(self.FOLDER_OUT, exist_ok=True) + + + + + + + + def _data_processing(self, + MR: bool = False, + SAVE_D: bool = False): + """ + This method pre-process the data before running the factorization. + If the memory saving option is active, the method ensures the correct splitting of the data matrix + into the required partitions. + If the mean removal is desired (MR: True), this method removes the time averaged column from the + data matrix D. + If neither the mean removal nor the memory saving options are active, this method is skipped. + + :param MR: bool + if True, the mean field is removed from the data matrix. + + :param SAVE_D: bool + if True, the D matrix is saved in the folder decided by the user + + + :return: D directly (as class attribute) if Memory saving is not active. + Otherwise, it returns None and the matrix is automatically saved on disk. + """ + + self.D = DataMatrix(self.D, self.FOLDER_OUT, MEMORY_SAVING=self.MEMORY_SAVING, + N_PARTITIONS=self.N_PARTITIONS, MR=MR, SAVE_D=SAVE_D) + + + return True + + def _correlation_matrix(self, + SAVE_K: bool = True): + """ + This method computes the time correlation matrix. Here the memory saving is + beneficial for large datasets. Since the matrix could be potentially heavy, it is automatically stored on disk to + minimize the usage of the RAM. This feature can be deactivated by setting SAVE_K = False. + In this case, the correlation matrix is returned to the main class. + + :param SAVE_K: bool + A flag deciding if the matrix will be stored in the disk (in FOLDER_OUT/MODULO_tmp) or not. + Default option is 'True'. This attribute is passed to the class + in order to decide if it has to be loaded from disk or not. + + + :return K: np.array + The correlation matrix D^T D (as class attribute) if Memory saving is not active. + Otherwise, it returns None and the matrix is automatically saved on disk. + + """ + + self.SAVE_K = SAVE_K + + self.K = CorrelationMatrix(self.N_T, self.N_PARTITIONS, self.MEMORY_SAVING, + self.FOLDER_OUT, self.SAVE_K, D=self.Dstar, weights=self.weights) + + return + + def _temporal_basis_POD(self, + SAVE_T_POD: bool = True): + """ + This method computes the temporal structure for the Proper Orthogonal Decomposition (POD) computation. + The theoretical background of the POD is briefly recalled here: + + https://youtu.be/8fhupzhAR_M + + The diagonalization of K is computed via Singular Value Decomposition (SVD). + A speedup is available if the user is on Linux machine, in which case MODULO + exploits the power of JAX and its Numpy implementation. + + For more on JAX: + + https://github.com/google/jax + https://jax.readthedocs.io/en/latest/jax.numpy.html + + If the user is on a Win machine, Linux OS can be used using + the Windows Subsystem for Linux. + + For more on WSL: + https://docs.microsoft.com/en-us/windows/wsl/install-win10 + + :param SAVE_T_POD: bool + Flag deciding if the results will be stored on the disk. + Default value is True, to limit the RAM's usage. + Note that this might cause a minor slowdown for the loading, + but the tradeoff seems worthy. + This attribute is passed to the MODULO class. + + + POD temporal basis are returned if MEMORY_SAVING is not active. Otherwise all the results are saved on disk. + + :return Psi_P: np.array + POD Psis + + :return Sigma_P: np.array + POD Sigmas. If needed, Lambdas can be easily computed recalling that: Sigma_P = np.sqrt(Lambda_P) + """ + + if self.MEMORY_SAVING: + K = np.load(self.FOLDER_OUT + "/correlation_matrix/k_matrix.npz")['K'] + SAVE_T_POD = True + else: + K = self.K + + Psi_P, Sigma_P = Temporal_basis_POD(K, SAVE_T_POD, + self.FOLDER_OUT, self.n_Modes, self.eig_solver) + + del K + return Psi_P, Sigma_P if not self.MEMORY_SAVING else None + + def _spatial_basis_POD(self, Psi_P, Sigma_P, + SAVE_SPATIAL_POD: bool = True): + """ + This method computes the spatial structure for the Proper Orthogonal Decomposition (POD) computation. + The theoretical background of the POD is briefly recalled here: + + https://youtu.be/8fhupzhAR_M + + :param Psi_P: np.array + POD temporal basis + :param Sigma_P: np.array + POD Sigmas + :param SAVE_SPATIAL_POD: bool + Flag deciding if the results will be stored on the disk. + Default value is True, to limit the RAM's usage. + Note that this might cause a minor slowdown for the loading, + but the tradeoff seems worthy. + This attribute is passed to the MODULO class. + + :return Phi_P: np.array + POD Phis + + """ + + self.SAVE_SPATIAL_POD = SAVE_SPATIAL_POD + + if self.MEMORY_SAVING: + '''Loading temporal basis from disk. They're already in memory otherwise.''' + Psi_P = np.load(self.FOLDER_OUT + 'POD/temporal_basis.npz')['Psis'] + Sigma_P = np.load(self.FOLDER_OUT + 'POD/temporal_basis.npz')['Sigmas'] + + Phi_P = Spatial_basis_POD(self.D, N_T=self.N_T, PSI_P=Psi_P, Sigma_P=Sigma_P, + MEMORY_SAVING=self.MEMORY_SAVING, FOLDER_OUT=self.FOLDER_OUT, + N_PARTITIONS=self.N_PARTITIONS, SAVE_SPATIAL_POD=SAVE_SPATIAL_POD) + + return Phi_P if not self.MEMORY_SAVING else None + + def _temporal_basis_mPOD(self, K, Nf, Ex, F_V, Keep, boundaries, MODE, dt, K_S=False): + """ + This function computes the temporal structures of each scale in the mPOD, as in step 4 of the algorithm + ref: Multi-Scale Proper Orthogonal Decomposition of Complex Fluid Flows - M. A. Mendez et al. + + :param K: np.array + Temporal correlation matrix + :param Nf: np.array + Order of the FIR filters that are used to isolate each of the scales + :param Ex: int + Extension at the boundaries of K to impose the boundary conditions (see boundaries) + It must be at least as Nf. + :param F_V: np.array + Frequency splitting vector, containing the frequencies of each scale (see article). + If the time axis is in seconds, these frequencies are in Hz. + :param Keep: np.array + Scale keep + :param boundaries: str -> {'nearest', 'reflect', 'wrap' or 'extrap'} + Define the boundary conditions for the filtering process, in order to avoid edge effects. + The available boundary conditions are the classic ones implemented for image processing: + nearest', 'reflect', 'wrap' or 'extrap'. See also https://docs.scipy.org/doc/scipy/reference/tutorial/ndimage.html + :param MODE: str -> {‘reduced’, ‘complete’, ‘r’, ‘raw’} + A QR factorization is used to enforce the orthonormality of the mPOD basis, to compensate + for the non-ideal frequency response of the filters. + The option MODE from np.linalg.qr carries out this operation. + + :return PSI_M: np.array + Multiscale POD temporal basis + + """ + + if self.MEMORY_SAVING: + K = np.load(self.FOLDER_OUT + "/correlation_matrix/k_matrix.npz")['K'] + + PSI_M = temporal_basis_mPOD(K=K, Nf=Nf, Ex=Ex, F_V=F_V, Keep=Keep, boundaries=boundaries, + MODE=MODE, dt=dt, FOLDER_OUT=self.FOLDER_OUT, + n_Modes=self.n_Modes, K_S=False, + MEMORY_SAVING=self.MEMORY_SAVING, SAT=self.SAT, eig_solver=self.eig_solver) + + + return PSI_M if not self.MEMORY_SAVING else None + + def _spatial_basis_mPOD(self, D, PSI_M, SAVE): + """ + This function implements the last step of the mPOD algorithm: + completing the decomposition. Here we project from psis, to get phis and sigmas + + :param D: np.array + data matrix + :param PSI_M: np.array + temporal basis for the mPOD. Remember that it is not possible to impose both basis matrices + phis and psis: given one of the two, the other is univocally determined. + :param SAVE: bool + if True, MODULO saves the results on disk. + + :return Phi_M: np.array + mPOD Phis (Matrix of spatial structures) + :return Psi_M: np.array + mPOD Psis (Matrix of temporal structures) + :return Sigma_M: np.array + mPOD Sigmas (vector of amplitudes, i.e. the diagonal of Sigma_M) + + """ + + Phi_M, Psi_M, Sigma_M = spatial_basis_mPOD(D, PSI_M, N_T=self.N_T, N_PARTITIONS=self.N_PARTITIONS, + N_S=self.N_S, MEMORY_SAVING=self.MEMORY_SAVING, + FOLDER_OUT=self.FOLDER_OUT, + SAVE=SAVE) + + return Phi_M, Psi_M, Sigma_M + + def compute_mPOD(self, Nf, Ex, F_V, Keep, SAT, boundaries, MODE, dt, SAVE=False): + """ + This function computes the temporal structures of each scale in the mPOD, as in step 4 of the algorithm + ref: Multi-Scale Proper Orthogonal Decomposition of Complex Fluid Flows - M. A. Mendez et al. + + :param K: np.array + Temporal correlation matrix + + :param Nf: np.array + Order of the FIR filters that are used to isolate each of the scales + + :param Ex: int + Extension at the boundaries of K to impose the boundary conditions (see boundaries) + It must be at least as Nf. + + :param F_V: np.array + Frequency splitting vector, containing the frequencies of each scale (see article). + If the time axis is in seconds, these frequencies are in Hz. + + :param Keep: np.array + Scale keep + + :param boundaries: str -> {'nearest', 'reflect', 'wrap' or 'extrap'} + Define the boundary conditions for the filtering process, in order to avoid edge effects. + The available boundary conditions are the classic ones implemented for image processing: + nearest', 'reflect', 'wrap' or 'extrap'. See also https://docs.scipy.org/doc/scipy/reference/tutorial/ndimage.html + + :param MODE: str -> {‘reduced’, ‘complete’, ‘r’, ‘raw’} + A QR factorization is used to enforce the orthonormality of the mPOD basis, to compensate + for the non-ideal frequency response of the filters. + The option MODE from np.linalg.qr carries out this operation. + + :param SAT: Maximum number of modes per scale. + Only used for mPOD (max number of modes per scale) + + :param dt: float + temporal step + + :return Phi_M: np.array + mPOD Phis (Matrix of spatial structures) + :return Psi_M: np.array + mPOD Psis (Matrix of temporal structures) + :return Sigma_M: np.array + mPOD Sigmas (vector of amplitudes, i.e. the diagonal of Sigma_M + + """ + + + print('Computing correlation matrix D matrix...') + self.K = CorrelationMatrix(self.N_T, self.N_PARTITIONS, + self.MEMORY_SAVING, + self.FOLDER_OUT, D=self.Dstar) + + if self.MEMORY_SAVING: + self.K = np.load(self.FOLDER_OUT + '/correlation_matrix/k_matrix.npz')['K'] + + + print("Computing Temporal Basis...") + + PSI_M = temporal_basis_mPOD(K=self.K, Nf=Nf, Ex=Ex, F_V=F_V, Keep=Keep, boundaries=boundaries, + MODE=MODE, dt=dt, FOLDER_OUT=self.FOLDER_OUT, + n_Modes=self.n_Modes,MEMORY_SAVING=self.MEMORY_SAVING, SAT=SAT, eig_solver=self.eig_solver) + + print("Done.") + + if hasattr(self, 'D'):# if self.D is available: + print('Computing Phi from D...') + Phi_M, Psi_M, Sigma_M = spatial_basis_mPOD(self.D, PSI_M, N_T=self.N_T, N_PARTITIONS=self.N_PARTITIONS, + N_S=self.N_S, MEMORY_SAVING=self.MEMORY_SAVING, + FOLDER_OUT=self.FOLDER_OUT, + SAVE=SAVE) + + else: # if not, the memory saving is on and D will not be used. We pass a dummy D + print('Computing Phi from partitions...') + Phi_M, Psi_M, Sigma_M = spatial_basis_mPOD(np.array([1]), PSI_M, N_T=self.N_T, N_PARTITIONS=self.N_PARTITIONS, + N_S=self.N_S, MEMORY_SAVING=self.MEMORY_SAVING, + FOLDER_OUT=self.FOLDER_OUT, + SAVE=SAVE) + + print("Done.") + + + + return Phi_M, Psi_M, Sigma_M + + def compute_POD_K(self, SAVE_T_POD: bool = True): + """ + This method computes the Proper Orthogonal Decomposition (POD) of a dataset + using the snapshot approach, i.e. working on the temporal correlation matrix. + The eig solver for K is defined in 'eig_solver' + The theoretical background of the POD is briefly recalled here: + + https://youtu.be/8fhupzhAR_M + + :return Psi_P: np.array + POD Psis + + :return Sigma_P: np.array + POD Sigmas. If needed, Lambdas can be easily computed recalling that: Sigma_P = np.sqrt(Lambda_P) + + :return Phi_P: np.array + POD Phis + """ + + + print('Computing correlation matrix D matrix...') + self.K = CorrelationMatrix(self.N_T, self.N_PARTITIONS, + self.MEMORY_SAVING, + self.FOLDER_OUT, D=self.Dstar, weights=self.weights) + + if self.MEMORY_SAVING: + self.K = np.load(self.FOLDER_OUT + '/correlation_matrix/k_matrix.npz')['K'] + + print("Computing Temporal Basis...") + Psi_P, Sigma_P = Temporal_basis_POD(self.K, SAVE_T_POD, + self.FOLDER_OUT, self.n_Modes, eig_solver=self.eig_solver) + print("Done.") + print("Computing Spatial Basis...") + + if hasattr(self, 'D'):# if self.D is available: + print('Computing Phi from D...') + Phi_P = Spatial_basis_POD(self.D, N_T=self.N_T, + PSI_P=Psi_P, + Sigma_P=Sigma_P, + MEMORY_SAVING=self.MEMORY_SAVING, + FOLDER_OUT=self.FOLDER_OUT, + N_PARTITIONS=self.N_PARTITIONS) + + else: # if not, the memory saving is on and D will not be used. We pass a dummy D + print('Computing Phi from partitions...') + Phi_P = Spatial_basis_POD(np.array([1]), N_T=self.N_T, + PSI_P=Psi_P, + Sigma_P=Sigma_P, + MEMORY_SAVING=self.MEMORY_SAVING, + FOLDER_OUT=self.FOLDER_OUT, + N_PARTITIONS=self.N_PARTITIONS) + + + print("Done.") + + return Phi_P, Psi_P, Sigma_P + + def compute_POD_svd(self, SAVE_T_POD: bool = True): + """ + This method computes the Proper Orthogonal Decomposition (POD) of a dataset + using the SVD decomposition. The svd solver is defined by 'svd_solver'. + Note that in this case, the memory saving option is of no help, since + the SVD must be performed over the entire dataset. + + https://youtu.be/8fhupzhAR_M + + :return Psi_P: np.array + POD Psis + + :return Sigma_P: np.array + POD Sigmas. If needed, Lambdas can be easily computed recalling that: Sigma_P = np.sqrt(Lambda_P) + + :return Phi_P: np.array + POD Phis + """ + # If Memory saving is active, we must load back the data. + # This process is memory demanding. Different SVD solver will handle this differently. + + if self.MEMORY_SAVING: + if self.N_T % self.N_PARTITIONS != 0: + tot_blocks_col = self.N_PARTITIONS + 1 + else: + tot_blocks_col = self.N_PARTITIONS + + # Prepare the D matrix again + D = np.zeros((self.N_S, self.N_T)) + R1 = 0 + + # print(' \n Reloading D from tmp...') + for k in tqdm(range(tot_blocks_col)): + di = np.load(self.FOLDER_OUT + f"/data_partitions/di_{k + 1}.npz")['di'] + R2 = R1 + np.shape(di)[1] + D[:, R1:R2] = di + R1 = R2 + + # Now that we have D back, we can proceed with the SVD approach + Phi_P, Psi_P, Sigma_P=switch_svds(D, self.n_Modes, self.svd_solver) + + + else: # self.MEMORY_SAVING: + Phi_P, Psi_P, Sigma_P=switch_svds(self.D, self.n_Modes, self.svd_solver) + + + return Phi_P, Psi_P, Sigma_P + + def compute_DMD_PIP(self, SAVE_T_DMD: bool = True, F_S=1): + """ + This method computes the Dynamic Mode Decomposition of the data + using the algorithm in https://arxiv.org/abs/1312.0041. + See also https://arxiv.org/abs/2001.01971 for more details. + + :return Phi_D: np.array + DMD Phis. As for the DFT, these are complex. + + :return Lambda_D: np.array + DMD Eigenvalues (of the reduced propagator). These are complex. + + :return freqs: np.array + Frequencies (in Hz, associated to the DMD modes) + + :return a0s: np.array + Initial Coefficients of the Modes + + """ + + # If Memory saving is active, we must load back the data + if self.MEMORY_SAVING: + if self.N_T % self.N_PARTITIONS != 0: + tot_blocks_col = self.N_PARTITIONS + 1 + else: + tot_blocks_col = self.N_PARTITIONS + + # Prepare the D matrix again + D = np.zeros((self.N_S, self.N_T)) + R1 = 0 + + # print(' \n Reloading D from tmp...') + for k in tqdm(range(tot_blocks_col)): + di = np.load(self.FOLDER_OUT + f"/data_partitions/di_{k + 1}.npz")['di'] + R2 = R1 + np.shape(di)[1] + D[:, R1:R2] = di + R1 = R2 + + + # Compute the DMD + Phi_D, Lambda, freqs, a0s = dmd_s(D[:, 0:self.N_T - 1], + D[:, 1:self.N_T], self.n_Modes, F_S,svd_solver=self.svd_solver) + + else: + Phi_D, Lambda, freqs, a0s = dmd_s(self.D[:, 0:self.N_T - 1], + self.D[:, 1:self.N_T], self.n_Modes, F_S, SAVE_T_DMD=SAVE_T_DMD, + svd_solver=self.svd_solver, FOLDER_OUT=self.FOLDER_OUT) + + return Phi_D, Lambda, freqs, a0s + + def compute_DFT(self, F_S, SAVE_DFT=False): + """ + This method computes the Discrete Fourier Transform of your data. + + :param F_S: float, + Sampling Frequency [Hz] + :param SAVE_DFT: bool, + If True, MODULO will save the output in self.FOLDER OUT/MODULO_tmp + + :return: Sorted_Freqs: np.array, + Sorted Frequencies + :return Phi_F: np.array, + DFT Phis + :return Sigma_F: np.array, + DFT Sigmas + """ + if self.D is None: + D = np.load(self.FOLDER_OUT + '/MODULO_tmp/data_matrix/database.npz')['D'] + SAVE_DFT = True + Sorted_Freqs, Phi_F, SIGMA_F = dft_fit(self.N_T, F_S, D, self.FOLDER_OUT, SAVE_DFT=SAVE_DFT) + + else: + Sorted_Freqs, Phi_F, SIGMA_F = dft_fit(self.N_T, F_S, self.D, self.FOLDER_OUT, SAVE_DFT=SAVE_DFT) + + return Sorted_Freqs, Phi_F, SIGMA_F + + def compute_SPOD_t(self, F_S, L_B=500, O_B=250, n_Modes=10, SAVE_SPOD=True): + """ + This method computes the Spectral POD of your data. This is the one by Towne et al + (https://www.cambridge.org/core/journals/journal-of-fluid-mechanics/article/abs/spectral-proper-orthogonal-decomposition-and-its-relationship-to-dynamic-mode-decomposition-and-resolvent-analysis/EC2A6DF76490A0B9EB208CC2CA037717) + + :param F_S: float, + Sampling Frequency [Hz] + :param L_B: float, + lenght of the chunks + :param O_B: float, + Overlapping between blocks in the chunk + :param n_Modes: float, + number of modes to be computed for each frequency + :param SAVE_SPOD: bool, + If True, MODULO will save the output in self.FOLDER OUT/MODULO_tmp + :return Psi_P_hat: np.array + Spectra of the SPOD Modes + :return Sigma_P: np.array + Amplitudes of the SPOD Modes. + :return Phi_P: np.array + SPOD Phis + :return freq: float + frequency bins for the Spectral POD + + + """ + if self.D is None: + D = np.load(self.FOLDER_OUT + '/MODULO_tmp/data_matrix/database.npz')['D'] + Phi_SP, Sigma_SP, Freqs_Pos = compute_SPOD_t(D, F_S, L_B=L_B, O_B=O_B, + n_Modes=n_Modes, SAVE_SPOD=SAVE_SPOD, + FOLDER_OUT=self.FOLDER_OUT,possible_svds=self.svd_solver) + else: + Phi_SP, Sigma_SP, Freqs_Pos = compute_SPOD_t(self.D, F_S, L_B=L_B, O_B=O_B, + n_Modes=n_Modes, SAVE_SPOD=SAVE_SPOD, + FOLDER_OUT=self.FOLDER_OUT,possible_svds=self.svd_solver) + + + return Phi_SP, Sigma_SP, Freqs_Pos + + # New Decomposition: SPOD f + + def compute_SPOD_s(self, F_S, N_O=100, f_c=0.3, n_Modes=10, SAVE_SPOD=True): + """ + This method computes the Spectral POD of your data. + This is the one by Sieber + et al (https://www.cambridge.org/core/journals/journal-of-fluid-mechanics/article/abs/spectral-proper-orthogonal-decomposition/DCD8A6EDEFD56F5A9715DBAD38BD461A) + + :param F_S: float, + Sampling Frequency [Hz] + :param N_o: float, + Semi-Order of the diagonal filter. + Note that the filter order will be 2 N_o +1 (to make sure it is odd) + :param f_c: float, + cut-off frequency of the diagonal filter + :param n_Modes: float, + number of modes to be computed + :param SAVE_SPOD: bool, + If True, MODULO will save the output in self.FOLDER OUT/MODULO_tmp + :return Psi_P: np.array + SPOD Psis + :return Sigma_P: np.array + SPOD Sigmas. + :return Phi_P: np.array + SPOD Phis + """ + + if self.D is None: + D = np.load(self.FOLDER_OUT + '/MODULO_tmp/data_matrix/database.npz')['D'] + + self.K = CorrelationMatrix(self.N_T, self.N_PARTITIONS, self.MEMORY_SAVING, + self.FOLDER_OUT, D=D) + + Phi_sP, Psi_sP, Sigma_sP = compute_SPOD_s(D, self.K, F_S, self.N_S, self.N_T, N_O, f_c, + n_Modes, SAVE_SPOD, self.FOLDER_OUT, self.MEMORY_SAVING, + self.N_PARTITIONS) + + else: + self.K = CorrelationMatrix(self.N_T, self.N_PARTITIONS, self.MEMORY_SAVING, + self.FOLDER_OUT, D=self.D) + + Phi_sP, Psi_sP, Sigma_sP = compute_SPOD_s(self.D, self.K, F_S, self.N_S, self.N_T, N_O, f_c, + n_Modes, SAVE_SPOD, self.FOLDER_OUT, self.MEMORY_SAVING, + self.N_PARTITIONS) + + + # if self.D is None: + # D = np.load(self.FOLDER_OUT + '/MODULO_tmp/data_matrix/database.npz')['D'] + # SAVE_SPOD = True + # # TODO : Lorenzo check this stuff + # else: + # D = self.D + # + # n_s = self.N_S # Repeat variable for debugging compatibility + # n_t = self.N_T + # + # print('Computing Correlation Matrix \n') + # + # # The first step is the same as the POD: we compute the correlation matrix + # K = CorrelationMatrix(self.N_T, self.N_PARTITIONS, self.MEMORY_SAVING, + # self.FOLDER_OUT, D=self.D) + # + # # 1. Initialize the extended + # K_e = np.zeros((n_t + 2 * N_o, n_t + 2 * N_o)) + # # From which we clearly know that: + # K_e[N_o:n_t + N_o, N_o:n_t + N_o] = K + # + # # 2. We fill the edges ( a bit of repetition but ok.. ) + # + # # Row-wise, Upper part + # for i in range(0, N_o): + # K_e[i, i:i + n_t] = K[0, :] + # + # # Row-wise, bottom part + # for i in range(N_o + n_t, n_t + 2 * N_o): + # K_e[i, i - n_t + 1:i + 1] = K[-1, :] + # + # # Column-wise, left part + # for j in range(0, N_o): + # K_e[j:j + n_t, j] = K[:, 0] + # + # # Column-wise, right part + # for j in range(N_o + n_t, 2 * N_o + n_t): + # K_e[j - n_t + 1:j + 1, j] = K[:, -1] + # + # # Now you create the diagonal kernel in 2D + # h_f = firwin(N_o, f_c) # Kernel in 1D + # # This is also something that must be put in a separate file: + # # To cancel the phase lag we make this non-causal with a symmetric + # # shift, hence with zero padding as equal as possible on both sides + # n_padd_l = round((n_t - N_o) / 2); + # n_padd_r = n_t - N_o - n_padd_l + # + # h_f_pad = np.pad(h_f, (n_padd_l, n_padd_r)) # symmetrically padded kernel in 1D + # h_f_2 = np.diag(h_f_pad) + # + # # Finally the filtered K is just + # K_F = signal.fftconvolve(K_e, h_f_2, mode='same')[N_o:n_t + N_o, N_o:n_t + N_o] + # # plt.plot(np.diag(K),'b--'); plt.plot(np.diag(K_F_e),'r') + # + # # From now on it's just POD: + # Psi_P, Sigma_P = Temporal_basis_POD(K_F, SAVE_SPOD, + # self.FOLDER_OUT, self.n_Modes) + # + # Phi_P = Spatial_basis_POD(self.D, N_T=self.N_T, PSI_P=Psi_P, Sigma_P=Sigma_P, + # MEMORY_SAVING=self.MEMORY_SAVING, FOLDER_OUT=self.FOLDER_OUT, + # N_PARTITIONS=self.N_PARTITIONS) + + return Phi_sP, Psi_sP, Sigma_sP + + + + def compute_kPOD(self, M_DIST=[1,10],k_m=0.1, cent='y', n_Modes=10,alpha=1e-6): + """ + This function implements the kernel PCA as described in + https://arxiv.org/pdf/2208.07746.pdf. + + :param M_DIST: array, + position of the two snapshots that will be considered to + estimate the minimal k. They should be the most different ones. + :param k_m: float, + minimum value for the kernelized correlation + :param alpha: float + regularization for K_zeta + :param cent: string, + if 'y', the matrix K is centered. + :param n_Modes: float, + number of modes to be computed + :param SAVE_SPOD: bool, + If True, MODULO will save the output in self.FOLDER OUT/MODULO_tmp + :return Psi_P: np.array + SPOD Psis + :return Sigma_P: np.array + SPOD Sigmas. + :return Phi_P: np.array + SPOD Phis + + """ + if self.D is None: + D = np.load(self.FOLDER_OUT + '/MODULO_tmp/data_matrix/database.npz')['D'] + else: + D = self.D + + # Compute Eucledean distances + i,j=M_DIST; n_s,n_t=np.shape(D) + M_ij=np.linalg.norm(D[:,i]-D[:,j])**2 + + gamma=-np.log(k_m)/M_ij + + K_zeta=pairwise_kernels(D.T,metric='rbf',gamma=gamma) + print('Kernel K ready') + + # Compute the Kernel Matrix + n_t=np.shape(D)[1] + # Center the Kernel Matrix (if cent is 'y'): + if cent =='y': + H=np.eye(n_t)-1/n_t*np.ones_like(K_zeta) + K_zeta=H@K_zeta@H.T + print('K_zeta centered') + # Diagonalize and Sort + lambdas,Psi_xi=linalg.eigh(K_zeta+alpha*np.eye(n_t),subset_by_index=[n_t-n_Modes, n_t-1]) + lambdas,Psi_xi=lambdas[::-1],Psi_xi[:,::-1] ; Sigma_xi=np.sqrt(lambdas); + print('K_zeta diagonalized') + # Encode + # Z_xi=np.diag(Sigma_xi)@Psi_xi.T + # We compute the spatial structures as projections of the data + # onto the Psi_xi! + R = Psi_xi.shape[1] + PHI_xi_SIGMA_xi = np.dot(D, (Psi_xi)) + # Initialize the output + PHI_xi = np.zeros((n_s, R)) + SIGMA_xi = np.zeros((R)) + + for i in tqdm(range(0, R)): + # print('Completing mPOD Mode ' + str(i)) + # Assign the norm as amplitude + SIGMA_xi[i] = np.linalg.norm(PHI_xi_SIGMA_xi[:, i]) + # Normalize the columns of C to get spatial modes + PHI_xi[:, i] = PHI_xi_SIGMA_xi[:, i] / SIGMA_xi[i] + + Indices = np.flipud(np.argsort(SIGMA_xi)) # find indices for sorting in decreasing order + Sorted_Sigmas = SIGMA_xi[Indices] # Sort all the sigmas + Phi_xi = PHI_xi[:, Indices] # Sorted Spatial Structures Matrix + Psi_xi = Psi_xi[:, Indices] # Sorted Temporal Structures Matrix + Sigma_xi = Sorted_Sigmas # Sorted Amplitude Matrix + print('Phi_xi computed') + return Phi_xi, Psi_xi, Sigma_xi + + diff --git a/modulo/utils/_data_matrix.py b/modulo/utils/_data_matrix.py new file mode 100644 index 0000000..6d4998b --- /dev/null +++ b/modulo/utils/_data_matrix.py @@ -0,0 +1,84 @@ +import math +import os +import numpy as np + + +def DataMatrix(database: np.array, FOLDER_OUT: str, + MEMORY_SAVING: bool = False, N_PARTITIONS: int =1, + MR: bool = False, SAVE_D: bool = False): + """ + This method performs pre-processing operations on the data matrix, D: + - Mean removing: if MR=True, the mean (per each column - i.e.: snapshot at time t_i) is removed; + - Splitting: if the MEMORY_SAVING=True the data matrix is splitted to optimize memory usage. Moreover, D is dumped + on disk and removed from the live memory. Finally, if in this condition, also the data type of the + matrix is self is changed: from float64 -> float32, with the same purpose. + + -------------------------------------------------------------------------------------------------------------------- + Parameters + ---------- + + :param database: np.array + data matrix D + :param FOLDER_OUT: str + folder in which the data (partitions and/or data matrix itself) will be eventually saved. + :param MEMORY_SAVING: bool, optional + If True, memory saving feature is activated. Passed through __init__ + :param N_PARTITIONS: int + In memory saving environment, this parameter refers to the number of partitions to be applied + to the data matrix. If the number indicated by the user is not a multiple of the N_T + i.e.: if (N_T % N_PARTITIONS) !=0 - then an additional partition is introduced, that contains + the remaining columns + :param MR: bool, optional + If True, it removes the mean (per column) from each snapshot + :param SAVE_D: bool, optional + If True, the matrix D is saved into memory. If the Memory Saving feature is active, this is performed + by default. + + Returns + ------- + :return: D if memory saving feature is active. Otherwise, D matrix is saved on disk and this method returns None + """ + N_S = int(np.shape(database)[0]) + N_T = int(np.shape(database)[1]) + + if MR: + '''Removing mean from data matrix''' + + print("Removing the mean from D ...") + D_MEAN = np.mean(database, 1) # Temporal average (along the columns) + D_Mr = database - np.array([D_MEAN, ] * N_T).transpose() # Mean Removed + print("Computing the mean-removed D ... ") + np.copyto(database, D_Mr) + del D_Mr + + if MEMORY_SAVING: + '''Converting D into float32, applying partitions and saving all.''' + SAVE_D = True + database = database.astype('float32', casting='same_kind') + os.makedirs(FOLDER_OUT + "/data_partitions/", exist_ok=True) + print("Memory Saving feature is active. Partitioning Data Matrix...") + if N_T % N_PARTITIONS != 0: + dim_col = math.floor(N_T / N_PARTITIONS) + + columns_to_part = dim_col * N_PARTITIONS + splitted_tmp = np.hsplit(database[:, :columns_to_part], N_PARTITIONS) + for ii in range(1, len(splitted_tmp) + 1): + np.savez(FOLDER_OUT + f"/data_partitions/di_{ii}", di=splitted_tmp[ii - 1]) + + np.savez(FOLDER_OUT + f"/data_partitions/di_{N_PARTITIONS + 1}", + di=database[:, columns_to_part:]) + else: + splitted_tmp = np.hsplit(database, N_PARTITIONS) + for ii in range(1, len(splitted_tmp) + 1): + np.savez(FOLDER_OUT + f"/data_partitions/di_{ii}", di=splitted_tmp[ii - 1]) + + print("\n Data Matrix has been successfully splitted. \n") + + if SAVE_D: + '''Saving data matrix in FOLDER_OUT''' + + os.makedirs(FOLDER_OUT + "./data_matrix", exist_ok=True) + print(f"Saving the matrix D in {FOLDER_OUT}") + np.savez(FOLDER_OUT + 'data_matrix/database', D=database, n_t=N_T, n_s=N_S) + + return database if not MEMORY_SAVING else None diff --git a/modulo/utils/_plots.py b/modulo/utils/_plots.py new file mode 100644 index 0000000..32dac62 --- /dev/null +++ b/modulo/utils/_plots.py @@ -0,0 +1,52 @@ +import numpy as np +import matplotlib.pyplot as plt +from matplotlib.patches import Rectangle + + +def plot_mesh(X, Y, x, y, X_mid, Y_mid): + """ + Plot original mesh and rectangular mesh over original mesh + :param X: + :param Y: + :param x: + :param y: + :param X_mid: + :param Y_mid: + :return: + """ + # Plot original mesh + fig, ax = plt.subplots() + ax.scatter(X, Y, marker='.', color='k') + ax.set_title('Original Mesh') + ax.set_xlabel('x') + ax.set_ylabel('y') + + # Plot rectangular mesh over original mesh + fig, ax = plt.subplots(figsize=(10, 8)) + im = ax.imshow(np.zeros_like(X), cmap='gray', extent=[x.min(), x.max(), y.min(), y.max()], origin='lower') + ax.scatter(X, Y, marker='.', color='k') + ax.plot(X_mid, Y_mid, '--k', lw=1) + ax.plot(X_mid.T, Y_mid.T, '--k', lw=1) + ax.set_title('Rectangular Mesh over Original Mesh') + ax.set_xlabel('x') + ax.set_ylabel('y') + fig.colorbar(im) + ax.set(xlim=(-1.2, 1.2), ylim=(0.9, 5.3)) + + return plt.show() + + +def plot_areas(X_mid, Y_mid, areas): + """ + Plot relative areas of sub-rectangles + :param X_mid: + :param Y_mid: + :param areas: + :return: + """ + fig, ax = plt.subplots() + ax.pcolormesh(X_mid, Y_mid, areas, cmap='viridis', shading='auto') + ax.set_xlabel('x') + ax.set_ylabel('y') + ax.set_title('Relative areas of sub-rectangles') + return plt.show() \ No newline at end of file diff --git a/modulo/utils/_utils.py b/modulo/utils/_utils.py new file mode 100644 index 0000000..143c5ce --- /dev/null +++ b/modulo/utils/_utils.py @@ -0,0 +1,340 @@ +import numpy as np +from scipy import signal +from scipy.sparse.linalg import svds, eigsh +from sklearn.decomposition import TruncatedSVD +from scipy.linalg import eigh +from sklearn.utils.extmath import randomized_svd + + +def Bound_EXT(S, Ex, boundaries): + """ + This function computes the extension of a signal for + filtering purposes + + :param S: The Input signal + :param Nf: The Size of the Kernel (must be an odd number!) + :param boundaries: The type of extension: + ‘reflect’ (d c b a | a b c d | d c b a) The input is extended by reflecting about the edge of the last pixel. + ‘nearest’ (a a a a | a b c d | d d d d) The input is extended by replicating the last pixel. + ‘wrap’ (a b c d | a b c d | a b c d) The input is extended by wrapping around to the opposite edge. + ‘extrap’ Extrapolation The input is extended via linear extrapolation. + + + """ + # We first perform a zero padding + # Ex=int((Nf-1)/2) # Extension on each size + size_Ext = 2 * Ex + len(S) # Compute the size of the extended signal + S_extend = np.zeros((int(size_Ext))) # Initialize extended signal + S_extend[Ex:int((size_Ext - Ex))] = S; # Assign the Signal on the zeroes + + if boundaries == "reflect": + LEFT = np.flip(S[0:Ex]) # Prepare the reflection on the left + RIGHT = np.flip(S[len(S) - Ex:len(S)]) # Prepare the reflectino on the right + S_extend[0:Ex] = LEFT + S_extend[len(S_extend) - Ex:len(S_extend)] = RIGHT + elif boundaries == "nearest": + LEFT = np.ones(Ex) * S[0] # Prepare the constant on the left + RIGHT = np.ones(Ex) * S[len(S) - 1] # Prepare the constant on the Right + S_extend[0:Ex] = LEFT + S_extend[len(S_extend) - Ex:len(S_extend)] = RIGHT + elif boundaries == "wrap": + LEFT = S[len(S) - Ex:len(S)] # Wrap on the Left + RIGHT = S[0:Ex] # Wrap on the Right + S_extend[0:Ex] = LEFT + S_extend[len(S_extend) - Ex:len(S_extend)] = RIGHT + elif boundaries == "extrap": + LEFT = np.ones(Ex) * S[0] # Prepare the constant on the left + RIGHT = np.ones(Ex) * S[len(S) - 1] # Prepare the constant on the Right + S_extend[0:Ex] = LEFT + S_extend[len(S_extend) - Ex:len(S_extend)] = RIGHT + print('Not active yet, replaced by nearest') + return S_extend + + +def conv_m(K, h, Ex, boundaries): + """ + This function computes the 2D convolution by perfoming 2 sets of 1D convolutions. + Moreover, we here use the fft with an appropriate extension + that avoids the periodicity condition. + + :param K: Matrix to be filtered + :param h: The 1D Kernel of the filter + :param boundaries: The type of extension: + ‘reflect’ (d c b a | a b c d | d c b a) The input is extended by reflecting about the edge of the last pixel. + ‘nearest’ (a a a a | a b c d | d d d d) The input is extended by replicating the last pixel. + ‘wrap’ (a b c d | a b c d | a b c d) The input is extended by wrapping around to the opposite edge. + ‘extrap’ Extrapolation The input is extended via linear extrapolation. + """ + # Filter along the raws + n_t = np.shape(K)[0] + # Ex=int(n_t/2) + K_F1 = np.zeros(np.shape(K)) + K_F2 = np.zeros(np.shape(K)) + # K_F=np.zeros(np.shape(K)) + for k in range(0, n_t): + S = K[:, k] + S_Ext = Bound_EXT(S, Ex, boundaries) + S_Filt = signal.fftconvolve(S_Ext, h, mode='valid') + # Compute where to take the signal + Ex1 = int((len(S_Filt) - len(S)) / 2) + # K_F1[k,:]=S_Filt[Ex:(len(S_Filt)-Ex)] + K_F1[:, k] = S_Filt[Ex1:(len(S_Filt) - Ex1)] + for k in range(0, n_t): + S = K_F1[k, :] + S_Ext = Bound_EXT(S, Ex, boundaries) + S_Filt = signal.fftconvolve(S_Ext, h, mode='valid') + # Compute where to take the signal + Ex1 = int((len(S_Filt) - len(S)) / 2) + # K_F2[:,k]=S_Filt[Ex:(len(S_Filt)-Ex)] + K_F2[k, :] = S_Filt[Ex1:(len(S_Filt) - Ex1)] + # K_F=K_F1+K_F2 + return K_F2 + + +def _loop_gemm(a, b, c=None, chunksize=100): + size_i = a.shape[0] + size_zip = a.shape[1] + + size_j = b.shape[1] + size_alt_zip = b.shape[0] + + if size_zip != size_alt_zip: + ValueError("Loop GEMM zip index is not of the same size for both tensors") + + if c is None: + c = np.zeros((size_i, size_j)) + + istart = 0 + for i in range(int(np.ceil(size_i / float(chunksize)))): + + left_slice = slice(istart, istart + chunksize) + left_view = a[left_slice] + + jstart = 0 + for j in range(int(np.ceil(size_j / float(chunksize)))): + right_slice = slice(jstart, jstart + chunksize) + right_view = b[:, right_slice] + + c[left_slice, right_slice] = np.dot(left_view, right_view) + jstart += chunksize + + istart += chunksize + + return c + + +def svds_RND(K, R_K): + """ + Quick and dirty implementation of randomized SVD + for computing eigenvalues of K. We follow same input/output structure + as for the svds in scipy + """ + svd = TruncatedSVD(R_K) + svd.fit_transform(K) + Psi_P = svd.components_.T + Lambda_P = svd.singular_values_ + return Psi_P, Lambda_P + + +def overlap(array, len_chunk, len_sep=1): + """ + Returns a matrix of all full overlapping chunks of the input `array`, with a chunk + length of `len_chunk` and a separation length of `len_sep`. Begins with the first full + chunk in the array. + This function is taken from https://stackoverflow.com/questions/38163366/split-list-into-separate-but-overlapping-chunks + it is designed to split an array with certain overlaps + + :param array: + :param len_chunk: + :param len_sep: + :return array_matrix: + """ + + n_arrays = int(np.ceil((array.size - len_chunk + 1) / len_sep)) + + array_matrix = np.tile(array, n_arrays).reshape(n_arrays, -1) + + columns = np.array(((len_sep * np.arange(0, n_arrays)).reshape(n_arrays, -1) + np.tile( + np.arange(0, len_chunk), n_arrays).reshape(n_arrays, -1)), dtype=np.intp) + + rows = np.array((np.arange(n_arrays).reshape(n_arrays, -1) + np.tile( + np.zeros(len_chunk), n_arrays).reshape(n_arrays, -1)), dtype=np.intp) + + return array_matrix[rows, columns] + + +# def switch_svds_K(A, n_modes, svd_solver): +# """ +# Utility function to switch between different svd solvers +# for the diagonalization of K. Being K symmetric and positive definite, +# Its eigenvalue decomposition is equivalent to an SVD. +# Thus we are using SVD solvers as eig solvers here. +# The options are the same used for switch_svds (which goes for the SVD of D) + +# -------------------------------------------------------------------------------------------------------------------- +# Parameters: +# ----------- +# :param A: np.array, +# Array of which compute the SVD +# :param n_modes: int, +# Number of modes to be computed. Note that if the `svd_numpy` method is chosen, the full matrix are +# computed, but then only the first n_modes are returned. Thus, it is not more computationally efficient. +# :param svd_solver: str, +# can be: +# 'svd_numpy'. +# This uses np.linalg.svd. +# It is the most accurate but the slowest and most expensive. +# It computes all the modes. + +# 'svd_sklearn_truncated' +# This uses the TruncatedSVD from scikitlearn. This uses either +# svds from scipy or randomized from sklearn depending on the size +# of the matrix. These are the two remaining options. +# To the merit of chosing this is to let sklearn take the decision as to +# what to use. + +# 'svd_scipy_sparse' +# This uses the svds from scipy. + +# 'svd_sklearn_randomized', +# This uses the randomized from sklearn. +# Returns +# -------- +# :return Psi_P, np.array (N_S x n_modes) +# :return Sigma_P, np.array (n_modes) +# """ +# if svd_solver.lower() == 'svd_sklearn_truncated': +# svd = TruncatedSVD(n_modes) +# svd.fit_transform(A) +# Psi_P = svd.components_.T +# Lambda_P = svd.singular_values_ +# Sigma_P = np.sqrt(Lambda_P) +# elif svd_solver.lower() == 'svd_numpy': +# Psi_P, Lambda_P, _ = np.linalg.svd(A) +# Sigma_P = np.sqrt(Lambda_P) +# Psi_P = Psi_P[:, :n_modes] +# Sigma_P = Sigma_P[:n_modes] +# elif svd_solver.lower() == 'svd_sklearn_randomized': +# Psi_P, Lambda_P = randomized_svd(A, n_modes) +# Sigma_P = np.sqrt(Lambda_P) +# elif svd_solver.lower() == 'svd_scipy_sparse': +# Psi_P, Lambda_P, _ = svds(A, k=n_modes) +# Sigma_P = np.sqrt(Lambda_P) + +# return Psi_P, Sigma_P + + +def switch_svds(A, n_modes, svd_solver='svd_sklearn_truncated'): + """ + Utility function to switch between different svd solvers + for the SVD of the snapshot matrix. This is a true SVD solver. + + -------------------------------------------------------------------------------------------------------------------- + Parameters: + ----------- + :param A: np.array, + Array of which compute the SVD + :param n_modes: int, + Number of modes to be computed. Note that if the `svd_numpy` method is chosen, the full matrix are + computed, but then only the first n_modes are returned. Thus, it is not more computationally efficient. + :param svd_solver: str, + can be: + 'svd_numpy'. + This uses np.linalg.svd. + It is the most accurate but the slowest and most expensive. + It computes all the modes. + + 'svd_sklearn_truncated' + This uses the TruncatedSVD from scikitlearn. This uses either + svds from scipy or randomized from sklearn depending on the size + of the matrix. These are the two remaining options. + The merit of chosing this is to let sklearn take the decision as to + what to use. One prefers to force the usage of any of those with the other two + options + + 'svd_scipy_sparse' + This uses the svds from scipy. + + 'svd_sklearn_randomized', + This uses the randomized from sklearn. + + Returns + -------- + :return Psi_P, np.array (N_S x n_modes) + :return Sigma_P, np.array (n_modes) + """ + if svd_solver.lower() == 'svd_sklearn_truncated': + svd = TruncatedSVD(n_modes) + X_transformed = svd.fit_transform(A) + Phi_P = X_transformed / svd.singular_values_ + Psi_P = svd.components_.T + Sigma_P = svd.singular_values_ + elif svd_solver.lower() == 'svd_numpy': + Phi_P, Sigma_P, Psi_P = np.linalg.svd(A) + Phi_P = Phi_P[:, 0:n_modes] + Psi_P = Psi_P.T[:, 0:n_modes] + Sigma_P = Sigma_P[0:n_modes] + elif svd_solver.lower() == 'svd_sklearn_randomized': + Phi_P, Sigma_P, Psi_P = randomized_svd(A, n_modes) + Psi_P = Psi_P.T + elif svd_solver.lower() == 'svd_scipy_sparse': + Phi_P, Sigma_P, Psi_P = svds(A, k=n_modes) + Psi_P = Psi_P.T + # It turns out that this does not rank them in decreasing order. + # Hence we do it manually: + idx = np.flip(np.argsort(Sigma_P)) + Sigma_P = Sigma_P[idx] + Phi_P = Phi_P[:, idx] + Psi_P = Psi_P[:, idx] + + return Phi_P, Psi_P, Sigma_P + + +def switch_eigs(A, n_modes, eig_solver): + """ + Utility function to switch between different eig solvers in a consistent way across different + methods of the package. + -------------------------------------------------------------------------------------------------------------------- + Parameters: + ----------- + :param A: np.array, + Array of which compute the eigenvalues + :param n_modes: int, + Number of modes to be computed. Note that if the `svd_numpy` method is chosen, the full matrix are + computed, but then only the first n_modes are returned. Thus, it is not more computationally efficient. + :param eig_solver: str, + can be: + 'svd_sklearn_randomized', + This uses svd truncated approach, which picks either randomized svd or scipy svds. + By default, it should pick mostly the first. + + 'eigsh' from scipy sparse. This is a compromise between the previous and the following. + + 'eigh' from scipy lin alg. This is the most precise, although a bit more expensive + + Returns + -------- + :return Psi_P, np.array (N_S x n_modes) + :return Sigma_P, np.array (n_modes) + """ + if eig_solver.lower() == 'svd_sklearn_randomized': + Psi_P, Lambda_P = svds_RND(A, n_modes) + elif eig_solver.lower() == 'eigh': + n = np.shape(A)[0] + Lambda_P, Psi_P = eigh(A, subset_by_index=[n - n_modes, n - 1]) + # It turns out that this does not rank them in decreasing order. + # Hence we do it manually: + idx = np.flip(np.argsort(Lambda_P)) + Lambda_P = Lambda_P[idx] + Psi_P = Psi_P[:, idx] + elif eig_solver.lower() == 'eigsh': + Lambda_P, Psi_P = eigsh(A, k=n_modes) + # It turns out that this does not rank them in decreasing order. + # Hence we do it manually: + idx = np.flip(np.argsort(Lambda_P)) + Lambda_P = Lambda_P[idx] + Psi_P = Psi_P[:, idx] + + Sigma_P = np.sqrt(Lambda_P) + + return Psi_P, Sigma_P \ No newline at end of file diff --git a/modulo/utils/others.py b/modulo/utils/others.py new file mode 100644 index 0000000..bf8553d --- /dev/null +++ b/modulo/utils/others.py @@ -0,0 +1,450 @@ +# -*- coding: utf-8 -*- +""" +Created on Mon Dec 30 20:33:42 2019 +@author: mendez +""" + +import numpy as np +import matplotlib.pyplot as plt + + +def Plot_Field_TEXT_JET(File): + """ + This function plots the vector field from the TR-PIV in Exercise 4. + + :param File: Name of the file to load + + """ + # We first perform a zero padding + Name = File + # Read data from a file + DATA = np.genfromtxt(Name) # Here we have the four colums + Dat = DATA[1:, :] # Here we remove the first raw, containing the header + nxny = Dat.shape[0] # is the to be doubled at the end we will have n_s=2 * n_x * n_y + n_s = 2 * nxny + ## 1. Reconstruct Mesh from file + X_S = Dat[:, 0]; + Y_S = Dat[:, 1]; + # Reshape also the velocity components + V_X = Dat[:, 2] # U component + V_Y = Dat[:, 3] # V component + # Put both components as fields in the grid + Xg, Yg, Vxg, Vyg, Magn = Plot_Field_JET(X_S, Y_S, V_X, V_Y, False, 2, 0.6) + # Show this particular step + fig, ax = plt.subplots(figsize=(8, 5)) # This creates the figure + # Or you can plot it as streamlines + plt.contourf(Xg, Yg, Magn) + # One possibility is to use quiver + STEPx = 2 + STEPy = 2 + plt.quiver(Xg[::STEPx, ::STEPy], Yg[::STEPx, ::STEPy], \ + Vxg[::STEPx, ::STEPy], -Vyg[::STEPx, ::STEPy], color='k') # Create a quiver (arrows) plot + + plt.rc('text', usetex=True) # This is Miguel's customization + plt.rc('font', family='serif') + plt.rc('xtick', labelsize=16) + plt.rc('ytick', labelsize=16) + ax.set_aspect('equal') # Set equal aspect ratio + ax.set_xlabel('$x[mm]$', fontsize=16) + ax.set_ylabel('$y[mm]$', fontsize=16) + ax.set_title('Velocity Field via TR-PIV', fontsize=18) + ax.set_xticks(np.arange(0, 40, 10)) + ax.set_yticks(np.arange(10, 30, 5)) + ax.set_xlim([0, 35]) + ax.set_ylim(10, 29) + ax.invert_yaxis() # Invert Axis for plotting purpose + plt.show() + Name[len(Name) - 12:len(Name)] + ' Plotted' + return n_s, Xg, Yg, Vxg, -Vyg, X_S, Y_S + + +def Plot_Field_JET(X_S, Y_S, V_X, V_Y, PLOT, Step, Scale): + # Number of n_X/n_Y from forward differences + GRAD_X = np.diff(X_S) + # GRAD_Y=np.diff(Y_S); + # Depending on the reshaping performed, one of the two will start with + # non-zero gradient. The other will have zero gradient only on the change. + IND_X = np.where(GRAD_X != 0); + DAT = IND_X[0]; + n_y = DAT[0] + 1; + nxny = X_S.shape[0] # is the to be doubled at the end we will have n_s=2 * n_x * n_y + # n_s=2*nxny + # Reshaping the grid from the data + n_x = (nxny // (n_y)) # Carefull with integer and float! + Xg = np.transpose(X_S.reshape((n_x, n_y))) + Yg = np.transpose(Y_S.reshape((n_x, n_y))) # This is now the mesh! 60x114. + Mod = np.sqrt(V_X ** 2 + V_Y ** 2) + Vxg = np.transpose(V_X.reshape((n_x, n_y))) + Vyg = np.transpose(V_Y.reshape((n_x, n_y))) + Magn = np.transpose(Mod.reshape((n_x, n_y))) + if PLOT: + # Show this particular step + # fig, ax = plt.subplots(figsize=(8, 5)) # This creates the figure + # Or you can plot it as streamlines + fig, ax = plt.subplots(figsize=(8, 5)) # This creates the figure + ax.contourf(Xg, Yg, Magn) + + # One possibility is to use quiver + STEPx = Step + STEPy = Step + + ax.quiver(Xg[::STEPx, ::STEPy], Yg[::STEPx, ::STEPy], \ + Vxg[::STEPx, ::STEPy], -Vyg[::STEPx, ::STEPy], color='k', + scale=Scale) # Create a quiver (arrows) plot + + plt.rc('text', usetex=True) # This is Miguel's customization + plt.rc('font', family='serif') + plt.rc('xtick', labelsize=16) + plt.rc('ytick', labelsize=16) + # ax.set_aspect('equal') # Set equal aspect ratio + # ax.set_xlabel('$x[mm]$',fontsize=16) + # ax.set_ylabel('$y[mm]$',fontsize=16) + # ax.set_title('Velocity Field via TR-PIV',fontsize=18) + # ax.set_xticks(np.arange(0,40,10)) + # ax.set_yticks(np.arange(10,30,5)) + # ax.set_xlim([0,35]) + # ax.set_ylim(10,29) + # ax.invert_yaxis() # Invert Axis for plotting purpose + # plt.show() + + return Xg, Yg, Vxg, Vyg, Magn + + +# Define the function to produce a gif file +def Animation_JET(Giff_NAME,D,X_S,Y_S,In,Fin,Step): + """ + The gif file is created from the provided data snapshot + """ + n_t=Fin-In + GRAD_X = np.diff(X_S) + IND_X = np.where(GRAD_X != 0); + DAT = IND_X[0]; + n_y = DAT[0] + 1; + nxny = X_S.shape[0] # is the to be doubled at the end we will have n_s=2 * n_x * n_y + # n_s=2*nxny + # Reshaping the grid from the data + n_x = (nxny // (n_y)) # Carefull with integer and float! + Xg = np.transpose(X_S.reshape((n_x, n_y))) + Yg = np.transpose(Y_S.reshape((n_x, n_y))) # This is now the mesh! 60x114. n_y, n_x = Yg.shape; + nxny = n_x * n_y + import os + # Create a temporary file to store the images in the GIF + Fol_Out = 'Gif_Images_temporary' + if not os.path.exists(Fol_Out): + os.mkdir(Fol_Out) + + # Prepare Animation of the analytical solution + for k in range(1,n_t,Step): + Dat = D[:, k+In] + V_X_m = Dat[0:nxny] + V_Y_m = Dat[nxny::] + # Put both components as fields in the grid + Mod = np.sqrt(V_X_m ** 2 + V_Y_m ** 2) + Vxg = np.transpose(V_X_m.reshape((n_x, n_y))) + Vyg = np.transpose(V_Y_m.reshape((n_x, n_y))) + Magn = np.transpose(Mod.reshape((n_x, n_y))) + fig, ax = plt.subplots(figsize=(8, 5)) # This creates the figure + # Or you can plot it as streamlines + plt.contourf(Xg, Yg, Magn) + # One possibility is to use quiver + STEPx = 2 + STEPy = 2 + plt.quiver(Xg[::STEPx, ::STEPy], Yg[::STEPx, ::STEPy], \ + Vxg[::STEPx, ::STEPy], -Vyg[::STEPx, ::STEPy], color='k') # Create a quiver (arrows) plot + + plt.rc('text', usetex=True) # This is Miguel's customization + plt.rc('font', family='serif') + plt.rc('xtick', labelsize=16) + plt.rc('ytick', labelsize=16) + ax.set_aspect('equal') # Set equal aspect ratio + ax.set_xlabel('$x[mm]$', fontsize=16) + ax.set_ylabel('$y[mm]$', fontsize=16) + #ax.set_title('Velocity Field via TR-PIV', fontsize=18) + ax.set_xticks(np.arange(0, 40, 10)) + ax.set_yticks(np.arange(10, 30, 5)) + ax.set_xlim([0, 35]) + ax.set_ylim(10, 29) + plt.clim(0, 6) + ax.invert_yaxis() # Invert Axis for plotting purpose + #plt.show() + + NameOUT = Fol_Out + os.sep + 'Im%03d' % (k) + '.png' + plt.savefig(NameOUT, dpi=100) + plt.close(fig) + print('Image n ' + str(k) + ' of ' + str(n_t)) + + # Assembly the GIF + import imageio # This used for the animation + + images = [] + + for k in range(1,n_t,Step): + MEX = 'Preparing Im ' + str(k) + print(MEX) + NameOUT = Fol_Out + os.sep + 'Im%03d' % (k) + '.png' + images.append(imageio.imread(NameOUT)) + + # Now we can assembly the video and clean the folder of png's (optional) + imageio.mimsave(Giff_NAME, images, duration=0.05) + import shutil # nice and powerfull tool to delete a folder and its content + + shutil.rmtree(Fol_Out) + return 'Gif Created' + + + + +def Plot_2D_CFD_Cyl(Xg,Yg,U,V,k=10,CL=16,Name=''): + # Make a 2D plot of the 2D cylinder test case in Openfoam. + n_x,n_y=np.shape(Xg) + U_g=U[:,k].reshape(n_y,n_x).T + V_g=V[:,k].reshape(n_y,n_x).T + # Prepare the plot + fig, ax = plt.subplots(figsize=(6, 3)) # This creates the figure + plt.contourf(Xg,Yg,np.sqrt(U_g**2+V_g**2),30) + # plt.quiver(Xg,Yg,U_g,V_g,scale=10000) + ax.set_aspect('equal') # Set equal aspect ratio + ax.set_xlabel('$x[mm]$',fontsize=13) + ax.set_ylabel('$y[mm]$',fontsize=13) + #ax.set_title('Tutorial 2: Cylinder Wake',fontsize=12) + ax.set_xticks(np.arange(-0.1,0.2,0.05)) + ax.set_yticks(np.arange(-0.1,0.1,0.05)) + ax.set_xlim([-0.05,0.2]) + ax.set_ylim(-0.05,0.05) + + if CL !=0: + plt.clim(0, CL) + + circle = plt.Circle((0,0),0.0075,fill=True,color='r',edgecolor='k',alpha=0.5) + plt.gcf().gca().add_artist(circle) + plt.tight_layout() + plt.show() + + if len(Name) !=0: + plt.savefig(Name, dpi=200) + plt.close(fig) + print('Image exported') + + return + + +def Animation_2D_CFD_Cyl(Giff_NAME,D,Xg,Yg,In,Fin,Step): + """ + The gif file is created from the provided data snapshot + """ + n_t=Fin-In + n_x,n_y=np.shape(Xg); nxny=n_x*n_y + U = D[0:nxny,:] + V = D[nxny::,] + + import os + # Create a temporary file to store the images in the GIF + Fol_Out = 'Gif_Images_temporary' + if not os.path.exists(Fol_Out): + os.mkdir(Fol_Out) + # Loop to produce the Gifs + for k in range(1,n_t,Step): + NameOUT = Fol_Out + os.sep + 'Im%03d' % (k) + '.png' + Plot_2D_CFD_Cyl(Xg,Yg,U,V,k=k+In,CL=16,Name=NameOUT) + + import imageio # This used for the animation + images = [] + + for k in range(1,n_t,Step): + MEX = 'Preparing Im ' + str(k) + print(MEX) + NameOUT = Fol_Out + os.sep + 'Im%03d' % (k) + '.png' + images.append(imageio.imread(NameOUT)) + + # Now we can assembly the video and clean the folder of png's (optional) + imageio.mimsave(Giff_NAME, images, duration=0.05) + import shutil # nice and powerfull tool to delete a folder and its content + + shutil.rmtree(Fol_Out) + return 'Gif Created' + + + +def Plot_Field_TEXT_Cylinder(File,Name_Mesh,Name_FIG): + """ + This function plots the vector field from the TR-PIV in Exercise 4. + + :param File: Name of the file to load + + """ + # We first perform a zero padding + Name=File + # Read data from a file + DATA = np.genfromtxt(Name) # Here we have the four colums + Dat=DATA[1:,:] # Here we remove the first raw, containing the header + nxny=Dat.shape[0] # is the to be doubled at the end we will have n_s=2 * n_x * n_y + n_s=2*nxny + ## 1. Reconstruct Mesh from file Name_Mesh + DATA_mesh=np.genfromtxt(Name_Mesh); + Dat_mesh=DATA_mesh[1:,:] + X_S=Dat_mesh[:,0]; + Y_S=Dat_mesh[:,1]; + # Reshape also the velocity components + V_X=Dat[:,0] # U component + V_Y=Dat[:,1] # V component + # Put both components as fields in the grid + Xg,Yg,Vxg,Vyg,Magn=Plot_Field_Cylinder(X_S,Y_S,V_X,V_Y,False,2,0.6) + # Show this particular step + fig, ax = plt.subplots(figsize=(5, 3)) # This creates the figure + # Or you can plot it as streamlines + CL=plt.contourf(Xg,Yg,Magn,levels=np.arange(0,18,2)) + # One possibility is to use quiver + STEPx=1; STEPy=1 + plt.quiver(Xg[::STEPx,::STEPy],Yg[::STEPx,::STEPy],\ + Vxg[::STEPx,::STEPy],Vyg[::STEPx,::STEPy],color='k') + plt.rc('text', usetex=True) + plt.rc('font', family='serif') + plt.rc('xtick',labelsize=12) + plt.rc('ytick',labelsize=12) + fig.colorbar(CL,pad=0.05,fraction=0.025) + ax.set_aspect('equal') # Set equal aspect ratio + ax.set_xlabel('$x[mm]$',fontsize=13) + ax.set_ylabel('$y[mm]$',fontsize=13) + #ax.set_title('Tutorial 2: Cylinder Wake',fontsize=12) + ax.set_xticks(np.arange(0,70,10)) + ax.set_yticks(np.arange(-10,11,10)) + ax.set_xlim([0,50]) + ax.set_ylim(-10,10) + circle = plt.Circle((0,0),2.5,fill=True,color='r',edgecolor='k',alpha=0.5) + plt.gcf().gca().add_artist(circle) + plt.tight_layout() + plt.savefig(Name_FIG, dpi=200) + plt.show() + print(Name_FIG+' printed') + return n_s, Xg, Yg, Vxg, Vyg, X_S, Y_S + + +def Plot_Field_Cylinder(X_S,Y_S,V_X,V_Y,PLOT,Step,Scale): + # Number of n_X/n_Y from forward differences + GRAD_X=np.diff(Y_S) + #GRAD_Y=np.diff(Y_S); + # Depending on the reshaping performed, one of the two will start with + # non-zero gradient. The other will have zero gradient only on the change. + IND_X=np.where(GRAD_X!=0); DAT=IND_X[0]; n_y=DAT[0]+1; + nxny=X_S.shape[0] # is the to be doubled at the end we will have n_s=2 * n_x * n_y + #n_s=2*nxny + # Reshaping the grid from the data + n_x=(nxny//(n_y)) # Carefull with integer and float! + Xg=np.transpose(X_S.reshape((n_x,n_y))) + Yg=np.transpose(Y_S.reshape((n_x,n_y))) # This is now the mesh! 60x114. + Mod=np.sqrt(V_X**2+V_Y**2) + Vxg=np.transpose(V_X.reshape((n_x,n_y))) + Vyg=np.transpose(V_Y.reshape((n_x,n_y))) + Magn=np.transpose(Mod.reshape((n_x,n_y))) + + if PLOT: + # One possibility is to use quiver + STEPx=Step; STEPy=Step + fig, ax = plt.subplots(figsize=(5, 3)) # This creates the figure + # Or you can plot it as streamlines + CL=plt.contourf(Xg,Yg,Magn,levels=np.arange(0,18,2)) + # One possibility is to use quiver + STEPx=1; STEPy=1 + plt.quiver(Xg[::STEPx,::STEPy],Yg[::STEPx,::STEPy],\ + Vxg[::STEPx,::STEPy],Vyg[::STEPx,::STEPy],color='k') + plt.rc('text', usetex=True) + plt.rc('font', family='serif') + plt.rc('xtick',labelsize=12) + plt.rc('ytick',labelsize=12) + fig.colorbar(CL,pad=0.05,fraction=0.025) + ax.set_aspect('equal') # Set equal aspect ratio + ax.set_xlabel('$x[mm]$',fontsize=13) + ax.set_ylabel('$y[mm]$',fontsize=13) + #ax.set_title('Tutorial 2: Cylinder Wake',fontsize=12) + ax.set_xticks(np.arange(0,70,10)) + ax.set_yticks(np.arange(-10,11,10)) + ax.set_xlim([0,50]) + ax.set_ylim(-10,10) + circle = plt.Circle((0,0),2.5,fill=True,color='r',edgecolor='k',alpha=0.5) + plt.gcf().gca().add_artist(circle) + plt.tight_layout() + plt.show() + + return Xg,Yg,Vxg,Vyg,Magn + +def Plot_Scalar_Field_Cylinder(X_S,Y_S,V_X,V_Y,Scalar,PLOT,Step,Scale): + # Number of n_X/n_Y from forward differences + GRAD_X=np.diff(Y_S) + #GRAD_Y=np.diff(Y_S); + # Depending on the reshaping performed, one of the two will start with + # non-zero gradient. The other will have zero gradient only on the change. + IND_X=np.where(GRAD_X!=0); DAT=IND_X[0]; n_y=DAT[0]+1; + nxny=X_S.shape[0] # is the to be doubled at the end we will have n_s=2 * n_x * n_y + #n_s=2*nxny + # Reshaping the grid from the data + n_x=(nxny//(n_y)) # Carefull with integer and float! + Xg=np.transpose(X_S.reshape((n_x,n_y))) + Yg=np.transpose(Y_S.reshape((n_x,n_y))) # This is now the mesh! 60x114. + Vxg=np.transpose(V_X.reshape((n_x,n_y))) + Vyg=np.transpose(V_Y.reshape((n_x,n_y))) + Magn=np.transpose(Scalar.reshape((n_x,n_y))) + + if PLOT: + # One possibility is to use quiver + STEPx=Step; STEPy=Step + fig, ax = plt.subplots(figsize=(5, 3)) # This creates the figure + # Or you can plot it as streamlines + CL=plt.contourf(Xg,Yg,Magn*100,levels=np.arange(0,70,10)) + # One possibility is to use quiver + STEPx=1; STEPy=1 + plt.quiver(Xg[::STEPx,::STEPy],Yg[::STEPx,::STEPy],\ + Vxg[::STEPx,::STEPy],Vyg[::STEPx,::STEPy],color='k') + plt.rc('text', usetex=True) + plt.rc('font', family='serif') + plt.rc('xtick',labelsize=12) + plt.rc('ytick',labelsize=12) + fig.colorbar(CL,pad=0.05,fraction=0.025) + ax.set_aspect('equal') # Set equal aspect ratio + ax.set_xlabel('$x[mm]$',fontsize=13) + ax.set_ylabel('$y[mm]$',fontsize=13) + #ax.set_title('Tutorial 2: Cylinder Wake',fontsize=12) + ax.set_xticks(np.arange(0,70,10)) + ax.set_yticks(np.arange(-10,11,10)) + ax.set_xlim([0,50]) + ax.set_ylim(-10,10) + circle = plt.Circle((0,0),2.5,fill=True,color='r',edgecolor='k',alpha=0.5) + plt.gcf().gca().add_artist(circle) + plt.tight_layout() + plt.show() + + return Xg,Yg,Vxg,Vyg,Magn + + + +def plot_grid_cylinder_flow(Xg,Yg,Vxg,Vyg): + STEPx=1; STEPy=1 + # This creates the figure + fig, ax = plt.subplots(figsize=(6, 3)) + Magn=np.sqrt(Vxg**2+Vyg**2) + # Plot Contour + #CL=plt.contourf(Xg,Yg,Magn,levels=np.linspace(0,np.max(Magn),5)) + CL=plt.contourf(Xg,Yg,Magn,20,cmap='viridis',alpha=0.95) + # One possibility is to use quiver + STEPx=1; STEPy=1 + plt.quiver(Xg[::STEPx,::STEPy],Yg[::STEPx,::STEPy],\ + Vxg[::STEPx,::STEPy],Vyg[::STEPx,::STEPy],color='k') + plt.rc('text', usetex=True) + plt.rc('font', family='serif') + plt.rc('xtick',labelsize=12) + plt.rc('ytick',labelsize=12) + #fig.colorbar(CL,pad=0.05,fraction=0.025) + ax.set_aspect('equal') # Set equal aspect ratio + ax.set_xlabel('$x[mm]$',fontsize=13) + ax.set_ylabel('$y[mm]$',fontsize=13) + #ax.set_title('Tutorial 2: Cylinder Wake',fontsize=12) + ax.set_xticks(np.arange(0,70,10)) + ax.set_yticks(np.arange(-10,11,10)) + ax.set_xlim([0,50]) + ax.set_ylim(-10,10) + circle = plt.Circle((0,0),2.5,fill=True,color='r',edgecolor='k',alpha=0.5) + plt.gcf().gca().add_artist(circle) + plt.tight_layout() + plt.show() + + + diff --git a/modulo/utils/read_db.py b/modulo/utils/read_db.py new file mode 100644 index 0000000..45b0524 --- /dev/null +++ b/modulo/utils/read_db.py @@ -0,0 +1,263 @@ +import numpy as np +import os +from tqdm import tqdm +import math + + +class ReadData: + """ + A MODULO helper class for input data. ReadData allows to load the data directly before using MODULO, and + hence assemblying the data matrix D from data, if needed. + + """ + + + def __init__(self): + pass + + + @classmethod + def _from_dat(cls, folder, filename, N, N_S, N_T, + h: int = 0, f: int = 0, + c: int = 0, N_PARTITIONS: int =1,MR: bool = False): + """ + This method imports data (in the specified format) and then assemblies the corresponding + data matrix, D. + + :param folder: str + Folder in which the data is stored + :param filename: str + Name of the files to be imported + :param N: int + Components to be analysed + :param N_S: int + Number of points in space + :param N_T: int + Components to be analysed + :param h: int + Lines to be skipped from header + :param f: int + Lines to be skipped from footer + :param c: int + Columns to be skipped (for example if the first c columns + contain the mesh grid.) + :param N_PARTITIONS: int + Number of partitions. if =1 , then the matrix will be built + and returned in output. if >1, then it means that the memory + saving is active. The code will break the data into blocks + and stored in a local directory + + :return: np.array + Assembled DataMarix + + """ + path, dirs, files = next(os.walk(folder)) + #N_T = len(files) # this is dangerous because the folder could contain + # also files that are related to the mesh. We give N_T as input! + + + if N_PARTITIONS==1: # If you have only one partition (one matrix! ) + D = np.zeros((N_S, N_T)) + + print("\n \n Importing data with no partitions... \n \n") + + if MR: + print("Mean removal activated") + D_MEAN=np.zeros(N_S) + + + for k in tqdm(range(0, N_T)): + Name = folder + os.sep + filename % (k + 1) + '.dat' # Name of the file to read + # Read data from a file + DATA = np.genfromtxt(Name, #usecols=np.arange(0, 2), + skip_header=h, skip_footer=f) # Here we have the two colums + #Dat = DATA[1:, :] # Here we remove the first raw, containing the header + for ii in range(c, N + c): + tmp = DATA[:, ii] + if ii == c: + V = np.copy(tmp) + else: + V = np.concatenate([V, tmp], axis = 0) + if MR: + D_MEAN += 1/N_T*V # Snapshot contribution to the mean + + D[:, k] = V # Reshape and assign + + if MR: + print("Removing the mean from D ...") + D_Mr = D - D_MEAN.reshape(-1,1) # Mean Removed + print("Computing the mean-removed D ... ") + np.copyto(D, D_Mr) + del D_Mr + + + elif N_PARTITIONS>1: # then we enter in the memory saving loop + # prepare the folder to store the parittions + os.makedirs("MODULO_tmp/data_partitions/", exist_ok=True) + print("Memory Saving feature is active. Partitioning Data Matrix...") + + dim_col = math.floor(N_T / N_PARTITIONS) + columns_to_part = dim_col * N_PARTITIONS # These are integer multiples of N_PARTITIONS + vec=np.arange(0,columns_to_part) + # This gets the blocks + splitted_tmp = np.hsplit(vec, N_PARTITIONS) + if columns_to_part != N_T : + print("WARNING: the last "+str(N_T-1-splitted_tmp[N_PARTITIONS-1][-1])+' snapshots are not considered') + + if MR: + print("Mean removal activated") + D_MEAN=np.zeros(N_S) + + for ii in range(1,len(splitted_tmp)+1): + count=0 + print('Working on block '+str(ii)+'/'+str(N_PARTITIONS)) + D=np.zeros((N_S,len(splitted_tmp[0]))) + i1=splitted_tmp[ii-1][0]; i2=splitted_tmp[ii-1][-1] # ranges + for k in tqdm(range(i1,i2+1)): + Name = folder + os.sep + filename % (k + 1) + '.dat' # Name of the file to read + DATA = np.genfromtxt(Name, #usecols=np.arange(0, 2), + skip_header=h, skip_footer=f) # Here we have the two colums + for nn in range(c, N + c): + tmp = DATA[:, nn] + if nn == c: + V = np.copy(tmp) + else: + V = np.concatenate([V, tmp], axis = 0) + + if MR: + D_MEAN += 1/N_T*V # Snapshot contribution to the mean + + + D[:, count] = V # Reshape and assign + count+=1 + np.savez(f"MODULO_tmp/data_partitions/di_{ii}", di=D) + print('Partition '+str(ii)+'/'+str(N_PARTITIONS)+' saved') + + if MR: + print('Reloading the data for removing the mean') + for ii in range(1,len(splitted_tmp)+1): + print(f"Mean centering block {ii}") + di = np.load(f"MODULO_tmp/data_partitions/di_{ii}.npz")['di'] + di_mr=di - D_MEAN.reshape(-1,1) # Mean Removed + np.savez(f"MODULO_tmp/data_partitions/di_{ii}", di=di_mr) + else: + raise TypeError("number of partitions not valid.") + + return D + + + + + + @classmethod + def from_xls(cls, filename, **kwargs): + """ + This class method builds the df from an excel file. + :param filename: str + filename (with path if needed) to the df file. + :return: constructor for the class. + + """ + ## TBD + return + + @classmethod + def _from_csv(cls, folder, filename, N, N_S, + h: int = 0, f: int = 0, + c: int = 0): + """ + This method imports data (in the specified format) and then assemblies the corresponding + data matrix, D. + + :param folder: str + Folder in which the data is stored + :param filename: str + Name of the files to be imported + :param N number of components: int + Components to be analysed + :param h: int + Lines to be skipped from header + :param f: int + Lines to be skipped from footer + :param c: int + Columns to be skipped + + :return: np.array + Assembled DataMarix + + """ + path, dirs, files = next(os.walk(folder)) + files = [f for f in files if f.endswith('.csv')] + N_T = len(files) + D = np.zeros((N_S, N_T)) + + print("\n \n Importing data... \n \n") + + for k in tqdm(range(0, N_T)): + Name = folder + files[k] #os.sep + filename % (k + 1) + '.csv' # Name of the file to read + # Read data from a file + DATA = np.genfromtxt(Name, # usecols=np.arange(0, 2), + skip_header=h, skip_footer=f) # Here we have the two colums + # Dat = DATA[1:, :] # Here we remove the first raw, containing the header + for ii in range(c, N + c): + tmp = DATA[:, ii] + if ii == c: + V = np.copy(tmp) + else: + V = np.concatenate([V, tmp], axis=0) + + D[:, k] = V # Reshape and assign + + return D + + @classmethod + def _from_txt(cls, folder, filename, N, N_S, + h: int = 0, f: int = 0, + c: int = 0): + """ + This method imports data (in the specified format) and then assemblies the corresponding + data matrix, D. + + :param folder: str + Folder in which the data is stored + :param filename: str + Name of the files to be imported + :param N number of components: int + Components to be analysed + :param h: int + Lines to be skipped from header + :param f: int + Lines to be skipped from footer + :param c: int + Columns to be skipped + + :return: np.array + Assembled DataMarix + + """ + path, dirs, files = next(os.walk(folder)) + N_T = len(files) + D = np.zeros((N_S, N_T)) + + print("\n \n Importing data... \n \n") + + for k in tqdm(range(0, N_T)): + Name = folder + os.sep + filename % (k + 1) + '.txt' # Name of the file to read + # Read data from a file + DATA = np.genfromtxt(Name, # usecols=np.arange(0, 2), + skip_header=h, skip_footer=f) # Here we have the two colums + # Dat = DATA[1:, :] # Here we remove the first raw, containing the header + for ii in range(c, N + c): + tmp = DATA[:, ii] + if ii == c: + V = np.copy(tmp) + else: + V = np.concatenate([V, tmp], axis=0) + + D[:, k] = V # Reshape and assign + + return D + + + +