Merge pull request #20 from scrm/update/v1.5.1

Update to scrm v1.5.1

DESCRIPTION

@@ -1,12 +1,14 @@
 Package: scrm
 Type: Package
 Title: Simulating the Evolution of Biological Sequences
-Version: 1.5.0-1.9001
-Date: 2015-05-13
-Authors@R: c( person('Paul', 'Staab', , '[email protected]', role=c('aut',
-    'cre', 'cph')), person('Zhu', 'Sha', role=c('aut', 'cph')), person('Dirk',
-    'Metzler', role=c('aut', 'cph', 'ths')), person('Gerton', 'Lunter',
-    role=c('aut', 'cph', 'ths')) )
+Version: 1.5.1-1
+Date: 2015-05-20
+Authors@R: c(
+    person('Paul', 'Staab', , '[email protected]', role=c('aut', 'cre', 'cph')),
+    person('Zhu', 'Sha', role=c('aut', 'cph')),
+    person('Dirk', 'Metzler', role=c('aut', 'cph', 'ths')),
+    person('Gerton', 'Lunter', role=c('aut', 'cph', 'ths'))
+    )
 Description: A coalescent simulator that allows the rapid simulation of
     biological sequences under neutral models of evolution. Different to
     other coalescent based simulations, scrm has an optional approximation

NAMESPACE

@@ -1,4 +1,4 @@
-# Generated by roxygen2 (4.1.0): do not edit by hand
+# Generated by roxygen2 (4.1.1): do not edit by hand
 
 export(scrm)
 importFrom(Rcpp,evalCpp)

NEWS

@@ -1,14 +1,15 @@
-scrm 1.5.0-1.9001
+scrm 1.5.1-1
 ============
 
+* Updated the sources of scrm to version 1.5.1.
 * Updated the citation information.
 
 
 scrm 1.5.0-1
 ============
 
