Script to run timecourse in parallel.

Author: joseph-hellerstein

data/badmodels.txt

@@ -48,6 +48,7 @@ BIOMD0000000133
 BIOMD0000000134
 BIOMD0000000135
 BIOMD0000000136
+BIOMD0000000137
 BIOMD0000000138
 BIOMD0000000139
 BIOMD0000000140
@@ -146,6 +147,7 @@ BIOMD0000000403
 BIOMD0000000404
 BIOMD0000000408
 BIOMD0000000416
+BIOMD0000000424
 BIOMD0000000429
 BIOMD0000000437
 BIOMD0000000450

docs/research_agenda.md

@@ -1,19 +1,23 @@
 # Research Agenda
 
-## To what extent is BioModels linear?
+## How well can BioModels be modelled by a system of linear differential equaitons?
 - [x] Fit first order monomial to all models and calculate $R^2$ for species time courses.
     - Minimum $R^2$ is model $R^2$
     - Individual $R^2$ are for species.
 - [x] Evaluate the density of coefficients in the Jacobian.
 - [ ] For the linear models, assess their dimensionality to see if dimension reduction is possible
 - [ ] Characterize the linear models based on the what is being model and possibly other characteristics.
-
-## What are the main reasons for nonlinear behavior?
-- [ ] Analyze the nonlinear models to determine which species are nonlinear and how/when the Jacobian changes to look at reactions.
 
   ## How robust is linearity to perturbations of initial values?
 - [ ] repeat the linear studies with perturbations of $\pm 5\%$, $\pm 19\%$, $\pm 20\%$, and $\pm 50\%$.
 - [ ] Can robustness be improved by training the regression on perturbation data?
 
+## What are the main reasons for nonlinear behavior?
+- [ ] Analyze the nonlinear models to determine which species are nonlinear and how/when the Jacobian changes to look at reactions.
+
 ## Are some nonlinear models piecewise linear?
-- [ ] Use a standard package for partitioning regressions to see if linearity can be achieved in segements of the time course.
+- [ ] Use a standard package for partitioning regressions to see if linearity can be achieved in segements of the time course.
+  - [ ] Normalize Jacobians
+  - [ ] k-means cluster with cluster distances with minimum cluster size.
+  - [ ] SystemDiscovery for each cluster 
+  - [ ] Prediction using: (a) ${\bf x} (t)$ using SystemDiscover of cluster at $t$ similarly for adjacent time points; (b) apply gaussian kernel on points; (c) denormalize

scripts/run_parallel_timecourse.sh

@@ -0,0 +1,47 @@
+#!/usr/bin/env bash
+# Run make_biomodels_timecourse.py in parallel across N workers.
+# Each worker receives an evenly divided slice of the 1078 BioModels.
+#
+# Usage: bash scripts/run_parallel_timecourse.sh <number_of_processes>
+
+#set -euo pipefail
+
+if [ $# -ne 1 ] || ! [[ $1 =~ ^[1-9][0-9]*$ ]]; then
+    echo "Usage: $0 <number_of_processes>" >&2
+    exit 1
+fi
+
+NUM_PROCESSES=$1
+TOTAL_MODELS=1078
+TOTAL_MODELS=20
+FIRST_MODEL=1
+LAST_MODEL=$(( FIRST_MODEL + TOTAL_MODELS - 1 ))
+
+# Ceiling division: chunk = ceil(TOTAL / N)
+CHUNK=$(( (TOTAL_MODELS + NUM_PROCESSES - 1) / NUM_PROCESSES ))
+
+SCRIPT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)"
+PROJECT_DIR="$(dirname "$SCRIPT_DIR")"
+
+source "$PROJECT_DIR/activate.sh"
+mkdir -p "$PROJECT_DIR/logs"
+
+echo "Launching $NUM_PROCESSES workers (~$CHUNK models each)..."
+
+for i in $(seq 0 $(( NUM_PROCESSES - 1 ))); do
+    FIRST=$(( FIRST_MODEL + i * CHUNK ))
+    LAST=$(( FIRST + CHUNK - 1 ))
+    if [ "$LAST" -gt "$LAST_MODEL" ]; then
+        LAST=$LAST_MODEL
+    fi
+    LOG="$PROJECT_DIR/logs/timecourse_${FIRST}_${LAST}.log"
+    echo "  Worker $(( i + 1 )): models $FIRST–$LAST  →  $LOG"
+    python "$SCRIPT_DIR/make_biomodels_timecourse.py" \
+        --first_model_num "$FIRST" \
+        --last_model_num  "$LAST"  \
+        > "$LOG" 2>&1 &
+done
+
+echo "All workers launched.  Waiting for completion..."
+wait
+echo "Done."

src/timecourse.py

@@ -76,9 +76,9 @@ def __eq__(self, other: object) -> bool:
                         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
+                        other.timecourse_df.values, equal_nan=True)) and
                 bool(np.allclose(self.jacobian_collection_arr,
-                        other.jacobian_collection_arr)))
+                        other.jacobian_collection_arr, equal_nan=True)))
 
     def _updateEndtime(self, end_time: Optional[float]=None):
         """Determine the end time and its source."""

tests/test_timecourse.py

@@ -18,6 +18,8 @@
 from timecourse import Timecourse  # type: ignore
 
 IGNORE_TESTS = False
+HAS_REAL_ZIP = os.path.isfile(cn.TIMECOURSE_ZIP_PATH)
+BIOMODEL_53 = "BIOMD0000000053"
 
 ANTIMONY_MODEL = """
 S1 -> S2; k1*S1
@@ -640,5 +642,17 @@ def test_subplot_ylabel_is_concentration(self) -> None:
             self.assertEqual(ax.get_ylabel(), "concentration")
 
 
+class TestTimecourseEqRealBiomodel(unittest.TestCase):
+    """Verifies __eq__ with a real BioModel timecourse loaded from the zip archive."""
+
+    @unittest.skipUnless(HAS_REAL_ZIP, "Real timecourse zip not found")
+    def test_biomodel53_equals_itself(self) -> None:
+        if IGNORE_TESTS:
+            return
+        from timecourse_iterator import TimecourseIterator  # type: ignore
+        tc = TimecourseIterator().getTimecourse(BIOMODEL_53)
+        self.assertEqual(tc, tc)
+
+
 if __name__ == "__main__":
     unittest.main()