Skip to content

Commit

Permalink
Add availability filter (#374)
Browse files Browse the repository at this point in the history 
* fix initialization of new variables

* add minAvailability and maxAvailability percentiles for timed human interventions

* add availability filter feature to timed cumulative and continuous deployments

* add test scenario for heterogeneity and availability based interventions

* cleanup

* improve error message for availability based interventions
  • Loading branch information
@acavelan
acavelan authored Aug 9, 2023
1 parent 4e6fabc commit 72baa9f
Show file tree
Hide file tree
Showing 10 changed files with 47,939 additions and 43 deletions.
10 changes: 10 additions & 0 deletions model/Host/Human.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -119,6 +119,11 @@ Human::Human(SimTime dateOfBirth) :
double iiFactor = infIncidence->getAvailabilityFactor(rng, 1.0);
perHostTransmission.initialise (rng, het.availabilityFactor * iiFactor);
clinicalModel = Clinical::ClinicalModel::createClinicalModel(het.treatmentSeekingFactor);

// Compute base availability
avail = 0.0;
for(size_t i = 0; i < Transmission::PerHostAnophParams::numSpecies(); ++i)
avail += perHostTransmission.anophEntoAvailability[i];
}

void Human::addToCohort(ComponentId id, SimTime duration )
Expand Down Expand Up @@ -173,6 +178,11 @@ void Human::updateCohortSet()
}
}

double Human::getAvailability() const
{
return avail;
}

