Line data    Source code

       1             : /* 
       2             : * Copyright (C) 2020-2025 German Aerospace Center (DLR-SC)
       3             : *
       4             : * Authors: René Schmieding, Julia Bicker
       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_SMM_MODEL_H
      22             : #define MIO_SMM_MODEL_H
      23             : 
      24             : #include "memilio/config.h"
      25             : #include "smm/parameters.h"
      26             : #include "memilio/compartments/compartmentalmodel.h"
      27             : #include "memilio/epidemiology/populations.h"
      28             : #include "memilio/geography/regions.h"
      29             : 
      30             : namespace mio
      31             : {
      32             : namespace smm
      33             : {
      34             : 
      35             : /**
      36             :  * @brief Stochastic Metapopulation Model.
      37             :  * @tparam regions Number of regions.
      38             :  * @tparam Status An infection state enum.
      39             :  */
      40             : template <size_t regions, class Status>
      41             : class Model
      42             :     : public mio::CompartmentalModel<ScalarType, Status, mio::Populations<ScalarType, mio::regions::Region, Status>,
      43             :                                      ParametersBase<Status>>
      44             : {
      45             :     using Base = mio::CompartmentalModel<ScalarType, Status, mio::Populations<ScalarType, mio::regions::Region, Status>,
      46             :                                          ParametersBase<Status>>;
      47             : 
      48             : public:
      49         104 :     Model()
      50         208 :         : Base(typename Base::Populations({static_cast<mio::regions::Region>(regions), Status::Count}, 0.),
      51         312 :                typename Base::ParameterSet())
      52             :     {
      53         104 :     }
      54             : 
      55             :     /**
      56             :      * @brief Calculate the current rate of the given adoption.
      57             :      * @param[in] rate An adoption rate from this model.
      58             :      * @param[in] x The current state of the model.
      59             :      * @return Current value of the adoption rate.
      60             :      */
      61          21 :     ScalarType evaluate(const AdoptionRate<Status>& rate, const Eigen::VectorXd& x) const
      62             :     {
      63          21 :         const auto& pop   = this->populations;
      64          21 :         const auto source = pop.get_flat_index({rate.region, rate.from});
      65             :         // determine order and calculate rate
      66          21 :         if (rate.influences.size() == 0) { // first order adoption
      67          16 :             return rate.factor * x[source];
      68             :         }
      69             :         else { // second order adoption
      70           5 :             ScalarType N = 0;
      71          35 :             for (size_t s = 0; s < static_cast<size_t>(Status::Count); ++s) {
      72          30 :                 N += x[pop.get_flat_index({rate.region, Status(s)})];
      73             :             }
      74             :             // accumulate influences
      75           5 :             ScalarType influences = 0.0;
      76          15 :             for (size_t i = 0; i < rate.influences.size(); i++) {
      77          10 :                 influences +=
      78          10 :                     rate.influences[i].factor * x[pop.get_flat_index({rate.region, rate.influences[i].status})];
      79             :             }
      80           5 :             return (N > 0) ? (rate.factor * x[source] * influences / N) : 0;
      81             :         }
      82             :     }
      83             : 
      84             :     /**
      85             :      * @brief Calculate the current rate of the given spatial transition.
      86             :      * @param[in] rate A transition rate from this model.
      87             :      * @param[in] x The current state of the model.
      88             :      * @return Current value of the transition rate.
      89             :      */
      90       18206 :     ScalarType evaluate(const TransitionRate<Status>& rate, const Eigen::VectorXd& x) const
      91             :     {
      92       18206 :         const auto source = this->populations.get_flat_index({rate.from, rate.status});
      93       18206 :         return rate.factor * x[source];
      94             :     }
      95             : 
      96             :     /**
      97             :     * Get the RandomNumberGenerator used by this Model for random events.
      98             :     * @return The random number generator.
      99             :     */
     100       18211 :     mio::RandomNumberGenerator& get_rng()
     101             :     {
     102       18211 :         return m_rng;
     103             :     }
     104             : 
     105             : private:
     106             :     mio::RandomNumberGenerator m_rng; ///< Model's random number generator.
     107             : };
     108             : 
     109             : } //namespace smm
     110             : 
     111             : } // namespace mio
     112             : 
     113             : #endif