Reduce coverage

src/Model.cpp

@@ -29,7 +29,7 @@ void Model::runScenarios(ScenariosList &scenarios, Population &popln,
                          double timestep, int index, int outputEndgame,
                          int outputEndgameDate, int reduceImpViaXml,
                          std::string randParamsfile, std::string RandomSeedFile,
-                         std::string opDir) {
+                         std::string RandomCovPropFile, std::string opDir) {
 
   std::cout << std::endl
             << "Index " << index << " running " << scenarios.getName()
@@ -59,6 +59,11 @@ void Model::runScenarios(ScenariosList &scenarios, Population &popln,
   // used for each set of parameters
   readSeedsFromFile(seeds, unsigned(replicates), RandomSeedFile);
 
+  std::vector<double> cov_props;
+  // we read in the entries of the random cov_props file. a different value will
+  // be used for each set of parameters
+  readCovPropFromFile(cov_props, unsigned(replicates), RandomCovPropFile);
+
   for (int rep = 0; rep < replicates; rep++) {
     // Read the seed from the seeds vector if it has been generated
     // othewise we set a random seed
@@ -70,6 +75,14 @@ void Model::runScenarios(ScenariosList &scenarios, Population &popln,
       rseed = std::chrono::system_clock::now().time_since_epoch().count();
       stats.set_seed(rseed);
     }
+    // Read the value to multiply the MDA's by if this has been supplied
+    // othewise we set this to 1
+    double cov_prop;
+    if (cov_props.size() > 0) {
+      cov_prop = cov_props[rep];
+    } else {
+      cov_prop = 1.0;
+    }
     getRandomParametersMultiplePerLine(rep + 1, k_vals, v_to_h_vals, aImp_vals,
                                        wPropMDA, unsigned(replicates),
                                        randParamsfile);
@@ -135,7 +148,7 @@ void Model::runScenarios(ScenariosList &scenarios, Population &popln,
       for (int y = 0; y < sc.getNumMonthsToSave(); y++)
         evolveAndSave(y, popln, vectors, worms, sc, currentOutput, rep, k_vals,
                       v_to_h_vals, popln.getUpdateParams(), outputEndgame,
-                      outputEndgameDate, reduceImpViaXml, opDir);
+                      outputEndgameDate, reduceImpViaXml, opDir, cov_prop);
 
       // done for this scenario, save the prevalence values for this replicate
       if (!_DEBUG)
@@ -202,7 +215,8 @@ void Model::evolveAndSave(int y, Population &popln, Vector &vectors,
                           int rep, std::vector<double> &k_vals,
                           std::vector<double> &v_to_h_vals, int updateParams,
                           int outputEndgame, int outputEndgameDate,
-                          int reduceImpViaXml, std::string opDir) {
+                          int reduceImpViaXml, std::string opDir,
+                          double cov_prop) {
 
   // advance to the next target month
   std::string folderName = opDir;
@@ -408,11 +422,15 @@ void Model::evolveAndSave(int y, Population &popln, Vector &vectors,
       }
     }
 
-    // std::cout << applyMDA->getMonth() << std::endl;
+    // std::cout << applyMDA->getMonth() << std::xendl;
     if (applyMDA) {
       MDAType = popln.returnMDAType(applyMDA);
       minAge = popln.returnMinAgeMDA(applyMDA);
-      double cov = applyMDA->getCoverage();
+      double coverage_multiplier =
+          multiplierForCoverage(t, cov_prop, popln.removeCoverageReduction,
+                                popln.removeCoverageReductionTime,
+                                popln.graduallyRemoveCoverageReduction);
+      double cov = applyMDA->getCoverage() * coverage_multiplier;
       double rho = applyMDA->getCompliance();
       // if we this is the first MDA, then we need to initialise the Probability
       // of treatment for each person
@@ -423,7 +441,7 @@ void Model::evolveAndSave(int y, Population &popln, Vector &vectors,
       }
       // if the MDA parameters have changed then we need to update the
       // probability of treatment for each person
-      if ((popln.prevCov != applyMDA->getCoverage()) ||
+      if ((popln.prevCov != applyMDA->getCoverage() * coverage_multiplier) ||
           (popln.prevRho != applyMDA->getCompliance())) {
         popln.checkForZeroPTreat(popln.prevCov, popln.prevRho);
         popln.editPTreat(cov, rho);
@@ -603,6 +621,26 @@ void Model::getRandomParameters(int index, std::vector<double> &k_vals,
   }
 }
 
+double Model::multiplierForCoverage(int t, double cov_prop,
+                                    int removeCoverageReduction,
+                                    int removeCoverageReductionTime,
+                                    int graduallyRemoveCoverageReduction) {
+
+  if (graduallyRemoveCoverageReduction == 1) {
+    if (t < removeCoverageReductionTime)
+      return ((1 - cov_prop)) * (double(t) / removeCoverageReductionTime) +
+             cov_prop;
+    else
+      return 1;
+  } else if (removeCoverageReduction == 1) {
+    if (t < removeCoverageReductionTime)
+      return cov_prop;
+    else
+      return 1;
+  } else
+    return cov_prop;
+}
+
 void Model::ProcessLine(const std::string &line, std::vector<double> &k_vals,
                         std::vector<double> &v_to_h_vals,
                         std::vector<double> &aImp_vals,
@@ -704,6 +742,40 @@ void Model::readSeedsFromFile(std::vector<unsigned long int> &seeds,
   }
 }
 
+void Model::readCovPropFromFile(std::vector<double> &cov_props,
+                                unsigned replicates, std::string fname) {
+  // We retrieve the random MDA shrinkage value from the input covproporation
+  // file. The line on which the value is on will correspond to the set of
+  // parameters on the same line of the input parameters file. If there is no
+  // covproporation file input we set it to 1.
+  //
+  if (fname.length() > 0) {
+    std::ifstream infile(fname, std::ios_base::in);
+    double cov_prop;
+    if (!infile.is_open()) {
+      std::cout << "Error in getting shrinkage. Cannot read file " << fname
+                << std::endl;
+      exit(1);
+    }
+    std::string line;
+    while (getline(infile, line)) {
+      std::istringstream linestream(line);
+      linestream >> cov_prop;
+      cov_props.push_back(cov_prop);
+    }
+    // if we haven't input enough values based on how many runs we want to make
+    // then we will abort the run as at some point we will run out of values and
+    // the simulations will crash.
+    if (cov_props.size() < replicates) {
+      std::cout << "Error in Model::runScenarios. " << fname
+                << " is too short for " << replicates << " replicates"
+                << std::endl;
+      exit(1);
+    }
+    infile.close();
+  }
+}
+
 std::vector<std::string> Model::printSeedName() const {
 
   // outputs "seed" as a column title for fitting output file

src/Model.hpp

@@ -30,11 +30,16 @@ class Model {
                     double timestep, int index, int outputEndgame,
                     int outputEndgameDate, int reduceImpViaXml,
                     std::string randParamsfile, std::string RandomSeedFile,
-                    std::string opDir);
+                    std::string RandomCovPropFile, std::string opDir);
   bool
   shouldReduceImportationViaPrevalance(int t, int reduceImpViaXml,
                                        int switchImportationReducingMethodTime);
 
+  double multiplierForCoverage(int t, double cov_prop,
+                               int removeCoverageReduction,
+                               int removeCoverageReductionTime,
+                               int graduallyRemoveCoverageReduction);
+
 protected:
   void burnIn(Population &popln, Vector &vectors, const Worm &worms,
               Output &currentOutput, PrevalenceEvent *pe);
@@ -43,7 +48,7 @@ class Model {
                      std::vector<double> &k_vals,
                      std::vector<double> &v_to_h_vals, int updateParams,
                      int outputEndgame, int outputEndgameDate,
-                     int reduceImpViaXml, std::string opDir);
+                     int reduceImpViaXml, std::string opDir, double cov_prop);
   void getRandomParameters(int index, std::vector<double> &k_vals,
                            std::vector<double> &v_to_h_vals,
                            std::vector<double> &aImp_vals,
@@ -55,6 +60,8 @@ class Model {
       unsigned replicates, std::string fname);
   void readSeedsFromFile(std::vector<unsigned long int> &seeds,
                          unsigned replicates, std::string fname);
+  void readCovPropFromFile(std::vector<double> &cov_props, unsigned replicates,
+                           std::string fname);
   void ProcessLine(const std::string &line, std::vector<double> &k_vals,
                    std::vector<double> &v_to_h_vals,
                    std::vector<double> &aImp_vals,

src/Population.cpp

@@ -143,6 +143,12 @@ Population::Population(TiXmlElement *xmlParameters) {
       neverTreatChangeCount = neverTreatChangeCount + 1;
     } else if (name == "switchImportationReducingMethodTime") {
       switchImportationReducingMethodTime = value;
+    } else if (name == "removeCoverageReduction") {
+      removeCoverageReduction = value;
+    } else if (name == "removeCoverageReductionTime") {
+      removeCoverageReductionTime = value;
+    } else if (name == "graduallyRemoveCoverageReduction") {
+      graduallyRemoveCoverageReduction = value;
     } else
       std::cout << "Unknown parameter " << name << " in Host parameter list."
                 << std::endl;

src/Population.hpp

@@ -90,7 +90,6 @@ class Population {
   void neverTreatToOriginal();
   void ICTestToOriginal();
   void updateKVal(double k_val);
-
   double getICTestChangeTime();
   void updateBedNetCoverage(BedNetEvent *bn);
   void updateImportationRate(double factor);
@@ -149,6 +148,16 @@ class Population {
       99999; // initialise the time to switch importation reduction to a large
              // arbitrary value, so that if we don't specify a different value
              // in the XML scenario file, then it will never be passed
+  int removeCoverageReduction =
+      0; // set to 0 meaning that if not read from file then we will not remove
+         // the reduction to MDA coverage
+  int removeCoverageReductionTime =
+      99999; // initialise the time to remove coverage reduction to a large
+             // arbitrary value, so that if we don't specify a different value
+             // in the XML scenario file, then it will never be passed
+  int graduallyRemoveCoverageReduction =
+      0; // set to 0 meaning that if not read from file then we will not
+         // gradually remove the reduction to MDA coverage
   std::vector<std::string> sensSpecChangeName;
   std::vector<std::string> neverTreatChangeName;
 

src/main.cpp

@@ -66,6 +66,7 @@ int main(int argc, char **argv) {
   std::string scenariosFile("");
   std::string opDir("");
   std::string RandomSeedFile("");
+  std::string CoverageReductionFile("");
 
   // initialize random seed value, whether the endgame output will be done
   // and whether the reduction in importation rate should be done via the
@@ -105,6 +106,8 @@ int main(int argc, char **argv) {
       opDir = argv[i + 1];
     else if (!strcmp(argv[i], "-g"))
       RandomSeedFile = argv[i + 1];
+    else if (!strcmp(argv[i], "-c"))
+      CoverageReductionFile = argv[i + 1];
     else if (!strcmp(argv[i], "-e"))
       outputEndgame = atoi(argv[i + 1]);
     else if (!strcmp(argv[i], "-D"))
@@ -191,7 +194,8 @@ int main(int argc, char **argv) {
   Model model;
   model.runScenarios(Scenarios, hostPopulation, vectors, worms, replicates, dt,
                      index, outputEndgame, outputEndgameDate, reduceImpViaXml,
-                     randParamsfile, RandomSeedFile, opDir);
+                     randParamsfile, RandomSeedFile, CoverageReductionFile,
+                     opDir);
 
   gettimeofday(&tv2, NULL);
   double timesofar = (double)(tv2.tv_usec - tv1.tv_usec) / 1000000.0 +

tests/test_model.cpp

@@ -37,4 +37,77 @@ TEST_CASE("Model", "[classic]") {
                      (t % 12 == 0)));
     }
   }
+
+  SECTION("Model::multiplierForCoverage") {
+
+    Model model;
+
+    SECTION("coverage multiplier should be unchanged when we don't remove "
+            "coverage reduction ") {
+      int removeCoverageReduction = 0;
+      int removeCoverageReductionTime = 276;
+      int graduallyRemoveCoverageReduction = 0;
+      double cov_prop = 0.6;
+      REQUIRE(model.multiplierForCoverage(0, cov_prop, removeCoverageReduction,
+                                          removeCoverageReductionTime,
+                                          graduallyRemoveCoverageReduction) ==
+              0.6);
+    }
+
+    SECTION("coverage multiplier should be 1 after removeCoverageReductionTime "
+            "when we remove "
+            "coverage reduction ") {
+      int removeCoverageReduction = 1;
+      int removeCoverageReductionTime = 276;
+      int graduallyRemoveCoverageReduction = 0;
+      double cov_prop = 0.6;
+      int t = 277;
+      REQUIRE(model.multiplierForCoverage(t, cov_prop, removeCoverageReduction,
+                                          removeCoverageReductionTime,
+                                          graduallyRemoveCoverageReduction) ==
+              1);
+      graduallyRemoveCoverageReduction = 1;
+      removeCoverageReduction = 0;
+      REQUIRE(model.multiplierForCoverage(t, cov_prop, removeCoverageReduction,
+                                          removeCoverageReductionTime,
+                                          graduallyRemoveCoverageReduction) ==
+              1);
+    }
+
+    SECTION("coverage multiplier should be linearly changed up until "
+            "removeCoverageReductionTime "
+            "when we gradually remove "
+            "coverage reduction ") {
+      int removeCoverageReduction = 0;
+      int removeCoverageReductionTime = 200;
+      int graduallyRemoveCoverageReduction = 1;
+      double cov_prop = 0.6;
+      int t = 0;
+      REQUIRE(model.multiplierForCoverage(t, cov_prop, removeCoverageReduction,
+                                          removeCoverageReductionTime,
+                                          graduallyRemoveCoverageReduction) ==
+              0.6);
+      graduallyRemoveCoverageReduction = 1;
+      t = 50;
+      REQUIRE(model.multiplierForCoverage(t, cov_prop, removeCoverageReduction,
+                                          removeCoverageReductionTime,
+                                          graduallyRemoveCoverageReduction) ==
+              0.7);
+      t = 100;
+      REQUIRE(model.multiplierForCoverage(t, cov_prop, removeCoverageReduction,
+                                          removeCoverageReductionTime,
+                                          graduallyRemoveCoverageReduction) ==
+              0.8);
+      t = 150;
+      REQUIRE(model.multiplierForCoverage(t, cov_prop, removeCoverageReduction,
+                                          removeCoverageReductionTime,
+                                          graduallyRemoveCoverageReduction) ==
+              0.9);
+      t = 200;
+      REQUIRE(model.multiplierForCoverage(t, cov_prop, removeCoverageReduction,
+                                          removeCoverageReductionTime,
+                                          graduallyRemoveCoverageReduction) ==
+              1);
+    }
+  }
 }