1213 PAIS for ABM

cpp/benchmarks/abm.cpp

@@ -76,7 +76,7 @@ mio::abm::Simulation<> make_simulation(size_t num_persons, std::initializer_list
                 mio::UniformIntDistribution<int>::get_instance()(prng, 1, int(mio::abm::InfectionState::Count) - 1));
             auto infection = mio::abm::Infection(prng, mio::abm::VirusVariant::Wildtype, person.get_age(),
                                                  model.parameters, mio::abm::TimePoint(0), state);
-            person.add_new_infection(std::move(infection));
+            person.add_new_infection(std::move(infection), prng, mio::abm::TimePoint(0), model.parameters);
         }
 
         //equal chance of (moderate) mask refusal and (moderate) mask eagerness

cpp/examples/CMakeLists.txt

@@ -116,6 +116,10 @@ add_executable(abm_minimal_example abm_minimal.cpp)
 target_link_libraries(abm_minimal_example PRIVATE memilio abm)
 target_compile_options(abm_minimal_example PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
 
+add_executable(abm_aims_halle abm_aims_halle.cpp)
+target_link_libraries(abm_aims_halle PRIVATE memilio abm)
+target_compile_options(abm_aims_halle PRIVATE ${MEMILIO_CXX_FLAGS_ENABLE_WARNING_ERRORS})
+
 if(MEMILIO_HAS_HDF5)
     add_executable(abm_parameter_study_example abm_parameter_study.cpp)
     target_link_libraries(abm_parameter_study_example PRIVATE memilio abm)

cpp/examples/abm_aims_visual.py

@@ -0,0 +1,13 @@
+import matplotlib.pyplot as plt
+import pandas as pd
+
+df = pd.read_csv("active_pais.txt", index_col=0, parse_dates=True)
+
+plt.figure(figsize=(10, 6))
+for column in df.columns:
+    plt.plot(df.index, df[column], label=column)
+plt.legend(["0-4", "5-14", "15-34", "35-59", "60-79", "80+"], title="Age Groups")
+plt.title("Active PAIS Over Time by Age Group")
+plt.xlabel("Time in Days since 01.01.1970")
+plt.ylabel("Number of Individuals")
+plt.show()

cpp/examples/abm_history_object.cpp

@@ -151,7 +151,8 @@ int main()
             (mio::abm::InfectionState)(rand() % ((uint32_t)mio::abm::InfectionState::Count - 1));
         if (infection_state != mio::abm::InfectionState::Susceptible)
             person.add_new_infection(mio::abm::Infection(rng, mio::abm::VirusVariant::Wildtype, person.get_age(),
-                                                         model.parameters, start_date, infection_state));
+                                                         model.parameters, start_date, infection_state),
+                                     rng, start_date, model.parameters);
     }
 
     // Assign locations to the people

cpp/examples/abm_minimal.cpp

@@ -124,7 +124,8 @@ int main()
         auto rng = mio::abm::PersonalRandomNumberGenerator(model.get_rng(), person);
         if (infection_state != mio::abm::InfectionState::Susceptible) {
             person.add_new_infection(mio::abm::Infection(rng, mio::abm::VirusVariant::Wildtype, person.get_age(),
-                                                         model.parameters, start_date, infection_state));
+                                                         model.parameters, start_date, infection_state),
+                                     rng, start_date, model.parameters);
         }
     }
 

cpp/examples/abm_parameter_study.cpp

@@ -135,7 +135,8 @@ mio::abm::Model make_model(const mio::RandomNumberGenerator& rng)
         auto person_rng = mio::abm::PersonalRandomNumberGenerator(model.get_rng(), person);
         if (infection_state != mio::abm::InfectionState::Susceptible) {
             person.add_new_infection(mio::abm::Infection(person_rng, mio::abm::VirusVariant::Wildtype, person.get_age(),
-                                                         model.parameters, start_date, infection_state));
+                                                         model.parameters, start_date, infection_state),
+                                     person_rng, start_date, model.parameters);
         }
     }
 

cpp/examples/graph_abm.cpp

