Creating paper plots and adding convenience methods

Author: joseph-hellerstein

notebooks/linearity_analysis.ipynb

@@ -0,0 +1,73 @@
+'''Creates the plots used in the paper.'''
+
+from src.biomodels_iterator import BiomodelsIterator
+import src.constants as cn
+from src.model import Model
+from src.score import Score
+from src.simulator import Simulator
+from src.system_discovery import SystemDiscovery, discoverNetwork
+from src.timecourse import Timecourse
+from src.timecourse_iterator import TimecourseIterator
+
+IS_PLOT = False
+
+import constants as cn
+if not IS_PLOT:
+    import matplotlib # type: ignore
+    matplotlib.use("PDF")  # Use non-interactive backend for testing
+import matplotlib.pyplot as plt # type: ignore
+import numpy as np  # type: ignore
+import os
+import pandas as pd # type: ignore
+from typing import Optional
+
+NUM_POINT = 1000
+
+############### Helper Functions####################
+
+def doPlot(model_num: int, poly_degree=1, threshold=0.001, species_names: Optional[list[str]] = None,
+        num_point:int=NUM_POINT):
+    start_time = 0
+    model = Model.makeBiomodel(model_num=model_num)
+    timecourse = Timecourse(model, start_time=start_time, num_point=num_point)
+    sdr = discoverNetwork(timecourse.timecourse_df, poly_degree=poly_degree, threshold=threshold,
+            plot_species_names=species_names, is_plot=IS_PLOT, subtitle=f"BioModel {model_num}",
+            is_plot_heatmap=False, is_print_equations=False, is_plot_comparisons=True, is_print_accuracy=False)
+    return sdr
+
+################################################
+# Linear Fits
+################################################
+############### Time course ####################
+sdr = doPlot(968, species_names=["SOCS1", "IL7IL7RJAK1"])
+sdr.fig.savefig(os.path.join(cn.PAPER_DIR, "linear_fit_968.pdf"), bbox_inches="tight", dpi=300) # type: ignore
+##
+sdr = doPlot(1004, species_names= ["IL6ext", "STAT3mRNA"])
+sdr.fig.savefig(os.path.join(cn.PAPER_DIR, "linear_fit_1004.pdf"), bbox_inches="tight", dpi=300) # type: ignore
+#
+############### CDFs ####################
+path = os.path.join(cn.DATA_DIR, "linear_predictor_scores-0.001.csv")
+score = Score.deserialize(path)
+fig = score.plotCDF(["min", "p10", "p50", "max"], xlabel= "model accuracy", is_plot_species=False,
+        is_plot_model=True, title=f"BioModel Models").fig
+if IS_PLOT:
+    plt.show()
+fig.savefig(os.path.join(cn.PAPER_DIR, "linear_fit_model_cdf.pdf"), bbox_inches="tight", dpi=300) # type: ignore
+#
+fig = score.plotCDF(["min", "p10", "p50", "max"], xlabel= "model accuracy", is_plot_species=True,
+        is_plot_model=False, title=f"BioModel Species").fig
+fig.savefig(os.path.join(cn.PAPER_DIR, "linear_fit_species_cdf.pdf"), bbox_inches="tight", dpi=300) # type: ignore
+if IS_PLOT:
+    plt.show()
+
+################################################
+# Perturbations
+################################################
+############### Time course ####################
+apr = SystemDiscovery.analyzePerturbations(1004, perturbations=[-50, -10, 0, 10, 50], frac_scatter_skip=0.05,
+        subtitle="BioModel 1004", plot_species_names=["IL6ext","IL6int"])
+apr.fig.savefig(os.path.join(cn.PAPER_DIR, "perturbation_fit_1004.pdf"), bbox_inches="tight", dpi=300) # type: ignore
+#
+apr = SystemDiscovery.analyzePerturbations(968, perturbations=[-50, -10, 0, 10, 50], frac_scatter_skip=0.05,
+        subtitle=f"BioModel 968", plot_species_names= ["SOCS1", "IL7IL7RJAK1"])
+apr.fig.savefig(os.path.join(cn.PAPER_DIR, "perturbation_fit_968.pdf"), bbox_inches="tight", dpi=300) # type: ignore

scripts/evaluate_monomial_models.py

