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 LCTSECIR_PARAMETERS_IO_H
      22             : #define LCTSECIR_PARAMETERS_IO_H
      23             : 
      24             : #include "memilio/config.h"
      25             : 
      26             : #ifdef MEMILIO_HAS_JSONCPP
      27             : 
      28             : #include "lct_secir/model.h"
      29             : #include "lct_secir/infection_state.h"
      30             : #include "lct_secir/parameters.h"
      31             : #include "memilio/io/epi_data.h"
      32             : #include "memilio/io/io.h"
      33             : #include "memilio/utils/date.h"
      34             : #include "memilio/utils/logging.h"
      35             : #include "memilio/utils/metaprogramming.h"
      36             : #include "memilio/math/eigen.h"
      37             : #include "memilio/math/floating_point.h"
      38             : #include "memilio/epidemiology/lct_populations.h"
      39             : 
      40             : #include <string>
      41             : namespace mio
      42             : {
      43             : namespace lsecir
      44             : {
      45             : 
      46             : namespace details
      47             : { // Use namespace to hide functions that are not intended to be used outside this file.
      48             : 
      49             : /**
      50             : * @brief Processes one entry of an RKI data set for the definition of an initial value vector for an LCT population.
      51             : *   
      52             : * Takes one entry of an RKI data vector and changes the value in populations accordingly.
      53             : * @param[out] populations The populations for which the inital data should be computed and set.
      54             : * @param[in] entry The entry of the RKI data set.
      55             : * @param[in] offset The offset between the date of the entry and the date for which the 
      56             : *   initial value vector is calculated.
      57             : * @param[in] staytimes A vector of the average time spent in each compartment defined in InfectionState
      58             : *    (for the group under consideration) for which the initial value vector is calculated.
      59             : * @param[in] inv_prob_SymptomsPerNoSymptoms The inverse of the probability InfectedSymptomsPerInfectedNoSymptoms
      60             : *    (for the group under consideration) for which the initial value vector is calculated.
      61             : * @param[in] severePerInfectedSymptoms Probability (for the group under consideration)
      62             : *    for which the initial value vector is calculated.
      63             : * @param[in] criticalPerSevere Probability (for the group under consideration)
      64             : *    for which the initial value vector is calculated.
      65             : * @param[in] scale_confirmed_cases Factor by which to scale the confirmed cases of RKI data to consider unreported cases.
      66             : * @tparam Populations is expected to be an LctPopulations defined in epidemiology/lct_populations. 
      67             : *   This defined the number of age groups and the number of subcompartments used.
      68             : * @tparam EntryType The type of the data entry of the RKI data.
      69             : * @tparam Group The age group of the entry the should be processed.
      70             : */
      71             : template <class Populations, class EntryType, size_t Group>
      72         194 : void process_entry(Populations& populations, const EntryType& entry, int offset,
      73             :                    const std::vector<ScalarType> staytimes, ScalarType inv_prob_SymptomsPerNoSymptoms,
      74             :                    ScalarType severePerInfectedSymptoms, ScalarType criticalPerSevere, ScalarType scale_confirmed_cases)
      75             : {
      76             :     using LctStateGroup      = type_at_index_t<Group, typename Populations::LctStatesGroups>;
      77         194 :     size_t first_index_group = populations.template get_first_index_of_group<Group>();
      78             : 
      79             :     // Add confirmed cases at date to compartment Recovered.
      80         194 :     if (offset == 0) {
      81          36 :         populations[first_index_group + LctStateGroup::template get_first_index<InfectionState::Recovered>()] +=
      82          18 :             scale_confirmed_cases * entry.num_confirmed;
      83             :     }
      84             : 
      85             :     // Compute initial values for compartment InfectedNoSymptoms.
      86         194 :     if (offset >= 0 && offset <= std::ceil(staytimes[(size_t)InfectionState::InfectedNoSymptoms])) {
      87             :         // Mean stay time in each subcompartment.
      88          54 :         ScalarType time_INS_per_subcomp =
      89          54 :             staytimes[(size_t)InfectionState::InfectedNoSymptoms] /
      90             :             (ScalarType)LctStateGroup::template get_num_subcompartments<InfectionState::InfectedNoSymptoms>();
      91             :         // Index of the last subcompartment of InfectedNoSymptoms.
      92          54 :         size_t last_index_INS =
      93          54 :             first_index_group + LctStateGroup::template get_first_index<InfectionState::InfectedNoSymptoms>() +
      94          54 :             LctStateGroup::template get_num_subcompartments<InfectionState::InfectedNoSymptoms>() - 1;
      95         198 :         for (int i = 0; i < (int)LctStateGroup::template get_num_subcompartments<InfectionState::InfectedNoSymptoms>();
      96             :              i++) {
      97         144 :             if (offset == std::floor(i * time_INS_per_subcomp)) {
      98          96 :                 populations[last_index_INS - i] -=
      99          48 :                     (1 - (i * time_INS_per_subcomp - std::floor(i * time_INS_per_subcomp))) *
     100          48 :                     inv_prob_SymptomsPerNoSymptoms * scale_confirmed_cases * entry.num_confirmed;
     101             :             }
     102         144 :             if (offset == std::ceil(i * time_INS_per_subcomp)) {
     103          96 :                 populations[last_index_INS - i] -= (i * time_INS_per_subcomp - std::floor(i * time_INS_per_subcomp)) *
     104          96 :                                                    inv_prob_SymptomsPerNoSymptoms * scale_confirmed_cases *
     105          48 :                                                    entry.num_confirmed;
     106             :             }
     107         144 :             if (offset == std::floor((i + 1) * time_INS_per_subcomp)) {
     108          96 :                 populations[last_index_INS - i] +=
     109          48 :                     (1 - ((i + 1) * time_INS_per_subcomp - std::floor((i + 1) * time_INS_per_subcomp))) *
     110          48 :                     inv_prob_SymptomsPerNoSymptoms * scale_confirmed_cases * entry.num_confirmed;
     111             :             }
     112         144 :             if (offset == std::ceil((i + 1) * time_INS_per_subcomp)) {
     113          96 :                 populations[last_index_INS - i] +=
     114          48 :                     ((i + 1) * time_INS_per_subcomp - std::floor((i + 1) * time_INS_per_subcomp)) *
     115          48 :                     inv_prob_SymptomsPerNoSymptoms * scale_confirmed_cases * entry.num_confirmed;
     116             :             }
     117             :         }
     118             :     }
     119             : 
     120             :     // Compute initial values for compartment Exposed.
     121         262 :     if (offset >= std::floor(staytimes[(size_t)InfectionState::InfectedNoSymptoms]) &&
     122          68 :         offset <= std::ceil(staytimes[(size_t)InfectionState::Exposed] +
     123             :                             staytimes[(size_t)InfectionState::InfectedNoSymptoms])) {
     124             :         // Mean stay time in each subcompartment.
     125          68 :         ScalarType time_Exposed_per_subcomp =
     126          68 :             staytimes[(size_t)InfectionState::Exposed] /
     127             :             (ScalarType)LctStateGroup::template get_num_subcompartments<InfectionState::Exposed>();
     128             :         // Index of the last subcompartment of Exposed.
     129         136 :         size_t last_index_Exposed = first_index_group +
     130          68 :                                     LctStateGroup::template get_first_index<InfectionState::Exposed>() +
     131          68 :                                     LctStateGroup::template get_num_subcompartments<InfectionState::Exposed>() - 1;
     132         192 :         for (int i = 0; i < (int)LctStateGroup::template get_num_subcompartments<InfectionState::Exposed>(); i++) {
     133         248 :             if (offset ==
     134         124 :                 std::floor(staytimes[(size_t)InfectionState::InfectedNoSymptoms] + i * time_Exposed_per_subcomp)) {
     135          33 :                 populations[last_index_Exposed - i] -=
     136          33 :                     (1 - (staytimes[(size_t)InfectionState::InfectedNoSymptoms] + i * time_Exposed_per_subcomp -
     137          33 :                           std::floor(staytimes[(size_t)InfectionState::InfectedNoSymptoms] +
     138          33 :                                      i * time_Exposed_per_subcomp))) *
     139          33 :                     inv_prob_SymptomsPerNoSymptoms * scale_confirmed_cases * entry.num_confirmed;
     140             :             }
     141         248 :             if (offset ==
     142         124 :                 std::ceil(staytimes[(size_t)InfectionState::InfectedNoSymptoms] + i * time_Exposed_per_subcomp)) {
     143          31 :                 populations[last_index_Exposed - i] -=
     144          31 :                     (staytimes[(size_t)InfectionState::InfectedNoSymptoms] + i * time_Exposed_per_subcomp -
     145          31 :                      std::floor(staytimes[(size_t)InfectionState::InfectedNoSymptoms] + i * time_Exposed_per_subcomp)) *
     146          31 :                     inv_prob_SymptomsPerNoSymptoms * scale_confirmed_cases * entry.num_confirmed;
     147             :             }
     148         124 :             if (offset == std::floor(staytimes[(size_t)InfectionState::InfectedNoSymptoms] +
     149         124 :                                      (i + 1) * time_Exposed_per_subcomp)) {
     150          31 :                 populations[last_index_Exposed - i] +=
     151          31 :                     (1 - (staytimes[(size_t)InfectionState::InfectedNoSymptoms] + (i + 1) * time_Exposed_per_subcomp -
     152          31 :                           std::floor(staytimes[(size_t)InfectionState::InfectedNoSymptoms] +
     153          31 :                                      (i + 1) * time_Exposed_per_subcomp))) *
     154          31 :                     inv_prob_SymptomsPerNoSymptoms * scale_confirmed_cases * entry.num_confirmed;
     155             :             }
     156         248 :             if (offset ==
     157         124 :                 std::ceil(staytimes[(size_t)InfectionState::InfectedNoSymptoms] + (i + 1) * time_Exposed_per_subcomp)) {
     158          29 :                 populations[last_index_Exposed - i] +=
     159          29 :                     (staytimes[(size_t)InfectionState::InfectedNoSymptoms] + (i + 1) * time_Exposed_per_subcomp -
     160          29 :                      std::floor(staytimes[(size_t)InfectionState::InfectedNoSymptoms] +
     161          29 :                                 (i + 1) * time_Exposed_per_subcomp)) *
     162          29 :                     inv_prob_SymptomsPerNoSymptoms * scale_confirmed_cases * entry.num_confirmed;
     163             :             }
     164             :         }
     165             :     }
     166             : 
     167             :     // Compute initial values for compartment InfectedSymptoms.
     168         194 :     if (offset >= std::floor(-staytimes[(size_t)InfectionState::InfectedSymptoms]) && offset <= 0) {
     169             :         // Mean stay time in each subcompartment.
     170          72 :         ScalarType time_IS_per_subcomp =
     171          72 :             staytimes[(size_t)InfectionState::InfectedSymptoms] /
     172             :             (ScalarType)LctStateGroup::template get_num_subcompartments<InfectionState::InfectedSymptoms>();
     173             :         // Index of the first subcompartment of InfectedSymptoms.
     174          72 :         size_t first_index_IS =
     175          72 :             first_index_group + LctStateGroup::template get_first_index<InfectionState::InfectedSymptoms>();
     176         204 :         for (int i = 0; i < (int)LctStateGroup::template get_num_subcompartments<InfectionState::InfectedSymptoms>();
     177             :              i++) {
     178         132 :             if (offset == std::floor(-time_IS_per_subcomp * (i + 1))) {
     179          66 :                 populations[first_index_IS + i] -=
     180          33 :                     (1 - (-(i + 1.) * time_IS_per_subcomp - std::floor(-(i + 1) * time_IS_per_subcomp))) *
     181          33 :                     scale_confirmed_cases * entry.num_confirmed;
     182             :             }
     183         132 :             if (offset == std::ceil(-time_IS_per_subcomp * (i + 1))) {
     184          66 :                 populations[first_index_IS + i] -=
     185          33 :                     (-(i + 1) * time_IS_per_subcomp - std::floor(-(i + 1) * time_IS_per_subcomp)) *
     186          33 :                     scale_confirmed_cases * entry.num_confirmed;
     187             :             }
     188         132 :             if (offset == std::floor(-time_IS_per_subcomp * i)) {
     189          66 :                 populations[first_index_IS + i] +=
     190          66 :                     (1 - (-i * time_IS_per_subcomp - std::floor(-i * time_IS_per_subcomp))) * scale_confirmed_cases *
     191          33 :                     entry.num_confirmed;
     192             :             }
     193         132 :             if (offset == std::ceil(-time_IS_per_subcomp * i)) {
     194          66 :                 populations[first_index_IS + i] += (-i * time_IS_per_subcomp - std::floor(-i * time_IS_per_subcomp)) *
     195          33 :                                                    scale_confirmed_cases * entry.num_confirmed;
     196             :             }
     197             :         }
     198             :     }
     199             : 
     200             :     // Compute initial values for compartment InfectedSevere.
     201         194 :     if (offset >= std::floor(-staytimes[(size_t)InfectionState::InfectedSymptoms] -
     202         370 :                              staytimes[(size_t)InfectionState::InfectedSevere]) &&
     203         176 :         offset <= std::ceil(-staytimes[(size_t)InfectionState::InfectedSymptoms])) {
     204             :         // Mean stay time in each subcompartment.
     205          72 :         ScalarType time_ISev_per_subcomp =
     206          72 :             staytimes[(size_t)InfectionState::InfectedSevere] /
     207             :             (ScalarType)LctStateGroup::template get_num_subcompartments<InfectionState::InfectedSevere>();
     208             :         // Index of the first subcompartment of InfectedSevere.
     209          72 :         size_t first_index_ISev =
     210          72 :             first_index_group + LctStateGroup::template get_first_index<InfectionState::InfectedSevere>();
     211         204 :         for (int i = 0; i < (int)LctStateGroup::template get_num_subcompartments<InfectionState::InfectedSevere>();
     212             :              i++) {
     213         264 :             if (offset ==
     214         132 :                 std::floor(-staytimes[(size_t)InfectionState::InfectedSymptoms] - time_ISev_per_subcomp * (i + 1))) {
     215          33 :                 populations[first_index_ISev + i] -=
     216          33 :                     severePerInfectedSymptoms *
     217          33 :                     (1 - (-staytimes[(size_t)InfectionState::InfectedSymptoms] - (i + 1) * time_ISev_per_subcomp -
     218          33 :                           std::floor(-staytimes[(size_t)InfectionState::InfectedSymptoms] -
     219          33 :                                      (i + 1) * time_ISev_per_subcomp))) *
     220          33 :                     scale_confirmed_cases * entry.num_confirmed;
     221             :             }
     222         264 :             if (offset ==
     223         132 :                 std::ceil(-staytimes[(size_t)InfectionState::InfectedSymptoms] - time_ISev_per_subcomp * (i + 1))) {
     224          33 :                 populations[first_index_ISev + i] -=
     225          33 :                     severePerInfectedSymptoms *
     226          33 :                     (-staytimes[(size_t)InfectionState::InfectedSymptoms] - (i + 1) * time_ISev_per_subcomp -
     227          33 :                      std::floor(-staytimes[(size_t)InfectionState::InfectedSymptoms] -
     228          33 :                                 (i + 1) * time_ISev_per_subcomp)) *
     229          33 :                     scale_confirmed_cases * entry.num_confirmed;
     230             :             }
     231         264 :             if (offset ==
     232         132 :                 std::floor(-staytimes[(size_t)InfectionState::InfectedSymptoms] - time_ISev_per_subcomp * i)) {
     233          33 :                 populations[first_index_ISev + i] +=
     234          33 :                     severePerInfectedSymptoms *
     235          33 :                     (1 -
     236          33 :                      (-staytimes[(size_t)InfectionState::InfectedSymptoms] - i * time_ISev_per_subcomp -
     237          33 :                       std::floor(-staytimes[(size_t)InfectionState::InfectedSymptoms] - i * time_ISev_per_subcomp))) *
     238          33 :                     scale_confirmed_cases * entry.num_confirmed;
     239             :             }
     240         132 :             if (offset == std::ceil(-staytimes[(size_t)InfectionState::InfectedSymptoms] - time_ISev_per_subcomp * i)) {
     241          33 :                 populations[first_index_ISev + i] +=
     242          33 :                     severePerInfectedSymptoms *
     243          33 :                     (-staytimes[(size_t)InfectionState::InfectedSymptoms] - i * time_ISev_per_subcomp -
     244          33 :                      std::floor(-staytimes[(size_t)InfectionState::InfectedSymptoms] - i * time_ISev_per_subcomp)) *
     245          33 :                     scale_confirmed_cases * entry.num_confirmed;
     246             :             }
     247             :         }
     248             :     }
     249             : 
     250             :     // Compute initial values for compartment InfectedCritical.
     251         388 :     if (offset >= std::floor(-staytimes[(size_t)InfectionState::InfectedSymptoms] -
     252         194 :                              staytimes[(size_t)InfectionState::InfectedSevere] -
     253         388 :                              staytimes[(size_t)InfectionState::InfectedCritical]) &&
     254         194 :         offset <= std::ceil(-staytimes[(size_t)InfectionState::InfectedSymptoms] -
     255             :                             staytimes[(size_t)InfectionState::InfectedSevere])) {
     256             :         // Mean stay time in each subcompartment.
     257          54 :         ScalarType time_ICri_per_subcomp =
     258          54 :             staytimes[(size_t)InfectionState::InfectedCritical] /
     259             :             (ScalarType)LctStateGroup::template get_num_subcompartments<InfectionState::InfectedCritical>();
     260             :         // Index of the first subcompartment of InfectedCritical.
     261          54 :         size_t first_index_ICri =
     262          54 :             first_index_group + LctStateGroup::template get_first_index<InfectionState::InfectedCritical>();
     263         141 :         for (int i = 0; i < (int)LctStateGroup::template get_num_subcompartments<InfectionState::InfectedCritical>();
     264             :              i++) {
     265         174 :             if (offset ==
     266          87 :                 std::floor(-staytimes[(size_t)InfectionState::InfectedSymptoms] -
     267          87 :                            staytimes[(size_t)InfectionState::InfectedSevere] - time_ICri_per_subcomp * (i + 1))) {
     268          29 :                 populations[first_index_ICri + i] -=
     269          58 :                     severePerInfectedSymptoms * criticalPerSevere *
     270          29 :                     (1 - (-staytimes[(size_t)InfectionState::InfectedSymptoms] -
     271          58 :                           staytimes[(size_t)InfectionState::InfectedSevere] - (i + 1) * time_ICri_per_subcomp -
     272          29 :                           std::floor(-staytimes[(size_t)InfectionState::InfectedSymptoms] -
     273             :                                      staytimes[(size_t)InfectionState::InfectedSevere] -
     274          29 :                                      (i + 1) * time_ICri_per_subcomp))) *
     275          29 :                     scale_confirmed_cases * entry.num_confirmed;
     276             :             }
     277         174 :             if (offset ==
     278          87 :                 std::ceil(-staytimes[(size_t)InfectionState::InfectedSymptoms] -
     279          87 :                           staytimes[(size_t)InfectionState::InfectedSevere] - time_ICri_per_subcomp * (i + 1))) {
     280          29 :                 populations[first_index_ICri + i] -=
     281          58 :                     severePerInfectedSymptoms * criticalPerSevere *
     282          29 :                     (-staytimes[(size_t)InfectionState::InfectedSymptoms] -
     283          58 :                      staytimes[(size_t)InfectionState::InfectedSevere] - (i + 1) * time_ICri_per_subcomp -
     284          29 :                      std::floor(-staytimes[(size_t)InfectionState::InfectedSymptoms] -
     285          29 :                                 staytimes[(size_t)InfectionState::InfectedSevere] - (i + 1) * time_ICri_per_subcomp)) *
     286          29 :                     scale_confirmed_cases * entry.num_confirmed;
     287             :             }
     288          87 :             if (offset == std::floor(-staytimes[(size_t)InfectionState::InfectedSymptoms] -
     289          87 :                                      staytimes[(size_t)InfectionState::InfectedSevere] - time_ICri_per_subcomp * i)) {
     290          29 :                 populations[first_index_ICri + i] +=
     291          58 :                     severePerInfectedSymptoms * criticalPerSevere *
     292          29 :                     (1 - (-staytimes[(size_t)InfectionState::InfectedSymptoms] -
     293          58 :                           staytimes[(size_t)InfectionState::InfectedSevere] - i * time_ICri_per_subcomp -
     294          29 :                           std::floor(-staytimes[(size_t)InfectionState::InfectedSymptoms] -
     295          29 :                                      staytimes[(size_t)InfectionState::InfectedSevere] - i * time_ICri_per_subcomp))) *
     296          29 :                     scale_confirmed_cases * entry.num_confirmed;
     297             :             }
     298          87 :             if (offset == std::ceil(-staytimes[(size_t)InfectionState::InfectedSymptoms] -
     299          87 :                                     staytimes[(size_t)InfectionState::InfectedSevere] - time_ICri_per_subcomp * i)) {
     300          29 :                 populations[first_index_ICri + i] +=
     301          58 :                     severePerInfectedSymptoms * criticalPerSevere *
     302          29 :                     (-staytimes[(size_t)InfectionState::InfectedSymptoms] -
     303          58 :                      staytimes[(size_t)InfectionState::InfectedSevere] - i * time_ICri_per_subcomp -
     304          29 :                      std::floor(-staytimes[(size_t)InfectionState::InfectedSymptoms] -
     305          29 :                                 staytimes[(size_t)InfectionState::InfectedSevere] - i * time_ICri_per_subcomp)) *
     306          29 :                     scale_confirmed_cases * entry.num_confirmed;
     307             :             }
     308             :         }
     309             :     }
     310             : 
     311             :     // Compute Dead.
     312         388 :     if (offset == std::floor(-staytimes[(size_t)InfectionState::InfectedSymptoms] -
     313         194 :                              staytimes[(size_t)InfectionState::InfectedSevere] -
     314             :                              staytimes[(size_t)InfectionState::InfectedCritical])) {
     315          18 :         populations[first_index_group + LctStateGroup::template get_first_index<InfectionState::Dead>()] +=
     316          18 :             (1 -
     317          18 :              (-staytimes[(size_t)InfectionState::InfectedSymptoms] - staytimes[(size_t)InfectionState::InfectedSevere] -
     318          18 :               staytimes[(size_t)InfectionState::InfectedCritical] -
     319          36 :               std::floor(-staytimes[(size_t)InfectionState::InfectedSymptoms] -
     320          18 :                          staytimes[(size_t)InfectionState::InfectedSevere] -
     321          18 :                          staytimes[(size_t)InfectionState::InfectedCritical]))) *
     322          18 :             entry.num_deaths;
     323             :     }
     324         388 :     if (offset == std::ceil(-staytimes[(size_t)InfectionState::InfectedSymptoms] -
     325         194 :                             staytimes[(size_t)InfectionState::InfectedSevere] -
     326             :                             staytimes[(size_t)InfectionState::InfectedCritical])) {
     327          18 :         populations[first_index_group + LctStateGroup::template get_first_index<InfectionState::Dead>()] +=
     328          18 :             (-staytimes[(size_t)InfectionState::InfectedSymptoms] - staytimes[(size_t)InfectionState::InfectedSevere] -
     329          18 :              staytimes[(size_t)InfectionState::InfectedCritical] -
     330          36 :              std::floor(-staytimes[(size_t)InfectionState::InfectedSymptoms] -
     331          18 :                         staytimes[(size_t)InfectionState::InfectedSevere] -
     332          18 :                         staytimes[(size_t)InfectionState::InfectedCritical])) *
     333          18 :             entry.num_deaths;
     334             :     }
     335         194 : }
     336             : 
     337             : /**
     338             : * @brief Computes an initialization vector for an LCT population with case data from RKI recursively for each age group
     339             : *    (or for one age group in the case without age resolution).
     340             : *   
     341             : * Please use the set_initial_data_from_confirmed_cases() function, which calls this function automatically!
     342             : * This function calculates a segment referring to the defined age group of the initial value vector with 
     343             : *   subcompartments using the rki_data and the parameters.
     344             : * The values for the whole initial value vector stored in populations are calculated recursively.
     345             : *
     346             : * @param[in] rki_data Vector with the RKI data.
     347             : * @param[out] populations The populations for which the inital data should be computed and set.
     348             : * @param[in] parameters The parameters that should be used to calculate the initial values. 
     349             : *   Probabilities and mean stay times are used.
     350             : * @param[in] date Date for which the initial values should be computed. date is the start date of the simulation.
     351             : * @param[in] total_population Total size of the population of Germany or of every age group. 
     352             : * @param[in] scale_confirmed_cases Factor(s for each age group) by which to scale the confirmed cases of the rki data 
     353             : *   to consider unreported cases.
     354             : * @tparam Populations is expected to be an LctPopulations defined in epidemiology/lct_populations. 
     355             : *   This defined the number of age groups and the number of subcompartments used.
     356             : * @tparam EntryType is expected to be ConfirmedCasesNoAgeEntry for data that is not age resolved and 
     357             : *   ConfirmedCasesDataEntry for age resolved data. See also epi_data.h.
     358             : * @tparam Group The age group for which the initial values should be calculated. The function is called recursively 
     359             : *   such that the initial values are calculated for every age group if Group is zero at the beginning.
     360             : * @returns Any io errors that happen during data processing.
     361             : */
     362             : template <class Populations, class EntryType, size_t Group = 0>
     363          18 : IOResult<void> set_initial_data_from_confirmed_cases_impl(Populations& populations,
     364             :                                                           const std::vector<EntryType>& rki_data,
     365             :                                                           const Parameters& parameters, Date date,
     366             :                                                           const std::vector<ScalarType>& total_population,
     367             :                                                           const std::vector<ScalarType>& scale_confirmed_cases)
     368             : {
     369             :     static_assert((Group < Populations::num_groups) && (Group >= 0), "The template parameter Group should be valid.");
     370             :     using LctStateGroup      = type_at_index_t<Group, typename Populations::LctStatesGroups>;
     371          18 :     size_t first_index_group = populations.template get_first_index_of_group<Group>();
     372             : 
     373             :     // Define variables for parameters.
     374          36 :     std::vector<ScalarType> staytimes((size_t)InfectionState::Count, -1.);
     375          18 :     staytimes[(size_t)InfectionState::Exposed]            = parameters.template get<TimeExposed>()[Group];
     376          18 :     staytimes[(size_t)InfectionState::InfectedNoSymptoms] = parameters.template get<TimeInfectedNoSymptoms>()[Group];
     377          18 :     staytimes[(size_t)InfectionState::InfectedSymptoms]   = parameters.template get<TimeInfectedSymptoms>()[Group];
     378          18 :     staytimes[(size_t)InfectionState::InfectedSevere]     = parameters.template get<TimeInfectedSevere>()[Group];
     379          18 :     staytimes[(size_t)InfectionState::InfectedCritical]   = parameters.template get<TimeInfectedCritical>()[Group];
     380          18 :     ScalarType inv_prob_SymptomsPerNoSymptoms =
     381          18 :         1 / (1 - parameters.template get<RecoveredPerInfectedNoSymptoms>()[Group]);
     382          18 :     ScalarType prob_SeverePerInfectedSymptoms = parameters.template get<SeverePerInfectedSymptoms>()[Group];
     383          18 :     ScalarType prob_CriticalPerSevere         = parameters.template get<CriticalPerSevere>()[Group];
     384             : 
     385          36 :     ScalarType min_offset_needed = std::floor(-staytimes[(size_t)InfectionState::InfectedSymptoms] -
     386          18 :                                               staytimes[(size_t)InfectionState::InfectedSevere] -
     387             :                                               staytimes[(size_t)InfectionState::InfectedCritical]);
     388             :     ScalarType max_offset_needed =
     389          18 :         std::ceil(staytimes[(size_t)InfectionState::Exposed] + staytimes[(size_t)InfectionState::InfectedNoSymptoms]);
     390             : 
     391          18 :     bool min_offset_needed_avail = false;
     392          18 :     bool max_offset_needed_avail = false;
     393             : 
     394             :     // Go through the entries of rki_data and check if the entry is age resolved and is referring to
     395             :     // the age_group Group in the case with age resolution. If the date is
     396             :     // needed for calculation, another function to handle the entry is called. Confirmed cases are scaled by scale_confirmed_cases.
     397         986 :     for (auto&& entry : rki_data) {
     398             :         if constexpr (std::is_same_v<EntryType, ConfirmedCasesDataEntry>) {
     399         924 :             if (!((size_t)entry.age_group == Group)) {
     400         770 :                 continue;
     401             :             }
     402             :         }
     403         198 :         int offset = get_offset_in_days(entry.date, date);
     404         198 :         if ((offset >= min_offset_needed) && (offset <= max_offset_needed)) {
     405         194 :             if (offset == max_offset_needed) {
     406          16 :                 max_offset_needed_avail = true;
     407             :             }
     408         194 :             if (offset == min_offset_needed) {
     409          18 :                 min_offset_needed_avail = true;
     410             :             }
     411         194 :             process_entry<Populations, EntryType, Group>(populations, entry, offset, staytimes,
     412             :                                                          inv_prob_SymptomsPerNoSymptoms, prob_SeverePerInfectedSymptoms,
     413             :                                                          prob_CriticalPerSevere, scale_confirmed_cases[Group]);
     414             :         }
     415             :     }
     416             : 
     417             :     // Compute Recovered.
     418          18 :     populations[first_index_group + LctStateGroup::template get_first_index<InfectionState::Recovered>()] -=
     419             :         populations.get_compartments()
     420          36 :             .segment(first_index_group + LctStateGroup::template get_first_index<InfectionState::InfectedSymptoms>(),
     421          18 :                      LctStateGroup::template get_num_subcompartments<InfectionState::InfectedSymptoms>() +
     422          18 :                          LctStateGroup::template get_num_subcompartments<InfectionState::InfectedSevere>() +
     423          18 :                          LctStateGroup::template get_num_subcompartments<InfectionState::InfectedCritical>())
     424          54 :             .sum();
     425          18 :     populations[first_index_group + LctStateGroup::template get_first_index<InfectionState::Recovered>()] -=
     426          18 :         populations[first_index_group + LctStateGroup::template get_first_index<InfectionState::Dead>()];
     427             : 
     428             :     // Compute Susceptibles.
     429          36 :     populations[first_index_group + LctStateGroup::template get_first_index<InfectionState::Susceptible>()] =
     430          18 :         total_population[Group] -
     431          36 :         populations.get_compartments().segment(first_index_group + 1, LctStateGroup::Count - 1).sum();
     432             : 
     433             :     // Check whether all necessary dates are available.
     434          18 :     if (!max_offset_needed_avail || !min_offset_needed_avail) {
     435           4 :         log_error(
     436             :             "Necessary range of dates needed to compute initial values does not exist in RKI data for group {:d}.",
     437             :             Group);
     438             :         return failure(StatusCode::OutOfRange,
     439           2 :                        "Necessary range of dates needed to compute initial values does not exist in RKI data.");
     440             :     }
     441             : 
     442             :     // Check if all values for populations are valid.
     443          70 :     for (size_t i = first_index_group; i < LctStateGroup::Count; i++) {
     444          56 :         if (floating_point_less((ScalarType)populations[i], 0., 1e-14)) {
     445           4 :             log_error("Something went wrong in the initialization of group {:d}. At least one entry is negative.",
     446             :                       Group);
     447             :             return failure(StatusCode::InvalidValue,
     448           2 :                            "Something went wrong in the initialization as at least one entry is negative.");
     449             :         }
     450             :     }
     451             : 
     452             :     if constexpr (Group + 1 < Populations::num_groups) {
     453             :         return set_initial_data_from_confirmed_cases_impl<Populations, EntryType, Group + 1>(
     454          10 :             populations, rki_data, parameters, date, total_population, scale_confirmed_cases);
     455             :     }
     456             :     else {
     457           8 :         return mio::success();
     458             :     }
     459          18 : }
     460             : } // namespace details
     461             : 
     462             : /**
     463             : * @brief Computes an initialization vector for an LCT population with case data from RKI.
     464             : *   
     465             : * Use just one group in the definition of the populations to not divide between age groups.
     466             : * Otherwise, the number of groups has to match the number of RKI age groups.
     467             : * The function calculates an initial value vector referring to an LCT population and updates the initial value vector
     468             : * in the populations class.
     469             : * For the computation expected stay times in the subcompartments defined in the parameters variable are used.
     470             : * To calculate the initial values, we assume for simplicity that individuals stay in the subcompartment 
     471             : * for exactly the expected time.
     472             : * The RKI data are linearly interpolated within one day to match the expected stay time in a subcompartment.
     473             : * The RKI data should contain data for each needed day with or without division of age groups, 
     474             : *   the completeness of the dates is not verified.
     475             : * Data can be downloaded e.g. with the file pycode/memilio-epidata/memilio/epidata/getCaseData.py, which creates files
     476             : * named e.g. cases_all_germany.json for no groups or cases_all_age.json with division in age groups or similar names.
     477             : * One should set impute_dates=True so that missing dates are imputed. 
     478             : * To read the data into a vector, use the functionality from epi_data.h.
     479             : * The data and the number of entries in the total_population and scale_confirmed_cases vectors have to match the 
     480             : *   number of groups used in Populations.
     481             : *
     482             : * @param[in] rki_data Vector with the RKI data.
     483             : * @param[out] populations The populations for which the inital data should be computed and set.
     484             : * @param[in] parameters The parameters that should be used to calculate the initial values. 
     485             : *   Probabilities and mean stay times are used.
     486             : * @param[in] date Date for which the initial values should be computed. date is the start date of the simulation.
     487             : * @param[in] total_population Total size of the population of Germany or of every age group. 
     488             : * @param[in] scale_confirmed_cases Factor(s for each age group) by which to scale the confirmed cases of the rki data 
     489             : *   to consider unreported cases.
     490             : * @tparam Populations is expected to be an LctPopulations defined in epidemiology/lct_populations. 
     491             : *   This defined the number of age groups and the number of subcompartments used.
     492             : * @tparam EntryType is expected to be ConfirmedCasesNoAgeEntry for data that is not age resolved and 
     493             : *   ConfirmedCasesDataEntry for age resolved data. See also epi_data.h.
     494             : * @returns Any io errors that happen during data processing.
     495             : */
     496             : template <class Populations, class EntryType>
     497          12 : IOResult<void> set_initial_data_from_confirmed_cases(const std::vector<EntryType>& rki_data, Populations& populations,
     498             :                                                      const Parameters& parameters, const Date date,
     499             :                                                      const std::vector<ScalarType>& total_population,
     500             :                                                      const std::vector<ScalarType>& scale_confirmed_cases)
     501             : { // Check if the inputs are matching.
     502          12 :     assert(total_population.size() == Populations::num_groups);
     503          12 :     assert(scale_confirmed_cases.size() == Populations::num_groups);
     504             :     if constexpr (Populations::num_groups > 1) {
     505             :         static_assert(std::is_same_v<EntryType, ConfirmedCasesDataEntry>);
     506           6 :         assert(ConfirmedCasesDataEntry::age_group_names.size() == Populations::num_groups);
     507             :     }
     508             :     else {
     509             :         static_assert(std::is_same_v<EntryType, ConfirmedCasesNoAgeEntry>);
     510             :     }
     511             :     // Check if RKI data vector is valid.
     512         387 :     auto max_date_entry = std::max_element(rki_data.begin(), rki_data.end(), [](auto&& a, auto&& b) {
     513         375 :         return a.date < b.date;
     514             :     });
     515          12 :     if (max_date_entry == rki_data.end()) {
     516           4 :         log_error("RKI data file is empty.");
     517           2 :         return failure(StatusCode::InvalidFileFormat, "RKI data is empty.");
     518             :     }
     519          10 :     auto max_date = max_date_entry->date;
     520          10 :     if (max_date < date) {
     521           4 :         log_error("Specified date does not exist in RKI data.");
     522           2 :         return failure(StatusCode::OutOfRange, "Specified date does not exist in RKI data.");
     523             :     }
     524             :     // Initially set populations to zero.
     525         378 :     for (size_t i = 0; i < populations.get_num_compartments(); i++) {
     526         370 :         populations[i] = 0;
     527             :     }
     528             :     return details::set_initial_data_from_confirmed_cases_impl<Populations, EntryType>(
     529           8 :         populations, rki_data, parameters, date, total_population, scale_confirmed_cases);
     530             : }
     531             : 
     532             : } // namespace lsecir
     533             : } // namespace mio
     534             : #endif // MEMILIO_HAS_JSONCPP
     535             : 
     536             : #endif // LCTSECIR_PARAMETERS_IO_H