fix: carry pre-bolus step size across dose ReInit (fixes #25)

cpp/include/qsp_sim_core/CVODEBase.h

@@ -74,6 +74,26 @@ class CVODEBase
 	//! / resetSolver. Wraps CVodeGetNumSteps for instrumentation.
 	long getNumSteps() const;
 
+	//! Per-step CVODE diagnostics for stiffness debugging. All counters are
+	//! cumulative since the last CVodeReInit (i.e. since the last
+	//! setupSamplingRun / resetSolver call). Take diffs across steps to get
+	//! per-step values. last_h / cur_h are the actually-taken last step and
+	//! the controller's predicted next step; last_order is the BDF order
+	//! used on the most recent successful step.
+	struct StepStats {
+		long nst;          // cumulative steps
+		long nfe;          // RHS evaluations
+		long netf;         // error-test failures (controller cut h)
+		long ncfn;         // nonlinear-solver convergence failures (Newton)
+		long nje;          // Jacobian evaluations
+		long nni;          // total Newton iterations
+		long nsetups;      // linear-solver setups
+		int  last_order;   // BDF order used on last step
+		double last_h;     // step size of last successful step
+		double cur_h;      // controller's predicted next step
+	};
+	StepStats getStepStats() const;
+
 	//! One-time setup before the first simOdeSample call: sets the CVODE
 	//! stop time to t_end_of_sim so CV_NORMAL stepping won't walk past it,
 	//! and runs any pending initial-assignment events at t0.
@@ -82,7 +102,17 @@ class CVODEBase
 	//! internal time when entering a sampling run after a prior
 	//! simOdeStep-based phase (e.g. evolve_to_diagnosis); if the fast path
 	//! runs from fresh ICs, leave it at the 0.0 default.
-	void setupSamplingRun(double tEndOfSim, double t0 = 0.0);
+	//!
+	//! `h0_hint` (optional, in solver time units): if > 0, passed to
+	//! CVodeSetInitStep so the post-ReInit step doesn't start from the
+	//! CVHin estimator. Use this when re-initing across a discontinuity
+	//! (e.g. a bolus event): the CVHin estimator looks at ||f(t,y)||_WRMS
+	//! at the new state, which can be huge if the perturbation jolts a
+	//! fast-binding subsystem, producing a sub-ULP h₀ that stalls
+	//! progress. Passing the pre-event h_cur (or any sane O(seconds)
+	//! value) restores MATLAB-parity behavior across events.
+	void setupSamplingRun(double tEndOfSim, double t0 = 0.0,
+	                      double h0_hint = 0.0);
 
 	//! examples of optional output
 	void PrintFinalStats(void *cvode_mem);

cpp/src/CVODEBase.cpp

@@ -542,11 +542,20 @@ void CVODEBase::resetSolver(realtype t0, realtype t1){
 	stale state, not the properly-initialized one.
 
 */
-void CVODEBase::setupSamplingRun(double tEndOfSim, double t0){
+void CVODEBase::setupSamplingRun(double tEndOfSim, double t0, double h0_hint){
 	restore_y();
 	resolveEvents(t0);
 	CVodeReInit(_cvode_mem, t0, _y);
 	CVodeSetStopTime(_cvode_mem, tEndOfSim);
+	if (h0_hint > 0.0) {
+		// Skip CVHin's heuristic h₀ pick. CVHin estimates h₀ from
+		// ||f(t0,y0)||_WRMS, which collapses to sub-ULP at solver times
+		// O(1e8) s when a bolus jolts a fast-binding subsystem — the
+		// step then advances t by less than one double-precision ULP
+		// and the outer loop's "did not advance" guard aborts the run.
+		// Passing the pre-event h_cur as the hint mirrors SimBiology.
+		CVodeSetInitStep(_cvode_mem, static_cast<realtype>(h0_hint));
+	}
 }
 
 /*! Advance CVODE from its current internal time to tEnd in CV_NORMAL mode.
@@ -600,6 +609,24 @@ long CVODEBase::getNumSteps() const {
 	}
 	return n_steps;
 }
+
+CVODEBase::StepStats CVODEBase::getStepStats() const {
+	StepStats s{};
+	CVodeGetNumSteps(_cvode_mem, &s.nst);
+	CVodeGetNumRhsEvals(_cvode_mem, &s.nfe);
+	CVodeGetNumErrTestFails(_cvode_mem, &s.netf);
+	CVodeGetNumNonlinSolvConvFails(_cvode_mem, &s.ncfn);
+	CVodeGetNumJacEvals(_cvode_mem, &s.nje);
+	CVodeGetNumNonlinSolvIters(_cvode_mem, &s.nni);
+	CVodeGetNumLinSolvSetups(_cvode_mem, &s.nsetups);
+	CVodeGetLastOrder(_cvode_mem, &s.last_order);
+	realtype hlast = 0.0, hcur = 0.0;
+	CVodeGetLastStep(_cvode_mem, &hlast);
+	CVodeGetCurrentStep(_cvode_mem, &hcur);
+	s.last_h = static_cast<double>(hlast);
+	s.cur_h = static_cast<double>(hcur);
+	return s;
+}
 // ----- Evolve-cache full-state serialization -----------------------------
 //
 // Pure data dump: ODE_system state vectors + delay-event queue + trigger

cpp/src/qsp_sim_main.cpp

@@ -47,6 +47,7 @@
  */
 #include <algorithm>
 #include <cstdint>
+#include <cstdlib>
 #include <cstring>
 #include <fstream>
 #include <iomanip>
@@ -878,6 +879,20 @@ int main(int argc, char* argv[]) {
         ? t_stop : std::min(dose_boundaries.front(), t_stop);
     ode.setupSamplingRun(seg_end, t);
 
+    // Per-step CVODE diagnostics: when QSP_TRACE_CVODE=1, log post-step
+    // stats every step (or every N steps if set to an integer > 1). Counters
+    // are cumulative since the last setupSamplingRun, so we keep the previous
+    // snapshot to emit deltas (errfail/convfail/jac evals per step). Diagnoses
+    // h-collapse: is the controller cutting h via error-test failures (stiff
+    // RHS at this tol) or via nonlinear-conv failures (bad Jacobian / linsol)?
+    const char* trace_env = std::getenv("QSP_TRACE_CVODE");
+    int trace_every = 0;
+    if (trace_env && *trace_env) {
+        trace_every = std::atoi(trace_env);
+        if (trace_every < 1) trace_every = 1;
+    }
+    CVODEBase::StepStats prev_stats{};
+
     // CV_ONE_STEP loop with cadence floor (D4). Each iteration advances by
     // one CVODE internal step, then emits a row if either:
     //   - elapsed solver time since last dump >= min_cadence_solver, or
@@ -890,11 +905,40 @@ int main(int argc, char* argv[]) {
         const double t_prev = t;
         t = ode.simOdeStepOne(seg_end);
 
+        if (trace_every > 0 && (step % trace_every) == 0) {
+            const auto s = ode.getStepStats();
+            std::cerr << "  [cvode] t=" << (t - t_offset) / time_factor
+                      << "d step=" << step
+                      << " nst=" << s.nst
+                      << " ord=" << s.last_order
+                      << " h_last=" << s.last_h
+                      << " h_cur=" << s.cur_h
+                      << " dnetf=" << (s.netf - prev_stats.netf)
+                      << " dncfn=" << (s.ncfn - prev_stats.ncfn)
+                      << " dnje=" << (s.nje - prev_stats.nje)
+                      << " dnfe=" << (s.nfe - prev_stats.nfe)
+                      << " dnni=" << (s.nni - prev_stats.nni)
+                      << std::endl;
+            prev_stats = s;
+        }
+
         const bool at_seg_end = !(t < seg_end);  // robust to float roundoff
         const bool at_t_stop = !(t < t_stop);
         const bool cadence_due = (t - t_last_dump) >= min_cadence_solver;
 
+        // Carry the pre-event step size across the dose discontinuity so the
+        // post-ReInit step doesn't fall into CVHin's degenerate h₀-pick when
+        // the bolus jolts a fast-binding subsystem (e.g. PD1/aPD1 synapse
+        // binding after a V_C.aPD1 bolus). See setupSamplingRun for context.
+        // Cap the hint at 1 s of solver time: pre-event h_cur can be O(1e4)s
+        // on smooth segments, which is far too aggressive for the post-bolus
+        // state (CVODE burns dozens of Newton iters cutting it back down).
+        // 1 s is well above any plausible ULP at solver times ≤ 1e10 s and
+        // is conservative enough that BDF converges in one Newton iter.
+        double h_carry = 0.0;
         if (at_seg_end && (next_dose_idx < dose_boundaries.size())) {
+            constexpr double H_CARRY_CAP = 1.0;  // seconds
+            h_carry = std::min(ode.getStepStats().cur_h, H_CARRY_CAP);
             // Apply bolus first so the dump captures post-dose state. Matches
             // MATLAB sbiosimulate dose_schedule semantics and the v2 path.
             apply_at(t);
@@ -918,7 +962,11 @@ int main(int argc, char* argv[]) {
             ++next_dose_idx;
             seg_end = (next_dose_idx < dose_boundaries.size())
                 ? std::min(dose_boundaries[next_dose_idx], t_stop) : t_stop;
-            ode.setupSamplingRun(seg_end, t);
+            ode.setupSamplingRun(seg_end, t, h_carry);
+            // CVODE counters are reset by CVodeReInit; baseline the
+            // per-step delta snapshot so the first post-dose log row
+            // reports clean per-step deltas instead of a negative jump.
+            prev_stats = CVODEBase::StepStats{};
         }
 
         // Guard against CV_ONE_STEP returning the same t when CVODE thinks

src/qsp_codegen/parity/__init__.py

@@ -154,6 +154,7 @@ def run_matlab_trajectories(
     stop_time: float = 365.0,
     scenario_yaml: Optional[Path] = None,
     evolve_function_name: Optional[str] = None,
+    output_times_csv: Optional[Path] = None,
     matlab_binary: str = "matlab",
     timeout: float = 600.0,
 ) -> Path:
@@ -183,6 +184,14 @@ def run_matlab_trajectories(
     ]
     if scenario_yaml is not None:
         cmd_parts.append(f"scenario_yaml='{scenario_yaml}';")
+    if output_times_csv is not None:
+        # When set, MATLAB reads this CSV (one float per line) and uses
+        # those values as its SolverOptions.OutputTimes — pinning MATLAB's
+        # sample grid to the C++ grid so the parity comparison is
+        # row-aligned with no interpolation. Without this, linear interp
+        # across a bolus discontinuity produces spurious post-bolus values
+        # at pre-bolus timestamps. Match how paper/benchmark/ does it.
+        cmd_parts.append(f"output_times_csv='{output_times_csv}';")
     if evolve_function_name is not None:
         # MATLAB scripts use the `matlab_evolve_function` handle variable.
         cmd_parts.append(f"matlab_evolve_function=@{evolve_function_name};")

src/qsp_codegen/parity/export_matlab_trajectories.m

@@ -54,7 +54,14 @@
     stop_time = 365;
 end
 cfg.StopTime = stop_time;
-cfg.SolverOptions.OutputTimes = (0:0.1:stop_time)';
+% If output_times_csv is provided, pin MATLAB's sample grid to the C++
+% grid so the parity comparison is row-aligned (no interpolation across
+% bolus discontinuities). Otherwise use the default 0.1 d grid.
+if exist('output_times_csv', 'var') && ~isempty(output_times_csv)
+    cfg.SolverOptions.OutputTimes = readmatrix(output_times_csv);
+else
+    cfg.SolverOptions.OutputTimes = (0:0.1:stop_time)';
+end
 
 % Optional: override model values with those from a param_all XML so MATLAB
 % simulates with the exact same ICs / parameters as the C++ run.
@@ -198,7 +205,11 @@
     end
     % Reinstate the scenario's output grid in case the evolve mutated it.
     cfg.StopTime = stop_time;
-    cfg.SolverOptions.OutputTimes = (0:0.1:stop_time)';
+    if exist('output_times_csv', 'var') && ~isempty(output_times_csv)
+        cfg.SolverOptions.OutputTimes = readmatrix(output_times_csv);
+    else
+        cfg.SolverOptions.OutputTimes = (0:0.1:stop_time)';
+    end
 end
 
 if isempty(dose_schedule)