@@ -1,14 +1,12 @@
 """Perturbation study: how do perturbed initial conditions affect SystemDiscovery R²?
 
-Runs only on models whose deg1_min >= MIN_R2 in evaluate_monomial_models-0.01.csv.
+Runs only on models with a explicitly specified end time.
 For each qualifying model, SystemDiscovery.analyzePerturbations is called with
 perturbation_value_fractions of -50%, -20%, -10%, -5%, 0%, +5%, +10%, +20%, +50%.
 R² is computed using the derivative method.
 
 Output CSV: data/perturbation_study.csv
 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
@@ -62,15 +60,16 @@ def main(is_initialize: bool = False) -> pd.DataFrame:
                 perturbation_species_fraction=SPECIES_FRACTION,
                 poly_degree=POLY_DEGREE,
                 is_plot=False,
-            )
+            ).df
         except Exception as exc:
             print(f"  [error] {item.model_name}: {exc}", file=sys.stderr)
             continue
 
+        analyze_df[cn.COL_THRESHOLD] = THRESHOLD
         current_df = pd.read_csv(OUTPUT_PATH) if os.path.isfile(OUTPUT_PATH) else pd.DataFrame()
         # Ensure analyze_df is a DataFrame before concatenating/writing.
         if isinstance(analyze_df, pd.Series):
-            analyze_df = pd.DataFrame([analyze_df.to_dict()])
+            analyze_df = pd.DataFrame([analyze_df.to_dict()]).df
         full_df = pd.concat([current_df, analyze_df], ignore_index=True) if len(current_df) > 0 else analyze_df
         full_df.to_csv(OUTPUT_PATH, index=False)
 

scripts/make_paper_plots.py

@@ -14,6 +14,7 @@
 TIMECOURSE_SERIALIZATION_DIR = os.path.join(SERIALIZATION_DIR,
         "timecourse")
 TIMECOURSE_ZIP_PATH = os.path.join(TIMECOURSE_SERIALIZATION_DIR, "timecourse.zip")
+PAPER_DIR = os.path.join(PROJECT_DIR, "paper")
 
 # Types
 TYPE_ROADRUNNER = "tellurium.roadrunner.extended_roadrunner.ExtendedRoadRunner"
@@ -43,6 +44,7 @@
 COL_NAMES = [COL_MODEL_NAME, COL_MAXCV, COL_ENDTIME, COL_ENDTIME_SOURCE]
 COL_PERTURBATION = "perturbation"  # Perturbation value fraction used in simulation
 COL_SYSTEM_ID = "system_id"  # Unique identifier for the system, e.g. model name or species name.
+COL_THRESHOLD = "threshold"  # Threshold used in SystemDiscovery
 COL_MEAN = "mean" # Mean value of the valid values used in calculating the statistics.
 COL_MIN = "min"  # Minimum value of the valid values used in calculating the statistics.
 COL_MAX = "max"  # Maximum value of the valid values used in calculating the statistics.

scripts/perturbation_study.py

@@ -271,25 +271,25 @@ def doPlot(value_arr: np.ndarray):
         plot_options = PlotOptions(**kwargs)
         legend = []
         if is_plot_model:
-            for metric_name in metric_names:
+            for mname in metric_names:
                 model_df = df[df[cn.COL_AGGREGATION_TYPE] == cn.COL_AGGREGATION_TYPE_MODEL]
-                if not model_df.empty and metric_name in model_df.columns:
-                    legend.append(f"{metric_name} (model)")
-                    value_arr = np.array(model_df[metric_name].values)
+                if not model_df.empty and mname in model_df.columns:
+                    legend.append(f"{mname} (model)")
+                    value_arr = np.array(model_df[mname].values)
                     doPlot(value_arr)
                     plotted_any = True
                 else:
-                    missing_metrics.append(f"model: {metric_name}")
+                    missing_metrics.append(f"model: {mname}")
         if is_plot_species:
-                for metric_name in metric_names:
+                for mname in metric_names:
                     species_df = df[df[cn.COL_AGGREGATION_TYPE] != cn.COL_AGGREGATION_TYPE_MODEL]
-                    if not species_df.empty and metric_name in species_df.columns:
-                        legend.append(f"{metric_name} (species)")
-                        value_arr = np.array(species_df[metric_name].values)
+                    if not species_df.empty and mname in species_df.columns:
+                        legend.append(f"{mname} (species)")
+                        value_arr = np.array(species_df[mname].values)
                         doPlot(value_arr)
                         plotted_any = True
                     else:
-                        missing_metrics.append(f"species: {metric_name}")
+                        missing_metrics.append(f"species: {mname}")
         if len(missing_metrics) > 0:
             warnings.warn(f"No data found for metrics {missing_metrics}; nothing was plotted.", UserWarning)
         if not plotted_any: