Compare to epicast

.gitignore

@@ -8,3 +8,4 @@ docs/amrex-doxygen-web.tag.xml
 .build
 src/version.h
 reports
+__pycache__

CMakeLists.txt

@@ -14,7 +14,7 @@ set(AMReX_PARTICLES ON)
 set(AMReX_FORTRAN OFF)
 set(AMReX_LINEAR_SOLVERS OFF)
 set(AMReX_AMRLEVEL OFF)
-set(AMReX_TINY_PROFILE ON)
+set(AMReX_TINY_PROFILE OFF)
 set(AMReX_SPACEDIM
     2
     CACHE INTERNAL "")
@@ -25,6 +25,8 @@ set(AMReX_PARTICLES_PRECISION
     SINGLE
     CACHE INTERNAL "")
 
+# add_compile_definitions(COMPARE_TO_EPICAST)
+
 #
 # Fetch amrex repo
 #
@@ -41,7 +43,8 @@ set(FETCHCONTENT_QUIET OFF) # Verbose ON
 FetchContent_Declare(
   amrex
   GIT_REPOSITORY https://github.com/AMReX-Codes/amrex.git/
-  GIT_TAG ${AMReX_GIT_BRANCH})
+  GIT_TAG ${AMReX_GIT_BRANCH}
+  GIT_SUBMODULES "")
 
 if(NOT ${amrex}_POPULATED)
   FetchContent_Populate(amrex)

data/.gitignore

@@ -2,3 +2,11 @@ Backtrace.*
 output-*
 out-*
 results
+EducationData/*
+!EducationData/schools_with_geoids.csv.gz
+LODES7
+US_2010_Census_BlockGroups
+US_2010_Census_States
+UrbanPop
+!UrbanPop/urbanpop_nm.bin.gz
+!UrbanPop/urbanpop_nm.idx

requirements.txt

@@ -1 +1,3 @@
 numpy~=1.15
+scipy
+matplotlib

src/AgentContainer.H

@@ -116,7 +116,7 @@ class AgentContainer : public amrex::ParticleContainer<0, 0, RealIdx::nattribs,
     amrex::GpuArray<amrex::Real, AgeGroups::total> m_symptomatic_withdraw_compliance_day_0 = {0.3_rt, 0.3_rt, 0.3_rt,
                                                                                               0.3_rt, 0.3_rt, 0.3_rt};
     /*!< Symptomatic withdrawal compliance rate, day 1 */
