updated sensitivity analysis

Author: matthiaskoenig

src/sbmlsim/sensitivity/analysis.py

@@ -105,21 +105,26 @@ class SensitivityAnalysis:
 
     def __init__(self,
                  sensitivity_simulation: SensitivitySimulation,
-                 parameters: list[SensitivityParameter]
+                 parameters: list[SensitivityParameter],
+                 results_path: Path,
                  ) -> None:
         """Create a sensitivity analysis for given parameter ids.
 
         Based on the results matrix the sensitivity is calculated.
         """
         self.sensitivity_simulation = sensitivity_simulation
 
+        # outputs to calculate sensitivity on; shape: (num_outputs,)
+        self.outputs: list[SensitivityOutput] = sensitivity_simulation.outputs
+        self.output_ids: list[str] = [q.uid for q in self.outputs]
+
         # parameters to vary; shape: (num_parameters,)
         self.parameters: list[SensitivityParameter] = parameters
         self.parameter_ids: list[str] = [p.uid for p in self.parameters]
 
-        # outputs to calculate sensitivity on; shape: (num_outputs,)
-        self.outputs: list[SensitivityOutput] = sensitivity_simulation.outputs
-        self.output_ids: list[str] = [q.uid for q in self.outputs]
+        # storage directory
+        self.results_path: Path = results_path
+        results_path.mkdir(parents=True, exist_ok=True)
 
         # parameter samples for sensitivity; shape: (num_samples x num_parameters)
         self.samples: Optional[xr.DataArray] = None
@@ -255,17 +260,19 @@ def plot_sensitivity(
         cluster_rows: bool = True,
         title: Optional[str] = None,
         cmap: str = "seismic",
+        fig_path: Optional[Path] = None,
         **kwargs
     ) -> None:
         df = self.sensitivity_df(key=key)
         heatmap(
             df=df,
-            parameter_labels={p.uid: p.name for p in self.parameters},
+            parameter_labels={p.uid: f"{p.uid}: {p.name}" for p in self.parameters},
             output_labels={q.uid: q.name for q in self.outputs},
             cutoff=cutoff,
             cluster_rows=cluster_rows,
             title=title,
             cmap=cmap,
+            fig_path=fig_path,
             **kwargs
         )
 
@@ -296,9 +303,11 @@ class LocalSensitivityAnalysis(SensitivityAnalysis):
     """
 
     def __init__(self, sensitivity_simulation: SensitivitySimulation,
-                 parameters: list[SensitivityParameter], difference: float = 0.01):
+                 parameters: list[SensitivityParameter],
+                 results_path: Path,
+                 difference: float = 0.01):
 
-        super().__init__(sensitivity_simulation, parameters)
+        super().__init__(sensitivity_simulation, parameters, results_path)
 
         self.difference: float = difference
 
@@ -345,6 +354,7 @@ def create_samples(self) -> None:
         samples[-1, :] = reference_values # reference
 
         self.samples = samples
+        console.print(self.samples)
 
     def calculate_sensitivity(self):
         """Calculate the two-sided local sensitivity matrix."""
@@ -395,9 +405,12 @@ class SobolSensitivityAnalysis(SensitivityAnalysis):
     def __init__(self,
                  sensitivity_simulation: SensitivitySimulation,
                  parameters: list[SensitivityParameter],
+                 N: int,
+                 results_path: Path,
                  ):
 
-        super().__init__(sensitivity_simulation, parameters)
+        super().__init__(sensitivity_simulation, parameters, results_path)
+        self.N: int = N
 
         # define the problem specification
         self.ssa_problem: ProblemSpec = ProblemSpec({
@@ -409,7 +422,7 @@ def __init__(self,
         # console.print(self.ssa_problem)
 
 
-    def create_samples(self, N: int=1024):
+    def create_samples(self) -> None:
         """Create samples for sobol.
 
         Generates model inputs using Saltelli's extension of the Sobol' sequence
