Finish relative time of first state inclusion

Author: DavidKerkmann

cpp/models/abm/analyze_result.h

@@ -188,7 +188,7 @@ std::vector<Model> ensemble_params_percentile(const std::vector<std::vector<Mode
                         return model.parameters.template get<ViralShedParameters>()[{virus_variant, age_group}];
                     },
                     [](auto& dist1, auto& dist2) {
-                        return dist1.virus_shed_alpha < dist2.virus_shed_alpha;
+                        return dist1.viral_shed_alpha < dist2.viral_shed_alpha;
                     });
                 param_percentile_dist(
                     node, std::vector<mio::AbstractParameterDistribution>(num_runs),

cpp/models/abm/infection.cpp

@@ -29,14 +29,14 @@ namespace abm
 {
 
 void Infection::initialize_viral_load(PersonalRandomNumberGenerator& rng, VirusVariant virus, AgeGroup age,
-                                      const Parameters& params, TimePoint init_date, ProtectionEvent latest_protection)
+                                      const Parameters& params, ProtectionEvent latest_protection)
 {
     auto& vl_params                   = params.get<ViralLoadDistributions>()[{virus, age}];
     ScalarType high_viral_load_factor = 1;
 
     if (latest_protection.type != ProtectionType::NoProtection) {
         high_viral_load_factor -= params.get<HighViralLoadProtectionFactor>()[{latest_protection.type, age, virus}](
-            init_date.days() - latest_protection.time.days());
+            m_viral_load.start_date.days() - latest_protection.time.days());
     }
 
     m_viral_load.peak     = vl_params.viral_load_peak.get(rng) * high_viral_load_factor;
@@ -57,15 +57,56 @@ void Infection::initialize_viral_shed(PersonalRandomNumberGenerator& rng, VirusV
     m_individual_viral_shed_factor = shedfactor_param.get(rng);
 }
 
+TimePoint Infection::shift_init_date(const InitialInfectionStateDistribution& init_state_dist,
+                                     PersonalRandomNumberGenerator& rng, VirusVariant virus, AgeGroup age,
+                                     TimePoint init_date)
+{
+    auto dist              = init_state_dist[{virus, age}];
+    TimeSpan time_in_state = days(dist.get(rng));
+    return init_date - time_in_state;
+}
+
 Infection::Infection(PersonalRandomNumberGenerator& rng, VirusVariant virus, AgeGroup age, const Parameters& params,
                      TimePoint init_date, InfectionState init_state, ProtectionEvent latest_protection, bool detected)
     : m_virus_variant(virus)
     , m_detected(detected)
 {
     assert(age.get() < params.get_num_groups());
-    m_viral_load.start_date = draw_infection_course(rng, age, params, init_date, init_state, latest_protection);
+    assert(init_state != InfectionState::Susceptible &&
+           "Initializatin of an Infection must happen with an InfectionState that is not Suscpetible.");
+
+    draw_infection_course_forward(rng, age, params, init_date, init_state, latest_protection);
+    m_viral_load.start_date = draw_infection_course_backward(rng, age, params, init_date, init_state);
+
+    initialize_viral_load(rng, virus, age, params, latest_protection);
+    initialize_viral_shed(rng, virus, age, params);
+}
+
+Infection::Infection(PersonalRandomNumberGenerator& rng, VirusVariant virus, AgeGroup age, const Parameters& params,
+                     TimePoint init_date, InfectionState init_state,
+                     const InitialInfectionStateDistribution& init_state_dist, ProtectionEvent latest_protection,
+                     bool detected)
+    : m_virus_variant(virus)
+    , m_detected(detected)
+{
+    assert(age.get() < params.get_num_groups());
+    assert(init_state != InfectionState::Susceptible &&
+           "Initializatin of an Infection must happen with an InfectionState that is not Suscpetible.");
+
+    // Draw the first transition and time that the agent has already spent in that state
+    StateTransition first_transition =
+        get_forward_transition(rng, age, params, init_state, init_date, latest_protection);
+    ScalarType relative_time_in_first_state = init_state_dist[{virus, age}].get(rng);
+
+    init_date -= first_transition.duration.multiply(relative_time_in_first_state);
+    m_infection_course.push_back({init_date, first_transition.from_state});
 
-    initialize_viral_load(rng, virus, age, params, init_date, latest_protection);
+    // Draw the rest of the infection course
+    draw_infection_course_forward(rng, age, params, init_date + first_transition.duration, first_transition.to_state,
+                                  latest_protection);
+    m_viral_load.start_date = draw_infection_course_backward(rng, age, params, init_date, init_state);
+
+    initialize_viral_load(rng, virus, age, params, latest_protection);
     initialize_viral_shed(rng, virus, age, params);
 }
 
@@ -131,9 +172,8 @@ TimePoint Infection::draw_infection_course(PersonalRandomNumberGenerator& rng, A
                                            TimePoint init_date, InfectionState init_state,
                                            ProtectionEvent latest_protection)
 {
-    TimePoint start_of_init_state =
-        draw_infection_course_forward(rng, age, params, init_date, init_state, latest_protection);
-    TimePoint start_date = draw_infection_course_backward(rng, age, params, start_of_init_state, init_state);
+    draw_infection_course_forward(rng, age, params, init_date, init_state, latest_protection);
+    TimePoint start_date = draw_infection_course_backward(rng, age, params, init_date, init_state);
     return start_date;
 }
 
@@ -145,10 +185,6 @@ StateTransition Infection::get_forward_transition(PersonalRandomNumberGenerator&
     StateTransition transition{current_state, current_state, TimeSpan{}};
 
     switch (current_state) {
-    case InfectionState::Susceptible:
-        transition.to_state = InfectionState::Exposed;
-        transition.duration = mio::abm::hours(0);
-        break;
 
     case InfectionState::Exposed:
         transition.to_state = InfectionState::InfectedNoSymptoms;
@@ -192,14 +228,14 @@ StateTransition Infection::get_forward_transition(PersonalRandomNumberGenerator&
         ScalarType critical_prob = params.get<CriticalPerInfectedSevere>()[{m_virus_variant, age}];
         ScalarType death_prob    = params.get<DeathsPerInfectedSevere>()[{m_virus_variant, age}];
 
-        if (p < critical_prob) {
-            transition.to_state = InfectionState::InfectedCritical;
-            transition.duration = days(params.get<TimeInfectedSevereToCritical>()[{m_virus_variant, age}].get(rng));
-        }
-        else if (p < critical_prob + death_prob) {
+        if (p < death_prob) {
             transition.to_state = InfectionState::Dead;
             transition.duration = days(params.get<TimeInfectedSevereToDead>()[{m_virus_variant, age}].get(rng));
         }
+        else if (p < critical_prob + death_prob) {
+            transition.to_state = InfectionState::InfectedCritical;
+            transition.duration = days(params.get<TimeInfectedSevereToCritical>()[{m_virus_variant, age}].get(rng));
+        }
         else {
             transition.to_state = InfectionState::Recovered;
             transition.duration = days(params.get<TimeInfectedSevereToRecovered>()[{m_virus_variant, age}].get(rng));
@@ -276,7 +312,7 @@ StateTransition Infection::get_recovered_backward_transition(PersonalRandomNumbe
 
     // Compute death probability to factor it out
     ScalarType p_death = calculate_death_probability(age, params);
-    assert(p_death == 1 && "Trying to create a recovered agent although the chance to die is 100%.");
+    assert(p_death < 1 && "Trying to create a recovered agent although the chance to die is 100%.");
     ScalarType inv_death = 1 / (1 - p_death);
 
     ScalarType symptoms_prob = params.get<SymptomsPerInfectedNoSymptoms>()[{m_virus_variant, age}];
@@ -343,40 +379,29 @@ ScalarType Infection::get_severity_protection_factor(const Parameters& params, P
         current_time.days() - latest_protection.time.days());
 }
 
-TimePoint Infection::draw_infection_course_forward(PersonalRandomNumberGenerator& rng, AgeGroup age,
-                                                   const Parameters& params, TimePoint init_date,
-                                                   InfectionState start_state, ProtectionEvent latest_protection)
+void Infection::draw_infection_course_forward(PersonalRandomNumberGenerator& rng, AgeGroup age,
+                                              const Parameters& params, TimePoint init_date, InfectionState start_state,
+                                              ProtectionEvent latest_protection)
 {
-    TimePoint start_of_init_state = init_date; // since there is no state transition from Recovered or Dead
-    bool init                     = true; // the random start time of the first state cannot be drawn
-    // thus, just use the init_date as the start for this
-    auto& init_state_dist = params.get<InitialInfectionStateDistributions>()[{m_virus_variant, age, start_state}];
-
     TimePoint current_time       = init_date;
     InfectionState current_state = start_state;
     m_infection_course.push_back({current_time, current_state});
 
     while (current_state != InfectionState::Recovered && current_state != InfectionState::Dead) {
         StateTransition transition =
             get_forward_transition(rng, age, params, current_state, current_time, latest_protection);
-        if (init && current_state != InfectionState::Susceptible) { // random init within first time period
-            ScalarType p = init_state_dist.get(rng);
-            start_of_init_state -= transition.duration.multiply(1 - p);
-            transition.duration = transition.duration.multiply(p);
-            init                = false;
-        }
+
         current_time += transition.duration;
         current_state = transition.to_state;
         m_infection_course.push_back({current_time, current_state});
     }
-    return start_of_init_state;
 }
 
 TimePoint Infection::draw_infection_course_backward(PersonalRandomNumberGenerator& rng, AgeGroup age,
-                                                    const Parameters& params, TimePoint start_of_init_state,
+                                                    const Parameters& params, TimePoint init_date,
                                                     InfectionState init_state)
 {
-    TimePoint current_time       = start_of_init_state;
+    TimePoint current_time       = init_date;
     InfectionState current_state = init_state;
 
     while (current_state != InfectionState::Exposed) {

cpp/models/abm/infection.h

@@ -67,6 +67,13 @@ struct ViralLoad {
     }
 };
 
+/**
+ * @brief Distributions of the relative time that people have been in their initial infection state at the beginning of the simulation.
+ * Values have to be within [0, 1].
+ * This makes it possible to draw from a user-defined distribution instead of drawing from a uniform distribution.
+ */
+using InitialInfectionStateDistribution = CustomIndexArray<AbstractParameterDistribution, VirusVariant, AgeGroup>;
+
 class Infection
 {
 public:
@@ -83,7 +90,23 @@ class Infection
      * @param[in] detected [Default: false] If the Infection is detected.     
      */
     Infection(PersonalRandomNumberGenerator& rng, VirusVariant virus, AgeGroup age, const Parameters& params,
-              TimePoint start_date, InfectionState start_state = InfectionState::Susceptible,
+              TimePoint start_date, InfectionState start_state = InfectionState::Exposed,
+              ProtectionEvent latest_protection = {ProtectionType::NoProtection, TimePoint(0)}, bool detected = false);
+
+    /**
+     * @brief Create an Infection for a single Person with a time spent in the given initial state that is drawn from the given distribution.
+     * @param[inout] rng PersonalRandomNumberGenerator of the Person.
+     * @param[in] virus Virus type of the Infection.
+     * @param[in] age AgeGroup to determine the ViralLoad course.
+     * @param[in] params Parameters of the Model.
+     * @param[in] init_date Date of initializing the Infection.
+     * @param[in] init_state #InfectionState at time of initializing the Infection.
+     * @param[in] init_state_dist Distribution to draw the relative time spent in the initial state from. Values have to be within [0, 1].
+     * @param[in] latest_protection The pair value of last ProtectionType (previous Infection/Vaccination) and TimePoint of that protection.
+     * @param[in] detected If the Infection is detected.
+     */
+    Infection(PersonalRandomNumberGenerator& rng, VirusVariant virus, AgeGroup age, const Parameters& params,
+              TimePoint init_date, InfectionState init_state, const InitialInfectionStateDistribution& init_state_dist,
               ProtectionEvent latest_protection = {ProtectionType::NoProtection, TimePoint(0)}, bool detected = false);
 
     /**
@@ -178,20 +201,16 @@ class Infection
      * InfectedCritical -> Recovered or InfectedCritical -> Dead,
      * until either Recoverd or Dead is reached.
      * The duration in each #InfectionState is taken from the respective parameter.
-     * The first transition has a random, uniformly distributed length from 0 to the respective parameter
-     * for persons that are initialized somewhere in the middle of the infection, i.e. not a newly infected
-     * person nor a Recovered/Dead person, to reflect uncertainty in the current state.
      * @param[inout] rng PersonalRandomNumberGenerator of the Person.
      * @param[in] age AgeGroup of the Person.
      * @param[in] params Parameters of the Model.
      * @param[in] init_date Date of initializing the Infection.
      * @param[in] init_state #InfectionState at time of initializing the Infection.
      * @param[in] latest_protection Latest protection against Infection, has an influence on transition probabilities.
-     * @return The starting date of the init_state.
      */
-    TimePoint draw_infection_course_forward(PersonalRandomNumberGenerator& rng, AgeGroup age, const Parameters& params,
-                                            TimePoint init_date, InfectionState init_state,
-                                            ProtectionEvent latest_protection);
+    void draw_infection_course_forward(PersonalRandomNumberGenerator& rng, AgeGroup age, const Parameters& params,
+                                       TimePoint init_date, InfectionState init_state,
+                                       ProtectionEvent latest_protection);
 
     /**
      * @brief Determine Infection course prior to the given #InfectionState start_state.
@@ -213,11 +232,10 @@ class Infection
      * @param[in] virus Virus type of the Infection.
      * @param[in] age AgeGroup of the Person.
      * @param[in] params Parameters of the Model.
-     * @param[in] init_date Date of initializing the Infection.
      * @param[in] latest_protection Latest protection against Infection.
      */
     void initialize_viral_load(PersonalRandomNumberGenerator& rng, VirusVariant virus, AgeGroup age,
-                               const Parameters& params, TimePoint init_date, ProtectionEvent latest_protection);
+                               const Parameters& params, ProtectionEvent latest_protection);
 
     /**
      * @brief Initialize the viral shed parameters and individual factor for the infection.
@@ -229,6 +247,18 @@ class Infection
     void initialize_viral_shed(PersonalRandomNumberGenerator& rng, VirusVariant virus, AgeGroup age,
                                const Parameters& params);
 
+    /**
+    * @brief Adjust start date based on initial infection state distribution
+    * @param[in] init_state_dist Distribution to draw the time spent in the initial state from.
+    * @param[inout] rng PersonalRandomNumberGenerator of the Person.
+    * @param[in] virus Virus type of the Infection.
+    * @param[in] age AgeGroup of the Person.
+    * @param[in] init_date Date of initializing the Infection.
+    */
+    TimePoint shift_init_date(const InitialInfectionStateDistribution& init_state_dist,
+                              PersonalRandomNumberGenerator& rng, VirusVariant virus, AgeGroup age,
+                              TimePoint init_date);
+
     /**
      * @brief Get the forward transition from a given infection state.
      * @param[inout] rng PersonalRandomNumberGenerator of the Person.

cpp/models/abm/model_functions.cpp

@@ -101,7 +101,7 @@ void interact(PersonalRandomNumberGenerator& personal_rng, Person& person, const
                                   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::Susceptible,
+                                                   mio::abm::InfectionState::Exposed,
                                                    person.get_latest_protection(t + dt / 2),
                                                    false)); // Starting time in second order approximation
             }

cpp/models/abm/parameters.h

@@ -277,23 +277,6 @@ struct DeathsPerInfectedCritical {
     }
 };
 
-/**
- * @brief Distributions of the time that people have been in their initial infection state at the beginning of the simulation.
- * This makes it possible to draw from a user-defined distribution instead of drawing from a uniform distribution.
- */
-struct InitialInfectionStateDistributions {
-    using Type = CustomIndexArray<AbstractParameterDistribution, VirusVariant, AgeGroup, InfectionState>;
-    static Type get_default(AgeGroup size)
-    {
-        return Type({VirusVariant::Count, size, InfectionState::Count},
-                    AbstractParameterDistribution(ParameterDistributionUniform(0., 1.)));
-    }
-    static std::string name()
-    {
-        return "InitialInfectionStateDistributions";
-    }
-};
-
 /**
  * @brief Parameters for the ViralLoad course. Default values taken as constant values from the average from
  * https://github.com/VirologyCharite/SARS-CoV-2-VL-paper/tree/main
@@ -718,13 +701,12 @@ using ParametersBase =
                  TimeInfectedSymptomsToSevere, TimeInfectedSymptomsToRecovered, TimeInfectedSevereToCritical,
                  TimeInfectedSevereToRecovered, TimeInfectedSevereToDead, TimeInfectedCriticalToDead,
                  TimeInfectedCriticalToRecovered, SymptomsPerInfectedNoSymptoms, SeverePerInfectedSymptoms,
-                 CriticalPerInfectedSevere, DeathsPerInfectedSevere, DeathsPerInfectedCritical,
-                 InitialInfectionStateDistributions, ViralLoadDistributions, ViralShedParameters, ViralShedFactor,
-                 InfectionRateFromViralShed, MaskProtection, AerosolTransmissionRates, LockdownDate, QuarantineDuration,
-                 QuarantineEffectiveness, SocialEventRate, BasicShoppingRate, WorkRatio, SchoolRatio,
-                 GotoWorkTimeMinimum, GotoWorkTimeMaximum, GotoSchoolTimeMinimum, GotoSchoolTimeMaximum,
-                 AgeGroupGotoSchool, AgeGroupGotoWork, InfectionProtectionFactor, SeverityProtectionFactor,
-                 HighViralLoadProtectionFactor, TestData>;
+                 CriticalPerInfectedSevere, DeathsPerInfectedSevere, DeathsPerInfectedCritical, ViralLoadDistributions,
+                 ViralShedParameters, ViralShedFactor, InfectionRateFromViralShed, MaskProtection,
+                 AerosolTransmissionRates, LockdownDate, QuarantineDuration, QuarantineEffectiveness, SocialEventRate,
+                 BasicShoppingRate, WorkRatio, SchoolRatio, GotoWorkTimeMinimum, GotoWorkTimeMaximum,
+                 GotoSchoolTimeMinimum, GotoSchoolTimeMaximum, AgeGroupGotoSchool, AgeGroupGotoWork,
+                 InfectionProtectionFactor, SeverityProtectionFactor, HighViralLoadProtectionFactor, TestData>;
 
 /**
  * @brief Maximum number of Person%s an infectious Person can infect at the respective Location.