Skip to content

Negative Compartments in stochastic models #1101

Open
@nijawa

Description

@nijawa

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

  • Attached labels, especially loc:: or model:: labels.
  • Linked to project

Metadata

Metadata

Labels

class::bugBugs found in the softwareloc::backendThis issue concerns the C++ backend implementation.model::sdeThis issue concerns any kind of stochastic differential equation-based model.

Type

No type

Projects

Status

Product Backlog

Milestone

No milestone

Relationships

None yet

Development

No branches or pull requests

Issue actions