Problems with predictions

05b_EvaluateGismo_GP.py

@@ -78,7 +78,10 @@ for model_name in models:
                 stdout=subprocess.PIPE,
                 text=True
             ).communicate()[0]
-            
+
+            print(f'sample: {sample}')
+            print(f'program_output: {program_output}')
+
             # Catch the errors
             matches = re.findall(pattern, program_output)
 

05c_EvaluateGismo_RBF.py

+import matplotlib.pyplot as plt
+import numpy as np
+import pandas as pd
+import subprocess
+import re
+import os
+import shutil
+
+from models.RBFInterpolator import RadialBasisRegressionModel
+
+from miscellaneous.DataPreparation import reconstruct_array, save_numpy_to_xml
+
+VELOCITY_SIZE = 1568
+PRESSURE_SIZE = 645
+
+path_to_executable = 'gismo/stokes_recreation'
+models = {
+    # 'radialBasisMod_Linear': RadialBasisRegressionModel(),
+    'radialBasisMod_thinplatespline': RadialBasisRegressionModel(),
+    # 'radialBasisMod_cubic': RadialBasisRegressionModel(),
+    # 'radialBasisMod_quintic': RadialBasisRegressionModel(),
+    # 'radialBasisMod_multiquadric': RadialBasisRegressionModel(),
+    # 'radialBasisMod_inversemultiquadric': RadialBasisRegressionModel(),
+    # 'radialBasisMod_inversequadratic': RadialBasisRegressionModel(),
+    # 'radialBasisMod_gaussian': RadialBasisRegressionModel(),
+}
+
+df = pd.read_excel(f'Data/test/parameter_input.xlsx').transpose()
+SAMPLES_TEST = df.shape[1]
+pattern = r'(L1|L2|Linf|H1)\s*\(?(rel\.)?\)?:\s*([\d.eE+-]+|inf)'
+
+for model_name in models:
+    print(f'Model: \t\t\t{model_name}')
+    # Create temporary folder for results
+    if not os.path.exists(f'temp_{model_name}'):
+        os.mkdir(f'temp_{model_name}')
+    
+    # Read POD modes
+    data_training = np.load(f'models/trained_models/{model_name}/MetaData.npz')
+    R = data_training['R']
+
+    # Create error vector
+    L1_p_vec, L2_p_vec, Linf_p_vec, H1_p_vec, L1_rel_p_vec, L2_rel_p_vec, Linf_rel_p_vec = [], [], [], [], [], [], []
+    L1_vx_vec, L2_vx_vec, Linf_vx_vec, H1_vx_vec, L1_rel_vx_vec, L2_rel_vx_vec, Linf_rel_vx_vec = [], [], [], [], [], [], []
+    L1_vy_vec, L2_vy_vec, Linf_vy_vec, H1_vy_vec, L1_rel_vy_vec, L2_rel_vy_vec, Linf_rel_vy_vec = [], [], [], [], [], [], []
+    
+    for R_ in R:
+        print(f'Singular value: \t{R_}')
+        # Load model
+        model = models[model_name]
+        model_velocity, model_pressure = model.load_model(f'{model_name}/velocity_{int(R_)}.pkl', f'{model_name}/pressure_{int(R_)}.pkl')
+
+        # Predict values
+        # Predicting test values
+        test_velocity_predict, test_pressure_predict = model.predict_test(model_velocity, model_pressure, rescale=True)
+
+        velocity_matrix_original_empty = np.empty((VELOCITY_SIZE, SAMPLES_TEST), np.float32)
+        pressure_matrix_original_empty = np.empty((PRESSURE_SIZE, SAMPLES_TEST), np.float32)
+
+        L1_p_vec_sample, L2_p_vec_sample, Linf_p_vec_sample, H1_p_vec_sample, L1_rel_p_vec_sample, L2_rel_p_vec_sample, Linf_rel_p_vec_sample = [], [], [], [], [], [], []
+        L1_vx_vec_sample, L2_vx_vec_sample, Linf_vx_vec_sample, H1_vx_vec_sample, L1_rel_vx_vec_sample, L2_rel_vx_vec_sample, Linf_rel_vx_vec_sample = [], [], [], [], [], [], []
+        L1_vy_vec_sample, L2_vy_vec_sample, Linf_vy_vec_sample, H1_vy_vec_sample, L1_rel_vy_vec_sample, L2_rel_vy_vec_sample, Linf_rel_vy_vec_sample = [], [], [], [], [], [], []
+
+        # Iterate over each testcase
+        for sample in range(SAMPLES_TEST):
+        # for sample in range(10):
+            save_numpy_to_xml(test_velocity_predict[:,sample], f'temp_{model_name}/predict_velocity_field_{sample}.xml', rows=VELOCITY_SIZE, cols='1')
+
+            save_numpy_to_xml(test_pressure_predict[:,sample], f'temp_{model_name}/predict_pressure_field_{sample}.xml', rows=PRESSURE_SIZE, cols='1')
+
+            # Get all the input files for the evaluation
+            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'
+
+            # Run evaluation with gismo script
+            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]
+
+            print(f'sample: {sample}')
+            print(f'program_output: {program_output}')
+
+            # Catch the errors
+            matches = re.findall(pattern, program_output)
+
+            results = []
+            # Parse the matches
+            # ToDo: Check if still necessary with new script?
+            for metric, rel, value in matches:
+                key = f"{metric} (rel.)" if rel else metric
+                if type(value) == str and value.lower() == 'inf':
+                    results.append(1e+10)
+                else:
+                    results.append(float(value))
+
+            # Collect the errors
+            # Pressure
+            L1_p_vec_sample.append(results[0])
+            L2_p_vec_sample.append(results[1])
+            Linf_p_vec_sample.append(results[2])
+            H1_p_vec_sample.append(results[3])
+            L1_rel_p_vec_sample.append(results[4])
+            L2_rel_p_vec_sample.append(results[5])
+            Linf_rel_p_vec_sample.append(results[6])
+
+            # Velocity (x)
+            # L1_vx_vec_sample.append(results[7])
+            # L2_vx_vec_sample.append(results[8])
+            # Linf_vx_vec_sample.append(results[9])
+            # H1_vx_vec_sample.append(results[10])
+            # L1_rel_vx_vec_sample.append(results[11])
+            # L2_rel_vx_vec_sample.append(results[12])
+            # Linf_rel_vx_vec_sample.append(results[13])
+
+            # Velocity (y)
+            # L1_vy_vec_sample.append(results[14])
+            # L2_vy_vec_sample.append(results[15])
+            # Linf_vy_vec_sample.append(results[16])
+            # H1_vy_vec_sample.append(results[17])
+            # L1_rel_vy_vec_sample.append(results[18])
+            # L2_rel_vy_vec_sample.append(results[19])
+            # Linf_rel_vy_vec_sample.append(results[20])
+
+        L1_p_vec.append(np.mean(L1_p_vec_sample))
+        L2_p_vec.append(np.mean(L2_p_vec_sample))
+        Linf_p_vec.append(np.mean(Linf_p_vec_sample))
+        H1_p_vec.append(np.mean(H1_p_vec_sample))
+        L1_rel_p_vec.append(np.mean(L1_rel_p_vec_sample))
+        L2_rel_p_vec.append(np.mean(L2_rel_p_vec_sample))
+        Linf_rel_p_vec.append(np.mean(Linf_rel_p_vec_sample))
+
+        # L1_vx_vec.append(np.mean(L1_vx_vec_sample))
+        # L2_vx_vec.append(np.mean(L2_vx_vec_sample))
+        # Linf_vx_vec.append(np.mean(Linf_vx_vec_sample))
+        # H1_vx_vec.append(np.mean(H1_vx_vec_sample))
+        # L1_rel_vx_vec.append(np.mean(L1_rel_vx_vec_sample))
+        # L2_rel_vx_vec.append(np.mean(L2_rel_vx_vec_sample))
+        # Linf_rel_vx_vec.append(np.mean(Linf_rel_vx_vec_sample))
+
+        # L1_vy_vec.append(np.mean(L1_vy_vec_sample))
+        # L2_vy_vec.append(np.mean(L2_vy_vec_sample))
+        # Linf_vy_vec.append(np.mean(Linf_vy_vec_sample))
+        # H1_vy_vec.append(np.mean(H1_vy_vec_sample))
+        # L1_rel_vy_vec.append(np.mean(L1_rel_vy_vec_sample))
+        # L2_rel_vy_vec.append(np.mean(L2_rel_vy_vec_sample))
+        # Linf_rel_vy_vec.append(np.mean(Linf_rel_vy_vec_sample))
+
+
+
+    # Save the errors
+    np.savez(f'models/trained_models/{model_name}/errors_recreated.npz', L1_p_vec=L1_p_vec, L2_p_vec=L2_p_vec, Linf_p_vec=Linf_p_vec, H1_p_vec=H1_p_vec, L1_rel_p_vec=L1_rel_p_vec, L2_rel_p_vec=L2_rel_p_vec, Linf_rel_p_vec=Linf_rel_p_vec, L1_vx_vec=L1_vx_vec, L2_vx_vec=L2_vx_vec, Linf_vx_vec=Linf_vx_vec, H1_vx_vec=H1_vx_vec, L1_rel_vx_vec=L1_rel_vx_vec, L2_rel_vx_vec=L2_rel_vx_vec, Linf_rel_vx_vec=Linf_rel_vx_vec, L1_vy_vec=L1_vy_vec, L2_vy_vec=L2_vy_vec, Linf_vy_vec=Linf_vy_vec, H1_vy_vec=H1_vy_vec, L1_rel_vy_vec=L1_rel_vy_vec, L2_rel_vy_vec=L2_rel_vy_vec, Linf_rel_vy_vec=Linf_rel_vy_vec, R=R)
+
+    # Remove temporary folder
+    if os.path.exists(f'temp_{model_name}'):
+        shutil.rmtree(f'temp_{model_name}')
\ No newline at end of file

