updated sensitivity analysis

Author: matthiaskoenig

src/sbmlsim/sensitivity/example/sensitivity_example.py

@@ -128,29 +128,39 @@ def _sensitivity_parameters() -> list[SensitivityParameter]:
         SobolSensitivityAnalysis,
         SamplingSensitivityAnalysis,
     )
-
     sensitivity_path = Path(__file__).parent / "results"
     console.print(SensitivityParameter.parameters_to_df(sensitivity_parameters))
 
-    SamplingSensitivityAnalysis.run_sensitivity_analysis(
-        results_path=sensitivity_path / "sampling",
-        sensitivity_simulation=sensitivity_simulation,
-        parameters=sensitivity_parameters,
-        groups=sensitivity_groups,
-        # cache_results=False,
-        # cache_sensitivity=False,
-        N=200,
-        seed=1234,
-    )
+    # SamplingSensitivityAnalysis.run_sensitivity_analysis(
+    #     results_path=sensitivity_path / "sampling",
+    #     sensitivity_simulation=sensitivity_simulation,
+    #     parameters=sensitivity_parameters,
+    #     groups=sensitivity_groups,
+    #     # cache_results=False,
+    #     # cache_sensitivity=False,
+    #     N=200,
+    #     seed=1234,
+    # )
+    #
+    # LocalSensitivityAnalysis.run_sensitivity_analysis(
+    #     results_path=sensitivity_path / "local",
+    #     sensitivity_simulation=sensitivity_simulation,
+    #     parameters=sensitivity_parameters,
+    #     groups=[sensitivity_groups[1]],
+    #     # cache_results=False,
+    #     # cache_sensitivity=False,
+    #     difference=0.01,
+    #     seed=1234,
+    # )
 
     SobolSensitivityAnalysis.run_sensitivity_analysis(
         results_path=sensitivity_path / "sobol",
         sensitivity_simulation=sensitivity_simulation,
         parameters=sensitivity_parameters,
-        groups=sensitivity_groups,
+        groups=[sensitivity_groups[1]],
         # cache_results=False,
         # cache_sensitivity=False,
-        # N=2048
+        # N=2048,
         N=8,
         seed=1234,
     )

src/sbmlsim/sensitivity/sensitivity_local.py

@@ -1,6 +1,9 @@
 from pathlib import Path
+from typing import Optional
 
+import numpy as np
 import xarray as xr
+from pymetadata.console import console
 
 from sbmlsim.sensitivity.analysis import SensitivitySimulation, AnalysisGroup, \
     SensitivityAnalysis