void Human::checkpoint(istream &stream)
{
perHostTransmission & stream;
Expand Down
6 changes: 6 additions & 0 deletions model/Host/Human.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -91,6 +91,9 @@ class Human {
/** Remove host from expired cohorts */
void updateCohortSet();

/** Returns the base of availability of this host **/
double getAvailability() const;

/** Checkpoint (read) */
void checkpoint(istream &stream);

Expand Down Expand Up @@ -131,6 +134,9 @@ class Human {
/** The state of he human. A human cannot be revived. */
bool dead = false;

/** Base availability **/
double avail = 1.0;

/** This lists sub-populations of which the human is a member together with
* expiry time.
*
Expand Down
2 changes: 1 addition & 1 deletion model/Transmission/PerHost.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -33,6 +33,7 @@ using namespace OM::util;

AgeGroupInterpolator PerHost::relAvailAge;
vector<PerHostAnophParams> PerHostAnophParams::params;
vector<double> PerHostAnophParams::entoAvailabilityPercentiles;

PerHostAnophParams::PerHostAnophParams (const scnXml::Mosq& mosq) {
const string &distr = mosq.getAvailability().getDistr();
Expand Down Expand Up @@ -95,7 +96,6 @@ void PerHost::deployComponent( LocalRng& rng, const HumanVectorInterventionCompo
activeComponents.push_back( params.makeHumanPart(rng) );
}


// Note: in the case an intervention is not present, we can use the approximation
// of Weibull decay over the time span now - sim::never()
// (easily large enough for conceivable Weibull params that the value is 0.0 when
Expand Down
59 changes: 51 additions & 8 deletions model/Transmission/PerHost.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,6 +28,8 @@
#include "util/DecayFunction.h"
#include "util/checkpoint_containers.h"

#include "util/random.h"

namespace OM {
namespace Transmission {

Expand Down Expand Up @@ -71,6 +73,45 @@ class PerHostAnophParams {
params[species].entoAvailability->scaleMean( entoAvailability );
}

inline static void calcAvailabilityPercentiles(){
int nSamples = 100000;
vector<double> samples(nSamples, 0.0);

// Using a different seed to make sure old results will not change
// Does it really matter?
util::MasterRng seed(0, 0);
seed.seed(0, 0);

util::LocalRng rng(seed);

// Sample
for(int i=0; i<nSamples; i++)
{
double avail = 0.0;
for(size_t s = 0; s < Transmission::PerHostAnophParams::numSpecies(); ++s)
avail += Transmission::PerHostAnophParams::get(s).entoAvailability->sample(rng);
samples[i] = avail;
}

// Sort samples
sort(samples.begin(), samples.end());

entoAvailabilityPercentiles.resize(100);

// Calc percetiles threshold values
for(int i=0; i<100; i++)
entoAvailabilityPercentiles[i] = samples[int(i*nSamples/100)];
entoAvailabilityPercentiles[0] = 0.0;
}

inline static double getEntoAvailabilityPercentile(int p){
if(p >= 0 && p < 100)
return entoAvailabilityPercentiles[p];
else if(p == 100)
return std::numeric_limits<double>::max();
else throw util::xml_scenario_error("Availability percentiles in intervention specifications can only be in [0 and 100]");
}

/** @brief Probabilities of finding a host and surviving a feeding cycle
*
* These parameters describe the mean and heterogeneity of α_i, P_B_i,
Expand All @@ -95,6 +136,8 @@ class PerHostAnophParams {
PerHostAnophParams (const scnXml::Mosq& mosq);

static vector<PerHostAnophParams> params;
static vector<double> entoAvailabilityPercentiles;

};

/**
Expand Down Expand Up @@ -302,14 +345,6 @@ class PerHost
// entoAvailability param stored in HostMosquitoInteraction.
double relativeAvailabilityHet;

private:
void checkpointIntervs( ostream& stream );
void checkpointIntervs( istream& stream );

vector<unique_ptr<PerHostInterventionData>> activeComponents;

static AgeGroupInterpolator relAvailAge;

/** Species availability rate of human to mosquitoes, including hetergeneity factor
* and base rate, but excluding age and intervention factors. */
vector<double> anophEntoAvailability;
Expand All @@ -322,6 +357,14 @@ class PerHost
* resting without dying, after biting the human (P_C_i * P_D_i) in the
* absense of interventions. */
vector<double> anophProbMosqResting;

private:
void checkpointIntervs( ostream& stream );
void checkpointIntervs( istream& stream );

vector<unique_ptr<PerHostInterventionData>> activeComponents;

static AgeGroupInterpolator relAvailAge;
};

}
Expand Down
69 changes: 49 additions & 20 deletions model/interventions/Deployments.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -198,23 +198,35 @@ class TimedHumanDeployment : public TimedDeployment, protected HumanDeploymentBa
TimedDeployment( date ),
HumanDeploymentBase( mass, intervention, subPop, complement ),
minAge( sim::fromYearsN( mass.getMinAge() ) ),
maxAge( sim::future() )
maxAge( sim::future() ),
minAvailability( 0 ),
maxAvailability( 100 )
{
if( mass.getMaxAge().present() )
maxAge = sim::fromYearsN( mass.getMaxAge().get() );

if( minAge < sim::zero() || maxAge < minAge ){
throw util::xml_scenario_error("timed intervention must have 0 <= minAge <= maxAge");
}
if( minAge < sim::zero() || maxAge < minAge )
throw util::xml_scenario_error("timed intervention must have 0 <= minAvailability <= maxAvailability <= 100");

if( mass.getMaxAvailability().present() )
maxAvailability = mass.getMaxAvailability().get();

if( mass.getMinAvailability().present() )
minAvailability = mass.getMinAvailability().get();

if( minAvailability < sim::zero() || maxAvailability < minAvailability )
throw util::xml_scenario_error("timed intervention must have 0 <= minAvailability <= maxAvailability <= 100");
}

virtual void deploy (vector<Host::Human> &population, Transmission::TransmissionModel& transmission) {
for(Human& human : population) {
SimTime age = human.age(sim::now());
if( age >= minAge && age < maxAge ){
if( subPop == ComponentId::wholePop() || (human.isInSubPop( subPop ) != complement) ){
if( human.rng.bernoulli( coverage ) ){
deployToHuman( human, mon::Deploy::TIMED );
if( human.getAvailability() >= Transmission::PerHostAnophParams::getEntoAvailabilityPercentile(minAvailability) && human.getAvailability() <= Transmission::PerHostAnophParams::getEntoAvailabilityPercentile(maxAvailability) ) {
if( subPop == ComponentId::wholePop() || (human.isInSubPop( subPop ) != complement) ){
if( human.rng.bernoulli( coverage ) ){
deployToHuman( human, mon::Deploy::TIMED );
}
}
}
}
Expand All @@ -233,6 +245,7 @@ class TimedHumanDeployment : public TimedDeployment, protected HumanDeploymentBa
protected:
// restrictions on deployment
SimTime minAge, maxAge;
SimTime minAvailability, maxAvailability;
};

/// Timed deployment of human-specific interventions in cumulative mode
Expand Down Expand Up @@ -263,10 +276,12 @@ class TimedCumulativeHumanDeployment : public TimedHumanDeployment {
for(Host::Human &human : population) {
SimTime age = human.age(sim::now());
if( age >= minAge && age < maxAge ){
if( subPop == ComponentId::wholePop() || (human.isInSubPop( subPop ) != complement) ){
total+=1;
if( !human.isInSubPop(cumCovInd) )
unprotected.push_back( &human );
if( human.getAvailability() >= Transmission::PerHostAnophParams::getEntoAvailabilityPercentile(minAvailability) && human.getAvailability() <= Transmission::PerHostAnophParams::getEntoAvailabilityPercentile(maxAvailability) ) {
if( subPop == ComponentId::wholePop() || (human.isInSubPop( subPop ) != complement) ){
total+=1;
if( !human.isInSubPop(cumCovInd) )
unprotected.push_back( &human );
}
}
}
}
Expand Down Expand Up @@ -515,7 +530,9 @@ class ContinuousHumanDeployment : protected HumanDeploymentBase {
ComponentId subPop, bool complement ) :
HumanDeploymentBase( elt, intervention, subPop, complement ),
begin( begin ), end( end ),
deployAge( sim::fromYearsN( elt.getTargetAgeYrs() ) )
deployAge( sim::fromYearsN( elt.getTargetAgeYrs() ) ),
minAvailability( 0 ),
maxAvailability( 100 )
{
if( begin < sim::startDate() || end < begin ){
throw util::xml_scenario_error("continuous intervention must have startDate <= begin <= end");
Expand All @@ -533,6 +550,14 @@ class ContinuousHumanDeployment : protected HumanDeploymentBase {
msg << sim::inYears(sim::maxHumanAge());
throw util::xml_scenario_error( msg.str() );
}
if( elt.getMaxAvailability().present() )
maxAvailability = elt.getMaxAvailability().get();

if( elt.getMinAvailability().present() )
minAvailability = elt.getMinAvailability().get();

if( minAvailability < sim::zero() || maxAvailability < minAvailability )
throw util::xml_scenario_error("timed intervention must have 0 <= minAvailability <= maxAvailability <= 100");
}

/** Apply filters and potentially deploy.
Expand All @@ -546,14 +571,17 @@ class ContinuousHumanDeployment : protected HumanDeploymentBase {
// human for now because remaining ones happen in the future
return false;
}else if( deployAge == age ){
auto now = sim::intervDate();
if( begin <= now && now < end &&
( subPop == ComponentId::wholePop() ||
(human.isInSubPop( subPop ) != complement)
) &&
human.rng.uniform_01() < coverage ) // RNG call should be last test
if( human.getAvailability() >= Transmission::PerHostAnophParams::getEntoAvailabilityPercentile(minAvailability) && human.getAvailability() <= Transmission::PerHostAnophParams::getEntoAvailabilityPercentile(maxAvailability) )
{
deployToHuman( human, mon::Deploy::CTS );
auto now = sim::intervDate();
if( begin <= now && now < end &&
( subPop == ComponentId::wholePop() ||
(human.isInSubPop( subPop ) != complement)
) &&
human.rng.uniform_01() < coverage ) // RNG call should be last test
{
deployToHuman( human, mon::Deploy::CTS );
}
}
}//else: for some reason, a deployment age was missed; ignore it
return true;
Expand All @@ -574,7 +602,8 @@ class ContinuousHumanDeployment : protected HumanDeploymentBase {
protected:
SimTime begin, end; // first time step active and first time step no-longer active
SimTime deployAge = sim::never();

SimTime minAvailability, maxAvailability;

friend ByDeployTime;
};

Expand Down
17 changes: 3 additions & 14 deletions model/main.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -224,20 +224,9 @@ int main(int argc, char* argv[])
Continuous.init(scenario->getMonitoring(), false);
population->createInitialHumans();
transmission->init2(population->humans);

// ofstream humanFile ("humans.txt");
// if (humanFile.is_open())
// {
// for(size_t i=0; i<population->getHumans()[0].perHostTransmission.speciesData.size(); i++)
// {
// for(const Host::Human &human : population->getHumans())
// {
// humanFile << human.perHostTransmission.speciesData[i].getEntoAvailability() << " ";
// }
// humanFile << endl;
// }
// humanFile.close();
// }

/** Calculate ento availability percentiles **/
Transmission::PerHostAnophParams::calcAvailabilityPercentiles();

/** Warm-up phase:
* Run the simulation using the equilibrium inoculation rates over one
Expand Down
23 changes: 23 additions & 0 deletions schema/interventions.xsd
View file
Original file line number Diff line number Diff line change
Expand Up @@ -968,6 +968,29 @@ Licence: GNU General Public Licence version 2 or later (see COPYING) -->
<xs:appinfo>name:Vaccine max cumulative doses;units:inoculations;min:0;</xs:appinfo>
</xs:annotation>
</xs:attribute>
<xs:attribute name="minAvailability" type="xs:int" use="optional">
<xs:annotation>
<xs:documentation>
Minimum ento availability percentile.

This option is meant to be used with heterogeneity of availability, which can be specified in the entomology section. Without heterogenity (default), all hosts have the same availability and this option will have no effect.

The percentile must be an integer value between 0 and 100. Percentile 99th represents individuals who are more available than 99% of the population. Percentile 0 represents the least available individuals. 100 is equivalent for infinity.
</xs:documentation>
<xs:appinfo>units:percent;min:0;name:Minimum ento availability;</xs:appinfo>
</xs:annotation>
</xs:attribute>
<xs:attribute name="maxAvailability" type="xs:int" use="optional">
<xs:annotation>
<xs:documentation>
Maximum ento availability percentile.

This option is meant to be used with heterogeneity of availability, which can be specified in the entomology section. Without heterogenity (default), all hosts have the same availability and this option will have no effect.

The percentile must be an integer value between 0 and 100. Percentile 99th represents individuals who are more available than 99% of the population. Percentile 0 represents the least available individuals. 100 is equivalent for infinity.</xs:documentation>
<xs:appinfo>units:percent;min:0;name:Maximum ento availability;</xs:appinfo>
</xs:annotation>
</xs:attribute>
</xs:complexType>
<xs:complexType name="MassDeployment">
<xs:complexContent>
Expand Down
1 change: 1 addition & 0 deletions test/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -54,6 +54,7 @@ set (OM_BOXTEST_NAMES
TBV
DecisionTree5DayDielmo
HetVecDailyEIR
AvailabilityHeterogeneity
)
# tests with broken checkpointing:
set (OM_BOXTEST_NC_NAMES)
Expand Down
Loading

0 comments on commit 72baa9f

Please sign in to comment.