Line data    Source code

       1             : /* 
       2             : * Copyright (C) 2020-2024 MEmilio
       3             : *
       4             : * Authors: Lena Ploetzke
       5             : *
       6             : * Contact: Martin J. Kuehn <Martin.Kuehn@DLR.de>
       7             : *
       8             : * Licensed under the Apache License, Version 2.0 (the "License");
       9             : * you may not use this file except in compliance with the License.
      10             : * You may obtain a copy of the License at
      11             : *
      12             : *     http://www.apache.org/licenses/LICENSE-2.0
      13             : *
      14             : * Unless required by applicable law or agreed to in writing, software
      15             : * distributed under the License is distributed on an "AS IS" BASIS,
      16             : * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
      17             : * See the License for the specific language governing permissions and
      18             : * limitations under the License.
      19             : */
      20             : 
      21             : #ifndef IDE_SEIR_H
      22             : #define IDE_SEIR_H
      23             : 
      24             : #include "memilio/math/eigen.h"
      25             : #include "memilio/utils/time_series.h"
      26             : #include "ide_seir/parameters.h"
      27             : #include "ide_seir/infection_state.h"
      28             : #include "memilio/utils/logging.h"
      29             : 
      30             : #include <iostream>
      31             : 
      32             : #include <vector>
      33             : 
      34             : namespace mio
      35             : {
      36             : namespace iseir
      37             : {
      38             : template <typename FP = double>
      39             : class Model
      40             : {
      41             :     using Pa  = ParametersBase<FP>;
      42             :     using Vec = TimeSeries<double>::Vector;
      43             : 
      44             : public:
      45             :     /**
      46             :         * @brief Create an IDE SEIR model.
      47             :         *
      48             :         * @param[in, out] init TimeSeries with the initial values of the number of susceptibles at associated initial times.
      49             :         *   The time steps in this vector should be equidistant and equal to the time step used for the simulation. 
      50             :         *   A certain history of time steps and values for susceptibles is needed. 
      51             :         *   A warning is displayed if the condition is violated.
      52             :         *   Co be more precise, the first time point needs to be smaller than -(k-1)*TimeStep with 
      53             :         *       k=ceil((InfectiousTime + LatencyTime)/TimeStep).
      54             :         *   The last time point in this vector should be a time 0.
      55             :         * @param[in] dt_init The size of the time step used for numerical simulation.
      56             :         * @param[in] N_init The population of the considered region. 
      57             :         */
      58           2 :     Model(TimeSeries<double>&& init, double dt_init, int N_init, const Pa& Parameterset_init = Pa())
      59           2 :         : parameters{Parameterset_init}
      60           2 :         , m_result{std::move(init)}
      61           2 :         , m_result_SEIR{TimeSeries<double>(4)}
      62           2 :         , m_dt{dt_init}
      63           2 :         , m_N{N_init}
      64             :     {
      65           2 :     }
      66             : 
      67             :     /**
      68             :         * @brief Simulate the evolution of infection numbers with the given IDE SEIR model.
      69             :         *
      70             :         * The simulation is performed by solving the underlying model equation numerically. 
      71             :         * Here, an integro-differential equation is to be solved. The model parameters and the initial data are used.
      72             :         *
      73             :         * @param[in] t_max Last simulation day. 
      74             :         *   If the last point of time of the initial TimeSeries was 0, the simulation will be executed for t_max days.
      75             :         * @return The result of the simulation, stored in a TimeSeries with simulation time and 
      76             :         *       associated number of susceptibles.
      77             :         */
      78           2 :     TimeSeries<double> const& simulate(int t_max)
      79             :     {
      80             : 
      81           2 :         m_l = (int)std::floor(parameters.template get<LatencyTime>() / m_dt);
      82           2 :         m_k =
      83           2 :             (int)std::ceil((parameters.template get<InfectiousTime>() + parameters.template get<LatencyTime>()) / m_dt);
      84           2 :         if (m_result.get_time(0) > -(m_k - 1) * m_dt) {
      85           0 :             log_warning("Constraint check: Initial data starts later than necessary. The simulation may be distorted. "
      86             :                         "Start the data at time {:.4f} at the latest.",
      87           0 :                         -(m_k - 1) * m_dt);
      88             :         }
      89             : 
      90         162 :         while (m_result.get_last_time() < t_max) {
      91         160 :             m_result.add_time_point(m_result.get_last_time() + m_dt);
      92         160 :             Eigen::Index idx = m_result.get_num_time_points();
      93             : 
      94             :             // R0t is the effective reproduction number at time t
      95         320 :             auto R0t1 = parameters.template get<ContactFrequency<double>>().get_cont_freq_mat().get_matrix_at(
      96         320 :                             m_result.get_time(idx - 2))(0, 0) *
      97         160 :                         parameters.template get<TransmissionRisk>() * parameters.template get<InfectiousTime>();
      98         320 :             auto R0t2 = parameters.template get<ContactFrequency<double>>().get_cont_freq_mat().get_matrix_at(
      99         320 :                             m_result.get_last_time())(0, 0) *
     100         160 :                         parameters.template get<TransmissionRisk>() * parameters.template get<InfectiousTime>();
     101             : 
     102         320 :             m_result.get_last_value() =
     103         320 :                 Vec::Constant(1, m_result[idx - 2][0] * exp(m_dt * (0.5 * 1 / m_N) *
     104         160 :                                                             (R0t1 * num_integration_inner_integral(idx - 2) +
     105         160 :                                                              R0t2 * num_integration_inner_integral(idx - 1))));
     106             :         }
     107           2 :         return m_result;
     108             :     }
     109             : 
     110             :     /**
     111             :         * @brief Calculate the distribution of the population in E, I and, R based on the calculated values for S.
     112             :         *
     113             :         * The values are calculated using the average latency and infection time, not using model equations. 
     114             :         * The simulated values of S are used for this purpose, so the simulate() function should be called beforehand.
     115             :         * 
     116             :         * @return The result of the calculation stored in an TimeSeries. The TimeSeries contains the simulation time and an
     117             :         *   associated Vector with values for S, E, I, and R.
     118             :         */
     119           2 :     TimeSeries<double> const& calculate_EIR()
     120             :     {
     121           2 :         Eigen::Index num_points = m_result.get_num_time_points();
     122           2 :         double S, E, I, R;
     123         206 :         for (Eigen::Index i = m_k; i < num_points; ++i) {
     124         204 :             S = m_result[i][Eigen::Index(InfectionState::S)];
     125         204 :             E = m_result[i - m_l][Eigen::Index(InfectionState::S)] - S;
     126         204 :             I = m_result[i - m_k][Eigen::Index(InfectionState::S)] - m_result[i - m_l][Eigen::Index(InfectionState::S)];
     127         204 :             R = m_N - S - E - I;
     128         204 :             m_result_SEIR.add_time_point(m_result.get_time(i), (Vec(4) << S, E, I, R).finished());
     129             :         }
     130           4 :         return m_result_SEIR;
     131             :     }
     132             : 
     133             :     // Used Parameters for the simulation.
     134             :     Pa parameters{};
     135             : 
     136             : private:
     137             :     /**
     138             :         * @brief Density of the generalized beta distribution used for the function f_{beta} of the IDE SEIR model.
     139             :         *
     140             :         * @param[in] tau evaluation point of the generalized Beta distribution.
     141             :         * @param[in] p parameter p of the generalized Beta distribution.
     142             :         * @param[in] q parameter q of the generalized Beta distribution.
     143             :         * @result Evaluation of the generalized beta distribution at the given evaluation point.
     144             :         */
     145       26560 :     double generalized_beta_distribution(double tau, double p = 3.0, double q = 10.0) const
     146             :     {
     147       52800 :         if ((parameters.template get<LatencyTime>() < tau) &&
     148       26240 :             (parameters.template get<InfectiousTime>() + parameters.template get<LatencyTime>() > tau)) {
     149       26240 :             return tgamma(p + q) * pow(tau - parameters.template get<LatencyTime>(), p - 1) *
     150       26240 :                    pow(parameters.template get<InfectiousTime>() + parameters.template get<LatencyTime>() - tau,
     151             :                        q - 1) /
     152       26240 :                    (tgamma(p) * tgamma(q) * pow(parameters.template get<InfectiousTime>(), p + q - 1));
     153             :         }
     154         320 :         return 0.0;
     155             :     }
     156             : 
     157             :     /**
     158             :         * @brief Numerical differentiation of one compartment using a central difference quotient.
     159             :         * 
     160             :         * @param[in] ts_ide TimeSeries with the time steps already calculated. 
     161             :         *       Used as function values in numerical differentiation.
     162             :         * @param[in] idx Time index at which the numerical differentiation should be performed.
     163             :         * @param[in] compartment Compartment for which the numerical differentiation is to be performed.
     164             :         * @return Numerically approximated derivative of the function belonging to the compartment at the point t[idx].
     165             :         */
     166       26560 :     double central_difference_quotient(TimeSeries<double> const& ts_ide, InfectionState compartment,
     167             :                                        Eigen::Index idx) const
     168             :     {
     169       26560 :         return (ts_ide[idx + 1][Eigen::Index(compartment)] - ts_ide[idx - 1][Eigen::Index(compartment)]) / (2 * m_dt);
     170             :     }
     171             : 
     172             :     /**
     173             :         * @brief Numerical integration of the inner integral of the integro-differential equation for the group S using
     174             :         *    a trapezoidal sum.
     175             :         * 
     176             :         * @param[in] idx Index of the point of time used in the inner integral.
     177             :         * @return Result of the numerical integration.
     178             :         */
     179         320 :     double num_integration_inner_integral(Eigen::Index idx) const
     180             :     {
     181         320 :         double res     = 0.5 * (generalized_beta_distribution(m_result.get_time(idx) - m_result.get_time(idx - m_k)) *
     182         320 :                                 central_difference_quotient(m_result, InfectionState::S, m_k) +
     183         320 :                             generalized_beta_distribution(m_result.get_time(idx) - m_result.get_time(idx - m_l)) *
     184         320 :                                 central_difference_quotient(m_result, InfectionState::S, idx - m_l));
     185         320 :         Eigen::Index i = idx - m_k + 1;
     186       26240 :         while (i <= idx - m_l - 2) {
     187       25920 :             res += (generalized_beta_distribution(m_result.get_time(idx) - m_result.get_time(i)) *
     188       25920 :                     central_difference_quotient(m_result, InfectionState::S, i));
     189       25920 :             ++i;
     190             :         }
     191         320 :         return res;
     192             :     }
     193             : 
     194             :     // TimeSeries containing points of time and the corresponding number of susceptibles.
     195             :     TimeSeries<double> m_result;
     196             :     // TimeSeries containing points of time and the corresponding number of susceptibles, exposed,
     197             :     // infected and recovered.
     198             :     TimeSeries<double> m_result_SEIR;
     199             : 
     200             :     // Timestep used for simulation.
     201             :     double m_dt{0};
     202             :     // Population of the considered region.
     203             :     int m_N{0};
     204             : 
     205             :     // Two Indices used for simulation.
     206             :     Eigen::Index m_k{0};
     207             :     Eigen::Index m_l{0};
     208             : };
     209             : } // namespace iseir
     210             : } // namespace mio
     211             : #endif