src/AlignmentTree.cpp

Go to the documentation of this file.
00001 #ifdef HAVE_CONFIG_H
00002 #include "config.h"
00003 #endif
00004 
00005 #include "PhyloTree.h"
00006 #include <sstream>
00007 #include <stack>
00008 using namespace std;
00009 
00010 typedef unsigned uint;
00011 
00012 PhyloTree::PhyloTree() : vector< TreeNode >() {
00013         weight = 0;
00014         root = 0;
00015 }
00016 
00017 PhyloTree::PhyloTree( const PhyloTree& pt ) :
00018 vector< TreeNode >( pt ),
00019 weight( pt.weight ),
00020 root( pt.root )
00021 {}
00022 
00023 PhyloTree& PhyloTree::operator=( const PhyloTree& pt )
00024 {
00025         vector< TreeNode >::operator=( pt );
00026         weight = pt.weight;
00027         root = pt.root;
00028         return *this;
00029 }
00030 
00031 PhyloTree::~PhyloTree()
00032 {}
00033 
00034 void PhyloTree::clear()
00035 {
00036         vector< TreeNode >::clear();
00037         weight = 0;
00038         root = 0;
00039 }
00040 
00041 
00046 void PhyloTree::readTree( istream& tree_file ){
00047         string line;
00048         clear();
00049         if( !getline( tree_file, line ) )
00050                 return;
00051 
00052         stringstream line_str( line );
00053 
00054         // look for a weight
00055         string::size_type open_bracket_pos = line.find( "[" );
00056         string::size_type bracket_pos = line.find( "]" );
00057         if( open_bracket_pos != string::npos && bracket_pos != string::npos && 
00058                 open_bracket_pos < bracket_pos && bracket_pos < line.find( "(" ) ){
00059                 // read in a weight
00060                 getline( line_str, line, '[' );
00061                 getline( line_str, line, ']' );
00062                 stringstream weight_str( line );
00063                 weight_str >> weight;
00064         }
00065         
00066         // ready to begin parsing the tree data.
00067         string tree_line;
00068         getline( line_str, tree_line, ';' );
00069         uint read_state = 0;    
00070         uint section_start = 0;
00071         stack< node_id_t > node_stack;
00072         stringstream blen_str;
00073         TreeNode new_node;
00074         new_node.distance = 0;  // default the distance to 0
00075         for( uint charI = 0; charI < tree_line.size(); charI++ ){
00076                 switch( tree_line[ charI ] ){
00077                         // if this is an open parens then simply create a new
00078                         // parent node and push it on the parent stack
00079                         case '(':
00080                                 if( node_stack.size() > 0 ){
00081                                         new_node.parents.clear();
00082                                         new_node.parents.push_back( node_stack.top() );
00083                                         (*this)[ node_stack.top() ].children.push_back( (*this).size() );
00084                                 }
00085                                 node_stack.push( (*this).size() );
00086                                 push_back( new_node );
00087                                 read_state = 1;
00088                                 section_start = charI + 1;
00089                                 break;
00090                         case ')':
00091                                 // read off a branch length
00092                                 blen_str.clear();
00093                                 blen_str.str( tree_line.substr( section_start, charI - section_start ) );
00094                                 blen_str >> (*this)[ node_stack.top() ].distance;
00095                                 if( read_state == 2 )
00096                                         node_stack.pop();
00097                                 section_start = charI + 1;
00098                                 // pop off the top of the node stack after its branch length is read:
00099                                 read_state = 2;
00100                                 break;
00101                         case ',':
00102                                 // read off a branch length
00103                                 blen_str.clear();
00104                                 blen_str.str( tree_line.substr( section_start, charI - section_start ) );
00105                                 blen_str >> (*this)[ node_stack.top() ].distance;
00106                                 if( read_state == 2 )
00107                                         node_stack.pop();
00108                                 section_start = charI + 1;
00109                                 read_state = 1; // indicates that we'll be creating a new node when we hit :
00110                                 break;
00111                         case ':':
00112                                 // read off a name, if possible
00113                                 if( read_state == 1 ){
00114                                         new_node.parents.clear();
00115                                         new_node.parents.push_back( node_stack.top() );
00116                                         (*this)[ node_stack.top() ].children.push_back( (*this).size() );
00117                                         node_stack.push( (*this).size() );
00118                                         push_back( new_node );
00119                                         read_state = 2; // pop this node after reading its branch length
00120                                 }
00121                                 (*this)[ node_stack.top() ].name = tree_line.substr( section_start, charI - section_start );
00122                                 section_start = charI + 1;
00123                                 break;
00124                         default:
00125                                 break;
00126                 }
00127         }
00128 
00129 }
00130 
00131 
00132 void PhyloTree::writeTree( ostream& os ) const{
00133         stack< node_id_t > node_stack;
00134         stack< uint > child_stack;
00135         node_stack.push( root );
00136         child_stack.push( 0 );
00137 
00138         if( (*this).weight != 0 )
00139                 os << "[" << weight << "]";
00140         os << "(";
00141 
00142         while( node_stack.size() > 0 ) {
00143                 if( (*this)[ node_stack.top() ].children.size() != 0 ){
00144                         // this is a parent node
00145                         // if we have scanned all its children then pop it
00146                         if( child_stack.top() == (*this)[ node_stack.top() ].children.size() ){
00147                                 os << ")";
00148                                 if( node_stack.size() > 1 )
00149                                         os << ":" << (*this)[ node_stack.top() ].distance;
00150                                 node_stack.pop();
00151                                 child_stack.pop();
00152                                 continue;
00153                         }
00154                         // try to recurse to its children
00155                         // if the child is a parent as well spit out a paren
00156                         node_id_t child = (*this)[ node_stack.top() ].children[ child_stack.top() ];
00157                         node_stack.push( child );
00158                         child_stack.top()++;
00159                         // print a comma to separate multiple children
00160                         if( child_stack.top() > 1 )
00161                                 os << ",";
00162                         if( (*this)[ child ].children.size() > 0 ){
00163                                 child_stack.push( 0 );
00164                                 os << "(";
00165                         }
00166                         continue;
00167                 }
00168                 
00169                 // this is a leaf node
00170                 os << (*this)[ node_stack.top() ].name << ":" << (*this)[ node_stack.top() ].distance;
00171                 
00172                 // pop the child
00173                 node_stack.pop();
00174         }
00175         os << ";" << endl;
00176 }
00177 
00178 
00179 double PhyloTree::getHeight() const
00180 {
00181         return getHeight( root );
00182 }
00183 double PhyloTree::getHeight( node_id_t nodeI ) const
00184 {
00185         if( (*this)[ nodeI ].children.size() == 0 )
00186                 return (*this)[ nodeI ].distance;
00187         return (*this)[ nodeI ].distance + getHeight( (*this)[ nodeI ].children[ 0 ] );
00188 }

Generated on Mon Aug 19 06:00:40 2013 for mauveAligner by doxygen 1.3.6