2026 Memilio paper

2026_Bicker_et_al_Memilio_paper/Fig5/README.md

@@ -0,0 +1,37 @@
+# Figure 5: Hybrid-OSECIR-dABM comparison and ABM-IDE-LCT-ODE comparison
+
+## Requirements
+
+- **MEmilio**: See `git_tag.cmake` for the required version.
+- **Python**: Required for plotting results and the Python benchmark script.
+
+## Files
+
+- **C++ Simulations**:
+    - `comparison_dabm_ode_hybrid.cpp`: Main simulation file for scenario comparison of dABM, OSECIR and temporal-hybrid (Fig. 5b and c).
+    - `abm_ide_lct_ode.cpp`: Main simulation file for scenario comparison of ABM, ODE-SECIR, LCT-SECIR and IDE-SECIR (Fig. 5d and e).
+- **Python Simulations**:
+    - `get_lognormal_and_erlang_parameters.py`: Returns for given mean and variance the corresponding Lognormal and Erlang parameters used for S2 in Fig. 5e.
+- **Plotting**:
+    - `plot_hybrid_application.py`: Generates Fig. 5b+c from outputs of comparison_dabm_ode_hybrid.cpp.
+    - `plot_model_comparison.py`: Generates Fig. 5e from outputs of abm_ide_lct_ode.cpp.
+    - `plot_pdfs.py`: Plots probability density function of distributions used for Fig.5e.
+
+## How to Run
+
+1.  **Build C++ Benchmarks**:
+    Configure and build the project using CMake, ensuring the `memilio` library is checked out at the commit specified in `git_tag.cmake`. For details, see the Memilio documentation.
+
+2. **Run C++ Simulations**:
+    Execute the compiled C++ simulation binary:
+    ```bash
+    ./bin/comparison_dabm_ode_hybrid
+    ./bin/abm_ide_lct_ode
+    ```
+
+4.  **Generate Plots**:
+    Use the Python scripts for visualization by setting the simulation output directory and executing e.g.:
+    ```bash
+    python3 plot_hybrid_application.py
+    ```
+    Note, that this script uses the general configuration file 'plottings_settins.py'.

2026_Bicker_et_al_Memilio_paper/Fig5/config.h

