cleanup sensitivity

Author: matthiaskoenig

pyproject.toml

@@ -50,7 +50,7 @@ dependencies = [
     # dependencies
     "dill>=0.4.0",
     "SAlib>=1.5.2",
-    "petab>=0.7.0",
+    "petab>=0.8.1",
     # "tables>=3.10.2",
     "statsmodels>=0.14.6",
     "typst>=0.14.5",

src/sbmlsim/sensitivity/__init__.py

@@ -1,17 +1,45 @@
-"""Module for sensitivity analysis.
+"""
+Sensitivity analysis framework for computational models.
+
+This package provides a unified framework for analyzing how uncertainty and
+variability in model parameters affect model outputs. It supports multiple
+complementary sensitivity analysis strategies, including local, sampling-based,
+and global methods, enabling both qualitative and quantitative assessment of
+parameter influence.
+
+Sensitivity analyses are performed by systematically perturbing or sampling
+model parameters, executing simulations, and evaluating changes in selected
+model outputs. The framework is designed for deterministic simulation models
+and integrates sampling, caching, statistical evaluation, and visualization
+within a consistent workflow.
 
 TODO implementation of alternative methods:
     - [ ] FAST
     - [ ] Morris
 
-FIXME: generate simple example => create tests for the sensitivity
-FIXME: general documentation
+FIXME: generate simple example
+FIXME: create unittests for the sensitivity
 FIXME: add a flag to control resources for parallelization (ncores)
 
 """
-
-from .sensitivity_sobol import SobolSensitivityAnalysis
-from .sensitivity_sampling import SamplingSensitivityAnalysis
+from .analysis import (
+    SensitivitySimulation,
+    SensitivityOutput,
+    AnalysisGroup,
+)
+from .parameters import (
+    SensitivityParameter,
+)
 from .sensitivity_local import LocalSensitivityAnalysis
+from .sensitivity_sampling import SamplingSensitivityAnalysis
+from .sensitivity_sobol import SobolSensitivityAnalysis
 
-
+__all__ = [
+    "SensitivityParameter",
+    "SensitivitySimulation",
+    "SensitivityOutput",
+    "AnalysisGroup",
+    "SobolSensitivityAnalysis",
+    "SamplingSensitivityAnalysis",
+    "LocalSensitivityAnalysis",
+]

src/sbmlsim/sensitivity/analysis.py

@@ -2,28 +2,33 @@
 
 
 """
+import multiprocessing
 import os
 import time
-import multiprocessing
 from dataclasses import dataclass
-from typing import Optional, Any
 from pathlib import Path
-from rich.progress import track
-from pymetadata.console import console
+from typing import Optional, Any
 
+import dill
 import numpy as np
 import pandas as pd
-import xarray as xr
-
-
 import roadrunner
-import dill
+import xarray as xr
+from pymetadata.console import console
+from rich.progress import track
 
 from sbmlsim.sensitivity.parameters import SensitivityParameter
-from sbmlsim.sensitivity.outputs import SensitivityOutput
 from sbmlsim.sensitivity.plots import heatmap
 
 
+@dataclass
+class SensitivityOutput:
+    """Output measurement for SensitivityAnalysis."""
+    uid: str
+    name: str
+    unit: Optional[str]
+
+
 @dataclass
 class AnalysisGroup:
     """Subgroup for analysis."""
@@ -59,8 +64,8 @@ def __init__(self, model_path: Path, selections: list[str],
         outputs_dict = {q.uid for q in self.outputs}
         for key in y:
             if key not in outputs_dict:
-                raise ValueError(f"Key '{key}' missing in outputs dictionary: '{outputs_dict}")
-
+                raise ValueError(
+                    f"Key '{key}' missing in outputs dictionary: '{outputs_dict}")
 
     @staticmethod
     def load_model(model_path: Path, selections: list[str]) -> roadrunner.RoadRunner:
@@ -72,7 +77,8 @@ def load_model(model_path: Path, selections: list[str]) -> roadrunner.RoadRunner
         return rr
 
     @staticmethod
-    def apply_changes(r: roadrunner.RoadRunner, changes: dict[str, float], reset_all: bool=True) -> None:
+    def apply_changes(r: roadrunner.RoadRunner, changes: dict[str, float],
+                      reset_all: bool = True) -> None:
         """Apply changes after possible reset of the model."""
         if reset_all:
             r.resetAll()
@@ -81,7 +87,8 @@ def apply_changes(r: roadrunner.RoadRunner, changes: dict[str, float], reset_all
             # print(f"{key=} {value=}")
             r.setValue(key, value)
 
-    def simulate(self, r: roadrunner.RoadRunner, changes: dict[str, float]) -> dict[str, float]:
+    def simulate(self, r: roadrunner.RoadRunner, changes: dict[str, float]) -> dict[
+        str, float]:
         """Run a model simulation and return scalar results dictionary."""
 
         raise NotImplemented
@@ -115,7 +122,7 @@ def __init__(self,
                  parameters: list[SensitivityParameter],
                  groups: list[AnalysisGroup],
                  results_path: Path,
-                 seed: Optional[int]=None,
+                 seed: Optional[int] = None,
                  ) -> None:
         """Create a sensitivity analysis for given parameter ids.
 
@@ -161,20 +168,15 @@ def __init__(self,
 
         # multiple sensitivities are stored
         # sensitivity matrix; shape: (num_parameters x num_outputs); could be multiple
-        self.sensitivity: dict[str, dict[str, xr.DataArray]] = {g.uid: {} for g in self.groups}
+        self.sensitivity: dict[str, dict[str, xr.DataArray]] = {g.uid: {} for g in
+                                                                self.groups}
 
     def samples_table(self) -> pd.DataFrame:
         return self._data_table(d=self.samples)
 
     def results_table(self) -> pd.DataFrame:
         return self._data_table(d=self.results)
 
-    # def sensitivity_tables(self) -> dict[str, pd.DataFrame]:
-    #
-    #     for group in self.groups:
-    #         for key in group.changes.keys():
-    #             return {k: self._data_table(d=d) for k, d in self.sensitivity.items()}
-
     def _data_table(self, d: dict[str, xr.DataArray]) -> pd.DataFrame:
         items = []
         for group in self.groups:
@@ -188,7 +190,8 @@ def _data_table(self, d: dict[str, xr.DataArray]) -> pd.DataFrame:
         return pd.DataFrame(items)
 
     def read_cache(self, cache_filename: str, cache: bool) -> Optional[Any]:
-        cache_path: Optional[Path] = self.results_path / cache_filename if cache_filename else None
+        cache_path: Optional[
+            Path] = self.results_path / cache_filename if cache_filename else None
         if cache and not cache_path:
             raise ValueError("Cache path is required for caching.")
 
@@ -202,7 +205,8 @@ def read_cache(self, cache_filename: str, cache: bool) -> Optional[Any]:
         return None
 
     def write_cache(self, data: Any, cache_filename: str, cache: bool) -> Optional[Any]:
-        cache_path: Optional[Path] = self.results_path / cache_filename if cache_filename else None
+        cache_path: Optional[
+            Path] = self.results_path / cache_filename if cache_filename else None
         if cache_path:
             with open(cache_path, 'wb') as f:
                 console.print(f"Simulated samples written to cache: '{cache_path}'")
@@ -247,7 +251,8 @@ def num_samples(self) -> int:
         samples = self.samples[self.group_ids[0]]
         return samples.shape[0]
 
-    def simulate_samples(self, cache_filename: Optional[str] = None, cache: bool = False) -> None:
+    def simulate_samples(self, cache_filename: Optional[str] = None,
+                         cache: bool = False) -> None:
         """Simulate all samples in parallel.
 
         :param cache_filename: Path to the cache path.
@@ -320,7 +325,8 @@ def split_into_chunks(items, n):
         # write to cache
         self.write_cache(data=self.results, cache_filename=cache_filename, cache=cache)
 
-    def calculate_sensitivity(self, cache_filename: Optional[str] = None, cache: bool = False):
+    def calculate_sensitivity(self, cache_filename: Optional[str] = None,
+                              cache: bool = False):
         """Calculate the sensitivity matrices."""
 
         raise NotImplemented
@@ -361,14 +367,14 @@ def plot_sensitivity(
         )
 
 
-
 def run_simulation(
     params_tuple
 ):
     """Pass all required arguments as parameter tuple."""
     sensitivity_simulation, r, chunked_changes = params_tuple
     outputs = []
-    for kc in track(range(len(chunked_changes)), description=f"Simulate samples PID={os.getpid()}"):
+    for kc in track(range(len(chunked_changes)),
+                    description=f"Simulate samples PID={os.getpid()}"):
         changes = chunked_changes[kc]
         # console.print(f"PID={os.getpid()} | k={kc}")
         Y = sensitivity_simulation.simulate(
@@ -378,6 +384,3 @@ def run_simulation(
         outputs.append(Y)
 
     return outputs
-
-
-

src/sbmlsim/sensitivity/example/sensitivity_example.py

@@ -3,7 +3,6 @@
 
 import numpy as np
 import roadrunner
-from roadrunner._roadrunner import NamedArray
 from pymetadata.console import console
 
 from sbmlsim.sensitivity.analysis import (
@@ -13,13 +12,11 @@
 )
 from sbmlsim.sensitivity.parameters import (
     SensitivityParameter,
-    ParameterType,
     parameters_for_sensitivity_analysis,
 )
 
 model_path: Path = Path(__file__).parent / "simple_chain.xml"
 
-
 # Subgroups to perform sensitivity analysis on
 sensitivity_groups: list[AnalysisGroup] = [
     AnalysisGroup(
@@ -48,22 +45,24 @@ class ExampleSensitivitySimulation(SensitivitySimulation):
     tend = 1000  #
     steps = 1000
 
-    def simulate(self, r: roadrunner.RoadRunner, changes: dict[str, float]) -> dict[str, float]:
+    def simulate(self, r: roadrunner.RoadRunner, changes: dict[str, float]) -> dict[
+        str, float]:
 
         # apply changes and simulate
         all_changes = {
             **self.changes_simulation,  # model
             **changes  # sensitivity
         }
         self.apply_changes(r, all_changes, reset_all=True)
-        # ensure tolerances
+
+        # ensure identical tolerances on all simulations
         r.integrator.setValue("absolute_tolerance", self.init_tolerances)
-        s: NamedArray = r.simulate(start=0, end=self.tend, steps=self.steps)
+        s = r.simulate(start=0, end=self.tend, steps=self.steps)
 
-        # pharmacokinetic parameters
+        # calculate outputs y (custom functions)
+        # this can be registered functions calculating scalars based on subsets of the
+        # timecourse vectors
         y: dict[str, float] = {}
-
-        # calculate outputs (custom functions)
         t = s["time"]
         for key in "S1", "S2", "S3":
             rr_key = f"[{key}]"
@@ -79,13 +78,8 @@ def simulate(self, r: roadrunner.RoadRunner, changes: dict[str, float]) -> dict[
 
 sensitivity_simulation = ExampleSensitivitySimulation(
     model_path=model_path,
-    selections=[
-        "time",
-        "[S1]",
-        "[S2]",
-        "[S3]",
-    ],
-    changes_simulation = {},
+    selections=["time", "[S1]", "[S2]", "[S3]"],
+    changes_simulation={},
     outputs=[
         SensitivityOutput(uid='[S1]_auc', name='[S1] AUC', unit=None),
         SensitivityOutput(uid='[S2]_tmax', name='[S2] time maximum', unit=None),
@@ -120,49 +114,46 @@ def _sensitivity_parameters() -> list[SensitivityParameter]:
 
 sensitivity_parameters = _sensitivity_parameters()
 
-
 if __name__ == "__main__":
-
     from sbmlsim.sensitivity import (
         LocalSensitivityAnalysis,
         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,
-    # )
-    #
-    # 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,
-    # )
+    SamplingSensitivityAnalysis.run_sensitivity_analysis(
+        results_path=sensitivity_path / "sampling",
+        sensitivity_simulation=sensitivity_simulation,
+        parameters=sensitivity_parameters,
+        groups=sensitivity_groups,
+        cache_results=True,
+        cache_sensitivity=True,
+        N=1000,
+        seed=1234,
+    )
+
+    LocalSensitivityAnalysis.run_sensitivity_analysis(
+        results_path=sensitivity_path / "local",
+        sensitivity_simulation=sensitivity_simulation,
+        parameters=sensitivity_parameters,
+        groups=[sensitivity_groups[1]],
+        cache_results=True,
+        cache_sensitivity=True,
+        difference=0.01,
+        seed=1234,
+    )
 
     SobolSensitivityAnalysis.run_sensitivity_analysis(
         results_path=sensitivity_path / "sobol",
         sensitivity_simulation=sensitivity_simulation,
         parameters=sensitivity_parameters,
         groups=[sensitivity_groups[1]],
-        # cache_results=False,
-        # cache_sensitivity=False,
-        # N=2048,
-        N=8,
+        cache_results=True,
+        cache_sensitivity=True,
+        N=2048,
+        # N=8,
         seed=1234,
     )
-
-