From 12526c7d1fbbd63bf791931dfe3f71023701878e Mon Sep 17 00:00:00 2001 From: Debojyoti Ghosh Date: Tue, 12 Nov 2024 09:05:06 -0800 Subject: [PATCH] Fix Disease Status (#88) * if -> else if so that an agent doesn't enter more than one if-block * renamed isHospitalized() and flag_hosp, etc for clarity * making sure latent, incubation, & infection periods are never negative; also ensuring latent period is not greater than the incubation+infection period * simplified the condition for latent period being smaller than incubation+infectious periods * not using std::max since Clang complains! * incubation period should be less than the latent+infectious periods! * fixed the if condition for when hospitalization starts --- src/AgentContainer.cpp | 16 +++++++---- src/AgentDefinitions.H | 4 +-- src/DiseaseStatus.H | 46 ++++++++++++++++---------------- src/HospitalModel.H | 10 +++---- src/InteractionModHome.H | 2 +- src/InteractionModHomeNborhood.H | 2 +- src/InteractionModSchool.H | 2 +- src/InteractionModWork.H | 2 +- src/InteractionModWorkNborhood.H | 2 +- 9 files changed, 46 insertions(+), 40 deletions(-) diff --git a/src/AgentContainer.cpp b/src/AgentContainer.cpp index 1ccf082..6804cf8 100644 --- a/src/AgentContainer.cpp +++ b/src/AgentContainer.cpp @@ -174,7 +174,7 @@ void AgentContainer::moveAgentsToWork () amrex::ParallelFor( np, [=] AMREX_GPU_DEVICE (int ip) noexcept { - if (!isHospitalized(ip, ptd)) { + if (!inHospital(ip, ptd)) { ParticleType& p = pstruct[ip]; p.pos(0) = (work_i_ptr[ip] + 0.5_prt)*dx[0]; p.pos(1) = (work_j_ptr[ip] + 0.5_prt)*dx[1]; @@ -222,7 +222,7 @@ void AgentContainer::moveAgentsToHome () amrex::ParallelFor( np, [=] AMREX_GPU_DEVICE (int ip) noexcept { - if (!isHospitalized(ip, ptd)) { + if (!inHospital(ip, ptd)) { ParticleType& p = pstruct[ip]; p.pos(0) = (home_i_ptr[ip] + 0.5_prt) * dx[0]; p.pos(1) = (home_j_ptr[ip] + 0.5_prt) * dx[1]; @@ -269,7 +269,7 @@ void AgentContainer::moveRandomTravel (const amrex::Real random_travel_prob) amrex::ParallelForRNG( np, [=] AMREX_GPU_DEVICE (int i, RandomEngine const& engine) noexcept { - if (!isHospitalized(i, ptd) && !withdrawn_ptr[i]) { + if (!inHospital(i, ptd) && !withdrawn_ptr[i]) { ParticleType& p = pstruct[i]; if (amrex::Random(engine) < random_travel_prob) { random_travel_ptr[i] = i; @@ -321,7 +321,7 @@ void AgentContainer::moveAirTravel (const iMultiFab& unit_mf, AirTravelFlow& air [=] AMREX_GPU_DEVICE (int i, RandomEngine const& engine) noexcept { int unit = unit_arr(home_i_ptr[i], home_j_ptr[i], 0); - if (!isHospitalized(i, ptd) && random_travel_ptr[i] <0 && air_travel_ptr[i] <0) { + if (!inHospital(i, ptd) && random_travel_ptr[i] <0 && air_travel_ptr[i] <0) { if (withdrawn_ptr[i] == 1) {return ;} if (amrex::Random(engine) < air_travel_prob_ptr[unit]) { ParticleType& p = pstruct[i]; @@ -548,7 +548,7 @@ void AgentContainer::updateStatus ( MFPtrVec& a_disease_stats /*!< Community-wis amrex::ParallelFor( np, [=] AMREX_GPU_DEVICE (int ip) noexcept { - if (isHospitalized(ip, ptd)) { + if (inHospital(ip, ptd)) { ParticleType& p = pstruct[ip]; p.pos(0) = (hosp_i_ptr[ip] + 0.5_prt)*dx[0]; p.pos(1) = (hosp_j_ptr[ip] + 0.5_prt)*dx[1]; @@ -674,6 +674,12 @@ void AgentContainer::infectAgents () latent_period_ptr[i] = amrex::RandomNormal(lparm->latent_length_mean, lparm->latent_length_std, engine); infectious_period_ptr[i] = amrex::RandomNormal(lparm->infectious_length_mean, lparm->infectious_length_std, engine); incubation_period_ptr[i] = amrex::RandomNormal(lparm->incubation_length_mean, lparm->incubation_length_std, engine); + if (latent_period_ptr[i] < 0) { latent_period_ptr[i] = 0.0_rt; } + if (infectious_period_ptr[i] < 0) { infectious_period_ptr[i] = 0.0_rt; } + if (incubation_period_ptr[i] < 0) { incubation_period_ptr[i] = 0.0_rt; } + if (incubation_period_ptr[i] > (infectious_period_ptr[i]+latent_period_ptr[i])) { + incubation_period_ptr[i] = std::floor(infectious_period_ptr[i]+latent_period_ptr[i]); + } return; } } diff --git a/src/AgentDefinitions.H b/src/AgentDefinitions.H index 1ad9d25..bbe3e92 100644 --- a/src/AgentDefinitions.H +++ b/src/AgentDefinitions.H @@ -208,8 +208,8 @@ bool notSusceptible ( const int a_idx, /*!< Agent index */ /*! \brief Is an agent hospitalized? */ template AMREX_GPU_DEVICE AMREX_FORCE_INLINE -bool isHospitalized ( const int a_idx, /*!< Agent index */ - const PTDType& a_ptd /*!< Particle tile data */ ) +bool inHospital ( const int a_idx, /*!< Agent index */ + const PTDType& a_ptd /*!< Particle tile data */ ) { return ( (a_ptd.m_idata[IntIdx::hosp_i][a_idx] >= 0) && (a_ptd.m_idata[IntIdx::hosp_j][a_idx] >= 0) ); diff --git a/src/DiseaseStatus.H b/src/DiseaseStatus.H index 91fafc1..b23ee0f 100644 --- a/src/DiseaseStatus.H +++ b/src/DiseaseStatus.H @@ -117,18 +117,18 @@ void DiseaseStatus::updateAgents(AC& a_agents, /*!< Agent contain int i_RT = IntIdx::nattribs; int r_RT = RealIdx::nattribs; - Gpu::DeviceVector flag_hosp, flag_ICU, flag_vent; - flag_hosp.resize(np); - flag_ICU.resize(np); - flag_vent.resize(np); - auto flag_hosp_ptr = flag_hosp.data(); - auto flag_ICU_ptr = flag_ICU.data(); - auto flag_vent_ptr = flag_vent.data(); + Gpu::DeviceVector marked_for_hosp, marked_for_ICU, marked_for_vent; + marked_for_hosp.resize(np); + marked_for_ICU.resize(np); + marked_for_vent.resize(np); + auto marked_for_hosp_ptr = marked_for_hosp.data(); + auto marked_for_ICU_ptr = marked_for_ICU.data(); + auto marked_for_vent_ptr = marked_for_vent.data(); ParallelFor( np, [=] AMREX_GPU_DEVICE (int i) noexcept { - flag_hosp_ptr[i] = 0; - flag_ICU_ptr[i] = 0; - flag_vent_ptr[i] = 0; + marked_for_hosp_ptr[i] = 0; + marked_for_ICU_ptr[i] = 0; + marked_for_vent_ptr[i] = 0; }); Gpu::synchronize(); @@ -175,7 +175,7 @@ void DiseaseStatus::updateAgents(AC& a_agents, /*!< Agent contain symptomatic_ptr[i] = SymptomStatus::presymptomatic; } } - if (counter_ptr[i] == Math::floor(incubation_period_ptr[i])) { + else if (counter_ptr[i] == Math::floor(incubation_period_ptr[i])) { if (symptomatic_ptr[i] != SymptomStatus::asymptomatic) { symptomatic_ptr[i] = SymptomStatus::symptomatic; } @@ -186,21 +186,21 @@ void DiseaseStatus::updateAgents(AC& a_agents, /*!< Agent contain } if (symptomatic_ptr[i] == SymptomStatus::symptomatic) { - int flag_ICU_i = 0, flag_vent_i = 0; + int marked_for_ICU_i = 0, marked_for_vent_i = 0; Real num_days = 0; disease_parm_d->check_hospitalization( num_days, - flag_ICU_i, - flag_vent_i, + marked_for_ICU_i, + marked_for_vent_i, age_group_ptr[i], u50frac, engine ); timer_ptr[i] = ParticleReal(num_days); - if (timer_ptr[i] > 0) { flag_hosp_ptr[i] = 1; } - if (flag_ICU_i) { flag_ICU_ptr[i] = 1; } - if (flag_vent_i) { flag_vent_ptr[i] = 1; } + if (timer_ptr[i] > 0) { marked_for_hosp_ptr[i] = 1; } + if (marked_for_ICU_i) { marked_for_ICU_ptr[i] = 1; } + if (marked_for_vent_i) { marked_for_vent_ptr[i] = 1; } } } - if (!isHospitalized(i,ptd)) { + else if (!inHospital(i,ptd)) { if (counter_ptr[i] >= (latent_period_ptr[i] + infectious_period_ptr[i])) { status_ptr[i] = Status::immune; counter_ptr[i] = amrex::RandomNormal(immune_length_mean, immune_length_std, engine); @@ -216,8 +216,8 @@ void DiseaseStatus::updateAgents(AC& a_agents, /*!< Agent contain ParallelFor( np, [=] AMREX_GPU_DEVICE (int i) noexcept { - if (flag_hosp_ptr[i] == 1) { - AMREX_ALWAYS_ASSERT(!isHospitalized(i, ptd)); + if (marked_for_hosp_ptr[i] == 1) { + AMREX_ALWAYS_ASSERT(!inHospital(i, ptd)); assign_hospital( i, hosp_i_ptr, hosp_j_ptr, ptd); } }); @@ -226,21 +226,21 @@ void DiseaseStatus::updateAgents(AC& a_agents, /*!< Agent contain auto ds_arr = (*a_stats[d])[mfi].array(); ParallelFor( np, [=] AMREX_GPU_DEVICE (int i) noexcept { - if (flag_hosp_ptr[i] == 1) { + if (marked_for_hosp_ptr[i] == 1) { Gpu::Atomic::AddNoRet( &ds_arr( home_i_ptr[i], home_j_ptr[i], 0, DiseaseStats::hospitalization ), 1.0_rt ); } - if (flag_ICU_ptr[i] == 1) { + if (marked_for_ICU_ptr[i] == 1) { Gpu::Atomic::AddNoRet( &ds_arr( home_i_ptr[i], home_j_ptr[i], 0, DiseaseStats::ICU ), 1.0_rt ); } - if (flag_vent_ptr[i] == 1) { + if (marked_for_vent_ptr[i] == 1) { Gpu::Atomic::AddNoRet( &ds_arr( home_i_ptr[i], home_j_ptr[i], 0, diff --git a/src/HospitalModel.H b/src/HospitalModel.H index 0223a60..c4d770e 100644 --- a/src/HospitalModel.H +++ b/src/HospitalModel.H @@ -85,13 +85,13 @@ void HospitalModel::treatAgents(PCType& a_agents, /*!< A int r_RT = RealIdx::nattribs; GpuArray status_ptrs, symptomatic_ptrs; - GpuArray counter_ptrs, timer_ptrs, latent_per_ptrs; + GpuArray counter_ptrs, timer_ptrs, incubation_per_ptrs; for (int d = 0; d < n_disease; d++) { status_ptrs[d] = soa.GetIntData(i_RT+i0(d)+IntIdxDisease::status).data(); symptomatic_ptrs[d] = soa.GetIntData(i_RT+i0(d)+IntIdxDisease::symptomatic).data(); counter_ptrs[d] = soa.GetRealData(r_RT+r0(d)+RealIdxDisease::disease_counter).data(); timer_ptrs[d] = soa.GetRealData(r_RT+r0(d)+RealIdxDisease::treatment_timer).data(); - latent_per_ptrs[d] = soa.GetRealData(r_RT+r0(d)+RealIdxDisease::latent_period).data(); + incubation_per_ptrs[d] = soa.GetRealData(r_RT+r0(d)+RealIdxDisease::incubation_period).data(); } Gpu::DeviceVector is_alive; @@ -132,11 +132,11 @@ void HospitalModel::treatAgents(PCType& a_agents, /*!< A [=] AMREX_GPU_DEVICE (int i, RandomEngine const& engine) noexcept { - if ( !isHospitalized(i, ptd) ) { + if ( !inHospital(i, ptd) ) { // agent is not in hospital return; } - if (counter_ptrs[d][i] == Math::ceil(latent_per_ptrs[d][i])) { + if (counter_ptrs[d][i] == Math::floor(incubation_per_ptrs[d][i])) { // agent just started treatment return; } @@ -213,7 +213,7 @@ void HospitalModel::treatAgents(PCType& a_agents, /*!< A ParallelFor( np, [=] AMREX_GPU_DEVICE (int i) noexcept { - if ( !isHospitalized(i, ptd) ) { return; } + if ( !inHospital(i, ptd) ) { return; } if (is_alive_ptr[i] == 0) { diff --git a/src/InteractionModHome.H b/src/InteractionModHome.H index c1b288c..13d4be1 100644 --- a/src/InteractionModHome.H +++ b/src/InteractionModHome.H @@ -59,7 +59,7 @@ template struct HomeCandidate { AMREX_GPU_HOST_DEVICE bool operator() (const int idx, const PTDType& ptd) const noexcept { - return !isHospitalized(idx, ptd) && ptd.m_idata[IntIdx::random_travel][idx] < 0 && ptd.m_idata[IntIdx::air_travel][idx] < 0; + return !inHospital(idx, ptd) && ptd.m_idata[IntIdx::random_travel][idx] < 0 && ptd.m_idata[IntIdx::air_travel][idx] < 0; } }; diff --git a/src/InteractionModHomeNborhood.H b/src/InteractionModHomeNborhood.H index 996d181..cb54397 100644 --- a/src/InteractionModHomeNborhood.H +++ b/src/InteractionModHomeNborhood.H @@ -47,7 +47,7 @@ struct HomeNborhoodCandidate { AMREX_GPU_HOST_DEVICE bool operator() (const int idx, const PTDType& ptd) const noexcept { // this is the only case where we allow random travelers to interact - return !isHospitalized(idx, ptd) && !ptd.m_idata[IntIdx::withdrawn][idx]; + return !inHospital(idx, ptd) && !ptd.m_idata[IntIdx::withdrawn][idx]; } }; diff --git a/src/InteractionModSchool.H b/src/InteractionModSchool.H index 889a3f1..f774a61 100644 --- a/src/InteractionModSchool.H +++ b/src/InteractionModSchool.H @@ -69,7 +69,7 @@ template struct SchoolCandidate { AMREX_GPU_HOST_DEVICE bool operator() (const int idx, const PTDType& ptd) const noexcept { - return !isHospitalized(idx, ptd) && + return !inHospital(idx, ptd) && ptd.m_idata[IntIdx::school][idx] > 0 && !ptd.m_idata[IntIdx::withdrawn][idx] && ptd.m_idata[IntIdx::air_travel][idx] < 0 && diff --git a/src/InteractionModWork.H b/src/InteractionModWork.H index 50275db..eb94c20 100644 --- a/src/InteractionModWork.H +++ b/src/InteractionModWork.H @@ -57,7 +57,7 @@ template struct WorkCandidate { AMREX_GPU_HOST_DEVICE bool operator() (const int idx, const PTDType& ptd) const noexcept { - return !isHospitalized(idx, ptd) && + return !inHospital(idx, ptd) && ptd.m_idata[IntIdx::work_i][idx] >= 0 && ptd.m_idata[IntIdx::workgroup][idx] > 0 && !ptd.m_idata[IntIdx::withdrawn][idx] && diff --git a/src/InteractionModWorkNborhood.H b/src/InteractionModWorkNborhood.H index 4f0a04d..425f202 100644 --- a/src/InteractionModWorkNborhood.H +++ b/src/InteractionModWorkNborhood.H @@ -52,7 +52,7 @@ template struct WorkNborhoodCandidate { AMREX_GPU_HOST_DEVICE bool operator() (const int idx, const PTDType& ptd) const noexcept { - return !isHospitalized(idx, ptd) && !ptd.m_idata[IntIdx::withdrawn][idx] && ptd.m_idata[IntIdx::random_travel][idx] < 0; + return !inHospital(idx, ptd) && !ptd.m_idata[IntIdx::withdrawn][idx] && ptd.m_idata[IntIdx::random_travel][idx] < 0; } };