simulation.hpp

 
// File: simulation.hpp                                                       //
// Project: respond                                                           //
// Created Date: 2026-02-05                                                   //
// Author: Matthew Carroll                                                    //
// -----                                                                      //
// Last Modified: 2026-10-06                                                  //
// Modified By: Matthew Carroll                                               //
// -----                                                                      //
// Copyright (c) 2026 Syndemics Lab at Boston Medical Center                  //
#ifndef RESPOND_SIMULATION_HPP_
#define RESPOND_SIMULATION_HPP_
 
#include <respond/constants.hpp>
#include <respond/eigen_config.hpp>
#include <respond/history.hpp>
#include <respond/logging.hpp>
#include <respond/model.hpp>
#include <respond/runtime_config.hpp>
 
#include <atomic>
#include <exception>
#include <map>
#include <memory>
#include <mutex>
#include <stdexcept>
#include <string>
#include <thread>
#include <utility>
#include <vector>
 
#include <Eigen/Dense>
 
namespace respond {
class Simulation {
public:
    class ModelSlotProxy {
    public:
        ModelSlotProxy(Simulation &owner, size_t idx)
            : _owner(&owner), _idx(idx) {}
 
        Model *operator->() { return &_owner->GetModelRefOrThrow(_idx); }
 
        Model &operator*() { return _owner->GetModelRefOrThrow(_idx); }
 
        operator Model &() { return _owner->GetModelRefOrThrow(_idx); }
 
        ModelSlotProxy &operator=(const ModelSlotProxy &other) {
            return *this = static_cast<const Model &>(
                       other._owner->GetModelRefOrThrow(other._idx));
        }
 
        ModelSlotProxy &operator=(const Model &model) {
            _owner->GetModelRefOrThrow(_idx);
            _owner->_models[_idx] = model.clone();
            return *this;
        }
 
        ModelSlotProxy &operator=(const std::unique_ptr<Model> &model) {
            if (!model) {
                LogError(_owner->_runtime_config.logging.logger_name,
                         "Cannot assign null model pointer to simulation "
                         "slot.");
                throw std::invalid_argument(
                    "Error attempting to assign null model pointer.");
            }
            _owner->GetModelRefOrThrow(_idx);
            _owner->_models[_idx] = model->clone();
            return *this;
        }
 
    private:
        Simulation *_owner;
        size_t _idx;
    };
 
    //
    // Rule of Five: Copy and Move Semantics
    //
 
    Simulation() : Simulation(RuntimeConfig{}) {}
 
    [[deprecated("Use Simulation(RuntimeConfig) instead")]]
    Simulation(const std::string &log_name)
        : Simulation(RuntimeConfig{
              ExecutionConfig{},
              LoggingConfig{log_name, RESPOND_DEFAULT_LOG_FILE, false}}) {}
 
    [[deprecated("Use Simulation(RuntimeConfig) instead")]]
    Simulation(const std::string &log_name, const std::string &log_filepath)
        : Simulation(
              RuntimeConfig{ExecutionConfig{},
                            LoggingConfig{log_name, log_filepath, false}}) {}
 
    [[deprecated("Use Simulation(RuntimeConfig) instead")]]
    Simulation(const std::string &log_name, const std::string &log_filepath,
               const ExecutionConfig &execution_config)
        : Simulation(RuntimeConfig{
              execution_config, LoggingConfig{log_name, log_filepath, false}}) {
    }
 
    explicit Simulation(const RuntimeConfig &runtime_config)
        : _runtime_config(runtime_config) {
        if (ConfigureLogger(_runtime_config.logging) ==
            CreationStatus::kError) {
            throw std::runtime_error(
                "Error attempting to initialize simulation logger.");
        }
    }
 
    ~Simulation() = default;
 
    Simulation(const Simulation &other) {
        _runtime_config = other._runtime_config;
        for (const auto &m : other._models) {
            _models.push_back(m->clone());
        }
        _duration = other._duration;
        _parameter_change_times = other._parameter_change_times;
        _stratify_entering_cohort = other._stratify_entering_cohort;
        _build_summary_stats = other._build_summary_stats;
        _save_state_history = other._save_state_history;
        _timesteps_to_report = other._timesteps_to_report;
        _pivot_long = other._pivot_long;
    }
 
    Simulation &operator=(const Simulation &other) {
        if (this != &other) {
            std::vector<std::unique_ptr<Model>> models;
            for (const auto &m : other._models) {
                models.push_back(m->clone());
            }
            _runtime_config = other._runtime_config;
            _models = std::move(models);
            _duration = other._duration;
            _parameter_change_times = other._parameter_change_times;
            _stratify_entering_cohort = other._stratify_entering_cohort;
            _build_summary_stats = other._build_summary_stats;
            _save_state_history = other._save_state_history;
            _timesteps_to_report = other._timesteps_to_report;
            _pivot_long = other._pivot_long;
        }
        return *this;
    }
 
