1+ """
2+ Creates ClusteredJacobianCollection objects for BioModels.
3+ """
4+ import src.constants as cn
5+ from trajectory import Trajectory # type: ignore
6+ from trajectory_collection import TrajectoryCollection # type: ignore
7+ from biomodels_iterator import BiomodelsIterator # type: ignore
8+
9+ import os
10+ import pandas as pd # type: ignore
11+ from typing import Dict, List, Optional, Tuple
12+
13+ import numpy as np # type: ignore
14+ import tellurium as te # type: ignore
15+ from src.l_roadrunner import LRoadrunner # type: ignore
16+
17+ OUTPUT_DATA_FILE = os.path.join(cn.DATA_DIR, "model_linearity_analysis_data.csv")
18+
19+
20+ class BiomodelsCluster:
21+ """Analyzes linearity of SBML or Antimony models by collecting Jacobians over time."""
22+
23+ def __init__(
24+ self,
25+ model_name: str,
26+ start_time: float = 0,
27+ end_time: float = np.nan,
28+ num_point: int = 100,
29+ diameter_metric: str = cn.DIAMETER_IVP,
30+ ) -> None:
31+ """
32+ Initialize a BiomodelsCluster with a model and simulation parameters.
33+
34+ Parameters
35+ ----------
36+ model_name : str
37+ Name of the model to analyze.
38+ start_time : float
39+ Simulation start time (default: 0).
40+ end_time : float, optional
41+ Simulation end time (default: 10).
42+ num_point : int
43+ Number of simulation timepoints (default: 100).
44+ """
45+ self.model_name = model_name
46+ self.start_time = start_time
47+ if np.isnan(end_time):
48+ end_time = LRoadrunner.endtime_dct.get(model_name, cn.END_TIME) # type: ignore
49+ self.end_time = end_time
50+ self.num_point = num_point
51+ self._diameter_metric = diameter_metric
52+ #
53+ self.sbml_str = self._getSbml()
54+ self.l_roadrunner = LRoadrunner(self.sbml_str, start_time=start_time,
55+ end_time=end_time, num_point=num_point)
56+ self._jacobian_collection = Trajectory(self.l_roadrunner,
57+ diameter_metric=diameter_metric)
58+
59+ @staticmethod
60+ def _report(text: str):
61+ """Print a report message if reporting is enabled."""
62+ print(text)
63+
64+ def _getSbml(self) -> str:
65+ """Search for an SBML file in the given directory and return its contents as a string."""
66+ dir_path = os.path.join(cn.BIOMODELS_DIR, self.model_name)
67+ sbml_files = os.listdir(dir_path)
68+ sbml_files = [f for f in sbml_files
69+ if f.endswith(".xml") and not "manifest" in f.lower()]
70+ if not sbml_files:
71+ raise FileNotFoundError("No SBML file found for model: " + self.model_name)
72+ with open(os.path.join(dir_path, sbml_files[0]), "r") as f:
73+ return f.read()
74+
75+ def cluster(self, n_cluster: int, is_sequential_partition: bool = True) -> TrajectoryCollection:
76+ """Partition the Jacobians into clusters and return a ClusteredJacobianCollection."""
77+ if is_sequential_partition:
78+ trajectory_collection = TrajectoryCollection(
79+ self._jacobian_collection.sequentialPartition(
80+ n_cluster=n_cluster))
81+ else:
82+ trajectory_collection = TrajectoryCollection(
83+ self._jacobian_collection.nonsequentialPartition(
84+ n_cluster=n_cluster))
85+ return trajectory_collection
86+
87+ @classmethod
88+ def clusterAnalysis(
89+ cls,
90+ directory: str = cn.BIOMODELS_DIR,
91+ output_data_file: str = OUTPUT_DATA_FILE,
92+ excluded_models: Optional[List[str]] = None,
93+ start_time: float = cn.START_TIME,
94+ num_point: int = cn.NUM_POINTS,
95+ n_cluster: int = 1,
96+ is_report: bool = True,
97+ is_sequential_partition: bool = True,
98+ diameter_metric: str = cn.DIAMETER_IVP,
99+ first_model_num: int = 0,
100+ last_model_num: int = int(1e9),
101+ ) -> pd.DataFrame:
102+ """
103+ For each model in BioModels, partition its Jacobians into n_cluster clusters and save
104+ the max CV of the clusters to a CSV.
105+ Two partitionation methods are available:
106+ k-means clustering (partitionJacobians) and sequential partitioning
107+
108+ Parameters
109+ ----------
110+ directory : str
111+ Path to the directory containing BioModel subdirectories. Defaults to cn.BIOMODELS_DIR.
112+ start_time : float
113+ The start time for the simulation.
114+ output_data_file : str
115+ Path to the CSV file where results will be saved.
116+ excluded_models : Optional[List[str]]
117+ List of model identifiers to exclude from processing.
118+ end_time : float
119+ The end time for the simulation. If NaN, it will be determined from LRoadrunner's endtime_dct or default to 10.
120+ num_point : int
121+ The number of time points to simulate.
122+ n_cluster : int
123+ Number of clusters of Jacobians for timepoints to use for k-means clustering.
124+ is_report : bool
125+ Whether to print report messages during processing.
126+ is_sequential_partition : bool
127+ Whether to use sequential partitioning instead of k-means clustering.
128+ diameter_metric : str
129+ The metric to use for calculating the diameter of each cluster.
130+ first_model_num : int
131+ The first model number to include (inclusive).
132+ last_model_num : int
133+ The last model number to include (inclusive).
134+
135+ Returns
136+ -------
137+ pd.DataFrame
138+ DataFrame containing the results for each model.
139+ """
140+ excluded_models = excluded_models if excluded_models is not None else []
141+ if excluded_models is None:
142+ excluded_models = []
143+ iterator = BiomodelsIterator(
144+ biomodels_dir=directory,
145+ excluded_models=excluded_models,
146+ existing_csv_path=output_data_file,
147+ is_report=is_report,
148+ first_model_num=first_model_num,
149+ last_model_num=last_model_num,
150+ )
151+ existing_df = iterator._existing_df
152+ ##
153+ def _write_csv(result_dct: Dict[str, float]) -> pd.DataFrame:
154+ """Write the given results to the output CSV, appending to existing data."""
155+ df = pd.DataFrame(result_dct)
156+ df = pd.concat([existing_df, df], ignore_index=False) if not existing_df.empty else df
157+ df.set_index(cn.COL_MODEL_NAME, inplace=True)
158+ df.to_csv(output_data_file, header=True, index=True)
159+ return df
160+ ##
161+ # Iterate over models and append results to CSV after each model is processed
162+ col_names = list(set(cn.COL_NAMES) - {cn.COL_ENDTIME_SOURCE})
163+ result_dct: dict = {c: [] for c in col_names}
164+ for item in iterator:
165+ try:
166+ biomodels_cluster = BiomodelsCluster(item.model_name,
167+ start_time=start_time, num_point=num_point,
168+ diameter_metric=diameter_metric)
169+ except Exception as e:
170+ print(e)
171+ continue
172+ cjc = biomodels_cluster.cluster(n_cluster=n_cluster, is_sequential_partition=is_sequential_partition)
173+ result_dct[cn.COL_MODEL_NAME].append(item.model_name)
174+ result_dct[cn.COL_MAXCV].append(cjc.max_cv)
175+ result_dct[cn.COL_ENDTIME].append(biomodels_cluster.end_time)
176+ result_df = _write_csv(result_dct)
177+ #
178+ result_df = _write_csv(result_dct)
179+ return result_df
0 commit comments