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_sir/infection_state.h"
      28             : #include "sde_sir/parameters.h"
      29             : 
      30             : namespace mio
      31             : {
      32             : namespace ssir
      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             : 
      42             : class Model : public FlowModel<ScalarType, InfectionState, Populations<ScalarType, InfectionState>, Parameters, Flows>
      43             : {
      44             :     using Base = FlowModel<ScalarType, InfectionState, mio::Populations<ScalarType, InfectionState>, Parameters, Flows>;
      45             : 
      46             : public:
      47           9 :     Model()
      48           9 :         : Base(Populations({InfectionState::Count}, 0.), ParameterSet())
      49             :     {
      50           9 :     }
      51             : 
      52          54 :     void get_flows(Eigen::Ref<const Vector<>> pop,
      53             :                    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             : 
      63             :         // Assuming that no person can change its InfectionState twice in a single time step,
      64             :         // take the minimum of the calculated flow and the source compartment, to ensure that
      65             :         // no compartment attains negative values.
      66             : 
      67          54 :         flows[get_flat_flow_index<InfectionState::Susceptible, InfectionState::Infected>()] = std::clamp(
      68         108 :             coeffStoI * y[(size_t)InfectionState::Susceptible] * pop[(size_t)InfectionState::Infected] +
      69          54 :                 sqrt(coeffStoI * y[(size_t)InfectionState::Susceptible] * pop[(size_t)InfectionState::Infected]) /
      70          54 :                     sqrt(step_size) * si,
      71         108 :             0.0, y[(size_t)InfectionState::Susceptible] / step_size);
      72             : 
      73          54 :         flows[get_flat_flow_index<InfectionState::Infected, InfectionState::Recovered>()] = std::clamp(
      74         108 :             (1.0 / params.get<TimeInfected>()) * y[(size_t)InfectionState::Infected] +
      75          54 :                 sqrt((1.0 / params.get<TimeInfected>()) * y[(size_t)InfectionState::Infected]) / sqrt(step_size) * ir,
      76         108 :             0.0, y[(size_t)InfectionState::Infected] / step_size);
      77          54 :     }
      78             : 
      79             :     ScalarType step_size; ///< A step size of the model with which the stochastic process is realized.
      80             :     mutable RandomNumberGenerator rng;
      81             : 
      82             : private:
      83             : };
      84             : 
      85             : } // namespace ssir
      86             : } // namespace mio
      87             : 
      88             : #endif // MIO_SDE_SIR_MODEL_H