#ifndef PHYCPP_TREE_TREE_HPP
#define PHYCPP_TREE_TREE_HPP

#include "taxa.hpp"
#include "node.hpp"
#include "edge.hpp"
#include "valueaccessor.hpp"
#include <cassert>

namespace PhyCpp {

namespace Detail {

	template <typename TaxonT>
	class TreeBase  {

	public:

		virtual void linkTaxa(BasicTaxa<TaxonT> * taxonGroup) = 0;
		virtual void unlinkTaxa() = 0;

	};

} // namespace detail

template <typename TaxonT, typename NodeT, typename EdgeT>
class BasicTree : public Detail::TreeBase<TaxonT>, public SimpleLabelledObject<std::string> {

	friend class BasicNode<TaxonT, NodeT, EdgeT>;

public:

	typedef BasicTaxon<TaxonT> 					Taxon;
	typedef BasicTaxa<TaxonT> 					Taxa;
	typedef BasicNode<TaxonT, NodeT, EdgeT> 	Node;
	typedef BasicEdge<TaxonT, NodeT, EdgeT> 	Edge;

	/** @name Construction */
	//@{
	BasicTree(Taxa * taxonGroup);
	BasicTree(const std::string & newickString, Taxa * taxonGroup, WhiteSpaceBehaviour wsb = TrimWhiteSpace);
	BasicTree(const BasicTree & other);
	BasicTree(const BasicTree & other, Taxa * taxonGroup);
	BasicTree(const Node * subtreeRoot);
	BasicTree(const Node * subtreeRoot, Taxa * taxonGroup);
	~BasicTree();
	BasicTree & operator =(const BasicTree & other);
	//@}

	/** @name Taxa */
	//@{
	inline Taxa * taxa() {return linked_taxa;}
	inline const Taxa * taxa() const {return linked_taxa;}
	inline unsigned int taxonCount() const {return taxonMap.size();}
	void linkTaxa(Taxa * taxonGroup);
	void unlinkTaxa();
	//@}

	/** @name Tree metrics*/
	//@{
	inline unsigned int size() const {return node_count;}
	double length() const;
	//@}

	/** @name Topology access */
	//@{
	Node * node(const std::string & taxonLabel);
	const Node * node(const std::string & taxonLabel) const;
	inline Node * root() {return (startNode->preordernext == endNode) ? 0 : startNode->preordernext;}
	inline const Node * root() const {return (startNode->preordernext == endNode) ? 0 : startNode->preordernext;}
	inline std::string toString() const {return (startNode->preordernext == 0) ? "" : startNode->preordernext->toString();}
	//@}

//	/** @name Topology modification */
//	//@{
	void clear();
	void removeNode(Node * n, bool retainUnbranchedParent = false, bool accumulateBranchLengths = true);
	void removeSubtree(Node * subtreeRoot, bool retainUnbranchedParent = false, bool accumulateBranchLengths = true);
	void resolvePolytomies(PolytomyResolutionMethod method);
	void setRoot(Node * n, bool retainOldUnbranchedRoot = false, bool accumulateBranchLengths = true);
	//@}

	/** @name Traversal */
	//@{
	typename Node::iterator begin(Node * subtreeRoot = 0);
	typename Node::iterator end(Node * subtreeRoot = 0);
	typename Node::const_iterator begin(const Node * subtreeRoot = 0) const;
	typename Node::const_iterator end(const Node * subtreeRoot = 0) const;
	typename Node::reverse_iterator rbegin(Node * subtreeRoot = 0);
	typename Node::reverse_iterator rend(Node * subtreeRoot = 0);
	typename Node::const_reverse_iterator rbegin(const Node * subtreeRoot = 0) const;
	typename Node::const_reverse_iterator rend(const Node * subtreeRoot = 0) const;
	//@}

private:

	void removeInternalEdgeNode(Node * n, bool accumulateBranchLengths);
	void copySubtree(Node * me, const Node * other);
	void processNewickToken(Node * n, std::string & token, WhiteSpaceBehaviour wsb) const;