-    amrex::GpuArray<amrex::Real, AgeGroups::total> m_symptomatic_withdraw_compliance_day_1 = {0.8_rt, 0.7_rt, 0.5_rt,
+    amrex::GpuArray<amrex::Real, AgeGroups::total> m_symptomatic_withdraw_compliance_day_1 = {0.8_rt, 0.6_rt, 0.5_rt,
                                                                                               0.5_rt, 0.5_rt, 0.5_rt};
     /*!< Symptomatic withdrawal compliance rate, day 2 */
     amrex::GpuArray<amrex::Real, AgeGroups::total> m_symptomatic_withdraw_compliance_day_2 = {0.9_rt, 0.8_rt, 0.7_rt,

src/AgentContainer.cpp

@@ -720,6 +720,7 @@ void AgentContainer::infectAgents (MFPtrVec& a_disease_stats /*!< Community-wise
                 auto ds_arr = (*a_disease_stats[d])[mfi].array();
 
                 auto status_ptr = soa.GetIntData(i_RT + i0(d) + IntIdxDisease::status).data();
+                auto symptomatic_ptr = soa.GetIntData(i_RT + i0(d) + IntIdxDisease::symptomatic).data();
 
                 auto prob_ptr = soa.GetRealData(r_RT + r0(d) + RealIdxDisease::prob).data();
                 auto counter_ptr = soa.GetRealData(r_RT + r0(d) + RealIdxDisease::disease_counter).data();
@@ -752,9 +753,9 @@ void AgentContainer::infectAgents (MFPtrVec& a_disease_stats /*!< Community-wise
                     prob_ptr[i] *= (1.0_prt - ci_arr[i]);
                     if (status_ptr[i] == Status::never || status_ptr[i] == Status::susceptible) {
                         if (amrex::Random(engine) < prob_ptr[i]) {
-                            setInfected(&(status_ptr[i]), &(counter_ptr[i]), &(latent_period_ptr[i]), &(infectious_period_ptr[i]),
-                                        &(incubation_period_ptr[i]), &(hospital_delay_ptr[i]), &(hospital_random_ptr[i]), engine,
-                                        lparm);
+                            setInfected(&(status_ptr[i]), &(symptomatic_ptr[i]), &(counter_ptr[i]), &(latent_period_ptr[i]),
+                                        &(infectious_period_ptr[i]), &(incubation_period_ptr[i]), &(hospital_delay_ptr[i]),
+                                        &(hospital_random_ptr[i]), engine, lparm);
                             Gpu::Atomic::AddNoRet(&ds_arr(home_i_ptr[i], home_j_ptr[i], 0, DiseaseStats::new_cases), 1.0_rt);
 
                             return;

src/AgentDefinitions.H

@@ -244,7 +244,9 @@ bool isNewlyAsymptomatic (const int a_idx,      /*!< Agent index */
     return status == Status::infected && disease_counter == incb_period && symptom_status == SymptomStatus::asymptomatic;
 }
 
-/*! \brief Did the agent just become presymptomatic this step? */
+/*! \brief Did the agent just transition out of presymptomatic this step (i.e., just became symptomatic or asymptomatic)?
+ *  This fires on the same step as isNewlySymptomatic / isNewlyAsymptomatic, so that
+ *  sum(NewP) == sum(NewS) + sum(NewA) over a complete run. */
 template <typename PTDType>
 AMREX_GPU_DEVICE AMREX_FORCE_INLINE
 bool isNewlyPresymptomatic (const int a_idx,      /*!< Agent index */
@@ -254,8 +256,9 @@ bool isNewlyPresymptomatic (const int a_idx,      /*!< Agent index */
     int symptom_status = a_ptd.m_runtime_idata[i0(a_d) + IntIdxDisease::symptomatic][a_idx];
     int disease_counter = (int)a_ptd.m_runtime_rdata[r0(a_d) + RealIdxDisease::disease_counter][a_idx];
     int incb_period = (int)amrex::Math::floor(a_ptd.m_runtime_rdata[r0(a_d) + RealIdxDisease::incubation_period][a_idx]);
-    // presymptomatic for one day only, just before becoming asymp/symp
-    return status == Status::infected && (disease_counter == incb_period - 1) && symptom_status == SymptomStatus::presymptomatic;
+    // fires on the same step as isNewlySymptomatic/isNewlyAsymptomatic (counter == incb_period, after transition)
+    return status == Status::infected && (disease_counter == incb_period) &&
+           (symptom_status == SymptomStatus::symptomatic || symptom_status == SymptomStatus::asymptomatic);
 }
 
 /*! \brief Is an agent asymptomatic (does not currently have symptoms and will not develop them)? */

src/DiseaseParm.H

@@ -47,25 +47,25 @@ struct DiseaseParm {
 
     // Transmission probabilities based on age group of receiver (0-4, 5-17, 18-29, 30-49, 50-64, 65+)
     /*! community transmission, both home and work */
-    Real xmit_comm[AgeGroups::total] = {Real(0.000015), Real(0.000055), Real(0.00015),
-                                        Real(0.00015),  Real(0.00015),  Real(0.00025)};
+    Real xmit_comm[AgeGroups::total] = {Real(0.000021), Real(0.000062), Real(0.000165),
+                                        Real(0.000165), Real(0.000165), Real(0.000247)};
     /*! neighborhood transmission, both home and work */
-    Real xmit_hood[AgeGroups::total] = {Real(0.000086), Real(0.00025), Real(0.00065), Real(0.00065), Real(0.00065), Real(0.001)};
+    Real xmit_hood[AgeGroups::total] = {Real(0.000082), Real(0.000247), Real(0.00066), Real(0.00066), Real(0.00066), Real(0.001)};
     /*! within household transmission, where transmitter is an adult */
     Real xmit_hh_adult[AgeGroups::total] = {Real(0.3), Real(0.3), Real(0.4), Real(0.4), Real(0.4), Real(0.4)};
     /*! within household transmission, where transmitter is a child */
     Real xmit_hh_child[AgeGroups::total] = {Real(0.6), Real(0.6), Real(0.3), Real(0.3), Real(0.3), Real(0.3)};
     /*! neighborhood cluster transmission, where transmitter is an adult */
-    Real xmit_nc_adult[AgeGroups::total] = {Real(0.0528), Real(0.0528), Real(0.066), Real(0.066), Real(0.066), Real(0.066)};
+    Real xmit_nc_adult[AgeGroups::total] = {Real(0.132), Real(0.132), Real(0.165), Real(0.165), Real(0.165), Real(0.165)};
     /*! neighborhood cluster transmission, where transmitter is a child */
-    Real xmit_nc_child[AgeGroups::total] = {Real(0.1), Real(0.1), Real(0.0528), Real(0.0528), Real(0.0528), Real(0.0528)};
+    Real xmit_nc_child[AgeGroups::total] = {Real(0.25), Real(0.25), Real(0.132), Real(0.132), Real(0.132), Real(0.132)};
     /// probabilities for school groups: none, college, high, middle, elementary, and daycare
     /*! child-to-child */
-    Real xmit_school[SchoolType::total] = {Real(0), Real(0.0315), Real(0.0315), Real(0.0375), Real(0.0435), Real(0.15)};
+    Real xmit_school[SchoolType::total] = {Real(0), Real(0.0002), Real(0.0024), Real(0.0027), Real(0.0088), Real(0.0125)};
     /*! for adult to child */
-    Real xmit_school_a2c[SchoolType::total] = {Real(0), Real(0.0315), Real(0.0315), Real(0.0375), Real(0.0435), Real(0.15)};
+    Real xmit_school_a2c[SchoolType::total] = {Real(0), Real(0.001), Real(0.01), Real(0.013), Real(0.043), Real(0.025)};
     /*! for child to adult */
-    Real xmit_school_c2a[SchoolType::total] = {Real(0), Real(0.0315), Real(0.0315), Real(0.0375), Real(0.0435), Real(0.15)};
+    Real xmit_school_c2a[SchoolType::total] = {Real(0), Real(0.001), Real(0.01), Real(0.013), Real(0.043), Real(0.025)};
 
     /*! Probabilities when school is closed */
     Real xmit_comm_SC[AgeGroups::total];
@@ -78,29 +78,33 @@ struct DiseaseParm {
     /*! probability for transmission within a workgroup, independent of age group currently */
     Real xmit_work = Real(0.0575);
 
-    Real p_trans = Real(0.20);                    /*!< probability of transimission given contact */
-    Real p_asymp = Real(0.30);                    /*!< fraction of cases that are asymptomatic */
-    Real asymp_relative_inf = Real(0.7);          /*!< relative infectiousness of asymptomatic individuals */
+    Real p_trans = Real(0.20);                 /*!< probability of transimission given contact */
+    Real p_asymp = Real(0.30);                 /*!< fraction of cases that are asymptomatic */
+    Real asymp_relative_inf = Real(0.7);       /*!< relative infectiousness of asymptomatic individuals */
 
-    Real vac_eff = Real(0.0);                     /*!< Vaccine efficacy */
+    Real vac_eff = Real(0.0);                  /*!< Vaccine efficacy */
 
-    Real child_compliance;                        /*!< Child compliance with masking ?? */
-    Real child_HH_closure;                        /*!< Multiplier for household contacts during school closure */
+    Real child_compliance;                     /*!< Child compliance with masking ?? */
+    Real child_HH_closure;                     /*!< Multiplier for household contacts during school closure */
 
-    Real immune_length_alpha = Real(540.0);       /*! alpha parameter for gamma distribution*/
-    Real immune_length_beta = Real(0.33);         /*! beta parameter for gamma distribution*/
+    Real immune_length_alpha = Real(1.0);      /*! alpha parameter for gamma distribution*/
+    Real immune_length_beta = Real(1.0);       /*! beta parameter for gamma distribution*/
+    Real immune_length_loc = Real(10000.0);    /*! location (shift) parameter for gamma distribution*/
 
-    Real latent_length_alpha = Real(5.2);         /*! alpha parameter for gamma distribution*/
-    Real latent_length_beta = Real(0.75);         /*! beta parameter for gamma distribution*/
+    Real latent_length_alpha = Real(6.0);      /*! alpha parameter for gamma distribution (unused - see latent_period_cdf)*/
+    Real latent_length_beta = Real(0.73);      /*! beta parameter for gamma distribution (unused - see latent_period_cdf)*/
 
-    Real infectious_length_alpha = Real(26.2);    /*! alpha parameter for gamma distribution*/
-    Real infectious_length_beta = Real(0.23);     /*! beta parameter for gamma distribution*/
+    Real infectious_length_alpha = Real(3.54); /*! alpha parameter for gamma distribution (unused - see infectious_period_cdf)*/
+    Real infectious_length_beta = Real(1.22);  /*! beta parameter for gamma distribution (unused - see infectious_period_cdf)*/
+    Real infectious_length_loc = Real(2.5);    /*! location (shift) parameter for gamma distribution */
 
-    Real incubation_length_alpha = Real(7.5);     /*! alpha parameter for gamma distribution*/
-    Real incubation_length_beta = Real(0.65);     /*! beta parameter for gamma distribution*/
+    Real incubation_length_alpha = Real(6.0);  /*! alpha parameter for gamma distribution (unused - see incubation_period_cdf)*/
+    Real incubation_length_beta = Real(0.73);  /*! beta parameter for gamma distribution (unused - see incubation_period_cdf)*/
+    Real incubation_length_loc = Real(1.0);    /*! location (shift) parameter for gamma distribution */
 
-    Real hospital_delay_length_alpha = Real(1.0); /*! alpha parameter for gamma distribution*/
-    Real hospital_delay_length_beta = Real(1.0);  /*! beta parameter for gamma distribution*/
+    Real hospital_delay_length_alpha = Real(0.1); /*! alpha parameter for gamma distribution*/
+    Real hospital_delay_length_beta = Real(0.1);  /*! beta parameter for gamma distribution*/
+    Real hospital_delay_length_loc = Real(1.0);   /*! location (shift) parameter for gamma distribution*/
 
     int m_hospital_stay_type = HospitalStayType::Constant;
 
@@ -187,6 +191,30 @@ struct DiseaseParm {
     }
 };
 
+#ifdef COMPARE_TO_EPICAST
+/*! Fixed cumulative distribution for latent period from Epicast (fig 6)
+ *  Entry i gives P(transition by day i). Days 0-7. */
+static constexpr int latent_period_cdf_size = 8;
+AMREX_GPU_DEVICE AMREX_GPU_CONSTANT
+constexpr Real latent_period_cdf[latent_period_cdf_size] = {Real(0.0),  Real(0.1),  Real(0.25), Real(0.55),
+                                                            Real(0.66), Real(0.78), Real(0.85), Real(1.0)};
+
+/*! Fixed cumulative distribution for incubation period from Epicast (fig 6)
+ *  Entry i gives P(transition by day i). Days 0-8. */
+static constexpr int incubation_period_cdf_size = 9;
+AMREX_GPU_DEVICE AMREX_GPU_CONSTANT
+constexpr Real incubation_period_cdf[incubation_period_cdf_size] = {Real(0.0),  Real(0.0),  Real(0.1),  Real(0.25), Real(0.55),
+                                                                    Real(0.66), Real(0.78), Real(0.85), Real(1.0)};
+
+/*! Fixed cumulative distribution for infectious period from Epicast (fig 6)
+ *  Entry i gives P(transition by day i). Days 0-10. */
+static constexpr int infectious_period_cdf_size = 11;
+AMREX_GPU_DEVICE AMREX_GPU_CONSTANT
+constexpr Real infectious_period_cdf[infectious_period_cdf_size] = {Real(0.0),  Real(0.0),  Real(0.0), Real(0.0),
+                                                                    Real(0.1),  Real(0.3),  Real(0.5), Real(0.7),
+                                                                    Real(0.85), Real(0.95), Real(1.0)};
+#endif
+
 /*! \brief Will a symptomatic agent be hospitalized? */
 template <typename PTDType>
 AMREX_GPU_DEVICE AMREX_FORCE_INLINE
@@ -198,30 +226,58 @@ bool willBeHospitalized (const int a_idx,      /*!< Agent index */
     return a_parm.checkHospitalization(draw, a_ptd.m_idata[IntIdx::age_group][a_idx]);
 }
 
+/*! \brief Sample a period (in days) from a fixed cumulative distribution array.
+ *
+ *  Draws a uniform random number and returns the first day index where the
+ *  cumulative probability exceeds the draw (inverse CDF / quantile method).
+ *  If the draw exceeds all entries the last valid day index is returned.
+ */
+AMREX_GPU_DEVICE AMREX_FORCE_INLINE
+int sampleFromCDF (const amrex::Real* cdf, int cdf_size, amrex::RandomEngine const& engine) {
+    amrex::Real u = amrex::Random(engine);
+    for (int i = 0; i < cdf_size; ++i) {
+        if (u <= cdf[i]) { return i; }
+    }
+    return cdf_size - 1;
+}
+
 /*! \brief Set this agent to infected status, and initialize disease periods. */
 AMREX_GPU_DEVICE AMREX_FORCE_INLINE
-void setInfected (int* status, amrex::ParticleReal* counter, amrex::ParticleReal* latent_period,
+void setInfected (int* status, int* symptomatic, amrex::ParticleReal* counter, amrex::ParticleReal* latent_period,
                   amrex::ParticleReal* infectious_period, amrex::ParticleReal* incubation_period,
                   amrex::ParticleReal* hospital_delay, amrex::ParticleReal* hospital_random, amrex::RandomEngine const& engine,
                   const DiseaseParm* lparm) {
     *status = Status::infected;
+    *symptomatic = SymptomStatus::presymptomatic;
     *counter = ParticleReal(0);
+
+#ifdef COMPARE_TO_EPICAST
+    // Sample latent, incubation, and infectious periods from fixed Epicast CDFs
+    *latent_period = static_cast<ParticleReal>(sampleFromCDF(latent_period_cdf, latent_period_cdf_size, engine));
+    *incubation_period = static_cast<ParticleReal>(sampleFromCDF(incubation_period_cdf, incubation_period_cdf_size, engine));
+    *infectious_period = static_cast<ParticleReal>(sampleFromCDF(infectious_period_cdf, infectious_period_cdf_size, engine));
+    *hospital_delay = amrex::Real(2);
+#else
     *latent_period = static_cast<ParticleReal>(amrex::RandomGamma(lparm->latent_length_alpha, lparm->latent_length_beta, engine));
     *incubation_period =
-            static_cast<ParticleReal>(amrex::RandomGamma(lparm->incubation_length_alpha, lparm->incubation_length_beta, engine));
+            static_cast<ParticleReal>(amrex::RandomGamma(lparm->incubation_length_alpha, lparm->incubation_length_beta, engine)) +
+            static_cast<ParticleReal>(lparm->incubation_length_loc);
     *infectious_period =
-            static_cast<ParticleReal>(amrex::RandomGamma(lparm->infectious_length_alpha, lparm->infectious_length_beta, engine));
-    *hospital_delay = static_cast<ParticleReal>(
-            amrex::RandomGamma(lparm->hospital_delay_length_alpha, lparm->hospital_delay_length_beta, engine));
-    *hospital_random = static_cast<ParticleReal>(amrex::Random(engine));
+            static_cast<ParticleReal>(amrex::RandomGamma(lparm->infectious_length_alpha, lparm->infectious_length_beta, engine)) +
+            static_cast<ParticleReal>(lparm->infectious_length_loc);
 
+    *hospital_delay = static_cast<ParticleReal>(
+                              amrex::RandomGamma(lparm->hospital_delay_length_alpha, lparm->hospital_delay_length_beta, engine)) +
+                      static_cast<ParticleReal>(lparm->hospital_delay_length_loc);
     if (*latent_period < 1) { *latent_period = amrex::Real(1); }
-    if (*infectious_period < 1) { *infectious_period = amrex::Real(1); }
     if (*incubation_period < 1) { *incubation_period = amrex::Real(1); }
-    if (*hospital_delay < 2) { *hospital_delay = amrex::Real(2); }
-    if (*incubation_period > (*infectious_period + *latent_period)) {
-        *incubation_period = std::floor(*infectious_period + *latent_period);
+    if (*infectious_period < 1) { *infectious_period = amrex::Real(1); }
+#endif
+    if (*incubation_period >= (*infectious_period + *latent_period)) {
+        *incubation_period = std::floor(*infectious_period + *latent_period) - amrex::Real(1);
     }
+
+    *hospital_random = static_cast<ParticleReal>(amrex::Random(engine));
 }
 
 template <typename T>

src/DiseaseParm.cpp

@@ -63,8 +63,13 @@ void DiseaseParm::readInputs (const std::string& a_pp_str /*!< Parmparse string
     pp.query("incubation_length_beta", incubation_length_beta);
     pp.query("hospital_delay_length_beta", hospital_delay_length_beta);
 
+    pp.query("infectious_length_loc", infectious_length_loc);
+    pp.query("incubation_length_loc", incubation_length_loc);
+    pp.query("hospital_delay_length_loc", hospital_delay_length_loc);
+
     pp.query("immune_length_alpha", immune_length_alpha);
     pp.query("immune_length_beta", immune_length_beta);
+    pp.query("immune_length_loc", immune_length_loc);
 
     std::string hospital_stay_type = m_hospital_stay_type == HospitalStayType::Constant ? "constant" : "random";
     pp.query("hospital_stay_type", hospital_stay_type);

src/DiseaseStatus.H

@@ -131,6 +131,7 @@ void DiseaseStatus<AC, ACT, ACTD, A>::updateAgents (AC& a_agents, /*!< Agent con
                 auto* disease_parm_h = a_agents.getDiseaseParameters_h(d);
                 auto immune_length_alpha = disease_parm_h->immune_length_alpha;
                 auto immune_length_beta = disease_parm_h->immune_length_beta;
+                auto immune_length_loc = disease_parm_h->immune_length_loc;
 
                 auto compliance_day_0 = a_agents.m_symptomatic_withdraw_compliance_day_0;
                 auto compliance_day_1 = a_agents.m_symptomatic_withdraw_compliance_day_1;
@@ -192,7 +193,8 @@ void DiseaseStatus<AC, ACT, ACTD, A>::updateAgents (AC& a_agents, /*!< Agent con
                     auto recovered_period = amrex::Math::floor(latent_period_ptr[i] + infectious_period_ptr[i]);
                     if (!marked_for_hosp_ptr[i] && !inHospital(i, ptd) && counter_ptr[i] == recovered_period) {
                         status_ptr[i] = Status::immune;
-                        counter_ptr[i] = static_cast<ParticleReal>(RandomGamma(immune_length_alpha, immune_length_beta, engine));
+                        counter_ptr[i] = static_cast<ParticleReal>(RandomGamma(immune_length_alpha, immune_length_beta, engine)) +
+                                         static_cast<ParticleReal>(immune_length_loc);
                         symptomatic_ptr[i] = SymptomStatus::presymptomatic;
                         withdrawn_ptr[i] = 0;
                     }