Serialization study

Author: joseph-hellerstein

data/perturbation_study.csv

@@ -3,7 +3,9 @@
 This study explores how perturbation initial values affects the quality of linear fits for models and for species.
 
 In @data/evaluate_monomial_models-0.01, we see that about 15% of the models have a deg1_min of at least 0.9 on their training data. This study compares these results with perturbations of initial values. In the studies, the training data are unperturbed time courses. Separately, time courses are constructed where initial values have been changed by $\pm 5\%$, $\pm 10\%$, $\pm 20\%$, $\pm 50\%$ with a perturbation species fraction of 1.0. A separate ``Timecourse`` is constructed for each perturbation (at total of 8 plus one for 0%), and then the model trained on unperturbed data is used to predict perturbed data. Please use and/or add capabilities to the class method ``SystemDiscovery.analyzePerturbations`` (formerly ``perturbationAnalysis``) to do this analysis. The result is a CSV file with the columns: model_name, threshold, r2_0, r2_-05, r2_-10, r2_-20, r2_-50, r2_+05, r2_+10, r2_+20, r2_+50. The output file path is new argument to ``analyzePerturbations``.
+
 * The arguments to ``analyzePerturbations``
+
   * model: Model
   * training_df: pd.DataFrame
   * threshold: float

docs/perturbation_study.md

@@ -1,23 +1,30 @@
 # Research Agenda
 
 ## How well can BioModels be modelled by a system of linear differential equaitons?
+
 - [x] Fit first order monomial to all models and calculate $R^2$ for species time courses.
-    - Minimum $R^2$ is model $R^2$
-    - Individual $R^2$ are for species.
-- [x] Evaluate the density of coefficients in the Jacobian.
+
+  - Minimum $R^2$ is model $R^2$
+  - Individual $R^2$ are for species.
+
+- [ ] Evaluate the density of coefficients in the Jacobian.
+  - [ ] Distributions don't look correct.
 - [ ] For the linear models, assess their dimensionality to see if dimension reduction is possible
 - [ ] Characterize the linear models based on the what is being model and possibly other characteristics.
 
   ## How robust is linearity to perturbations of initial values?
+
 - [ ] repeat the linear studies with perturbations of $\pm 5\%$, $\pm 19\%$, $\pm 20\%$, and $\pm 50\%$.
 - [ ] Can robustness be improved by training the regression on perturbation data?
 
 ## What are the main reasons for nonlinear behavior?
+
 - [ ] Analyze the nonlinear models to determine which species are nonlinear and how/when the Jacobian changes to look at reactions.
 
 ## Are some nonlinear models piecewise linear?
+
 - [ ] Use a standard package for partitioning regressions to see if linearity can be achieved in segements of the time course.
   - [ ] Normalize Jacobians
   - [ ] k-means cluster with cluster distances with minimum cluster size.
-  - [ ] SystemDiscovery for each cluster 
-  - [ ] Prediction using: (a) ${\bf x} (t)$ using SystemDiscover of cluster at $t$ similarly for adjacent time points; (b) apply gaussian kernel on points; (c) denormalize
+  - [ ] SystemDiscovery for each cluster
+  - [ ] Prediction using: (a) ${\bf x} (t)$ using SystemDiscover of cluster at $t$ similarly for adjacent time points; (b) apply gaussian kernel on points; (c) denormalize

docs/research_agenda.md

@@ -6,14 +6,15 @@
 R² is computed using the derivative method.
 
 Output CSV: data/perturbation_study.csv