	Taxon * setTaxon(Node * n, Taxon * t);
	Taxon * setTaxon(Node * n, const std::string & label);

	Taxa * linked_taxa;
	Node * startNode;
	Node * endNode;
	unsigned int node_count;
	std::map<Taxon*, Node*> taxonMap;

};

} // namespace phycpp

template <typename TaxonT, typename NodeT, typename EdgeT>
PhyCpp::BasicTree<TaxonT, NodeT, EdgeT>::BasicTree(Taxa * taxonGroup)
: SimpleLabelledObject<std::string>(),
  linked_taxa(0),
  startNode(new Node(this, -1)),
  endNode(new Node(this, -1)),
  node_count(0) {
	assert(taxonGroup != 0);
	linkTaxa(taxonGroup);
}

template <typename TaxonT, typename NodeT, typename EdgeT>
PhyCpp::BasicTree<TaxonT, NodeT, EdgeT>::BasicTree(const std::string & newickString, Taxa * taxonGroup, WhiteSpaceBehaviour wsb)
: SimpleLabelledObject<std::string>(),
  linked_taxa(0),
  startNode(new Node(this, -1)),
  endNode(new Node(this, -1)),
  node_count(0) {
	assert(taxonGroup != 0);
	linkTaxa(taxonGroup);

	std::string::const_iterator it = newickString.begin();
	std::string::const_iterator ite = newickString.end();

	Node * root_node = new Node(this, -1);
	startNode->preordernext = root_node;
	root_node->postordernext = startNode;

	Node * activeNode = 0;
	std::string label;
	char lastToken = ' ';
	while (true) {
		if (*it == '[') {
			while (it != ite && *it != ']')
				++it;
		}
		else if (*it == ';' || it == ite) {
			if (activeNode != 0)
				processNewickToken(activeNode, label, wsb);
			break;
		}
		else if (*it == '(') {
			if (activeNode == 0) {
				activeNode = root_node;
			} else {
				activeNode = activeNode->appendChild();
			}
		} else if (*it == ')') {
			Node * n;
			if (lastToken != ')')
				n = activeNode->appendChild();
			else {
				n = activeNode->leftchild;
				while (n->rightsibling != 0)
					n = n->rightsibling;
			}
			processNewickToken(n, label, wsb);
			if (activeNode != root_node)
				activeNode = activeNode->parentnode;
		} else if (*it == ',') {
			Node * n;
			if (lastToken != ')')
				n = activeNode->appendChild();
			else {
				n = activeNode->leftchild;
				while (n->rightsibling != 0)
					n = n->rightsibling;
			}
			processNewickToken(n, label, wsb);
		} else {
			label += *it;
		}
		if (*it == ')' || *it == '(' || *it == ',')
			lastToken = *it;
		++it;
	};
}

template <typename TaxonT, typename NodeT, typename EdgeT>
void PhyCpp::BasicTree<TaxonT, NodeT, EdgeT>::processNewickToken(Node * n, std::string & token, WhiteSpaceBehaviour wsb) const {

	//todo : pos = token.size()-1 (i.e. "A:")  currently set len as 0
	size_t pos = token.find(':');
	if (pos == std::string::npos) {
		if (wsb == TrimWhiteSpace)
			n->setTaxon(trimWhiteSpace(token));
		else if (wsb == StripWhiteSpace)
			n->setTaxon(stripWhiteSpace(token));
		else
			n->setTaxon(token);
	} else if (pos > 0) {
		if (wsb == TrimWhiteSpace)
			n->setTaxon(trimWhiteSpace(token.substr(0, pos)));
		else if (wsb == StripWhiteSpace)
			n->setTaxon(stripWhiteSpace(token.substr(0, pos)));
		else
			n->setTaxon(token.substr(0, pos));
		n->parentedge->setLength(strtod(token.substr(pos+1, token.length() - pos+1).data(), 0));
	} else {
		n->parentedge->setLength(strtod(token.substr(1, token.length() - 1).data(), 0));
	}
	token = "";
}

