libMems/PhyloTree.h

Go to the documentation of this file.

#ifdef HAVE_CONFIG_H
#include "config.h"
#endif

#ifndef __PhyloTree_h__
#define __PhyloTree_h__

#include <vector>
#include <string>
#include <iostream>
#include <sstream>
#include <stack>

//typedef unsigned int node_id_t;
typedef size_t node_id_t;
class TreeNode 
{
public:
        TreeNode() : distance(0) {};
        std::string name;       
        double distance;        
        std::vector< node_id_t > parents;       
        std::vector< node_id_t > children;      
};

template< class T >
class PhyloTree 
{
public:
        PhyloTree();
        PhyloTree( const PhyloTree<T>& pt );
        PhyloTree<T>& operator=( const PhyloTree<T>& pt );
        ~PhyloTree();
        double weight;  
        node_id_t root; 
        std::vector< T > nodes; 
        void clear();
        void readTree( std::istream& tree_file );
        void writeTree( std::ostream& os ) const;
        double getHeight() const;
        double getHeight( node_id_t nodeI ) const;

        T& operator[]( const unsigned i ){ return nodes[i]; }
        const T& operator[]( const unsigned i ) const{ return nodes[i]; }
        size_t size() const{ return nodes.size(); }
        void push_back( T& t ){ nodes.push_back(t); }
        T& back() { return nodes.back(); }
        const T& back() const{ return nodes.back(); }
        void resize( const unsigned s ){ nodes.resize(s); }


        void swap( PhyloTree<T>& other )
        {
                std::swap( weight, other.weight );
                std::swap( root, other.root );
                nodes.swap( other.nodes );
        }
protected:
};


template< class T >
PhyloTree<T>::PhyloTree()
{
        weight = 0;
        root = 0;
}

template< class T >
PhyloTree<T>::PhyloTree( const PhyloTree<T>& pt ) :
nodes( pt.nodes ),
weight( pt.weight ),
root( pt.root )
{}

template< class T >
PhyloTree<T>& PhyloTree<T>::operator=( const PhyloTree<T>& pt )
{
        nodes = pt.nodes;
        weight = pt.weight;
        root = pt.root;
        return *this;
}

template< class T >
PhyloTree<T>::~PhyloTree()
{}

template< class T >
void PhyloTree<T>::clear()
{
        nodes.clear();
        weight = 0;
        root = 0;
}


template< class T >
void PhyloTree<T>::readTree( std::istream& tree_file )
{
        std::string line;
        clear();
        if( !std::getline( tree_file, line ) )
                return;
        // look for either a ; or a matched number of parenthesis, if
        // not found then read another line
        while(true){
                int paren_count = 0;
                for( size_t charI = 0; charI < line.size(); charI++ )
                {
                        if( line[charI] == '(' )
                                paren_count++;
                        if( line[charI] == ')' )
                                paren_count--;
                }
                if( paren_count == 0 )
                        break;
                if( paren_count != 0 ){
                        std::string another_line;
                        if( !std::getline( tree_file, another_line ) )
                                return;
                        line += another_line;
                }
        }

        std::stringstream line_str( line );

        // look for a weight
        std::string::size_type open_bracket_pos = line.find( "[" );
        std::string::size_type bracket_pos = line.find( "]" );
        if( open_bracket_pos != std::string::npos && bracket_pos != std::string::npos && 
                open_bracket_pos < bracket_pos && bracket_pos < line.find( "(" ) ){
                // read in a weight
                getline( line_str, line, '[' );
                getline( line_str, line, ']' );
                std::stringstream weight_str( line );
                weight_str >> weight;
        }
        
        // ready to begin parsing the tree data.
        std::string tree_line;
        std::getline( line_str, tree_line, ';' );
        uint read_state = 0;    
        uint section_start = 0;
        std::stack< node_id_t > node_stack;
        std::stringstream blen_str;
        T new_node;
        new_node.distance = 0;  // default the distance to 0
        bool already_read_name = false;
        bool blen_found = false;
        for( uint charI = 0; charI < tree_line.size(); charI++ ){
                switch( tree_line[ charI ] ){
                        // if this is an open parens then simply create a new
                        // parent node and push it on the parent stack
                        case '(':
                                if( node_stack.size() > 0 ){
                                        new_node.parents.clear();
                                        new_node.parents.push_back( node_stack.top() );
                                        (*this)[ node_stack.top() ].children.push_back( (node_id_t)(*this).size() );
                                }
                                node_stack.push( (node_id_t)(*this).size() );
                                nodes.push_back( new_node );
                                read_state = 1;
                                section_start = charI + 1;
                                break;
                        case ')':
                                if( blen_found )
                                {
                                        // read off a branch length
                                        blen_str.clear();
                                        blen_str.str( tree_line.substr( section_start, charI - section_start ) );
                                        blen_str >> (*this)[ node_stack.top() ].distance;
                                }else{
                                        // read off a name, if possible
                                        if( read_state == 1 ){
                                                new_node.parents.clear();
                                                new_node.parents.push_back( node_stack.top() );
                                                (*this)[ node_stack.top() ].children.push_back( (node_id_t)(*this).size() );
                                                node_stack.push( (node_id_t)(*this).size() );
                                                nodes.push_back( new_node );
                                                read_state = 2; // pop this node after reading its branch length
                                        }
                                        (*this)[ node_stack.top() ].name = tree_line.substr( section_start, charI - section_start );
                                }
                                if( read_state == 2 )
                                        node_stack.pop();
                                section_start = charI + 1;
                                blen_found = false;

                                // pop off the top of the node stack
                                read_state = 2;
                                break;
                        case ',':
                                if( blen_found ){
                                        // read off a branch length
                                        blen_str.clear();
                                        blen_str.str( tree_line.substr( section_start, charI - section_start ) );
                                        blen_str >> (*this)[ node_stack.top() ].distance;
                                }else{
                                        // read off a name, if possible
                                        if( read_state == 1 ){
                                                new_node.parents.clear();
                                                new_node.parents.push_back( node_stack.top() );
                                                (*this)[ node_stack.top() ].children.push_back( (node_id_t)(*this).size() );
                                                node_stack.push( (node_id_t)(*this).size() );
                                                nodes.push_back( new_node );
                                                read_state = 2; // pop this node after reading its name
                                        }
                                        (*this)[ node_stack.top() ].name = tree_line.substr( section_start, charI - section_start );
                                }
                                if( read_state == 2 )
                                        node_stack.pop();
                                section_start = charI + 1;
                                read_state = 1; // indicates that we'll be creating a new node when we hit :
                                blen_found = false;
                                break;
                        case ':':
                                // read off a name, if possible
                                if( read_state == 1 ){
                                        new_node.parents.clear();
                                        new_node.parents.push_back( node_stack.top() );
                                        (*this)[ node_stack.top() ].children.push_back( (node_id_t)(*this).size() );
                                        node_stack.push( (node_id_t)(*this).size() );
                                        nodes.push_back( new_node );
                                        read_state = 2; // pop this node after reading its branch length
                                }
                                (*this)[ node_stack.top() ].name = tree_line.substr( section_start, charI - section_start );
                                section_start = charI + 1;
                                blen_found = true;
                                break;
                        default:
                                break;
                }
        }

}


