Line data    Source code

       1             : /* 
       2             : * Copyright (C) 2020-2024 MEmilio
       3             : *
       4             : * Authors: Daniel Abele
       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             : #include "abm/time.h"
      21             : #include "memilio/utils/random_number_generator.h"
      22             : #include <algorithm>
      23             : #include <array>
      24             : #include <numeric>
      25             : 
      26             : namespace mio
      27             : {
      28             : namespace abm
      29             : {
      30             : 
      31             : /**
      32             :  * @brief Select a random transition from a list of possible transitions from the current state to others.
      33             :  * Each transition is represented by the new state and the probability of the transition, e.g. 
      34             :  * a pair {1, 0.5} is the transition to state 1 with rate 0.5.
      35             :  * Transition rates are not probabilities but the parameters of an exponential distribution.
      36             :  * One of the transitions happens if x < dt, where x is a sample from the exponential distribution Exp(S), 
      37             :  * S begin the sum of all rates. Which transition happens is determined by sampling from a discrete distribution
      38             :  * with the rates as weights. It's also possible that no transition happens in this time step.
      39             :  * In this case the current state is returned.
      40             :  * @tparam RNG Type that satisfies the UniformRandomBitGenerator concept.
      41             :  * @tparam T Type that represents the states.
      42             :  * @tparam NumTransitions Number of possible transitions.
      43             :  * @param[inout] rng RandomNumberGenerator.
      44             :  * @param[in] current_state Current state before transition.
      45             :  * @param[in] dt Length of the time step.
      46             :  * @param[in] transitions Array of pairs of new states and their rates (probabilities).
      47             :  * @return New state from the list if transition happens, current_state otherwise.
      48             :  */
      49             : template <class RNG, class T, size_t NumTransitions>
      50        2367 : T random_transition(RNG& rng, T current_state, TimeSpan dt, const std::pair<T, double> (&transitions)[NumTransitions])
      51             : {
      52        4806 :     assert(std::all_of(std::begin(transitions), std::end(transitions),
      53             :                        [](auto& p) {
      54             :                            return p.second >= 0.0;
      55             :                        }) &&
      56             :            "transition rates must be non-negative");
      57             : 
      58             :     //check if any transition happens using exponential distribution with the sum of all transition rates
      59        2826 :     auto sum = std::accumulate(std::begin(transitions), std::end(transitions), 0.0, [](auto&& a, auto&& t) {
      60        2439 :         return a + t.second;
      61             :     });
      62        2367 :     if (sum <= 0) { //no transitions or all transitions have rate zero
      63        1854 :         return current_state;
      64             :     }
      65         513 :     auto v = ExponentialDistribution<double>::get_instance()(rng, sum);
      66         513 :     if (v < dt.days()) {
      67             :         //pick one of the possible transitions using discrete distribution
      68          81 :         std::array<double, NumTransitions> rates;
      69         171 :         std::transform(std::begin(transitions), std::end(transitions), rates.begin(), [](auto&& t) {
      70         117 :             return t.second;
      71             :         });
      72          81 :         auto random_idx = DiscreteDistribution<size_t>::get_instance()(rng, rates);
      73          81 :         return transitions[random_idx].first;
      74             :     }
      75             : 
      76         432 :     return current_state;
      77             : }
      78             : 
      79             : } // namespace abm
      80             : } // namespace mio