Line data    Source code

       1             : /* 
       2             : * Copyright (C) 2020-2024 MEmilio
       3             : *
       4             : * Authors: Daniel Abele, Martin J. Kuehn
       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             : #ifndef ODESECIR_PARAMETER_SPACE_H
      21             : #define ODESECIR_PARAMETER_SPACE_H
      22             : 
      23             : #include "memilio/mobility/metapopulation_mobility_instant.h"
      24             : #include "memilio/utils/memory.h"
      25             : #include "memilio/utils/logging.h"
      26             : #include "memilio/utils/parameter_distributions.h"
      27             : #include "ode_secir/model.h"
      28             : 
      29             : #include <assert.h>
      30             : #include <string>
      31             : #include <vector>
      32             : #include <random>
      33             : #include <memory>
      34             : 
      35             : namespace mio
      36             : {
      37             : namespace osecir
      38             : {
      39             : /** Sets alls SecirParams parameters normally distributed, 
      40             :  *  using the current value and a given standard deviation
      41             :  * @param[inout] params SecirParams including contact patterns for alle age groups
      42             :  * @param[in] t0 start time
      43             :  * @param[in] tmax end time
      44             :  * @param[in] dev_rel maximum relative deviation from particular value(s) given in params
      45             :  */
      46             : template <typename FP = double>
      47          41 : void set_params_distributions_normal(Model<FP>& model, double t0, double tmax, double dev_rel)
      48             : {
      49        6137 :     auto set_distribution = [dev_rel](UncertainValue<FP>& v, double min_val = 0.001) {
      50       10160 :         v.set_distribution(ParameterDistributionNormal(
      51             :             //add add limits for nonsense big values. Also mscv has a problem with a few doubles so this fixes it
      52        6096 :             std::min(std::max(min_val, (1 - dev_rel * 2.6) * v), 0.1 * std::numeric_limits<double>::max()),
      53        6096 :             std::min(std::max(min_val, (1 + dev_rel * 2.6) * v), 0.5 * std::numeric_limits<double>::max()),
      54        4064 :             std::min(std::max(min_val, double(v)), 0.3 * std::numeric_limits<double>::max()),
      55        4064 :             std::min(std::max(min_val, dev_rel * v), std::numeric_limits<double>::max())));
      56             :     };
      57             : 
      58          41 :     set_distribution(model.parameters.template get<Seasonality<FP>>(), 0.0);
      59          41 :     set_distribution(model.parameters.template get<ICUCapacity<FP>>());
      60          41 :     set_distribution(model.parameters.template get<TestAndTraceCapacity<FP>>());
      61          41 :     set_distribution(model.parameters.template get<DynamicNPIsImplementationDelay<FP>>());
      62             : 
      63             :     // populations
      64         128 :     for (auto i = AgeGroup(0); i < model.parameters.get_num_groups(); i++) {
      65         957 :         for (auto j = Index<InfectionState>(0); j < Index<InfectionState>(InfectionState::Count); j++) {
      66             : 
      67             :             // don't touch S and D
      68         870 :             if (j == InfectionState::Susceptible || j == InfectionState::Dead) {
      69         174 :                 continue;
      70             :             }
      71             : 
      72             :             // variably sized groups
      73         696 :             set_distribution(model.populations[{i, j}], 0.0);
      74             :         }
      75             :     }
      76             : 
      77             :     // times
      78         128 :     for (auto i = AgeGroup(0); i < model.parameters.get_num_groups(); i++) {
      79             : 
      80          87 :         set_distribution(model.parameters.template get<TimeExposed<FP>>()[i]);
      81          87 :         set_distribution(model.parameters.template get<TimeInfectedNoSymptoms<FP>>()[i]);
      82          87 :         set_distribution(model.parameters.template get<TimeInfectedSymptoms<FP>>()[i]);
      83          87 :         set_distribution(model.parameters.template get<TimeInfectedSevere<FP>>()[i]);
      84          87 :         set_distribution(model.parameters.template get<TimeInfectedCritical<FP>>()[i]);
      85             : 
      86          87 :         set_distribution(model.parameters.template get<TransmissionProbabilityOnContact<FP>>()[i]);
      87          87 :         set_distribution(model.parameters.template get<RelativeTransmissionNoSymptoms<FP>>()[i]);
      88          87 :         set_distribution(model.parameters.template get<RecoveredPerInfectedNoSymptoms<FP>>()[i]);
      89          87 :         set_distribution(model.parameters.template get<RiskOfInfectionFromSymptomatic<FP>>()[i]);
      90          87 :         set_distribution(model.parameters.template get<MaxRiskOfInfectionFromSymptomatic<FP>>()[i]);
      91          87 :         set_distribution(model.parameters.template get<DeathsPerCritical<FP>>()[i]);
      92          87 :         set_distribution(model.parameters.template get<SeverePerInfectedSymptoms<FP>>()[i]);
      93          87 :         set_distribution(model.parameters.template get<CriticalPerSevere<FP>>()[i]);
      94             :     }
      95             : 
      96             :     // dampings
      97          41 :     auto matrices = std::vector<size_t>();
      98          82 :     for (size_t i = 0; i < model.parameters.template get<ContactPatterns<FP>>().get_cont_freq_mat().get_num_matrices();
      99             :          ++i) {
     100          41 :         matrices.push_back(i);
     101             :     }
     102          41 :     auto groups = Eigen::VectorXd::Constant(Eigen::Index(model.parameters.get_num_groups().get()), 1.0);
     103          82 :     model.parameters.template get<ContactPatterns<FP>>().get_dampings().emplace_back(
     104          82 :         mio::UncertainValue(0.5), mio::DampingLevel(0), mio::DampingType(0), mio::SimulationTime(t0 + (tmax - t0) / 2),
     105             :         matrices, groups);
     106          41 :     set_distribution(model.parameters.template get<ContactPatterns<FP>>().get_dampings()[0].get_value(), 0.0);
     107          82 : }
     108             : 
     109             : /**
     110             :  * draws a sample from the specified distributions for all parameters related to the demographics, e.g. population.
     111             :  * @param[inout] model Model including contact patterns for alle age groups
     112             :  */
     113             : template <typename FP = double>
     114         138 : void draw_sample_demographics(Model<FP>& model)
     115             : {
     116         138 :     model.parameters.template get<ICUCapacity<FP>>().draw_sample();
     117         138 :     model.parameters.template get<TestAndTraceCapacity<FP>>().draw_sample();
     118             : 
     119         354 :     for (auto i = AgeGroup(0); i < model.parameters.get_num_groups(); i++) {
     120         216 :         double group_total = model.populations.get_group_total(i);
     121             : 
     122         216 :         model.populations[{i, InfectionState::Exposed}].draw_sample();
     123         216 :         model.populations[{i, InfectionState::InfectedNoSymptoms}].draw_sample();
     124         216 :         model.populations[{i, InfectionState::InfectedSymptoms}].draw_sample();
     125         216 :         model.populations[{i, InfectionState::InfectedSevere}].draw_sample();
     126         216 :         model.populations[{i, InfectionState::InfectedCritical}].draw_sample();
     127         216 :         model.populations[{i, InfectionState::Recovered}].draw_sample();
     128             : 
     129             :         // no sampling for dead and total numbers
     130             :         // [...]
     131             : 
     132         216 :         model.populations.template set_difference_from_group_total<AgeGroup>({i, InfectionState::Susceptible},
     133             :                                                                              group_total);
     134         432 :         model.populations.template set_difference_from_group_total<AgeGroup>({i, InfectionState::Susceptible},
     135         216 :                                                                              model.populations.get_group_total(i));
     136             :     }
     137         138 : }
     138             : 
     139             : /**
     140             :  * draws a sample from the specified distributions for all parameters related to the infection.
     141             :  * @param[inout] model Model including contact patterns for alle age groups
     142             :  */
     143             : template <typename FP = double>
     144         123 : void draw_sample_infection(Model<FP>& model)
     145             : {
     146         123 :     model.parameters.template get<Seasonality<FP>>().draw_sample();
     147             : 
     148             :     //not age dependent
     149         123 :     model.parameters.template get<TimeExposed<FP>>()[AgeGroup(0)].draw_sample();
     150         123 :     model.parameters.template get<TimeInfectedNoSymptoms<FP>>()[AgeGroup(0)].draw_sample();
     151         123 :     model.parameters.template get<TimeInfectedSymptoms<FP>>()[AgeGroup(0)].draw_sample();
     152         123 :     model.parameters.template get<RelativeTransmissionNoSymptoms<FP>>()[AgeGroup(0)].draw_sample();
     153         123 :     model.parameters.template get<RiskOfInfectionFromSymptomatic<FP>>()[AgeGroup(0)].draw_sample();
     154         123 :     model.parameters.template get<MaxRiskOfInfectionFromSymptomatic<FP>>()[AgeGroup(0)].draw_sample();
     155             : 
     156         294 :     for (auto i = AgeGroup(0); i < model.parameters.get_num_groups(); i++) {
     157             :         //not age dependent
     158         171 :         model.parameters.template get<TimeExposed<FP>>()[i] =
     159         342 :             model.parameters.template get<TimeExposed<FP>>()[AgeGroup(0)];
     160         171 :         model.parameters.template get<TimeInfectedNoSymptoms<FP>>()[i] =
     161         342 :             model.parameters.template get<TimeInfectedNoSymptoms<FP>>()[AgeGroup(0)];
     162         171 :         model.parameters.template get<RelativeTransmissionNoSymptoms<FP>>()[i] =
     163         342 :             model.parameters.template get<RelativeTransmissionNoSymptoms<FP>>()[AgeGroup(0)];
     164         171 :         model.parameters.template get<RiskOfInfectionFromSymptomatic<FP>>()[i] =
     165         342 :             model.parameters.template get<RiskOfInfectionFromSymptomatic<FP>>()[AgeGroup(0)];
     166         171 :         model.parameters.template get<MaxRiskOfInfectionFromSymptomatic<FP>>()[i] =
     167         342 :             model.parameters.template get<MaxRiskOfInfectionFromSymptomatic<FP>>()[AgeGroup(0)];
     168             : 
     169             :         //age dependent
     170         171 :         model.parameters.template get<TimeInfectedSymptoms<FP>>()[i].draw_sample();
     171         171 :         model.parameters.template get<TimeInfectedSevere<FP>>()[i].draw_sample();
     172         171 :         model.parameters.template get<TimeInfectedCritical<FP>>()[i].draw_sample();
     173             : 
     174         171 :         model.parameters.template get<TransmissionProbabilityOnContact<FP>>()[i].draw_sample();
     175         171 :         model.parameters.template get<RecoveredPerInfectedNoSymptoms<FP>>()[i].draw_sample();
     176         171 :         model.parameters.template get<DeathsPerCritical<FP>>()[i].draw_sample();
     177         171 :         model.parameters.template get<SeverePerInfectedSymptoms<FP>>()[i].draw_sample();
     178         171 :         model.parameters.template get<CriticalPerSevere<FP>>()[i].draw_sample();
     179             :     }
     180         123 : }
     181             : 
     182             : /** Draws a sample from Model parameter distributions and stores sample values
     183             :  * as SecirParams parameter values (cf. UncertainValue and SecirParams classes).
     184             :  * @param[inout] model Model including contact patterns for alle age groups.
     185             :  */
     186             : template <typename FP = double>
     187          99 : void draw_sample(Model<FP>& model)
     188             : {
     189          99 :     draw_sample_infection(model);
     190          99 :     draw_sample_demographics(model);
     191          99 :     model.parameters.template get<ContactPatterns<FP>>().draw_sample();
     192          99 :     model.apply_constraints();
     193          99 : }
     194             : 
     195             : template <typename FP = double>
     196          24 : Graph<Model<FP>, MobilityParameters<FP>> draw_sample(Graph<Model<FP>, MobilityParameters<FP>>& graph)
     197             : {
     198          24 :     Graph<Model<FP>, MobilityParameters<FP>> sampled_graph;
     199             : 
     200             :     //sample global parameters
     201          24 :     auto& shared_params_model = graph.nodes()[0].property;
     202          24 :     draw_sample_infection(shared_params_model);
     203          24 :     auto& shared_contacts = shared_params_model.parameters.template get<ContactPatterns<FP>>();
     204          24 :     shared_contacts.draw_sample_dampings();
     205          24 :     auto& shared_dynamic_npis = shared_params_model.parameters.template get<DynamicNPIsInfectedSymptoms<FP>>();
     206          24 :     shared_dynamic_npis.draw_sample();
     207          24 :     auto& shared_dynamic_npis_delay = shared_params_model.parameters.template get<DynamicNPIsImplementationDelay<FP>>();
     208          24 :     shared_dynamic_npis_delay.draw_sample();
     209             : 
     210         102 :     for (auto& params_node : graph.nodes()) {
     211          39 :         auto& node_model = params_node.property;
     212             : 
     213             :         //sample local parameters
     214          39 :         draw_sample_demographics(params_node.property);
     215             : 
     216             :         //copy global parameters
     217             :         //save demographic parameters so they aren't overwritten
     218          39 :         auto local_icu_capacity = node_model.parameters.template get<ICUCapacity<FP>>();
     219          39 :         auto local_tnt_capacity = node_model.parameters.template get<TestAndTraceCapacity<FP>>();
     220          39 :         auto local_holidays     = node_model.parameters.template get<ContactPatterns<FP>>().get_school_holidays();
     221          39 :         node_model.parameters   = shared_params_model.parameters;
     222          39 :         node_model.parameters.template get<ICUCapacity<FP>>()                           = local_icu_capacity;
     223          39 :         node_model.parameters.template get<TestAndTraceCapacity<FP>>()                  = local_tnt_capacity;
     224          39 :         node_model.parameters.template get<ContactPatterns<FP>>().get_school_holidays() = local_holidays;
     225             : 
     226          39 :         node_model.parameters.template get<ContactPatterns<FP>>().make_matrix();
     227          39 :         node_model.apply_constraints();
     228             : 
     229          39 :         sampled_graph.add_node(params_node.id, node_model);
     230             :     }
     231             : 
     232          54 :     for (auto& edge : graph.edges()) {
     233          15 :         auto edge_params = edge.property;
     234          45 :         apply_dampings(edge_params.get_coefficients(), shared_contacts.get_dampings(), [&edge_params](auto& v) {
     235          30 :             return make_mobility_damping_vector(edge_params.get_coefficients().get_shape(), v);
     236             :         });
     237          15 :         edge_params.set_dynamic_npis_infected(shared_dynamic_npis);
     238          15 :         sampled_graph.add_edge(edge.start_node_idx, edge.end_node_idx, edge_params);
     239             :     }
     240             : 
     241          48 :     return sampled_graph;
     242          24 : }
     243             : 
     244             : } // namespace osecir
     245             : } // namespace mio
     246             : 
     247             : #endif // ODESECIR_PARAMETER_SPACE_H