Line data    Source code

       1             : /* 
       2             : * Copyright (C) 2020-2025 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 Eigen::VectorX<ScalarType>> pop, Eigen::Ref<const Eigen::VectorX<ScalarType>> y,
      53             :                    ScalarType t, Eigen::Ref<Eigen::VectorX<ScalarType>> flows) const
      54             :     {
      55          54 :         auto& params         = this->parameters;
      56         108 :         ScalarType coeffStoI = params.get<ContactPatterns>().get_matrix_at(t)(0, 0) *
      57         162 :                                params.get<TransmissionProbabilityOnContact>() / populations.get_total();
      58             : 
      59          54 :         ScalarType si = mio::DistributionAdapter<std::normal_distribution<ScalarType>>::get_instance()(rng, 0.0, 1.0);
      60          54 :         ScalarType ir = mio::DistributionAdapter<std::normal_distribution<ScalarType>>::get_instance()(rng, 0.0, 1.0);
      61             : 
      62             :         // Assuming that no person can change its InfectionState twice in a single time step,
      63             :         // take the minimum of the calculated flow and the source compartment, to ensure that
      64             :         // no compartment attains negative values.
      65             : 
      66          54 :         flows[get_flat_flow_index<InfectionState::Susceptible, InfectionState::Infected>()] = std::clamp(
      67         108 :             coeffStoI * y[(size_t)InfectionState::Susceptible] * pop[(size_t)InfectionState::Infected] +
      68          54 :                 sqrt(coeffStoI * y[(size_t)InfectionState::Susceptible] * pop[(size_t)InfectionState::Infected]) /
      69          54 :                     sqrt(step_size) * si,
      70         108 :             0.0, y[(size_t)InfectionState::Susceptible] / step_size);
      71             : 
      72          54 :         flows[get_flat_flow_index<InfectionState::Infected, InfectionState::Recovered>()] = std::clamp(
      73         108 :             (1.0 / params.get<TimeInfected>()) * y[(size_t)InfectionState::Infected] +
      74          54 :                 sqrt((1.0 / params.get<TimeInfected>()) * y[(size_t)InfectionState::Infected]) / sqrt(step_size) * ir,
      75         108 :             0.0, y[(size_t)InfectionState::Infected] / step_size);
      76          54 :     }
      77             : 
      78             :     ScalarType step_size; ///< A step size of the model with which the stochastic process is realized.
      79             :     mutable RandomNumberGenerator rng;
      80             : 
      81             : private:
      82             : };
      83             : 
      84             : } // namespace ssir
      85             : } // namespace mio
      86             : 
      87             : #endif // MIO_SDE_SIR_MODEL_H