template <typename TaxonT, typename NodeT, typename EdgeT>
void PhyCpp::BasicTree<TaxonT, NodeT, EdgeT>::copySubtree(BasicNode<TaxonT, NodeT, EdgeT> * me, const BasicNode<TaxonT, NodeT, EdgeT> * other) {
	if (other == 0)
		return;
	assert(me->leftchild == 0);
	if (other->leftchild != 0) {
		Node * k = other->leftchild;
		Node * j = me->appendChild(k->parentedge->length());
		if (k->tax != 0)
			j->setTaxon(k->tax->label());
		copySubtree(j, k);
		while (k->rightsibling != 0) {
			k = k->rightsibling;
			j = me->appendChild(k->parentedge->length());
			if (k->tax != 0)
				j->setTaxon(k->tax->label());
			copySubtree(j, k);
		}
	}
}

template <typename TaxonT, typename NodeT, typename EdgeT>
PhyCpp::BasicTree<TaxonT, NodeT, EdgeT>::BasicTree(const BasicTree<TaxonT, NodeT, EdgeT> & other)
: SimpleLabelledObject<std::string>(other.label()),
  linked_taxa(0),
  startNode(new Node(this, -1)),
  endNode(new Node(this, -1)),
  node_count(0) {
	linkTaxa(other.linked_taxa);
	if (other.node_count > 0) {
		Node * m = new Node(this, -1);
		startNode->preordernext = m;
		m->postordernext = startNode;
		Node * o = other.startNode->preordernext;
		m->parentedge->setLength(o->parentedge->length());
		if (o->tax != 0)
		m->setTaxon(o->tax->label());
		copySubtree(m, o);
	}
}

template <typename TaxonT, typename NodeT, typename EdgeT>
PhyCpp::BasicTree<TaxonT, NodeT, EdgeT>::BasicTree(const BasicTree<TaxonT, NodeT, EdgeT> & other, BasicTaxa<TaxonT> * taxonGroup)
: SimpleLabelledObject<std::string>(other.label()),
  linked_taxa(0),
  startNode(new Node(this, -1)),
  endNode(new Node(this, -1)),
  node_count(0) {
	assert(taxonGroup != 0);
	linkTaxa(taxonGroup);
	if (other.node_count > 0) {
		Node * m = new Node(this, -1);
		startNode->preordernext = m;
		m->postordernext = startNode;
		Node * o = other.startNode->preordernext;
		m->parentedge->setLength(o->parentedge->length());
		if (o->tax != 0)
		m->setTaxon(o->tax->label());
		copySubtree(m, o);
	}
}

template <typename TaxonT, typename NodeT, typename EdgeT>
PhyCpp::BasicTree<TaxonT, NodeT, EdgeT>::BasicTree(const BasicNode<TaxonT, NodeT, EdgeT> * subtreeRoot)
: SimpleLabelledObject<std::string>(),
  linked_taxa(0),
  startNode(new Node(this, -1)),
  endNode(new Node(this, -1)),
  node_count(0) {
	linkTaxa(subtreeRoot->tree_ptr->linked_taxa);
	Node * m = new Node(this, -1);
	startNode->preordernext = m;
	m->postordernext = startNode;
	m->parentedge->setLength(subtreeRoot->parentedge->length());
	if (subtreeRoot->tax != 0)
		m->setTaxon(subtreeRoot->tax->label());
	copySubtree(m, subtreeRoot);
}

template <typename TaxonT, typename NodeT, typename EdgeT>
PhyCpp::BasicTree<TaxonT, NodeT, EdgeT>::BasicTree(const BasicNode<TaxonT, NodeT, EdgeT> * subtreeRoot, Taxa * taxonGroup)
: SimpleLabelledObject<std::string>(),
  linked_taxa(0),
  startNode(new Node(this, -1)),
  endNode(new Node(this, -1)),
  node_count(0) {
	assert(taxonGroup != 0);
	linkTaxa(taxonGroup);
	Node * m = new Node(this, -1);
	startNode->preordernext = m;
	m->postordernext = startNode;
	m->parentedge->setLength(subtreeRoot->parentedge->length());
	if (subtreeRoot->tax != 0)
		m->setTaxon(subtreeRoot->tax->label());
	copySubtree(m, subtreeRoot);
}

