From d4cb5ee1c15881901d0c5e1422fcfee753d277a0 Mon Sep 17 00:00:00 2001 From: lenaploetzke <70579874+lenaploetzke@users.noreply.github.com> Date: Fri, 16 Feb 2024 12:11:49 +0100 Subject: [PATCH] 769 new initialization for ide model (#770) Added a different method to initialize the ide model. Is comparable to the method used for the ODE model for compartment Recovered. Co-authored-by: annawendler <106674756+annawendler@users.noreply.github.com> --- cpp/examples/ide_secir.cpp | 10 +-- cpp/models/ide_secir/model.cpp | 124 ++++++++++++++++----------- cpp/models/ide_secir/model.h | 30 +++++-- cpp/tests/test_ide_secir.cpp | 148 ++++++++++++++++++++++++++++----- 4 files changed, 230 insertions(+), 82 deletions(-) diff --git a/cpp/examples/ide_secir.cpp b/cpp/examples/ide_secir.cpp index 685bed5a56..517fac5da5 100644 --- a/cpp/examples/ide_secir.cpp +++ b/cpp/examples/ide_secir.cpp @@ -32,10 +32,10 @@ int main() { using Vec = mio::TimeSeries::Vector; - ScalarType tmax = 10; - ScalarType N = 10000; - ScalarType Dead_before = 12; - ScalarType dt = 1; + ScalarType tmax = 10; + ScalarType N = 10000; + ScalarType deaths = 13.10462213; + ScalarType dt = 1; int num_transitions = (int)mio::isecir::InfectionTransition::Count; @@ -63,7 +63,7 @@ int main() } // Initialize model. - mio::isecir::Model model(std::move(init), N, Dead_before); + mio::isecir::Model model(std::move(init), N, deaths); // model.m_populations.get_last_value()[(Eigen::Index)mio::isecir::InfectionState::Susceptible] = 1000; // model.m_populations.get_last_value()[(Eigen::Index)mio::isecir::InfectionState::Recovered] = 0; diff --git a/cpp/models/ide_secir/model.cpp b/cpp/models/ide_secir/model.cpp index 196e6dd775..d9e33f3ab9 100644 --- a/cpp/models/ide_secir/model.cpp +++ b/cpp/models/ide_secir/model.cpp @@ -29,64 +29,33 @@ namespace mio namespace isecir { -Model::Model(TimeSeries&& init, ScalarType N_init, ScalarType Dead_before, +Model::Model(TimeSeries&& init, ScalarType N_init, ScalarType deaths, ScalarType total_confirmed_cases, const ParameterSet& Parameterset_init) : parameters{Parameterset_init} , m_transitions{std::move(init)} , m_populations{TimeSeries(Eigen::Index(InfectionState::Count))} , m_N{N_init} - , m_deaths_before{Dead_before} + , m_total_confirmed_cases{total_confirmed_cases} { + m_deaths_before = + deaths - m_transitions.get_last_value()[Eigen::Index(InfectionTransition::InfectedCriticalToDead)]; m_populations.add_time_point( 0, TimeSeries::Vector::Constant((int)InfectionState::Count, 0)); - m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Dead)] = - m_deaths_before + m_transitions.get_last_value()[Eigen::Index(InfectionTransition::InfectedCriticalToDead)]; + m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Dead)] = deaths; } void Model::initialize(ScalarType dt) { - // compute Susceptibles at time 0 and m_forceofinfection at time -m_dt as initial values for discretization scheme - // use m_forceofinfection at -m_dt to be consistent with further calculations of S (see compute_susceptibles()), - // where also the value of m_forceofinfection for the previous timestep is used - update_forceofinfection(dt, true); - if (m_forceofinfection > 0) { - m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Susceptible)] = - m_transitions.get_last_value()[Eigen::Index(InfectionTransition::SusceptibleToExposed)] / - (dt * m_forceofinfection); - - //calculate other compartment sizes for t=0 + if (m_total_confirmed_cases > 1e-12) { + m_initialization_method = 1; other_compartments_current_timestep(dt); - //R; need an initial value for R, therefore do not calculate via compute_recovered() + // The scheme of the ODE model for initialization is applied here. m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Recovered)] = - m_N - m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Susceptible)] - - m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Exposed)] - - m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedNoSymptoms)] - - m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedSymptoms)] - + m_total_confirmed_cases - m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedSymptoms)] - m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedSevere)] - m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedCritical)] - m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Dead)]; - } - else if (m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Susceptible)] > 1e-12) { - //take initialized value for Susceptibles if value can't be calculated via the standard formula - //calculate other compartment sizes for t=0 - other_compartments_current_timestep(dt); - - //R; need an initial value for R, therefore do not calculate via compute_recovered() - m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Recovered)] = - m_N - m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Susceptible)] - - m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Exposed)] - - m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedNoSymptoms)] - - m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedSymptoms)] - - m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedSevere)] - - m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedCritical)] - - m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Dead)]; - } - else if (m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Recovered)] > 1e-12) { - //if value for Recovered is initialized and standard method is not applicable, calculate Susceptibles via other compartments - //determining other compartment sizes is not dependent of Susceptibles(0), just of the transitions of the past. - //calculate other compartment sizes for t=0 - other_compartments_current_timestep(dt); m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Susceptible)] = m_N - m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Exposed)] - @@ -98,10 +67,69 @@ void Model::initialize(ScalarType dt) m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Dead)]; } else { - log_error("Error occured while initializing compartments: Force of infection is evaluated to 0 and neither " - "Susceptibles nor Recovered for time 0 were set. One of them should be larger 0."); - } + // compute Susceptibles at time 0 and m_forceofinfection at time -dt as initial values for discretization scheme + // use m_forceofinfection at -dt to be consistent with further calculations of S (see compute_susceptibles()), + // where also the value of m_forceofinfection for the previous timestep is used + update_forceofinfection(dt, true); + if (m_forceofinfection > 1e-12) { + m_initialization_method = 2; + m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Susceptible)] = + m_transitions.get_last_value()[Eigen::Index(InfectionTransition::SusceptibleToExposed)] / + (dt * m_forceofinfection); + + //calculate other compartment sizes for t=0 + other_compartments_current_timestep(dt); + + //R; need an initial value for R, therefore do not calculate via compute_recovered() + m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Recovered)] = + m_N - m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Susceptible)] - + m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Exposed)] - + m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedNoSymptoms)] - + m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedSymptoms)] - + m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedSevere)] - + m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedCritical)] - + m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Dead)]; + } + else if (m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Susceptible)] > 1e-12) { + //take initialized value for Susceptibles if value can't be calculated via the standard formula + m_initialization_method = 3; + //calculate other compartment sizes for t=0 + other_compartments_current_timestep(dt); + + //R; need an initial value for R, therefore do not calculate via compute_recovered() + m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Recovered)] = + m_N - m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Susceptible)] - + m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Exposed)] - + m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedNoSymptoms)] - + m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedSymptoms)] - + m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedSevere)] - + m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedCritical)] - + m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Dead)]; + } + else if (m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Recovered)] > 1e-12) { + //if value for Recovered is initialized and standard method is not applicable, calculate Susceptibles via other compartments + //determining other compartment sizes is not dependent of Susceptibles(0), just of the transitions of the past. + //calculate other compartment sizes for t=0 + m_initialization_method = 4; + other_compartments_current_timestep(dt); + + m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Susceptible)] = + m_N - m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Exposed)] - + m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedNoSymptoms)] - + m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedSymptoms)] - + m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedSevere)] - + m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::InfectedCritical)] - + m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Recovered)] - + m_populations[Eigen::Index(0)][Eigen::Index(InfectionState::Dead)]; + } + else { + m_initialization_method = -1; + log_error("Error occured while initializing compartments: Force of infection is evaluated to 0 and neither " + "Susceptibles nor Recovered or total confirmed cases for time 0 were set. One of them should be " + "larger 0."); + } + } // compute m_forceofinfection at time 0 needed for further simulation update_forceofinfection(dt); } @@ -118,11 +146,11 @@ void Model::compute_susceptibles(ScalarType dt) void Model::compute_flow(int idx_InfectionTransitions, Eigen::Index idx_IncomingFlow, ScalarType dt) { ScalarType sum = 0; - /* In order to satisfy TransitionDistribution(m_dt*i) = 0 for all i >= k, k is determined by the maximum support of the distribution. + /* In order to satisfy TransitionDistribution(dt*i) = 0 for all i >= k, k is determined by the maximum support of the distribution. Since we are using a backwards difference scheme to compute the derivative, we have that the - derivative of TransitionDistribution(m_dt*i) = 0 for all i >= k+1. + derivative of TransitionDistribution(dt*i) = 0 for all i >= k+1. - Hence calc_time_index goes until std::ceil(support_max/m_dt) since for std::ceil(support_max/m_dt)+1 all terms are already zero. + Hence calc_time_index goes until std::ceil(support_max/dt) since for std::ceil(support_max/dt)+1 all terms are already zero. This needs to be adjusted if we are changing the finite difference scheme */ Eigen::Index calc_time_index = (Eigen::Index)std::ceil( @@ -131,7 +159,7 @@ void Model::compute_flow(int idx_InfectionTransitions, Eigen::Index idx_Incoming Eigen::Index num_time_points = m_transitions.get_num_time_points(); for (Eigen::Index i = num_time_points - 1 - calc_time_index; i < num_time_points - 1; i++) { - // (num_time_points - 1 - i)* m_dt is the time, the individuals has already spent in this state. + // (num_time_points - 1 - i)* dt is the time, the individuals has already spent in this state. ScalarType state_age = (num_time_points - 1 - i) * dt; @@ -214,7 +242,7 @@ void Model::update_forceofinfection(ScalarType dt, bool initialization) ScalarType deaths; if (initialization) { - // determine m_forceofinfection at time -m_dt which is the penultimate timepoint in m_transitions + // determine m_forceofinfection at time -dt which is the penultimate timepoint in m_transitions num_time_points = m_transitions.get_num_time_points() - 1; current_time = -dt; deaths = m_deaths_before; diff --git a/cpp/models/ide_secir/model.h b/cpp/models/ide_secir/model.h index cfcf3c9822..951ddcea74 100644 --- a/cpp/models/ide_secir/model.h +++ b/cpp/models/ide_secir/model.h @@ -39,17 +39,17 @@ class Model * @brief Constructor to create an IDE SECIR model. * * @param[in, out] init TimeSeries with the initial values of the number of individuals, - * which transit within one timestep dt_init from one compartment to another. + * which transit within one timestep dt from one compartment to another. * Possible transitions are specified in as #InfectionTransition%s. - * Considered points of times should have the distance dt_init and the last time point should be 0. + * Considered points of times should have the distance dt and the last time point should be 0. * The time history must reach a certain point in the past so that the simulation can be performed. * A warning is displayed if the condition is violated. - * @param[in] dt_init The size of the time step used for numerical simulation. * @param[in] N_init The population of the considered region. - * @param[in] Dead_before The total number of deaths at the time point - dt_init. + * @param[in] deaths The total number of deaths at the time zero. + * @param[in] total_confirmed_cases Total confirmed cases at time t0 can be set if it should be used for initialisation. * @param[in, out] Parameterset_init Used Parameters for simulation. */ - Model(TimeSeries&& init, ScalarType N_init, ScalarType Dead_before, + Model(TimeSeries&& init, ScalarType N_init, ScalarType deaths, ScalarType total_confirmed_cases = 0, const ParameterSet& Parameterset_init = ParameterSet()); /** @@ -207,6 +207,21 @@ class Model m_tol = new_tol; } + /** + * @brief Specifies a number associated with the method used for initialization. + * + * @returns 0 if the initialization method has not yet been selected, + * 1 if the method using the total number of confirmed cases at time 0 is used, + * 2 if the force of infection method is used, + * 3 if the initialization is calculated using a prior set value for S, + * 4 if the initialization is calculated using a prior set value for R and + * -1 if the initialization was not possible using any of the methods. + */ + int get_initialization_method() + { + return m_initialization_method; + } + ParameterSet parameters{}; ///< ParameterSet of Model Parameters. /* Attention: m_populations and m_transitions do not necessarily have the same number of time points due to the initialization part. */ TimeSeries @@ -217,8 +232,11 @@ class Model private: ScalarType m_forceofinfection{0}; ///< Force of infection term needed for numerical scheme. ScalarType m_N{0}; ///< Total population size of the considered region. - ScalarType m_deaths_before{0}; ///< Deaths before start of simulation (at time -m_dt). + ScalarType m_deaths_before{0}; ///< Total number of deaths at the time point - dt. + ScalarType m_total_confirmed_cases{0}; ///< Total number of confirmed cases at time t0. ScalarType m_tol{1e-10}; ///< Tolerance used to calculate the maximum support of the TransitionDistributions. + int m_initialization_method{ + 0}; ///< Gives the index of the method used for the initialization of the model. See also get_initialization_method() for the number code. }; } // namespace isecir diff --git a/cpp/tests/test_ide_secir.cpp b/cpp/tests/test_ide_secir.cpp index 23bd2d935d..e5c523bd0b 100755 --- a/cpp/tests/test_ide_secir.cpp +++ b/cpp/tests/test_ide_secir.cpp @@ -41,8 +41,8 @@ class ModelTestIdeSecir : public testing::Test using Vec = mio::TimeSeries::Vector; //Set initial conditions - ScalarType N = 10000; - ScalarType Dead_before = 12; + ScalarType N = 10000; + ScalarType deaths = 13.10462213; int num_transitions = (int)mio::isecir::InfectionTransition::Count; @@ -65,7 +65,7 @@ class ModelTestIdeSecir : public testing::Test } // Initialize model - model = new mio::isecir::Model(std::move(init), N, Dead_before); + model = new mio::isecir::Model(std::move(init), N, deaths); // Set working parameters. mio::SmootherCosine smoothcos(2.0); @@ -163,10 +163,10 @@ TEST(IdeSecir, checkSimulationFunctions) { using Vec = mio::TimeSeries::Vector; - ScalarType tmax = 0.5; - ScalarType N = 10000; - ScalarType Dead_before = 10; - ScalarType dt = 0.5; + ScalarType tmax = 0.5; + ScalarType N = 10000; + ScalarType deaths = 10; + ScalarType dt = 0.5; int num_transitions = (int)mio::isecir::InfectionTransition::Count; @@ -192,7 +192,7 @@ TEST(IdeSecir, checkSimulationFunctions) } // Initialize model. - mio::isecir::Model model(std::move(init), N, Dead_before); + mio::isecir::Model model(std::move(init), N, deaths); // Set working parameters. // In this example, SmootherCosine with parameter 1 (and thus with a maximum support of 1) @@ -246,15 +246,117 @@ TEST(IdeSecir, checkSimulationFunctions) } } +// Check if the model uses the correct method for initialization using the function get_initialization_method(). +TEST(IdeSecir, checkInitializations) +{ + using Vec = mio::TimeSeries::Vector; + + ScalarType tmax = 1; + ScalarType N = 10000; + ScalarType deaths = 13.10462213; + ScalarType dt = 1; + + int num_transitions = (int)mio::isecir::InfectionTransition::Count; + + // Create TimeSeries with num_transitions elements where transitions needed for simulation will be stored. + mio::TimeSeries init(num_transitions); + // Add initial time point to time series. + init.add_time_point(-10, Vec::Constant(num_transitions, 3.0)); + // Add further time points until time 0. + while (init.get_last_time() < 0) { + init.add_time_point(init.get_last_time() + dt, Vec::Constant(num_transitions, 3.0)); + } + + // --- Case with total_confirmed_cases. + mio::TimeSeries init_copy1(init); + mio::isecir::Model model1(std::move(init_copy1), N, deaths, 1000); + + // Check that the initialization method is not already set. + EXPECT_EQ(0, model1.get_initialization_method()); + + // Carry out simulation. + mio::isecir::Simulation sim1(model1, 0, dt); + sim1.advance(tmax); + + // Verify that the expected initialization method was used. + EXPECT_EQ(1, sim1.get_model().get_initialization_method()); + + // --- Case with forceofinfection. + mio::TimeSeries init_copy2(init); + mio::isecir::Model model2(std::move(init_copy2), N, deaths, 0); + + // Carry out simulation. + mio::isecir::Simulation sim2(model2, 0, dt); + sim2.advance(tmax); + + // Verify that the expected initialization method was used. + EXPECT_EQ(2, sim2.get_model().get_initialization_method()); + + // --- Case with S. + /* !! For the other tests, the contact rate is set to 0 so that the force of infection is zero. + The forceofinfection initialization method is therefore not used for these tests.*/ + mio::isecir::Parameters parameters; + mio::ContactMatrixGroup contact_matrix = mio::ContactMatrixGroup(1, 1); + contact_matrix[0] = mio::ContactMatrix(Eigen::MatrixXd::Constant(1, 1, 0)); + parameters.get() = mio::UncertainContactMatrix(contact_matrix); + + mio::TimeSeries init_copy3(init); + mio::isecir::Model model3(std::move(init_copy3), N, deaths, 0, std::move(parameters)); + + model3.m_populations.get_last_value()[(Eigen::Index)mio::isecir::InfectionState::Susceptible] = 5000; + model3.m_populations.get_last_value()[(Eigen::Index)mio::isecir::InfectionState::Recovered] = 0; + + // Carry out simulation. + mio::isecir::Simulation sim3(model3, 0, dt); + sim3.advance(tmax); + + // Verify that the expected initialization method was used. + EXPECT_EQ(3, sim3.get_model().get_initialization_method()); + + // --- Case with R. + mio::TimeSeries init_copy4(init); + mio::isecir::Model model4(std::move(init_copy4), N, deaths, 0, std::move(parameters)); + + model4.m_populations.get_last_value()[(Eigen::Index)mio::isecir::InfectionState::Susceptible] = 0; + model4.m_populations.get_last_value()[(Eigen::Index)mio::isecir::InfectionState::Recovered] = 1000; + + // Carry out simulation. + mio::isecir::Simulation sim4(model4, 0, dt); + sim4.advance(tmax); + + // Verify that the expected initialization method was used. + EXPECT_EQ(4, sim4.get_model().get_initialization_method()); + + // --- Case without fitting initialization method. + // Deactivate temporarily log output for next test. Error is expected here. + mio::set_log_level(mio::LogLevel::off); + + // Here we do not need a copy of init as this is the last use of the vector. We can apply move directly. + mio::isecir::Model model5(std::move(init), N, deaths, 0, std::move(parameters)); + + model5.m_populations.get_last_value()[(Eigen::Index)mio::isecir::InfectionState::Susceptible] = 0; + model5.m_populations.get_last_value()[(Eigen::Index)mio::isecir::InfectionState::Recovered] = 0; + + // Carry out simulation. + mio::isecir::Simulation sim5(model5, 0, dt); + sim5.advance(tmax); + + // Verify that initialization was not possible with one of the models methods. + EXPECT_EQ(-1, sim5.get_model().get_initialization_method()); + + // Reactive log output. + mio::set_log_level(mio::LogLevel::warn); +} + // a) Test if check_constraints() function correctly reports wrongly set parameters. // b) Test if check_constraints() does not complain if parameters are set within correct ranges. TEST(IdeSecir, testValueConstraints) { using Vec = mio::TimeSeries::Vector; - ScalarType N = 10000; - ScalarType Dead_before = 10; - ScalarType dt = 1; + ScalarType N = 10000; + ScalarType deaths = 10; + ScalarType dt = 1; int num_transitions = (int)mio::isecir::InfectionTransition::Count; @@ -280,7 +382,7 @@ TEST(IdeSecir, testValueConstraints) } // Initialize a model. - mio::isecir::Model model(std::move(init), N, Dead_before); + mio::isecir::Model model(std::move(init), N, deaths); // Deactivate temporarily log output for next tests. mio::set_log_level(mio::LogLevel::off); @@ -408,10 +510,10 @@ TEST(IdeSecir, checkProportionRecoveredDeath) { using Vec = mio::TimeSeries::Vector; - ScalarType tmax = 30; - ScalarType N = 10000; - ScalarType Dead_before = 10; - ScalarType dt = 1; + ScalarType tmax = 30; + ScalarType N = 10000; + ScalarType deaths = 10; + ScalarType dt = 1; int num_transitions = (int)mio::isecir::InfectionTransition::Count; @@ -437,7 +539,7 @@ TEST(IdeSecir, checkProportionRecoveredDeath) } // Initialize model. - mio::isecir::Model model(std::move(init), N, Dead_before); + mio::isecir::Model model(std::move(init), N, deaths); // Set working parameters. // All TransitionDistribution%s are ExponentialDecay functions. @@ -496,10 +598,10 @@ TEST(IdeSecir, compareEquilibria) { using Vec = mio::TimeSeries::Vector; - ScalarType tmax = 20; - ScalarType N = 10000; - ScalarType Dead_before = 10; - ScalarType dt = 1; + ScalarType tmax = 20; + ScalarType N = 10000; + ScalarType deaths = 10; + ScalarType dt = 1; int num_transitions = (int)mio::isecir::InfectionTransition::Count; @@ -527,8 +629,8 @@ TEST(IdeSecir, compareEquilibria) mio::TimeSeries init2(init); // Initialize two models. - mio::isecir::Model model(std::move(init), N, Dead_before); - mio::isecir::Model model2(std::move(init2), N, Dead_before); + mio::isecir::Model model(std::move(init), N, deaths); + mio::isecir::Model model2(std::move(init2), N, deaths); // Set working parameters. // Here the maximum support for the TransitionDistribution%s is set differently for each model