pratikunterwegs
commented
4 months ago A few thoughts on this issue:
- Introducing large numbers of vaccinations converts an otherwise non-stiff system into a potentially stiff system (depending on other compartmental transitions), which poses an issue for the simple RK4 stepper used in {epidemics} models.
- I haven't been able to find a lot of standard guidance on what to do for stiff systems. Some steppers are held to work well with stiff systems including the Rosenbrock 4 stepper in Boost (and this is used in Epidemics.jl).
- The Boost Rosenbrock 4 stepper and some others only appear to work with Boost's uBLAS array type, and not with Eigen arrays which we use in {epidemics} and {finalsize}. I could move to using uBLAS arrays in {epidemics} but the library is old and its own documentation points to Eigen as a potential replacement.
- Furthermore, the Rosenbrock stepper (and stiff-system steppers in Boost generally) seems to require a pre-calculated Jacobian matrix to be efficient. This becomes challenging when interventions, time-dependence, and vaccinations are added to a scenario as the Jacobian would, if I understand correctly, need to be recalculated in each timestep.
The workable solution is probably to manually change the number of vaccinations so that the $S \rightarrow V$ flow is always $\leq S$.
Happy to hear solutions on this that could work for this situation.
This issue is a subset of #224, where simulating vaccination leads to negative numbers of susceptibles as per-day vaccination counts may be greater than available susceptibles towards the end of an outbreak.
One solution would be the use of a stepper that is more suitable for stiff ODE systems. It should especially be able to handle > 1% instantaneous changes as these are expected in vaccination scenarios. Issue #224 shows that integration fails and produces negative values for the current fixed size stepper used. Options from Boost are listed here, and some guidance on related functionality is here.
An alternative solution is to build on PR #202 and add a simple
if elsecheck which sets vaccinations to the lesser of either available doses or susceptibles.