Line data    Source code

       1             : /*
       2             : * Copyright (C) 2020-2024 MEmilio
       3             : *
       4             : * Authors: Lena Ploetzke, Anna Wendler
       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             : #include "ide_secir/parameters_io.h"
      22             : #include "memilio/config.h"
      23             : 
      24             : #ifdef MEMILIO_HAS_JSONCPP
      25             : 
      26             : #include "ide_secir/model.h"
      27             : #include "ide_secir/infection_state.h"
      28             : #include "memilio/math/eigen.h"
      29             : #include "memilio/io/epi_data.h"
      30             : #include "memilio/io/io.h"
      31             : #include "memilio/utils/date.h"
      32             : 
      33             : #include <string>
      34             : #include <cmath>
      35             : 
      36             : namespace mio
      37             : {
      38             : namespace isecir
      39             : {
      40             : 
      41           5 : IOResult<void> set_initial_flows(Model& model, ScalarType dt, std::string const& path, Date date,
      42             :                                  ScalarType scale_confirmed_cases)
      43             : {
      44             :     //--- Preparations ---
      45             :     // Try to get RKI data from path.
      46           5 :     BOOST_OUTCOME_TRY(auto&& rki_data, mio::read_confirmed_cases_noage(path));
      47           5 :     auto max_date_entry = std::max_element(rki_data.begin(), rki_data.end(), [](auto&& a, auto&& b) {
      48          56 :         return a.date < b.date;
      49             :     });
      50           5 :     if (max_date_entry == rki_data.end()) {
      51           2 :         log_error("RKI data file is empty.");
      52           1 :         return failure(StatusCode::InvalidFileFormat, path + ", file is empty.");
      53             :     }
      54           4 :     auto max_date = max_date_entry->date;
      55           4 :     if (max_date < date) {
      56           2 :         log_error("Specified date does not exist in RKI data.");
      57           1 :         return failure(StatusCode::OutOfRange, path + ", specified date does not exist in RKI data.");
      58             :     }
      59           3 :     auto min_date_entry = std::min_element(rki_data.begin(), rki_data.end(), [](auto&& a, auto&& b) {
      60          42 :         return a.date < b.date;
      61             :     });
      62           3 :     auto min_date       = min_date_entry->date;
      63             : 
      64             :     // Get (global) support_max to determine how many flows in the past we have to compute.
      65           3 :     ScalarType global_support_max         = model.get_global_support_max(dt);
      66           3 :     Eigen::Index global_support_max_index = Eigen::Index(std::ceil(global_support_max / dt));
      67             : 
      68             :     // m_transitions should be empty at the beginning.
      69           3 :     if (model.m_transitions.get_num_time_points() > 0) {
      70           2 :         model.m_transitions = TimeSeries<ScalarType>((int)InfectionTransition::Count);
      71             :     }
      72           3 :     if (model.m_populations.get_time(0) != 0) {
      73           1 :         model.m_populations.remove_last_time_point();
      74           2 :         model.m_populations.add_time_point<Eigen::VectorXd>(
      75           2 :             0, TimeSeries<ScalarType>::Vector::Constant((int)InfectionState::Count, 0));
      76             :     }
      77             : 
      78             :     // Set the Dead compartment to 0 so that RKI data can be added correctly.
      79           3 :     model.m_populations[0][Eigen::Index(InfectionState::Dead)] = 0;
      80             : 
      81             :     // Define variables for the mean of transitions used.
      82             :     ScalarType mean_ExposedToInfectedNoSymptoms =
      83           3 :         model.parameters.get<TransitionDistributions>()[Eigen::Index(InfectionTransition::ExposedToInfectedNoSymptoms)]
      84           3 :             .get_mean(dt);
      85             :     ScalarType mean_InfectedNoSymptomsToInfectedSymptoms =
      86             :         model.parameters
      87           3 :             .get<TransitionDistributions>()[(Eigen::Index)InfectionTransition::InfectedNoSymptomsToInfectedSymptoms]
      88           3 :             .get_mean(dt);
      89             :     ScalarType mean_InfectedSymptomsToInfectedSevere =
      90             :         model.parameters
      91           3 :             .get<TransitionDistributions>()[(Eigen::Index)InfectionTransition::InfectedSymptomsToInfectedSevere]
      92           3 :             .get_mean();
      93             :     ScalarType mean_InfectedSevereToInfectedCritical =
      94             :         model.parameters
      95           3 :             .get<TransitionDistributions>()[(Eigen::Index)InfectionTransition::InfectedSevereToInfectedCritical]
      96           3 :             .get_mean();
      97             :     ScalarType mean_InfectedCriticalToDead =
      98           3 :         model.parameters.get<TransitionDistributions>()[(Eigen::Index)InfectionTransition::InfectedCriticalToDead]
      99           3 :             .get_mean();
     100             : 
     101             :     // The first time we need is -4 * global_support_max.
     102           3 :     Eigen::Index start_shift = 4 * global_support_max_index;
     103             :     // The last time needed is dependent on the mean stay time in the Exposed compartment and
     104             :     // the mean stay time of asymptomatic individuals in InfectedNoSymptoms.
     105             :     Eigen::Index last_time_index_needed =
     106           3 :         Eigen::Index(std::ceil((mean_ExposedToInfectedNoSymptoms + mean_InfectedNoSymptomsToInfectedSymptoms) / dt));
     107             :     // Create TimeSeries with zeros. The index of time zero is start_shift.
     108          66 :     for (Eigen::Index i = -start_shift; i <= last_time_index_needed; i++) {
     109             :         // Add time point.
     110         126 :         model.m_transitions.add_time_point(
     111         126 :             i * dt, TimeSeries<ScalarType>::Vector::Constant((int)InfectionTransition::Count, 0.));
     112             :     }
     113             : 
     114             :     //--- Calculate the flow InfectedNoSymptomsToInfectedSymptoms using the RKI data and store in the m_transitions object.---
     115           3 :     ScalarType min_offset_needed = std::ceil(
     116           3 :         model.m_transitions.get_time(0) -
     117             :         1); // Need -1 if first time point is integer and just the floor value if not, therefore use ceil and -1
     118           3 :     ScalarType max_offset_needed = std::ceil(model.m_transitions.get_last_time());
     119             : 
     120           3 :     bool min_offset_needed_avail = false;
     121           3 :     bool max_offset_needed_avail = false;
     122             : 
     123             :     // Go through the entries of rki_data and check if date is needed for calculation. Confirmed cases are scaled.
     124             :     // Define dummy variables to store the first and the last index of the TimeSeries where the considered entry of rki_data is potentially needed.
     125           3 :     Eigen::Index idx_needed_first = 0;
     126           3 :     Eigen::Index idx_needed_last  = 0;
     127           3 :     ScalarType time_idx           = 0;
     128          48 :     for (auto&& entry : rki_data) {
     129          45 :         int offset = get_offset_in_days(entry.date, date);
     130          45 :         if ((offset >= min_offset_needed) && (offset <= max_offset_needed)) {
     131          25 :             if (offset == min_offset_needed) {
     132           2 :                 min_offset_needed_avail = true;
     133             :             }
     134          25 :             if (offset == max_offset_needed) {
     135           2 :                 max_offset_needed_avail = true;
     136             :             }
     137             :             // Smallest index for which the entry is needed.
     138          25 :             idx_needed_first =
     139          25 :                 Eigen::Index(std::max(std::floor((offset - model.m_transitions.get_time(0) - 1) / dt), 0.));
     140             :             // Biggest index for which the entry is needed.
     141          25 :             idx_needed_last = Eigen::Index(std::min(std::ceil((offset - model.m_transitions.get_time(0) + 1) / dt),
     142          50 :                                                     double(model.m_transitions.get_num_time_points() - 1)));
     143             : 
     144         134 :             for (Eigen::Index i = idx_needed_first; i <= idx_needed_last; i++) {
     145         109 :                 time_idx = model.m_transitions.get_time(i);
     146         109 :                 if (offset == int(std::floor(time_idx))) {
     147          88 :                     model.m_transitions[i][Eigen::Index(InfectionTransition::InfectedNoSymptomsToInfectedSymptoms)] +=
     148          44 :                         (1 - (time_idx - std::floor(time_idx))) * scale_confirmed_cases * entry.num_confirmed;
     149             :                 }
     150         109 :                 if (offset == int(std::ceil(time_idx))) {
     151          88 :                     model.m_transitions[i][Eigen::Index(InfectionTransition::InfectedNoSymptomsToInfectedSymptoms)] +=
     152          44 :                         (time_idx - std::floor(time_idx)) * scale_confirmed_cases * entry.num_confirmed;
     153             :                 }
     154         109 :                 if (offset == int(std::floor(time_idx - dt))) {
     155          88 :                     model.m_transitions[i][Eigen::Index(InfectionTransition::InfectedNoSymptomsToInfectedSymptoms)] -=
     156          44 :                         (1 - (time_idx - dt - std::floor(time_idx - dt))) * scale_confirmed_cases * entry.num_confirmed;
     157             :                 }
     158         109 :                 if (offset == int(std::ceil(time_idx - dt))) {
     159          88 :                     model.m_transitions[i][Eigen::Index(InfectionTransition::InfectedNoSymptomsToInfectedSymptoms)] -=
     160          44 :                         (time_idx - dt - std::floor(time_idx - dt)) * scale_confirmed_cases * entry.num_confirmed;
     161             :                 }
     162             :             }
     163             : 
     164             :             // Compute Dead by shifting RKI data according to mean stay times.
     165             :             // This is done because the RKI reports death with the date of positive test instead of the date of deaths.
     166          25 :             if (offset == std::floor(-mean_InfectedSymptomsToInfectedSevere - mean_InfectedSevereToInfectedCritical -
     167             :                                      mean_InfectedCriticalToDead)) {
     168           4 :                 model.m_populations[0][Eigen::Index(InfectionState::Dead)] +=
     169           2 :                     (1 - (-mean_InfectedSymptomsToInfectedSevere - mean_InfectedSevereToInfectedCritical -
     170           2 :                           mean_InfectedCriticalToDead -
     171           2 :                           std::floor(-mean_InfectedSymptomsToInfectedSevere - mean_InfectedSevereToInfectedCritical -
     172           2 :                                      mean_InfectedCriticalToDead))) *
     173           2 :                     entry.num_deaths;
     174             :             }
     175          25 :             if (offset == std::ceil(-mean_InfectedSymptomsToInfectedSevere - mean_InfectedSevereToInfectedCritical -
     176             :                                     mean_InfectedCriticalToDead)) {
     177           4 :                 model.m_populations[0][Eigen::Index(InfectionState::Dead)] +=
     178           2 :                     (-mean_InfectedSymptomsToInfectedSevere - mean_InfectedSevereToInfectedCritical -
     179           2 :                      mean_InfectedCriticalToDead -
     180           2 :                      std::floor(-mean_InfectedSymptomsToInfectedSevere - mean_InfectedSevereToInfectedCritical -
     181           2 :                                 mean_InfectedCriticalToDead)) *
     182           2 :                     entry.num_deaths;
     183             :             }
     184             : 
     185          25 :             if (offset == 0) {
     186           3 :                 model.m_total_confirmed_cases = scale_confirmed_cases * entry.num_confirmed;
     187             :             }
     188             :         }
     189             :     }
     190             : 
     191           3 :     if (!max_offset_needed_avail) {
     192           2 :         log_error("Necessary range of dates needed to compute initial values does not exist in RKI data.");
     193           1 :         return failure(StatusCode::OutOfRange, path + ", necessary range of dates does not exist in RKI data.");
     194             :     }
     195             : 
     196           2 :     if (!min_offset_needed_avail) {
     197           1 :         std::string min_date_string =
     198           7 :             std::to_string(min_date.day) + "." + std::to_string(min_date.month) + "." + std::to_string(min_date.year);
     199             :         // Get first date that is needed.
     200           1 :         mio::Date min_offset_date          = offset_date_by_days(date, int(min_offset_needed));
     201           6 :         std::string min_offset_date_string = std::to_string(min_offset_date.day) + "." +
     202           5 :                                              std::to_string(min_offset_date.month) + "." +
     203           3 :                                              std::to_string(min_offset_date.year);
     204           4 :         log_warning("RKI data is needed from " + min_offset_date_string +
     205           3 :                     " to compute initial values. RKI data is only available from " + min_date_string +
     206             :                     ". Missing dates were set to 0.");
     207           1 :     }
     208             : 
     209             :     //--- Calculate the flows "after" InfectedNoSymptomsToInfectedSymptoms. ---
     210             :     // Set m_support_max_vector and m_derivative_vector in the model which is needed for the following computations.
     211           2 :     model.set_transitiondistributions_support_max(dt);
     212           2 :     model.set_transitiondistributions_derivative(dt);
     213             :     // Compute flow InfectedSymptomsToInfectedSevere for -3 * global_support_max, ..., 0.
     214          28 :     for (Eigen::Index i = -3 * global_support_max_index; i <= 0; i++) {
     215          26 :         model.compute_flow(Eigen::Index(InfectionTransition::InfectedSymptomsToInfectedSevere),
     216             :                            Eigen::Index(InfectionTransition::InfectedNoSymptomsToInfectedSymptoms), dt,
     217             :                            i + start_shift);
     218             :     }
     219             :     // Compute flow InfectedSevereToInfectedCritical for -2 * global_support_max, ..., 0.
     220          20 :     for (Eigen::Index i = -2 * global_support_max_index; i <= 0; i++) {
     221          18 :         model.compute_flow(Eigen::Index(InfectionTransition::InfectedSevereToInfectedCritical),
     222             :                            Eigen::Index(InfectionTransition::InfectedSymptomsToInfectedSevere), dt, i + start_shift);
     223             :     }
     224             :     // Compute flows from InfectedSymptoms, InfectedSevere and InfectedCritical to Recovered and
     225             :     // flow InfectedCriticalToDead for -global_support_max, ..., 0.
     226          12 :     for (Eigen::Index i = -global_support_max_index; i <= 0; i++) {
     227             :         // Compute flow InfectedSymptomsToRecovered.
     228          10 :         model.compute_flow(Eigen::Index(InfectionTransition::InfectedSymptomsToRecovered),
     229             :                            Eigen::Index(InfectionTransition::InfectedNoSymptomsToInfectedSymptoms), dt,
     230             :                            i + start_shift);
     231             :         // Compute flow InfectedSevereToRecovered.
     232          10 :         model.compute_flow(Eigen::Index(InfectionTransition::InfectedSevereToRecovered),
     233             :                            Eigen::Index(InfectionTransition::InfectedSymptomsToInfectedSevere), dt, i + start_shift);
     234             :         // Compute flow InfectedCriticalToRecovered.
     235          10 :         model.compute_flow(Eigen::Index(InfectionTransition::InfectedCriticalToRecovered),
     236             :                            Eigen::Index(InfectionTransition::InfectedSevereToInfectedCritical), dt, i + start_shift);
     237             :         // Compute flow InfectedCriticalToDead.
     238          10 :         model.compute_flow(Eigen::Index(InfectionTransition::InfectedCriticalToDead),
     239             :                            Eigen::Index(InfectionTransition::InfectedSevereToInfectedCritical), dt, i + start_shift);
     240             :     }
     241             : 
     242             :     //--- Calculate the remaining flows. ---
     243             :     // Compute flow ExposedToInfectedNoSymptoms for -2 * global_support_max, ..., 0.
     244             :     // Use mean value of the TransitionDistribution InfectedNoSymptomsToInfectedSymptoms for the calculation.
     245           2 :     Eigen::Index index_shift_mean = Eigen::Index(std::round(mean_InfectedNoSymptomsToInfectedSymptoms / dt));
     246          20 :     for (Eigen::Index i = -2 * global_support_max_index; i <= 0; i++) {
     247          18 :         model.m_transitions[i + start_shift][Eigen::Index(InfectionTransition::ExposedToInfectedNoSymptoms)] =
     248          36 :             (1 / model.parameters.get<TransitionProbabilities>()[Eigen::Index(
     249          36 :                      InfectionTransition::InfectedNoSymptomsToInfectedSymptoms)]) *
     250          18 :             model.m_transitions[i + start_shift + index_shift_mean]
     251          36 :                                [Eigen::Index(InfectionTransition::InfectedNoSymptomsToInfectedSymptoms)];
     252             :     }
     253             : 
     254             :     // Compute flow SusceptibleToExposed for -global_support_max, ..., 0.
     255             :     // Use mean values of the TransitionDistribution ExposedToInfectedNoSymptoms and of the TransitionDistribution InfectedNoSymptomsToInfectedSymptoms for the calculation.
     256           2 :     index_shift_mean =
     257           2 :         Eigen::Index(std::round((mean_ExposedToInfectedNoSymptoms + mean_InfectedNoSymptomsToInfectedSymptoms) / dt));
     258          12 :     for (Eigen::Index i = -global_support_max_index; i <= 0; i++) {
     259          10 :         model.m_transitions[i + start_shift][Eigen::Index(InfectionTransition::SusceptibleToExposed)] =
     260          20 :             (1 / model.parameters.get<TransitionProbabilities>()[Eigen::Index(
     261          20 :                      InfectionTransition::InfectedNoSymptomsToInfectedSymptoms)]) *
     262          10 :             model.m_transitions[i + start_shift + index_shift_mean]
     263          20 :                                [Eigen::Index(InfectionTransition::InfectedNoSymptomsToInfectedSymptoms)];
     264             :     }
     265             : 
     266             :     // InfectedNoSymptomsToRecovered for -global_support_max, ..., 0.
     267             :     // If we previously calculated the transition ExposedToInfectedNoSymptoms, we can calculate this transition using the standard formula.
     268          12 :     for (Eigen::Index i = -global_support_max_index; i <= 0; i++) {
     269          10 :         model.compute_flow(Eigen::Index(InfectionTransition::InfectedNoSymptomsToRecovered),
     270             :                            Eigen::Index(InfectionTransition::ExposedToInfectedNoSymptoms), dt, i + start_shift);
     271             :     }
     272             : 
     273             :     // At the end of the calculation, delete all time points that are not required for the simulation.
     274           2 :     auto transition_copy(model.m_transitions);
     275           2 :     model.m_transitions = TimeSeries<ScalarType>((int)InfectionTransition::Count);
     276          12 :     for (Eigen::Index i = -global_support_max_index; i <= 0; i++) {
     277          10 :         model.m_transitions.add_time_point(i * dt, transition_copy.get_value(i + start_shift));
     278             :     }
     279             : 
     280           4 :     return mio::success();
     281           5 : }
     282             : 
     283             : } // namespace isecir
     284             : } // namespace mio
     285             : 
     286             : #endif // MEMILIO_HAS_JSONCPP