SciCompMod / memilio

Modular spatio-temporal models for epidemic and pandemic simulations
https://scicompmod.github.io/memilio/
Apache License 2.0
55 stars 16 forks source link

Negative Compartments in stochastic models #1101

Open nijawa opened 3 months ago

nijawa commented 3 months ago

Bug description

The current implementation of the stochastic compartment models can result in negative compartments. This is problematic for multiple reasons:

  1. Negative Compartments do not make sense
  2. The stochastic part tries to take the square root of a negative number
  3. Clamping bounds for flows do not match so the program aborts

Version

Windows

To reproduce

  1. in "memilio\cpp\examples\sde_seirvv.cpp" set model.populations[{mio::sseirvv::InfectionState::InfectedV2}] = 0.1;

  2. 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

nijawa commented 3 months ago

Screenshot 2024-08-08 191321