lct_secir/model.h

Go to the documentation of this file.

/*
* Copyright (C) 2020-2026 MEmilio
*
* Authors: Lena Ploetzke
*
* 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 LCT_SECIR_MODEL_H
#define LCT_SECIR_MODEL_H
 
#include "lct_secir/parameters.h"
#include "lct_secir/infection_state.h"
#include "memilio/compartments/compartmental_model.h"
#include "memilio/epidemiology/lct_populations.h"
#include "memilio/epidemiology/lct_infection_state.h"
#include "memilio/config.h"
#include "memilio/utils/time_series.h"
#include "memilio/utils/logging.h"
#include "memilio/math/eigen.h"
#include "memilio/utils/type_list.h"
#include "memilio/utils/metaprogramming.h"
 
#include <numbers>
 
namespace mio
{
namespace lsecir
{
template <typename FP, class... LctStates>
class Model : public CompartmentalModel<FP, InfectionState, LctPopulations<FP, LctStates...>, Parameters<FP>>
{
public:
    using LctStatesGroups = TypeList<LctStates...>;
    using Base            = CompartmentalModel<FP, InfectionState, LctPopulations<FP, LctStates...>, Parameters<FP>>;
    using typename Base::ParameterSet;
    using typename Base::Populations;
    static size_t constexpr num_groups = sizeof...(LctStates);
    static_assert(num_groups >= 1, "The number of LctStates provided should be at least one.");
 
    Model()
        : Base(Populations(), ParameterSet(num_groups))
    {
    }
 
    Model(const Populations& pop, const ParameterSet& params)
        : Base(pop, params)
    {
    }
 
    void get_derivatives(Eigen::Ref<const Eigen::VectorX<FP>> pop, Eigen::Ref<const Eigen::VectorX<FP>> y, FP t,
                         Eigen::Ref<Eigen::VectorX<FP>> dydt) const override
    {
        // Vectors are sorted such that we first have all InfectionState%s for AgeGroup 0,
        // afterwards all for AgeGroup 1 and so on.
        dydt.setZero();
        get_derivatives_impl(pop, y, t, dydt);
    }
 
    TimeSeries<FP> calculate_compartments(const TimeSeries<FP>& subcompartments_ts) const
    {
        Eigen::Index count_InfStates  = (Eigen::Index)InfectionState::Count;
        Eigen::Index num_compartments = count_InfStates * num_groups;
        TimeSeries<FP> compartments_ts(num_compartments);
        if (!(this->populations.get_num_compartments() == (size_t)subcompartments_ts.get_num_elements())) {
            log_error("Result does not match InfectionState of the Model.");
            Eigen::VectorX<FP> wrong_size = Eigen::VectorX<FP>::Constant(num_compartments, -1);
            compartments_ts.add_time_point(-1, wrong_size);
            return compartments_ts;
        }
        Eigen::VectorX<FP> compartments(num_compartments);
        for (Eigen::Index timepoint = 0; timepoint < subcompartments_ts.get_num_time_points(); ++timepoint) {
            compress_vector(subcompartments_ts[timepoint], compartments);
            compartments_ts.add_time_point(subcompartments_ts.get_time(timepoint), compartments);
        }
 
        return compartments_ts;
    }
 
    bool check_constraints() const
    {
 
        return (check_constraints_impl() || this->parameters.check_constraints() ||
                this->populations.check_constraints());
    }
 
private:
    template <size_t Group = 0>
    void compress_vector(const Eigen::VectorX<FP>& subcompartments, Eigen::VectorX<FP>& compartments) const
    {
        static_assert((Group < num_groups) && (Group >= 0), "The template parameter Group should be valid.");
        using LctStateGroup = type_at_index_t<Group, LctStates...>;
 
        // Define first index of the group Group in a vector including all compartments without a resolution
        // in subcompartments.
        Eigen::Index count_InfStates         = (Eigen::Index)InfectionState::Count;
        Eigen::Index first_index_group_comps = Group * count_InfStates;
 
        // Use function from the LctState of the Group to calculate the vector without subcompartments
        // using the corresponding vector with subcompartments.
        compartments.segment(first_index_group_comps, count_InfStates) =
            LctStateGroup::calculate_compartments(subcompartments.segment(
                this->populations.template get_first_index_of_group<Group>(), LctStateGroup::Count));
 
        // Function call for next group if applicable.
        if constexpr (Group + 1 < num_groups) {
            compress_vector<Group + 1>(subcompartments, compartments);
        }
    }
 
    template <size_t Group = 0>
    void get_derivatives_impl(Eigen::Ref<const Eigen::VectorX<FP>> pop, Eigen::Ref<const Eigen::VectorX<FP>> y, FP t,
                              Eigen::Ref<Eigen::VectorX<FP>> dydt) const
    {
        static_assert((Group < num_groups) && (Group >= 0), "The template parameter Group should be valid.");
        using LctStateGroup = type_at_index_t<Group, LctStates...>;
 
        size_t first_index_group = this->populations.template get_first_index_of_group<Group>();
        const auto& params       = this->parameters;
        FP flow                  = 0;
 
        // Indices of first subcompartment of the InfectionState for the group in the vectors.
        size_t Ei_first_index = first_index_group + LctStateGroup::template get_first_index<InfectionState::Exposed>();
        size_t INSi_first_index =
            first_index_group + LctStateGroup::template get_first_index<InfectionState::InfectedNoSymptoms>();
        size_t ISyi_first_index =
            first_index_group + LctStateGroup::template get_first_index<InfectionState::InfectedSymptoms>();
        size_t ISevi_first_index =
            first_index_group + LctStateGroup::template get_first_index<InfectionState::InfectedSevere>();
        size_t ICri_first_index =
            first_index_group + LctStateGroup::template get_first_index<InfectionState::InfectedCritical>();
        size_t Ri = first_index_group + LctStateGroup::template get_first_index<InfectionState::Recovered>();
        size_t Di = first_index_group + LctStateGroup::template get_first_index<InfectionState::Dead>();
 
        // Calculate derivative of the Susceptible compartment.
        interact<Group, 0>(pop, y, t, dydt);
 
        // Calculate derivative of the Exposed compartment.
        dydt[Ei_first_index] = -dydt[first_index_group];
        for (size_t subcomp = 0; subcomp < LctStateGroup::template get_num_subcompartments<InfectionState::Exposed>();
             subcomp++) {
            // Variable flow stores the value of the flow from one subcompartment to the next one.
            // Ei_first_index + subcomp is always the index of a (sub-)compartment of Exposed and Ei_first_index
            // + subcomp + 1 can also be the index of the first (sub-)compartment of InfectedNoSymptoms.
            flow = (FP)LctStateGroup::template get_num_subcompartments<InfectionState::Exposed>() *
                   (1 / params.template get<TimeExposed<FP>>()[Group]) * y[Ei_first_index + subcomp];
            // Subtract flow from dydt[Ei_first_index + subcomp] and add to next subcompartment.
            dydt[Ei_first_index + subcomp] -= flow;
            dydt[Ei_first_index + subcomp + 1] = flow;
        }
 
        // Calculate derivative of the InfectedNoSymptoms compartment.
        for (size_t subcomp = 0;
             subcomp < LctStateGroup::template get_num_subcompartments<InfectionState::InfectedNoSymptoms>();
             subcomp++) {
            flow = (FP)LctStateGroup::template get_num_subcompartments<InfectionState::InfectedNoSymptoms>() *
                   (1 / params.template get<TimeInfectedNoSymptoms<FP>>()[Group]) * y[INSi_first_index + subcomp];
            dydt[INSi_first_index + subcomp] -= flow;
            dydt[INSi_first_index + subcomp + 1] = flow;
        }
 
        // Calculate derivative of the InfectedSymptoms compartment.
        // Flow from last (sub-) compartment of C must be split between
        // the first subcompartment of InfectedSymptoms and Recovered.
        dydt[Ri] = dydt[ISyi_first_index] * params.template get<RecoveredPerInfectedNoSymptoms<FP>>()[Group];
        dydt[ISyi_first_index] =
            dydt[ISyi_first_index] * (1 - params.template get<RecoveredPerInfectedNoSymptoms<FP>>()[Group]);
        for (size_t subcomp = 0;
             subcomp < LctStateGroup::template get_num_subcompartments<InfectionState::InfectedSymptoms>(); subcomp++) {
            flow = (FP)LctStateGroup::template get_num_subcompartments<InfectionState::InfectedSymptoms>() *
                   (1 / params.template get<TimeInfectedSymptoms<FP>>()[Group]) * y[ISyi_first_index + subcomp];
            dydt[ISyi_first_index + subcomp] -= flow;
            dydt[ISyi_first_index + subcomp + 1] = flow;
        }
 
        // Calculate derivative of the InfectedSevere compartment.
        dydt[Ri] += dydt[ISevi_first_index] * (1 - params.template get<SeverePerInfectedSymptoms<FP>>()[Group]);
        dydt[ISevi_first_index] = dydt[ISevi_first_index] * params.template get<SeverePerInfectedSymptoms<FP>>()[Group];
        for (size_t subcomp = 0;
             subcomp < LctStateGroup::template get_num_subcompartments<InfectionState::InfectedSevere>(); subcomp++) {
            flow = (FP)LctStateGroup::template get_num_subcompartments<InfectionState::InfectedSevere>() *
                   (1 / params.template get<TimeInfectedSevere<FP>>()[Group]) * y[ISevi_first_index + subcomp];
            dydt[ISevi_first_index + subcomp] -= flow;
            dydt[ISevi_first_index + subcomp + 1] = flow;
        }
 
        // Calculate derivative of the InfectedCritical compartment.
        dydt[Ri] += dydt[ICri_first_index] * (1 - params.template get<CriticalPerSevere<FP>>()[Group]);
        dydt[ICri_first_index] = dydt[ICri_first_index] * params.template get<CriticalPerSevere<FP>>()[Group];
        for (size_t subcomp = 0;
             subcomp < LctStateGroup::template get_num_subcompartments<InfectionState::InfectedCritical>() - 1;
             subcomp++) {
            flow = (FP)LctStateGroup::template get_num_subcompartments<InfectionState::InfectedCritical>() *
                   (1 / params.template get<TimeInfectedCritical<FP>>()[Group]) * y[ICri_first_index + subcomp];
            dydt[ICri_first_index + subcomp] -= flow;
            dydt[ICri_first_index + subcomp + 1] = flow;
        }
        // Last flow from InfectedCritical has to be divided between Recovered and Dead.
        // Must be calculated separately in order not to overwrite the already calculated values ​​for Recovered.
        flow = (FP)LctStateGroup::template get_num_subcompartments<InfectionState::InfectedCritical>() *
               (1 / params.template get<TimeInfectedCritical<FP>>()[Group]) * y[Ri - 1];
        dydt[Ri - 1] -= flow;
        dydt[Ri] = dydt[Ri] + (1 - params.template get<DeathsPerCritical<FP>>()[Group]) * flow;
        dydt[Di] = params.template get<DeathsPerCritical<FP>>()[Group] * flow;
 
        // Function call for next group if applicable.
        if constexpr (Group + 1 < num_groups) {
            get_derivatives_impl<Group + 1>(pop, y, t, dydt);
        }
    }
 
    template <size_t Group1, size_t Group2 = 0>
    void interact(Eigen::Ref<const Eigen::VectorX<FP>> pop, Eigen::Ref<const Eigen::VectorX<FP>> y, FP t,
                  Eigen::Ref<Eigen::VectorX<FP>> dydt) const
    {
        static_assert((Group1 < num_groups) && (Group1 >= 0) && (Group2 < num_groups) && (Group2 >= 0),
                      "The template parameters Group1 & Group2 should be valid.");
        using LctStateGroup2    = type_at_index_t<Group2, LctStates...>;
        size_t Si_1             = this->populations.template get_first_index_of_group<Group1>();
        FP infectedNoSymptoms_2 = 0;
        FP infectedSymptoms_2   = 0;
        auto params             = this->parameters;
 
        size_t first_index_group2 = this->populations.template get_first_index_of_group<Group2>();
 
        // Calculate sum of all subcompartments for InfectedNoSymptoms of Group2.
        infectedNoSymptoms_2 =
            pop.segment(first_index_group2 +
                            LctStateGroup2::template get_first_index<InfectionState::InfectedNoSymptoms>(),
                        LctStateGroup2::template get_num_subcompartments<InfectionState::InfectedNoSymptoms>())
                .sum();
        // Calculate sum of all subcompartments for InfectedSymptoms of Group2.
        infectedSymptoms_2 =
            pop.segment(first_index_group2 +
                            LctStateGroup2::template get_first_index<InfectionState::InfectedSymptoms>(),
                        LctStateGroup2::template get_num_subcompartments<InfectionState::InfectedSymptoms>())
                .sum();
        // Size of the Subpopulation Group2 without dead people.
        FP N_2          = pop.segment(first_index_group2, LctStateGroup2::Count - 1).sum();
        const FP divN_2 = (N_2 < Limits<FP>::zero_tolerance()) ? 0.0 : 1.0 / N_2;
        FP season_val   = 1 + params.template get<Seasonality<FP>>() *
                                std::sin(std::numbers::pi_v<ScalarType> *
                                         ((params.template get<StartDay<FP>>() + t) / 182.5 + 0.5));
        dydt[Si_1] += -y[Si_1] * divN_2 * season_val *
                      params.template get<TransmissionProbabilityOnContact<FP>>()[Group1] *
                      params.template get<ContactPatterns<FP>>().get_cont_freq_mat().get_matrix_at(
                          SimulationTime<FP>(t))(static_cast<Eigen::Index>(Group1), static_cast<Eigen::Index>(Group2)) *
                      (params.template get<RelativeTransmissionNoSymptoms<FP>>()[Group2] * infectedNoSymptoms_2 +
                       params.template get<RiskOfInfectionFromSymptomatic<FP>>()[Group2] * infectedSymptoms_2);
        // Function call for next interacting group if applicable.
        if constexpr (Group2 + 1 < num_groups) {
            interact<Group1, Group2 + 1>(pop, y, t, dydt);
        }
    }
 
    template <size_t Group = 0>
    bool check_constraints_impl() const
    {
        static_assert((Group < num_groups) && (Group >= 0), "The template parameter Group should be valid.");
        using LctStateGroup = type_at_index_t<Group, LctStates...>;
 
        if (LctStateGroup::template get_num_subcompartments<InfectionState::Susceptible>() != 1) {
            log_warning("Constraint check: The number of subcompartments for Susceptibles of group {} should be one!",
                        Group);
            return true;
        }
        if (LctStateGroup::template get_num_subcompartments<InfectionState::Recovered>() != 1) {
            log_warning("Constraint check: The number of subcompartments for Recovered of group {} should be one!",
                        Group);
            return true;
        }
        if (LctStateGroup::template get_num_subcompartments<InfectionState::Dead>() != 1) {
            log_warning("Constraint check: The number of subcompartments for Dead of group {} should be one!", Group);
            return true;
        }
 
        if constexpr (Group == num_groups - 1) {
            return false;
        }
        else {
            return check_constraints_impl<Group + 1>();
        }
    }
};
 
} // namespace lsecir
} // namespace mio
 
#endif // LCTSECIR_MODEL_H