@@ -419,15 +432,14 @@ def create_samples(self, N: int=1024):
         """
 
         # libsa samples based on definition
-        ssa_samples = saltelli.sample(self.ssa_problem, N=N, calc_second_order=True)
+        ssa_samples = saltelli.sample(self.ssa_problem, N=self.N, calc_second_order=True)
         self.ssa_problem.set_samples(ssa_samples)
 
         # (num_samples x num_outputs)
         #  total model evaluations are (2d+2) * N for d input factors
-        num_samples = (2 * self.num_parameters + 2) * N
+        num_samples = (2 * self.num_parameters + 2) * self.N
 
         self.samples = xr.DataArray(
-            # np.full((num_samples, self.num_parameters), np.nan),
             ssa_samples,
             dims=["sample", "parameter"],
             coords={"sample": range(num_samples),
@@ -436,7 +448,7 @@ def create_samples(self, N: int=1024):
         )
 
 
-    def calculate_sensitivity(self):
+    def calculate_sensitivity(self) -> None:
         """Calculate the sensitivity matrices."""
 
         Y = self.results.values
@@ -465,16 +477,38 @@ def calculate_sensitivity(self):
             Si = SALib.analyze.sobol.analyze(
                 self.ssa_problem, Yo,
                 calc_second_order=True,
-                print_to_console=True,
+                print_to_console=False,
+                n_processors=4,
             )
-            console.print("S1")
-            console.print(Si["S1"])
             for key in sensitivity_keys:
                 self.sensitivity[key][:, ko] = Si[key]
 
 
-    def plot(self):
-        Si.plot()
-        from matplotlib import pyplot as plt
-        plt.show()
+    def plot_sobol_indices(
+        self,
+        fig_path: Path,
+        ):
+        """Barplots for the Sobol indices.
+
+        """
+        parameter_labels: dict[str, str] = {p.uid: f"{p.uid}: {p.name}" for p in self.parameters}
+        output_labels: dict[str, str] = {q.uid: q.name for q in self.outputs}
+
+
+        for ko, output in enumerate(self.outputs):
+            S1 = self.sensitivity["S1"][:, ko]
+            ST = self.sensitivity["ST"][:, ko]
+            S1_conf = self.sensitivity["S1_conf"][:, ko]
+            ST_conf = self.sensitivity["ST_conf"][:, ko]
+            console.print(S1)
+            console.print(type(S1))
+
+            break
+
+
+
+
+
+
+
 

src/sbmlsim/sensitivity/plots.py

@@ -3,6 +3,8 @@
 FIXME: use patchcollection
 https://stackoverflow.com/questions/59381273/heatmap-with-circles-indicating-size-of-population
 """
+from pathlib import Path
+
 import xarray as xr
 from typing import Optional
 
@@ -26,6 +28,7 @@ def heatmap(
     vcenter: float = 0.0,
     vmin: float = -2.0,
     vmax: float = 2.0,
+    fig_path: Optional[Path] = None,
 ):
     """Creates heatmap of model sensitivity"""
 
@@ -63,12 +66,12 @@ def calculate_subset(df, cutoff=0.01) -> pd.DataFrame:
     # parameters
     yticklabels = [pid for pid in df_subset.index]
     if parameter_labels:
-        yticklabels = [f"{pid}: {parameter_labels[pid]}" for pid in yticklabels]
+        yticklabels = [parameter_labels[pid] for pid in yticklabels]
 
     n_outputs = df_subset.shape[1]
     n_parameters = df_subset.shape[0]
-    figsize = (int(n_outputs/n_parameters*15), 15)
-
+    # (width, height)
+    figsize = (15, int(n_parameters/n_outputs*15))
 
     # plot heatmap
     cg = sns.clustermap(
@@ -79,7 +82,6 @@ def calculate_subset(df, cutoff=0.01) -> pd.DataFrame:
         xticklabels=xticklabels,
         yticklabels=yticklabels,
         cmap=cmap,
-        # cbar_pos=(0.0, 0.0, 0.6, 0.05), #  (left, bottom, width, height),
         cbar_pos=(0.0, 0.4, 0.03, 0.2),  # (left, bottom, width, height),
         cbar_kws={
             "orientation": "vertical",
@@ -100,13 +102,20 @@ def calculate_subset(df, cutoff=0.01) -> pd.DataFrame:
         horizontalalignment="right",
         size=20,
     )
-    label_fontsize=13
-    plt.setp(cg.ax_heatmap.get_yticklabels(), size=label_fontsize)
-    plt.setp(cg.ax_heatmap.get_xticklabels(), size=label_fontsize)
+    label_fontsize=15
+    plt.setp(cg.ax_heatmap.get_yticklabels(), size=label_fontsize, weight="bold")
+    plt.setp(cg.ax_heatmap.get_xticklabels(), size=label_fontsize, weight="bold")
     cg.ax_cbar.tick_params(labelsize=label_fontsize)
     cg.ax_row_dendrogram.set_visible(False)
     cg.ax_col_dendrogram.set_visible(False)
 
     if title:
-        plt.suptitle(title)
+        plt.suptitle(title, fontsize=40, fontweight="bold")
+
+    if fig_path:
+        plt.savefig(fig_path, dpi=300, bbox_inches="tight")
+    plt.show()
+
 
+def sobol_barplot():
+    pass