First test complete implementation of TrajectoryCollection.

Author: joseph-hellerstein

scripts/calculate_linear_prediction_scores.py

@@ -130,7 +130,7 @@ def processModels(first_model_num: int, last_model_num: int, process_index: int,
                 # Handle large Jacobians
                 trajectory = Trajectory.makeBiomodel(model_name=model_name)
                 prediction_df = trajectory.predict(is_adjust_fitted_jacobian=True)
-            score.addTestResult(trajectory.timecourse, prediction_df, description=model_name)
+            score.addTestResult(trajectory.timecourse_df, prediction_df, description=model_name)
         except Exception as e:
             print(f"Error occurred while processing model {model_name}: {e}")
             continue

src/l_roadrunner.py

@@ -92,10 +92,10 @@ def __init__(self, roadrunner_specification: str,
         self._start_time = start_time
         self.num_point = num_point
         self._end_time: float = end_time if end_time is not None else np.nan
-        self._sedml_str = sedml_str
+        self.sedml_str = sedml_str
         self.end_time_source: Optional[str] = None
         self._species_names: List[str] = []
-        self._timecourse = pd.DataFrame()
+        self._timecourse_df = pd.DataFrame()
         self.model_name = model_name
 
     @classmethod
@@ -176,7 +176,7 @@ def start_time(self, start_time: float) -> None:
         self._start_time = start_time
 
     @property
-    def timecourse(self) -> pd.DataFrame:
+    def timecourse_df(self) -> pd.DataFrame:
         """
         Simulate the model and return the time course as a DataFrame.
 
@@ -186,12 +186,12 @@ def timecourse(self) -> pd.DataFrame:
             A DataFrame containing the time course of the simulation, with columns for time and each floating species.
             The index are the time values.
         """
-        if self._timecourse.empty:
+        if self._timecourse_df.empty:
             result_arr = self.simulate(is_with_timepoints=True)
             df = pd.DataFrame(result_arr, columns=["time"] + self.species_names)
             df = df.set_index("time")
-            self._timecourse = df
-        return self._timecourse
+            self._timecourse_df = df
+        return self._timecourse_df
 
     def getInitialValues(self) -> np.ndarray:
         """
@@ -208,7 +208,6 @@ def getForcingInputs(self) -> np.ndarray:
         This is used in the linear predictor to extrapolate from the initial state.
         """
         rr = self.getRoadrunner()
-        #_ = rr.simulate(self.start_time, self.end_time, 2)
         _ = rr.simulate(self.end_time, self.end_time*1.001, 2)
         jacobian_arr = np.array(rr.getFullJacobian())
         f_arr = np.array(rr.getRatesOfChange())
@@ -403,17 +402,17 @@ def _calculateEndtimeSBML(self) -> float:
         float or np.nan
             The end time specified in the SBML string, or np.nan if no valid specification is found.
         """
-        if not self._sedml_str:
+        if not self.sedml_str:
             self._msg("No SED-ML string provided, skipping SED-ML end time extraction.")
             return self._end_time
         #
-        end_time_match = re.search(r'outputEndTime\s*=\s*"([0-9.]+)"', self._sedml_str)
+        end_time_match = re.search(r'outputEndTime\s*=\s*"([0-9.]+)"', self.sedml_str)
         if end_time_match:
             end_time_match = re.search(r'outputEndTime\s*=\s*"([0-9.]+)"',
-                    self._sedml_str)
+                    self.sedml_str)
             if end_time_match:
                 end_time = float(end_time_match.group(1))
-                if not DEFAULT_END_TIME_STR in self._sedml_str and end_time > 0:
+                if not DEFAULT_END_TIME_STR in self.sedml_str and end_time > 0:
                     self._end_time = end_time
                 elif end_time != SBML_DEFAULT_END_TIME:
                     self._end_time = end_time
@@ -423,13 +422,6 @@ def _calculateEndtimeSBML(self) -> float:
                 self._msg("No valid outputEndTime attribute found in SED-ML string.")
         else:
             self._msg("No outputEndTime attribute found in SED-ML string.")
-        """ # Validate the end time by simulating to it and checking for errors.
-        if self._end_time is not np.nan:
-            try:
-                rr = self.getRoadrunner()
-                rr.simulate(self.start_time, self._end_time, 2)
-            except Exception:
-                self._end_time = np.nan """
         #
         return self._end_time
 
@@ -460,8 +452,11 @@ def simulate(self, is_with_timepoints: bool=False) -> np.ndarray:
         idx = 1
         if is_with_timepoints:
             idx = 0
+        rr = self.getRoadrunner()
+        if self.start_time > 0:
+            result_arr = rr.simulate(0, self.start_time, 2)
         try:
-            result_arr = self.getRoadrunner().simulate(self.start_time, self.end_time, self.num_point)
+            result_arr = rr.simulate(self.start_time, self.end_time, self.num_point)
         except Exception as e:
             print(f"Error occurred while simulating: {e}")
             return np.array([]).reshape(0, 0)
@@ -485,25 +480,27 @@ def makeJacobians(self)->Tuple[np.ndarray, np.ndarray]:
         ValueError
             If the model has no floating species after reset.
         """
-        rr = self.getRoadrunner()
-        if len(rr.getFloatingSpeciesIds()) == 0:
+        if len(self.species_names) == 0:
             raise ValueError("Model has no floating species; cannot compute Jacobian.")
         try:
-            result_arr = rr.simulate(self.start_time, self.end_time, self.num_point)
+            result_arr = self.simulate(is_with_timepoints=True)
         except RuntimeError as e:
             raise ValueError(f"CVODE failed during initial simulation: {e}") from e
-        times_arr = np.array(result_arr["time"])  # copy before reset invalidates buffer
+        time_arr = np.array(result_arr[:, 0])  # copy before reset invalidates buffer
 
+        rr = self.getRoadrunner()
         rr.reset()
+        if self.start_time > 0:
+            _ = rr.simulate(0, self.start_time, 2)
         jacobians = []
         valid_times = []
-        for i, t in enumerate(times_arr):
+        for i, t in enumerate(time_arr):
             try:
                 if i == 0:
                     t2 = self.start_time + 1e-10
                     rr.simulate(self.start_time, t2, 2)
                 else:
-                    rr.simulate(times_arr[i - 1], t, 2)
+                    rr.simulate(time_arr[i - 1], t, 2)
             except RuntimeError as e:
                 raise ValueError(f"CVODE failed at t={t}: {e}") from e
             try:
@@ -516,7 +513,7 @@ def makeJacobians(self)->Tuple[np.ndarray, np.ndarray]:
             valid_times.append(t)
         if not jacobians:
             raise ValueError("No valid Jacobians could be computed (all timepoints failed, likely due to assignment-rule species).")
-        _ = self.timecourse  # Cache the timecourse DataFrame for later use in plotting, to avoid simulating again. This is done after the Jacobian collection loop to avoid state corruption that can cause getFullJacobian to segfault even after reset().
+        _ = self.timecourse_df  # Cache the timecourse DataFrame for later use in plotting, to avoid simulating again. This is done after the Jacobian collection loop to avoid state corruption that can cause getFullJacobian to segfault even after reset().
         return np.array(jacobians), np.array(valid_times)
 
     def _calculateEndtimeJacobian(self) -> float:

src/plot_options.py

@@ -3,37 +3,47 @@
 import matplotlib.pyplot as plt # type: ignore
 from typing import Optional, Tuple
 
+AX = plt.gca()
+FIG = plt.gcf()
+
 class PlotOptions(object):
 
     def __init__(self,
-            ax: Optional[plt.Axes] = None,   # type: ignore
+            ax=AX, 
+            fig = FIG,   # type: ignore
             title: Optional[str] = None,
-            xlabel: Optional[str] = None,
-            ylabel: Optional[str] = None,
+            xlabel: str = "time",
+            ylabel: str = "concentration",
             legend: bool = True,
             xlim: Optional[Tuple[float, float]] = None,
             ylim: Optional[Tuple[float, float]] = None,
             ):
-        if ax is None:
-            _, ax = plt.subplots()
+        if (AX == ax) and (FIG == fig):
+            fig, ax = plt.subplots()
         self.ax = ax
+        self.fig = fig
         self.title = title
         self.xlabel = xlabel
         self.ylabel = ylabel
         self.legend = legend
         self.xlim = xlim
         self.ylim = ylim
+
+    def to_dict(self):
+        return self.__dict__
 
     def apply(self):
         if self.title is not None:
-            self.ax.set_title(self.title)
+            self.ax.set_title(self.title) # type: ignore
         if self.xlabel is not None:
-            self.ax.set_xlabel(self.xlabel)
+            self.ax.set_xlabel(self.xlabel)  # type: ignore
         if self.ylabel is not None:
-            self.ax.set_ylabel(self.ylabel)
-        if self.legend:
-            self.ax.legend()
+            self.ax.set_ylabel(self.ylabel) # type: ignore
+        if isinstance(self.legend, bool) and self.legend:
+            self.ax.legend() # type: ignore
+        if isinstance(self.legend, list):
+            self.ax.legend(self.legend) # type: ignore
         if self.xlim is not None:
-            self.ax.set_xlim(self.xlim)
+            self.ax.set_xlim(self.xlim) # type: ignore
         if self.ylim is not None:
-            self.ax.set_ylim(self.ylim)
+            self.ax.set_ylim(self.ylim)  # type: ignore

src/trajectory.py

@@ -4,6 +4,7 @@
 from src.l_roadrunner import LRoadrunner, NULL_L_ROADRUNNER  # type: ignore
 import src.utils as utils
 from src.plot_options import PlotOptions  # type: ignore
+from src.score import Score  # type: ignore
 
 import collections
 import matplotlib.axes as maxes  # type: ignore
@@ -169,10 +170,6 @@ def max_cv(self) -> float:
             cv_arr[~np.isfinite(cv_arr)] = 0.0
         return np.max(cv_arr)
 
-    @property
-    def timecourse(self) -> pd.DataFrame:
-        return self.l_roadrunner.timecourse
-
     # ------------------------------------------------------------------
     # Public methods (alphabetical)
     # ------------------------------------------------------------------
@@ -277,7 +274,7 @@ def _calculateResiduals(params: Parameters, ispecies:int) -> np.ndarray:
             for i in range(self.num_species):
                 jacobian_arr[ispecies, i] = params[f'd{i}'].value
             prediction_arr = self._predict(jacobian_arr=jacobian_arr)[:, ispecies]
-            residual_arr = self.timecourse.iloc[:, ispecies].values - prediction_arr
+            residual_arr = self.timecourse_df.iloc[:, ispecies].values - prediction_arr
             return residual_arr[1:]  # Exclude timepoint 0 to avoid issues with initial state
         ##
         for ispecies, _ in enumerate(self.l_roadrunner.species_names):
@@ -305,6 +302,11 @@ def fromArrays(cls, jacobian_arr: np.ndarray, timepoint_arr: np.ndarray,
         jc._initialize(l_roadrunner, eigenvalues_collection_arr,
                 eigenvector_collection_arr, **kwargs)
         jc._sortArrays()
+        if l_roadrunner is not NULL_L_ROADRUNNER:
+            try:
+                jc.timecourse_df = l_roadrunner.timecourse_df.loc[jc.timepoint_arr]
+            except KeyError:
+                pass  # timepoint_arr not a subset of timecourse_df; keep full timecourse
         return jc
 
     def getCost(self, istart: int, iend: int) -> float:
@@ -557,23 +559,19 @@ def plot(self,
         plt.show()
         return PlotInfo(top_ax=ax1, bottom_ax=ax2, fig=fig)
 
-    def plotPredictions(self,
-            ax: Optional[plt.Axes] = None,   # type: ignore
-            ylim: Optional[Tuple[float, float]]=None,
-            xlim: Optional[Tuple[float, float]]=None,
-            model_name: str = "",
-            legend: bool = True,
-            ) -> PlotOptions:
+    def plotPrediction(self, **kwargs) -> PlotOptions:
         """
         Plot the predicted timecourse of simulation species concentrations.
         The first plot shows how the Jacobian changes over time relative to the centroid.
         The second plot shows the dynamics of the model's species concentrations
-        over time.
+        over time. Does not plot the first value since this is the initial state and not a prediction.
 
         Parameters
         ----------
         ax : Optional[plt.Axes]
             An optional matplotlib Axes
+        title: str
+            The title for the plot
         fig : Optional[plt.Figure]
             An optional matplotlib Figure object to use. If None, a new figure will be created.
         is_legend : bool
@@ -586,29 +584,29 @@ def plotPredictions(self,
             The model name
         """
         if hasattr(self.l_roadrunner, "getRoadrunner"):
-            roadrunner = self.l_roadrunner.getRoadrunner()
-            species_ids = roadrunner.getFloatingSpeciesIds()
-            data_arr = self.l_roadrunner.simulate(is_with_timepoints=True)
-            species_data = data_arr[:, 1:]  # Exclude time column
-            species_times = data_arr[:, 0]  # Extract time column
+            species_ids = self.l_roadrunner.species_names
+            species_times = self.timecourse_df.index.values
+            species_data = self.timecourse_df.values
         else:
             raise ValueError("Cannot plot species timecourse has a NULL LRoadrunner instance.")
-        prediction_df = self.predict()
+        pred_df = self.predict()
+        # Extract model_name before passing kwargs to PlotOptions (not a PlotOptions param)
+        model_name = kwargs.pop("model_name", "")
+        if model_name and "title" not in kwargs:
+            kwargs["title"] = f"{model_name}: Species Timecourse"
         # Timecourse plot
-        plot_options = PlotOptions(ax=ax,
-                legend=legend,
-                ylim=ylim,
-                xlim=xlim, title=f"{model_name}: Species Timecourse")
-        ax = plot_options.ax
+        plt_opt = PlotOptions(**kwargs)
+        ax = plt_opt.ax
         colors = [sns.color_palette("tab10")[i % 10] for i in range(len(species_ids))]
+        # Do separate loops so that legend works out correctly
+        for i, species_id in enumerate(species_ids):
+            ax.plot(species_times, species_data[:, i], # type: ignore
+                    label=species_id, color=colors[i], alpha=0.7)
         for i, species_id in enumerate(species_ids):
-            ax.plot(species_times, species_data[:, i], label=species_id, color=colors[i], alpha=0.7)
-            ax.scatter(species_times, prediction_df[species_id], s=8, alpha=0.7, color=colors[i])
-        ax.set_xlabel("Time")
-        ax.set_ylabel("Concentration")
-        ax.set_title(f"{model_name}: Species Timecourse")
-        plot_options.apply()
-        return plot_options
+            ax.scatter(species_times[1:], pred_df[species_id].values[1:],  # type: ignore
+                    s=8, alpha=0.7, color=colors[i])
+        plt_opt.apply()
+        return plt_opt
 
     def predict(self, **kwargs) -> pd.DataFrame:
         """
@@ -761,6 +759,7 @@ def _initialize(self, l_roadrunner: LRoadrunner,
         self._fitted_jacobian_arr = cn.NULL_ARRAY
         self._num_fit = num_fit
         self._jacobian_selection = jacobian_selection
+        self.timecourse_df = self.l_roadrunner.timecourse_df.copy()
 
     @staticmethod
     def _ivp(_: float, x: np.ndarray, jacobian_arr: np.ndarray) -> np.ndarray:
@@ -819,12 +818,12 @@ def ode(t: float, x: np.ndarray) -> np.ndarray:
         n_time = len(self.timepoint_arr)
         n_species = len(forcing_input_arr)
         result_arr = np.full((n_time, n_species), np.nan)
-        result_arr[0] = self.l_roadrunner.timecourse.iloc[0, :].values # type: ignore
+        result_arr[0] = self.l_roadrunner.timecourse_df.loc[self.timepoint_arr[0]].values # type: ignore
         try:
             with warnings.catch_warnings():
                 warnings.simplefilter("ignore", RuntimeWarning)
                 for itime, timepoint in enumerate(self.timepoint_arr[:-1]):
-                    initial_state_arr = self.l_roadrunner.timecourse.loc[timepoint].values # type: ignore
+                    initial_state_arr = self.l_roadrunner.timecourse_df.loc[timepoint].values # type: ignore
                     sol = solve_ivp(ode,
                             (timepoint, self.timepoint_arr[itime+1]),
                             initial_state_arr,
@@ -833,7 +832,7 @@ def ode(t: float, x: np.ndarray) -> np.ndarray:
                     if sol.success and sol.y.shape == (n_species, 1):
                         result_arr[itime+1] = sol.y.T
             return result_arr
-        except Exception:
+        except Exception as e:
             return np.full((n_time, n_species), np.nan)
 
     def _sortArrays(self) -> None: