working on sensitivity

Author: matthiaskoenig

src/sbmlsim/sensitivity/global_sensitivity.py

@@ -19,6 +19,8 @@
 - [ ] parallelization ? (benchmark)
 
 """
+from typing import Optional
+
 import SALib
 from SALib import ProblemSpec
 from SALib.sample import saltelli
@@ -35,7 +37,6 @@
 from roadrunner._roadrunner import NamedArray
 
 
-
 @dataclass
 class SensitivitySimulation:
     """Base class for sensitivity calculation.
@@ -49,8 +50,9 @@ class SensitivitySimulation:
     selections: list[str]
     rr: roadrunner.RoadRunner = None
     outputs: list[str] = None
+    changes_simulation: dict[str, float] = None
 
-    def __init__(self, model_path: Path, selections: list[str]):
+    def __init__(self, model_path: Path, selections: list[str], changes_simulation: dict[str, float]):
         self.model_path = model_path
         self.selections = selections
         self.rr: roadrunner.RoadRunner = roadrunner.RoadRunner(str(model_path))
@@ -59,6 +61,9 @@ def __init__(self, model_path: Path, selections: list[str]):
         integrator.setSetting("variable_step_size", True)
         # state = rr.saveStateS()
 
+        # store the simulation changes
+        self.changes_simulation = changes_simulation
+
         # get the outputs from the simulation
         y = self.simulate(changes={})
         self.outputs = list(y.keys())
@@ -71,6 +76,15 @@ def simulate(self, changes: dict[str, float]) -> dict[str, float]:
         """
         raise NotImplemented
 
+    def parameter_values(self, parameters: list[str], changes: dict[str, float]) -> dict[str, float]:
+        """Get the parameter values for a given set of changes."""
+        self.apply_changes(changes, reset_all=True)
+        values: dict[str, float] = {}
+        for pid in parameters:
+            values[pid] = self.rr.getValue(pid)
+        return values
+
+
     def plot(self) -> None:
         """Plots the model simulation for debugging."""
         raise NotImplemented
@@ -84,9 +98,116 @@ def apply_changes(self, changes: dict[str, float], reset_all: bool=True) -> None
             self.rr.setValue(key, value)
 
 
+@dataclass
+class SensitivityAnalysis:
+    """Parent class for all sensitivity analysis.
+
+    TODO: additional metadata for the outputs and the parameters; i.e. name, units, bounds, ....
+    """
+
+    sensitivity_simulation: SensitivitySimulation
+
+    def __init__(self, sensitivity_simulation: SensitivitySimulation,
+                 parameters: list[str]) -> None:
+        """Create a sensitivity analysis for given parameter ids.
+
+        Based on the results matrix the sensitivity is calculated.
+        """
+        self.sensitivity_simulation = sensitivity_simulation
+
+        # parameters to vary; shape: (num_parameters,)
+        self.parameters: list[str] = parameters
+        # outputs to calculate sensitivity on; shape: (num_outputs,)
+        self.outputs: list[str] = sensitivity_simulation.outputs
+        # parameter samples for sensitivity; shape: (num_samples x num_parameters)
+        self.samples: Optional[np.ndarray] = None
+        # outputs for given samples; shape: (num_samples x num_outputs)
+        self.results: Optional[np.ndarray] = None
+        # sensitivity matrix; shape: (num_parameters x num_outputs)
+        self.sensitivity_results: Optional[np.ndarray] = None
+
+    @property
+    def num_parameters(self) -> int:
+        return len(self.parameters)
+
+    @property
+    def num_outputs(self) -> int:
+        return len(self.outputs)
+
+    def create_samples(self) -> None:
+        """Create and set parameter samples."""
+
+        raise NotImplemented
+
+    def num_samples(self) -> int:
+        """Number of samples.
+
+        Requires that samples have been created.
+        """
+        return self.samples.shape[0]
+
+    def simulate_samples(self) -> None:
+        """Simulate all samples."""
+        self.samples = np.zeros(shape=(self.num_samples, self.num_parameters))
+        self.outputs = np.zeros(shape=(self.num_samples, self.num_outputs))
+
+        for k in range(self.num_samples()):
+            changes = dict(zip(self.parameters, self.samples[k, :]))
+            outputs = self.sensitivity_simulation.simulate(changes=changes)
+            self.outputs[k, :] = outputs
+
+    def calculate_sensitivity(self):
+        """Calculate the sensitivity matrix."""
+
+        raise NotImplemented
+
+@dataclass
+class LocalSensitivityAnalysis(SensitivityAnalysis):
+    """Local sensitivity analysis based on local differences."""
+
+    difference: float
+    sensitivity: np.ndarray = None
+
+    def __init__(self, sensitivity_simulation: SensitivitySimulation,
+                 parameters: list[str], difference: float = 0.1):
+
+        self.sensitivity = np.zeros(shape=(self.num_parameters, self.num_outputs))
+        self.difference = difference
+        self.samples = self.create_samples()
+
+    @property
+    def num_samples(self) -> int:
+        """Number of parameter samples to simulate."""
+        return 2 * self.num_parameters
+
+    def create_samples(self) -> np.ndarray:
+
+        for key, value in p_ref.items():
+            values = np.ones(shape=(2 * num_pars,)) * value.magnitude
+            # change parameters in correct position
+            values[index] = value.magnitude * (1.0 + difference)
+            values[index + num_pars] = value.magnitude * (1.0 - difference)
+            changes[key] = Q_(values, value.units)
+            index += 1
+
+    def calculate_sensitivity(self):
+
+        pass
+
+    def plot_sensitivity(self):
+
+        pass
+
+
+@dataclass
+class SamplingSensitivityAnalysis(SensitivityAnalysis):
+    """Sample from provided parameter distributions."""
+
+    # TODO: implement
+    pass
 
 @dataclass
-class SBMLSensitivityAnalysis:
+class GlobalSobolSensitivityAnalysis:
     """Parent class for sensitivity analysis."""
 
     sensitivity_simulation: SensitivitySimulation