Add option for driver effects

heavywatal · Mar 2, 2016 · 8caab26 · 8caab26
1 parent 9cddf4d
commit 8caab26
Show file tree

Hide file tree

Showing 2 changed files with 33 additions and 20 deletions.
diff --git a/cell.cpp b/cell.cpp
@@ -14,6 +14,7 @@
 
 double Cell::MUTATION_RATE_ = 1e-1;
 double Cell::MUTATION_SIGMA_ = 0.0;
+int Cell::DRIVER_EFFECTS_ = 1;
 double Cell::DRIVER_FRACTION_ = 0.0;
 double Cell::BIRTH_RATE_ = 1.0;
 double Cell::DEATH_RATE_ = 0.0;
@@ -27,22 +28,25 @@ std::vector<size_t> Cell::MUTANT_IDS_;
 //! Program options
 /*! @return Program options description
 
-    Command line option | Symbol       | Variable
-    --------------------| ------------ | -------------------------
-    `-u,--mutation`     | \f$\mu\f$    | Cell::MUTATION_RATE_
-    `-s,--sigma`        | \f$\sigma\f$ | Cell::MUTATION_SIGMA_
-    `-b,--birth`        | \f$b\f$      | Cell::BIRTH_RATE_
-    `-d,--death`        | \f$d\f$      | Cell::DEATH_RATE_
-    `-m,--migration`    | \f$m\f$      | Cell::MIGRATION_RATE_
-    `-k,--shape`        | \f$k\f$      | Cell::GAMMA_SHAPE_
+    Command line option | Symbol         | Variable
+    --------------------| -------------- | -------------------------
+    `-u,--mutation`     | \f$\mu\f$      | Cell::MUTATION_RATE_
+    `-s,--sigma`        | \f$\sigma\f$   | Cell::MUTATION_SIGMA_
+    `-e,--effect`       |                | Cell::DRIVER_EFFECTS_
+    `-f,--fraction`     |                | Cell::DRIVER_FRACTION_
+    `-b,--birth`        | \f$\beta_0\f$  | Cell::BIRTH_RATE_
+    `-d,--death`        | \f$\delta_0\f$ | Cell::DEATH_RATE_
+    `-m,--migration`    | \f$m\f$        | Cell::MIGRATION_RATE_
+    `-k,--shape`        | \f$k\f$        | Cell::GAMMA_SHAPE_
 */
 boost::program_options::options_description& Cell::opt_description() {
     namespace po = boost::program_options;
     static po::options_description desc{"Cell"};
     desc.add_options()
         ("mutation,u", po::value<double>(&MUTATION_RATE_)->default_value(MUTATION_RATE_))
         ("sigma,s", po::value<double>(&MUTATION_SIGMA_)->default_value(MUTATION_SIGMA_))
-        ("driver,f", po::value<double>(&DRIVER_FRACTION_)->default_value(DRIVER_FRACTION_))
+        ("effect,e", po::value<int>(&DRIVER_EFFECTS_)->default_value(DRIVER_EFFECTS_))
+        ("fraction,f", po::value<double>(&DRIVER_FRACTION_)->default_value(DRIVER_FRACTION_))
         ("birth,b", po::value<double>(&BIRTH_RATE_)->default_value(BIRTH_RATE_))
         ("death,d", po::value<double>(&DEATH_RATE_)->default_value(DEATH_RATE_))
         ("migration,m", po::value<double>(&MIGRATION_RATE_)->default_value(MIGRATION_RATE_))
@@ -65,7 +69,12 @@ void Cell::mutate() {
     sites_.push_back(MUTATION_EFFECTS_.size());
     MUTATION_EFFECTS_.push_back(effect);
     MUTANT_IDS_.push_back(id());
-    fitness_ *= (effect += 1.0);
+    if (DRIVER_EFFECTS_ & 0b01) {
+        birth_rate_ *= (effect += 1.0);
+    }
+    if (DRIVER_EFFECTS_ & 0b10) {
+        death_rate_ *= (1.0 - effect);
+    }
 }
 
 double Cell::delta_time() {
@@ -110,7 +119,8 @@ std::string Cell::header(const size_t dimensions, const std::string& sep) {
         << "birth" << sep << "death" << sep;
     std::vector<std::string> axes{"x", "y", "z"};
     axes.resize(dimensions);
-    wtl::ost_join(oss, axes, sep) << sep << "sites" << sep << "fitness\n";
+    wtl::ost_join(oss, axes, sep) << sep << "sites" << sep
+        << "beta" << sep << "delta\n";
     return oss.str();
 }
 
@@ -119,7 +129,7 @@ std::ostream& Cell::write(std::ostream& ost, const std::string& sep) const {
         << time_of_birth_ << sep << time_of_death_ << sep;
     wtl::ost_join(ost, coord(), sep) << sep;
     wtl::ost_join(ost, sites(), ":") << sep
-        << fitness() << "\n";
+        << birth_rate_ << death_rate_ << "\n";
     return ost;
 }
 

diff --git a/cell.h b/cell.h
@@ -34,7 +34,8 @@ class Cell {
     Cell(const std::vector<int>& v): coord_(v), id_(++ID_TAIL_), ancestor_(id_) {}
     //! Copy constructor
     Cell(const Cell& other):
-        coord_(other.coord_), sites_(other.sites_), fitness_(other.fitness_),
+        coord_(other.coord_), sites_(other.sites_),
+        birth_rate_(other.birth_rate_), death_rate_(other.death_rate_),
         id_(++ID_TAIL_), mother_(other.id_), ancestor_(other.ancestor_) {}
     //! Destructor
     ~Cell() {--ID_TAIL_;}
@@ -71,8 +72,6 @@ class Cell {
     //! convert site positions to 01 vector
     std::vector<size_t> haplotype(std::vector<size_t> segsites) const;
 
-    double fitness() const {return fitness_;}
-
     static const std::vector<double>& MUTATION_EFFECTS() {return MUTATION_EFFECTS_;}
     static const std::vector<size_t>& MUTANT_IDS() {return MUTANT_IDS_;}
 
@@ -86,9 +85,9 @@ class Cell {
 
   private:
     //! finite per capita rates
-    double birth_rate() const {return BIRTH_RATE_ * fitness_;}
-    double death_rate() const {return DEATH_RATE_;}
-    double growth_rate() const {return 1.0 + birth_rate() - death_rate();}
+    double birth_rate() const {return birth_rate_;}
+    double death_rate() const {return death_rate_;}
+    double growth_rate() const {return 1.0 + birth_rate_ - death_rate_;}
     //! instantaneous rate for time increment
     double instantaneous_event_rate() const {
         const double lambda = growth_rate();
@@ -100,6 +99,9 @@ class Cell {
 
     static double MUTATION_SIGMA_;
 
+    //! 1: birth, 2: death, 3: both
+    static int DRIVER_EFFECTS_;
+
     static double DRIVER_FRACTION_;
 
     static double BIRTH_RATE_;
@@ -120,8 +122,9 @@ class Cell {
     std::vector<int> coord_;
     //! Mutated sites (infinite-site model)
     std::vector<size_t> sites_;
-    //! Sum of mutation effects + 1
-    double fitness_ = 1.0;
+
+    double birth_rate_ = BIRTH_RATE_;
+    double death_rate_ = DEATH_RATE_;
 
     //! Extra data
     size_t id_ = 0;