parameter_studies.h

Go to the documentation of this file.

/*
* Copyright (C) 2020-2026 MEmilio
*
* Authors: Daniel Abele, Martin J. Kuehn
*
* 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.
*/
#ifndef MIO_COMPARTMENTS_PARAMETER_STUDIES_H
#define MIO_COMPARTMENTS_PARAMETER_STUDIES_H
 
#include "memilio/io/binary_serializer.h"
#include "memilio/io/io.h"
#include "memilio/utils/logging.h"
#include "memilio/utils/miompi.h"
#include "memilio/utils/random_number_generator.h"
#include "memilio/mobility/metapopulation_mobility_instant.h"
 
#include <cassert>
#include <cmath>
#include <cstdint>
#include <type_traits>
#include <utility>
#include <vector>
 
namespace mio
{
 
template <class ParameterType, typename TimeType, typename StepType = TimeType>
class ParameterStudy
{
public:
    using Parameters = ParameterType;
    using Time       = TimeType;
    using Step       = StepType;
 
private:
    template <class CreateSimulationFunction>
        requires std::is_invocable_v<CreateSimulationFunction, Parameters, Time, Step, size_t>
    using SimulationT = std::decay_t<std::invoke_result_t<CreateSimulationFunction, Parameters, Time, Step, size_t>>;
    template <class CreateSimulationFunction, class ProcessSimulationResultFunction>
        requires std::is_invocable_v<ProcessSimulationResultFunction, SimulationT<CreateSimulationFunction>, size_t>
    using ProcessedResultT = std::decay_t<
        std::invoke_result_t<ProcessSimulationResultFunction, SimulationT<CreateSimulationFunction>, size_t>>;
    template <class CreateSimulationFunction, class ProcessSimulationResultFunction>
    using EnsembleResultT =
        std::conditional_t<std::is_void_v<ProcessedResultT<CreateSimulationFunction, ProcessSimulationResultFunction>>,
                           void,
                           std::vector<ProcessedResultT<CreateSimulationFunction, ProcessSimulationResultFunction>>>;
 
public:
    // TODO: replacement for "set_params_distributions_normal"? Maybe a special ctor for UncertainParameterSet?
 
    ParameterStudy(const Parameters& parameters, Time t0, Time tmax, Step dt, size_t num_runs)
        : m_parameters(parameters)
        , m_num_runs(num_runs)
        , m_t0{t0}
        , m_tmax{tmax}
        , m_dt(dt)
    {
    }
 
    template <class CreateSimulationFunction, class ProcessSimulationResultFunction>
    EnsembleResultT<CreateSimulationFunction, ProcessSimulationResultFunction>
    run_serial(CreateSimulationFunction&& create_simulation,
               ProcessSimulationResultFunction&& process_simulation_result)
    {
        return run_impl(0, m_num_runs, std::forward<CreateSimulationFunction>(create_simulation),
                        std::forward<ProcessSimulationResultFunction>(process_simulation_result));
    }
 