template <typename TaxonT, typename NodeT, typename EdgeT>
PhyCpp::BasicTree<TaxonT, NodeT, EdgeT>::~BasicTree() {
	clear();
	unlinkTaxa();
	delete startNode;
	delete endNode;
}

template <typename TaxonT, typename NodeT, typename EdgeT>
PhyCpp::BasicTree<TaxonT, NodeT, EdgeT> & PhyCpp::BasicTree<TaxonT, NodeT, EdgeT>::operator =(const BasicTree & other) {
	clear();
	setLabel(other.label());
	if (other.node_count > 0) {
		Node * m = new Node(this, -1);
		startNode->preordernext = m;
		m->postordernext = startNode;
		Node * o = other.startNode->preordernext;
		m->parentedge->setLength(o->parentedge->length());
		if (o->tax != 0)
		m->setTaxon(o->tax->label());
		copySubtree(m, o);
	}
	return *this;
}

template <typename TaxonT, typename NodeT, typename EdgeT>
void PhyCpp::BasicTree<TaxonT, NodeT, EdgeT>::clear() {
	if (startNode->preordernext != 0) {
		startNode->preordernext->destroy();
		delete startNode->preordernext;
		startNode->preordernext = endNode->postordernext = 0;
		taxonMap.clear();
	}
}

template <typename TaxonT, typename NodeT, typename EdgeT>
PhyCpp::BasicTaxon<TaxonT> * PhyCpp::BasicTree<TaxonT, NodeT, EdgeT>::setTaxon(BasicNode<TaxonT, NodeT, EdgeT> * n, BasicTaxon<TaxonT> * t) {
	if (n->tree_ptr != this)
		return 0;
	if (t != 0) {
		std::pair<typename std::map<BasicTaxon<TaxonT>*, BasicNode<TaxonT, NodeT, EdgeT>* >::iterator, bool> z = taxonMap.insert(std::make_pair(t, n));
		if (z.second) {
			if (n->tax != 0)
				taxonMap.erase(taxonMap.find(n->tax));
				n->tax = t;
				return t;
		} else {
			return 0;
		}
	} else {
		taxonMap.erase(n->tax);
		n->tax = 0;
		return 0;
	}
}

template <typename TaxonT, typename NodeT, typename EdgeT>
PhyCpp::BasicTaxon<TaxonT> * PhyCpp::BasicTree<TaxonT, NodeT, EdgeT>::setTaxon(BasicNode<TaxonT, NodeT, EdgeT> * n, const std::string & label) {
	if (n->tree_ptr != this)
		return 0;
	Taxon * t = linked_taxa->insert(label);
	return (t != 0) ? setTaxon(n, t) : 0;
}

template <typename TaxonT, typename NodeT, typename EdgeT>
PhyCpp::BasicNode<TaxonT, NodeT, EdgeT> * PhyCpp::BasicTree<TaxonT, NodeT, EdgeT>::node(const std::string & taxonLabel) {
	typename std::map<Taxon*, Node*>::iterator it = taxonMap.begin();
	while (it != taxonMap.end()) {
		if (it->first->label() == taxonLabel)
			return it->second;
		++it;
	}
	return 0;
}

template <typename TaxonT, typename NodeT, typename EdgeT>
const PhyCpp::BasicNode<TaxonT, NodeT, EdgeT> * PhyCpp::BasicTree<TaxonT, NodeT, EdgeT>::node(const std::string & taxonLabel) const {
	typename std::map<Taxon*, Node*>::iterator it = taxonMap.begin();
	while (it != taxonMap.end()) {
		if (it->first->label() == taxonLabel)
			return it->second;
		++it;
	}
	return 0;
}

