Add R0 tracking

Author: emilydolson

emsdk

@@ -59,7 +59,7 @@ EMP_BUILD_CONFIG(SymConfigBase,
     VALUE(WRITE_CURRENT_INTERACTION_COUNTS, bool, 0, "Should the world write the count of only-currently-present interactions? (0 for no, 1 for yes)"),
     VALUE(PHYLOGENY_SNAPSHOT_INTERVAL, int, 10001, "How often to output phylogeny snapshots"),
     VALUE(NUM_PHYLO_BINS, size_t, 5, "How many bins should organisms be separated into if phylogeny is on?"),
-    VALUE(PHYLOGENY_TAXON_TYPE, std::string, "interaction-value-binned", "What are phylogeny taxa based on? Options: interaction-value-binned, interaction-value-exact, tag, individual"),
+    VALUE(PHYLOGENY_TAXON_TYPE, std::string, "interaction-value-binned", "What are phylogeny taxa based on? Options: interaction-value-binned, interaction-value-exact, tag, individual, horizontal-clade"),
     VALUE(STORE_EXTINCT, bool, 0, "Should extinct taxa be stored? (0 for no, 1 for yes)"),
 
     GROUP(MUTATION, "Mutation"),
@@ -95,6 +95,9 @@ EMP_BUILD_CONFIG(SymConfigBase,
     VALUE(TAG_MUTATION_SIZE, double, 0.01, "What is the probability that any given position in the bitstring tag flips during mutation?"),
     VALUE(WRITE_TAG_MATRIX, bool, 0, "At the end of the experiment, should a similarity matrix of all persisting tags be generated?"),
     VALUE(TAG_MATRIX_SAMPLE_PROPORTION, double, 0.1, "What proportion of positions in the world should be sampled to produce the tag matrix from?"),
-    VALUE(HOST_STARTING_TAGS_ONE_PROB, double, 0, "What probability should initializing bits in tags have of being 1s? Hosted symbionts will be assigned their host's tag. (0 for basic, all-0 only tags)")
+    VALUE(HOST_STARTING_TAGS_ONE_PROB, double, 0, "What probability should initializing bits in tags have of being 1s? Hosted symbionts will be assigned their host's tag. (0 for basic, all-0 only tags)"),
+
+    GROUP(ANTIBIOTIC_RESISTANCE, "Settings for antibiotic resistance"),
+    VALUE(R0_WINDOW, double, 100, "How many updates should be used to calculate R0?")
 )
 #endif

source/ConfigSetup.h

@@ -29,6 +29,8 @@ class Organism {
     return !(*this == other);
   }
 
+  uint64_t taxon_id;
+
   virtual std::string const GetName() const {
     std::cout << "GetName called from Organism" << std::endl;
     throw "Organism method called!";

source/Organism.h

@@ -21,6 +21,7 @@ void SymWorld::CreateDataFiles() {
   SetupTransmissionFile(my_config->FILE_PATH()+"TransmissionRates"+my_config->FILE_NAME()+file_ending).SetTimingRepeat(TIMING_REPEAT);
   SetupSymDiversityFile(my_config->FILE_PATH()+"SymDiversity"+my_config->FILE_NAME()+file_ending).SetTimingRepeat(TIMING_REPEAT);
   SetupReproHistFile(my_config->FILE_PATH() + "ReproHist" + my_config->FILE_NAME() + file_ending).SetTimingRepeat(TIMING_REPEAT);
+  SetupAntibioticResistanceFile(my_config->FILE_PATH() + "AntibioticResistance" + my_config->FILE_NAME() + file_ending).SetTimingRepeat(TIMING_REPEAT);  
   if (my_config->FREE_LIVING_SYMS() == 1) {
     SetupFreeLivingSymFile(my_config->FILE_PATH()+"FreeLivingSyms_"+my_config->FILE_NAME()+file_ending).SetTimingRepeat(TIMING_REPEAT);
   }
@@ -228,6 +229,50 @@ emp::DataFile& SymWorld::SetupReproHistFile(const std::string& filename) {
   return file;
 }
 
+emp::DataFile & SymWorld::SetupAntibioticResistanceFile(const std::string & filename) {
+  auto & file = SetupFile(filename);
+  // auto & node = GetR0();
+
+  file.AddVar(update, "update", "Update");
+
+  // node.Reset();
+  file.AddFun((std::function<double()>) [this](){
+    sym_sys->RemoveAncestorsBefore(GetUpdate() - GetConfig()->R0_WINDOW());
+    sym_sys->RemoveOutsideBefore(GetUpdate() - GetConfig()->R0_WINDOW());
+
+    double numerator = 0;
+    double denominator = 0;
+    // This is going to bias R0 to be high, since we're not counting horizontal transmissions that
+    // quickly died out
+    for (auto sym_taxon : sym_sys->GetActive()) {
+      denominator += sym_taxon->GetNumOffEver();
+      numerator += 1;
+    }
+    // std::cout << numerator << " " << denominator << std::endl;
+    for (auto sym_taxon : sym_sys->GetAncestors()) {
+      denominator += sym_taxon->GetNumOffEver();
+      numerator += 1;
+    }
+    // std::cout << numerator << " " << denominator << std::endl;
+    for (auto sym_taxon : sym_sys->GetOutside()) {
+      denominator += sym_taxon->GetNumOffEver();
+      numerator += 1;
+    }
+    // std::cout << numerator << " " << denominator << std::endl <<std::endl;
+    return numerator / denominator;
+  }, "R0", "R0");
+  // node.AddDatum(numerator / denominator);
+
+  // file.AddMean(node, "R0", "R0");
+  // file.AddMax(node, "max_intval", "Maximum symbiont interaction value");
+  // file.AddMin(node, "min_intval", "Minimum symbiont interaction value");
+  // file.AddTotal(node1, "count", "Total number of symbionts");
+
+  file.PrintHeaderKeys();
+
+  return file;
+}
+
 
 /**
  * Input: The address of the string representing the suffixes for the files to be created.
@@ -1413,4 +1458,20 @@ emp::DataMonitor<double>& SymWorld::GetHostTagPermissiveness() {
     return *data_node_symbiont_tag_shannon;
   }
 
+  /**
+   * Input: None
+   *
+   * Output: The DataMonitor<double>& that has the information representing
+   * the symbiont R0 this update.
+   *
+   * Purpose: To retrieve the data node that is tracking the
+   * R0.
+   */
+  emp::DataMonitor<double>& SymWorld::GetR0() {
+    if (!data_node_R0) {
+      data_node_R0.New();
+    }      
+    return *data_node_R0;
+  }
+
 #endif

source/default_mode/DataNodes.h

@@ -46,7 +46,7 @@ class SymWorld : public emp::World<Organism> {
   enum class SPATIAL_STRUCT_MODE { WELL_MIXED, GRID, LOAD };
   static const std::unordered_map<std::string, SPATIAL_STRUCT_MODE> spatial_struct_mode_cfg_mapping;
 
-  enum class PHYLO_TAXON_TYPE { INTERACTION_VALUE_BINNED, INTERACTION_VALUE_EXACT, TAG, INDIVIDUAL };
+  enum class PHYLO_TAXON_TYPE { INTERACTION_VALUE_BINNED, INTERACTION_VALUE_EXACT, TAG, INDIVIDUAL, HORIZONTAL_CLADE };
   static const std::unordered_map<std::string, PHYLO_TAXON_TYPE> phylo_taxon_type_cfg_mapping;
 
   enum class TAG_METRIC_TYPE { HAMMING, STREAK, HASH };
@@ -190,6 +190,7 @@ class SymWorld : public emp::World<Organism> {
   emp::Ptr<emp::DataMonitor<int>> data_node_freesymcount;
   emp::Ptr<emp::DataMonitor<int>> data_node_hostedsymcount;
   emp::Ptr<emp::DataMonitor<int>> data_node_uninf_hosts;
+  emp::Ptr<emp::DataMonitor<double>> data_node_R0;
   emp::Ptr<emp::DataMonitor<double, emp::data::Histogram>> data_node_attempts_horiztrans;
   emp::Ptr<emp::DataMonitor<double, emp::data::Histogram>> data_node_tagfail_horiztrans;
   emp::Ptr<emp::DataMonitor<double, emp::data::Histogram>> data_node_sizefail_horiztrans;
@@ -343,6 +344,7 @@ class SymWorld : public emp::World<Organism> {
     if (data_node_successes_horiztrans) data_node_successes_horiztrans.Delete();
     if (data_node_attempts_verttrans) data_node_attempts_verttrans.Delete();
     if (data_node_successes_verttrans) data_node_successes_verttrans.Delete();
+    if (data_node_R0) data_node_R0.Delete();
 
     for (size_t i = 0; i < sym_pop.size(); i++) { //host population deletion is handled by empirical world destructor
       if (sym_pop[i]) {
@@ -1040,6 +1042,7 @@ class SymWorld : public emp::World<Organism> {
   emp::DataFile& SetupTransmissionFile(const std::string& filename);
   emp::DataFile& SetupTagDistFile(const std::string& filename);
   emp::DataFile& SetupSymDiversityFile(const std::string& filename);
+  emp::DataFile& SetupAntibioticResistanceFile(const std::string& filename);
   virtual void SetupTransmissionFileColumns(emp::DataFile& file);
   virtual void SetupHostFileColumns(emp::DataFile& file);
   emp::DataMonitor<int>& GetHostCountDataNode();
@@ -1064,6 +1067,7 @@ class SymWorld : public emp::World<Organism> {
   emp::DataMonitor<double>& GetHostTagShannonDiversity();
   emp::DataMonitor<int>& GetSymbiontTagRichness();
   emp::DataMonitor<double>& GetSymbiontTagShannonDiversity();
+  emp::DataMonitor<double>& GetR0();
   emp::DataMonitor<double, emp::data::Histogram>& GetTagDistanceDataNode();
   emp::DataMonitor<double, emp::data::Histogram>& GetHostIntValDataNode();
   emp::DataMonitor<double, emp::data::Histogram>& GetSymIntValDataNode();
@@ -1503,7 +1507,8 @@ const std::unordered_map<
     {"interaction-value-binned", PHYLO_TAXON_TYPE::INTERACTION_VALUE_BINNED},
     {"interaction-value-exact", PHYLO_TAXON_TYPE::INTERACTION_VALUE_EXACT},
     {"tag", PHYLO_TAXON_TYPE::TAG},
-    {"individual", PHYLO_TAXON_TYPE::INDIVIDUAL}
+    {"individual", PHYLO_TAXON_TYPE::INDIVIDUAL},
+    {"horizontal-clade", PHYLO_TAXON_TYPE::HORIZONTAL_CLADE}
 };
 
 const std::unordered_map<

source/default_mode/SymWorld.h

@@ -123,6 +123,13 @@ class Symbiont: public Organism {
    */
   emp::WorldPosition location;
 
+  /**
+   *
+   * Purpose: Janky approach to tracking R0
+   */
+  bool next_reproduction_is_horizontal = false;
+  uint64_t taxon_id;
+
 public:
   /**
    * The constructor for symbiont
@@ -698,6 +705,12 @@ class Symbiont: public Organism {
     sym_baby->Mutate();
     sym_baby->SetReproCount(reproductions + 1);
     if(my_config->PHYLOGENY() == 1) {
+      if (next_reproduction_is_horizontal) {
+        sym_baby->taxon_id = my_world->GetSymSys()->GetNextID();
+      } else {
+        sym_baby->taxon_id = my_taxon->GetID();
+      }
+
       my_world->AddSymToSystematic(sym_baby, my_taxon);
       //baby's taxon will be set in AddSymToSystematic
     }
@@ -860,7 +873,9 @@ class Symbiont: public Organism {
    */
   void IndependentReproduction(emp::WorldPosition location) {
     if (AttemptIndependentReproduction(location)) {
+      next_reproduction_is_horizontal = true;
       emp::Ptr<Organism> sym_baby = Reproduce();
+      next_reproduction_is_horizontal = false;
       if (my_config->TAG_MATCHING() || my_config->FREE_HT_FAILURE()) sym_baby->SetPoints(0);
       emp::WorldPosition new_pos = my_world->SymDoBirth(sym_baby, location);
 

source/default_mode/Symbiont.h

@@ -214,6 +214,15 @@ void SymWorld::SetupPhylogenyTracking() {
         return (long unsigned) sym_sys->GetNextID();
       };
       break;
+    case PHYLO_TAXON_TYPE::HORIZONTAL_CLADE:
+      calc_host_info_fun = [&](Organism& org) {
+        return 0;
+      };
+
+      calc_sym_info_fun = [&](Organism& org) {
+        return org.taxon_id;
+      };
+      break;
     default:
       emp_error("Unimplemented phylogeny taxon type.");
       break;