template< class T >
void PhyloTree<T>::writeTree( std::ostream& os ) const{
        std::stack< node_id_t > node_stack;
        std::stack< uint > child_stack;
        node_stack.push( root );
        child_stack.push( 0 );
        bool write_branch_lengths = false;
        for( size_t nodeI = 0; nodeI < this->size(); nodeI++ )
        {
                if( (*this)[nodeI].distance != 0 )
                {
                        write_branch_lengths = true;
                        break;
                }
        }

        if( (*this).weight != 0 )
                os << "[" << weight << "]";
        os << "(";

        while( node_stack.size() > 0 ) {
                if( (*this)[ node_stack.top() ].children.size() != 0 ){
                        // this is a parent node
                        // if we have scanned all its children then pop it
                        if( child_stack.top() == (*this)[ node_stack.top() ].children.size() ){
                                os << ")";
                                if( node_stack.size() > 1 && write_branch_lengths )
                                        os << ":" << (*this)[ node_stack.top() ].distance;
                                node_stack.pop();
                                child_stack.pop();
                                continue;
                        }
                        // try to recurse to its children
                        // if the child is a parent as well spit out a paren
                        node_id_t child = (*this)[ node_stack.top() ].children[ child_stack.top() ];
                        node_stack.push( child );
                        child_stack.top()++;
                        // print a comma to separate multiple children
                        if( child_stack.top() > 1 )
                                os << ",";
                        if( (*this)[ child ].children.size() > 0 ){
                                child_stack.push( 0 );
                                os << "(";
                        }
                        continue;
                }
                
                // this is a leaf node
                os << (*this)[ node_stack.top() ].name;
                if( write_branch_lengths )
                        os << ":" << (*this)[ node_stack.top() ].distance;
                
                // pop the child
                node_stack.pop();
        }
        os << ";" << std::endl;
}


template< class T >
double PhyloTree<T>::getHeight() const
{
        return getHeight( root );
}

template< class T >
double PhyloTree<T>::getHeight( node_id_t nodeI ) const
{
        if( (*this)[ nodeI ].children.size() == 0 )
                return (*this)[ nodeI ].distance;
        return (*this)[ nodeI ].distance + getHeight( (*this)[ nodeI ].children[ 0 ] );
}


template< class TreeType >
void getDescendants( TreeType& alignment_tree, node_id_t node, std::vector< node_id_t >& descendants )
{
        // do a depth first search
        std::stack< node_id_t > node_stack;
        node_stack.push( node );
        descendants.clear();
        while( node_stack.size() > 0 )
        {
                node_id_t cur_node = node_stack.top();
                node_stack.pop();
                if( alignment_tree[cur_node].children.size() > 0 )
                {
                        node_stack.push(alignment_tree[cur_node].children[0]);
                        node_stack.push(alignment_tree[cur_node].children[1]);
                }
                descendants.push_back(cur_node);
        }
}


template< class TreeType >
void getLeaves( TreeType& tree, node_id_t node, std::vector< node_id_t >& leaves )
{
        // do a depth first search
        std::stack< node_id_t > node_stack;
        node_stack.push( node );
        leaves.clear();
        while( node_stack.size() > 0 )
        {
                node_id_t cur_node = node_stack.top();
                node_stack.pop();
                if( tree[cur_node].children.size() > 0 )
                {
                        node_stack.push(tree[cur_node].children[0]);
                        node_stack.push(tree[cur_node].children[1]);
                }else
                        leaves.push_back(cur_node);
        }
}

namespace std {

template< class T > inline
void swap( PhyloTree<T>& a, PhyloTree<T>& b )
{
        a.swap(b);
}

template<> inline void swap( PhyloTree<TreeNode>& a, PhyloTree<TreeNode>& b){ a.swap(b); }
}

#endif // __PhyloTree_h__

---