Fixes to Timecourse

Author: joseph-hellerstein

scripts/make_biomodels_timecourse.py

@@ -19,6 +19,7 @@ def main(
         first_model_num: int = 0,
         last_model_num: int = int(1e9),
         excluded_models: List[str] = EXCLUDED_MODELS,
+        is_initialize: bool = True, # Ignore existing serialized Timecourse when initializing (for testing).
 ) -> None:
     '''Serialize timecourses for all BioModels.
 
@@ -41,16 +42,21 @@ def main(
             last_model_num=last_model_num):
         model_name = item.model_name
         if not item.sbml_paths:
+            print(f"Skipping {model_name} (no SBML files)")
             continue
         pkl_path = os.path.join(cn.TIMECOURSE_SERIALIZATION_DIR,
                 f"{model_name}_timecourse.pkl")
-        if os.path.isfile(pkl_path):
+        if os.path.isfile(pkl_path) and (not is_initialize):
+            print(f"Skipping {model_name} (already serialized)")
             continue
         try:
             model = Model.makeBiomodel(model_name)
             timecourse = Timecourse(model=model, end_time=item.end_time)
             _ = timecourse.jacobian_collection_arr  # Force calculations
             path = timecourse.serialize()
+            serialized_timecourse = Timecourse.deserialize(path=path)
+            if not serialized_timecourse == timecourse:
+                raise ValueError(f"Deserialized timecourse does not match original for {model_name}.")  
         except Exception as e:
             print(f"Error processing {model_name}: {e}")
 

src/timecourse.py

@@ -26,7 +26,7 @@ def __init__(self, model: Model,
         end_time: Optional[float] = None,
         num_point: int = cn.NUM_POINTS,
         timecourse_df: pd.DataFrame = pd.DataFrame(),
-        jacobian_collection_arr: np.ndarray = np.array([])
+        jacobian_collection_arr: np.ndarray = np.array([]),
         ) -> None:
         """ 
         Parameters
@@ -39,6 +39,10 @@ def __init__(self, model: Model,
             Time to end the simulation.
         num_points : int
             Number of time points to simulate.
+        timecourse_df : pd.DataFrame
+            Optional pre-computed timecourse DataFrame (index: time, columns: species).
+        jacobian_collection_arr : np.ndarray
+            Optional pre-computed Jacobian collection (shape: [num_time_points, num_species, num_species]).
         """
         self.model = model
         self.start_time = start_time
@@ -48,6 +52,20 @@ def __init__(self, model: Model,
         self._timecourse_df = timecourse_df
         self._jacobian_collection_arr = jacobian_collection_arr
 
+    def __eq__(self, other: object) -> bool:
+        if not isinstance(other, Timecourse):
+            return NotImplemented
+        return (self.model == other.model and
+                bool(np.isclose(self.start_time, other.start_time)) and
+                (self.end_time == other.end_time
+                        if self.end_time is None or other.end_time is None
+                        else bool(np.isclose(self.end_time, other.end_time))) and
+                self.num_point == other.num_point and
+                bool(np.allclose(self.timecourse_df.values,
+                        other.timecourse_df.values)) and
+                bool(np.allclose(self.jacobian_collection_arr,
+                        other.jacobian_collection_arr)))
+
     def _updateEndtime(self, end_time: Optional[float]=None)->float | None:
         """Determine the end time and its source."""
         if end_time is not None:
@@ -83,6 +101,14 @@ def jacobian_collection_arr(self) -> np.ndarray:
             self._jacobian_collection_arr = simulation_result.jacobian_collection_arr
             self._timecourse_df = simulation_result.timecourse_df
         return self._jacobian_collection_arr
+    
+    def _checkSpeciesNames(self, names: List[str]) -> None:
+        """Check that the species names in the simulation result match the model."""
+        result_species = list(names)
+        if result_species != self.model.species_names:
+            raise ValueError(
+                    f"Simulation species {result_species} do not match "
+                    f"model species {self.model.species_names}.")
 
     def _simulate(self, is_jacobian_collection: bool = False) -> SimulationResult:
         """Create a Trajectory by running a simulation.
@@ -112,12 +138,18 @@ def _simulate(self, is_jacobian_collection: bool = False) -> SimulationResult:
         if self.start_time > 0:
             rr.simulate(0, self.start_time, 2)
         try:
-            result_arr = np.array(rr.simulate(self.start_time,
-                    self.end_time, self.num_point))
+            rr_result = rr.simulate(self.start_time, self.end_time, self.num_point)
         except Exception as e:
             raise ValueError(f"Simulation failed: {e}")
+        # Check column order before converting to ndarray (colnames lost after np.array).
+        # Skip the leading 'time' column and strip brackets from species names.
+        result_species = [
+                c[1:-1] if c.startswith("[") and c.endswith("]") else c
+                for c in rr_result.colnames[1:]  # type: ignore
+        ]
+        self._checkSpeciesNames(result_species)
+        result_arr = np.array(rr_result)
         timepoint_arr = result_arr[:, 0]
-        # FIXME: Use species names in NamedArray and sort by model.species_names --- IGNORE ---
         timecourse_df = pd.DataFrame(
                 result_arr[:, 1:],
                 index=timepoint_arr,
@@ -136,7 +168,10 @@ def _simulate(self, is_jacobian_collection: bool = False) -> SimulationResult:
                     rr.simulate(self.start_time, self.start_time + 1e-10, 2)
                 else:
                     rr.simulate(timepoint_arr[i - 1], t, 2)
-                jacobian_arr = np.array(rr.getFullJacobian()).copy()
+                jacobian_arr = rr.getFullJacobian()
+                self._checkSpeciesNames(jacobian_arr.rownames)
+                self._checkSpeciesNames(jacobian_arr.colnames)
+                jacobian_arr = np.array(jacobian_arr).copy()
                 if np.all(np.isclose(jacobian_arr, 0.0)):
                     raise ValueError(
                             f"Jacobian at t={t} is all zeros; model may be degenerate.")
@@ -158,8 +193,7 @@ def serialize(self) -> str:
         """
         if not self.model.model_name:
             raise ValueError("Model must have a name to serialize Timecourse.")
-        path = os.path.join(cn.TIMECOURSE_SERIALIZATION_DIR,
-                f"{self.model.model_name}_timecourse.pkl")
+        path = self.makeBiomodelSerializePath(self.model.model_name)
         dct = {
             "model": self.model,
             "start_time": self.start_time,
@@ -171,17 +205,38 @@ def serialize(self) -> str:
             pickle.dump(dct, f)
         return path
 
+    @staticmethod
+    def makeBiomodelSerializePath(model_name: str) -> str:
+        """
+        Get the expected path for a serialized Timecourse of a BioModel.
+
+        Parameters:
+            model_name (str): The name of the BioModel.
+        """
+        return os.path.join(cn.TIMECOURSE_SERIALIZATION_DIR, f"{model_name}_timecourse.pkl")
+
     @classmethod
-    def deserialize(cls, path: str) -> 'Timecourse':
+    def deserialize(cls, path: str = "", model_name: str = "") -> 'Timecourse':
         """
         Deserialize a Timecourse from a file
+        At least one of `path` or `model_name` must be provided.
+        If both are provided, `path` takes precedence.
 
         Parameters:
             path (str): The path to the serialized file.
+            model_name (str): The name of the BioModel (used if path is not specified).
 
         Returns:
             Timecourse: The deserialized Timecourse object.
         """
+        if not path and not model_name:
+            raise ValueError("At least one of `path` or `model_name` must be provided.")
+        if not path:
+            path = cls.makeBiomodelSerializePath(model_name)
+        # Check if the file exists
+        if not os.path.isfile(path):
+            raise FileNotFoundError(f"No serialized Timecourse found at {path}")    
+        # Deserialize
         with open(path, 'rb') as f:
             dct = pickle.load(f)
         return cls(

tests/test_timecourse.py

@@ -80,6 +80,46 @@ def test_jacobian_collection_arr_stored(self) -> None:
         self.assertGreater(tc._jacobian_collection_arr.size, 0)
 
 
+class TestTimecourseJacobianCollectionArr(unittest.TestCase):
+    """Tests for Timecourse.jacobian_collection_arr property."""
+
+    def test_returns_ndarray(self) -> None:
+        if IGNORE_TESTS:
+            return
+        tc = _makeTimecourse()
+        self.assertIsInstance(tc.jacobian_collection_arr, np.ndarray)
+
+    def test_shape(self) -> None:
+        if IGNORE_TESTS:
+            return
+        tc = _makeTimecourse()
+        self.assertEqual(tc.jacobian_collection_arr.shape,
+                (NUM_POINT, NUM_SPECIES, NUM_SPECIES))
+
+    def test_returns_prepopulated_array(self) -> None:
+        if IGNORE_TESTS:
+            return
+        tc = _makeTimecourse()
+        np.testing.assert_array_equal(
+                tc.jacobian_collection_arr, tc._jacobian_collection_arr)
+
+    def test_cached_on_second_access(self) -> None:
+        if IGNORE_TESTS:
+            return
+        tc = _makeTimecourse()
+        first = tc.jacobian_collection_arr
+        second = tc.jacobian_collection_arr
+        self.assertIs(first, second)
+
+    def test_simulation_populates_timecourse_df(self) -> None:
+        if IGNORE_TESTS:
+            return
+        tc = Timecourse(model=_makeModel(), end_time=10.0)
+        self.assertTrue(tc._timecourse_df.empty)
+        _ = tc.jacobian_collection_arr
+        self.assertFalse(tc._timecourse_df.empty)
+
+
 class TestTimecourseSerialize(unittest.TestCase):
     """Tests for Timecourse.serialize."""
 
@@ -121,6 +161,71 @@ def test_raises_without_model_name(self) -> None:
             tc.serialize()
 
 
+class TestTimecourseEq(unittest.TestCase):
+    """Tests for Timecourse.__eq__."""
+
+    def test_equal_timecourses(self) -> None:
+        if IGNORE_TESTS:
+            return
+        tc1 = _makeTimecourse()
+        tc2 = _makeTimecourse()
+        self.assertEqual(tc1, tc2)
+
+    def test_different_model_not_equal(self) -> None:
+        if IGNORE_TESTS:
+            return
+        tc1 = _makeTimecourse()
+        tc2 = _makeTimecourse()
+        tc2.model = Model(ANTIMONY_MODEL, model_name="other")
+        self.assertNotEqual(tc1, tc2)
+
+    def test_different_start_time_not_equal(self) -> None:
+        if IGNORE_TESTS:
+            return
+        tc1 = _makeTimecourse()
+        tc2 = _makeTimecourse()
+        tc2.start_time = 99.0
+        self.assertNotEqual(tc1, tc2)
+
+    def test_different_end_time_not_equal(self) -> None:
+        if IGNORE_TESTS:
+            return
+        tc1 = _makeTimecourse()
+        tc2 = _makeTimecourse()
+        tc2.end_time = 99.0
+        self.assertNotEqual(tc1, tc2)
+
+    def test_different_num_point_not_equal(self) -> None:
+        if IGNORE_TESTS:
+            return
+        tc1 = _makeTimecourse()
+        tc2 = _makeTimecourse()
+        tc2.num_point = 999
+        self.assertNotEqual(tc1, tc2)
+
+    def test_different_timecourse_df_not_equal(self) -> None:
+        if IGNORE_TESTS:
+            return
+        tc1 = _makeTimecourse()
+        tc2 = _makeTimecourse()
+        tc2._timecourse_df = tc2._timecourse_df * 2
+        self.assertNotEqual(tc1, tc2)
+
+    def test_different_jacobian_not_equal(self) -> None:
+        if IGNORE_TESTS:
+            return
+        tc1 = _makeTimecourse()
+        tc2 = _makeTimecourse()
+        tc2._jacobian_collection_arr = tc2._jacobian_collection_arr * 2
+        self.assertNotEqual(tc1, tc2)
+
+    def test_not_equal_to_non_timecourse(self) -> None:
+        if IGNORE_TESTS:
+            return
+        tc = _makeTimecourse()
+        self.assertIs(tc.__eq__("not a timecourse"), NotImplemented)
+
+
 class TestTimecourseRoundtrip(unittest.TestCase):
     """Roundtrip tests for Timecourse.serialize / Timecourse.deserialize."""
 
@@ -138,44 +243,17 @@ def _serializeAndDeserialize(self) -> tuple:
         restored = Timecourse.deserialize(path)
         return original, restored
 
-    def test_timecourse_df_roundtrip(self) -> None:
-        if IGNORE_TESTS:
-            return
-        original, restored = self._serializeAndDeserialize()
-        pd.testing.assert_frame_equal(
-                original._timecourse_df, restored._timecourse_df)
-
-    def test_jacobian_collection_arr_roundtrip(self) -> None:
+    def test_roundtrip_equal(self) -> None:
         if IGNORE_TESTS:
             return
         original, restored = self._serializeAndDeserialize()
-        np.testing.assert_array_equal(
-                original._jacobian_collection_arr,
-                restored._jacobian_collection_arr)
+        self.assertEqual(original, restored)
 
-    def test_start_time_roundtrip(self) -> None:
+    def test_deserialize_raises_with_empty_path(self) -> None:
         if IGNORE_TESTS:
             return
-        original, restored = self._serializeAndDeserialize()
-        self.assertEqual(original.start_time, restored.start_time)
-
-    def test_end_time_roundtrip(self) -> None:
-        if IGNORE_TESTS:
-            return
-        original, restored = self._serializeAndDeserialize()
-        self.assertEqual(original.end_time, restored.end_time)
-
-    def test_num_points_roundtrip(self) -> None:
-        if IGNORE_TESTS:
-            return
-        original, restored = self._serializeAndDeserialize()
-        self.assertEqual(original.num_point, restored.num_point)
-
-    def test_model_name_roundtrip(self) -> None:
-        if IGNORE_TESTS:
-            return
-        original, restored = self._serializeAndDeserialize()
-        self.assertEqual(original.model.model_name, restored.model.model_name)
+        with self.assertRaises(Exception):
+            Timecourse.deserialize("")
 
 
 if __name__ == "__main__":