Skip to content

Commit 45bd126

Browse files
Better cli for perturbation_study.py
1 parent 6d274a8 commit 45bd126

10 files changed

Lines changed: 784 additions & 4343 deletions

data/perturbation_study.csv

Lines changed: 0 additions & 3691 deletions
This file was deleted.

docs/perturbation_study.md

Lines changed: 68 additions & 23 deletions
Original file line numberDiff line numberDiff line change
@@ -1,25 +1,70 @@
11
# Perturbation Study
22

3-
This study explores how perturbation initial values affects the quality of linear fits for models and for species.
4-
5-
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``.
6-
7-
* The arguments to ``analyzePerturbations``
8-
9-
* model: Model
10-
* training_df: pd.DataFrame
11-
* threshold: float
12-
* perturbations: list[float]
13-
* perturbation_species_fraction: float = 1.0
14-
* figsize: tuple[float, float] | None = None
15-
* poly_degree: int = 1
16-
* frac_keep: float = 0.2
17-
* is_plot: bool = True
18-
* ``analyzePerturbations`` returns a pandas.Series
19-
* @scripts/perturbation_study has a default path for CSV of ``perturbation_study.csv``
20-
* Calculate $R^2$ using the "derivative" method.
21-
* Revise and/or rewrite @perturbation_study.py.
22-
* As needed, make use of existing classes and functions in @src/. Use code in @scripts/ as a guide for implementation, but do not import these modules.
23-
* The $R_2$ value for a model is the minimum $R^2$ value for all species in the model.
24-
* The poly_degree is 1.
25-
* Do $R^2$ clamping so that $0 \leq R^2 \leq 1$.
3+
This study explores how perturbing initial species values affects the quality of linear fits for models and for individual species.
4+
5+
In `data/evaluate_monomial_models-0.01`, about 15% of the models have a `deg1_min` of at least 0.9 on their training data. This study compares those unperturbed results with perturbations applied to initial values. In each run, the training timecourse is simulated unperturbed. Separately, new timecourses are constructed where all species' initial values are changed by a signed fractional amount: ±5%, ±10%, ±20%, ±50% (plus an unperturbed 0% reference), with `perturbation_species_fraction=1.0`. For each perturbation level a separate ``Timecourse`` is simulated, and then the SINDy model trained on the unperturbed training data is used to predict each perturbed timecourse. Use the class method ``SystemDiscovery.analyzePerturbations`` (formerly ``perturbationAnalysis``) for this analysis. The per-model results are written as CSV files in `data/`.
6+
7+
## Parameters of ``analyzePerturbations``
8+
9+
| Parameter | Type | Default | Description |
10+
|---|---|---|---|
11+
| ``model`` | Model \| int || Simulates ground-truth timecourses for each perturbation. Accepts a model object or an integer BioModel number. |
12+
| ``training_df`` | pd.DataFrame | NULL_DF | Unperturbed timecourse used to fit the SINDy model. |
13+
| ``threshold`` | float | 0.001 | STLSQ sparsity threshold. |
14+
| ``poly_degree`` | int | 1 | Degree of the polynomial library (default linear). |
15+
| ``perturbations`` | list[float] | [-0.5, -0.2, -0.1, 0.0, 0.1, 0.2, 0.5] | Signed fractional perturbation values applied to initial species concentrations. |
16+
| ``col_percentile`` | str | "p10" | Column of the accuracy DataFrame used for per-perturbation summaries in plots. |
17+
| ``perturbation_species_fraction`` | float | 1.0 | Fraction of non-zero species whose initial values are perturbed when simulating each perturbed timecourse. |
18+
| ``frac_scatter_skip`` | float | 0.2 | Scatter-plot density: `num_skip = max(1, int(n_points * frac_scatter_skip))`. |
19+
| ``figsize`` | tuple[float, float] \| None | None | Figure size in inches for the trajectory comparison plot; auto-sized when None. |
20+
| ``plot_species_names`` | list[str] \| None | None | Species to include in the per-perturbation trajectory figure. Defaults to all species. |
21+
| ``subtitle`` | str | "Perturbation Analysis" | Title rendered on the output figure. |
22+
| ``is_analyze_model`` | bool | True | Include model-level rows (`aggregation_type='model'`) in the returned DataFrame. |
23+
| ``is_analyze_species`` | bool | True | Include per-species rows (one row per species name) in the returned DataFrame. |
24+
| ``is_plot`` | bool | True | Whether to show a trajectory comparison figure alongside the accuracy metrics. |
25+
26+
## What is calculated
27+
28+
For each perturbed timecourse, the fitted SINDy model predicts concentrations at every timepoint. Accuracy is computed pointwise as:
29+
30+
```
31+
accuracy = max(0, 1 - abs(prediction - actual) / abs(actual))
32+
```
33+
34+
Values are clipped to `[0, 1]`. Timepoints where `actual` is zero or non-finite receive a sentinel accuracy of `-1` and are excluded from aggregation. For each species × timepoint pair the Accuracy score aggregates across time via `StatisticCalculator`, producing percentile columns (`mean`, `min`, `max`, `count`, `invalid_count`, `p05`, `p10`, `p20`, `p25`, `p30`, `p50`, `p80`, `p90`, `p95`, `p99`).
35+
36+
## Output structure
37+
38+
``analyzePerturbations`` returns an ``AnalyzePerturbationsResult`` named tuple containing:
39+
40+
- ``.df`` — a DataFrame with one row per perturbation value at each requested aggregation level.
41+
- ``.fig`` — the trajectory comparison figure (or `None` if `is_plot=False`).
42+
43+
### Rows in the DataFrame
44+
45+
| Column | Source | Notes |
46+
|---|---|---|
47+
| ``perturbation`` | set from input list | The signed fractional perturbation value for this row; preserved exactly per-row. |
48+
| ``fraction_species_perturbable`` | `perturbation_species_fraction` argument | Same for every row in the result. |
49+
| ``system_id`` | `model.model_name` | Model name, repeated on each row. |
50+
| ``aggregation_type`` | computed | `'model'` for aggregated rows; species name string for per-species rows. |
51+
| ``mean``, ``min``, ``max``, ``count``, ``invalid_count``, ``p05````p99`` | `StatisticCalculator` on the Accuracy score | Aggregated across valid (non-sentinel) timepoints within that aggregation level. |
52+
53+
### Model-level vs species-level rows
54+
55+
- With **both** flags True (the default), the result contains both:
56+
- One model-level row per perturbation (`aggregation_type='model'`); statistics are aggregated across all non-zero species at each timepoint.
57+
- One per-species row per perturbation (`aggregation_type=<species name>`).
58+
- Set `is_analyze_model=False` or `is_analyze_species=False` to drop the corresponding level from the output entirely.
59+
60+
### CSV outputs in ``scripts/perturbation_study.py``
61+
62+
The script writes three variant CSVs, chosen by which aggregation levels are requested:
63+
64+
| Mode | Path pattern |
65+
|---|---|
66+
| model + species (default) | `data/perturbation_study-model_species{THRESHOLD}.csv` |
67+
| model only | `data/perturbation_study-model{THRESHOLD}.csv` |
68+
| species only | `data/perturbation_study-species{THRESHOLD}.csv` |
69+
70+
Each row in a CSV file is one perturbation × aggregation-level combination (i.e. not one row per model). The script skips models listed in its `EXCLUDES` set and resumes from an existing CSV on subsequent runs by checking `system_id`.

