diff --git a/cpp/include/qsp_sim_core/CVODEBase.h b/cpp/include/qsp_sim_core/CVODEBase.h index ed7d0fa..7c6da46 100644 --- a/cpp/include/qsp_sim_core/CVODEBase.h +++ b/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); diff --git a/cpp/src/CVODEBase.cpp b/cpp/src/CVODEBase.cpp index 72cd7e2..51a1e6c 100644 --- a/cpp/src/CVODEBase.cpp +++ b/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(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(hlast); + s.cur_h = static_cast(hcur); + return s; +} // ----- Evolve-cache full-state serialization ----------------------------- // // Pure data dump: ODE_system state vectors + delay-event queue + trigger diff --git a/cpp/src/qsp_sim_main.cpp b/cpp/src/qsp_sim_main.cpp index da66b04..3866606 100644 --- a/cpp/src/qsp_sim_main.cpp +++ b/cpp/src/qsp_sim_main.cpp @@ -47,6 +47,7 @@ */ #include #include +#include #include #include #include @@ -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 diff --git a/src/qsp_codegen/parity/__init__.py b/src/qsp_codegen/parity/__init__.py index d275fbc..e9255c1 100644 --- a/src/qsp_codegen/parity/__init__.py +++ b/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};") diff --git a/src/qsp_codegen/parity/export_matlab_trajectories.m b/src/qsp_codegen/parity/export_matlab_trajectories.m index 53a64f3..ff071f4 100644 --- a/src/qsp_codegen/parity/export_matlab_trajectories.m +++ b/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)