Merge pull request #21 from scrm/update/v1.6.0

Update to v1.6.0

DESCRIPTION

@@ -1,8 +1,8 @@
 Package: scrm
 Type: Package
 Title: Simulating the Evolution of Biological Sequences
-Version: 1.5.1-1
-Date: 2015-05-20
+Version: 1.6.0-1
+Date: 2015-06-05
 Authors@R: c(
     person('Paul', 'Staab', , '[email protected]', role=c('aut', 'cre', 'cph')),
     person('Zhu', 'Sha', role=c('aut', 'cph')),

NEWS

@@ -1,6 +1,12 @@
 scrm 1.5.1-1
 ============
 
+* Updated the sources of scrm to version 1.6.0.
+
+
+scrm 1.5.1-1
+============
+
 * Updated the sources of scrm to version 1.5.1.
 * Updated the citation information.
 

src/Makevars

@@ -15,4 +15,4 @@ OBJECTS = $(OBJECTS.scrmr) $(OBJECTS.scrm) $(OBJECTS.random) \
           $(OBJECTS.sumstats)
 
 CXX_STD = CXX11
-PKG_CPPFLAGS = -DVERSION="\"1.5.0\"" -DRBUILD
+PKG_CPPFLAGS = -DVERSION="\"R\"" -DRBUILD

src/Makevars.win

@@ -15,4 +15,4 @@ OBJECTS = $(OBJECTS.scrmr) $(OBJECTS.scrm) $(OBJECTS.random) \
           $(OBJECTS.sumstats)
 
 CXX_STD = CXX11
-PKG_CPPFLAGS = -DVERSION="\"1.5.0\"" -DRBUILD
+PKG_CPPFLAGS = -DVERSION="\"R\"" -DRBUILD

src/scrm/forest-debug.cc

@@ -24,11 +24,16 @@
 
 void Forest::createExampleTree() {
   this->clear();
+  this->writable_model()->disable_approximation();
   // Only set the number of samples to 4, but keep rest of the model
   this->writable_model()->sample_times_.clear();
   this->writable_model()->sample_populations_.clear();
   this->writable_model()->addSampleSizes(0.0, std::vector<size_t>(1, 4));
 
+  this->rec_bases_.push_back(5.0);  
+  this->current_rec_ = 1;
+  //this->rec_bases_.push_back(105.0);  
+
   Node* leaf1 = nodes()->createNode(0, 1);
   Node* leaf2 = nodes()->createNode(0, 2);
   Node* leaf3 = nodes()->createNode(0, 3);
@@ -57,9 +62,9 @@ void Forest::createExampleTree() {
 
   // Add a non-local tree
   Node* nl_node = nodes()->createNode(4); 
-  nl_node->make_nonlocal(5);
+  nl_node->make_nonlocal(current_rec_);
   Node* nl_root = nodes()->createNode(6);
-  nl_root->make_nonlocal(5);
+  nl_root->make_nonlocal(current_rec_);
 
   nl_node->set_parent(nl_root);
   nl_root->set_first_child(nl_node);
@@ -72,8 +77,6 @@ void Forest::createExampleTree() {
   updateAbove(leaf3);
   updateAbove(leaf4);
 
-  this->set_current_base(5);  
-  this->set_next_base(105);
   this->set_sample_size(4);
 
   this->contemporaries_ = ContemporariesContainer(model().population_number(), 

src/scrm/forest.cc

@@ -40,7 +40,11 @@ void Forest::initialize(Model* model,
 
   this->set_model(model);
   this->set_random_generator(rg);
-  this->set_current_base(1);
+
+  current_rec_ = 0;
+  rec_bases_ = std::vector<double>(1, -1);
+  rec_bases_.reserve(1000);
+
   this->set_sample_size(0);
 
   this->coalescence_finished_ = true;
@@ -75,9 +79,8 @@ Forest::Forest(const Forest &current_forest) {
 
   // Copy state information
   this->set_sample_size(current_forest.sample_size());
-  this->set_current_base(current_forest.current_base());
-  this->set_next_base(current_forest.next_base());
-  this->segment_count_ = current_forest.segment_count_;
+  this->rec_bases_ = current_forest.rec_bases_;
+  this->current_rec_ = current_forest.current_rec_;
 
   // Copy the nodes
   this->nodes_ = NodeContainer(*current_forest.getNodes());
@@ -126,7 +129,7 @@ Node* Forest::cut(const TreePoint &cut_point) {
   if ( !cut_point.base_node()->local() )
     new_leaf->make_nonlocal(cut_point.base_node()->last_update());
   else
-    new_leaf->make_nonlocal(current_base());
+    new_leaf->make_nonlocal(current_rec());
   assert( !new_leaf->local() );
 
   new_leaf->set_population(cut_point.base_node()->population());
@@ -197,7 +200,7 @@ void Forest::updateAbove(Node* node, bool above_local_root,
   // Fast forward above local root because this part is fairly straight forward
   if (above_local_root) {
     // Assure that everything is non-local
-    if (node->local()) node->make_nonlocal(current_base());
+    if (node->local()) node->make_nonlocal(current_rec());
 
     // Update the primary root if needed
     if ( node->is_root() ) {
@@ -210,7 +213,7 @@ void Forest::updateAbove(Node* node, bool above_local_root,
     return;
   }
 
-  node->set_last_change(current_base());
+  node->set_last_change(current_rec());
 
   // Calculate new values for samples_below and length_below for the current
   // node
@@ -237,10 +240,10 @@ void Forest::updateAbove(Node* node, bool above_local_root,
   // Update whether the node is local or not 
   if (!invariants_only) {
     if (samples_below == 0) {
-      if ( node->local() ) node->make_nonlocal(current_base());
+      if ( node->local() ) node->make_nonlocal(current_rec());
     }
     else if ( samples_below == sample_size() ) {
-      if ( node->local() ) node->make_nonlocal(current_base());
+      if ( node->local() ) node->make_nonlocal(current_rec());
 
       // Are we the local root?
       if (node->countChildren() == 2 && 
@@ -278,9 +281,12 @@ void Forest::updateAbove(Node* node, bool above_local_root,
  */
 void Forest::buildInitialTree() {
   dout << "===== BUILDING INITIAL TREE =====" << std::endl;
-  assert( this->nodes()->size() == 0 );
-  this->set_current_base(0.0);
-  this->segment_count_ = 1;
+  assert(this->nodes()->size() == 0);
+  assert(this->segment_count() == 0);
+  assert(this->rec_bases_.size() == 1);
+  assert(this->model().getCurrentSequencePosition() == 0.0);
+  this->set_next_base(0.0);
+  ++current_rec_;
 
   dout << "* Adding first node... ";
   Node* first_node = nodes()->createNode(model().sample_time(0), 1);
@@ -314,6 +320,7 @@ void Forest::buildInitialTree() {
     assert(this->printNodes());
   }
   this->sampleNextBase();
+  dout << "Next Sequence position: " << this->next_base() << std::endl;
   this->calcSegmentSumStats();
 }
 
@@ -404,23 +411,26 @@ TreePoint Forest::samplePoint(Node* node, double length_left) {
  * @ingroup group_pf_update
  */
 void Forest::sampleNextGenealogy() {
-  this->set_current_base(next_base_);
+  ++current_rec_; // Move to next recombination;
 
-  if (current_base_ == model().getCurrentSequencePosition()) {
+  if (current_base() == model().getCurrentSequencePosition()) {
     // Don't implement a recombination if we are just here because rates changed
     dout << std::endl << "Position: " << this->current_base() << ": Changing rates." << std::endl;
     this->sampleNextBase();
     this->calcSegmentSumStats();
+    dout << "Next Position: " << this->next_base() << std::endl;
     return;
   }
 
+  assert( current_base() > model().getCurrentSequencePosition() );
+  assert( current_base() < model().getNextSequencePosition() );
+
   assert( tmp_event_time_ >= 0 );
   this->contemporaries_.buffer(tmp_event_time_);
 
-  ++segment_count_;
   dout << std::endl << "===== BUILDING NEXT GENEALOGY =====" << std::endl;
-  dout << "Segment Nr: " << segment_count_ << std::endl;
-  dout << "Sequence position: " << this->current_base() << std::endl;
+  dout << "Segment Nr: " << current_rec() << " | "
+       << "Sequence position: " << this->current_base() << std::endl;
   assert( this->current_base() <= this->model().loci_length() );
 
   // Sample the recombination point
@@ -444,6 +454,7 @@ void Forest::sampleNextGenealogy() {
   assert( this->printNodes() );
 
   this->sampleNextBase();
+  dout << "Next Sequence position: " << this->next_base() << std::endl;
   this->calcSegmentSumStats();
 }
 
@@ -466,7 +477,7 @@ void Forest::sampleCoalescences(Node *start_node) {
   coalescence_finished_ = false;
 
   // This assertion needs an exception for building the initial tree
-  assert ( segment_count_ == 1 || 
+  assert ( current_rec() == 1 || 
            active_node(1)->in_sample() || 
            start_node->height() <= active_node(1)->height() );
 
@@ -609,6 +620,10 @@ void Forest::calcRates(const TimeInterval &ti) {
     // recombining
     rates_[0] += calcRecombinationRate(active_node(1));
   }
+
+  assert(rates_[0] >= 0);
+  assert(rates_[1] >= 0);
+  assert(rates_[2] >= 0);
 }  
 
 
@@ -763,24 +778,27 @@ double Forest::calcPwCoalescenceRate(const size_t pop, const TimeInterval &ti) c
 
 
 double Forest::calcRecombinationRate(Node const* node) const {
-  assert( !node->local() );
-  if (node->last_update() >= model().getCurrentSequencePosition()) {
+  assert(!node->local());
+  assert(this->current_base() >= model().getCurrentSequencePosition());
+  double last_update_pos = get_rec_base(node->last_update());
+
+  if (last_update_pos >= model().getCurrentSequencePosition()) {
     // Rec rate is constant for the relevant sequence part
-    return ( model().recombination_rate() * (this->current_base() - node->last_update()) );
+    return (this->current_base() - last_update_pos) * model().recombination_rate();
   } else {
     // Rec rate may change. Accumulate the total rate.
 
     double rate = model().recombination_rate() * 
         (this->current_base() - model().getCurrentSequencePosition()); 
     size_t idx = model().get_position_index() - 1;
 
-    while (model().change_position(idx) > node->last_update()) {
+    while (model().change_position(idx) > last_update_pos) {
       assert(idx > 0);
       rate += model().recombination_rate(idx) * (model().change_position(idx+1)-model().change_position(idx));
       --idx;
     }
 
-    rate += model().recombination_rate(idx) * (model().change_position(idx+1)-node->last_update());
+    rate += model().recombination_rate(idx) * (model().change_position(idx+1)-last_update_pos);
     return rate;
   }
 }
@@ -1128,7 +1146,7 @@ Node* Forest::possiblyMoveUpwards(Node* node, const TimeInterval &time_interval)
 
 bool Forest::pruneNodeIfNeeded(Node* node, const bool prune_orphans) {
   assert(node != NULL);
-  if (model().exact_window_length() == -1) return false;
+  if (!model().has_approximation()) return false;
   if (node->in_sample()) return false;
 
   if (node->is_root()) {
@@ -1211,11 +1229,16 @@ void Forest::clear() {
   nodes()->clear();
 
   // Reset Position & Segment Counts
-  this->set_current_base(0.0);
-  this->segment_count_ = 0;
+  this->rec_bases_.clear();
+  this->set_next_base(-1.0);
+  this->current_rec_ = 0;
 
   // Clear Summary Statistics
   this->clearSumStats();
+
+  // Reset Model
+  writable_model()->resetTime();
+  writable_model()->resetSequencePosition();
 }
 
 
@@ -1264,7 +1287,7 @@ Node* Forest::readNewickNode( std::string &in_str, std::string::iterator &it, si
     } else if ( (*(it)) == ';' ) {
       dout <<" Node " << node << " closed " << std::endl;
       this->nodes()->add( node );
-      node->make_nonlocal( this->current_base() );
+      node->make_nonlocal(current_rec());
       return node;
     } else {
       continue;
@@ -1276,8 +1299,7 @@ Node* Forest::readNewickNode( std::string &in_str, std::string::iterator &it, si
 
 
 void Forest::readNewick( std::string &in_str ){
-  this->set_current_base(0.0);
-  this->segment_count_ = 1;
+  this->current_rec_ = 1;
   std::string::iterator it = in_str.begin();
   (void)this->readNewickNode( in_str, it );
   this->set_local_root( this->nodes()->last() );
@@ -1289,8 +1311,24 @@ void Forest::readNewick( std::string &in_str ){
   }
   assert(this->printNodes());
   assert(this->printTree());
-dout << "contemporaries_.size()"<<contemporaries_.size(0) <<std::endl;
+  dout << "contemporaries_.size()"<<contemporaries_.size(0) <<std::endl;
   this->sampleNextBase();
   this->calcSegmentSumStats();
   this->tmp_event_time_ = this->local_root()->height();
 }
+
+
+// Must be called AFTER the tree was modified.
+void Forest::sampleNextBase() {
+  double length = random_generator()->sampleExpoLimit(getLocalTreeLength() * model().recombination_rate(),
+                                                      model().getNextSequencePosition() - current_base());
+  if (length == -1) {
+    // No recombination until the model changes
+    set_next_base(model().getNextSequencePosition());
+    if (next_base() < model().loci_length()) writable_model()->increaseSequencePosition();
+  } else {
+    // Recombination in the sequence segment
+    set_next_base(current_base() + length);
+  }
+} 
+