Bug description
The current implementation of the stochastic compartment models can result in negative compartments. This is problematic for multiple reasons:
- Negative Compartments do not make sense
- The stochastic part tries to take the square root of a negative number
- Clamping bounds for flows do not match so the program aborts
Version
Windows
To reproduce
-
in "memilio\cpp\examples\sde_seirvv.cpp" set
model.populations[{mio::sseirvv::InfectionState::InfectedV2}] = 0.1;
-
in "memilio\cpp\models\sde_seirvv\model.h" set
std::initializer_list<uint32_t> seeds = {14159265u, 35897932u};
rng.seed(seeds);
(3. Use your favorite way to watch compartment values, for example the following before flow calculation)
printf("%f %f %f %f %f %f %f %f %f %f\n", y[(size_t)InfectionState::Susceptible], y[(size_t)InfectionState::ExposedV1],
y[(size_t)InfectionState::InfectedV1], y[(size_t)InfectionState::RecoveredV1],
y[(size_t)InfectionState::ExposedV2], y[(size_t)InfectionState::InfectedV2],
y[(size_t)InfectionState::RecoveredV2], y[(size_t)InfectionState::ExposedV1V2],
y[(size_t)InfectionState::InfectedV1V2], y[(size_t)InfectionState::RecoveredV1V2]);
Relevant log output
No response
Add any relevant information, e.g. used compiler, screenshots.
For some reason I cannot put a screenshot of the output here.
This should be an easy fix. From my tests, the error seems to originate from numerical errors in the Euler-Maruyama scheme together with the flow clamping:
compartment * inv_step_size * step_size != compartment
I suggest adding an error bound to the upper clamping bound. I suggest a multiplicative error term (so replace compartment * inv_step_size with compartment * inv_step_size * eps with eps being something like 1-1e-10 but I am not sure what a proper value would be here). An additive error term might result in illegal clamping bounds again.
Checklist
Bug description
The current implementation of the stochastic compartment models can result in negative compartments. This is problematic for multiple reasons:
Version
Windows
To reproduce
in "memilio\cpp\examples\sde_seirvv.cpp" set
model.populations[{mio::sseirvv::InfectionState::InfectedV2}] = 0.1;in "memilio\cpp\models\sde_seirvv\model.h" set
(3. Use your favorite way to watch compartment values, for example the following before flow calculation)
Relevant log output
No response
Add any relevant information, e.g. used compiler, screenshots.
Checklist