work on sensitivity

Author: matthiaskoenig

src/sbmlsim/sensitivity/analysis.py

@@ -46,7 +46,7 @@
 
 from sbmlsim.sensitivity.parameters import SensitivityParameter
 from sbmlsim.sensitivity.outputs import SensitivityOutput
-
+import pandas as pd
 
 @dataclass
 class SensitivitySimulation:
@@ -69,7 +69,7 @@ def __init__(self, model_path: Path, selections: list[str], changes_simulation:
         self.rr: roadrunner.RoadRunner = roadrunner.RoadRunner(str(model_path))
         self.rr.selections = self.selections
         integrator: roadrunner.Integrator = self.rr.integrator
-        integrator.setSetting("variable_step_size", True)
+        # integrator.setSetting("variable_step_size", True)
         # state = rr.saveStateS()
 
         # store the simulation changes
@@ -148,6 +148,7 @@ def __init__(self, sensitivity_simulation: SensitivitySimulation,
         self.results: Optional[xr.DataArray] = None
         # sensitivity matrix; shape: (num_parameters x num_outputs); could be multiple
         self.sensitivity: Optional[xr.DataArray] = None
+        self.sensitivity_normalized: Optional[xr.DataArray] = None
 
 
     @property
@@ -192,25 +193,22 @@ def simulate_samples(self) -> None:
 
     def calculate_sensitivity(self):
         """Calculate the sensitivity matrix."""
+        pass
 
-        self.sensitivity = xr.DataArray(
-            np.full((self.num_parameters, self.num_outputs), np.nan),
-            dims=["parameter", "output"],
-            coords={"parameter": [p.uid for p in self.parameters],
-                    "output": self.outputs},
-            name="sensitivity"
-        )
 
 
 @dataclass
 class LocalSensitivityAnalysis(SensitivityAnalysis):
-    """Local sensitivity analysis based on local differences."""
+    """Local sensitivity analysis based on local differences.
+
+    param difference: change for calculation of local sensitivity (0.01 = 1% change)
+    """
 
     difference: float
     sensitivity: np.ndarray = None
 
     def __init__(self, sensitivity_simulation: SensitivitySimulation,
-                 parameters: list[SensitivityParameter], difference: float = 0.1):
+                 parameters: list[SensitivityParameter], difference: float = 0.01):
 
         super().__init__(sensitivity_simulation, parameters)
         self.sensitivity = np.zeros(shape=(self.num_parameters, self.num_outputs))
@@ -260,34 +258,62 @@ def calculate_sensitivity(self):
         """Calculate the two-sided local sensitivity matrix."""
 
         # num_parameters x num_outputs
-        super().calculate_sensitivity()
+        # empty sensitivity
+        self.sensitivity = xr.DataArray(
+            np.full((self.num_parameters, self.num_outputs), np.nan),
+            dims=["parameter", "output"],
+            coords={"parameter": [p.uid for p in self.parameters],
+                    "output": self.outputs},
+            name="sensitivity"
+        )
+        self.sensitivity_normalized = xr.DataArray(
+            np.full((self.num_parameters, self.num_outputs), np.nan),
+            dims=["parameter", "output"],
+            coords={"parameter": [p.uid for p in self.parameters],
+                    "output": self.outputs},
+            name="sensitivity"
+        )
 
         for kp, p in enumerate(self.parameters):
             pid = self.parameters[kp].uid
+            p_ref = self.samples[-1, kp]
+            p_up = self.samples[2*kp, kp]
+            p_down = self.samples[2 * kp + 1, kp]
+
             for ko, oid in enumerate(self.outputs):
                 # num_samples x num_outputs
-                value_ref = self.results[-1, ko]
-                value_up = self.results[2*kp, ko]
-                value_down = self.results[2 * kp + 1, ko]
+                q_ref = self.results[-1, ko]
+                q_up = self.results[2*kp, ko]
+                q_down = self.results[2 * kp + 1, ko]
 
-                # midpoint method, two-sided sensitivity
-                self.sensitivity[kp, ko] = (value_up - value_down) / (2.0 * value_ref)
+                # two-sided sensitivity
+                self.sensitivity[kp, ko] = (q_up - q_down) / (p_up - p_down)
+                # normalized: relative change in output per relative change in parameter
+                self.sensitivity_normalized[kp, ko] = self.sensitivity[kp, ko] * p_ref/q_ref
 
-
-    def plot_sensitivity(self):
-        from sbmlsim.sensitivity.plots import heatmap
-        import pandas as pd
-
-        df = pd.DataFrame(
-            self.sensitivity.values,
+    @property
+    def sensitivity_df(self) -> pd.DataFrame:
+        """Convert sensitivity information to dataframe."""
+        return pd.DataFrame(
+            self.sensitivity_normalized.values,
             columns=self.sensitivity.coords["output"],
             index=self.sensitivity.coords["parameter"]
         )
-        console.print(df)
-        heatmap(df, cutoff=0.01)
 
+    def plot_sensitivity(self):
+        df = self.sensitivity_df
+        self.plot_sensitivity_df(df)
+
+    @staticmethod
+    def plot_sensitivity_df(df: pd.DataFrame, cutoff=0.1, cluster_rows: bool = True):
+        from sbmlsim.sensitivity.plots import heatmap
+        console.print(df)
 
+        # TODO: labels of parameters
+        # TODO: labels of outputs
+        # TODO: better position of colorbar
 
+        heatmap(df, cutoff=cutoff, cluster_rows=False)
 
 
 

src/sbmlsim/sensitivity/parameters.py

@@ -8,6 +8,7 @@
 from sbmlutils.console import console
 from sbmlutils.factory import ValueWithUnit
 
+import roadrunner
 
 
 @dataclass
@@ -27,6 +28,7 @@ def parameters_for_sensitivity_analysis(
     sbml_path: Path,
     exclude_ids: Optional[set[str]] = None,
     exclude_na: bool = True,
+    exclude_zero: bool = True,
 ) -> list[SensitivityParameter]:
     """Retrieve parameters from model for the sensitivity analysis.
 
@@ -50,7 +52,7 @@ def parameter_from_sbase(sbase: libsbml.SBase) -> SensitivityParameter:
 
         return parameter
 
-
+    r: roadrunner.RoadRunner = roadrunner.RoadRunner(str(sbml_path))
     doc: libsbml.SBMLDocument = libsbml.readSBMLFromFile(str(sbml_path))
     sbml_model: libsbml.Model = doc.getModel()
 
@@ -64,6 +66,8 @@ def parameter_from_sbase(sbase: libsbml.SBase) -> SensitivityParameter:
         if p.getConstant() is True:
             if exclude_na and np.isnan(p.getValue()):
                 exclude_ids.add(sid)
+            if exclude_zero and np.isclose(r.getValue(sid), 0.0):
+                exclude_ids.add(sid)
             parameters.append(parameter_from_sbase(p))
 
     # constant compartments
@@ -73,6 +77,8 @@ def parameter_from_sbase(sbase: libsbml.SBase) -> SensitivityParameter:
         if c.getConstant() is True:
             if exclude_na and np.isnan(c.getSize()):
                 exclude_ids.add(sid)
+            if exclude_zero and np.isclose(r.getValue(sid), 0.0):
+                exclude_ids.add(sid)
             parameters.append(parameter_from_sbase(c))
 
     # constant species or boundaryCondition == True
@@ -87,6 +93,11 @@ def parameter_from_sbase(sbase: libsbml.SBase) -> SensitivityParameter:
                 exclude_ids.add(sid)
             elif s.isSetInitialConcentration() and np.isnan(s.getInitialConcentration()):
                 exclude_ids.add(sid)
+        if exclude_zero:
+            if s.isSetInitialAmount() and np.isclose(s.getInitialAmount(), 0.0):
+                exclude_ids.add(sid)
+            elif s.isSetInitialConcentration() and np.isclose(s.getInitialConcentration(), 0.0):
+                exclude_ids.add(sid)
 
         if s.getConstant() is True or s.getBoundaryCondition() is True:
             parameters.append(parameter_from_sbase(s))

src/sbmlsim/sensitivity/plots.py

@@ -11,7 +11,14 @@
 import numpy as np
 import pandas as pd
 
-def heatmap(df: pd.DataFrame, cutoff: float=0.01, annotate_values=True, transpose: bool=False):
+def heatmap(
+    df: pd.DataFrame,
+    cutoff: float=0.01,
+    annotate_values=True,
+    cluster_rows: bool = True, # cluster parameters
+    cluster_cols: bool = False, # cluster outputs
+    transpose: bool=False
+):
     """Creates heatmap of model sensitivity"""
 
     def calculate_mask(df, cutoff=0.01):
@@ -24,7 +31,7 @@ def calculate_mask(df, cutoff=0.01):
                 mask[index] = False
         return pd.DataFrame(data=mask, columns=df.columns, index=df.index)
 
-    def calculate_subset(df, cutoff=0.01):
+    def calculate_subset(df, cutoff=0.01) -> pd.DataFrame:
         """Calculates subset of data frame consisting of rows where at least
         one value is above cutoff."""
         return df[(df.abs() >= cutoff).any(axis=1)]
@@ -40,23 +47,35 @@ def calculate_subset(df, cutoff=0.01):
     yticklabels = [pid for pid in df_subset.index]
     xticklabels = [pid for pid in df_subset.columns]
 
+    n_outputs = df_subset.shape[1]
+    n_parameters = df_subset.shape[0]
+    figsize = (7, int(n_parameters / n_outputs * 7)/2)
+
+    colorbar_range = 2.0
+
     # plot heatmap
     ax = sns.clustermap(
         df_subset,
         center=0,
-        vmin=-0.2,
-        vmax=0.2,
+        vmin=-colorbar_range,
+        vmax=colorbar_range,
         xticklabels=xticklabels,
         yticklabels=yticklabels,
         cmap="seismic",
-        cbar_pos=(0.05, 0.25, 0.03, 0.4),
+        # 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",
+            "label": "sensitivity"
+        },
         annot=annotate_values,
         fmt="1.2f",
-        annot_kws={"size": 13},
+        annot_kws={"size": 11},
         mask=df_subset_mask,
         col_cluster=False,
+        row_cluster=True,
         method="single",
-        figsize=(20, 20),
+        figsize=figsize,
     )
     plt.setp(
         ax.ax_heatmap.get_xticklabels(),