@@ -204,7 +204,8 @@ int main()
         auto rng = mio::abm::PersonalRandomNumberGenerator(model1.get_rng(), person);
         if (infection_state != mio::abm::InfectionState::Susceptible) {
             person.add_new_infection(mio::abm::Infection(rng, mio::abm::VirusVariant::Wildtype, person.get_age(),
-                                                         model1.parameters, start_date, infection_state));
+                                                         model1.parameters, start_date, infection_state),
+                                     rng, start_date, model1.parameters);
         }
         person.set_assigned_location(mio::abm::LocationType::SocialEvent, event_m1, model1.get_id());
         person.set_assigned_location(mio::abm::LocationType::BasicsShop, shop_m1, model1.get_id());
@@ -234,7 +235,8 @@ int main()
         auto rng = mio::abm::PersonalRandomNumberGenerator(model2.get_rng(), person);
         if (infection_state != mio::abm::InfectionState::Susceptible) {
             person.add_new_infection(mio::abm::Infection(rng, mio::abm::VirusVariant::Wildtype, person.get_age(),
-                                                         model2.parameters, start_date, infection_state));
+                                                         model2.parameters, start_date, infection_state),
+                                     rng, start_date, model2.parameters);
         }
         person.set_assigned_location(mio::abm::LocationType::SocialEvent, event_m2, model2.get_id());
         person.set_assigned_location(mio::abm::LocationType::BasicsShop, shop_m2, model2.get_id());

cpp/models/abm/CMakeLists.txt

