Reverse BFS order

Author: melkebir

result/sims/neighbor_joining/results_m5.txt

@@ -11,3 +11,32 @@ mS,8,Neighbor joining,8
 mS,9,Neighbor joining,3
 S,0,Neighbor joining,11
 S,12,Neighbor joining,9
+S,15,Neighbor joining,10
+S,17,Neighbor joining,6
+S,2,Neighbor joining,8
+S,23,Neighbor joining,11
+S,3,Neighbor joining,9
+S,4,Neighbor joining,11
+S,5,Neighbor joining,12
+S,9,Neighbor joining,12
+M,132,Neighbor joining,14
+M,134,Neighbor joining,16
+M,141,Neighbor joining,12
+M,152,Neighbor joining,17
+M,209,Neighbor joining,16
+M,45,Neighbor joining,13
+M,5,Neighbor joining,10
+M,50,Neighbor joining,10
+M,7,Neighbor joining,14
+M,73,Neighbor joining,15
+R,17,Neighbor joining,15
+R,2057,Neighbor joining,8
+R,2093,Neighbor joining,13
+R,2155,Neighbor joining,16
+R,2438,Neighbor joining,8
+R,2574,Neighbor joining,8
+R,2627,Neighbor joining,12
+R,2992,Neighbor joining,16
+R,3178,Neighbor joining,10
+R,325,Neighbor joining,7
+R,442,Neighbor joining,11

src/simulation/celltree.cpp

@@ -7,6 +7,7 @@
 
 #include "celltree.h"
 #include "beta_distribution.hpp"
+#include <lemon/bfs.h>
 
 CellTree::CellTree(double K,
                    double migrationRate,
@@ -1389,6 +1390,25 @@ void CellTree::constructSampledCloneTree()
     }
   }
 
+  // Find right order of arcs, reverse BFS order
+  IntNodeMap dist(_pCloneTT->tree(), -1);
+  lemon::bfs(_pCloneTT->tree()).distMap(dist).run(_pCloneTT->root());
+  int max_dist = lemon::mapMaxValue(_pCloneTT->tree(), dist);
+  ArcVector arcs;
+  while (max_dist > 0)
+  {
+    for (NodeIt v_i(_pCloneTT->tree()); v_i != lemon::INVALID; ++v_i)
+    {
+      if (dist[v_i] == max_dist)
+      {
+        Arc inArc = InArcIt(_pCloneTT->tree(), v_i);
+        assert(inArc != lemon::INVALID);
+        arcs.push_back(inArc);
+      }
+    }
+    --max_dist;
+  }
+  
   // construct vertex labeling
   for (ArcIt aa(_G); aa != lemon::INVALID; ++aa)
   {
@@ -1398,7 +1418,7 @@ void CellTree::constructSampledCloneTree()
     const int t = _anatomicalSiteMap[v_t];
 
     // find matching migration edge in T
-    for (ArcIt a(_pCloneTT->tree()); a != lemon::INVALID; ++a)
+    for (Arc a : arcs)
     {
       Node v_i = _pCloneTT->tree().source(a);
       Node v_j = _pCloneTT->tree().target(a);