template <typename TaxonT, typename NodeT, typename EdgeT>
void PhyCpp::BasicTree<TaxonT, NodeT, EdgeT>::linkTaxa(BasicTaxa<TaxonT> * taxonGroup) {
	if (linked_taxa == taxonGroup)
		return;
	if (taxonGroup == 0)
		unlinkTaxa();
	else if (linked_taxa != 0) {
		std::map<Node*, std::string> oldconnections;
		typename std::map<Taxon*, Node*>::iterator it = taxonMap.begin();
		while (it != taxonMap.end()) {
			oldconnections[it->second] = it->first->label();
			++it;
		}
		unlinkTaxa();
		linked_taxa = taxonGroup;
		taxonGroup->trees.push_back(this);
		typename std::map<Node*, std::string>::iterator k = oldconnections.begin();
		while (k != oldconnections.end()) {
			k->first->setTaxon(k->second);
			++k;
		}
	} else {
		linked_taxa = taxonGroup;
		taxonGroup->trees.push_back(this);
	}
}

template <typename TaxonT, typename NodeT, typename EdgeT>
void PhyCpp::BasicTree<TaxonT, NodeT, EdgeT>::unlinkTaxa() {
	if (linked_taxa == 0)
		return;
	typename std::vector<Detail::TreeBase<TaxonT>*>::iterator it = std::find(linked_taxa->trees.begin(), linked_taxa->trees.end(), this);
	if (it != linked_taxa->trees.end()) {
		linked_taxa->trees.erase(it);
		linked_taxa = 0;
		typename std::map<Taxon*, Node*>::iterator k = taxonMap.begin();
		while (k != taxonMap.end()) {
			k->second->tax = 0;
			++k;
		}
		taxonMap.clear();
	}
}

template <typename TaxonT, typename NodeT, typename EdgeT>
typename PhyCpp::BasicNode<TaxonT, NodeT, EdgeT>::iterator PhyCpp::BasicTree<TaxonT, NodeT, EdgeT>::begin(Node * subtreeRoot) {
	return (subtreeRoot == 0 || subtreeRoot == startNode->preordernext) ? startNode->preordernext : subtreeRoot;
}

template <typename TaxonT, typename NodeT, typename EdgeT>
typename PhyCpp::BasicNode<TaxonT, NodeT, EdgeT>::iterator PhyCpp::BasicTree<TaxonT, NodeT, EdgeT>::end(Node * subtreeRoot) {
	if (subtreeRoot == 0 || subtreeRoot == startNode->preordernext) {
		return endNode;
	} else {
		while (subtreeRoot != 0 && subtreeRoot->rightsibling == 0)
			subtreeRoot = subtreeRoot->parentnode;
		return (subtreeRoot == 0) ? endNode : subtreeRoot->rightsibling;
	}
}

template <typename TaxonT, typename NodeT, typename EdgeT>
typename PhyCpp::BasicNode<TaxonT, NodeT, EdgeT>::const_iterator PhyCpp::BasicTree<TaxonT, NodeT, EdgeT>::begin(const Node * subtreeRoot) const {
	return (subtreeRoot == 0 || subtreeRoot == startNode->preordernext) ? startNode->preordernext : subtreeRoot;
}

template <typename TaxonT, typename NodeT, typename EdgeT>
typename PhyCpp::BasicNode<TaxonT, NodeT, EdgeT>::const_iterator PhyCpp::BasicTree<TaxonT, NodeT, EdgeT>::end(const Node * subtreeRoot) const {
	if (subtreeRoot == 0 || subtreeRoot == startNode->preordernext) {
		return endNode;
	} else {
		while (subtreeRoot != 0 && subtreeRoot->rightsibling == 0)
			subtreeRoot = subtreeRoot->parentnode;
		return (subtreeRoot == 0) ? endNode : subtreeRoot->rightsibling;
	}
}

template <typename TaxonT, typename NodeT, typename EdgeT>
typename PhyCpp::BasicNode<TaxonT, NodeT, EdgeT>::reverse_iterator PhyCpp::BasicTree<TaxonT, NodeT, EdgeT>::rbegin(Node * subtreeRoot) {
	return (subtreeRoot == 0 || subtreeRoot == startNode->preordernext || *end(subtreeRoot) == endNode) ? endNode->postordernext : end(subtreeRoot)->postordernext;
}

