Repository navigation
fix: carry pre-bolus step size across dose ReInit (fixes #25) - #26
Merged
Merged
Conversation
CVODE's initial-step estimator (CVHin) picks h₀ from ||f(t,y)||_WRMS at the post-bolus state. When a bolus jolts a fast-binding subsystem (e.g. nivolumab → V_C.aPD1, which couples into the PD1/aPD1 synapse binding cross-section), the norm is huge and CVHin returns a sub-ULP h₀. At solver times O(1e8) s, the resulting step doesn't change t in double precision and the CV_ONE_STEP loop's "did not advance" guard aborts — even though the corrector succeeded with zero error-test and Newton convergence failures. setupSamplingRun now accepts an optional h0_hint that bypasses CVHin via CVodeSetInitStep. The sampling loop captures pre-event h_cur and passes std::min(h_cur, 1.0 s) as the hint — well above ULP, conservative enough that BDF converges without retries on most discontinuities. This mirrors how SimBiology carries h-hints across events. Also adds: - CVODEBase::getStepStats() exposing per-step CVODE counters (order, h_last, h_cur, error-test fails, Newton conv fails, Jac/RHS evals). Opt-in via QSP_TRACE_CVODE env var in qsp_sim_main for diagnosing stiff-step issues without rebuilding. - Optional output_times_csv plumbing through parity infra (run_matlab_trajectories + export_matlab_trajectories.m). Pinning MATLAB's grid to the C++ grid eliminates spurious "post-bolus value at pre-bolus timestamp" failures from linear interpolation across a dose discontinuity. Pattern lifted from paper/benchmark/. Verified against pdac-build gvax_nivo_neoadjuvant_zheng2022 scenario: the day-7 nivo bolus no longer collapses h to 1e-9 s, and grid-aligned MATLAB/C++ parity passes at rtol=5e-2, atol=1e-6 over 30 d.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Summary
Fixes #25 — CVODE collapse at the day-7 nivolumab bolus in the
gvax_nivo_neoadjuvant_zheng2022scenario.Root cause was not stiffness or a missing ReInit (both were red herrings in the original triage). Per-step CVODE diagnostics show the post-bolus step had zero error-test failures and zero Newton convergence failures — CVODE succeeded, but with
h = 1.4e-9 sat solver timet ≈ 1.6e8 s, which is below the double-precision ULP. The CV_ONE_STEP loop's "did not advance" guard then aborts the run, even though CVODE was about to growhto8e-6 son the next call.The bad
hcame from CVODE's initial-step estimator (CVHin), which picksh₀from||f(t,y)||_WRMSat the post-bolus state. The nivo perturbation toV_C.aPD1couples into the fast PD1/aPD1 synapse binding cross-section, the norm explodes, andh₀ → sub-ULP.Fix
setupSamplingRungains an optionalh0_hintthat bypasses CVHin viaCVodeSetInitStep.qsp_sim_main.cppcaptures pre-eventh_curand passesmin(h_cur, 1.0 s)as the hint. The 1 s cap is conservative enough that BDF converges without Newton retries on most discontinuities; still well above ULP at any realistic solver time.This mirrors how SimBiology carries h-hints across events, which is why MATLAB never tripped this.
Also included
CVODEBase::getStepStats()— per-step CVODE diagnostics (order,h_last,h_cur, error-test fails, Newton conv fails, Jac/RHS evals). Opt-in viaQSP_TRACE_CVODEenv var inqsp_sim. This was the diagnostic that distinguished "CVODE struggling" from "CVODE took a sub-ULP step and the outer loop bailed."output_times_csvplumbed throughrun_matlab_trajectoriesandexport_matlab_trajectories.m. When provided, MATLAB pins its grid to the C++ grid — eliminating spurious "post-bolus value at pre-bolus timestamp" parity failures from linear interpolation across a discontinuity. Pattern lifted frompaper/benchmark/.Verification
Against pdac-build's
tests/test_cpp_matlab_parity.py::test_sbi_scenario_paritywithoutput_times_csvwired up:gvax_nivo— PASS (was crashing at t=7d)gvax_neoadjuvant— PASSbaseline_no_treatment/clinical_progression— pre-existing MATLAB bug inschedule_dosingfor drugless scenarios; unrelated.Follow-ups (not in this PR)
output_times_csvthrough its parity test.scripts/evolve_to_diagnosis.mstill uses the legacyspqsp_resource()lookup forhealthy_state.yaml— part of issue Rename CancerVCT namespace and scrub legacy SPQSP_PDAC references from runtime #19's SPQSP_PDAC scrub.schedule_dosingbug for drugless scenarios.H_PD1_APCshows a ~8.7% transient ~10 min post-bolus ongvax_nivo; still inside the documented 5–20% envelope but worth a closer look.