@@ -0,0 +1,163 @@
+/* 
+* Copyright (C) 2020-2025 MEmilio
+*
+* Authors: Julia Bicker, Anna Wendler
+*
+* Contact: Martin J. Kuehn <Martin.Kuehn@DLR.de>
+*
+* Licensed under the Apache License, Version 2.0 (the "License");
+* you may not use this file except in compliance with the License.
+* You may obtain a copy of the License at
+*
+*     http://www.apache.org/licenses/LICENSE-2.0
+*
+* Unless required by applicable law or agreed to in writing, software
+* distributed under the License is distributed on an "AS IS" BASIS,
+* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+* See the License for the specific language governing permissions and
+* limitations under the License.
+*/
+
+// Define parameters for the simulation.
+#include "abm/location.h"
+#include "abm/model.h"
+#include "abm/simulation.h"
+#include "abm/time.h"
+#include "abm/virus_variant.h"
+#include "memilio/epidemiology/age_group.h"
+#include "memilio/io/history.h"
+#include "memilio/utils/compiler_diagnostics.h"
+#include <algorithm>
+#include <numeric>
+#include <vector>
+namespace params
+{
+const size_t num_age_groups = 6;
+// Define age group sizes and resulting total population of Germany.
+const std::vector<double> age_group_sizes_germany = {3969138.0, 7508662, 18921292, 28666166, 18153339, 5936434};
+const ScalarType total_population_germany =
+    std::accumulate(age_group_sizes_germany.begin(), age_group_sizes_germany.end(), 0.);
+// Define scaling factor so that we can scale population down for faster simulation.
+const size_t scaling_factor       = 1000;
+const ScalarType total_population = total_population_germany / scaling_factor;
+
+// Define transition probabilities per age group.
+const ScalarType infectedSymptomsPerInfectedNoSymptoms[] = {0.75, 0.75, 0.8, 0.8, 0.8, 0.8};
+const ScalarType severePerInfectedSymptoms[]             = {0.0075, 0.0075, 0.019, 0.0615, 0.165, 0.225};
+const ScalarType criticalPerSevere[]                     = {0.075, 0.075, 0.075, 0.15, 0.3, 0.4};
+const ScalarType deathsPerCritical[]                     = {0.05, 0.05, 0.14, 0.14, 0.4, 0.6};
+
+// Define lognormal parameters. For each transition, we need shape and scale.
+// These are given in this order below. The transition distributions are the same for all age groups.
+const ScalarType lognorm_EtINS[]     = {1.45139714, 0.32459285}; //{0.32459285, 4.26907484};
+const ScalarType lognorm_INStISy[]   = {-0.1609284, 0.7158751}; //{0.7158751, 0.85135303};
+const ScalarType lognorm_INStR[]     = {2.04912923, 0.24622068}; //{0.24622068, 7.76114};
+const ScalarType lognorm_ISytISev[]  = {1.66755725, 0.66258947}; //{0.66258947, 5.29920733};
+const ScalarType lognorm_ISytR[]     = {2.04912923, 0.24622068}; //{0.24622068, 7.76114};
+const ScalarType lognorm_ISevtICri[] = {-0.10536052, 1.01076765}; //{1.01076765, 0.9};
+const ScalarType lognorm_ISevtR[]    = {2.8387344, 0.33816427}; //{0.33816427, 17.09411753};
+const ScalarType lognorm_ICritD[]    = {2.27856645, 0.42819924}; //{0.42819924, 9.76267505};
+const ScalarType lognorm_ICritR[]    = {2.8387344, 0.33816427}; //{0.33816427, 17.09411753};
+
+// Define mean stay times per age group. Note that these are different per age group as the transition probabilities
+// differ between age groups.
+const ScalarType timeExposed[]            = {4.5, 4.5, 4.5, 4.5, 4.5, 4.5};
+const ScalarType timeInfectedNoSymptoms[] = {2.825, 2.825, 2.48, 2.48, 2.48, 2.48};
+const ScalarType timeInfectedSymptoms[]   = {7.9895, 7.9895, 7.9734, 7.9139, 7.9139, 7.685};
+const ScalarType timeInfectedSevere[]     = {16.855, 16.855, 16.855, 15.61, 13.12, 11.46};
+const ScalarType timeInfectedCritical[]   = {17.73, 17.73, 17.064, 17.064, 15.14, 13.66};
+
+// Define number of subcompartments for each compartment per age group. Note that these are different per age group as
+// the transition probabilities differ between age groups.
+// These values are used as a template argument and thus have to be constexpr.
+constexpr size_t n_subcomps_E[]    = {9, 9, 9, 9, 9, 9};
+constexpr size_t n_subcomps_INS[]  = {5, 5, 4, 4, 4, 4};
+constexpr size_t n_subcomps_ISy[]  = {16, 16, 16, 15, 15, 13};
+constexpr size_t n_subcomps_ISev[] = {8, 8, 8, 7, 6, 5};
+constexpr size_t n_subcomps_ICri[] = {8, 8, 8, 8, 7, 6};
+
+// Define epidemiological parameters.
+const ScalarType relativeTransmissionNoSymptoms = 1;
+const ScalarType riskOfInfectionFromSymptomatic = 1.;
+const ScalarType seasonality                    = 0.; // this implies that we don't have any seasonal effects in EBMs
+const ScalarType scale_confirmed_cases          = 1.;
+
+// Define tolerance for computing the support max for IDE model and initialization of LCT model.
+const ScalarType tol_supp_max = 1e-8;
+
+// Define simulation parameters.
+const ScalarType t0        = 0.;
+const ScalarType init_tmax = 20.;
+const ScalarType tmax      = 120.;
+const ScalarType dt        = 1. / 24.; // corresponds to hours that are used as time step in ABM simulation
+
+} // namespace params
+
+namespace ABMparams
+{
+
+// Distribution of household sizes from 1 to 5
+const ScalarType household_size_distribution[] = {0.415, 0.342, 0.118, 0.091, 0.034};
+// Location parameters. All location sizes are sampled from a normal distribution which is capped at the given minimum value.
+// Mean size and stddev of workplaces
+const size_t workplace_size[] = {15, 25};
+// Minimum workplace size
+const size_t min_workplace_size = 5;
+// Mean size and stddev of schools
+const size_t school_size[] = {42, 3};
+// Minimum school size
+const size_t min_school_size = 15;
+// Mean size and stddev of events
+const size_t event_size[] = {42, 941};
+// Minimum event size
+const size_t min_event_size = 10;
+// Mean size and stddev of shops
+const size_t shop_size[] = {90, 1000};
+// Minimum shop size
+const size_t min_shop_size = 30;
+
+} // namespace ABMparams
+
+namespace ABMLoggers
+{
+struct LogTimePoint : mio::LogAlways {
+    using Type = mio::abm::TimePoint;
+    static Type log(const mio::abm::Simulation<mio::abm::Model>& sim)
+    {
+        return sim.get_time();
+    }
+};
+
+struct LogExposureRate : mio::LogAlways {
+    using Type = std::vector<double>;
+    static Type log(const mio::abm::Simulation<mio::abm::Model>& sim)
+    {
+        std::vector<double> cum_rate(params::num_age_groups);
+        auto& rates = sim.get_model().get_contact_exposure_rates();
+        if (rates.size() > 0) {
+            size_t num_locs = 0;
+            for (auto& loc : sim.get_model().get_locations()) {
+                bool counted_loc = false;
+                for (size_t age = 0; age < params::num_age_groups; ++age) {
+                    auto& rate = rates[loc.get_id().get()]
+                                      [{mio::abm::CellIndex(0), mio::abm::VirusVariant::Wildtype, mio::AgeGroup(age)}];
+                    cum_rate[age] += rate;
+                    if (rate > 0 && !counted_loc) {
+                        num_locs += 1;
+                        counted_loc = true;
+                    }
+                }
+            }
+
+            if (num_locs < 1) {
+                num_locs = 1;
+            }
+            // Divide rate by number of locations
+            std::transform(cum_rate.begin(), cum_rate.end(), cum_rate.begin(), [num_locs](double x) {
+                return x / num_locs;
+            });
+        }
+        return cum_rate;
+    }
+};
+} // namespace ABMLoggers