    Simulation(Simulation &&other) noexcept
        : _runtime_config(std::move(other._runtime_config)),
          _models(std::move(other._models)), _duration(other._duration),
          _parameter_change_times(std::move(other._parameter_change_times)),
          _stratify_entering_cohort(other._stratify_entering_cohort),
          _build_summary_stats(other._build_summary_stats),
          _save_state_history(other._save_state_history),
          _timesteps_to_report(std::move(other._timesteps_to_report)),
          _pivot_long(other._pivot_long) {}
 
    Simulation &operator=(Simulation &&other) noexcept {
        if (this != &other) {
            _runtime_config = std::move(other._runtime_config);
            _models = std::move(other._models);
            _duration = other._duration;
            _parameter_change_times = std::move(other._parameter_change_times);
            _stratify_entering_cohort = other._stratify_entering_cohort;
            _build_summary_stats = other._build_summary_stats;
            _save_state_history = other._save_state_history;
            _timesteps_to_report = std::move(other._timesteps_to_report);
            _pivot_long = other._pivot_long;
        }
        return *this;
    }
 
    //
    // Simulation Behavior Methods: Model Management
    //
 
    std::unique_ptr<Model> CreateNewModel(const std::string &model_name) {
        _models.push_back(Model::Create(model_name, _runtime_config));
        return _models.back()->clone();
    }
 
    void ClearModels() { _models.clear(); }
 
    void AddModel(const std::unique_ptr<Model> &model) {
        if (!model) {
            LogError(_runtime_config.logging.logger_name,
                     "Cannot add a null model to the simulation.");
            throw std::invalid_argument(
                "Error attempting to add a null model to simulation.");
        }
        _models.push_back(model->clone());
    }
 
    void Run(int duration = -1) {
        if (duration == 0 || duration < -1) {
            throw std::invalid_argument(
                "Simulation duration must be positive or -1.");
        }
        if (_models.empty()) {
            LogError(_runtime_config.logging.logger_name,
                     "Cannot run a simulation with no models.");
            throw std::invalid_argument(
                "Error attempting to run simulation with no models.");
        }
        std::lock_guard<std::mutex> execution_lock(
            detail::GetEigenExecutionMutex());
        if (duration > 0) {
            _duration = duration;
        }
        const auto &execution = _runtime_config.execution;
        const unsigned int thread_limit = std::thread::hardware_concurrency();
        const unsigned int eigen_threads =
            thread_limit == 0 ? execution.eigen_threads
                              : std::min(execution.eigen_threads, thread_limit);
        Eigen::setNbThreads(eigen_threads);
 
        LogInfo(_runtime_config.logging.logger_name,
                "Running simulation for duration of " +
                    std::to_string(_duration) + " timesteps.");
        const auto run_model = [this](const std::unique_ptr<Model> &model) {
            model->SetFinalTimestep(_duration);
            model->RunTimesteps();
        };
        if (!execution.run_models_concurrently || _models.size() < 2) {
            for (const auto &model : _models) {
                run_model(model);
            }
            return;
        }
 
        if (execution.eigen_threads > 1) {
            throw std::invalid_argument(
                "Concurrent model execution requires eigen_threads == 1.");
        }
 
        unsigned int worker_limit = execution.total_threads;
        if (worker_limit == 0) {
            worker_limit = std::thread::hardware_concurrency();
            if (worker_limit == 0) {
                worker_limit = 1;
            }
        }
        const auto worker_count =
            std::min<size_t>(worker_limit, _models.size());
        if (worker_count <= 1) {
            for (const auto &model : _models) {
                run_model(model);
            }
            return;
        }
 
        std::atomic<size_t> next_model{0};
        std::exception_ptr first_exception;
        std::mutex exception_mutex;
        std::vector<std::thread> workers;
        workers.reserve(worker_count);
 
        for (size_t worker = 0; worker < worker_count; ++worker) {
            workers.emplace_back([&]() {
                while (true) {
                    const size_t index = next_model.fetch_add(1);
                    if (index >= _models.size()) {
                        return;
                    }
                    try {
                        run_model(_models[index]);
                    } catch (...) {
                        std::lock_guard<std::mutex> lock(exception_mutex);
                        if (!first_exception) {
                            first_exception = std::current_exception();
                        }
                        return;
                    }
                }
            });
        }
 
        for (auto &worker : workers) {
            worker.join();
        }
        if (first_exception) {
            std::rethrow_exception(first_exception);
        }
    }
 
