Add elimination model: track imported, introduced and indigenous infe…

…ctions and corresponding new outputs (#384)

* add InfectionOrigin to Infection and WithinHost classes

* simplify probTransmissionToMosquito()

* split m_y_lag into m_y_lag_i and am_y_lag_l

* simplify probTransmissionToMosauito

* split EIR survey calculation from updateKappa()

* merge probTransmissionToMosquito and probTransGenotype

* remove old code

* calculate probTransmission from imported and local infections separetly

* calculate P_dif from imported and local infections separately

* compute P_dif_i, P_dif_l, O_v_i, O_v_l, S_v_i, S_v_l

* add EIR_i, EIR_l

* compute nNewInfs using two poissons

* add introduced infections in WithinHost models

* complete imported -> introduced -> indigenous cycle

* cleanup comments

* update expected output of tests using imported infections due to additional rng call

* remove warning for unused variable

* add Introduced and Indigenous EIR outputs

* add new output measures related to elimination

This commit adds "Imported", "Introduced" and "Indigenous" for the output measures:
nInfect
nPatent
totalInfs
totalPatentInf
SimulatedEIR (only Introduced and Indigenous)

The corresponding output numbers are calculated as follows:
Imported: +1000
Introduced: +2000
Indigenous: +3000

For example "nPatent" is as follows:
nPatent: 3
nPatent_Imported: 1003
nPatent_Introduced: 2003
nPatent_Indigenous: 3003

* add elimination support for nUncomp

model/Clinical/Episode.cpp

@@ -22,6 +22,7 @@
 #include "Clinical/Episode.h"
 #include "Clinical/ClinicalModel.h"
 #include "Host/Human.h"
+#include "Host/WithinHost/WHInterface.h"
 
 namespace OM {
 namespace Clinical {
@@ -40,6 +41,8 @@ void Episode::flush() {
 void Episode::update (const Host::Human& human, Episode::State newState)
 {
     if( time + ClinicalModel::hsMemory() <= sim::ts0() ){
+        infectionType = human.withinHostModel->getInfectionType();
+
         report ();
 
         time = sim::ts0();
@@ -66,6 +69,12 @@ void Episode::report () {
             }
         } else { // UC or UC2
             mon::reportMSACI( mon::MHE_UNCOMPLICATED_EPISODES, surveyPeriod, ageGroup, cohortSet, 1 );
+            if(infectionType == WithinHost::InfectionOrigin::Indigenous)
+                mon::reportMSACI( mon::MHE_UNCOMPLICATED_EPISODES_INDIGENOUS, surveyPeriod, ageGroup, cohortSet, 1 );
+            else if(infectionType == WithinHost::InfectionOrigin::Introduced)
+                mon::reportMSACI( mon::MHE_UNCOMPLICATED_EPISODES_INTRODUCED, surveyPeriod, ageGroup, cohortSet, 1 );
+            else
+                mon::reportMSACI( mon::MHE_UNCOMPLICATED_EPISODES_IMPORTED, surveyPeriod, ageGroup, cohortSet, 1 );
         }
 
         // Report outcomes of malarial fevers

model/Clinical/Episode.h

@@ -24,6 +24,7 @@
 
 #include "Global.h"
 #include "Host/WithinHost/Pathogenesis/State.h"
+#include "Host/WithinHost/Infection/Infection.h"
 #include "mon/AgeGroup.h"
 #include "mon/info.h"
 #include <ostream>
@@ -112,6 +113,8 @@ class Episode{
     /// Descriptor of state, containing reporting info. Not all information will
     /// be reported (e.g. indirect deaths are reported independantly).
     Episode::State state;
+
+    WithinHost::InfectionOrigin infectionType = WithinHost::InfectionOrigin::Indigenous;
 
 private:
     /** Report a clinical episode.

model/Host/Human.cpp

@@ -214,8 +214,6 @@ void Human::checkpoint(ostream &stream)
 }
 
 // -----  Non-static functions: per-time-step update  -----
-vector<double> EIR_per_genotype;        // cache (not thread safe)
-
 void summarize(Human &human, bool surveyOnlyNewEp) {
     if( surveyOnlyNewEp && human.clinicalModel->isExistingCase() ){
         // This modifies the denominator to treat the health-system-memory
@@ -257,14 +255,30 @@ void update(Human &human, Transmission::TransmissionModel& transmission)
     double age1 = sim::inYears(human.age(sim::ts1()));
 
     // age1 used only in PerHost::relativeAvailabilityAge(); difference to age0 should be minor
-    double EIR = transmission.getEIR( human, age0, age1, EIR_per_genotype );
+    vector<double> EIR_per_genotype_i, EIR_per_genotype_l;
+    transmission.getEIR(human, age0, age1, EIR_per_genotype_i, EIR_per_genotype_l);
+
+    double EIR_i = util::vectors::sum(EIR_per_genotype_i);
+    double EIR_l = util::vectors::sum(EIR_per_genotype_l);
+    double EIR = EIR_i + EIR_l;
 
-    int nNewInfs = human.infIncidence->numNewInfections( human, EIR );
+    double x = 0.0;
 
+    if(EIR > 0)
+        x = EIR_i / EIR;
+
+    double expectedNNewInfs = human.infIncidence->expectedNumNewInfections( human, EIR);
+    double expectedNNewInfs_i = x * expectedNNewInfs;
+    double expectedNNewInfs_l = (1.0-x) * expectedNNewInfs;
+
+    // One more rng call if imported infection
+    int nNewInfs_l = human.infIncidence->numNewInfections( human, expectedNNewInfs_l);
+    int nNewInfs_i = human.infIncidence->numNewInfections( human, expectedNNewInfs_i);
+
     // age1 used when medicating drugs (small effect) and in immunity model (which was parameterised for it)
-    human.withinHostModel->update(human, human.rng, nNewInfs, EIR_per_genotype, age1);
-    human.infIncidence->reportNumNewInfections(human, nNewInfs);
-    
+    human.withinHostModel->update(human, human.rng, nNewInfs_i, nNewInfs_l, EIR_per_genotype_i, EIR_per_genotype_l, age1);
+    human.infIncidence->reportNumNewInfections(human, nNewInfs_i+nNewInfs_l);
+
     // age1 used to get case fatality and sequelae probabilities, determine pathogenesis
     human.clinicalModel->update( human, age1, age0 == sim::zero() );
     human.clinicalModel->updateInfantDeaths( age0 );

model/Host/ImportedInfections.h

@@ -28,6 +28,7 @@
 #include "util/UnitParse.h"
 #include "Host/Human.h"
 #include "Host/WithinHost/WHInterface.h"
+#include "Host/WithinHost/Infection/Infection.h"
 
 namespace OM {
     class Population;
@@ -115,7 +116,7 @@ namespace Host {
             if( rateNow > 0.0 ){
                 for(Human& human : population){
                     if(human.rng.bernoulli( rateNow )){
-                        human.withinHostModel->importInfection(human.rng);
+                        human.withinHostModel->importInfection(human.rng, WithinHost::InfectionOrigin::Imported);
                     }
                 }
             }

model/Host/InfectionIncidenceModel.cpp

@@ -214,30 +214,36 @@ double InfectionIncidenceModel::susceptibility () {
   }
 }
 
-int InfectionIncidenceModel::numNewInfections (Human& human, double effectiveEIR) {
-  double expectedNumInfections = getModelExpectedInfections (human.rng, effectiveEIR, human.perHostTransmission);
-  // error check (should be OK if kappa is checked, for nonVector model):
-  if( !(std::isfinite)(effectiveEIR) ){
-    ostringstream out;
-    out << "effectiveEIR is not finite: " << effectiveEIR << endl;
-    throw TRACED_EXCEPTION (out.str(), util::Error::EffectiveEIR);
-  }
-
-  // Only to be consistent with old simulation runs when set to false
-  // Setting this option to true will only affect reporting
-  if(opt_vaccine_genotype == false)
-      //Introduce the effect of vaccination. Note that this does not affect cumEIR.
-      expectedNumInfections *= human.vaccine.getFactor( interventions::Vaccine::PEV );
-
-  //Update pre-erythrocytic immunity
-  m_cumulativeEIRa+=effectiveEIR;
+double InfectionIncidenceModel::expectedNumNewInfections(OM::Host::Human& human, double effectiveEIR)
+{
+    double expectedNumInfections = getModelExpectedInfections (human.rng, effectiveEIR, human.perHostTransmission);
+    // error check (should be OK if kappa is checked, for nonVector model):
+    if( !(std::isfinite)(effectiveEIR) ){
+        ostringstream out;
+        out << "effectiveEIR is not finite: " << effectiveEIR << endl;
+        throw TRACED_EXCEPTION (out.str(), util::Error::EffectiveEIR);
+    }
 
-  m_pInfected = 1.0 - exp(-expectedNumInfections) * (1.0-m_pInfected);
-  if (m_pInfected < 0.0)
-    m_pInfected = 0.0;
-  else if (m_pInfected > 1.0)
-    m_pInfected = 1.0;
+    // Only to be consistent with old simulation runs when set to false
+    // Setting this option to true will only affect reporting
+    if(opt_vaccine_genotype == false)
+        //Introduce the effect of vaccination. Note that this does not affect cumEIR.
+        expectedNumInfections *= human.vaccine.getFactor( interventions::Vaccine::PEV );
 
+    //Update pre-erythrocytic immunity
+    m_cumulativeEIRa+=effectiveEIR;
+
+    m_pInfected = 1.0 - exp(-expectedNumInfections) * (1.0-m_pInfected);
+    if (m_pInfected < 0.0)
+        m_pInfected = 0.0;
+    else if (m_pInfected > 1.0)
+        m_pInfected = 1.0;
+
+    return expectedNumInfections;
+}
+
+int InfectionIncidenceModel::numNewInfections (Human& human, double expectedNumInfections)
+{ 
   if (expectedNumInfections > 0.0000001){
     int n = human.rng.poisson(expectedNumInfections);
     if( n > WithinHost::WHInterface::MAX_INFECTIONS ){

model/Host/InfectionIncidenceModel.h

@@ -91,7 +91,7 @@ class InfectionIncidenceModel
   /// Output _pinfected to the summary
   void summarize (const Host::Human& human);
 
-  /** Calculate the number of new infections to introduce.
+    /** Calculate the expected number of new infections to introduce.
    * 
    * Firstly converts the EIR into an expected number of infections (what
    * getExpectedNumberOfInfections() used to do):
@@ -106,7 +106,15 @@ class InfectionIncidenceModel
    * @param effectiveEIR    EIR adjusted for availability (age, het,
    *        outside transmission due to hospitalisation)
    */
-  int numNewInfections(OM::Host::Human& human, double effectiveEIR);
+  double expectedNumNewInfections(OM::Host::Human& human, double effectiveEIR);
+
+  /** Calculate the number of new infections to introduce.
+   * 
+   * @param human   Reference to human
+   * @param expectedNumNewInfections   Expected number of infections
+   * to be introduced
+   */
+  int numNewInfections(OM::Host::Human& human, double expectedNumNewInfections);
 
   /** Report the number of new infections
    *