template <typename TaxonT, typename NodeT, typename EdgeT>
typename PhyCpp::BasicNode<TaxonT, NodeT, EdgeT>::reverse_iterator PhyCpp::BasicTree<TaxonT, NodeT, EdgeT>::rend(Node * subtreeRoot) {
	return (subtreeRoot == 0 || subtreeRoot == startNode->preordernext) ? startNode : subtreeRoot->postordernext;
}

template <typename TaxonT, typename NodeT, typename EdgeT>
typename PhyCpp::BasicNode<TaxonT, NodeT, EdgeT>::const_reverse_iterator PhyCpp::BasicTree<TaxonT, NodeT, EdgeT>::rbegin(const Node * subtreeRoot) const {
	return (subtreeRoot == 0 || subtreeRoot == startNode->preordernext || *end(subtreeRoot) == endNode) ? endNode->postordernext : end(subtreeRoot)->postordernext;
}

template <typename TaxonT, typename NodeT, typename EdgeT>
typename PhyCpp::BasicNode<TaxonT, NodeT, EdgeT>::const_reverse_iterator PhyCpp::BasicTree<TaxonT, NodeT, EdgeT>::rend(const Node * subtreeRoot) const {
	return (subtreeRoot == 0 || subtreeRoot == startNode->preordernext) ? startNode : subtreeRoot->postordernext;
}

template <typename TaxonT, typename NodeT, typename EdgeT>
void PhyCpp::BasicTree<TaxonT, NodeT, EdgeT>::removeSubtree(Node * subtreeRoot, bool retainUnbranchedParent, bool accumulateBranchLengths) {
	if (subtreeRoot == startNode->preordernext) {
		clear();
		return;
	} else {
		Node * n = subtreeRoot->parentnode;
		if (!retainUnbranchedParent && n == startNode->preordernext && (n->leftchild->rightsibling != 0 && n->leftchild->rightsibling->rightsibling == 0)) {
			Node::setPreorderRelation(0, n->leftchild);
			bool recalc = (subtreeRoot != n->leftchild);
			if (!recalc)
				Node::setPreorderRelation(startNode, n->leftchild->rightsibling);
			subtreeRoot->destroy();
			delete subtreeRoot;
			if (recalc)
				Node::setPreorderRelation(n->rightSubtreeExtent(), 0);
			delete n;
			n = startNode->preordernext;
		} else if (n == startNode->preordernext && startNode->preordernext->leftchild->rightsibling == 0) {
			clear();
			return;
		} else {
			bool recalc = (n->leftchild != subtreeRoot);
			if (!recalc) {
				n->leftchild = subtreeRoot->rightsibling;
				Node::setPreorderRelation(n, subtreeRoot->subtreePreorderNext());
			}
			else {
				Node * L = subtreeRoot->leftSibling();
				L->rightsibling = subtreeRoot->rightsibling;
				Node::setPreorderRelation(L->rightSubtreeExtent(), subtreeRoot->subtreePreorderNext());
			}
			subtreeRoot->destroy();
			delete subtreeRoot;
		}
		if (!retainUnbranchedParent && n->leftchild == 0)
			removeSubtree(n);
		 if (n->leftchild != 0 && n->leftchild->rightsibling == 0 && !retainUnbranchedParent)
			removeNode(n, retainUnbranchedParent, accumulateBranchLengths);
	}
}

