Line data    Source code

       1             : /* 
       2             : * Copyright (C) 2020-2024 MEmilio
       3             : *
       4             : * Authors: Nils Wassmuth, Rene Schmieding, Martin J. Kuehn
       5             : *
       6             : * Contact: Martin J. Kuehn <Martin.Kuehn@DLR.de>
       7             : *
       8             : * Licensed under the Apache License, Version 2.0 (the "License");
       9             : * you may not use this file except in compliance with the License.
      10             : * You may obtain a copy of the License at
      11             : *
      12             : *     http://www.apache.org/licenses/LICENSE-2.0
      13             : *
      14             : * Unless required by applicable law or agreed to in writing, software
      15             : * distributed under the License is distributed on an "AS IS" BASIS,
      16             : * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
      17             : * See the License for the specific language governing permissions and
      18             : * limitations under the License.
      19             : */
      20             : 
      21             : #ifndef MIO_SDE_SIR_MODEL_H
      22             : #define MIO_SDE_SIR_MODEL_H
      23             : 
      24             : #include "memilio/compartments/flow_model.h"
      25             : #include "memilio/epidemiology/populations.h"
      26             : #include "memilio/utils/random_number_generator.h"
      27             : #include "sde_sirs/infection_state.h"
      28             : #include "sde_sirs/parameters.h"
      29             : 
      30             : namespace mio
      31             : {
      32             : namespace ssirs
      33             : {
      34             : 
      35             : /********************
      36             :  * define the model *
      37             :  ********************/
      38             : 
      39             : using Flows = TypeList<Flow<InfectionState::Susceptible, InfectionState::Infected>,
      40             :                        Flow<InfectionState::Infected, InfectionState::Recovered>,
      41             :                        Flow<InfectionState::Recovered, InfectionState::Susceptible>>;
      42             : 
      43             : class Model : public FlowModel<ScalarType, InfectionState, Populations<ScalarType, InfectionState>, Parameters, Flows>
      44             : {
      45             :     using Base = FlowModel<ScalarType, InfectionState, mio::Populations<ScalarType, InfectionState>, Parameters, Flows>;
      46             : 
      47             : public:
      48           9 :     Model()
      49           9 :         : Base(Populations({InfectionState::Count}, 0.), ParameterSet())
      50             :     {
      51           9 :     }
      52             : 
      53          54 :     void get_flows(Eigen::Ref<const Vector<>> pop, Eigen::Ref<const Vector<>> y, ScalarType t,
      54             :                    Eigen::Ref<Vector<>> flows) const
      55             :     {
      56          54 :         auto& params         = this->parameters;
      57         108 :         ScalarType coeffStoI = params.get<ContactPatterns>().get_matrix_at(t)(0, 0) *
      58         162 :                                params.get<TransmissionProbabilityOnContact>() / populations.get_total();
      59             : 
      60          54 :         ScalarType si = mio::DistributionAdapter<std::normal_distribution<ScalarType>>::get_instance()(rng, 0.0, 1.0);
      61          54 :         ScalarType ir = mio::DistributionAdapter<std::normal_distribution<ScalarType>>::get_instance()(rng, 0.0, 1.0);
      62          54 :         ScalarType rs = mio::DistributionAdapter<std::normal_distribution<ScalarType>>::get_instance()(rng, 0.0, 1.0);
      63             : 
      64          54 :         const ScalarType inv_sqrt_dt = 1 / sqrt(step_size);
      65             : 
      66             :         // Assuming that no person can change its InfectionState twice in a single time step,
      67             :         // take the minimum of the calculated flow and the source compartment, to ensure that
      68             :         // no compartment attains negative values.
      69             : 
      70          54 :         flows[get_flat_flow_index<InfectionState::Susceptible, InfectionState::Infected>()] = std::clamp(
      71         108 :             coeffStoI * y[(size_t)InfectionState::Susceptible] * pop[(size_t)InfectionState::Infected] +
      72          54 :                 sqrt(coeffStoI * y[(size_t)InfectionState::Susceptible] * pop[(size_t)InfectionState::Infected]) *
      73          54 :                     inv_sqrt_dt * si,
      74         108 :             0.0, y[(size_t)InfectionState::Susceptible] / step_size);
      75             : 
      76          54 :         flows[get_flat_flow_index<InfectionState::Infected, InfectionState::Recovered>()] = std::clamp(
      77         108 :             (1.0 / params.get<TimeInfected>()) * y[(size_t)InfectionState::Infected] +
      78          54 :                 sqrt((1.0 / params.get<TimeInfected>()) * y[(size_t)InfectionState::Infected]) * inv_sqrt_dt * ir,
      79         108 :             0.0, y[(size_t)InfectionState::Infected] / step_size);
      80             : 
      81          54 :         flows[get_flat_flow_index<InfectionState::Recovered, InfectionState::Susceptible>()] = std::clamp(
      82         108 :             (1.0 / params.get<TimeImmune>()) * y[(size_t)InfectionState::Recovered] +
      83          54 :                 sqrt((1.0 / params.get<TimeImmune>()) * y[(size_t)InfectionState::Recovered]) * inv_sqrt_dt * rs,
      84         108 :             0.0, y[(size_t)InfectionState::Recovered] / step_size);
      85          54 :     }
      86             : 
      87             :     ScalarType step_size; ///< A step size of the model with which the stochastic process is realized.
      88             :     mutable RandomNumberGenerator rng;
      89             : 
      90             : private:
      91             : };
      92             : 
      93             : } // namespace ssirs
      94             : } // namespace mio
      95             : 
      96             : #endif // MIO_SDE_SIRS_MODEL_H