Skip to content

Commit

Permalink
Fix Disease Status (#88)
Browse files Browse the repository at this point in the history
* 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
  • Loading branch information
debog authored Nov 12, 2024
1 parent 5d83e0a commit 12526c7
Show file tree
Hide file tree
Showing 9 changed files with 46 additions and 40 deletions.
16 changes: 11 additions & 5 deletions src/AgentContainer.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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];
Expand Down Expand Up @@ -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];
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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];
Expand Down Expand Up @@ -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];
Expand Down Expand Up @@ -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;
}
}
Expand Down
4 changes: 2 additions & 2 deletions src/AgentDefinitions.H
Original file line number Diff line number Diff line change
Expand Up @@ -208,8 +208,8 @@ bool notSusceptible ( const int a_idx, /*!< Agent index */
/*! \brief Is an agent hospitalized? */
template <typename PTDType>
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) );
Expand Down
46 changes: 23 additions & 23 deletions src/DiseaseStatus.H
Original file line number Diff line number Diff line change
Expand Up @@ -117,18 +117,18 @@ void DiseaseStatus<AC,ACT,ACTD,A>::updateAgents(AC& a_agents, /*!< Agent contain
int i_RT = IntIdx::nattribs;
int r_RT = RealIdx::nattribs;

Gpu::DeviceVector<int> 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<int> 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();

Expand Down Expand Up @@ -175,7 +175,7 @@ void DiseaseStatus<AC,ACT,ACTD,A>::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;
}
Expand All @@ -186,21 +186,21 @@ void DiseaseStatus<AC,ACT,ACTD,A>::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);
Expand All @@ -216,8 +216,8 @@ void DiseaseStatus<AC,ACT,ACTD,A>::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);
}
});
Expand All @@ -226,21 +226,21 @@ void DiseaseStatus<AC,ACT,ACTD,A>::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,
Expand Down
10 changes: 5 additions & 5 deletions src/HospitalModel.H
Original file line number Diff line number Diff line change
Expand Up @@ -85,13 +85,13 @@ void HospitalModel<PCType, PTDType, PType>::treatAgents(PCType& a_agents, /*!< A
int r_RT = RealIdx::nattribs;

GpuArray<int*,ExaEpi::max_num_diseases> status_ptrs, symptomatic_ptrs;
GpuArray<ParticleReal*,ExaEpi::max_num_diseases> counter_ptrs, timer_ptrs, latent_per_ptrs;
GpuArray<ParticleReal*,ExaEpi::max_num_diseases> 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<int> is_alive;
Expand Down Expand Up @@ -132,11 +132,11 @@ void HospitalModel<PCType, PTDType, PType>::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;
}
Expand Down Expand Up @@ -213,7 +213,7 @@ void HospitalModel<PCType, PTDType, PType>::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) {

Expand Down
2 changes: 1 addition & 1 deletion src/InteractionModHome.H
Original file line number Diff line number Diff line change
Expand Up @@ -59,7 +59,7 @@ template <typename PTDType>
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;
}
};

Expand Down
2 changes: 1 addition & 1 deletion src/InteractionModHomeNborhood.H
Original file line number Diff line number Diff line change
Expand Up @@ -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];
}
};

Expand Down
2 changes: 1 addition & 1 deletion src/InteractionModSchool.H
Original file line number Diff line number Diff line change
Expand Up @@ -69,7 +69,7 @@ template <typename PTDType>
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 &&
Expand Down
2 changes: 1 addition & 1 deletion src/InteractionModWork.H
Original file line number Diff line number Diff line change
Expand Up @@ -57,7 +57,7 @@ template <typename PTDType>
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] &&
Expand Down
2 changes: 1 addition & 1 deletion src/InteractionModWorkNborhood.H
Original file line number Diff line number Diff line change
Expand Up @@ -52,7 +52,7 @@ template <typename PTDType>
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;
}
};

Expand Down

0 comments on commit 12526c7

Please sign in to comment.