-* Updated scrm to version 1.5.0
-* Added citation information
+* Updated the sources of scrm to version 1.5.0.
+* Added citation information.
 * Replaced callback to warning with Rf_warning (#12). Thanks to J.J. Allaire for
   the patch.
 
@@ -17,7 +18,7 @@ scrm 1.3-3
 ==========
 
 * Added dependency to R >= 3.1.0.
-* Updated scrm source to version 1.3.1.
+* Updated the sources of scrm to version 1.3.1.
 * Added title to vignettes.
 
 

man/scrm-package.Rd

@@ -1,4 +1,4 @@
-% Generated by roxygen2 (4.1.0): do not edit by hand
+% Generated by roxygen2 (4.1.1): do not edit by hand
 % Please edit documentation in R/scrm.R
 \docType{package}
 \name{scrm-package}

man/scrm.Rd

@@ -1,4 +1,4 @@
-% Generated by roxygen2 (4.1.0): do not edit by hand
+% Generated by roxygen2 (4.1.1): do not edit by hand
 % Please edit documentation in R/RcppExports.R
 \name{scrm}
 \alias{scrm}

src/scrm/forest-debug.cc

@@ -496,12 +496,12 @@ std::vector<Node const*> Forest::determinePositions() const {
   }
 
   bool Forest::printNodes() const {
-    std::cout << std::setw(10) << std::right << "Node";
-    std::cout << std::setw(10) << std::right << "Height";
+    std::cout << std::setw(15) << std::right << "Node";
+    std::cout << std::setw(15) << std::right << "Height";
     std::cout << std::setw(6) << std::right << "label";
-    std::cout << std::setw(10) << std::right << "Parent";
-    std::cout << std::setw(10) << std::right << "1th_child";
-    std::cout << std::setw(10) << std::right << "2nd_child";
+    std::cout << std::setw(15) << std::right << "Parent";
+    std::cout << std::setw(15) << std::right << "1th_child";
+    std::cout << std::setw(15) << std::right << "2nd_child";
     std::cout << std::setw(6) << std::right << "local";
     std::cout << std::setw(6) << std::right << "pop";
     std::cout << std::setw(10) << std::right << "l_upd";
@@ -510,14 +510,14 @@ std::vector<Node const*> Forest::determinePositions() const {
     std::cout << std::endl;
 
     for(size_t i = 0; i < this->getNodes()->size(); ++i) {
-      std::cout << std::setw(10) << std::right << this->getNodes()->get(i);
-      std::cout << std::setw(10) << std::right << this->getNodes()->get(i)->height();
+      std::cout << std::setw(15) << std::right << this->getNodes()->get(i);
+      std::cout << std::setw(15) << std::right << this->getNodes()->get(i)->height();
       std::cout << std::setw(6) << std::right << this->getNodes()->get(i)->label();
       if (!getNodes()->get(i)->is_root()) 
-        std::cout << std::setw(10) << std::right << this->getNodes()->get(i)->parent();
-      else std::cout << std::setw(10) << std::right << 0;
-      std::cout << std::setw(10) << std::right << this->getNodes()->get(i)->first_child();
-      std::cout << std::setw(10) << std::right << this->getNodes()->get(i)->second_child();
+        std::cout << std::setw(15) << std::right << this->getNodes()->get(i)->parent();
+      else std::cout << std::setw(15) << std::right << 0;
+      std::cout << std::setw(15) << std::right << this->getNodes()->get(i)->first_child();
+      std::cout << std::setw(15) << std::right << this->getNodes()->get(i)->second_child();
       std::cout << std::setw(6) << std::right << this->getNodes()->get(i)->local();
       std::cout << std::setw(6) << std::right << this->getNodes()->get(i)->population();
       std::cout << std::setw(10) << std::right << this->getNodes()->get(i)->last_update();

src/scrm/forest.cc

@@ -414,7 +414,6 @@ void Forest::sampleNextGenealogy() {
     return;
   }
 
-  dout << tmp_event_time_ << std::endl;
   assert( tmp_event_time_ >= 0 );
   this->contemporaries_.buffer(tmp_event_time_);
 
@@ -444,9 +443,6 @@ void Forest::sampleNextGenealogy() {
   assert( this->printTree() );
   assert( this->printNodes() );
 
-  // Rarely prune all possible nodes 
-  if (segment_count_ % 100000 == 0 && model().exact_window_length() != -1) this->doCompletePruning();
-
   this->sampleNextBase();
   this->calcSegmentSumStats();
 }
@@ -1206,14 +1202,6 @@ void Forest::printLocusSumStats(std::ostream &output) const {
   }
 }
 
-void Forest::doCompletePruning() {
-  dout << "Pruning Tree:" << std::endl;
-  for (NodeIterator ni = nodes()->iterator(nodes()->first()->next()); ni.good(); ++ni) {
-    pruneNodeIfNeeded((*ni)->previous());
-  }
-  pruneNodeIfNeeded(nodes()->last());
-}
-
 void Forest::clear() {
   // Clear roots tracking
   set_local_root(NULL);

src/scrm/forest.h

@@ -245,7 +245,6 @@ class Forest
   }
 
   bool pruneNodeIfNeeded(Node* node, const bool prune_orphans = true);
-  void doCompletePruning();
 
   // Calculation of Rates
   double calcCoalescenceRate(const size_t pop, const TimeInterval &ti) const;

src/scrm/model.cc

@@ -496,17 +496,17 @@ void Model::addMigrationRates(double time, const std::vector<double> &mig_rates,
  */
 void Model::addSymmetricMigration(const double time, const double mig_rate, 
                                   const bool &time_scaled, const bool &rate_scaled) {
-  std::vector<double> mig_rates = std::vector<double>(population_number()*population_number(), mig_rate);
-  this->addMigrationRates(time, mig_rates, time_scaled, rate_scaled);
-}
+    std::vector<double> mig_rates = std::vector<double>(population_number()*population_number(), mig_rate);
+    this->addMigrationRates(time, mig_rates, time_scaled, rate_scaled);
+  }
 
 
-void Model::addSingleMigrationEvent(const double time, const size_t source_pop, 
-                                    const size_t sink_pop, const double fraction,
-                                    const bool &time_scaled) {
-
-  size_t position = addChangeTime(time, time_scaled);
-  size_t popnr = population_number();
+  void Model::addSingleMigrationEvent(const double time, const size_t source_pop, 
+                                      const size_t sink_pop, const double fraction,
+                                      const bool &time_scaled) {
+    
+    size_t position = addChangeTime(time, time_scaled);
+    size_t popnr = population_number();
 
   if ( time < 0.0 ) throw std::invalid_argument("Single migration event: Negative time");
   if ( source_pop >= population_number() ) throw std::invalid_argument("Single migration event: Unknown population");
@@ -564,7 +564,7 @@ std::ostream& operator<<(std::ostream& os, Model& model) {
       for (size_t j = 0; j < n_pops; ++j) {
         if (model.single_mig_pop(i, j) != 0) {
           os << " " << model.single_mig_pop(i, j) * 100 << "\% of pop " 
-             << i << " move to pop " << j << std::endl;
+             << i + 1 << " move to pop " << j + 1 << std::endl;
         }
       }
     }
@@ -782,3 +782,59 @@ void Model::addPopToVectorList(std::vector<std::vector<double>*> &vector_list) {
     (*it)->push_back(nan("value to replace"));
   }
 }
+
+
+/**
+ * @brief Sets the recombination rate
+ *
+ * @param rate The recombination rate per sequence length per time unit. 
+ * The sequence length can be either per locus or per base pair and the time
+ * can be given in units of 4N0 generations ("scaled") or per generation. 
+ *
+ * @param loci_length The length of every loci.
+ * @param per_locus Set to TRUE, if the rate is given per_locus, and to FALSE
+ * if the rate is per base pair.
+ * @param scaled Set to TRUE is the rate is scaled with 4N0, or to FALSE if
+ * it isn't
+ */
+void Model::setRecombinationRate(double rate,  
+                                 const bool &per_locus, 
+                                 const bool &scaled,
+                                 const double seq_position) { 
+
+  if (rate < 0.0) { 
+    throw std::invalid_argument("Recombination rate must be non-negative");
+  }
+
+  if (scaled) rate /= 4.0 * default_pop_size; 
+  if (per_locus) {
+    if (loci_length() <= 1) {
+      throw std::invalid_argument("Locus length must be at least two");
+    }
+    rate /= loci_length()-1;
+  }
+
+  recombination_rates_[addChangePosition(seq_position)] = rate; 
+}
+
+
+/**
+ * @brief Sets the mutation rate
+ *
+ * @param rate The mutation rate per locus, either scaled as theta=4N0*u or
+ * unscaled as u.
+ * @param per_locus TRUE if the rate is per locus, FALSE if per base.
+ * @param scaled Set this to TRUE if you give the mutation rate in scaled
+ * units (e.g. as theta rather than as u).
+ */
+void Model::setMutationRate(double rate, const bool &per_locus, const bool &scaled, 
+                            const double seq_position) { 
+  if (scaled) rate /= 4.0 * default_pop_size; 
+
+  size_t idx = addChangePosition(seq_position);
+  if (per_locus) {
+    mutation_rates_.at(idx) = rate / loci_length();
+  } else {
+    mutation_rates_.at(idx) = rate;
+  }
+}

src/scrm/model.h

@@ -496,59 +496,5 @@ std::ostream& operator<<(std::ostream& os, const std::vector<T> &vec) {
 }
 
 
-/**
- * @brief Sets the recombination rate
- *
- * @param rate The recombination rate per sequence length per time unit. 
- * The sequence length can be either per locus or per base pair and the time
- * can be given in units of 4N0 generations ("scaled") or per generation. 
- *
- * @param loci_length The length of every loci.
- * @param per_locus Set to TRUE, if the rate is given per_locus, and to FALSE
- * if the rate is per base pair.
- * @param scaled Set to TRUE is the rate is scaled with 4N0, or to FALSE if
- * it isn't
- */
-inline void Model::setRecombinationRate(double rate,  
-                                        const bool &per_locus, 
-                                        const bool &scaled,
-                                        const double seq_position) { 
-
-  if (rate < 0.0) { 
-    throw std::invalid_argument("Recombination rate must be non-negative");
-  }
-
-  if (scaled) rate /= 4.0 * default_pop_size; 
-  if (per_locus) {
-    if (loci_length() <= 1) {
-      throw std::invalid_argument("Locus length must be at least two");
-    }
-    rate /= loci_length()-1;
-  }
-
-  recombination_rates_[addChangePosition(seq_position)] = rate; 
-}
-
-
-/**
- * @brief Sets the mutation rate
- *
- * @param rate The mutation rate per locus, either scaled as theta=4N0*u or
- * unscaled as u.
- * @param per_locus TRUE if the rate is per locus, FALSE if per base.
- * @param scaled Set this to TRUE if you give the mutation rate in scaled
- * units (e.g. as theta rather than as u).
- */
-inline void Model::setMutationRate(double rate, const bool &per_locus, const bool &scaled, 
-                                   const double seq_position) { 
-  if (scaled) rate /= 4.0 * default_pop_size; 
-
-  size_t idx = addChangePosition(seq_position);
-  if (per_locus) {
-    mutation_rates_.at(idx) = rate / loci_length();
-  } else {
-    mutation_rates_.at(idx) = rate;
-  }
-}
 
 #endif

src/scrm/node.cc

@@ -61,7 +61,13 @@ void Node::change_child(Node* from, Node* to) {
   }
   else if ( this->second_child() == from )
     this->set_second_child(to);
-  else throw std::invalid_argument("Can't find child node to replace");
+  else {
+    dout << "Error when changing child of " << this << " form "
+         << from << " to " << to << std::endl;
+    dout << "Children are " << this->first_child() << " and "
+         << this->second_child() << std::endl;
+    throw std::invalid_argument("Can't find child node to replace");
+  }
 }
 
 void Node::remove_child(Node* child) {

vignettes/scrm-Arguments.Rmd

@@ -1,5 +1,5 @@
 ---
-title: "scrm's Arguments"
+title: "Command Line Arguments"
 author: "Paul Staab"
 date: "`r Sys.Date()`"
 output: rmarkdown::html_vignette
@@ -107,10 +107,13 @@ Viewed backwards in time, this moves a fraction _1-p_ of the linages in populati
 
 ## Sequence specific parameters
 The following commands change the model parameters from at a sequence position _s_.
+You should still set the initial rate with `-r` or `-t`, respectively, and then use
+the commands prefixed with `s` for all changes. Note that `-r` also takes the total
+length of the sequence as second argument, while `-sr` just has the rate as argument.
 
 * `-sr <s> <R>`: Set the recombination rate to _R_ starting at position _s_.
 * `-st <s> <$\theta$>`: Set the mutation rate to $\theta$ starting at position _s_. 
 
 
 [1]: http://home.uchicago.edu/~rhudson1/source/mksamples.html
-[2]: https://webshare.uchicago.edu/users/rhudson1/Public/ms.folder/msdoc.pdf?ticket=t_DunQ7c99
+[2]: https://webshare.uchicago.edu/users/rhudson1/Public/ms.folder/msdoc.pdf?ticket=t_DunQ7c99