Added post-processing scripts

Author: brunosorban

docs/notebooks/monte_carlo_analysis/mc_example.py

@@ -187,5 +187,5 @@
 )
 
 test_dispersion.simulate(
-    number_of_simulations=1000, append=False, light_mode=True, parallel=True
+    number_of_simulations=10, append=False, light_mode=False, parallel=True
 )

rocketpy/simulation/monte_carlo.py

@@ -938,9 +938,16 @@ def dict_to_h5(h5_file, path, dic):
         """
         for key, item in dic.items():
             if isinstance(
-                item, (np.ndarray, np.int64, np.float64, str, bytes, int, float)
+                item, (np.int64, np.float64, int, float)
             ):
-                h5_file[path + key] = item
+                data = np.array([[item]])
+                h5_file.create_dataset(path + key, data=data, shape=data.shape, dtype=data.dtype)
+            elif isinstance(item, np.ndarray):
+                if len(item.shape) < 2:
+                    item = item.reshape(-1, 1) # Ensure it is a column vector
+                h5_file.create_dataset(path + key, data=item, shape=item.shape, dtype=item.dtype)
+            elif isinstance(item, (str, bytes)):
+                h5_file.create_dataset(path + key, data=item)
             elif isinstance(item, Function):
                 raise TypeError(
                     "Function objects should be preprocessed before saving."

rocketpy/simulation/sim_config/serializer.py

@@ -25,7 +25,7 @@ def function_serializer(function_object: Function, t_range=None):
         else:
             raise ValueError("t_range must be specified for callable functions")
 
-    return func.source
+    return func.get_source()
 
     # with open(csv_path, "w+", newline="") as f:
     #     writer = csv.writer(f)

rocketpy/stochastic/post_processing/__init__.py

@@ -0,0 +1,52 @@
+from pathlib import Path
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+from rocketpy.stochastic.post_processing.stochastic_cache import \
+    SimulationCache
+
+
+def compute_apogee(file_name, batch_path):
+    cache = SimulationCache(
+        file_name,
+        batch_path,
+    )
+    apogee = cache.read_outputs('apogee')
+
+    mean_apogee = float(np.nanmean(apogee, axis=0))
+
+    return apogee, mean_apogee
+
+
+def plot_apogee(batch_path, apogee, mean_apogee, save=False, show=True):
+    # Histogram
+    fig, ax = plt.subplots()
+    ax.hist(apogee.flatten(), bins=50)
+    ax.axvline(
+        mean_apogee,
+        color='black',
+        linestyle='--',
+        linewidth=2,
+        label=f'Mean: {mean_apogee:.2f} m',
+    )
+    ax.set_title('Mean Apogee Distribution')
+    ax.set_xlabel('Apogee [m]')
+    ax.set_ylabel('Frequency')
+    ax.legend()
+
+    if save:
+        plt.savefig(batch_path / 'mean_apogee_distribution.png')
+    if show:
+        plt.show()
+
+
+def run(file_name, batch_path, save, show):
+    apogee, mean_apogee = compute_apogee(file_name, batch_path)
+    plot_apogee(batch_path, apogee, mean_apogee, save=save, show=show)
+
+
+if __name__ == '__main__':
+    batch_path = Path("mc_simulations/")
+    file_name = 'monte_carlo_class_example'
+    run(file_name, batch_path, save=True, show=True)

rocketpy/stochastic/post_processing/scripts/__init__.py

@@ -0,0 +1,87 @@
+from pathlib import Path
+
+import matplotlib.pyplot as plt
+import numpy as np
+from matplotlib.patches import Ellipse
+
+from rocketpy.stochastic.post_processing.stochastic_cache import \
+    SimulationCache
+
+# 1-3 sigma
+lower_percentiles = [0.16, 0.03, 0.003]
+upper_percentiles = [0.84, 0.97, 0.997]
+
+
+# Define function to calculate eigen values
+def eigsorted(cov):
+    vals, vecs = np.linalg.eigh(cov)
+    order = vals.argsort()[::-1]
+    return vals[order], vecs[:, order]
+
+
+def compute_impact(file_name, batch_path, save, show):
+    cache = SimulationCache(
+        file_name,
+        batch_path,
+    )
+    x_impact = cache.read_outputs('x_impact')
+    y_impact = cache.read_outputs('y_impact')
+
+    x_mean_impact = np.nanmean(x_impact, axis=0)
+    y_mean_impact = np.nanmean(y_impact, axis=0)
+
+    # Calculate error ellipses for impact
+    impact_cov = np.cov(x_impact.flatten(), y_impact.flatten())
+    impact_vals, impactVecs = eigsorted(impact_cov)
+    impact_theta = np.degrees(np.arctan2(*impactVecs[:, 0][::-1]))
+    impact_w, impactH = 2 * np.sqrt(impact_vals)
+
+    fig, ax = plt.subplots()
+    ax.scatter(x_impact, y_impact, c='blue')
+    ax.scatter(
+        x_mean_impact,
+        y_mean_impact,
+        marker='x',
+        c='red',
+        label='Mean Impact Point',
+    )
+
+    # Draw error ellipses for impact
+    impact_ellipses = []
+    for j in [1, 2, 3]:
+        impactEll = Ellipse(
+            xy=(np.mean(x_impact), np.mean(y_impact)),
+            width=impact_w * j,
+            height=impactH * j,
+            angle=impact_theta,
+            color="black",
+        )
+        impactEll.set_facecolor((0, 0, 1, 0.2))
+        impact_ellipses.append(impactEll)
+        ax.add_artist(impactEll)
+
+    ax.set_xlabel('X Impact Point [m]')
+    ax.set_ylabel('Y Impact Point [m]')
+    ax.set_title('Impact Point Distribution')
+    ax.set_aspect('equal')
+    ax.legend()
+    ax.grid()
+
+    if save:
+        plt.savefig(batch_path / 'mean_impact_distribution.png')
+
+    if show:
+        plt.show()
+    plt.show()
+
+
+def run(file_name, batch_path, save, show):
+    compute_impact(file_name, batch_path, save, show)
+
+
+if __name__ == '__main__':
+    # import easygui
+
+    batch_path = Path("mc_simulations/")
+    file_name = 'monte_carlo_class_example'
+    run(file_name, batch_path, save=True, show=True)

rocketpy/stochastic/post_processing/scripts/stochastic_apogee.py

@@ -0,0 +1,42 @@
+from pathlib import Path
+
+import matplotlib.pyplot as plt
+
+from rocketpy.stochastic.post_processing.stochastic_cache import \
+    SimulationCache
+
+
+def compute_mach(file_name, batch_path, save, show):
+    cache = SimulationCache(
+        file_name,
+        batch_path,
+    )
+    mach = cache.read_outputs('mach_number')
+
+    fig, ax = plt.subplots()
+    for i in range(len(mach)):
+        ax.plot(mach[i, :, 0], mach[i, :, 1], c='blue')
+    ax.set_xlabel('Time [s]')
+    ax.set_ylabel('Mach [-]')
+    ax.set_title('Mach number Distribution')
+    ax.legend()
+    ax.grid()
+
+    if save:
+        plt.savefig(batch_path / 'mach_distribution.png')
+
+    if show:
+        plt.show()
+    plt.show()
+
+
+def run(file_name, batch_path, save, show):
+    compute_mach(file_name, batch_path, save, show)
+
+
+if __name__ == '__main__':
+    # import easygui
+
+    batch_path = Path("mc_simulations/")
+    file_name = 'monte_carlo_class_example'
+    run(file_name, batch_path, save=True, show=True)

rocketpy/stochastic/post_processing/scripts/stochastic_impact_point.py

@@ -0,0 +1,46 @@
+from pathlib import Path
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+from rocketpy.stochastic.post_processing.stochastic_cache import \
+    SimulationCache
+
+
+def compute_maxQ(file_name, batch_path, save, show):
+    cache = SimulationCache(
+        file_name,
+        batch_path,
+    )
+    dyn_press = cache.read_outputs('dynamic_pressure') / 1000
+    maxQarg = np.nanargmax(dyn_press[:, :, 1], axis=1)
+    maxQ = dyn_press[np.arange(len(dyn_press)), maxQarg]
+
+    fig, ax = plt.subplots()
+    for i in range(len(dyn_press)):
+        ax.plot(dyn_press[i, :, 0], dyn_press[i, :, 1], c='blue')
+    ax.scatter(maxQ[:, 0], maxQ[:, 1], c='red', label='Max Q')
+    ax.set_xlabel('Time [s]')
+    ax.set_ylabel('Dynamic Pressure [kPa]')
+    ax.set_title('Max Q Distribution')
+    ax.legend()
+    ax.grid()
+
+    if save:
+        plt.savefig(batch_path / 'maxQ_distribution.png')
+
+    if show:
+        plt.show()
+    plt.show()
+
+
+def run(file_name, batch_path, save, show):
+    compute_maxQ(file_name, batch_path, save, show)
+
+
+if __name__ == '__main__':
+    # import easygui
+
+    batch_path = Path("mc_simulations/")
+    file_name = 'monte_carlo_class_example'
+    run(file_name, batch_path, save=True, show=True)

rocketpy/stochastic/post_processing/scripts/stochastic_mach.py

@@ -0,0 +1,60 @@
+from functools import lru_cache
+from pathlib import Path
+
+import h5py
+import numpy as np
+
+
+class SimulationCache:
+    def __init__(self, file_name, batch_path):
+        self.file_name = file_name
+        self.batch_path = Path(batch_path)
+        self.sim_number = self._get_sim_number()
+
+    def _get_sim_number(self):
+        with h5py.File(self.batch_path / (self.file_name + ".inputs.h5"), 'r') as f:
+            return len(f.keys())
+
+    @lru_cache(maxsize=32)
+    def read_inputs(self, var_name):
+        data = []
+
+        with h5py.File(self.batch_path / (self.file_name + ".inputs.h5"), 'r') as f:
+            for i in range(self.sim_number):
+                data.append(np.array(f[str(i)][var_name]))
+
+        # Add nans at the end of the data and convert to numpy array
+        max_len = max([len(d) for d in data])
+        data = [np.concatenate([d, np.full(max_len - len(d), np.nan)]) for d in data]
+
+        return np.array(data)
+
+    @lru_cache(maxsize=32)
+    def read_outputs(self, var_name):
+        data = []
+
+        with h5py.File(self.batch_path / (self.file_name + ".outputs.h5"), 'r') as f:
+            for i in range(self.sim_number):
+                data.append(np.array(f[str(i)][var_name]))
+
+        # Add nans at the end of the data and convert to numpy array
+        max_len = max([len(d) for d in data])
+        data = [
+            np.concatenate([d, np.full((max_len - len(d), d.shape[1]), np.nan)])
+            for d in data
+        ]
+
+        return np.array(data)
+
+
+if __name__ == '__main__':
+    import easygui
+
+    batch_path = easygui.diropenbox()
+    cache = SimulationCache(
+        'monte_carlo_class_example',
+        batch_path,
+    )
+    data = cache.read_outputs('ax')
+    print(data)
+    print(data.shape)