Abm paper test bs

cpp/CMakeLists.txt

@@ -6,7 +6,7 @@ option(MEMILIO_BUILD_TESTS "Build memilio unit tests." ON)
 option(MEMILIO_BUILD_EXAMPLES "Build memilio examples." ON)
 option(MEMILIO_BUILD_MODELS "Build memilio models." ON)
 option(MEMILIO_BUILD_SIMULATIONS "Build memilio simulations that were used for scientific articles." ON)
-option(MEMILIO_BUILD_BENCHMARKS "Build memilio benchmarks with google benchmark." OFF)
+option(MEMILIO_BUILD_BENCHMARKS "Build memilio benchmarks with google benchmark." ON)
 option(MEMILIO_USE_BUNDLED_SPDLOG "Use spdlog bundled with epi" ON)
 option(MEMILIO_USE_BUNDLED_EIGEN "Use eigen bundled with epi" ON)
 option(MEMILIO_USE_BUNDLED_BOOST "Use boost bundled with epi (only for epi-io)" ON)
@@ -15,8 +15,8 @@ option(MEMILIO_SANITIZE_ADDRESS "Enable address sanitizer." OFF)
 option(MEMILIO_SANITIZE_UNDEFINED "Enable undefined behavior sanitizer." OFF)
 option(MEMILIO_BUILD_SHARED_LIBS "Build memilio as a shared library." ON)
 option(MEMILIO_BUILD_STATIC_LIBS "Build memilio as a static library." ON)
-option(MEMILIO_ENABLE_MPI "Build memilio with MPI." OFF)
-option(MEMILIO_ENABLE_OPENMP "Enable Multithreading with OpenMP." OFF)
+option(MEMILIO_ENABLE_MPI "Build memilio with MPI." ON)
+option(MEMILIO_ENABLE_OPENMP "Enable Multithreading with OpenMP." ON)
 
 mark_as_advanced(MEMILIO_USE_BUNDLED_SPDLOG MEMILIO_SANITIZE_ADDRESS MEMILIO_SANITIZE_UNDEFINED)
 

cpp/benchmarks/abm.cpp

@@ -1,4 +1,5 @@
 #include "abm/simulation.h"
+#include "memilio/math/interpolation.h"
 #include "memilio/utils/stl_util.h"
 #include "benchmark/benchmark.h"
 