07_VisualizeErrorsRecreated.py

@@ -22,7 +22,7 @@ models = [
     # 'gaussianProcessMod_Matern',
     # 'gaussianProcessMod_RationalQuadratic',
     # 'gaussianProcessMod_ExpSineSquared',
-    #     'radialBasisMod_Linear',
+    'radialBasisMod_Linear',
     # 'radialBasisMod_thinplatespline',
     # 'radialBasisMod_cubic',
     # 'radialBasisMod_quintic',
@@ -39,7 +39,7 @@ label_names = [
     # 'Gaussian Process Matern',
     # 'Gaussian Process RationalQuadratic',
     # 'Gaussian Process ExpSineSquared',
-    # 'Radial Basis Linear',
+    'Radial Basis Linear',
     # 'Radial Basis ThinPlateSpline',
     # 'Radial Basis Cubic',
     # 'Radial Basis Quintic',
@@ -52,16 +52,22 @@ label_names = [
 fig_pressure, axs_pressure = plt.subplots(1, 1, figsize=(3.54*4, 3.54*2), layout='constrained')
 fig_velx, axs_velx = plt.subplots(1, 1, figsize=(3.54*4, 3.54*2), layout='constrained')
 fig_vely, axs_vely = plt.subplots(1, 1, figsize=(3.54*4, 3.54*2), layout='constrained')
+L2_rel_p_vec = []
+L2_rel_vx_vec = []
+L2_rel_vy_vec = []
 for model in models:
     data = np.load(f'models/trained_models/{model}/errors_recreated.npz')
     R = data['R']
-    L2_rel_p_vec = data['L2_rel_p_vec']
-    L2_rel_vx_vec = data['L2_rel_vx_vec']
-    L2_rel_vy_vec = data['L2_rel_vy_vec']
+    L2_rel_p = data['L2_rel_p_vec']
+    L2_rel_p_vec.append(L2_rel_p)
+    print(f'model: {model}')
+    print(f'L2_rel_p_vec.shape: {L2_rel_p_vec.shape}')
+    # L2_rel_vx_vec = data['L2_rel_vx_vec']
+    # L2_rel_vy_vec = data['L2_rel_vy_vec']
 
     axs_pressure.bar(label_names, L2_rel_p_vec)
-    axs_velx.bar(label_names, L2_rel_vx_vec)
-    axs_vely.bar(label_names, L2_rel_vy_vec)
+    # axs_velx.bar(label_names, L2_rel_vx_vec)
+    # axs_vely.bar(label_names, L2_rel_vy_vec)
     
 for axs in [axs_pressure, axs_velx, axs_vely]:
     axs.set_ylabel('L2_rel')