Gillespie timedep#1334
Merged
Merged
Conversation
…y on time. Fixes #1318. The direct method samples the waiting time from a propensity it treats as constant until the next reaction, which is only correct when the rate laws are time-homogeneous. When a rate depends explicitly on time (directly, or through a time-dependent assignment or rate rule), the SSA mean drifts off the ODE: a rate that is zero at t=0 froze the first reporting interval, and coarse output intervals trailed the true mean by about one interval. Now we check once, on the first step, whether any rate varies with time (the detection adds no random draws, so time-homogeneous models keep the original code path and produce bit-identical results). When a rate does depend on time, the waiting time tau satisfies integral over [t, t+tau] of s(u) du = -log(r1), where s is the total propensity re-evaluated as time advances; we accumulate that integral with an adaptive trapezoidal sub-step and select the reaction using the propensities at the firing time. This is the direct-method form of the integrated-propensity representation in D. F. Anderson, J. Chem. Phys. 127, 214107 (2007). The representation is exact; this implementation realizes it up to the quadrature error of the sub-step, so it is exact in the limit of small steps and converges as the tolerance is tightened, not bit-exact at a finite tolerance. Adds regression test TimeDependentPropensityMatchesODE, an immigration-death process with a time-dependent immigration rate and a closed-form ODE mean. The time-homogeneous path is byte-for-byte identical to before, and the stochastic test suite is unaffected (its rate laws are all time-homogeneous; the CSymbolTime models use time only in event triggers). The time-dependence check is a value-based heuristic; inspecting the rate-law expression trees for the time symbol would be more robust, which is noted in the code.
…rules The timeDependentRates comment claimed the fix handles time entering through a rate rule. It does not: Gillespie freezes the state vector between reaction events and never integrates a rate-rule variable, so setTime() alone does not advance it and it is neither detected nor tracked by the propensity integral. Drop the rate-rule claim and state the limitation explicitly. Comment-only; no behavior change.
Same fix as the header comment, in the duplicate that sat in the .cpp. A rate-rule variable is part of the frozen state, so it does not vary across the detection probes; time dependence entering through a rate rule is not detected (and not handled) here. Comment-only.
Gillespie 1318 timedep
* Actually check for the 'time' csymbol instead of integrating. * Make the comments less voluminous. * Update dependencies, since VS updated on github.
The latest MSVC seems to only use <random>, not <tr1/random>. Or something like that.
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
2 participants
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.
Slightly altered fix for #1318.