2026_Bicker_et_al_Memilio_paper/Fig5/get_erlang_parameters.py

@@ -0,0 +1,130 @@
+#############################################################################
+# Copyright (C) 2020-2025 MEmilio
+#
+# Authors: Anna Wendler
+#
+# Contact: Martin J. Kuehn <Martin.Kuehn@DLR.de>
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+#############################################################################
+import numpy as np
+
+
+def compute_optimal_num_subcompartments(mean, std):
+
+    optimal_num_subcompartments = mean**2/std**2
+
+    return optimal_num_subcompartments
+
+
+def get_weighted_value(prob_1, value_1, value_2):
+
+    weighted_value = prob_1*value_1 + (1-prob_1)*value_2
+
+    return weighted_value
+
+
+def main():
+
+    # Mean and standard deviation of lognormal distributions per transiton from Covasim paper.
+    mean_and_std_per_transition = {"ExposedToInfectedNoSymptoms": [4.5, 1.5],
+                                   "InfectedNoSymptomsToInfectedSymptoms": [1.1, 0.9],
+                                   "InfectedNoSymptomsToRecovered": [8.0, 2.0],
+                                   "InfectedSymptomsToInfectedSevere": [6.6, 4.9],
+                                   "InfectedSymptomsToRecovered": [8.0, 2.0],
+                                   "InfectedSevereToInfectedCritical": [1.5, 2.0],
+                                   "InfectedSevereToRecovered": [18.1, 6.3],
+                                   "InfectedCriticalToDead": [10.7, 4.8],
+                                   "InfectedCriticalToRecovered": [18.1, 6.3]}
+
+    # Transition probabilities from Assessment paper, see LCT paper.
+    # Store only the probabilities that are not considering the transition to the Recovered compartment as this is
+    # determined by the given probabilities.
+    transition_probabilities_non_recovered = {"ExposedToInfectedNoSymptoms": [1., 1., 1., 1., 1., 1.],
+                                              "InfectedNoSymptomsToInfectedSymptoms": [0.75, 0.75, 0.8, 0.8, 0.8, 0.8],
+                                              "InfectedSymptomsToInfectedSevere": [0.0075, 0.0075, 0.019, 0.0615, 0.0615, 0.225],
+                                              "InfectedSevereToInfectedCritical": [0.075, 0.075, 0.075, 0.15, 0.3, 0.4],
+                                              "InfectedCriticalToDead": [0.05, 0.05, 0.14, 0.14, 0.4, 0.6]}
+
+    compartments = ["Exposed", "InfectedNoSymptoms",
+                    "InfectedSymptoms", "InfectedSevere", "InfectedCritical"]
+
+    # Get the number of subcompartment for the LCT model.
+    # To get the numvber of subcompartments, we first calculate the optimal number of subcompartments that is not
+    # necessarily an integer. In the LCT model we cannot consider different distributions for e.g. the transition
+    # for InfectedNoSymptomsToInfectedSymptoms and InfectedNoSymptomsToRecovered but only have one transition for
+    # leaving the InfectedNoSymptoms compartments. This is why we weigh the optimal number of subcompartments
+    # of both transitions with the corresponding transition probability.
+    #
+    # For the same reason, we weigh the mean of the corresponding transitions from Covasim by the corresponding
+    # transition probabilities.
+    # Note that the mean is the same for the LCT and the ODE model, only the number of subcompartments differs and is
+    # always one for the ODE model.
+
+    print("Number of subcompartments for the LCT and mean for the LCT and ODE models:")
+    print()
+
+    # Compute the optimal number of subcompartments for each transition.
+    optimal_num_subcompartments = []
+    for i in mean_and_std_per_transition.keys():
+        optimal_num_subcompartment = compute_optimal_num_subcompartments(
+            mean_and_std_per_transition[i][0], mean_and_std_per_transition[i][1])
+        optimal_num_subcompartments.append(optimal_num_subcompartment)
+
+    for i, transition in enumerate(transition_probabilities_non_recovered.keys()):
+
+        # For the transition ExposedToInfectedNoSymptoms we do not need to weigh anything as there is only one
+        # possible transition from Exposed.
+        if transition == "ExposedToInfectedNoSymptoms":
+            weighted_num_subcompartments_per_age = []
+            weighted_mean_per_age = []
+
+            for age in range(len(transition_probabilities_non_recovered[transition])):
+                weighted_num_subcompartments_per_age.append(round(
+                    optimal_num_subcompartments[i]))
+                weighted_mean_per_age.append(
+                    mean_and_std_per_transition[transition][0])
+
+        else:
+            weight = transition_probabilities_non_recovered[transition]
+            # Get the indices of the two transitions that we need to weigh.
+            first_index = 2*i-1
+            second_index = 2*i
+            # print(
+            #     f"Weighing transitions {keys_list[first_index]} and {keys_list[second_index]} with weight {weight}.")
+
+            # Get the corresponding keys of the transitions.
+            keys_list = list(mean_and_std_per_transition.keys())
+            first_key = keys_list[first_index]
+            second_key = keys_list[second_index]
+
+            # Weigh the optimal number of subcompartments of both transitions with the corresponding transition probability.
+            weighted_num_subcompartments_per_age = []
+            weighted_mean_per_age = []
+            for age in range(len(transition_probabilities_non_recovered[transition])):
+                # Round to integer.
+                weighted_num_subcompartments_per_age.append(round(get_weighted_value(
+                    weight[age], optimal_num_subcompartments[first_index], optimal_num_subcompartments[second_index])))
+                # Round to 8 digits.
+                weighted_mean_per_age.append(round(get_weighted_value(
+                    weight[age], mean_and_std_per_transition[first_key][0], mean_and_std_per_transition[second_key][0]), 8))
+
+        print(compartments[i])
+        print(
+            f"Num subcompartments: {weighted_num_subcompartments_per_age}")
+        print(f"Mean: {weighted_mean_per_age}")
+        print()
+
+
+if __name__ == '__main__':
+    main()

2026_Bicker_et_al_Memilio_paper/Fig5/git_tag.cmake

@@ -0,0 +1 @@
+set(GIT_TAG_MEMILIO "91b2665bc02b01d1332961f9e05a0e5c00fb12b0")