run sensitivity analysis

Author: matthiaskoenig

src/sbmlsim/sensitivity/__init__.py

@@ -12,14 +12,6 @@
 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:
-    - [ ] Morris
-
-FIXME: generate simple example
-FIXME: create unittests for the sensitivity
-FIXME: add a flag to control resources for parallelization (ncores)
-
 """
 from .analysis import (
     SensitivityAnalysis,
@@ -30,10 +22,10 @@
 from .parameters import (
     SensitivityParameter,
 )
+from .sensitivity_fast import FASTSensitivityAnalysis
 from .sensitivity_local import LocalSensitivityAnalysis
 from .sensitivity_sampling import SamplingSensitivityAnalysis
 from .sensitivity_sobol import SobolSensitivityAnalysis
-from .sensitivity_fast import FASTSensitivityAnalysis
 
 __all__ = [
     "SensitivityParameter",

src/sbmlsim/sensitivity/analysis.py

@@ -123,6 +123,8 @@ def __init__(self,
                  groups: list[AnalysisGroup],
                  results_path: Path,
                  seed: Optional[int] = None,
+                 n_cores: Optional[int] = None,
+                 cache_results: bool = False,
                  ) -> None:
         """Create a sensitivity analysis for given parameter ids.
 
@@ -160,6 +162,15 @@ def __init__(self,
         if seed is not None:
             np.random.seed(seed)
 
+        # caching
+        self.cache_results: bool = cache_results
+        self.prefix: str = self.__class__.__name__
+
+        # handle compute resources
+        if not n_cores:
+            n_cores = int(round(0.9 * multiprocessing.cpu_count()))
+        self.n_cores = n_cores
+
         # parameter samples for sensitivity; shape: (num_samples x num_parameters)
         self.samples: dict[str, Optional[xr.DataArray]] = {}
 
@@ -171,47 +182,6 @@ def __init__(self,
         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 _data_table(self, d: dict[str, xr.DataArray]) -> pd.DataFrame:
-        items = []
-        for group in self.groups:
-            da: xr.DataArray = d[group.uid]
-            item = {
-                'group': group.uid,
-                # 'group_name': group.name,
-                **da.sizes,
-            }
-            items.append(item)
-        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
-        if cache and not cache_path:
-            raise ValueError("Cache path is required for caching.")
-
-        # retrieve from cache
-        if cache and cache_path.exists():
-            with open(cache_path, 'rb') as f:
-                data = dill.load(f)
-                console.print(f"Simulated samples loaded from cache: '{cache_path}'")
-                return data
-
-        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
-        if cache_path:
-            with open(cache_path, 'wb') as f:
-                console.print(f"Simulated samples written to cache: '{cache_path}'")
-                dill.dump(data, f)
-
     @property
     def output_ids(self) -> list[str]:
         return [o.uid for o in self.outputs]
@@ -236,6 +206,30 @@ def num_outputs(self) -> int:
     def num_groups(self) -> int:
         return len(self.groups)
 
+    def execute(self):
+        """Execute the sensitivity analysis."""
+        console.rule(
+            f"{self.__class__.__name__}",
+            style="blue bold",
+            align="center",
+        )
+        console.rule("Samples", style="white")
+        self.create_samples()
+        console.print(self.samples_table())
+
+        console.rule("Results", style="white")
+        self.simulate_samples(
+            cache_filename=f"{self.prefix}_results.pkl",
+            cache=self.cache_results,
+        )
+        console.print(self.results_table())
+
+        console.rule("Sensitivity", style="white")
+        self.calculate_sensitivity(
+            cache_filename=f"{self.prefix}_sensitivity.pkl",
+            cache=self.cache_results,
+        )
+
     def create_samples(self) -> None:
         """Create and set parameter samples."""
 
@@ -283,8 +277,6 @@ def simulate_samples(self, cache_filename: Optional[str] = None,
             )
 
             # number of cores
-            n_cores = multiprocessing.cpu_count()
-
             samples = self.samples[group.uid]
 
             # create chunk of samples for core
@@ -305,13 +297,13 @@ def split_into_chunks(items, n):
                 return chunks, chunked_samples
 
             items = list(range(self.num_samples))
-            chunks, chunked_samples = split_into_chunks(items, n_cores)
+            chunks, chunked_samples = split_into_chunks(items, self.n_cores)
 
             # parameters for multiprocessing
             sa_sim = self.sensitivity_simulation
-            rrs = [(sa_sim, r, chunked_samples[i]) for i in range(n_cores)]
+            rrs = [(sa_sim, r, chunked_samples[i]) for i in range(self.n_cores)]
 
-            with multiprocessing.Pool(processes=n_cores) as pool:
+            with multiprocessing.Pool(processes=self.n_cores) as pool:
                 outputs_list: list = pool.map(run_simulation, rrs)
 
             for kc, chunk in enumerate(chunks):
@@ -331,6 +323,47 @@ def calculate_sensitivity(self, cache_filename: Optional[str] = None,
 
         raise NotImplemented
 
+    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 _data_table(self, d: dict[str, xr.DataArray]) -> pd.DataFrame:
+        items = []
+        for group in self.groups:
+            da: xr.DataArray = d[group.uid]
+            item = {
+                'group': group.uid,
+                # 'group_name': group.name,
+                **da.sizes,
+            }
+            items.append(item)
+        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
+        if cache and not cache_path:
+            raise ValueError("Cache path is required for caching.")
+
+        # retrieve from cache
+        if cache and cache_path.exists():
+            with open(cache_path, 'rb') as f:
+                data = dill.load(f)
+                console.print(f"Simulated samples loaded from cache: '{cache_path}'")
+                return data
+
+        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
+        if cache_path:
+            with open(cache_path, 'wb') as f:
+                console.print(f"Simulated samples written to cache: '{cache_path}'")
+                dill.dump(data, f)
+
     def sensitivity_df(self, group_id: str, key: str) -> pd.DataFrame:
         """Convert sensitivity information to dataframes."""
 
@@ -341,6 +374,10 @@ def sensitivity_df(self, group_id: str, key: str) -> pd.DataFrame:
             index=sensitivity.coords["parameter"]
         )
 
+    def plot(self):
+        """Should be implemented by subclass."""
+        console.rule("Plotting", style="white")
+
     def plot_sensitivity(
         self,
         group_id: str,

src/sbmlsim/sensitivity/example/sensitivity_example.py

@@ -112,58 +112,63 @@ def _sensitivity_parameters() -> list[SensitivityParameter]:
 sensitivity_parameters = _sensitivity_parameters()
 
 if __name__ == "__main__":
+    import multiprocessing
     from sbmlsim.sensitivity import (
-        LocalSensitivityAnalysis,
         SobolSensitivityAnalysis,
+        LocalSensitivityAnalysis,
         SamplingSensitivityAnalysis,
         FASTSensitivityAnalysis,
     )
 
     sensitivity_path = Path(__file__).parent / "results"
     console.print(SensitivityParameter.parameters_to_df(sensitivity_parameters))
-    cache = False
 
-    SamplingSensitivityAnalysis.run_sensitivity_analysis(
-        results_path=sensitivity_path / "sampling",
+    settings = {
+        "cache_results": False,
+        "n_cores": int(round(0.9 * multiprocessing.cpu_count())),
+        "seed": 1234
+    }
+
+    sa_sampling = SamplingSensitivityAnalysis(
         sensitivity_simulation=sensitivity_simulation,
         parameters=sensitivity_parameters,
         groups=sensitivity_groups,
-        cache_results=cache,
-        cache_sensitivity=cache,
+        results_path=sensitivity_path / "sampling",
         N=1000,
-        seed=1234,
+        **settings,
     )
+    sa_sampling.execute()
+    sa_sampling.plot()
 
-    LocalSensitivityAnalysis.run_sensitivity_analysis(
-        results_path=sensitivity_path / "local",
+    sa_local = LocalSensitivityAnalysis(
         sensitivity_simulation=sensitivity_simulation,
         parameters=sensitivity_parameters,
-        groups=[sensitivity_groups[1]],
-        cache_results=cache,
-        cache_sensitivity=cache,
+        groups=sensitivity_groups,
+        results_path=sensitivity_path / "local",
         difference=0.01,
-        seed=1234,
+        **settings,
     )
+    sa_local.execute()
+    sa_local.plot()
 
-    SobolSensitivityAnalysis.run_sensitivity_analysis(
-        results_path=sensitivity_path / "sobol",
+    sa_sobol = SobolSensitivityAnalysis(
         sensitivity_simulation=sensitivity_simulation,
         parameters=sensitivity_parameters,
         groups=[sensitivity_groups[1]],
-        cache_results=cache,
-        cache_sensitivity=cache,
+        results_path=sensitivity_path / "sobol",
         N=4096,
-        # N=8,
-        seed=1234,
+        **settings,
     )
+    sa_sobol.execute()
+    sa_sobol.plot()
 
-    FASTSensitivityAnalysis.run_sensitivity_analysis(
-        results_path=sensitivity_path / "fast",
+    sa_fast = FASTSensitivityAnalysis(
         sensitivity_simulation=sensitivity_simulation,
         parameters=sensitivity_parameters,
-        groups=[sensitivity_groups[1]],
-        cache_results=cache,
-        cache_sensitivity=cache,
+        groups=sensitivity_groups,
+        results_path=sensitivity_path / "fast",
         N=1000,
-        seed=1234,
+        **settings,
     )
+    sa_fast.execute()
+    sa_fast.plot()

src/sbmlsim/sensitivity/plots.py

@@ -142,6 +142,7 @@ def plot_S1_ST_indices(
                 ymin=np.min([-0.05, ymin]),
             )
 
+
 def S1_ST_barplot(
     S1, ST, S1_conf, ST_conf,
     parameter_labels: dict[str, str],