    //
    // Getters and Setters for Transitions and Metadata
    //
 
    std::vector<std::unique_ptr<Model>> GetModels() const {
        std::vector<std::unique_ptr<Model>> _models_copy;
        for (const auto &model : _models) {
            _models_copy.push_back(model->clone());
        }
        return _models_copy;
    }
 
    ModelSlotProxy operator[](size_t idx) {
        GetModelRefOrThrow(idx);
        return ModelSlotProxy(*this, idx);
    }
 
    const Model &operator[](size_t idx) const {
        return GetModelRefOrThrow(idx);
    }
 
    std::map<size_t, std::string> GetModelIndexNameMap() const {
        std::map<size_t, std::string> ret;
        for (size_t i = 0; i < _models.size(); ++i) {
            ret[i] = _models[i]->GetName();
        }
        return ret;
    }
 
    std::unique_ptr<Model> GetModel(int idx) const {
        if (_models.empty()) {
            LogError(_runtime_config.logging.logger_name,
                     "No models available in GetModel.");
            throw std::out_of_range("Error attempting to GetModel: no models.");
        }
        if (idx < -1 || idx >= static_cast<int>(_models.size())) {
            LogError(_runtime_config.logging.logger_name,
                     "Index out of range in GetModel: " + std::to_string(idx));
            throw std::out_of_range("Error attempting to GetModel by index.");
        }
        // Return the last model if idx is -1
        if (idx == -1) {
            return _models.back()->clone();
        }
        return _models[idx]->clone();
    }
 
    std::vector<std::string> GetModelNames() const {
        std::vector<std::string> ret;
        for (auto &m : _models) {
            ret.push_back(m->GetName());
        }
        return ret;
    }
 
    const std::map<std::string, History> &GetModelHistory(size_t idx) const {
        if (idx >= _models.size()) {
            LogError(_runtime_config.logging.logger_name,
                     "Index out of range in GetModelHistory: " +
                         std::to_string(idx));
            throw std::out_of_range(
                "Error attempting to GetModelHistory by index.");
        }
        return _models[idx]->GetHistories();
    }
 
    const std::vector<std::string> GetModelHistoryNames(size_t idx) const {
        if (idx >= _models.size()) {
            LogError(_runtime_config.logging.logger_name,
                     "Index out of range in GetModelHistoryNames: " +
                         std::to_string(idx));
            throw std::out_of_range(
                "Error attempting to GetModelHistoryNames by index.");
        }
        std::vector<std::string> ret;
        for (const auto &kv : _models[idx]->GetHistories()) {
            ret.push_back(kv.first);
        }
        return ret;
    }
 
    void SetDuration(int duration) {
        if (duration <= 0) {
            throw std::invalid_argument(
                "Simulation duration must be positive.");
        }
        _duration = duration;
    }
 
    const ExecutionConfig &GetExecutionConfig() const {
        return _runtime_config.execution;
    }
 
    void SetExecutionConfig(const ExecutionConfig &execution_config) {
        _runtime_config.execution = execution_config;
    }
 
    const RuntimeConfig &GetRuntimeConfig() const { return _runtime_config; }
 
    void SetRuntimeConfig(const RuntimeConfig &runtime_config) {
        if (ConfigureLogger(runtime_config.logging) == CreationStatus::kError) {
            LogError(_runtime_config.logging.logger_name,
                     "Unable to apply simulation runtime logging config.");
            throw std::invalid_argument(
                "Error attempting to apply simulation runtime logging "
                "configuration.");
        }
        _runtime_config = runtime_config;
    }
 
private:
    Model &GetModelRefOrThrow(size_t idx) {
        if (idx >= _models.size()) {
            LogError(_runtime_config.logging.logger_name,
                     "Index out of range in model access: " +
                         std::to_string(idx));
            throw std::out_of_range("Error attempting to access model by "
                                    "index.");
        }
        return *_models[idx];
    }
 
    const Model &GetModelRefOrThrow(size_t idx) const {
        if (idx >= _models.size()) {
            LogError(_runtime_config.logging.logger_name,
                     "Index out of range in model access: " +
                         std::to_string(idx));
            throw std::out_of_range("Error attempting to access model by "
                                    "index.");
        }
        return *_models[idx];
    }
 
    RuntimeConfig _runtime_config;
    std::vector<std::unique_ptr<Model>> _models;
 
    int _duration = 1; // Default simulation duration in timesteps
    std::vector<int> _parameter_change_times;
    bool _stratify_entering_cohort = false;
 
    bool _build_summary_stats = false;
    bool _save_state_history = false;
    std::vector<int> _timesteps_to_report;
    bool _pivot_long = false;
};
} // namespace respond
 
#endif // RESPOND_SIMULATION_HPP_