    template <class CreateSimulationFunction>
    std::vector<SimulationT<CreateSimulationFunction>> run_serial(CreateSimulationFunction&& create_simulation)
    {
        return run_serial(
            std::forward<CreateSimulationFunction>(create_simulation),
            [](SimulationT<CreateSimulationFunction>&& sim, size_t) -> SimulationT<CreateSimulationFunction>&& {
                return std::move(sim);
            });
    }
    template <class CreateSimulationFunction, class ProcessSimulationResultFunction>
    EnsembleResultT<CreateSimulationFunction, ProcessSimulationResultFunction>
    run(CreateSimulationFunction&& create_simulation, ProcessSimulationResultFunction&& process_simulation_result)
    {
        using ResultT = EnsembleResultT<CreateSimulationFunction, ProcessSimulationResultFunction>;
        int num_procs = 1, rank = 0;
 
#ifdef MEMILIO_ENABLE_MPI
        MPI_Comm_size(mpi::get_world(), &num_procs);
        MPI_Comm_rank(mpi::get_world(), &rank);
#endif
 
        //The ParameterDistributions used for sampling parameters use thread_local_rng()
        //So we set our own RNG to be used.
        //Assume that sampling uses the thread_local_rng() and isn't multithreaded
        m_rng.synchronize();
 
        std::vector<size_t> run_distribution = distribute_runs(m_num_runs, num_procs);
        size_t start_run_idx = std::accumulate(run_distribution.begin(), run_distribution.begin() + rank, size_t(0));
        size_t end_run_idx   = start_run_idx + run_distribution[rank];
 
        if constexpr (std::is_void_v<ResultT>) {
            // if the processor returns nothing, there is nothing to synchronize
            run_impl(start_run_idx, end_run_idx, std::forward<CreateSimulationFunction>(create_simulation),
                     std::forward<ProcessSimulationResultFunction>(process_simulation_result));
            return;
        }
        else {
            auto ensemble_result =
                run_impl(start_run_idx, end_run_idx, std::forward<CreateSimulationFunction>(create_simulation),
                         std::forward<ProcessSimulationResultFunction>(process_simulation_result));
 
#ifdef MEMILIO_ENABLE_MPI
            //gather results
            if (rank == 0) {
                for (int src_rank = 1; src_rank < num_procs; ++src_rank) {
                    int bytes_size;
                    MPI_Recv(&bytes_size, 1, MPI_INT, src_rank, 0, mpi::get_world(), MPI_STATUS_IGNORE);
                    ByteStream bytes(bytes_size);
                    MPI_Recv(bytes.data(), bytes.data_size(), MPI_BYTE, src_rank, 0, mpi::get_world(),
                             MPI_STATUS_IGNORE);
 
                    IOResult<ResultT> src_ensemble_results = deserialize_binary(bytes, Tag<ResultT>{});
                    if (!src_ensemble_results) {
                        log_error("Error receiving ensemble results from rank {}.", src_rank);
                    }
                    std::copy(src_ensemble_results.value().begin(), src_ensemble_results.value().end(),
                              std::back_inserter(ensemble_result));
                }
            }
            else {
                ByteStream bytes = serialize_binary(ensemble_result);
                int bytes_size   = int(bytes.data_size());
                MPI_Send(&bytes_size, 1, MPI_INT, 0, 0, mpi::get_world());
                MPI_Send(bytes.data(), bytes.data_size(), MPI_BYTE, 0, 0, mpi::get_world());
                ensemble_result.clear(); //only return root process
            }
#endif
 
            return ensemble_result;
        }
    }
 
    size_t get_num_runs() const
    {
        return m_num_runs;
    }
 
    Time get_tmax() const
    {
        return m_tmax;
    }
 
    Time get_t0() const
    {
        return m_t0;
    }
 
    Time get_dt() const
    {
        return m_dt;
    }
 
    const Parameters& get_parameters() const
    {
        return m_parameters;
    }
    Parameters& get_parameters()
    {
        return m_parameters;
    }
    RandomNumberGenerator& get_rng()
    {
        return m_rng;
    }
 
private:
    template <class T>
    inline auto make_ensemble_result()
    {
        if constexpr (std::is_void_v<T>) {
            return char(0); // placeholder, as we cannot instanciate void
        }
        else {
            T result;
            result.reserve(m_num_runs);
            return result;
        }
    }
 
