Evalute Error over Singular-Values for gismo recreate

08_EvaluteReconstructedGismoSVD.py

+import matplotlib.pyplot as plt
+import numpy as np
+import pandas as pd
+import os
+import shutil
+import subprocess 
+import re
+
+# from models.RegressionModels import RegressionModels
+from models.LinearRegression import LinearRegressionModel
+from models.GaussianProcessRegressor import GaussianProcessRegressionModel
+from models.RBFInterpolator import RadialBasisRegressionModel
+
+from miscellaneous.DataPreparation import reconstruct_array, save_numpy_to_xml
+
+path_to_executable = 'stokes_recreation/stokes_recreation'
+paraview_export = False
+df = pd.read_excel(f'Data/test/parameter_input.xlsx').transpose()
+N_FILES_STAGE = df.shape[1]
+pattern = r'L1\s*:\s*([\d.]+)|L2\s*:\s*([\d.]+)|Linf\s*:\s*([\d.]+)|H1\s*:\s*([\d.]+)|L1\s+\(rel\.\):\s*([\d.]+)|L2\s+\(rel\.\):\s*([\d.]+)|Linf+\(rel\.\):\s*([\d.]+)'
+
+VELOCITY_SIZE = 3840
+PRESSURE_SIZE = 1080
+
+models = {
+    'linregMod': LinearRegressionModel(),
+    'gaussianProcessMod': GaussianProcessRegressionModel(),
+    'radialBasisMod': RadialBasisRegressionModel(),
+}
+
+SINGULAR_VALUES = np.linspace(1,21,21)
+
+# Iterate over the models
+for model_name in models:
+    print(f'Model: \t\t\t{model_name}')
+    if not os.path.exists(f'temp_{model_name}'):
+        os.mkdir(f'temp_{model_name}')
+
+    # iterate over the singular values
+    L1_svd_p, L2_svd_p, Linf_svd_p, H1_svd_p, L1_rel_svd_p, L2_rel_svd_p, Linf_rel_svd_p = [], [], [], [], [], [], []
+    L1_svd_v, L2_svd_v, Linf_svd_v, H1_svd_v, L1_rel_svd_v, L2_rel_svd_v, Linf_rel_svd_v = [], [], [], [], [], [], []
+    for singular_value in SINGULAR_VALUES:
+        print(f'Singular value: \t{singular_value}')
+        # Load model
+        model = models[model_name]
+        model_velocity, model_pressure = model.load_model(f'{model_name}/velocity_{int(singular_value)}.pkl', f'{model_name}/pressure_{int(singular_value)}.pkl')
+
+        # Predicting test values
+        test_velocity_predict, test_pressure_predict = model.predict_test(model_velocity, model_pressure, rescale=True)
+
+        # Reconstruct the arrays
+        # Get dimensions
+        samples = test_velocity_predict.shape[1]
+        velocity_matrix_original_empty = np.empty((VELOCITY_SIZE, samples), np.float32)
+        pressure_matrix_original_empty = np.empty((PRESSURE_SIZE, samples), np.float32)
+        velocity_removed_rows_identifier = np.loadtxt('Data/train/matrices/velocity_removed_rows_identifier.csv')
+        pressure_removed_rows_identifier = np.loadtxt('Data/train/matrices/pressure_removed_rows_identifier.csv')
+
+        # Convert the removed_rows_identifier arrays to 1D arrays of integers
+        velocity_removed_rows_identifier = velocity_removed_rows_identifier.astype(int)
+        pressure_removed_rows_identifier = pressure_removed_rows_identifier.astype(int)
+
+        # Reconstructing the arrays
+        test_velocity_predict_reconstructed = reconstruct_array(test_velocity_predict, velocity_removed_rows_identifier, velocity_matrix_original_empty.shape)
+        test_pressure_predict_reconstructed = reconstruct_array(test_pressure_predict, pressure_removed_rows_identifier, pressure_matrix_original_empty.shape)
+
+        # Save each sample of the reconstructed arrays separately
+        L1_singularvalue_p, L2_singularvalue_p, Linf_singularvalue_p, H1_singularvalue_p, L1_rel_singularvalue_p, L2_rel_singularvalue_p, Linf_rel_singularvalue_p = [], [], [], [], [], [], []
+        L1_singularvalue_v, L2_singularvalue_v, Linf_singularvalue_v, H1_singularvalue_v, L1_rel_singularvalue_v, L2_rel_singularvalue_v, Linf_rel_singularvalue_v = [], [], [], [], [], [], []
+
+        for sample in range(samples):
+        # for sample in range(10):
+            save_numpy_to_xml(test_velocity_predict_reconstructed[:,sample], f'temp_{model_name}/predict_velocity_field_{sample}.xml', rows='3840', cols='1')
+
+            save_numpy_to_xml(test_pressure_predict_reconstructed[:,sample], f'temp_{model_name}/predict_pressure_field_{sample}.xml', rows='1080', cols='1')
+
+            geometry_file = f'Data/test/microstructure/index_{sample}.xml'
+            velocity_solution = f'Data/test/velocities/velocity_field_{sample}.xml'
+            pressure_solution = f'Data/test/pressure/pressure_field_{sample}.xml'
+            velocity_predicted = f'temp_{model_name}/predict_velocity_field_{sample}.xml'
+            pressure_predicted = f'temp_{model_name}/predict_pressure_field_{sample}.xml'
+
+            program_output = subprocess.Popen(
+            [f'./{path_to_executable}', f'-f {geometry_file}', f'-v {velocity_solution}', f'-p {pressure_solution}', f'-w {velocity_predicted}', f'-q {pressure_predicted}', '-e 1', '--no-plot'],
+            stdout=subprocess.PIPE,
+            text=True
+            ).communicate()[0]
+
+            matches = re.findall(pattern, program_output)
+
+            # Extract error values into variables
+            L1_p_, L2_p_, Linf_p_, H1_p_, L1_rel_p_, L2_rel_p_, Linf_rel_p_, L1_v_, L2_v_, Linf_v_, H1_v_, L1_rel_v_, L2_rel_v_, Linf_rel_v_ = (float(value) for match in matches for value in match if value)
+
+            L1_singularvalue_p.append(L1_p_)
+            L2_singularvalue_p.append(L2_p_)
+            Linf_singularvalue_p.append(Linf_p_)
+            H1_singularvalue_p.append(H1_p_)
+            L1_rel_singularvalue_p.append(L1_rel_p_)
+            L2_rel_singularvalue_p.append(L2_rel_p_)
+            Linf_rel_singularvalue_p.append(Linf_rel_p_)
+            L1_singularvalue_v.append(L1_v_)
+            L2_singularvalue_v.append(L2_v_)
+            Linf_singularvalue_v.append(Linf_v_)
+            H1_singularvalue_v.append(H1_v_)
+            L1_rel_singularvalue_v.append(L1_rel_v_)
+            L2_rel_singularvalue_v.append(L2_rel_v_)
+            Linf_rel_singularvalue_v.append(Linf_rel_v_)
+        for arr in [L1_singularvalue_p, L2_singularvalue_p, Linf_singularvalue_p, H1_singularvalue_p, L1_rel_singularvalue_p, L2_rel_singularvalue_p, Linf_rel_singularvalue_p, L1_singularvalue_v, L2_singularvalue_v, Linf_singularvalue_v, H1_singularvalue_v, L1_rel_singularvalue_v, L2_rel_singularvalue_v, Linf_rel_singularvalue_v]:
+            arr = np.array(arr)
+
+        L1_svd_p.append(np.mean(L1_singularvalue_p))
+        L2_svd_p.append(np.mean(L2_singularvalue_p))
+        Linf_svd_p.append(np.mean(Linf_singularvalue_p))
+        H1_svd_p.append(np.mean(H1_singularvalue_p))
+        L1_rel_svd_p.append(np.mean(L1_rel_singularvalue_p))
+        L2_rel_svd_p.append(np.mean(L2_rel_singularvalue_p))
+        Linf_rel_svd_p.append(np.mean(Linf_rel_singularvalue_p))
+        L1_svd_v.append(np.mean(L1_singularvalue_v))
+        L2_svd_v.append(np.mean(L2_singularvalue_v))
+        Linf_svd_v.append(np.mean(Linf_singularvalue_v))
+        H1_svd_v.append(np.mean(H1_singularvalue_v))
+        L1_rel_svd_v.append(np.mean(L1_rel_singularvalue_v))
+        L2_rel_svd_v.append(np.mean(L2_rel_singularvalue_v))
+        Linf_rel_svd_v.append(np.mean(Linf_rel_singularvalue_v))
+
+    # Save the errors
+    Errors = {
+        'L1_svd_p': L1_svd_p,
+        'L2_svd_p': L2_svd_p,
+        'Linf_svd_p': Linf_svd_p,
+        'H1_svd_p': H1_svd_p,
+        'L1_rel_svd_p': L1_rel_svd_p,
+        'L2_rel_svd_p': L2_rel_svd_p,
+        'Linf_rel_svd_p': Linf_rel_svd_p,
+        'L1_svd_v': L1_svd_v,
+        'L2_svd_v': L2_svd_v,
+        'Linf_svd_v': Linf_svd_v,
+        'H1_svd_v': H1_svd_v,
+        'L1_rel_svd_v': L1_rel_svd_v,
+        'L2_rel_svd_v': L2_rel_svd_v,
+        'Linf_rel_svd_v': Linf_rel_svd_v
+    }
+    df = pd.DataFrame(Errors)
+    df['SingularValue'] = SINGULAR_VALUES
+    df.to_excel(f'Results/{model_name}/ErrorsSVD.xlsx')
+    # os.rmdir(f'temp_{model_name}')
+    shutil.rmtree(f'temp_{model_name}')
+
+    # Visualize errors over singular values
+    norm_names = ['$L_1$ (rel)', '$L_2$ (rel)', '$L_\infty$ (rel)']
+    norm_velocities = [L1_rel_svd_v, L2_rel_svd_v, Linf_rel_svd_v]
+    norm_pressures = [L1_rel_svd_p, L2_rel_svd_p, Linf_rel_svd_p]
+    fig, axs = plt.subplots(nrows=3, ncols=2, figsize=(3.54*3,3.54*3))
+    fig.suptitle(f'Pressure - {model_name} - Errors recreated with Gismo')
+    
+    for i, norm_name in enumerate(norm_names):
+        axs[i,0].semilogy(SINGULAR_VALUES, norm_pressures[i])
+        axs[i,0].set_title(f'{norm_name} pressure')
+        axs[i,0].set_xlabel('Singular value')
+        axs[i,0].set_ylabel(f'{norm_name}')
+        axs[i,1].semilogy(SINGULAR_VALUES, norm_velocities[i])
+        axs[i,1].set_title(f'{norm_name} velocity')
+        axs[i,1].set_xlabel('Singular value')
+        axs[i,1].set_ylabel(f'{norm_name}')
+
+    fig.tight_layout()
+    fig.savefig(f'Results/{model_name}/recreated_errors_svd.jpg')