-Columns: model_name, threshold, r2_-50, r2_-20, r2_-10, r2_-05, r2_0,
-         r2_+05, r2_+10, r2_+20, r2_+50
+Columns: model_name, threshold, and for each perturbation level
+         (r2_-50, r2_-20, r2_-10, r2_-05, r2_0, r2_+05, r2_+10, r2_+20, r2_+50):
+         three columns _min, _med, _max (clamped derivative R² across species).
 """
 
 import os
 import sys
 
-import pandas as pd
+import pandas as pd  # type: ignore
 
 import src.constants as cn
 from src.system_discovery import SystemDiscovery
@@ -26,15 +27,20 @@
 
 SOURCE_PATH = os.path.join(cn.DATA_DIR, "evaluate_monomial_models-0.01.csv")
 OUTPUT_PATH = os.path.join(cn.DATA_DIR, "perturbation_study.csv")
-MIN_DEG1 = 0.9
+MIN_R2 = 0.8
+COL_DEG1_MODEL_MAX = "deg1_max"
+
+EXCLUDES = [
+            "BIOMD0000000718",
+            ]
 
 
 def main(is_initialize: bool = False) -> pd.DataFrame:
     source_df = pd.read_csv(SOURCE_PATH)
     model_names: set[str] = set(
-        source_df.loc[source_df["deg1_min"] >= MIN_DEG1, cn.COL_MODEL_NAME]
+        source_df.loc[source_df[COL_DEG1_MODEL_MAX] >= MIN_R2, cn.COL_MODEL_NAME]
     )
-    print(f"Models with deg1_min >= {MIN_DEG1}: {len(model_names)}")
+    print(f"Models ({len(source_df)}) with deg1_min >= {MIN_R2}: {len(model_names)}")
 
     if not is_initialize and os.path.isfile(OUTPUT_PATH):
         print(f"Loading existing results from {OUTPUT_PATH}...")
@@ -53,6 +59,9 @@ def main(is_initialize: bool = False) -> pd.DataFrame:
         if item.model_name in already_done:
             print(f"Skipping {item.model_name} (already processed)", flush=True)
             continue
+        if item.model_name in EXCLUDES:
+            print(f"Skipping {item.model_name} (excluded)", flush=True)
+            continue
         print(f"Processing {item.model_name}...", flush=True)
         try:
             series = SystemDiscovery.analyzePerturbations(

notebooks/linearity_analysis.ipynb

@@ -48,6 +48,7 @@
 import pysindy as ps # type: ignore
 from pysindy.feature_library import PolynomialLibrary # type: ignore
 from scipy.integrate import solve_ivp # type: ignore
+import sys
 from typing import Literal, cast
 import warnings
 
@@ -170,7 +171,6 @@ def __init__(
         self.species_cols = species_cols
         self._X_list: list[np.ndarray] = [d[species_cols].to_numpy(dtype=float) for d in dfs]
         self._time_list: list[np.ndarray] = [d.index.to_numpy(dtype=float) for d in dfs]
-        # First trajectory used for simulation and plotting
         self.time_arr: np.ndarray = self._time_list[0]
         self.X: np.ndarray = self._X_list[0]
         #
@@ -301,17 +301,17 @@ def analyzePerturbations(
         poly_degree: int = 1,
         frac_keep: float = 0.2,
         is_plot: bool = True,
-    ) -> "pd.Series":
+    ) -> pd.Series:
         """Fit on training_df and evaluate derivative R² at each perturbation level.
 
         For each value in *perturbations* a fresh Timecourse is simulated from
         *model* with that ``perturbation_value_fraction``.  The SINDy model
         (fitted on the unperturbed *training_df*) is evaluated against each
         perturbed timecourse using the derivative R² method.  The reported R²
-        per perturbation is the minimum clamped R² across all species.
+        per perturbation are the min, median, and max clamped R² across all species.
 
         When *is_plot* is True, a trajectory comparison figure is also shown
-        (simulation R² shown in the legend for visual context).
+        (derivative) R² shown in the legend for visual context).
 
         Parameters
         ----------
