Allow branch length updates

ComparativeGenomicsToolkit · Oct 1, 2013 · 26e7923 · 26e7923
1 parent 0483cc1
commit 26e7923
Show file tree

Hide file tree

Showing 4 changed files with 38 additions and 2 deletions.
diff --git a/api/hdf5_impl/hdf5Alignment.cpp b/api/hdf5_impl/hdf5Alignment.cpp
@@ -402,6 +402,27 @@ string HDF5Alignment::getParentName(const string& name) const
   return stTree_getLabel(parent);
 }
 
+
+void HDF5Alignment::updateBranchLength(const string& parentName,
+                                       const string& childName,
+                                       double length)
+{
+  map<string, stTree*>::iterator findIt = _nodeMap.find(childName);
+  if (findIt == _nodeMap.end())
+  {
+    throw hal_exception(string("node ") + childName + " not found");
+  }
+  stTree* node = findIt->second;
+  stTree* parent = stTree_getParent(node);
+  if (parent == NULL || parentName != stTree_getLabel(parent))
+  {
+    throw hal_exception(string("edge ") + parentName + "--" + childName + 
+                        " not found");
+  }
+  stTree_setBranchLength(node, length);
+  _dirty = true;
+}
+
 double HDF5Alignment::getBranchLength(const string& parentName,
                                       const string& childName) const
 {

diff --git a/api/hdf5_impl/hdf5Alignment.h b/api/hdf5_impl/hdf5Alignment.h
@@ -55,6 +55,10 @@ class HDF5Alignment : public Alignment
 
    std::string getParentName(const std::string& name) const;
 
+   void updateBranchLength(const std::string& parentName,
+                           const std::string& childName,
+                           double length);
+
    double getBranchLength(const std::string& parentName,
                           const std::string& childName) const;   
 

diff --git a/api/inc/halAlignment.h b/api/inc/halAlignment.h
@@ -91,7 +91,14 @@ class Alignment
     * @param childName name of child genome */
    virtual double getBranchLength(const std::string& parentName,
                                   const std::string& childName) const = 0;
-
+
+   /** Change the branch length between two genomes in the phylogeny 
+    * @param parentName name of parent genome
+    * @param childName name of child genome */   
+   virtual void updateBranchLength(const std::string& parentName,
+                                   const std::string& childName,
+                                   double length) = 0;
+
    /** Get names of child genomes in the phylogeny (empty vector for leaves)
     * @param name Name of genome */
    virtual std::vector<std::string> 

diff --git a/api/tests/halAlignmentTest.cpp b/api/tests/halAlignmentTest.cpp
@@ -103,14 +103,18 @@ void AlignmentTestTrees::createCallBack(hal::AlignmentPtr alignment)
   alignment->addLeafGenome("Leaf2", "Root", 5.1);
   alignment->addLeafGenome("Leaf3", "Root", 6.1);
   alignment->addLeafGenome("Leaf4", "Root", 7.1);
+  alignment->updateBranchLength("Root", "Leaf1", 3.0);
+  alignment->updateBranchLength("Root", "Leaf2", 6.1);
+  alignment->updateBranchLength("Root", "Leaf2", 5.1);
 }
 
 void AlignmentTestTrees::checkCallBack(hal::AlignmentConstPtr alignment)
 {
   CuAssertTrue(_testCase, alignment->getRootName() == "NewRoot");
   CuAssertTrue(_testCase, alignment->getNewickTree() == 
-        "((Leaf:10,Leaf1:4.1,Leaf2:5.1,Leaf3:6.1,Leaf4:7.1)Root:15)NewRoot;");
+        "((Leaf:10,Leaf1:3,Leaf2:5.1,Leaf3:6.1,Leaf4:7.1)Root:15)NewRoot;");
   CuAssertTrue(_testCase, alignment->getBranchLength("Root", "Leaf") == 10.0);
+  CuAssertTrue(_testCase, alignment->getBranchLength("Root", "Leaf1") == 3.0);
   vector<string> children = alignment->getChildNames("Root");
   CuAssertTrue(_testCase, children.size() == 5);