notebooks/perturbation_analysis.ipynb

Lines changed: 176 additions & 255 deletions
Large diffs are not rendered by default.

scripts/make_paper_plots.py

Lines changed: 15 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -23,6 +23,12 @@
2323

2424
NUM_POINT = 1000
2525

26+
################ Data #####################
27+
path = os.path.join(cn.DATA_DIR, "perturbation_study-0.001.csv")
28+
df = pd.read_csv(path)
29+
DF_P_DCT = {p: df[p] for p in ['min', 'p10', 'p50']}
30+
#
31+
2632
############### Helper Functions####################
2733

2834
def doPlot(model_num: int, poly_degree=1, threshold=0.001, species_names: Optional[list[str]] = None,
@@ -35,6 +41,13 @@ def doPlot(model_num: int, poly_degree=1, threshold=0.001, species_names: Option
3541
is_plot_heatmap=False, is_print_equations=False, is_plot_comparisons=True, is_print_accuracy=False)
3642
return sdr
3743

44+
def plotPerturbation(col: str):
45+
dff = DF_P_DCT[col]
46+
dff["aggregation_type"] = "model"
47+
score = Score.deserialize()
48+
score.score_df = dff
49+
score.plotCDF([-0.5, -0.1, 0, 0.10, 0.50], title="CDF: " + col)
50+
3851
################################################
3952
# Linear Fits
4053
################################################
@@ -70,4 +83,5 @@ def doPlot(model_num: int, poly_degree=1, threshold=0.001, species_names: Option
7083
#
7184
apr = SystemDiscovery.analyzePerturbations(968, perturbations=[-50, -10, 0, 10, 50], frac_scatter_skip=0.05,
7285
subtitle=f"BioModel 968", plot_species_names= ["SOCS1", "IL7IL7RJAK1"])
73-
apr.fig.savefig(os.path.join(cn.PAPER_DIR, "perturbation_fit_968.pdf"), bbox_inches="tight", dpi=300) # type: ignore
86+
apr.fig.savefig(os.path.join(cn.PAPER_DIR, "perturbation_fit_968.pdf"), bbox_inches="tight", dpi=300) # type: ignore
87+

scripts/perturbation_study.py

Lines changed: 106 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -9,6 +9,7 @@
99
Columns: model_name, threshold, and for each perturbation level
1010
"""
1111

12+
import argparse
1213
import os
1314
import sys
1415

@@ -18,24 +19,107 @@
1819
from src.system_discovery import SystemDiscovery
1920
from src.timecourse_iterator import TimecourseIterator
2021

21-
THRESHOLD = 0.001
22+
DEFAULT_THRESHOLD = 0.001
2223
POLY_DEGREE = 1
2324
SPECIES_FRACTION = 1.0
2425
PERTURBATIONS: list[float] = [-0.50, -0.20, -0.10, -0.05, 0.00, 0.05, 0.10, 0.20, 0.50]
2526

26-
OUTPUT_PATH = os.path.join(cn.DATA_DIR, f"perturbation_study-{THRESHOLD}.csv")
2727
MIN_R2 = 0.8
2828
COL_DEG1_MODEL_MAX = "deg1_max"
2929

3030
EXCLUDES = [
31-
"BIOMD0000000718",
31+
"BIOMD0000000338", # Got many errors
32+
"BIOMD0000000339", # Got many errors
33+
"BIOMD0000000378", # Got many errors
34+
"BIOMD0000000531", # Got many errors
35+
"BIOMD0000000532", # Got many errors
36+
"BIOMD0000000555", # Got many errors
37+
"BIOMD0000000559", # Got many errors
38+
"BIOMD0000000561", # Got many errors
39+
"BIOMD0000000570", # Got many errors
40+
"BIOMD0000000572", # Got many errors
41+
"BIOMD0000000627", # Got many errors
42+
"BIOMD0000000673", # Got many errors
43+
"BIOMD0000000711", # Got many errors
44+
"BIOMD0000000718", # Got many errors
45+
"BIOMD0000000721", # Got many errors
46+
"BIOMD0000000734", # Got many errors
47+
"BIOMD0000000763", # Got many errors
48+
"BIOMD0000000787", # Got many errors
49+
"BIOMD0000000809", # Got many errors
50+
"BIOMD0000000810", # Got many errors
51+
"BIOMD0000000834", # Got many errors
52+
"BIOMD0000000856", # Got many errors
53+
"BIOMD0000000864", # Got many errors
54+
"BIOMD0000000876", # Got many errors
55+
"BIOMD0000000879", # Got many errors
56+
"BIOMD0000000923", # Got many errors
57+
"BIOMD0000000943", # Got many errors
58+
"BIOMD0000000961", # Got many errors
59+
"BIOMD0000000972", # Got many errors
60+
"BIOMD0000000989", # Got many errors
61+
"BIOMD0000000990", # Got many errors
62+
"BIOMD0000001019", # Got many errors
63+
"BIOMD0000001020", # Got many errors
64+
"BIOMD0000001027", # Got many errors
3265
]
3366

3467

35-
def main(is_initialize: bool = False) -> pd.DataFrame:
36-
if not is_initialize and os.path.isfile(OUTPUT_PATH):
37-
print(f"Loading existing results from {OUTPUT_PATH}...")
38-
initial_df = pd.read_csv(OUTPUT_PATH)
68+
def _build_parser() -> argparse.ArgumentParser:
69+
parser = argparse.ArgumentParser(
70+
description="Perturbation study: how do perturbed initial conditions affect SystemDiscovery R²?",
71+
)
72+
parser.add_argument(
73+
"--threshold",
74+
type=float,
75+
default=DEFAULT_THRESHOLD,
76+
help=f"Feature selection threshold (default: {DEFAULT_THRESHOLD})",
77+
)
78+
parser.add_argument(
79+
"--is-analyze-model",
80+
dest="is_analyze_model",
81+
action="store_true",
82+
default=True,
83+
help="Analyze model-level results (default: True)",
84+
)
85+
parser.add_argument(
86+
"--no-is-analyze-model",
87+
dest="is_analyze_model",
88+
action="store_false",
89+
help="Disable model-level analysis",
90+
)
91+
parser.add_argument(
92+
"--is-analyze-species",
93+
dest="is_analyze_species",
94+
action="store_true",
95+
default=True,
96+
help="Analyze species-level results (default: True)",
97+
)
98+
parser.add_argument(
99+
"--no-is-analyze-species",
100+
dest="is_analyze_species",
101+
action="store_false",
102+
help="Disable species-level analysis",
103+
)
104+
return parser
105+
106+
107+
def main(is_initialize: bool = False, is_analyze_model: bool = True,
108+
is_analyze_species: bool = True, threshold: float = DEFAULT_THRESHOLD) -> pd.DataFrame:
109+
if is_analyze_model and is_analyze_species:
110+
print("Analyzing both model-level and species-level results...")
111+
output_path = os.path.join(cn.DATA_DIR, f"perturbation_study-model_species{threshold}.csv")
112+
elif is_analyze_model and not is_analyze_species:
113+
print("Analyzing model-level results only...")
114+
output_path = os.path.join(cn.DATA_DIR, f"perturbation_study-model{threshold}.csv")
115+
elif not is_analyze_model and is_analyze_species:
116+
print("Analyzing species-level results only...")
117+
output_path = os.path.join(cn.DATA_DIR, f"perturbation_study-species{threshold}.csv")
118+
else:
119+
raise ValueError("At least one of is_analyze_model or is_analyze_species must be True.")
120+
if not is_initialize and os.path.isfile(output_path):
121+
print(f"Loading existing results from {output_path}...")
122+
initial_df = pd.read_csv(output_path)
39123
else:
40124
initial_df = pd.DataFrame()
41125

@@ -55,31 +139,38 @@ def main(is_initialize: bool = False) -> pd.DataFrame:
55139
analyze_df = SystemDiscovery.analyzePerturbations(
56140
model=item.timecourse.model,
57141
training_df=item.timecourse.timecourse_df,
58-
threshold=THRESHOLD,
142+
threshold=threshold,
59143
perturbations=PERTURBATIONS,
60144
perturbation_species_fraction=SPECIES_FRACTION,
61145
poly_degree=POLY_DEGREE,
62146
is_plot=False,
147+
is_analyze_model=is_analyze_model,
148+
is_analyze_species=is_analyze_species,
63149
).df
64150
except Exception as exc:
65151
print(f" [error] {item.model_name}: {exc}", file=sys.stderr)
66152
continue
67153

68-
analyze_df[cn.COL_THRESHOLD] = THRESHOLD
69-
current_df = pd.read_csv(OUTPUT_PATH) if os.path.isfile(OUTPUT_PATH) else pd.DataFrame()
154+
analyze_df[cn.COL_THRESHOLD] = threshold
155+
current_df = pd.read_csv(output_path) if os.path.isfile(output_path) else pd.DataFrame()
70156
# Ensure analyze_df is a DataFrame before concatenating/writing.
71157
if isinstance(analyze_df, pd.Series):
72158
analyze_df = pd.DataFrame([analyze_df.to_dict()]).df
73159
full_df = pd.concat([current_df, analyze_df], ignore_index=True) if len(current_df) > 0 else analyze_df
74-
full_df.to_csv(OUTPUT_PATH, index=False)
160+
full_df.to_csv(output_path, index=False)
75161

76-
if os.path.isfile(OUTPUT_PATH):
77-
full_df = pd.read_csv(OUTPUT_PATH)
162+
if os.path.isfile(output_path):
163+
full_df = pd.read_csv(output_path)
78164
else:
79165
full_df = initial_df
80-
print(f"\nDone. {len(full_df)} rows in {OUTPUT_PATH}")
166+
print(f"\nDone. {len(full_df)} rows in {output_path}")
81167
return full_df
82168

83169

84170
if __name__ == "__main__":
85-
main()
171+
args = _build_parser().parse_args()
172+
main(
173+
is_analyze_model=args.is_analyze_model,
174+
is_analyze_species=args.is_analyze_species,
175+
threshold=args.threshold,
176+
)

src/score.py

Lines changed: 6 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -14,7 +14,7 @@
1414
import matplotlib.pyplot as plt # type: ignore
1515
import numpy as np # type: ignore
1616
import pandas as pd # type: ignore
17-
from typing import List # type: ignore
17+
from typing import List, Optional # type: ignore
1818
import warnings # type: ignore
1919

2020

@@ -25,7 +25,7 @@ class Score:
2525
"""
2626
SERIALIZATION_PATH = "score.csv" # Default path for persistence, can be overridden in subclasses.
2727

28-
def __init__(self, serialization_path: str = "", is_persist: bool = True,
28+
def __init__(self, serialization_path: Optional[str] = None, is_persist: bool = True,
2929
col_percentile: str = cn.COL_P10) -> None:
3030
"""
3131
Parameters
@@ -38,15 +38,15 @@ def __init__(self, serialization_path: str = "", is_persist: bool = True,
3838
The percentile to compute (e.g., 'p10', 'p50', 'p90'). Default is 'p10'.
3939
"""
4040
self._is_persist = is_persist
41-
if len(serialization_path) == 0:
41+
if serialization_path is None:
4242
serialization_path = self.SERIALIZATION_PATH
4343
self._serialization_path = serialization_path
4444
self._col_percentile = col_percentile
4545
#
4646
self.score_df = pd.DataFrame()
4747

4848
@classmethod
49-
def deserialize(cls, serialization_path: str) -> 'Score':
49+
def deserialize(cls, serialization_path: Optional[str]=None) -> 'Score':
5050
"""Loads a previously serialized score DataFrame from CSV.
5151
Uses a defulat of cn.COL_P10 for the percentile column, but this can be changed later if needed.
5252
@@ -61,7 +61,8 @@ def deserialize(cls, serialization_path: str) -> 'Score':
6161
The deserialized score object.
6262
"""
6363
score = cls(serialization_path=serialization_path, is_persist=False)
64-
score.score_df = pd.read_csv(serialization_path)
64+
if serialization_path is not None:
65+
score.score_df = pd.read_csv(serialization_path)
6566
score._col_percentile = cn.COL_P10 # Default to p10; can be changed later if needed.
6667
return score
6768

0 commit comments

Comments
 (0)