@@ -337,9 +337,9 @@ def analyzePerturbations(
         Returns
         -------
         pd.Series
-            Index: ``model_name``, ``threshold``, and one ``r2_*`` key per
-            perturbation value.  R² values are the minimum clamped derivative
-            R² across species, in [0, 1].
+            Index: ``model_name``, ``threshold``, and three keys per
+            perturbation value (``r2_*_min``, ``r2_*_med``, ``r2_*_max``).
+            R² values are clamped derivative R² across species, in [0, 1].
         """
         import src.constants as cn  # avoid circular at module level
 
@@ -368,19 +368,22 @@ def analyzePerturbations(
                 )
                 test_df = tc.timecourse_df
                 r2_dict = disc.calculateRsq(method="derivative", test_df=test_df)
-                result[col_name] = cls._normalize_rsq(
-                    float(np.min(list(r2_dict.values())))
-                )
+                r2_clamped = {k: cls._normalize_rsq(v) for k, v in r2_dict.items()}
+                vals = list(r2_clamped.values())
+                result[f"{col_name}_min"] = float(np.min(vals))
+                result[f"{col_name}_med"] = float(np.median(vals))
+                result[f"{col_name}_max"] = float(np.max(vals))
                 if is_plot:
                     try:
                         pred_df: pd.DataFrame | None = disc.predict(test_df)
                     except Exception:
                         pred_df = None
-                    r2_sim = disc.calculateRsq(method="simulation", test_df=test_df)
-                    plot_records.append((p, test_df, pred_df, r2_sim))
+                    plot_records.append((p, test_df, pred_df, r2_clamped))
             except Exception as exc:
                 print(f"  [p={p}] {model.model_name}: {exc}", file=sys.stderr)
-                result[col_name] = float("nan")
+                result[f"{col_name}_min"] = float("nan")
+                result[f"{col_name}_med"] = float("nan")
+                result[f"{col_name}_max"] = float("nan")
 
         if is_plot and plot_records:
             n = len(disc.species_names)
@@ -532,7 +535,7 @@ def plotResult(
             pred_df = None
             prediction_ok = False
 
-        r2_vals = self.calculateRsq(method="simulation", test_df=test_df)
+        r2_vals = self.calculateRsq(method="derivative", test_df=test_df)
 
         num_skip_point = max(1, len(time_arr) // num_true_point)
         for idx, name in enumerate(self.species_names):
@@ -545,7 +548,7 @@ def plotResult(
             r2 = r2_vals.get(name, float("nan"))
             title = f"{name}"
             if not np.isnan(r2):
-                title += f"   R²={r2:.4f}"  # may be negative; clamping removed
+                title += f"   R²={r2:.4f}"
             low_y = X[:, idx].min()
             high_y = X[:, idx].max()
             if np.isclose(low_y, high_y):
@@ -678,17 +681,16 @@ def minR2(self, test_df: pd.DataFrame = NULL_DF) -> float:
         Parameters
         ----------
         test_df : pd.DataFrame, optional
-            If provided, R² is evaluated against this DataFrame (simulation
-            from its initial conditions vs observed).  When omitted, the
-            training data are used.
+            If provided, R² is evaluated against this DataFrame (derivative)
+            When omitted, the training data are used.
 
         Returns
         -------
         float
             Minimum R² across species, clamped to [0, 1].
         """
         self._require_fitted()
-        r2_raw = self.calculateRsq(method="simulation", test_df=test_df)
+        r2_raw = self.calculateRsq(method="derivative", test_df=test_df)
         return self._normalize_rsq(float(np.min(list(r2_raw.values()))))
 
     def score(self) -> ScoreInfo:
@@ -967,7 +969,7 @@ def discoverNetwork(
     print()
 
     try:
-        r2_sim = disc.calculateRsq(method="simulation", test_df=test_df)
+        r2_sim = disc.calculateRsq(method="derivative", test_df=test_df)
         print("R² on simulated trajectories per species:")
         for name, val in r2_sim.items():
             print(f"  {name}: {val:.6f}")

scripts/perturbation_study.py

@@ -16,7 +16,7 @@
 sys.path.insert(0, os.path.join(cn.PROJECT_DIR, "scripts"))
 
 from perturbation_study import (  # type: ignore  # noqa: E402
-    MIN_DEG1,
+    MIN_R2,
     PERTURBATIONS,
     THRESHOLD,
     main,