@@ -136,45 +139,61 @@ def calculate_sensitivity(self, cache_filename: Optional[str] = None, cache: boo
         # write to cache
         self.write_cache(data=self.sensitivity, cache_filename=cache_filename, cache=cache)
 
-
-def local_sensitivity_analysis():
-    """Local sensitivity analysis"""
-    console.rule("LOCAL SENSITIVITY ANALYSIS", style="blue bold", align="center")
-
-    sensitivity_simulation = CanagliflozinSensitivitySimulation.sensitivity_simulation()
-    parameters = sensitivity_simulation.sensitivity_parameters()
-
-    sa = LocalSensitivityAnalysis(
-        sensitivity_simulation=sensitivity_simulation,
-        parameters=sensitivity_parameters,
-        groups=sensitivity_groups,
-        results_path=RESULTS_PATH / "sensitivity",
-        seed=1234,
-        difference=0.01,  # 1% change
-    )
-
-    console.rule("Samples", style="white")
-    sa.create_samples()
-
-    console.rule("Results", style="white")
-    sa.simulate_samples()
-    console.print(sa.results)
-
-    console.rule("Sensitivity", style="white")
-    sa.calculate_sensitivity()
-    console.print(sa.sensitivity)
-
-    console.rule("Plotting", style="white")
-    for kg, group in enumerate(sa.groups):
-        sa.plot_sensitivity(
-            group_id=group.uid,
-            sensitivity_key="normalized",
-            # title=f"{group.name}",
-            cutoff=0.05,
-            cluster_rows = False,
-            cmap = "seismic",
-            vcenter=0.0,
-            vmin=-2.0,
-            vmax=2.0,
-            fig_path=sa.results_path / f"local_sensitivity_{kg:>02}_{group.uid}_{sa.difference}.png",
+    @staticmethod
+    def run_sensitivity_analysis(
+            results_path: Path,
+            sensitivity_simulation: SensitivitySimulation,
+            parameters: list[SensitivityParameter],
+            groups: list[AnalysisGroup],
+            seed: int,
+            difference: float = 0.01,
+            cache_results: bool = False,
+            cache_sensitivity: bool = False,
+    ) -> None:
+        """Local sensitivity analysis.
+
+        :param sensitivity_simulation: Sensitivity simulation.
+        :param parameters: Sensitivity parameters.
+        :param groups: Sensitivity groups.
+        :param difference: relative change of parameters.
+        :param seed: Random seed.
+        """
+        console.rule("LOCAL SENSITIVITY ANALYSIS", style="blue bold", align="center")
+        if cache_sensitivity and not cache_results:
+            # sensitivity must be recalculated for new results
+            cache_sensitivity = False
+
+        sa = LocalSensitivityAnalysis(
+            sensitivity_simulation=sensitivity_simulation,
+            parameters=parameters,
+            groups=groups,
+            results_path=results_path,
+            seed=1234,
+            difference=difference,
         )
+        console.rule("Samples", style="white")
+        sa.create_samples()
+        console.print(sa.samples_table())
+
+        console.rule("Results", style="white")
+        sa.simulate_samples(cache_filename=f"local_results_difference{sa.difference}.pkl", cache=cache_results)
+        console.print(sa.results_table())
+
+        console.rule("Sensitivity", style="white")
+        sa.calculate_sensitivity(cache_filename=f"local_sensitivity_difference{sa.difference}.pkl", cache=cache_sensitivity)
+        # console.print(sa.sensitivity_tables())
+
+        console.rule("Plotting", style="white")
+        for kg, group in enumerate(sa.groups):
+            sa.plot_sensitivity(
+                group_id=group.uid,
+                sensitivity_key="normalized",
+                # title=f"{group.name}",
+                cutoff=0.05,
+                cluster_rows=False,
+                cmap="seismic",
+                vcenter=0.0,
+                vmin=-2.0,
+                vmax=2.0,
+                fig_path=sa.results_path / f"local_sensitivity_{kg:>02}_{group.uid}_{sa.difference}.png",
+            )

src/sbmlsim/sensitivity/sensitivity_sobol.py

@@ -64,7 +64,6 @@ def create_samples(self) -> None:
         The Sobol' sequence is a popular quasi-random low-discrepancy sequence used
         to generate uniform samples of parameter space.
         """
-        console.rule("Samples", style="white")
         # (num_samples x num_outputs)
         #  total model evaluations are (2d+2) * N for d input factors
         num_samples = (2 * self.num_parameters + 2) * self.N
@@ -81,7 +80,6 @@ def create_samples(self) -> None:
                         "parameter": self.parameter_ids},
                 name="samples"
             )
-        console.print(self.samples)
 
 
     def calculate_sensitivity(self, cache_filename: Optional[str] = None, cache: bool = False):
@@ -128,40 +126,6 @@ def calculate_sensitivity(self, cache_filename: Optional[str] = None, cache: boo
         # write to cache
         self.write_cache(data=self.sensitivity, cache_filename=cache_filename, cache=cache)
 
-    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}
-        parameter_labels: dict[str, str] = {p.uid: p.uid for p in self.parameters}
-        output_labels: dict[str, str] = {q.uid: q.name for q in self.outputs}
-
-        for group in self.groups:
-            gid = group.uid
-            ymax = self.sensitivity[gid]["ST"].max(dim=None)
-            ymin = self.sensitivity[gid]["S1"].min(dim=None)
-
-            for ko, output in enumerate(self.outputs):
-                # f_path = fig_path.parent / f"FigS{ko+22}_{fig_path.stem}_{ko:>03}_{output.uid}{fig_path.suffix}"
-                f_path = fig_path.parent / f"{fig_path.stem}_{ko:>03}_{output.uid}{fig_path.suffix}"
-
-                S1 = self.sensitivity[gid]["S1"][:, ko]
-                ST = self.sensitivity[gid]["ST"][:, ko]
-                S1_conf = self.sensitivity[gid]["S1_conf"][:, ko]
-                ST_conf = self.sensitivity[gid]["ST_conf"][:, ko]
-                sobol_barplot(
-                    S1=S1,
-                    ST=ST,
-                    S1_conf=S1_conf,
-                    ST_conf=ST_conf,
-                    title=f"{output_labels[output.uid]} ({group.name})",
-                    fig_path=f_path,
-                    parameter_labels=parameter_labels,
-                    ymax=np.max([1.05, ymax]),
-                    ymin=np.min([-0.05, ymin]),
-                )
-
     @staticmethod
     def run_sensitivity_analysis(
         results_path: Path,
@@ -192,9 +156,8 @@ def run_sensitivity_analysis(
             parameters=parameters,
             groups=groups,
             results_path=results_path,
-            # N=4096,
-            N=16,
-            seed=1234,
+            N=N,
+            seed=seed,
         )
 
         console.rule("Samples", style="white")
@@ -232,46 +195,79 @@ def run_sensitivity_analysis(
                 fig_path=sa.results_path / f"sobol_sensitivity_N{sa.N}_{kg:>02}_{group.uid}.png",
             )
 
+    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}
+        parameter_labels: dict[str, str] = {p.uid: p.uid for p in self.parameters}
+        output_labels: dict[str, str] = {q.uid: q.name for q in self.outputs}
 
+        for group in self.groups:
+            gid = group.uid
+            ymax = self.sensitivity[gid]["ST"].max(dim=None)
+            ymin = self.sensitivity[gid]["S1"].min(dim=None)
 
-def sobol_barplot(
-    S1, ST, S1_conf, ST_conf,
-    parameter_labels: dict[str, str],
-    fig_path: Optional[Path] = None,
-    title: Optional[str] = None,
-    ymax: float = 1.1,
-    ymin: float = -0.1,
-):
-    # width
-    figsize = (15, 3)
-    label_fontsize = 15
-
-    categories: list[str] = list(parameter_labels.values())
-    f, ax = plt.subplots(figsize=figsize)
-
-    ax.bar(categories, ST, label='ST',
-           color="tab:orange",
-           alpha=1.0,
-           edgecolor="black",
-           yerr=ST_conf, capsize=5
-           )
-
-    ax.bar(categories, S1, label='S1', color="tab:blue",
-           edgecolor="black", yerr=S1_conf, capsize=5)
-
-
-    # ax.set_xlabel('Parameter', fontsize=label_fontsize, fontweight="bold")
-    ax.set_ylabel('Sobol Index', fontsize=label_fontsize, fontweight="bold")
-    ax.set_ylim(bottom=ymin, top=ymax)
-    ax.grid(True, axis="y")
-    ax.tick_params(axis='x', labelrotation=90)
-    # ax.tick_params(axis='x', labelweight='bold')
-    ax.legend()
+            for ko, output in enumerate(self.outputs):
+                # f_path = fig_path.parent / f"FigS{ko+22}_{fig_path.stem}_{ko:>03}_{output.uid}{fig_path.suffix}"
+                f_path = fig_path.parent / f"{fig_path.stem}_{ko:>03}_{output.uid}{fig_path.suffix}"
 
-    if title:
-        plt.suptitle(title, fontsize=20, fontweight="bold")
+                S1 = self.sensitivity[gid]["S1"][:, ko]
+                ST = self.sensitivity[gid]["ST"][:, ko]
+                S1_conf = self.sensitivity[gid]["S1_conf"][:, ko]
+                ST_conf = self.sensitivity[gid]["ST_conf"][:, ko]
+                self.sobol_barplot(
+                    S1=S1,
+                    ST=ST,
+                    S1_conf=S1_conf,
+                    ST_conf=ST_conf,
+                    title=f"{output_labels[output.uid]} ({group.name})",
+                    fig_path=f_path,
+                    parameter_labels=parameter_labels,
+                    ymax=np.max([1.05, ymax]),
+                    ymin=np.min([-0.05, ymin]),
+                )
 
-    if fig_path:
-        plt.savefig(fig_path, dpi=300, bbox_inches="tight")
-    plt.show()
+    @staticmethod
+    def sobol_barplot(
+        S1, ST, S1_conf, ST_conf,
+        parameter_labels: dict[str, str],
+        fig_path: Optional[Path] = None,
+        title: Optional[str] = None,
+        ymax: float = 1.1,
+        ymin: float = -0.1,
+    ):
+        # width
+        figsize = (15, 3)
+        label_fontsize = 15
+
+        categories: list[str] = list(parameter_labels.values())
+        f, ax = plt.subplots(figsize=figsize)
+
+        ax.bar(categories, ST, label='ST',
+               color="tab:orange",
+               alpha=1.0,
+               edgecolor="black",
+               yerr=ST_conf, capsize=5
+               )
+
+        ax.bar(categories, S1, label='S1', color="tab:blue",
+               edgecolor="black", yerr=S1_conf, capsize=5)
+
+
+        # ax.set_xlabel('Parameter', fontsize=label_fontsize, fontweight="bold")
+        ax.set_ylabel('Sobol Index', fontsize=label_fontsize, fontweight="bold")
+        ax.set_ylim(bottom=ymin, top=ymax)
+        ax.grid(True, axis="y")
+        ax.tick_params(axis='x', labelrotation=90)
+        # ax.tick_params(axis='x', labelweight='bold')
+        ax.legend()
+
+        if title:
+            plt.suptitle(title, fontsize=20, fontweight="bold")
+
+        if fig_path:
+            plt.savefig(fig_path, dpi=300, bbox_inches="tight")
+        plt.show()