@@ -35,6 +35,9 @@ add_library(abm
     mask.h
     mask.cpp
     common_abm_loggers.h
+    sex.h
+    pais.h
+    pais.cpp
 )
 target_link_libraries(abm PUBLIC memilio)
 target_include_directories(abm PUBLIC

cpp/models/abm/common_abm_loggers.h

@@ -178,7 +178,6 @@ struct LogInfectionState : mio::LogAlways {
         Eigen::VectorX<ScalarType> sum =
             Eigen::VectorX<ScalarType>::Zero(Eigen::Index(mio::abm::InfectionState::Count));
         auto curr_time = sim.get_time();
-        PRAGMA_OMP(for)
         for (auto& location : sim.get_model().get_locations()) {
             for (uint32_t inf_state = 0; inf_state < (int)mio::abm::InfectionState::Count; inf_state++) {
                 sum[inf_state] += sim.get_model().get_subpopulation(location.get_id(), curr_time,
@@ -189,6 +188,30 @@ struct LogInfectionState : mio::LogAlways {
     }
 };
 
+/**
+ * @brief Looger to log the TimeSeries of the number of Person%s that have an active PAIS.
+ */
+struct LogPAIS : mio::LogAlways {
+    using Type = std::pair<mio::abm::TimePoint, ScalarType>;
+    /** 
+     * @brief Log the TimeSeries of the number of Person%s that have an active PAIS.
+     * @param[in] sim The simulation of the abm.
+     * @return A pair of the TimePoint and the TimeSeries of the number of Person%s that have an active PAIS.
+     */
+    static Type log(const mio::abm::Simulation<>& sim)
+    {
+        ScalarType sum = 0;
+        auto curr_time = sim.get_time();
+        for (auto& person : sim.get_model().get_persons()) {
+            auto person_id = person.get_id();
+            if (sim.get_model().get_person(person_id).has_active_pais(curr_time)) {
+                sum++;
+            }
+        }
+        return std::make_pair(curr_time, sum);
+    }
+};
+
 /**
 * @brief This is like the DataWriterToMemory, but it only logs time series data.
 * @tparam Loggers The loggers that are used to log data. The loggers must return a touple with a TimePoint and a value.

cpp/models/abm/infection.cpp

@@ -144,6 +144,26 @@ InfectionState Infection::get_infection_state(TimePoint t) const
     return std::prev(it)->second;
 }
 
+std::pair<TimePoint, InfectionState> Infection::get_highest_infection_state() const
+{
+    if (m_infection_course.back().second == InfectionState::Dead) {
+        return m_infection_course.back();
+    }
+    else {
+        return m_infection_course[m_infection_course.size() - 2];
+    }
+}
+
+TimePoint Infection::get_infection_state_start_date(InfectionState state) const
+{
+    for (const auto& [time_point, inf_state] : m_infection_course) {
+        if (inf_state == state) {
+            return time_point;
+        }
+    }
+    return TimePoint(-1); // invalid TimePoint
+}
+
 void Infection::set_detected()
 {
     m_detected = true;

cpp/models/abm/infection.h

@@ -144,6 +144,22 @@ class Infection
      */
     InfectionState get_infection_state(TimePoint t) const;
 
+    /**
+     * @brief Get the highest #InfectionState of the Infection and the TimePoint when that state is reached.
+      * The highest #InfectionState is the state with the most severe symptoms that is reached during the Infection.
+      * For example, if a Person goes through the states InfectedNoSymptoms -> InfectedSymptoms -> Recovered, the highest #InfectionState is InfectedSymptoms.
+      * If a Person goes through the states InfectedNoSymptoms -> InfectedSymptoms -> InfectedSevere -> Recovered, the highest #InfectionState is InfectedSevere.
+      * @return A pair of the highest #InfectionState and the TimePoint when that state is reached.
+     */
+    std::pair<TimePoint, InfectionState> get_highest_infection_state() const;
+
+    /**
+     * @brief Get the start date of a specific #InfectionState.
+     * @param[in] state #InfectionState for which the start date is queried.
+     * @return The start date of the given #InfectionState. If the Person does not reach that state during the Infection, an invalid TimePoint is returned.
+     */
+    TimePoint get_infection_state_start_date(InfectionState state) const;
+
     /**
      * @brief Set the Infection to detected.
      */

cpp/models/abm/model.cpp

@@ -52,10 +52,10 @@ LocationId Model::add_location(LocationType type, uint32_t num_cells)
     return id;
 }
 
-PersonId Model::add_person(const LocationId id, AgeGroup age)
+PersonId Model::add_person(const LocationId id, AgeGroup age, Sex sex)
 {
     PersonId person_id = (static_cast<int64_t>(m_id)) << 32 | static_cast<uint32_t>(m_persons.size());
-    return add_person(Person(m_rng, get_location(id).get_type(), id, m_id, age, person_id));
+    return add_person(Person(m_rng, get_location(id).get_type(), id, m_id, age, sex, person_id));
 }
 
 PersonId Model::add_person(Person&& person)

cpp/models/abm/model.h

@@ -20,6 +20,7 @@
 #ifndef MIO_ABM_MODEL_H
 #define MIO_ABM_MODEL_H
 
+#include "abm/sex.h"
 #include "abm/infection_state.h"
 #include "abm/model_functions.h"
 #include "abm/location_type.h"
@@ -203,9 +204,10 @@ class Model
      * @brief Add a Person to the Model.
      * @param[in] id The LocationID of the initial Location of the Person.
      * @param[in] age AgeGroup of the person.
+     * @param[in] sex Sex of the person.
      * @return Id of the newly created Person.
      */
-    PersonId add_person(const LocationId id, AgeGroup age);
+    PersonId add_person(const LocationId id, AgeGroup age, Sex sex = Sex::Male);
 
     /**
      * @brief Adds a copy of a given Person to the Model.

cpp/models/abm/model_functions.cpp

@@ -100,10 +100,10 @@ void interact(PersonalRandomNumberGenerator& personal_rng, Person& person, const
                 random_transition(personal_rng, VirusVariant::Count, dt,
                                   local_indiv_expected_trans); // use VirusVariant::Count for no virus submission
             if (virus != VirusVariant::Count) {
-                person.add_new_infection(Infection(personal_rng, virus, age_receiver, global_parameters, t + dt / 2,
-                                                   mio::abm::InfectionState::Exposed,
-                                                   person.get_latest_protection(t + dt / 2),
-                                                   false)); // Starting time in second order approximation
+                person.add_new_infection(
+                    Infection(personal_rng, virus, age_receiver, global_parameters, t + dt / 2,
+                              mio::abm::InfectionState::Exposed, person.get_latest_protection(t + dt / 2), false),
+                    personal_rng, t + dt / 2, global_parameters); // Starting time in second order approximation
             }
         }
     }

cpp/models/abm/pais.cpp

@@ -0,0 +1,115 @@
+/*
+* Copyright (C) 2020-2026 MEmilio
+*
+* Authors: David Kerkmann
+*
+* 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.
+*/
+
+#include "abm/pais.h"
+#include "abm/person.h"
+#include "abm/random_events.h"
+
+namespace mio
+{
+namespace abm
+{
+
+void PAIS::update_severity(const Parameters& params, PersonalRandomNumberGenerator& rng, TimePoint t, TimeSpan dt)
+{
+    if (severity.empty() || t > severity.back().first) {
+        return; // only update if the last update was before t
+    }
+    std::pair<PAISState, ScalarType> transmission_probs[static_cast<uint32_t>(PAISState::Count)];
+
+    for (auto&& v : enum_members<PAISState>()) {
+        transmission_probs[static_cast<uint32_t>(v)] = {
+            v, params.get<PAISTransitionMatrix>()(static_cast<Eigen::Index>(severity.back().second),
+                                                  static_cast<Eigen::Index>(v))};
+    }
+    auto severity_new = random_transition(rng, severity.back().second, dt, transmission_probs);
+    if (severity_new != severity.back().second) {
+        this->severity.push_back({t, severity_new}); // only update if there is a change in severity
+    }
+}
+
+void PAIS::init_or_refresh(const Parameters& params, Person& p, PersonalRandomNumberGenerator& rng,
+                           const Infection& inf, TimePoint t)
+{
+    // get highest InfectionState of the new infection
+    auto highest_state = inf.get_highest_infection_state();
+    if (highest_state.second != InfectionState::Dead) {
+        // if the Person already had an active PAIS and gets a reinfection, refresh the PAIS status
+        if (get_severity(t) != PAISState::Count) {
+            add_new_severity(t, highest_state);
+        }
+        else {
+            // base probability of developing PAIS based on age, sex, virus variant and number of vaccinations
+            ScalarType pais_prob = params.get<PAISProbability>()[{inf.get_virus_variant(), p.get_age(), p.get_sex(),
+                                                                  get_vaccination_class(p.get_vaccinations().size())}];
+            // increase probability of developing PAIS if the Person had a severe acute infection or worse
+            if (highest_state.second == InfectionState::InfectedSevere ||
+                highest_state.second == InfectionState::InfectedCritical) {
+                pais_prob *= params.get<PAISProbabilitySeverityFactor>()[{
+                    inf.get_virus_variant(), get_vaccination_class(p.get_vaccinations().size())}];
+            }
+            // reduce probability of developing PAIS if the Person has not had PAIS after an earlier infection
+            if (!p.get_infections().empty() && get_severity(t) == PAISState::Count) {
+                pais_prob *= params.get<PAISProtectionAtSecondInfection>()[{
+                    inf.get_virus_variant(), get_vaccination_class(p.get_vaccinations().size())}];
+            }
+
+            auto& uniform_dist = UniformDistribution<ScalarType>::get_instance();
+            if (uniform_dist(rng) < pais_prob) {
+                TimePoint time_recovered = inf.get_infection_state_start_date(InfectionState::Recovered);
+                add_new_severity(time_recovered, highest_state);
+            }
+        }
+    }
+}
+
+void PAIS::add_new_severity(TimePoint t, std::pair<TimePoint, InfectionState> highest_state)
+{
+    PAISState severity_new;
+    if (highest_state.second == InfectionState::InfectedSevere ||
+        highest_state.second == InfectionState::InfectedCritical) {
+        severity_new = PAISState::Severe;
+    }
+    else {
+        severity_new = PAISState::Medium;
+    }
+    if (severity.empty() || (t > severity.back().first && severity_new != severity.back().second)) {
+        this->severity.push_back({t, severity_new});
+    }
+}
+
+PAISState PAIS::get_severity(TimePoint t) const
+{
+    if (severity.empty()) {
+        return PAISState::Count;
+    }
+    if (t < severity[0].first) {
+        return PAISState::Count;
+    }
+
+    auto it = std::upper_bound(severity.begin(), severity.end(), t,
+                               [](const TimePoint& s, const std::pair<TimePoint, PAISState>& state) {
+                                   return state.first > s;
+                               });
+    return std::prev(it)->second;
+}
+
+} // namespace abm
+} // namespace mio