Line data    Source code

       1             : /* 
       2             : * Copyright (C) 2020-2025 MEmilio
       3             : *
       4             : * Authors: Henrik Zunker, Wadim Koslow, Daniel Abele, Martin J. Kühn
       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 MIO_ODE_SECIRTS_ANALYZE_RESULT_H
      21             : #define MIO_ODE_SECIRTS_ANALYZE_RESULT_H
      22             : 
      23             : #include "ode_secirts/model.h"
      24             : #include "memilio/data/analyze_result.h"
      25             : #include "ode_secirts/parameters.h"
      26             : 
      27             : #include <functional>
      28             : #include <vector>
      29             : 
      30             : namespace mio
      31             : {
      32             : namespace osecirts
      33             : {
      34             : /**
      35             :     * @brief computes the p percentile of the parameters for each node.
      36             :     * @param ensemble_result graph of multiple simulation runs
      37             :     * @param p percentile value in open interval (0, 1)
      38             :     * @return p percentile of the parameters over all runs
      39             :     */
      40             : template <class Model>
      41           9 : std::vector<Model> ensemble_params_percentile(const std::vector<std::vector<Model>>& ensemble_params, double p)
      42             : {
      43           9 :     assert(p > 0.0 && p < 1.0 && "Invalid percentile value.");
      44             : 
      45           9 :     auto num_runs   = ensemble_params.size();
      46           9 :     auto num_nodes  = ensemble_params[0].size();
      47           9 :     auto num_groups = (size_t)ensemble_params[0][0].parameters.get_num_groups();
      48           9 :     auto num_days   = ensemble_params[0][0]
      49           9 :                         .parameters.template get<DailyPartialVaccinations<double>>()
      50           9 :                         .template size<mio::SimulationDay>();
      51             : 
      52          18 :     std::vector<double> single_element_ensemble(num_runs);
      53             : 
      54             :     // lambda function that calculates the percentile of a single parameter
      55          27 :     std::vector<Model> percentile(num_nodes, Model((int)num_groups));
      56     5082921 :     auto param_percentil = [&ensemble_params, p, num_runs, &percentile](auto n, auto get_param) mutable {
      57      412128 :         std::vector<double> single_element(num_runs);
      58     1648512 :         for (size_t run = 0; run < num_runs; run++) {
      59     2747520 :             auto const& params  = ensemble_params[run][n];
      60     1373760 :             single_element[run] = get_param(params);
      61             :         }
      62      137376 :         std::sort(single_element.begin(), single_element.end());
      63      274752 :         auto& new_params = get_param(percentile[n]);
      64      274752 :         new_params       = single_element[static_cast<size_t>(num_runs * p)];
      65      137376 :     };
      66             : 
      67          18 :     for (size_t node = 0; node < num_nodes; node++) {
      68           9 :         percentile[node].parameters.template get<DailyPartialVaccinations<double>>().resize(num_days);
      69           9 :         percentile[node].parameters.template get<DailyFullVaccinations<double>>().resize(num_days);
      70           9 :         percentile[node].parameters.template get<DailyBoosterVaccinations<double>>().resize(num_days);
      71             : 
      72          54 :         for (auto i = AgeGroup(0); i < AgeGroup(num_groups); i++) {
      73             :             //Population
      74        1350 :             for (auto compart = Index<InfectionState>(0); compart < InfectionState::Count; ++compart) {
      75        1305 :                 param_percentil(
      76       14355 :                     node, [ compart, i ](auto&& model) -> auto& {
      77       14355 :                         return model.populations[{i, compart}];
      78             :                     });
      79             :             }
      80             :             // times
      81          45 :             param_percentil(
      82         495 :                 node, [i](auto&& model) -> auto& { return model.parameters.template get<TimeExposed<double>>()[i]; });
      83          45 :             param_percentil(
      84         495 :                 node, [i](auto&& model) -> auto& {
      85         495 :                     return model.parameters.template get<TimeInfectedNoSymptoms<double>>()[i];
      86             :                 });
      87          45 :             param_percentil(
      88         495 :                 node, [i](auto&& model) -> auto& {
      89         495 :                     return model.parameters.template get<TimeInfectedSymptoms<double>>()[i];
      90             :                 });
      91          45 :             param_percentil(
      92         495 :                 node, [i](auto&& model) -> auto& {
      93         495 :                     return model.parameters.template get<TimeInfectedSevere<double>>()[i];
      94             :                 });
      95          45 :             param_percentil(
      96         495 :                 node, [i](auto&& model) -> auto& {
      97         495 :                     return model.parameters.template get<TimeInfectedCritical<double>>()[i];
      98             :                 });
      99             :             //probs
     100          45 :             param_percentil(
     101         495 :                 node, [i](auto&& model) -> auto& {
     102         495 :                     return model.parameters.template get<TransmissionProbabilityOnContact<double>>()[i];
     103             :                 });
     104          45 :             param_percentil(
     105         495 :                 node, [i](auto&& model) -> auto& {
     106         495 :                     return model.parameters.template get<RelativeTransmissionNoSymptoms<double>>()[i];
     107             :                 });
     108          45 :             param_percentil(
     109         495 :                 node, [i](auto&& model) -> auto& {
     110         495 :                     return model.parameters.template get<RiskOfInfectionFromSymptomatic<double>>()[i];
     111             :                 });
     112          45 :             param_percentil(
     113         495 :                 node, [i](auto&& model) -> auto& {
     114         495 :                     return model.parameters.template get<MaxRiskOfInfectionFromSymptomatic<double>>()[i];
     115             :                 });
     116          45 :             param_percentil(
     117         495 :                 node, [i](auto&& model) -> auto& {
     118         495 :                     return model.parameters.template get<RecoveredPerInfectedNoSymptoms<double>>()[i];
     119             :                 });
     120          45 :             param_percentil(
     121         495 :                 node, [i](auto&& model) -> auto& {
     122         495 :                     return model.parameters.template get<SeverePerInfectedSymptoms<double>>()[i];
     123             :                 });
     124          45 :             param_percentil(
     125         495 :                 node, [i](auto&& model) -> auto& {
     126         495 :                     return model.parameters.template get<CriticalPerSevere<double>>()[i];
     127             :                 });
     128          45 :             param_percentil(
     129         495 :                 node, [i](auto&& model) -> auto& {
     130         495 :                     return model.parameters.template get<DeathsPerCritical<double>>()[i];
     131             :                 });
     132             :             //vaccinations
     133          45 :             param_percentil(
     134         495 :                 node, [i](auto&& model) -> auto& {
     135         495 :                     return model.parameters.template get<ReducExposedPartialImmunity<double>>()[i];
     136             :                 });
     137          45 :             param_percentil(
     138         495 :                 node, [i](auto&& model) -> auto& {
     139         495 :                     return model.parameters.template get<ReducExposedImprovedImmunity<double>>()[i];
     140             :                 });
     141          45 :             param_percentil(
     142         495 :                 node, [i](auto&& model) -> auto& {
     143         495 :                     return model.parameters.template get<ReducInfectedSymptomsPartialImmunity<double>>()[i];
     144             :                 });
     145          45 :             param_percentil(
     146         495 :                 node, [i](auto&& model) -> auto& {
     147         495 :                     return model.parameters.template get<ReducInfectedSymptomsImprovedImmunity<double>>()[i];
     148             :                 });
     149          45 :             param_percentil(
     150         495 :                 node, [i](auto&& model) -> auto& {
     151         495 :                     return model.parameters.template get<ReducInfectedSevereCriticalDeadPartialImmunity<double>>()[i];
     152             :                 });
     153          45 :             param_percentil(
     154         495 :                 node, [i](auto&& model) -> auto& {
     155         495 :                     return model.parameters.template get<ReducInfectedSevereCriticalDeadImprovedImmunity<double>>()[i];
     156             :                 });
     157          45 :             param_percentil(
     158         495 :                 node, [i](auto&& model) -> auto& {
     159         495 :                     return model.parameters.template get<ReducTimeInfectedMild<double>>()[i];
     160             :                 });
     161          45 :             param_percentil(
     162         495 :                 node, [i](auto&& model) -> auto& {
     163         495 :                     return model.parameters.template get<DaysUntilEffectivePartialVaccination<double>>()[i];
     164             :                 });
     165          45 :             param_percentil(
     166         495 :                 node, [i](auto&& model) -> auto& {
     167         495 :                     return model.parameters.template get<DaysUntilEffectiveImprovedVaccination<double>>()[i];
     168             :                 });
     169             : 
     170       45045 :             for (auto day = SimulationDay(0); day < num_days; ++day) {
     171       45000 :                 param_percentil(
     172      495000 :                     node, [ i, day ](auto&& model) -> auto& {
     173      495000 :                         return model.parameters.template get<DailyPartialVaccinations<double>>()[{i, day}];
     174             :                     });
     175       45000 :                 param_percentil(
     176      495000 :                     node, [ i, day ](auto&& model) -> auto& {
     177      495000 :                         return model.parameters.template get<DailyFullVaccinations<double>>()[{i, day}];
     178             :                     });
     179       45000 :                 param_percentil(
     180      495000 :                     node, [ i, day ](auto&& model) -> auto& {
     181      495000 :                         return model.parameters.template get<DailyBoosterVaccinations<double>>()[{i, day}];
     182             :                     });
     183             :             }
     184             :             //virus variants
     185          45 :             param_percentil(
     186         495 :                 node, [i](auto&& model) -> auto& {
     187         495 :                     return model.parameters.template get<InfectiousnessNewVariant<double>>()[i];
     188             :                 });
     189             :         }
     190             :         // group independent params
     191           9 :         param_percentil(
     192          99 :             node, [](auto&& model) -> auto& { return model.parameters.template get<Seasonality<double>>(); });
     193           9 :         param_percentil(
     194          99 :             node, [](auto&& model) -> auto& { return model.parameters.template get<TestAndTraceCapacity<double>>(); });
     195           9 :         param_percentil(
     196          99 :             node, [](auto&& model) -> auto& { return model.parameters.template get<ICUCapacity<double>>(); });
     197           9 :         param_percentil(
     198          99 :             node, [](auto&& model) -> auto& {
     199          99 :                 return model.parameters.template get<DynamicNPIsImplementationDelay<double>>();
     200             :             });
     201             : 
     202          99 :         for (size_t run = 0; run < num_runs; run++) {
     203             : 
     204          90 :             auto const& params = ensemble_params[run][node];
     205          90 :             single_element_ensemble[run] =
     206          90 :                 params.parameters.template get<ICUCapacity<double>>() * params.populations.get_total();
     207             :         }
     208           9 :         std::sort(single_element_ensemble.begin(), single_element_ensemble.end());
     209          18 :         percentile[node].parameters.template set<ICUCapacity<double>>(
     210           9 :             single_element_ensemble[static_cast<size_t>(num_runs * p)]);
     211             :     }
     212          18 :     return percentile;
     213           9 : }
     214             : 
     215             : } // namespace osecirts
     216             : } // namespace mio
     217             : 
     218             : #endif //MIO_ODE_SECIRTS_ANALYZE_RESULT_H