    template <class CreateSimulationFunction, class ProcessSimulationResultFunction>
    EnsembleResultT<CreateSimulationFunction, ProcessSimulationResultFunction>
    run_impl(size_t start_run_idx, size_t end_run_idx, CreateSimulationFunction&& create_simulation,
             ProcessSimulationResultFunction&& process_simulation_result)
    {
        using ResultT = EnsembleResultT<CreateSimulationFunction, ProcessSimulationResultFunction>;
        assert(start_run_idx <= end_run_idx);
        // this bool (and all code that it is used in) enables using void functions
        constexpr bool should_gather_results = !std::is_void_v<ResultT>;
 
        // Note that this overwrites seed and counter of thread_local_rng, but it does not replace it.
        thread_local_rng() = m_rng;
 
        // the result is not used or returned if the processing function returns nothing (i.e. void)
        [[maybe_unused]] auto ensemble_result = make_ensemble_result<ResultT>();
 
        for (size_t run_idx = start_run_idx; run_idx < end_run_idx; run_idx++) {
            log(LogLevel::info, "ParameterStudies: run {}", run_idx);
 
            //prepare rng for this run by setting the counter to the right offset
            //Add the old counter so that this call of run() produces different results
            //from the previous call
            Counter<uint64_t> run_rng_counter =
                m_rng.get_counter() +
                rng_totalsequence_counter<uint64_t>(static_cast<uint32_t>(run_idx), Counter<uint32_t>(0));
            thread_local_rng().set_counter(run_rng_counter);
 
            //sample
            SimulationT<CreateSimulationFunction> sim =
                create_simulation(std::as_const(m_parameters), std::as_const(m_t0), std::as_const(m_dt), run_idx);
 
            // the create_counter is only used for debug logging, hence it is unused in Release builds
            [[maybe_unused]] const uint64_t create_counter = (thread_local_rng().get_counter() - run_rng_counter).get();
            log_debug("ParameterStudy: Generated {} random numbers creating simulation #{}.", create_counter, run_idx);
 
            //perform run
            sim.advance(m_tmax);
 
            log_debug("ParameterStudy: Generated {} random numbers running simulation #{}.",
                      run_rng_counter.get() - create_counter, run_idx);
 
            //handle result
            if constexpr (should_gather_results) {
                ensemble_result.emplace_back(process_simulation_result(std::move(sim), run_idx));
            }
            else {
                process_simulation_result(std::move(sim), run_idx);
            }
        }
 
        //Set the counter of our RNG so that future calls of run() produce different parameters.
        m_rng.set_counter(m_rng.get_counter() + rng_totalsequence_counter<uint64_t>(m_num_runs, Counter<uint32_t>(0)));
 
        if constexpr (should_gather_results) {
            return ensemble_result;
        }
        // else: there is nothing to return, because process_simulation_result returns void
    }
 
    static std::vector<size_t> distribute_runs(size_t num_runs, int num_procs)
    {
        assert(num_procs > 0);
        //evenly distribute runs
        //lower processes do one more run if runs are not evenly distributable
        size_t num_runs_local = num_runs / num_procs; //integer division!
        size_t remainder      = num_runs % num_procs;
 
        std::vector<size_t> run_distribution(num_procs);
        std::fill(run_distribution.begin(), run_distribution.begin() + remainder, num_runs_local + 1);
        std::fill(run_distribution.begin() + remainder, run_distribution.end(), num_runs_local);
 
        return run_distribution;
    }
 
    ParameterType m_parameters; 
    size_t m_num_runs; 
    Time m_t0, m_tmax; 
    Step m_dt; 
    RandomNumberGenerator m_rng; 
};
 
template <typename FP, class Sim>
auto make_sampled_graph_simulation(const Graph<typename Sim::Model, MobilityParameters<FP>>& sampled_graph, FP t0,
                                   FP dt_node_sim, FP dt_graph_sim)
{
    using SimGraph = Graph<SimulationNode<FP, Sim>, MobilityEdge<FP>>;
 
    SimGraph sim_graph;
 
    for (auto&& node : sampled_graph.nodes()) {
        sim_graph.add_node(node.id, node.property, t0, dt_node_sim);
    }
    for (auto&& edge : sampled_graph.edges()) {
        sim_graph.add_edge(edge.start_node_idx, edge.end_node_idx, edge.property);
    }
 
    return make_mobility_sim<FP, Sim>(t0, dt_graph_sim, std::move(sim_graph));
}
 
} // namespace mio
 
#endif // MIO_COMPARTMENTS_PARAMETER_STUDIES_H