@@ -35,7 +36,7 @@ mio::abm::Simulation make_simulation(size_t num_persons, std::initializer_list<u
          {mio::abm::LocationType::School, mio::abm::LocationType::Work, mio::abm::LocationType::SocialEvent,
           mio::abm::LocationType::BasicsShop, mio::abm::LocationType::Hospital, mio::abm::LocationType::ICU}) {
 
-        const auto num_locs = std::max(size_t(1), num_persons / 2'000);
+        const auto num_locs = std::max(size_t(1), num_persons / 100);
         std::vector<mio::abm::LocationId> locs(num_locs);
         std::generate(locs.begin(), locs.end(), [&] {
             return world.add_location(loc_type);
@@ -108,6 +109,13 @@ mio::abm::Simulation make_simulation(size_t num_persons, std::initializer_list<u
         mio::abm::TestingScheme(random_criteria(), mio::abm::days(3), mio::abm::TimePoint(0),
                                 mio::abm::TimePoint(0) + mio::abm::days(10), {}, 0.5));
 
+    // Necessary to set a parameter for the world
+    world.parameters.get<mio::abm::HighViralLoadProtectionFactor>() = [](ScalarType days) -> ScalarType {
+        return mio::linear_interpolation_of_data_set<ScalarType, ScalarType>(
+            {{0, 0.863}, {1, 0.969}, {7, 0.029}, {10, 0.002}, {14, 0.0014}, {21, 0}}, days);
+    };
+    world.parameters.get<mio::abm::InfectionRateFromViralShed>() = 100;
+
     return mio::abm::Simulation(mio::abm::TimePoint(0), std::move(world));
 }
 
@@ -118,6 +126,15 @@ mio::abm::Simulation make_simulation(size_t num_persons, std::initializer_list<u
  */
 void abm_benchmark(benchmark::State& state, size_t num_persons, std::initializer_list<uint32_t> seeds)
 {
+    int tid = -1;
+    std::cout << "Number of persons = " << num_persons << "\n";
+#pragma omp parallel private(tid) // Start of parallel region: forks threads
+    {
+        tid = omp_get_thread_num(); // default is number of CPUs on machine
+        if (tid == 0) {
+            printf("Number of threads = %d\n", omp_get_num_threads());
+        }
+    } // ** end of the the parallel: joins threads
     mio::set_log_level(mio::LogLevel::warn);
 
     for (auto&& _ : state) {
@@ -126,7 +143,7 @@ void abm_benchmark(benchmark::State& state, size_t num_persons, std::initializer
         state.ResumeTiming();
 
         //simulated time should be long enough to have full infection runs and migration to every location
-        auto final_time = sim.get_time() + mio::abm::days(10);
+        auto final_time = sim.get_time() + mio::abm::days(5);
         sim.advance(final_time);
 
         //debug output can be enabled to check for unexpected results (e.g. infections dieing out)
@@ -152,8 +169,16 @@ void abm_benchmark(benchmark::State& state, size_t num_persons, std::initializer
 //have to be adjusted to get the benchmark back to normal.
 //For small sizes (e.g. 10k) extreme cases are too likely, i.e. infections die out
 //or overwhelm everything, so we don't benchmark these. Results should be mostly transferrable.
-BENCHMARK_CAPTURE(abm_benchmark, abm_benchmark_50k, 50000, {14159265u, 35897932u})->Unit(benchmark::kMillisecond);
-BENCHMARK_CAPTURE(abm_benchmark, abm_benchmark_100k, 100000, {38462643u, 38327950u})->Unit(benchmark::kMillisecond);
-BENCHMARK_CAPTURE(abm_benchmark, abm_benchmark_200k, 200000, {28841971u, 69399375u})->Unit(benchmark::kMillisecond);
-
+BENCHMARK_CAPTURE(abm_benchmark, abm_benchmark_25k_pt, 25000, {14159265u, 35897932u})->Unit(benchmark::kMillisecond);
+BENCHMARK_CAPTURE(abm_benchmark, abm_benchmark_50k_pt, 50000, {14159265u, 35897932u})->Unit(benchmark::kMillisecond);
+BENCHMARK_CAPTURE(abm_benchmark, abm_benchmark_100k_pt, 100000, {14159265u, 35897932u})->Unit(benchmark::kMillisecond);
+BENCHMARK_CAPTURE(abm_benchmark, abm_benchmark_200k_pt, 200000, {14159265u, 35897932u})->Unit(benchmark::kMillisecond);
+BENCHMARK_CAPTURE(abm_benchmark, abm_benchmark_400k_pt, 400000, {14159265u, 35897932u})->Unit(benchmark::kMillisecond);
+BENCHMARK_CAPTURE(abm_benchmark, abm_benchmark_800k_pt, 800000, {14159265u, 35897932u})->Unit(benchmark::kMillisecond);
+BENCHMARK_CAPTURE(abm_benchmark, abm_benchmark_1600k_pt, 1600000, {14159265u, 35897932u})
+    ->Unit(benchmark::kMillisecond);
+BENCHMARK_CAPTURE(abm_benchmark, abm_benchmark_3200k_pt, 3200000, {14159265u, 35897932u})
+    ->Unit(benchmark::kMillisecond);
+BENCHMARK_CAPTURE(abm_benchmark, abm_benchmark_3200k_pt, 6400000, {14159265u, 35897932u})
+    ->Unit(benchmark::kMillisecond);
 BENCHMARK_MAIN();

cpp/examples/abm_history_object.cpp

@@ -125,10 +125,11 @@ int main()
     auto probability           = 0.5;
     auto start_date            = mio::abm::TimePoint(0);
     auto end_date              = mio::abm::TimePoint(0) + mio::abm::days(30);
-    auto test_type             = mio::abm::AntigenTest();
+    auto test_type             = mio::abm::TestType::Antigen;
+    auto test_parameters       = world.parameters.get<mio::abm::TestData>()[test_type];
     auto testing_criteria_work = mio::abm::TestingCriteria();
-    auto testing_scheme_work =
-        mio::abm::TestingScheme(testing_criteria_work, testing_min_time, start_date, end_date, test_type, probability);
+    auto testing_scheme_work   = mio::abm::TestingScheme(testing_criteria_work, testing_min_time, start_date, end_date,
+                                                         test_parameters, probability);
     world.get_testing_strategy().add_testing_scheme(mio::abm::LocationType::Work, testing_scheme_work);
 
     // Assign infection state to each person.

cpp/examples/abm_minimal.cpp

@@ -41,10 +41,10 @@ int main()
     world.parameters.get<mio::abm::IncubationPeriod>() = 4.;
 
     // Set the age group the can go to school is AgeGroup(1) (i.e. 5-14)
-    world.parameters.get<mio::abm::AgeGroupGotoSchool>() = false;
+    world.parameters.get<mio::abm::AgeGroupGotoSchool>()                    = false;
     world.parameters.get<mio::abm::AgeGroupGotoSchool>()[age_group_5_to_14] = true;
     // Set the age group the can go to work is AgeGroup(2) and AgeGroup(3) (i.e. 15-34 and 35-59)
-    world.parameters.get<mio::abm::AgeGroupGotoWork>() = false;
+    world.parameters.get<mio::abm::AgeGroupGotoWork>()                     = false;
     world.parameters.get<mio::abm::AgeGroupGotoWork>()[age_group_15_to_34] = true;
     world.parameters.get<mio::abm::AgeGroupGotoWork>()[age_group_35_to_59] = true;
 
@@ -110,10 +110,11 @@ int main()
     auto probability           = 0.5;
     auto start_date            = mio::abm::TimePoint(0);
     auto end_date              = mio::abm::TimePoint(0) + mio::abm::days(10);
-    auto test_type             = mio::abm::AntigenTest();
+    auto test_type             = mio::abm::TestType::Antigen;
+    auto test_parameters       = world.parameters.get<mio::abm::TestData>()[test_type];
     auto testing_criteria_work = mio::abm::TestingCriteria();
-    auto testing_scheme_work =
-        mio::abm::TestingScheme(testing_criteria_work, testing_min_time, start_date, end_date, test_type, probability);
+    auto testing_scheme_work   = mio::abm::TestingScheme(testing_criteria_work, testing_min_time, start_date, end_date,
+                                                         test_parameters, probability);
     world.get_testing_strategy().add_testing_scheme(mio::abm::LocationType::Work, testing_scheme_work);
 
     // Assign infection state to each person.

cpp/memilio/epidemiology/damping.h

@@ -519,7 +519,7 @@ class Dampings
     static void inclusive_exclusive_sum_rec(Iter b, Iter e, Matrix& sum)
     {
         if (b != e) {
-            sum = sum + std::get<Matrix>(*b) - (sum.array() * std::get<Matrix>(*b).array()).matrix();
+            sum = (sum + std::get<Matrix>(*b) - (sum.array() * std::get<Matrix>(*b).array()).matrix()).eval();
             inclusive_exclusive_sum_rec(++b, e, sum);
         }
     }

cpp/memilio/io/epi_data.h

@@ -456,7 +456,7 @@ class VaccinationDataEntry
 public:
     static std::vector<const char*> age_group_names;
 
-    double num_vaccinations_completed;
+    double num_vaccinations_completed, num_vaccinations_partially;
     Date date;
     AgeGroup age_group;
     boost::optional<regions::StateId> state_id;
@@ -468,14 +468,16 @@ class VaccinationDataEntry
     {
         auto obj                        = io.expect_object("VaccinationDataEntry");
         auto num_vaccinations_completed = obj.expect_element("Vacc_completed", Tag<double>{});
+        auto num_vaccinations_partially = obj.expect_element("Vacc_partially", Tag<double>{});
         auto date                       = obj.expect_element("Date", Tag<StringDate>{});
         auto age_group_str              = obj.expect_element("Age_RKI", Tag<std::string>{});
         auto state_id                   = obj.expect_optional("ID_State", Tag<regions::StateId>{});
         auto county_id                  = obj.expect_optional("ID_County", Tag<regions::CountyId>{});
         auto district_id                = obj.expect_optional("ID_District", Tag<regions::DistrictId>{});
         return mio::apply(
             io,
-            [](auto nf, auto d, auto&& a_str, auto sid, auto cid, auto did) -> IOResult<VaccinationDataEntry> {
+            [](auto nf, auto nf2, auto d, auto&& a_str, auto sid, auto cid,
+               auto did) -> IOResult<VaccinationDataEntry> {
                 auto it = std::find(age_group_names.begin(), age_group_names.end(), a_str);
                 auto a  = AgeGroup(0);
                 if (it != age_group_names.end()) {
@@ -484,9 +486,10 @@ class VaccinationDataEntry
                 else {
                     return failure(StatusCode::InvalidValue, "Invalid vaccination data age group.");
                 }
-                return success(VaccinationDataEntry{nf, d, a, sid, cid, did});
+                return success(VaccinationDataEntry{nf, nf2, d, a, sid, cid, did});
             },
-            num_vaccinations_completed, date, age_group_str, state_id, county_id, district_id);
+            num_vaccinations_completed, num_vaccinations_partially, date, age_group_str, state_id, county_id,
+            district_id);
     }
 };
 

cpp/memilio/io/result_io.h

@@ -135,11 +135,15 @@ IOResult<void> save_result_with_params(const std::vector<TimeSeries<double>>& re
 template <class Model>
 IOResult<void> save_results(const std::vector<std::vector<TimeSeries<double>>>& ensemble_results,
                             const std::vector<std::vector<Model>>& ensemble_params, const std::vector<int>& county_ids,
-                            const fs::path& result_dir, bool save_single_runs = true, bool save_percentiles = true)
+                            const fs::path& result_dir, bool save_single_runs = true, bool save_percentiles = true,
+                            bool just_one_ag = false)
 {
     //save results and sum of results over nodes
     auto ensemble_result_sum = sum_nodes(ensemble_results);
     auto num_groups          = (int)(size_t)ensemble_params[0][0].parameters.get_num_groups();
+    if (just_one_ag) {
+        num_groups = 1;
+    }
     if (save_single_runs) {
         for (size_t i = 0; i < ensemble_result_sum.size(); ++i) {
             BOOST_OUTCOME_TRY(save_result(ensemble_result_sum[i], {0}, num_groups,

cpp/memilio/utils/random_number_generator.h

@@ -182,7 +182,9 @@ static_assert(sizeof(Key<uint32_t>) == sizeof(uint32_t), "Empty Base Optimizatio
 * @tparam an unsigned integer type that determines the size of the counter, i.e., the length of the random sequence.
 */
 template <class T>
-struct MEMILIO_ENABLE_EBO Counter : TypeSafe<T, Counter<T>>, OperatorComparison<Counter<T>>, OperatorAdditionSubtraction<Counter<T>> {
+struct MEMILIO_ENABLE_EBO Counter : TypeSafe<T, Counter<T>>,
+                                    OperatorComparison<Counter<T>>,
+                                    OperatorAdditionSubtraction<Counter<T>> {
     static_assert(std::is_unsigned<T>::value, "Underlying Integer type must be unsigned.");
     using TypeSafe<T, Counter<T>>::TypeSafe;
 };
@@ -243,12 +245,13 @@ Counter<UIntC> rng_totalsequence_counter(UIntN subsequence_idx, CounterS counter
     static_assert(N_BITS <= C_BITS, "Subsequence index must not be bigger than total sequence counter.");
     static_assert(N_BITS <= sizeof(UIntN) * BITS_PER_BYTE, "Subsequence index must be at least N bits");
 
-    assert(UIntC(subsequence_idx) <= (UIntC(1) << N_BITS) && "Subsequence index is too large."); //(1 << N) is the same as (2^N)
+    assert(UIntC(subsequence_idx) <= (UIntC(1) << N_BITS) &&
+           "Subsequence index is too large."); //(1 << N) is the same as (2^N)
 
     //N high bits: subsequence idx
     //S low bits: subsequence counter
     //=> C = N + S bits: total sequence counter
-    //example: 
+    //example:
     //subsequence index uint32_t(181) = 0x000000B5
     //subsequence counter uint32_t(41309) = 0x0000A15D
     //total sequence counter = 0x000000B50000A15D
@@ -271,7 +274,7 @@ Counter<UIntC> rng_totalsequence_counter(UIntN subsequence_idx, CounterS counter
 template <class UIntS, class CounterC>
 Counter<UIntS> rng_subsequence_counter(CounterC counter)
 {
-    using UIntC = typename CounterC::ValueType;
+    using UIntC                = typename CounterC::ValueType;
     static const UIntC C_BYTES = sizeof(UIntC);
     static const UIntC S_BYTES = sizeof(UIntS);
 
@@ -673,6 +676,20 @@ using UniformDistribution = DistributionAdapter<std::uniform_real_distribution<R
 template <class Int>
 using PoissonDistribution = DistributionAdapter<std::poisson_distribution<Int>>;
 
+/**
+ * adapted lognormal_distribution.
+ * @see DistributionAdapter
+ */
+template <class Real>
+using LogNormalDistribution = DistributionAdapter<std::lognormal_distribution<Real>>;
+
+/**
+ * adapted gamma_distribution.
+ * @see DistributionAdapter
+ */
+template <class Real>
+using GammaDistribution = DistributionAdapter<std::gamma_distribution<Real>>;
+
 } // namespace mio
 
 #endif

cpp/models/abm/abm.h

@@ -38,5 +38,6 @@
 #include "abm/vaccine.h"
 #include "abm/household.h"
 #include "abm/lockdown_rules.h"
+#include "abm/analyze_result.h"
 
 #endif