template <typename TaxonT, typename NodeT, typename EdgeT>
void PhyCpp::BasicTree<TaxonT, NodeT, EdgeT>::removeNode(Node * n, bool retainUnbranchedParent, bool accumulateBranchLengths) {
	if (n->leftchild == 0) {																			// a tip
		removeSubtree(n, retainUnbranchedParent);
	} else if (n == startNode->preordernext) {															// the root
		if (n->leftchild == 0)
			clear();
		else
			setRoot(n->leftchild, false, accumulateBranchLengths);
	} else if (n->leftchild != 0 && n->leftchild->rightsibling == 0 && n->parentnode != 0) {			// an unbranched internal node
		if (accumulateBranchLengths) {
			if (n->leftchild->parentedge->length() >= 0 && n->parentedge->length() >= 0)
				n->leftchild->parentedge->setLength(n->leftchild->parentedge->length() + n->parentedge->length());
			else if (n->parentedge->length() >= 0)
				n->leftchild->parentedge->setLength(n->parentedge->length());
		}
		n->leftchild->parentnode = n->parentnode;
		if (n->parentnode->leftchild == n) {
			n->parentnode->leftchild = n->leftchild;
			n->leftchild->rightsibling = n->rightsibling;
			n->leftchild->parentnode = n->parentnode;
		} else {
			n->leftSibling()->rightsibling = n->leftchild;
			n->leftchild->parentnode = n->parentnode;
			n->leftchild->rightsibling = n->rightsibling;
		}
		n->postordernext->preordernext = n->preordernext;
		n->preordernext->postordernext = n->postordernext;
		if (n->postordernext == 0) {
			n->preordernext->preordernext = endNode;
			endNode->postordernext = n->preordernext->preordernext;
		}
		Node * k = n->parentnode;
		delete n;
		if (!retainUnbranchedParent && k != 0 && k->leftchild != 0 && k->leftchild->rightsibling == 0 && k->parentnode != 0)
			removeNode(n->parentnode);
	} else {																							// a branched internal node
		Node * parent = n->parentnode;
		Node * leftsib = n->leftSibling();
		Node * rightsib = n->rightsibling;
		Node * leftchild = n->leftchild;
		double len = n->parentedge->length();
		n->parentnode = n->rightsibling = 0;
		delete n;
		if (leftsib) {
			leftsib->rightsibling = leftchild;
			Node::setPreorderRelation(leftsib->rightSubtreeExtent(), leftchild);
		} else {
			parent->leftchild = leftchild;
			Node::setPreorderRelation(parent, leftchild);
		}
		n = leftchild;
		while (n->rightsibling) {
			n->parentnode = parent;
			if (accumulateBranchLengths)
				n->parentedge->setLength(n->parentedge->length() + len);
			n = n->rightsibling;
		}
		n->parentnode = parent;
		if (accumulateBranchLengths)
			n->parentedge->setLength(n->parentedge->length() + len);
		if (rightsib) {
			n->rightsibling = rightsib;
			Node::setPreorderRelation(n->rightSubtreeExtent(), rightsib);
		} else {
			Node::setPreorderRelation(n->rightSubtreeExtent(), n->subtreePreorderNext());
		}
	}
}

template <typename TaxonT, typename NodeT, typename EdgeT>
void PhyCpp::BasicTree<TaxonT, NodeT, EdgeT>::setRoot(Node * n, bool retainOldUnbranchedRoot, bool accumulateBranchLengths) {
	if (n == 0 || n->tree_ptr != this || node_count == 0) return;
	Node * k = startNode->preordernext;
	n->reroot((Node*)0);
	n->recalculateTraversalOrder();
	if (!retainOldUnbranchedRoot && k->parentnode != 0 && k->leftchild != 0 && k->leftchild->rightsibling == 0) {
		this->removeNode(k, retainOldUnbranchedRoot, accumulateBranchLengths);
	}
}

template <typename TaxonT, typename NodeT, typename EdgeT>
void PhyCpp::BasicTree<TaxonT, NodeT, EdgeT>::resolvePolytomies(PolytomyResolutionMethod method) {
	typename Node::iterator it = begin();
	while (it != end()) {
		it.node()->resolvePolytomy(method);
		++it;
	}
}

template <typename TaxonT, typename NodeT, typename EdgeT>
double PhyCpp::BasicTree<TaxonT, NodeT, EdgeT>::length() const {
	typename Node::const_iterator it = begin();
	double result = 0;
	while (it != end()) {
		if (it->parentEdge()->length() > 0)
			result += it->parentEdge()->length();
		++it;
	}
	return result;
}

#endif /* PHYCPP_TREE_TREE_HPP */

