Skip to content

Commit

Permalink
changes undone, only time independent
Browse files Browse the repository at this point in the history
  • Loading branch information
@johapau
johapau committed Nov 20, 2023
1 parent e33b287 commit 200c9b7
Show file tree
Hide file tree
Showing 5,351 changed files with 916,336 additions and 248 deletions.
The diff you're trying to view is too large. We only load the first 3000 changed files.
248 changes: 0 additions & 248 deletions cpp/models/ode_secir/model.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,15 +30,6 @@
#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 @@ -532,245 +523,6 @@ 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.
Expand Down
Binary file modified cpp/thirdparty/boost_1_75_0.tar.gz
View file
Binary file not shown.
Binary file added cpp/thirdparty/boost_1_75_0.tar.gz.old
View file
Binary file not shown.
Loading

0 comments on commit 200c9b7

Please sign in to comment.