Skip to content

Commit

Permalink
work in progress
Browse files Browse the repository at this point in the history
  • Loading branch information
@johapau
johapau committed Nov 6, 2023
1 parent b0b451c commit e33b287
Showing 1 changed file with 250 additions and 0 deletions.
250 changes: 250 additions & 0 deletions cpp/models/ode_secir/model.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,11 +24,21 @@
#include "memilio/compartments/simulation.h"
#include "memilio/config.h"
#include "memilio/epidemiology/populations.h"
#include "memilio/io/io.h"
#include "ode_secir/infection_state.h"
#include "memilio/math/interpolation.h"
#include "ode_secir/parameters.h"
#include "memilio/math/smoother.h"
#include "memilio/math/eigen_util.h"
#include "unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h"
#include <boost/math/tools/roots.hpp>
#include <boost/math/special_functions/next.hpp> // For float_distance.
#include <functional>
#include <limits>
#include <tuple> // for std::tuple and std::make_tuple.
#include <boost/math/special_functions/cbrt.hpp> // For boost::math::cbrt.
#include <boost/math/quadrature/gauss_kronrod.hpp>
#include <boost/math/differentiation/finite_difference.hpp>

namespace mio
{
Expand Down Expand Up @@ -317,6 +327,7 @@ double get_infections_relative(const Simulation<Base>& sim, double /*t*/, const
*@tparam Base simulation type that uses a secir compartment model. see Simulation.
*@returns The computed reproduction number at the provided index time.
*/

template <class Base>
IOResult<ScalarType> get_reproduction_number(size_t t_idx, const Simulation<Base>& sim)
{
Expand Down Expand Up @@ -522,11 +533,250 @@ IOResult<ScalarType> get_reproduction_number(size_t t_idx, const Simulation<Base
return mio::success(eigen_vals_abs.maxCoeff());
}

IOResult<ScalarType> calc_for_lambda(Eigen::MatrixXd F, Eigen::MatrixXd V, ScalarType lambda)
{
//return rho(W(t,0,lambda))
Eigen::MatrixXd M = 365.0 * (-V + (1 / lambda) * F);
//compute evolution operator
Eigen::MatrixXd evolution_op = M.exp();
//Compute the largest eigenvalue in absolute value
Eigen::ComplexEigenSolver<Eigen::MatrixXd> ces;

ces.compute(M);
const Eigen::VectorXcd eigen_vals = ces.eigenvalues();

Eigen::VectorXd eigen_vals_abs;
eigen_vals_abs.resize(eigen_vals.size());

for (int i = 0; i < eigen_vals.size(); i++) {
eigen_vals_abs[i] = std::abs(eigen_vals[i]);
}
return mio::success(eigen_vals_abs.maxCoeff());
}

template <class Base>
IOResult<ScalarType> get_reproduction_number_with_seasonality(size_t t_idx, const Simulation<Base>& sim)
{

if (!(t_idx < static_cast<size_t>(sim.get_result().get_num_time_points()))) {
return mio::failure(mio::StatusCode::OutOfRange, "t_idx is not a valid index for the TimeSeries");
}

auto const& params = sim.get_model().parameters;
const size_t num_groups = (size_t)sim.get_model().parameters.get_num_groups();
size_t num_infected_compartments = 5;
size_t total_infected_compartments = num_infected_compartments * num_groups;

Eigen::MatrixXd F(total_infected_compartments, total_infected_compartments);
Eigen::MatrixXd V(total_infected_compartments, total_infected_compartments);
F = Eigen::MatrixXd::Zero(total_infected_compartments,
total_infected_compartments); //Initialize matrices F and V with zeroes
V = Eigen::MatrixXd::Zero(total_infected_compartments, total_infected_compartments);

auto test_and_trace_required = 0.0;
auto icu_occupancy = 0.0;
for (auto i = AgeGroup(0); i < (mio::AgeGroup)num_groups; ++i) {
auto rateINS = 0.5 / (params.template get<IncubationTime>()[i] - params.template get<SerialInterval>()[i]);
test_and_trace_required +=
(1 - params.template get<RecoveredPerInfectedNoSymptoms>()[i]) * rateINS *
sim.get_result().get_value(
t_idx)[sim.get_model().populations.get_flat_index({i, InfectionState::InfectedNoSymptoms})];
icu_occupancy += sim.get_result().get_value(
t_idx)[sim.get_model().populations.get_flat_index({i, InfectionState::InfectedCritical})];
}

ContactMatrixGroup const& contact_matrix = sim.get_model().parameters.template get<ContactPatterns>();

Eigen::MatrixXd cont_freq_eff(num_groups, num_groups);
Eigen::VectorXd divN(num_groups);
Eigen::VectorXd riskFromInfectedSymptomatic(num_groups);
Eigen::MatrixXd riskFromInfectedSymptomatic_derivatives(num_groups, num_groups);
Eigen::VectorXd rateINS(num_groups);

for (mio::AgeGroup k = 0; k < (mio::AgeGroup)num_groups; k++) {
double temp = sim.get_result().get_value(
t_idx)[sim.get_model().populations.get_flat_index({k, InfectionState::Susceptible})] +
sim.get_result().get_value(
t_idx)[sim.get_model().populations.get_flat_index({k, InfectionState::Exposed})] +
sim.get_result().get_value(
t_idx)[sim.get_model().populations.get_flat_index({k, InfectionState::InfectedNoSymptoms})] +
sim.get_result().get_value(
t_idx)[sim.get_model().populations.get_flat_index({k, InfectionState::InfectedSymptoms})] +
sim.get_result().get_value(
t_idx)[sim.get_model().populations.get_flat_index({k, InfectionState::InfectedSevere})] +
sim.get_result().get_value(
t_idx)[sim.get_model().populations.get_flat_index({k, InfectionState::InfectedCritical})] +
sim.get_result().get_value(
t_idx)[sim.get_model().populations.get_flat_index({k, InfectionState::Recovered})];
if (temp == 0) {
temp = 1;
}
divN[(size_t)k] = 1 / temp;

riskFromInfectedSymptomatic[(size_t)k] = smoother_cosine(
test_and_trace_required, params.template get<TestAndTraceCapacity>(),
(params.template get<TestAndTraceCapacity>()) * 5, params.template get<RiskOfInfectionFromSymptomatic>()[k],
params.template get<MaxRiskOfInfectionFromSymptomatic>()[k]);

rateINS[(size_t)k] =
0.5 / (params.template get<IncubationTime>()[k] - params.template get<SerialInterval>()[(mio::AgeGroup)k]);

for (mio::AgeGroup l = 0; l < (mio::AgeGroup)num_groups; l++) {
if (test_and_trace_required < params.template get<TestAndTraceCapacity>() ||
test_and_trace_required > 5 * params.template get<TestAndTraceCapacity>()) {
riskFromInfectedSymptomatic_derivatives((size_t)k, (size_t)l) = 0;
}
else {
riskFromInfectedSymptomatic_derivatives((size_t)k, (size_t)l) =
-0.5 * 3.14159265358979323846 *
(params.template get<MaxRiskOfInfectionFromSymptomatic>()[k] -
params.template get<RiskOfInfectionFromSymptomatic>()[k]) /
(4 * params.template get<TestAndTraceCapacity>()) *
(1 - params.template get<RecoveredPerInfectedNoSymptoms>()[l]) * rateINS[(size_t)l] *
std::sin(3.14159265358979323846 / (4 * params.template get<TestAndTraceCapacity>()) *
(test_and_trace_required - params.template get<TestAndTraceCapacity>()));
}
}

for (Eigen::Index l = 0; l < (Eigen::Index)num_groups; l++) {
cont_freq_eff(l, (size_t)k) =
contact_matrix.get_matrix_at(static_cast<double>(t_idx))( //Discard season_val on purpose
static_cast<Eigen::Index>((size_t)l), static_cast<Eigen::Index>((size_t)k));
}
}

//Initialize the matrix F
for (size_t i = 0; i < num_groups; i++) {

for (size_t j = 0; j < num_groups; j++) {

double temp = 0;
for (Eigen::Index k = 0; k < (Eigen::Index)num_groups; k++) {
temp += cont_freq_eff(i, k) *
sim.get_result().get_value(t_idx)[sim.get_model().populations.get_flat_index(
{(mio::AgeGroup)k, InfectionState::InfectedSymptoms})] *
riskFromInfectedSymptomatic_derivatives(k, j) * divN[k];
}

F(i, j + num_groups) =
sim.get_result().get_value(t_idx)[sim.get_model().populations.get_flat_index(
{(mio::AgeGroup)i, InfectionState::Susceptible})] *
params.template get<TransmissionProbabilityOnContact>()[(mio::AgeGroup)i] *
(cont_freq_eff(i, j) * params.template get<RelativeTransmissionNoSymptoms>()[(mio::AgeGroup)j] *
divN[(size_t)j] +
temp);
}

for (size_t j = 0; j < num_groups; j++) {
F(i, j + 2 * num_groups) = sim.get_result().get_value(t_idx)[sim.get_model().populations.get_flat_index(
{(mio::AgeGroup)i, InfectionState::Susceptible})] *
params.template get<TransmissionProbabilityOnContact>()[(mio::AgeGroup)i] *
cont_freq_eff(i, j) * riskFromInfectedSymptomatic[(size_t)j] * divN[(size_t)j];
}
}

//Initialize the matrix V
for (Eigen::Index i = 0; i < (Eigen::Index)num_groups; i++) {

double rateE = 1.0 / (2 * params.template get<SerialInterval>()[(mio::AgeGroup)i] -
params.template get<IncubationTime>()[(mio::AgeGroup)i]);

double criticalPerSevereAdjusted = smoother_cosine(
icu_occupancy, 0.90 * params.template get<ICUCapacity>(), params.template get<ICUCapacity>(),
params.template get<CriticalPerSevere>()[(mio::AgeGroup)i], 0);

V(i, i) = rateE;
V(i + num_groups, i) = -rateE;
V(i + num_groups, i + num_groups) = rateINS[i];
V(i + 2 * num_groups, i + num_groups) =
-(1 - params.template get<RecoveredPerInfectedNoSymptoms>()[(mio::AgeGroup)i]) * rateINS[i];
V(i + 2 * num_groups, i + 2 * num_groups) = (1 / params.template get<TimeInfectedSymptoms>()[(mio::AgeGroup)i]);
V(i + 3 * num_groups, i + 2 * num_groups) =
-params.template get<SeverePerInfectedSymptoms>()[(mio::AgeGroup)i] /
params.template get<TimeInfectedSymptoms>()[(mio::AgeGroup)i];
V(i + 3 * num_groups, i + 3 * num_groups) = 1 / (params.template get<TimeInfectedSevere>()[(mio::AgeGroup)i]);
V(i + 4 * num_groups, i + 3 * num_groups) =
-criticalPerSevereAdjusted / (params.template get<TimeInfectedSevere>()[(mio::AgeGroup)i]);

V(i + 4 * num_groups, i + 4 * num_groups) = 1 / (params.template get<TimeInfectedCritical>()[(mio::AgeGroup)i]);

if (!(icu_occupancy < 0.9 * params.template get<ICUCapacity>() ||
icu_occupancy > (double)(params.template get<ICUCapacity>()))) {
for (size_t j = 0; j < num_groups; j++) {
V(i + 4 * num_groups, j + 4 * num_groups) -=
sim.get_result().get_value(t_idx)[sim.get_model().populations.get_flat_index(
{(mio::AgeGroup)i, InfectionState::InfectedSevere})] /
params.template get<TimeInfectedSevere>()[(mio::AgeGroup)i] * 5 *
params.template get<CriticalPerSevere>()[(mio::AgeGroup)i] * 3.141592653589793 /
(params.template get<ICUCapacity>()) *
std::sin(3.141592653589793 / (0.1 * params.template get<ICUCapacity>()) *
(icu_occupancy - 0.9 * params.template get<ICUCapacity>()));
}
}
}

//Set up the iteration
ScalarType start_guess = get_reproduction_number(t_idx, sim);
boost::uintmax_t maxit = 30;
std::function<std::vector<double>(ScalarType)> calc_val = [F, V](ScalarType lambda) {
//return spectral radius and derivative of spectral radius
//return rho(W(t,0,lambda))
Eigen::MatrixXd M = 365.0 * (-V + (1 / lambda) * F);
//compute evolution operator
Eigen::MatrixXd evolution_op = M.exp();
//Compute the largest eigenvalue in absolute value
Eigen::ComplexEigenSolver<Eigen::MatrixXd> ces;

ces.compute(M);
const Eigen::VectorXcd eigen_vals = ces.eigenvalues();
const Eigen::MatrixXcd eigen_vecs = ces.eigenvectors();

Eigen::VectorXcd eigen_vals_abs;
eigen_vals_abs.resize(eigen_vals.size());

for (int i = 0; i < eigen_vals.size(); i++) {
eigen_vals_abs[i] = std::abs(eigen_vals[i]);
}
auto it_max_ev =
std::find(eigen_vals_abs.data(), eigen_vals_abs.data() + eigen_vals_abs.size(), eigen_vals_abs.maxCoeff());

auto idx_max_ev = std::distance(eigen_vals_abs.data(), it_max_ev);
Eigen::VectorXcd max_eigen_vec = eigen_vecs.col(idx_max_ev);

//Now calculate derivative of exp(A(t)) evaluated at max_eigen_vec at lambda
Eigen::MatrixXcd MatDeriv(F.rows(), F.cols());
for (Eigen::Index i = 0; i < F.rows(); i++) {
for (Eigen::Index j = 0; j < F.cols(); j++) {
std::function<std::complex<double>(std::complex<double>)> integrand = [F, evolution_op, lambda, i,
j](ScalarType alpha) {
Eigen::MatrixXcd result(F.rows(), F.cols());
result = -(1 / (lambda * lambda)) * (alpha * evolution_op).exp() * F *
((1 - alpha) * evolution_op).exp();
return result(i, j);
};
MatDeriv(i, j) =
boost::math::quadrature::gauss_kronrod<std::complex<ScalarType>, 100>::integrate(integrand, 0, 1);
}
}

Eigen::VectorXcd term1(F.rows());
term1 = MatDeriv * max_eigen_vec;

//Calculate derivative of largest eigenvector numerically

Eigen::VectorXcd dx_max_eigen_vec(F.rows());
Eigen::VectorXd term2 = evolution_op * dx_max_eigen_vec;

std::function<ScalarType(ScalarType)> deriv_max_eigen_vec = return std::make_tuple(*(it_max_ev), );
};
}

/**
*@brief Computes the reproduction number for all time points of the Model output obtained by the Simulation.
*@param sim The Model Simulation.
*@returns vector containing all reproduction numbers
*/

template <class Base>
Eigen::VectorXd get_reproduction_numbers(const Simulation<Base>& sim)
{
Expand Down

0 comments on commit e33b287

Please sign in to comment.