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_D_ABM_SIMULATION_H
      22             : #define MIO_D_ABM_SIMULATION_H
      23             : 
      24             : #include "d_abm/model.h"
      25             : #include "memilio/config.h"
      26             : #include "memilio/compartments/simulation.h"
      27             : #include "memilio/utils/random_number_generator.h"
      28             : 
      29             : namespace mio
      30             : {
      31             : 
      32             : namespace dabm
      33             : {
      34             : 
      35             : /**
      36             :  * @brief A specialized Simulation for mio::dabm::Model.
      37             :  * @tparam Implementation A class implementing all functions and types marked with the using keyword in Model, see dabm::Model.
      38             :  */
      39             : template <class Implementation>
      40             : class Simulation
      41             : {
      42             :     using Status = typename Implementation::Status;
      43             : 
      44             : public:
      45             :     using Model = dabm::Model<Implementation>;
      46             : 
      47             :     /**
      48             :      * @brief Set up the simulation for a diffusive ABM.
      49             :      * @param[in] model An instance of mio::dabm::Model.
      50             :      * @param[in] t0 Start time.
      51             :      * @param[in] dt Step size of integration.
      52             :      */
      53           7 :     Simulation(const Model& model, ScalarType t0 = 0, ScalarType dt = 0.1)
      54           7 :         : m_t0(t0)
      55           7 :         , m_dt(dt)
      56           7 :         , m_model(std::make_unique<Model>(model))
      57           7 :         , m_result(t0, m_model->time_point())
      58             :     {
      59           7 :         assert(dt > 0);
      60           7 :         m_current_rates.reserve(m_model->populations.size());
      61           7 :         m_current_events.reserve(m_model->populations.size());
      62           7 :     }
      63             : 
      64             :     /**
      65             :      * @brief Advance simulation to tmax.
      66             :      * This function performs a temporal Gillespie.
      67             :      * @param tmax Next stopping point of simulation.
      68             :      */
      69           7 :     void advance(const ScalarType t_max)
      70             :     {
      71             :         // draw time until an adoption takes place
      72           7 :         ScalarType waiting_time = mio::ExponentialDistribution<ScalarType>::get_instance()(m_model->get_rng(), 1.0);
      73        2107 :         while (m_t0 < t_max) {
      74        2100 :             ScalarType dt             = std::min({m_dt, t_max - m_t0});
      75        2100 :             ScalarType remaining_time = dt;
      76        2107 :             while (remaining_time > 0) {
      77        2107 :                 compute_current_rates_and_events(); // lambda_k (aka f-hat(N))
      78             :                 ScalarType cumulative_adoption_rate =
      79        2107 :                     std::accumulate(m_current_rates.begin(), m_current_rates.end(), 0.0); // Lambda
      80             :                 // status update
      81        2107 :                 if (waiting_time < cumulative_adoption_rate * remaining_time) {
      82             :                     // draw which adoption takes place
      83             :                     const size_t event_id =
      84           7 :                         mio::DiscreteDistribution<size_t>::get_instance()(m_model->get_rng(), m_current_rates);
      85           7 :                     Event& event = m_current_events[event_id];
      86             :                     // perform adoption
      87           7 :                     m_model->adopt(event.agent, event.new_status);
      88             :                     // draw new waiting time
      89           7 :                     remaining_time -= waiting_time / cumulative_adoption_rate;
      90           7 :                     waiting_time = mio::ExponentialDistribution<ScalarType>::get_instance()(m_model->get_rng(), 1.0);
      91             :                 }
      92             :                 else {
      93             :                     // no event, decrease waiting time
      94        2100 :                     waiting_time -= cumulative_adoption_rate * remaining_time;
      95        2100 :                     break;
      96             :                 }
      97             :             }
      98             :             // position update
      99        8400 :             for (auto& agent : m_model->populations) {
     100        6300 :                 m_model->move(m_t0, dt, agent);
     101             :             }
     102        2100 :             m_t0 += dt;
     103             :             // store result
     104        2100 :             m_result.add_time_point(m_t0, m_model->time_point());
     105             :         }
     106           7 :     }
     107             : 
     108             :     /**
     109             :      * @brief Returns the final simulation result.
     110             :      * @return A TimeSeries to represent the final simulation result.
     111             :      */
     112          70 :     TimeSeries<ScalarType>& get_result()
     113             :     {
     114          70 :         return m_result;
     115             :     }
     116             :     const TimeSeries<ScalarType>& get_result() const
     117             :     {
     118             :         return m_result;
     119             :     }
     120             : 
     121             :     /**
     122             :      * @brief Returns the model used in the simulation.
     123             :      */
     124             :     const Model& get_model() const
     125             :     {
     126             :         return *m_model;
     127             :     }
     128             :     Model& get_model()
     129             :     {
     130             :         return *m_model;
     131             :     }
     132             : 
     133             : private:
     134             :     /**
     135             :      * @brief Struct defining an adoption event for an agent and target infection state.
     136             :      */
     137             :     struct Event {
     138             :         typename Model::Agent& agent;
     139             :         Status new_status;
     140             :     };
     141             : 
     142             :     /**
     143             :      * @brief Calculate current values for m_current_rates and m_current_events.
     144             :      */
     145        2107 :     inline void compute_current_rates_and_events()
     146             :     {
     147        2107 :         m_current_rates.clear();
     148        2107 :         m_current_events.clear();
     149             :         // compute rate for each (agent, status) combination
     150        8428 :         for (auto& agent : m_model->populations) {
     151       44247 :             for (int s = 0; s < static_cast<int>(Status::Count); s++) {
     152       37926 :                 Status new_status = static_cast<Status>(s);
     153             :                 // check if an adoption from the agents status is possible
     154       37926 :                 auto adoption_rate = m_model->adoption_rate(agent, new_status);
     155       37926 :                 if (adoption_rate > 0) {
     156             :                     // add rate and corresponding event
     157          42 :                     m_current_rates.push_back(adoption_rate);
     158          42 :                     m_current_events.push_back({agent, new_status});
     159             :                 }
     160             :             }
     161             :         }
     162        2107 :     }
     163             : 
     164             :     ScalarType m_t0, m_dt; ///< Start time of the simulation and integration step size.
     165             :     std::unique_ptr<Model> m_model; ///< Pointer to the model used in the simulation.
     166             :     std::vector<ScalarType> m_current_rates; ///< Current adoption rates.
     167             :     std::vector<Event> m_current_events; ///< Contains an event corresponding to each rate in m_current_rates.
     168             :     mio::TimeSeries<ScalarType> m_result; ///< Result time series.
     169             : };
     170             : } //namespace dabm
     171             : } // namespace mio
     172             : 
     173             : #endif