random_events.h

Go to the documentation of this file.

/* 
* Copyright (C) 2020-2026 MEmilio
*
* Authors: Daniel Abele
*
* Contact: Martin J. Kuehn <Martin.Kuehn@DLR.de>
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
*     http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
#include "abm/time.h"
#include "memilio/utils/random_number_generator.h"
#include <algorithm>
#include <array>
#include <numeric>
 
namespace mio
{
namespace abm
{
 
template <class RNG, class T, size_t NumTransitions>
T random_transition(RNG& rng, T current_state, TimeSpan dt,
                    const std::pair<T, ScalarType> (&transitions)[NumTransitions])
{
    assert(std::all_of(std::begin(transitions), std::end(transitions),
                       [](auto& p) {
                           return p.second >= 0.0;
                       }) &&
           "transition rates must be non-negative");
 
    //check if any transition happens using exponential distribution with the sum of all transition rates
    auto sum = std::accumulate(std::begin(transitions), std::end(transitions), 0.0, [](auto&& a, auto&& t) {
        return a + t.second;
    });
    if (sum <= 0) { //no transitions or all transitions have rate zero
        return current_state;
    }
    auto v = ExponentialDistribution<ScalarType>::get_instance()(rng, sum);
    if (v < dt.days()) {
        //pick one of the possible transitions using discrete distribution
        std::array<ScalarType, NumTransitions> rates;
        std::transform(std::begin(transitions), std::end(transitions), rates.begin(), [](auto&& t) {
            return t.second;
        });
        auto random_idx = DiscreteDistribution<size_t>::get_instance()(rng, rates);
        return transitions[random_idx].first;
    }
 
    return current_state;
}
 
} // namespace abm
} // namespace mio