Go back to Richel Bilderbeek's homepage.

Go back to Richel Bilderbeek's C++ page.

 

 

 

 

 

(C++) Newick

 

A Newick is a way to write down a phylogeny as a std::string. This page shows how to check this std::string and how to store a Newick more efficiently.

 

Because personally, I only work with Newicks in the form '(((A,B),(C,D)),E)', these algorithms will so as well.

 

 

 

 

 

Newick code snippets

 

  1. CalcComplexity
  2. CalcNumOfCombinations
  3. CalcNumOfSymmetries
  4. CheckNewick
  5. CreateInvalidNewicks
  6. CreateValidNewicks
  7. DumbNewickToString
  8. GetRootBranches
  9. GetSimplerNewicks
  10. InspectInvalidNewick
  11. IsBinaryNewick
  12. IsNewick
  13. NewickToString
  14. BinaryNewickVector
  15. StringToNewick

 

 

 

 

 

newick.h

 

//---------------------------------------------------------------------------
/*
Newick, Newick functions
Copyright (C) 2010-2011 Richel Bilderbeek

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program.If not, see <http://www.gnu.org/licenses/>.
*/
//---------------------------------------------------------------------------
//From http://www.richelbilderbeek.nl/CppNewick.htm
//---------------------------------------------------------------------------
#ifndef NEWICK_H
#define NEWICK_H
//---------------------------------------------------------------------------
#ifdef _WIN32
//See http://www.richelbilderbeek.nl/CppCompileErrorSwprintfHasNotBeenDeclared.htm
#undef __STRICT_ANSI__
#endif
//---------------------------------------------------------------------------
#include <cmath>
#include <string>
#include <vector>
#include <boost/foreach.hpp>
#include <boost/lexical_cast.hpp>
#include <boost/tuple/tuple.hpp>
//---------------------------------------------------------------------------
#include "BigIntegerLibrary.hh"
#include "newickcpp98.h"
#include "newickstorage.h"
#include "trace.h"
//---------------------------------------------------------------------------
///namespace Newick contains general Newick functions,
///not using an Newick class
namespace Newick
{
  enum { bracket_open  = -1 };
  enum { bracket_close = -2 };
  enum { comma         = -3 };
  enum { new_line      = -4 };
  enum { null          = -5 };

  ///fuzzy_equal_to is a predicate to test two doubles for equality
  ///with a certain tolerance. A tolerance of 0.0 denotes that
  ///an exact match is requested. Note that the value of 0.0 cannot
  ///be compared fuzzily.
  //From http://www.richelbilderbeek.nl/CppFuzzy_equal_to.htm
  struct fuzzy_equal_to : public std::binary_function<double,double,bool>
  {
    fuzzy_equal_to(const double tolerance = 0.01)
      : m_tolerance(tolerance)
    {
      assert(tolerance >= 0.0);
    }
    bool operator()(const double lhs, const double rhs) const
    {
      const double d = std::fabs(m_tolerance * lhs);
      return rhs > lhs - d
          && rhs < lhs + d;
    }
    const double m_tolerance;
  };
  //---------------------------------------------------------------------------
  ///CreateVector creates a std::vector from three arguments
  ///From http://www.richelbilderbeek.nl/CppCreateVector.htm

  template <class T>
  const std::vector<T> CreateVector(const T& a, const T& b, const T& c)
  {
    std::vector<T> v;
    v.reserve(3);
    v.push_back(a);
    v.push_back(b);
    v.push_back(c);
    return v;
  }
  //---------------------------------------------------------------------------
  bool AllAboutEqual(const std::vector<double>& v,const double tolerance = 0.01);
  const BigInteger CalcComplexity(const std::vector<int>& v);
  const BigInteger CalcNumOfCombinationsBinary(const std::vector<int>& v);
  const BigInteger CalcNumOfSymmetriesBinary(std::vector<int> v);
  double CalcDenominator(const std::vector<int>& v,const double theta);
  double CalcProbabilitySimpleNewick(const std::vector<int>& v,const double theta);
  void CheckNewick(const std::string& s);
  void CheckNewick(const std::vector<int>& v);
  const std::vector<std::string> CreateInvalidNewicks();
  const std::string CreateRandomNewick(const int n,const int max);
  const std::vector<int> CreateRandomBinaryNewickVector(const int n,const int max);
  const std::vector<std::string> CreateValidBinaryNewicks();
  const std::vector<std::string> CreateValidNewicks();
  const std::vector<std::string> CreateValidTrinaryNewicks();
  const std::vector<std::string> CreateValidUnaryNewicks();
  const std::string DumbNewickToString(const std::vector<int>& v);
  const std::vector<int> Factorial(const std::vector<int>& v_original);
  const BigInteger FactorialBigInt(const int n);
  int Factorial(const int n);
  int FindPosAfter(const std::vector<int>& v,const int index,const int value);
  int FindPosBefore(const std::vector<int>& v,const int index,const int value);
  const std::vector<int> GetDepth(const std::vector<int>& n);
  const std::vector<int> GetFactorialTerms(const int n);
  std::vector<boost::tuple<std::string,double,double> > GetKnownProbabilities();
  int GetLeafMaxArity(const std::vector<int>& n);
  const std::vector<std::vector<int> >
    GetRootBranches(const std::vector<int>& n);
  const std::pair<std::vector<int>,std::vector<int> >
    GetRootBranchesBinary(const std::vector<int>& n);
  const std::vector<std::vector<int> > GetSimplerBinaryNewicks(
    const std::vector<int>& n);
  const std::vector<std::vector<int> > GetSimplerNewicks(
    const std::vector<int>& n);
  const std::vector<std::pair<std::vector<int>,int> >
    GetSimplerNewicksFrequencyPairs(
    const std::vector<int>& n);
  const std::vector<std::pair<std::vector<int>,int> >
    GetSimplerBinaryNewicksFrequencyPairs(
    const std::vector<int>& n);
  void InspectInvalidNewick(std::ostream& os, const std::vector<int>& v);
  bool IsBinaryNewick(std::vector<int> v);
  bool IsTrinaryNewick(std::vector<int> v);
  bool IsUnaryNewick(const std::vector<int>& v);
  bool IsNewick(const std::string& s);
  bool IsNewick(const std::vector<int>& v);
  bool IsSimple(const std::vector<int>& v);
  const std::string NewickToString(const std::vector<int>& v);
  const std::vector<int> ReplaceLeave(const std::vector<int>& newick,const int value);
  const std::vector<int> StringToNewick(const std::string& newick);
  const std::vector<int> Surround(const std::vector<int>& newick);
  void Test();
  char ValueToChar(const int value);

  template <class NewickType>
  double CalculateProbability(
    const NewickType& n,
    const double theta,
    NewickStorage<NewickType>& storage)
{
  //#define TRACE_NEWICK_CALCULATEPROBABILITY
  while(1)
  {
    try
    {
      //Is n already known?
      {
        const double p = storage.Find(n);
        if (p!=0.0)
        {
          return p;
        }
      }

      //Check for simple phylogeny
      if (n.IsSimple())
      {
        const double p = n.CalcProbabilitySimpleNewick(theta);
        storage.Store(n,p);
        return p;
      }
      //Complex
      //Generate other Newicks and their coefficients
      std::vector<double> coefficients;
      std::vector<NewickType> newicks;
      {
        const double d = n.CalcDenominator(theta);
        #ifdef TRACE_NEWICK_CALCULATEPROBABILITY
        TRACE("Denominator for "
          + n.ToStr()
          + " = "
          + boost::lexical_cast<std::string>(d));
        #endif
        typedef std::pair<std::vector<int>,int> NewickFrequencyPair;
        const std::vector<NewickFrequencyPair> newick_freqs
          = Newick::GetSimplerNewicksFrequencyPairs(n.Peek());
        BOOST_FOREACH(const NewickFrequencyPair& p,newick_freqs)
        {
          const int frequency = p.second;
          assert(frequency > 0);
          if (frequency == 1)
          {
            newicks.push_back(p.first);
            coefficients.push_back(theta / d);
          }
          else
          {
            const double f_d = static_cast<double>(frequency);
            newicks.push_back(p.first);
            coefficients.push_back( (f_d*(f_d-1.0)) / d);
          }
          #ifdef TRACE_NEWICK_CALCULATEPROBABILITY
          TRACE(std::string("BinaryNewickVector ")
            + Newick::NewickToString(p.first)
            + " has coefficient "
            + boost::lexical_cast<std::string>(coefficients.back()))
          #endif
        }
      }
      //Ask help about these new Newicks
      {
        const int sz = newicks.size();
        assert(newicks.size() == coefficients.size() );
        double p = 0.0;
        for (int i=0; i!=sz; ++i)
        {
          //Recursive function call
          p+=(coefficients[i] * CalculateProbability(newicks[i],theta,storage));
        }
        storage.Store(n,p);
        return p;
      }
    }
    catch (std::bad_alloc& e)
    {
      storage.CleanUp();
      std::cerr << "std::bad_alloc\n";
    }
    catch (std::exception& e)
    {
      storage.CleanUp();
      std::cerr << "std::exception";
    }
    catch (...)
    {
      storage.CleanUp();
      std::cerr << "Unknown exception\n";
    }
  }
}
//---------------------------------------------------------------------------

  /*
  template <class NewickType>
  double CalculateProbabilityMessy(
    const NewickType& n,
    const double theta,
    NewickStorage<NewickType>& storage)
  {
  //#define TRACE_CALCULATEPROBABILITYMESSY
  while(1)
  {
    try
    {

      //Is n already known?
      {
        const double p = storage.Find(n);
        if (p!=0.0)
        {
          return p;
        }
      }

      //Check for simple phylogeny
      if (n.IsSimple())
      {
        const double p = n.CalcProbabilitySimpleNewick(theta);
        storage.Store(n,p);
        return p;
      }

      //Complex
      //Generate other Newicks and their coefficients
      std::vector<NewickType> newicks;
      std::vector<double> coefficients;
      {
        const double d = n.CalcDenominator(theta);
        #ifdef TRACE_CALCULATEPROBABILITYMESSY
        TRACE("Denominator for "
          + n.ToStr()
          + " = "
          + boost::lexical_cast<std::string>(d));
        #endif
        const std::vector<int>& v = n.Peek();
        const int sz = v.size();
        for (int i=0; i!=sz; ++i)
        {
          const int x = v[i];
          if (x < 0) //x is not a digit
          {
            continue;
          }
          if (x == 1)
          {
            NewickType next_newick = n.TermIsOne(i);
            if (!next_newick.Empty())
            {
              newicks.push_back(next_newick);
              coefficients.push_back(theta / d);
            }
          }
          else
          {
            NewickType next_newick = n.TermIsNotOne(i);
            if (!next_newick.Empty())
            {
              newicks.push_back(next_newick);
              const double x_d = static_cast<double>(x);
              coefficients.push_back( (x_d*(x_d-1.0)) / d);
            }
          }
        }
      }
      //Ask help about these new Newicks
      {
        const int sz = newicks.size();
        assert(newicks.size() == coefficients.size() );
        double p = 0.0;
        for (int i=0; i!=sz; ++i)
        {
          //Recursive function call
          p+=(coefficients[i] * CalculateProbabilityMessy(newicks[i],theta,storage));
        }
        storage.Store(n,p);
        return p;
      }
    }
    catch (std::bad_alloc& e)
    {
      storage.CleanUp();
      std::cerr << "std::bad_alloc\n";
    }
    catch (std::exception& e)
    {
      storage.CleanUp();
      std::cerr << "std::exception";
    }
    catch (...)
    {
      storage.CleanUp();
      std::cerr << "Unknown exception\n";
    }

  }
}
*/
//#endif
//---------------------------------------------------------------------------

};
//---------------------------------------------------------------------------
#endif // NEWICK_H


 

 

 

 

 

newick.cpp

 

//---------------------------------------------------------------------------
/*
Newick, Newick functions
Copyright (C) 2010-2011 Richel Bilderbeek

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program.If not, see <http://www.gnu.org/licenses/>.
*/
//---------------------------------------------------------------------------
//From http://www.richelbilderbeek.nl/CppNewick.htm
//---------------------------------------------------------------------------
#ifdef _WIN32
//See http://www.richelbilderbeek.nl/CppCompileErrorSwprintfHasNotBeenDeclared.htm
#undef __STRICT_ANSI__
#endif
//---------------------------------------------------------------------------
#include <algorithm>
#include <cassert>
#include <deque>
#include <iostream>
#include <functional>
#include <map>
#include <numeric>
#include <stdexcept>
#include <string>
#include <sstream>
#include <vector>
//---------------------------------------------------------------------------
#include <boost/foreach.hpp>
#include <boost/lexical_cast.hpp>
#include <boost/numeric/conversion/cast.hpp>
//---------------------------------------------------------------------------
#include "BigIntegerLibrary.hh"
//---------------------------------------------------------------------------
#include "newick.h"
#include "newickcpp98.h"
#include "trace.h"
//---------------------------------------------------------------------------
namespace Newick {
//---------------------------------------------------------------------------
//From http://www.richelbilderbeek.nl/CppAccumulate_if.htm
template
  <
  typename InputIterator,
  typename ElementType,
  typename Predicate
  >
const ElementType Accumulate_if(
  InputIterator first,
  const InputIterator last,
  ElementType init,
  const Predicate predicate)
{
  for (; first != last; ++first)
    if (predicate(*first)) init += *first;
  return init;
}
//---------------------------------------------------------------------------
///Copy_if was dropped from the standard library by accident.
///From http://www.richelbilderbeek.nl/CppCopy_if.htm
template<typename In, typename Out, typename Pred>
Out Copy_if(In first, In last, Out res, Pred Pr)
{
  while (first != last)
  {
    if (Pr(*first))
      *res++ = *first;
    ++first;
  }
  return res;
}
//---------------------------------------------------------------------------
//From http://www.richelbilderbeek.nl/CppFunctorIncrease.htm
struct Increase : public std::unary_function<void,int>
{
  explicit Increase(const int& initValue = 0) : mValue(initValue) {}
  void operator()(int& anything)
  {
    anything = mValue;
    ++mValue;
  }
  private:
  int mValue;
};
//---------------------------------------------------------------------------
// Free functions
//---------------------------------------------------------------------------
///AllAboutEqual tests if all values in a std::vector are about equal.
///From http://www.richelbilderbeek.nl/CppAllAboutEqual.htm
bool AllAboutEqual(
  const std::vector<double>& v,
  const double tolerance)
{
  assert(!v.empty());
  return std::count_if(
    v.begin(),
    v.end(),
    std::bind2nd(fuzzy_equal_to(tolerance),v[0]))
    == boost::numeric_cast<int>(v.size());

}
//---------------------------------------------------------------------------
///CalcComplexity calculates the complexity of a Newick.
///From http://www.richelbilderbeek.nl/CppCalcComplexity.htm
const BigInteger CalcComplexity(const std::vector<int>& v)
{
  if (v.empty()) return 0;
  //assert(IsNewick(v));
  BigInteger complexity = 1;
  int n_frequencies = 0;
  const int sz = v.size();
  for (int i=0; i!=sz; ++i)
  {
    const int x = v[i];
    if (x < 0) continue; //Ignore if x is not a number
    ++n_frequencies;
    complexity*=x;
  }
  complexity*=n_frequencies;
  return complexity;
}
//---------------------------------------------------------------------------
double CalcDenominator(const std::vector<int>& v,const double theta)
{
  int sum_above_zero = 0;
  int sum_above_one  = 0;
  BOOST_FOREACH(const int& i,v)
  {
    if (i > 0) sum_above_zero+= i;
    if (i > 1) sum_above_one += i;
  }
  const double d
    = boost::numeric_cast<double>(
      sum_above_zero * (sum_above_zero - 1))
    + (boost::numeric_cast<double>(sum_above_one)
       * theta);
  return d;
}
//---------------------------------------------------------------------------
///CalcNumOfCombinations returns the number of combinations a Newick can have.
///
///The number of possible combinations equals
///     !(v0 + v1 + v2 + etc)
/// N = -------------------------- / 2^number_of_symmetries
///     !v0 * !v1 * !v2 * etc
///
///      n
///   = --- / 2^number_of_symmetries
///      d
///
/// where v denotes an element in vector v in range [1,-> >
/// where v0 denotes the first element in vector v
/// and where !v0 denotes the factorial of v0
///     {factorial(!SUM(v)) product terms}
/// N = --------------------------------------------------
///     {product terms} + { number_symmetries times a '2'}
///
///     numerator_terms
/// N = --------------------------------------------------
///     denominator_terms with appended number_symmetries times a '2'
///
///From http://www.richelbilderbeek.nl/CppCalcNumOfCombinationsBinary.htm
const BigInteger CalcNumOfCombinationsBinary(const std::vector<int>& v)
{
  assert(Newick::IsNewick(v));

  //Get all positives
  std::vector<BigInteger> positives;
  Copy_if(
    v.begin(),v.end(),
    std::back_inserter(positives),
    std::bind2nd(std::greater<BigInteger>(),0));

  //Obtain numerator = (SUM(x))!
  const int sum_values = Accumulate_if(v.begin(),v.end(),0,std::bind2nd(std::greater<int>(),0));

  //std::clog << "sum_values:" << sum_values << '\n';
  BigInteger numerator = FactorialBigInt(sum_values);
  //std::clog << "Numerator:" << numerator << '\n';

  //Obtain factorialated positives
  BigInteger denominator = 1;
  BOOST_FOREACH(const int& i, v)
  {
    if (i<=0) continue;
    const BigInteger i_temp = FactorialBigInt(i);
    denominator*=i_temp;
    //std::clog << "Denominator:" << denominator << '\n';
  }

  //Obtain number_of_symmetries
  const BigInteger number_of_symmetries = CalcNumOfSymmetriesBinary(v);

  //Add number_of_symmetries times a 2 to denominator terms
  for(BigInteger i=0; i!=number_of_symmetries; ++i)
  {
    denominator*=2;
    //std::clog << "Denominator:" << denominator << '\n';
  }

  //Return the division
  numerator/=denominator;
  return numerator;
}
//---------------------------------------------------------------------------
///CalcNumOfSymmetries calculates the number of symmetries in a Newick.
///From http://www.richelbilderbeek.nl/CppCalcNumOfSymmetriesBinary.htm
const BigInteger CalcNumOfSymmetriesBinary(std::vector<int> v)
{
  assert(IsNewick(v));
  assert(IsBinaryNewick(v));
  if (v.size() == 3) return 0;
  if (v.size() == 4) return (v[1] > 0 && v[1]==v[2] ? 1 : 0);

  const int n_reserved
    = *std::max_element(v.begin(),v.end())
    + std::count_if(v.begin(), v.end(), std::bind2nd(std::greater<int>(),0));

  BigInteger n_symmetries = 0;
  int id = n_reserved + 1;

  std::map<std::pair<int,int>,int> ids;

  while (1)
  {
    //std::copy(v.begin(),v.end(),std::ostream_iterator<int>(std::clog," ")); std::clog << '\n';
    //Count number of symmetries
    assert(!v.empty());
    const std::size_t sz = v.size();
    assert(sz >= 2);
    const std::size_t j = sz - 1;
    for (std::size_t i = 0; i!=j; ++i)
    {
      if (v[i] > 0 && v[i]==v[i+1]) ++n_symmetries;
    }
    //Collect all leaves and store new leaves
    std::vector<std::pair<int,int> > leaves;
    for (std::size_t i = 0; i!=j; ++i)
    {
      if (v[i] > 0 && v[i+1] > 0)
      {
        //Keep pair sorted
        const std::pair<int,int> p
          = (v[i] <= v[i+1]
          ? std::make_pair(v[i+0],v[i+1])
          : std::make_pair(v[i+1],v[i+0]));
        //If this leaf is new, store it
        if (ids.find(p)==ids.end())
        {
          ids[p] = id;
          ++id;
        }
      }
    }
    //Generalize all leaves
    for (std::size_t i = 0; i < v.size()-1; ++i)
    {
      assert(v.size()>2);
      if (v[i] > 0 && v[i+1] > 0)
      {
        //Keep pair sorted
        const std::pair<int,int> p
          = (v[i] <= v[i+1]
          ? std::make_pair(v[i+0],v[i+1])
          : std::make_pair(v[i+1],v[i+0]));
        //If this leaf is new, store it
        assert(ids.find(p)!=ids.end() && "Leaf should have been stored already");
        assert(i > 0);
        std::vector<int> v_new;
        std::copy(v.begin(),v.begin() + i - 1,std::back_inserter(v_new));
        const int id = ids[p];
        v_new.push_back(id);
        std::copy(v.begin() + i + 3,v.end(),std::back_inserter(v_new));
        v = v_new;
        i = (i-1 > 0 ? i-1 : 0);
        //Check if there are more leaves to be generalized
        if (v.size()<=4)
        {
          //Check if the last (X,Y) is symmetrical...
          return n_symmetries + (v[1] > 0 && v[1]==v[2] ? 1 : 0);
          break;
        }
      }
    }
  }
}
//---------------------------------------------------------------------------
///CalcProbabilitySimpleNewick returns the probability of
///a Newick for a value of theta
///using the Ewens formula
double CalcProbabilitySimpleNewick(
  const std::vector<int>& v,
  const double theta)
{
  assert(IsNewick(v));
  assert(IsSimple(v));
  const int sz = v.size();

  int n=0;
  int k=0;

  double probability = 1.0;

  for (int i=0; i!=sz; ++i)
  {
    if (v[i]>0)
    {
      const int ni = v[i];
      ++k;
      ++n;
      for (int p=1; p!=ni; ++p, ++n)
      {
        probability *= (static_cast<double>(p)
          / ( static_cast<double>(n) + theta));
      }
      probability /= ( static_cast<double>(n) + theta);
    }
  }
  probability *= (static_cast<double>(n)+theta)
    * std::pow(theta,static_cast<double>(k-1));
  return probability;
}
//---------------------------------------------------------------------------
///CheckNewick checks if a std::string is a valid Newick.
///If this std::string is not a valid Newick,
///CheckNewick throws an exception with a detailed description
///From http://www.richelbilderbeek.nl/CppCheckNewick.htm
void CheckNewick(const std::string& s)
{
  #ifndef NDEBUG
  //std::clog
  //  << "Researching Newick string: '"
  //  << s << "'\n";
  #endif

  if (s.size()<3)
    throw std::invalid_argument(
      "The Newick std::string must have a size of at least three characters");
  if (s[0]!='(')
    throw std::invalid_argument(
      "The Newick std::string must start with an opening bracket ('(').");
  if (s[s.size()-1]!=')')
    throw std::invalid_argument(
      "The Newick std::string must end with a closing bracket (')').");
  if (std::count(s.begin(),s.end(),'(')
    !=std::count(s.begin(),s.end(),')'))
    throw std::invalid_argument(
       "The Newick std::string must have as much opening "
       "as closing brackets #1");
  if (s.find("(0")!=std::string::npos)
    throw std::invalid_argument(
      "A std::string Newick frequency cannot be or "
      "start with a zero (#1)");
  if (s.find(",0")!=std::string::npos)
    throw std::invalid_argument(
      "A std::string Newick frequency cannot be or "
      "start with a zero (#2)");
  if (s.find("()")!=std::string::npos)
    throw std::invalid_argument(
      "The Newick std::string cannot have "
      "a consecutive opening and closing bracket");
  if (s.find(",,")!=std::string::npos)
    throw std::invalid_argument(
      "A Newick std::string can have no consecutive comma's");
  if (s.find("(,")!=std::string::npos)
    throw std::invalid_argument(
      "A Newick std::string cannot have comma "
      "directly after an opening bracket");
  if (s.find(",)")!=std::string::npos)
    throw std::invalid_argument(
      "A Newick std::string cannot have comma "
      "directly before a closing bracket");

  std::string s_copy = s;
  while(s_copy.size()>2) //Find a leaf and cut it until the string is empty
  {
    #ifndef NDEBUG
    //std::clog
    //  << "Researching Newick string leaf: '"
    //  << s_copy
    //  << "'\n";
    #endif
    //Find a leaf
    //Find index i (starting opening bracket) and j (closing bracket)
    const std::size_t sz = s_copy.size();
    std::size_t i = 0;
    std::size_t j = 0;
    for (i=0 ; i!=sz; ++i) //Index of opening bracket
    {
      if (s_copy[i]!='(') continue;

      assert(s_copy[i]=='(');
      assert(i+1 < s_copy.size());
      for (j=i+1; j!=sz; ++j)
      {
        if (s_copy[j]=='(') { j = 0; break; }
        if (s_copy[j]!=')') continue;
        break;
      }

      if (j ==  0) continue; //j cannot be 0 after previous for loop
      if (j == sz)
        throw std::invalid_argument(
          "The Newick std::string must have as much opening as closing brackets #2");
      break;
    }
    if (s_copy[i]!='(')
      throw std::invalid_argument(
        "The Newick std::string must have as much opening as closing brackets #3");
    //Indices i and j found
    //Is range between i and j valid?
    if (s_copy[i]!='(')
      throw std::logic_error(
        "Bilderbikkel incorrectly assumes that s_copy[i]=='('");
    if (s_copy[j]!=')')
      throw std::logic_error(
        "Bilderbikkel incorrectly assumes that s_copy[j]==')'");
    //Check the range
    for (size_t k=i+1; k!=j; ++k)
    {
      if ( s_copy[k]!='0'
        && s_copy[k]!='1'
        && s_copy[k]!='2'
        && s_copy[k]!='3'
        && s_copy[k]!='4'
        && s_copy[k]!='5'
        && s_copy[k]!='6'
        && s_copy[k]!='7'
        && s_copy[k]!='8'
        && s_copy[k]!='9'
        && s_copy[k]!='0'
        && s_copy[k]!=',')
      {
        std::stringstream err_msg;
        err_msg << "Invalid non-number character in input: '" << s_copy[k] << "'";
        throw std::invalid_argument(err_msg.str().c_str());
      }
    }
    if (i > 0 && s_copy[i-1]=='(' &&
      j +1 < sz && s_copy[j + 1] == ')')
    {
      throw std::invalid_argument(
        "Newicks must not have the form ((X))");
    }
    //Check if there is a comma somewhere between brackets
    if (i > 0 //< (1) is valid, (1,(2)) not, ((1),2) not
      && std::find(
        s_copy.begin()+i,s_copy.begin()+j,',')
          == s_copy.begin()+j)
    {
      throw std::invalid_argument(
        "The Newick std::string cannot have the sequence "
        "of an opening bracket, a value and a closing bracket "
        "as a \'complex\' leaf");
    }
    //Range is assumed valid
    //Cut the leaf (turns '(1,2)' to (9) )
    assert(s_copy[i]=='(');
    assert(s_copy[j]==')');
    const std::string s_new_1 = s_copy.substr(0,i);
    const std::string s_new_2 = s_copy.substr(j+1,sz-j-1); //OK
    const std::string s_new =  s_new_1 + "9" + s_new_2;
    s_copy = s_new;
  }
}
//---------------------------------------------------------------------------
///CheckNewick checks if a std::vector<int> is a valid Newick.
///If this std::vector<int> is not a valid Newick,
///CheckNewick throws an exception with a detailed description
///From http://www.richelbilderbeek.nl/CppCheckNewick.htm
void CheckNewick(const std::vector<int>& v)
{
  #ifndef NDEBUG
  //std::clog << "Researching newick: '"
  //  << DumbNewickToString(v) << '\n';
  #endif

  if (v.size()<3)
    throw std::invalid_argument(
      "The Newick std::vector<int> must have a size of at least three characters");
  if (v[0]!=bracket_open)
    throw std::invalid_argument(
      "The Newick std::vector<int> must start with an opening bracket ('(').");
  if (v[v.size()-1]!=bracket_close)
    throw std::invalid_argument(
      "The Newick std::vector<int> must end with a closing bracket (')').");
  if (std::count(v.begin(),v.end(),static_cast<int>(bracket_open))
    !=std::count(v.begin(),v.end(),static_cast<int>(bracket_close)))
    throw std::invalid_argument(
       "The Newick std::string must have as much opening "
       "as closing brackets #1");
  if (std::count(v.begin(),v.end(),0))
    throw std::invalid_argument(
      "A std::vector<int> Newick frequency cannot be "
      "zero");


  std::vector<int> v_copy = v;
  while(v_copy.size()>2) //Find a leaf and cut it until the string is empty
  {
    #ifndef NDEBUG
    //std::clog << "Researching leaf: '";
    //std::copy(v_copy.begin(),v_copy.end(),std::ostream_iterator<int>(std::clog," "));
    //std::clog << "'\n";
    #endif
    //Find a leaf
    //Find index i (starting opening bracket) and j (closing bracket)
    const std::size_t sz = v_copy.size();
    std::size_t i = 0;
    std::size_t j = 0;
    for (i=0 ; i!=sz; ++i) //Index of opening bracket
    {
      if (v_copy[i]!=bracket_open) continue;

      assert(v_copy[i]==bracket_open);

      assert(i+1 < v_copy.size());

      if (v_copy[i+1]==bracket_close)
        throw std::invalid_argument(
          "The Newick std::vector<int> cannot have "
          "a consecutive opening and closing bracket");

      for (j=i+1; j!=sz; ++j)
      {
        if (v_copy[j]==bracket_open) { j = 0; break; }
        if (v_copy[j]!=bracket_close) continue;
        break;
      }
      if (i + 2 == j && j < sz - 1) //< (1) is valid, (1,(2)) not, ((1),2) not
      {
        throw std::invalid_argument(
          "The Newick std::vector<int> cannot have the sequence"
          "of an opening bracket, a value and a closing bracket"
          "as a \'complex\' leaf");
      }

      if (j ==  0) continue; //j cannot be 0 after previous for loop
      if (j == sz)
        throw std::invalid_argument(
          "The Newick std::vector<int> must have as much opening "
          "as closing brackets #2");
      break;
    }
    if (v_copy[i]!=bracket_open)
      throw std::invalid_argument(
        "The Newick std::vector<int> must have as much opening "
        "as closing brackets #3");
    //Indices i and j found
    //Is range between i and j valid?
    if (v_copy[i]!=bracket_open)
      throw std::logic_error(
        "Bilderbikkel incorrectly assumes that s_copy[i]=='('");
    if (v_copy[j]!=bracket_close)
      throw std::logic_error(
        "Bilderbikkel incorrectly assumes that s_copy[j]==')'");
    //Check the range
    for (size_t k=i+1; k!=j; ++k)
    {
      if (v_copy[k] < 0)
      {
        std::ostringstream err_msg;
        err_msg << "Invalid non-number in input: '" << v_copy[k] << "'";
        throw std::invalid_argument(err_msg.str().c_str());
      }
    }
    //Range is assumed valid
    //Cut the leaf
    //Changes '(1,2)' to '(999)'
    assert(v_copy[i]==bracket_open);
    assert(v_copy[j]==bracket_close);
    std::vector<int> v_new(v_copy.begin(),v_copy.begin() + i);
    v_new.push_back(999);
    std::copy(v_copy.begin() + j + 1, v_copy.end(),std::back_inserter(v_new));
    v_copy = v_new;
  }
}
//---------------------------------------------------------------------------
///CreateInvalidNewicks creates std::strings
///that cannot and must not be converted to a Newick
///From http://www.richelbilderbeek.nl/CppCreateInvalidNewicks.htm
const std::vector<std::string> CreateInvalidNewicks()
{
  std::vector<std::string> v;
  v.push_back("");
  v.push_back("(");
  v.push_back(")");
  v.push_back("1");
  v.push_back("1234");
  v.push_back(")1234(");
  v.push_back("()1234()");
  v.push_back("(1234,)");
  v.push_back("(,1234,)");
  v.push_back("()");
  v.push_back("(0)");
  v.push_back("(-)");
  v.push_back("(-1)");
  v.push_back("(0,0)");
  v.push_back("(1,0)");
  v.push_back("(0,1)");
  v.push_back("(0,(1,1))");
  v.push_back("(1,(0,1))");
  v.push_back("(1,(1,0))");
  v.push_back("((0,1),1)");
  v.push_back("((1,0),1)");
  v.push_back("((1,1),0)");
  v.push_back("((2))");
  v.push_back("(1,(2,3)");
  v.push_back("(1,(2))");
  v.push_back("(1,((3)))");
  v.push_back("(11,(22,33)");
  v.push_back("(22,33),33)");
  v.push_back("1,2");
  v.push_back("(1,1),");
  v.push_back("(2,2),");
  v.push_back("((2,2),2),");
  v.push_back(",(1,1)");
  v.push_back(",(2,2)");
  v.push_back(",((2,2),2)");
  v.push_back(",(1,1),");
  v.push_back(",(2,2),");
  v.push_back("(2,(1,1)),");
  v.push_back(",((2,2),2),");
  v.push_back("(1,2");
  v.push_back("(-1,2");
  v.push_back("(1,-2");
  v.push_back("(0,-2");
  v.push_back("(-0,2");
  v.push_back("(1,,2)");
  v.push_back("(1,2))");
  v.push_back("(1,(2),3)");
  v.push_back("((1,2),(3,4))()");
  return v;
}
//---------------------------------------------------------------------------
///CreateRandomNewick creates an unsorted Newick string,
///with n values, with each value e [0,max>.
///From http://www.richelbilderbeek.nl/CppCreateRandomNewick.htm
const std::string CreateRandomNewick(const int n,const int max)
{
  const std::vector<int> v = CreateRandomBinaryNewickVector(n,max);
  return NewickToString(v);
}
//---------------------------------------------------------------------------
///CreateRandomBinaryNewickVector creates an unsorted Newick
///std::vector<int>, with n values, with each value e [0,max>.
///From http://www.richelbilderbeek.nl/CppCreateRandomBinaryNewickVector.htm
const std::vector<int> CreateRandomBinaryNewickVector(const int n,const int max)
{
  assert(n>0);
  assert(max>1);

  std::vector<int> v;
  v.reserve(2 + (n*2));

  v.push_back(bracket_open);

  v.push_back(1 + (std::rand() % (max-1) ));
  if (n==1)
  {
    v.push_back(bracket_close);
    return v;
  }
  v.push_back(1 + (std::rand() % (max-1) ));

  v.push_back(bracket_close); //??? IntVector format has no trailing bracket

  if (n==2)
  {
    return v;
  }

  for (int i=2; i!=n; ++i)
  {
    if ((std::rand() >> 4) % 2)
    {
      //Append
      std::vector<int> new_v;
      new_v.reserve(2 + v.size());
      new_v.push_back(bracket_open);

      std::copy(v.begin(),v.end(),std::back_inserter(new_v));

      new_v.push_back(1 + (std::rand() % (max-1)));
      new_v.push_back(bracket_close);
      std::swap(v,new_v);
    }
    else
    {
      //Prepend
      std::vector<int> new_v;
      new_v.reserve(2 + v.size());
      new_v.push_back(bracket_open);
      new_v.push_back(1 + (std::rand() % (max-1)));

      std::copy(v.begin(),v.end(),std::back_inserter(new_v));

      new_v.push_back(bracket_close);
      std::swap(v,new_v);
    }
    assert(std::count(v.begin(),v.end(),static_cast<int>(bracket_open ))
        == std::count(v.begin(),v.end(),static_cast<int>(bracket_close)));
  }
  assert(IsNewick(v));
  return v;
}
//---------------------------------------------------------------------------
///CreateValidBinaryNewicks creates std::strings
///that can be converted to a BinaryNewickVector.
///From http://www.richelbilderbeek.nl/CppCreateValidBinaryNewicks.htm
const std::vector<std::string> CreateValidBinaryNewicks()
{
  std::vector<std::string> v;
  v.push_back("(1,2)");
  v.push_back("(11,22)");
  v.push_back("(1,(1,1))");
  v.push_back("(1,(1,2))");
  v.push_back("(1,(2,1))");
  v.push_back("(1,(2,2))");
  v.push_back("(1,(2,3))");
  v.push_back("(2,(1,1))");
  v.push_back("(2,(1,2))");
  v.push_back("(2,(2,1))");
  v.push_back("(2,(2,2))");
  v.push_back("(4,(2,3))");
  v.push_back("((2,3),4)");
  v.push_back("(2,((2,3),4))");
  v.push_back("(11,(22,33))");
  v.push_back("((22,33),33)");
  v.push_back("((1,2),(3,4))");
  v.push_back("(((1,2),(3,4)),5)");
  v.push_back("(1,((2,3),(4,5)))");
  v.push_back("((11,2),(3,44))");
  v.push_back("(((1,22),(33,4)),(55,6))");
  return v;
}
//---------------------------------------------------------------------------
///CreateValidTrinaryNewicks creates std::strings
///that can be converted to a TrinaryNewickVector.
///From http://www.richelbilderbeek.nl/CppCreateValidTinaryNewicks.htm
const std::vector<std::string> CreateValidTrinaryNewicks()
{
  //#ifdef __GXX_EXPERIMENTAL_CXX0X__
  ///\note
  ///The tests below must be put back in again once
  #ifdef DEBUG_TEMP_REMOVE_2738236826438
  return
    {
      "(1,1,1)",
      "(1,2,3)",
      "((1,1),1,1)",
      "(1,(1,1),1)",
      "(1,1,(1,1))",
      "(1,(2,3,4))",
      "(1,2,(3,4))",
      "(1,2,(3,4,5))",
      "((1,2,3),4,5)",
      "(11,22,33)",
      "(11,(22,33,44))",
      "(11,22,(33,44))",
      "(11,22,(33,44,55))",
      "((11,22,33),44,55)",
      "((1,2),(3,4),(5,6))",
      "((1,2,3),(4,5),(6,7))",
      "((1,2),(3,4,5),(6,7))",
      "((1,2),(3,4),(5,6,7))",
      "((1,2,3),(4,5),(6,7))",
      "((1,2),(3,4,5),(6,7))",
      "((1,2),(3,4),(5,6,7))",
      "((1,2,3),(4,5,6),(7,8))",
      "((1,2),(3,4,5),(6,7,8))",
      "((1,2,3),(4,5),(6,7,8))",
      "((1,2,3),(4,5,6),(7,8,9))",
      "((11,22,33),(44,55,66),(77,88,99))"
    };
  #else
    return NewickCpp98::CreateValidTrinaryNewicks();
  #endif
}
//---------------------------------------------------------------------------
///CreateValidNewicks creates std::strings
///that are valid newicks.
///From http://www.richelbilderbeek.nl/CppCreateValidNewicks.htm
const std::vector<std::string> CreateValidNewicks()
{
  std::vector<std::string> v;
  {
    std::vector<std::string> w = CreateValidUnaryNewicks();
    std::copy(w.begin(),w.end(),std::back_inserter(v));
  }
  {
    std::vector<std::string> w = CreateValidBinaryNewicks();
    std::copy(w.begin(),w.end(),std::back_inserter(v));
  }
  {
    std::vector<std::string> w = CreateValidTrinaryNewicks();
    std::copy(w.begin(),w.end(),std::back_inserter(v));
  }
  v.push_back("(1,2,3,4)");
  v.push_back("(1,2,3,4,5)");
  v.push_back("(1,2,3,4,5,6)");
  v.push_back("(1,2,3,4,5,6,7)");
  v.push_back("(1,2,3,4,5,6,7,8)");
  v.push_back("((1,2,3,4,5,6,7,8),(1,2,3,4,5,6,7,8),(1,2,3,4,5,6,7,8))");
  return v;
}
//---------------------------------------------------------------------------
///CreateValidBinaryNewicks creates std::strings
///that can be converted to a BinaryNewickVector.
///From http://www.richelbilderbeek.nl/CppCreateValidBinaryNewicks.htm
const std::vector<std::string> CreateValidUnaryNewicks()
{
  std::vector<std::string> v;
  v.push_back("(1)");
  v.push_back("(9)");
  v.push_back("(123)");
  return v;
}
//---------------------------------------------------------------------------
///DumbNewickToString converts a Newick std::vector<int> to a
///standard-format std::string without error checking.
///From http://www.richelbilderbeek.nl/CppDumbNewickToString.htm
const std::string DumbNewickToString(const std::vector<int>& v)
{
  std::string s;
  s.reserve(2 * v.size()); //Just a guess
  const int sz = v.size();
  for (int i=0; i!=sz; ++i)
  {
    const int x = v[i];
    if (x >= 0)
    {
      s+=boost::lexical_cast<std::string>(x);
      const int next = v[i+1];
      if (next > 0 || next == bracket_open)
      {
        s+=",";
      }
    }
    else if (x==bracket_open)
    {
      s+="(";
    }
    else if (x==bracket_close)
    {
      s+=")";
      //Final closing bracket?
      if (i+1==sz) break;
      const int next = v[i+1];
      if (next > 0 || next == bracket_open)
      {
        s+=",";
      }
    }
    else
    {
      s+="x"; //Unknown character
    }
  }
  return s;
}
//---------------------------------------------------------------------------
///Factorial calculates the factorial of all std::vector elements.
///From http://www.richelbilderbeek.nl/CppFactorial.htm
const std::vector<int> Factorial(const std::vector<int>& v_original)
{
  std::vector<int> v(v_original);
  std::transform(v.begin(),v.end(),v.begin(),std::ptr_fun<int,int>(Factorial));
  return v;
}
//---------------------------------------------------------------------------
///Factorial calculates the factorial of a value.
///From http://www.richelbilderbeek.nl/CppFactorial.htm
int Factorial(const int n)
{
  assert(n>=0);
  int result = 1;
  for (int i=1; i<=n; ++i)
  {
    result*=i;
  }
  return result;
}
//---------------------------------------------------------------------------
///FactorialBigInt returns the factorial of an integer
///as a BigInteger.
///From http://www.richelbilderbeek.nl/CppFactorialBigInt.htm
const BigInteger FactorialBigInt(const int n)
{
  assert(n>=0);
  BigInteger result = 1;
  for (int i=1; i<=n; ++i)
  {
    result*=i;
  }
  return result;
}
//---------------------------------------------------------------------------
int FindPosAfter(const std::vector<int>& v,const int x, const int index)
{
  const int sz = v.size();
  for (int i=index+1; i!=sz; ++i)
  {
    if (v[i]==x) return i;
  }
  return sz;
}
//---------------------------------------------------------------------------
int FindPosBefore(const std::vector<int>& v,const int x, const int index)
{
  for (int i=index-1; i!=-1; --i)
  {
    if (v[i]==x) return i;
  }
  return -1;
}
//---------------------------------------------------------------------------
///GetDepth returns the depth of each Newick element
///Example #1
///(1,2,3)
///01 1 10
///Example #2
///((1,2),(3,(4,5)))
///000 00 00 00 0000 <- depth layer 0 (comma's are skipped)
///.11 11 11 11 111. <- depth layer 1
///... .. .. .2 22.. <- depth layer 2
///011 11 11 22 2210 <- result of GetDepth
const std::vector<int> GetDepth(const std::vector<int>& n)
{
  assert(IsNewick(n));
  std::vector<int> v;
  int depth = -1;
  BOOST_FOREACH(const int& x,n)
  {
    if (x == Newick::bracket_open) ++depth;
    v.push_back(depth);
    if (x == Newick::bracket_close) --depth;
  }
  assert(n.size() == v.size());
  return v;
}
//---------------------------------------------------------------------------
///GetFactorialTerms returns all terms from a factorial.
///For example, 4! return {1,2,3,4}
///From http://www.richelbilderbeek.nl/CppGetFactorialTerms.htm
const std::vector<int> GetFactorialTerms(const int n)
{
  std::vector<int> v(n);
  std::for_each(v.begin(), v.end(),Increase(1));
  assert(std::count(v.begin(),v.end(),0)==0);
  return v;
}
//---------------------------------------------------------------------------
std::vector<boost::tuple<std::string,double,double> > GetKnownProbabilities()
{
  //#ifdef __GXX_EXPERIMENTAL_CXX0X__
  ///\note
  ///The tests below must be put back in again once
  #ifdef DEBUG_TEMP_REMOVE_2738236826438
  const std::vector<boost::tuple<std::string,double,double> > v
    = {
    //Sum equals 1
    { "(1)"  , 10.0, 1.0000000 },
    //Sum equals 2
    { "(2)"  , 10.0, 0.0909091 },
    { "(1,1)", 10.0, 0.9090909 },
    //Sum equals 3
    { "(3)"      , 10.0, 0.0151515 },
    { "(1,2)"    , 10.0, 0.0757576 },
    { "(2,1)"    , 10.0, 0.0757576 },
    { "(1,(1,1))", 10.0, 0.2525253 },
    { "((1,1),1)", 10.0, 0.2525253 },
    //Trinary
    { "(1,1,1)"  , 10.0, 0.7575758 },
    //Sum equals 4
    { "(4)"      , 10.0, 0.0034965 },
    { "(1,3)"    , 10.0, 0.0116550 },
    { "(2,2)"    , 10.0, 0.0058275 },
    { "(3,1)"    , 10.0, 0.0116550 },
    { "(1,(1,2))", 10.0, 0.0194250 },
    { "(1,(2,1))", 10.0, 0.0194250 },
    { "(2,(1,1))", 10.0, 0.0194250 },
    { "((1,2),1)", 10.0, 0.0194250 },
    { "((2,1),1)", 10.0, 0.0194250 },
    { "((1,1),2)", 10.0, 0.0194250 },
    //Trinary
    { "(1,1,2)"    , 10.0, 0.0582751 },
    { "(1,2,1)"    , 10.0, 0.0582751 },
    { "(2,1,1)"    , 10.0, 0.0582751 },
    { "(1,1,(1,1))", 10.0, 0.1295001 }, //(1)(confirmed)
    { "(1,(1,1),1)", 10.0, 0.1295001 },
    { "((1,1),1,1)", 10.0, 0.1295001 },
    { "(1,(1,1,1))", 10.0, 0.0971251 }, //(2)(confirmed)
    { "((1,1,1),1)", 10.0, 0.0971251 },
    //Quadrary
    { "(1,1,1,1)", 10.0, 0.5827505 },
    //Sum equals 5
    { "(1,4)"        , 10.0, 0.0024975 },
    { "(2,3)"        , 10.0, 0.0008325 },
    { "(3,2)"        , 10.0, 0.0008325 },
    { "(4,1)"        , 10.0, 0.0024975 },
    { "(1,(1,3))"    , 10.0, 0.0028305 },
    { "(1,(2,2))"    , 10.0, 0.0012950 },
    { "(1,(3,1))"    , 10.0, 0.0028305 },
    { "(2,(1,2))"    , 10.0, 0.0014338 },
    { "(2,(2,1))"    , 10.0, 0.0014338 },
    { "(3,(1,1))"    , 10.0, 0.0026640 },
    //Trinary
    { "(1,1,(1,2))"  , 10.0, 0.0092731 }, //(3)(confirmed)
    { "(1,1,(2,1))"  , 10.0, 0.0092731 },
    { "(1,1,(1,1,1))", 10.0, 0.0348263 }, //(4)(confirmed)
    { "(1,(1,1,1),1)", 10.0, 0.0348263 },
    { "((1,1,1),1,1)", 10.0, 0.0348263 },
    { "(2,(1,1,1))"  , 10.0, 0.0070069 }, //(5)(confirmed)
    { "((1,1,1),2)"  , 10.0, 0.0070069 },
    { "(1,1,1,(1,1))", 10.0, 0.0692918 }, //(6)(confirmed)
    { "(1,2,(1,1))"  , 10.0, 0.0092223 }, //(7)(confirmed)
    { "(2,1,(1,1))"  , 10.0, 0.0092223 },
    { "(1,(1,1),2)"  , 10.0, 0.0092223 },
    { "(2,(1,1),1)"  , 10.0, 0.0092223 },
    { "((1,1),1,2)"  , 10.0, 0.0092223 },
    { "((1,1),2,1)"  , 10.0, 0.0092223 },
    { "(1,(1,1,2))"  , 10.0, 0.0069190 }, //(9)(confirmed)
    { "(1,(1,2,1))"  , 10.0, 0.0069190 },
    { "(1,(2,1,1))"  , 10.0, 0.0069190 },
    { "((1,1,2),1)"  , 10.0, 0.0069190 },
    { "((1,2,1),1)"  , 10.0, 0.0069190 },
    { "((2,1,1),1)"  , 10.0, 0.0069190 },
    //Quadrary
    { "(1,(1,1,1,1))", 10.0, 0.0415140 }, //(8)(confirmed)
    //Pentary
    { "(1,1,1,1,1)"  , 10.0, 0.4162504 },
    //Sum equals 6
    { "(1,5)"        , 10.0, 0.0006660 },
    { "(2,4)"        , 10.0, 0.0001665 },
    { "(3,3)"        , 10.0, 0.0001110 },
    { "(1,(1,4))"    , 10.0, 0.0005804 },
    { "(1,(2,3))"    , 10.0, 0.0001679 },
    { "(1,(3,2))"    , 10.0, 0.0001679 },
    { "(1,(4,1))"    , 10.0, 0.0005804 },
    { "(2,(1,3))"    , 10.0, 0.0001991 },
    { "(2,(2,2))"    , 10.0, 0.0000925 },
    { "(2,(3,1))"    , 10.0, 0.0001991 },
    { "(3,(1,2))"    , 10.0, 0.0001880 },
    { "(3,(2,1))"    , 10.0, 0.0001880 },
    { "(4,(1,1))"    , 10.0, 0.0005043 },
    //Trinary
    { "(1,1,(1,3))"  , 10.0, 0.0012712 },
    { "(1,1,(2,2))"  , 10.0, 0.0005563 },
    { "(1,1,(3,1))"  , 10.0, 0.0012712 },
    { "(1,(1,3),1)"  , 10.0, 0.0012712 },
    { "(1,(2,2),1)"  , 10.0, 0.0005563 },
    { "(1,(3,1),1)"  , 10.0, 0.0012712 },
    { "((1,3),1,1)"  , 10.0, 0.0012712 },
    { "((2,2),1,1)"  , 10.0, 0.0005563 },
    { "((3,1),1,1)"  , 10.0, 0.0012712 },
    { "(1,2,(1,2))"  , 10.0, 0.0006346 },
    { "(2,1,(1,2))"  , 10.0, 0.0006346 },
    { "(1,2,(2,1))"  , 10.0, 0.0006346 },
    { "(2,1,(2,1))"  , 10.0, 0.0006346 },
    { "(1,(2,1),2)"  , 10.0, 0.0006346 },
    { "(1,(1,2),2)"  , 10.0, 0.0006346 },
    { "(2,(2,1),1)"  , 10.0, 0.0006346 },
    { "(2,(1,2),1)"  , 10.0, 0.0006346 },
    { "(1,3,(1,1))"  , 10.0, 0.0011913 },
    { "(1,(1,1),3)"  , 10.0, 0.0011913 },
    { "((1,1),1,3)"  , 10.0, 0.0011913 },
    { "(3,(1,1),1)"  , 10.0, 0.0011913 },
    { "((1,1),3,1)"  , 10.0, 0.0011913 },
    { "(1,1,(1,1,2))", 10.0, 0.0023165 },
    { "(1,1,(1,2,1))", 10.0, 0.0023165 },
    { "(1,1,(2,1,1))", 10.0, 0.0023165 },
    { "(1,(1,1,2),1)", 10.0, 0.0023165 },
    { "(1,(1,2,1),1)", 10.0, 0.0023165 },
    { "(1,(2,1,1),1)", 10.0, 0.0023165 },
    { "((1,1,2),1,1)", 10.0, 0.0023165 },
    { "((1,2,1),1,1)", 10.0, 0.0023165 },
    { "((2,1,1),1,1)", 10.0, 0.0023165 },
    { "(1,2,(1,1,1))", 10.0, 0.0023323 },
    { "(2,1,(1,1,1))", 10.0, 0.0023323 },
    { "(1,(1,1,1),2)", 10.0, 0.0023323 },
    { "(2,(1,1,1),1)", 10.0, 0.0023323 },
    { "((1,1,1),1,2)", 10.0, 0.0023323 },
    { "((1,1,1),2,1)", 10.0, 0.0023323 },
    //Quadrary
    { "(1,(1,1,1,2))"  , 10.0, 0.0027574 },
    { "(1,(1,1,2,1))"  , 10.0, 0.0027574 },
    { "(1,(1,2,1,1))"  , 10.0, 0.0027574 },
    { "(1,(2,1,1,1))"  , 10.0, 0.0027574 },
    { "((1,1,1,2),1)"  , 10.0, 0.0027574 },
    { "((1,1,2,1),1)"  , 10.0, 0.0027574 },
    { "((1,2,1,1),1)"  , 10.0, 0.0027574 },
    { "((2,1,1,1),1)"  , 10.0, 0.0027574 },
    { "(2,(1,1,1,1))"  , 10.0, 0.0028154 },
    { "((1,1,1,1),2)"  , 10.0, 0.0028154 },
    //Pentary
    { "(1,(1,1,1,1,1))", 10.0, 0.0183824 }, //(7)
    { "((1,1,1,1,1),1)", 10.0, 0.0183824 },
    //Hexary
    { "(1,1,1,1,1,1)"  , 10.0, 0.2775003 }
  };
  return v;
  #else
    return NewickCpp98::GetKnownProbabilities();
  #endif
}
//---------------------------------------------------------------------------
int GetLeafMaxArity(const std::vector<int>& n)
{
  assert(IsNewick(n));
  const int size = boost::numeric_cast<int>(n.size());
  if (IsSimple(n)) return size - 2;

  int max = 0;
  for (int from = 0; from!=size; ++from)
  {
    if (n[from] != Newick::bracket_open) continue;
    for (int to = from+1; to!=size; ++to)
    {
      if (n[to] == Newick::bracket_open) break;
      if (n[to]  > 0) continue;
      if (n[to] == Newick::bracket_close)
      {
        assert(from < to);
        max = to - from - 1;
        break;
      }
    }
  }
  return max;
}
//---------------------------------------------------------------------------
///GetRootBranches obtains the root branches from a non-unary Newick.
///Examples:
///(1,2)               -> { 1     , 2             }
///(1,2,3)             -> { 1     , 2     , 3     }
///((1,1),(2,2),(3,3)) -> { (1,1) , (2,2) , (3,3) }
///From http://www.richelbilderbeek.nl/CppGetRootBranchesBinary.htm
const std::vector<std::vector<int> >
  GetRootBranches(const std::vector<int>& n)
{
  //#ifdef __GXX_EXPERIMENTAL_CXX0X__
  ///\note
  ///The tests below must be put back in again once
  #ifdef DEBUG_TEMP_REMOVE_2738236826438

  //#define DEBUG_GETROOTBRANCHES
  #ifdef  DEBUG_GETROOTBRANCHES
  TRACE_FUNC();
  TRACE(Newick::NewickToString(n));
  #endif

  assert(IsNewick(n));
  assert(!IsUnaryNewick(n));

  const int size = boost::numeric_cast<int>(n.size());
  std::vector<std::vector<int> > v;

  if (IsSimple(n))
  {
    for (int i=1; i!=size-1; ++i) //Skip brackets
    {
      v.push_back(
        CreateVector(
          static_cast<int>(bracket_open),
          n[i],
          static_cast<int>(bracket_close)));
    }
    assert(IsNewick(v.back()));
    assert(v.size() > 1);
    return v;
  }
  //Complex newick
  assert(!IsSimple(n));
  const std::vector<int> depth = GetDepth(n);

  assert(depth.size() == n.size());
  //Search for open and closing brackets in depth 1
  for (int i=0; i!=size; ++i)
  {
    if (depth[i] == 0 && n[i] > 0)
    {
      //C++0x initialization list
      v.push_back(
        {
          Newick::bracket_open,
          n[i],
          Newick::bracket_close
        }
      );
      assert(Newick::IsNewick(v.back()));
      continue;
    }
    if (depth[i] != 1 || n[i]!=Newick::bracket_open) continue;
    for (int j=i+1; j!=size; ++j)
    {
      if (depth[j] != 1 || n[j]!=Newick::bracket_close) continue;
      std::vector<int> w;
      w.push_back(Newick::bracket_open);
      std::copy(n.begin() + i + 1,n.begin() + j,std::back_inserter(w));
      w.push_back(Newick::bracket_close);
      assert(IsNewick(w));
      v.push_back(w);
      //Set from index i after current end
      i = j;
      break;
    }
  }
  assert(v.size() > 1);
  return v;
  #else
  return NewickCpp98::GetRootBranches(n);
  #endif
}
//---------------------------------------------------------------------------
///GetRootBranchesBinary obtains the two root branches from a binary Newick.
///Examples:
///(1,2)                 -> { 1             , 2     }
///(1,(2,3))             -> { 1             , (2,3) }
///((1,2),(3,4))         -> { (1,2)         , (3,4) }
///(((1,2),(3,4)),(5,6)) -> { ((1,2),(3,4)) , (5,6) }
///From http://www.richelbilderbeek.nl/CppGetRootBranchesBinary.htm
const std::pair<std::vector<int>,std::vector<int> >
  GetRootBranchesBinary(const std::vector<int>& n)
{
  assert(IsNewick(n));
  assert(IsBinaryNewick(n) && "Only binary Newicks can have two roots");

  assert(n[0] == bracket_open);
  assert(n[n.size()-1] == bracket_close);
  const int sz = boost::numeric_cast<int>(n.size());
  //Return the answer directly is Newick consists
  //out of two values only
  if (sz==4)
  {
    assert(n[1] > 0);
    assert(n[2] > 0);
    return std::make_pair(
      CreateVector(
        static_cast<int>(bracket_open),
        n[1],
        static_cast<int>(bracket_close)),
      CreateVector(
        static_cast<int>(bracket_open),
        n[2],
        static_cast<int>(bracket_close)));
  }
  typedef std::vector<int>::const_iterator Iterator;
  const Iterator j = n.end() - 1;
  for (Iterator i = n.begin() + 2; i!=j; ++i)
  {
    std::vector<int> lhs(n.begin() + 1,i);
    std::vector<int> rhs(i,n.end() - 1);
    if ( lhs.front() != Newick::bracket_open
      || lhs.back()  != Newick::bracket_close)
    {
      lhs = Surround(lhs);
    }
    if ( rhs.front() != Newick::bracket_open
      || rhs.back()  != Newick::bracket_close)
    {
      rhs = Surround(rhs);
    }
    //std::clog
    //  << NewickToString(lhs)
    //  << " and "
    //  << NewickToString(rhs)
    //  << '\n';
    if (IsNewick(lhs) && IsNewick(rhs))
    {
      return std::make_pair(lhs,rhs);
    }
  }
  assert(!"Should not get here");
  throw std::logic_error("Should not get here");
}
//---------------------------------------------------------------------------
///GetRootBranches obtains the two root branches from a binary Newick.
///Examples:
///(1,2)                 -> { 1             , 2     }
///(1,(2,3))             -> { 1             , (2,3) }
///((1,2),(3,4))         -> { (1,2)         , (3,4) }
///(((1,2),(3,4)),(5,6)) -> { ((1,2),(3,4)) , (5,6) }
///From http://www.richelbilderbeek.nl/CppGetRootBranches.htm
const std::pair<std::vector<int>,std::vector<int> >
  GetRootBranchesOld(const std::vector<int>& n)
{
  assert(IsNewick(n));
  assert(IsBinaryNewick(n));
  assert(n.size() > 3
    && "A Newick must have at least two values");
  assert(n[0] == bracket_open);
  assert(n[n.size()-1] == bracket_close);
  const int sz = boost::numeric_cast<int>(n.size());
  //Return the answer directly is Newick consists
  //out of two values only
  if (sz==4)
  {
    assert(n[1] > 0);
    assert(n[2] > 0);
    return std::make_pair(
      CreateVector(
        static_cast<int>(bracket_open),
        n[1],
        static_cast<int>(bracket_close)),
      CreateVector(
        static_cast<int>(bracket_open),
        n[2],
        static_cast<int>(bracket_close)));
  }
  //Check the kind of Newick
  //1) (x(yz))          return {    x       ,(yz)}
  //2) ((xy)z)          return { (xy)       , (z)}
  //3) ((ab)(cd))       return { (ab)       ,(cd)}
  //4) (((ab)(cd))(ef)) return { (((ab)(cd)),(ef))
  if (n[1]>0)
  {
    //1) (x(yz))    return {    x,(yz) }
    return std::make_pair(
      CreateVector(
        static_cast<int>(bracket_open)
        ,n[1],
        static_cast<int>(bracket_close)),
      std::vector<int>(n.begin()+2,n.end()-1));
  }
  else if (n[n.size()-2]>0)
  {
    //2) ((xy)z)    return { (xy), (z) }
    return std::make_pair(
      std::vector<int>(n.begin()+1,n.end()-2),
      CreateVector(
        static_cast<int>(bracket_open),
        n[sz-2],
        static_cast<int>(bracket_close)));
  }
  else
  {
    //3) ((ab)(cd)) return { (ab),(cd) }
    //Find ')('
    int i = 0;
    for (; i!=sz-1; ++i)
    {
      if ( n[i  ]==bracket_close
        && n[i+1]==bracket_open)
      break;
    }
    assert(i!=sz-1 && "Sequence )( must be found");
    return std::make_pair(
      std::vector<int>(n.begin()+1,n.begin()+i+1),
      std::vector<int>(n.begin()+i+1,n.end()-1));
  }
}
//---------------------------------------------------------------------------
///GetSimplerBinaryNewicks creates simpler, derived Newicks from a binary Newick.
///From http://www.richelbilderbeek.nl/CppGetSimplerBinaryNewicks.htm
const std::vector<std::vector<int> >
  GetSimplerBinaryNewicks(const std::vector<int>& n)
{
  assert(IsNewick(n));
  assert(IsUnaryNewick(n) || IsBinaryNewick(n));

  std::vector<std::vector<int> > newicks;

  //If newick is simple (by counting the number of opening brackets),
  //there are no simpler Newicks
  if (std::count( n.begin(),n.end(),
    static_cast<int>(bracket_open))==1)
  {
    //Simple newicks do not need to be simplified
    assert(n.size()==3 || n.size() == 4);
    assert(n[0]==bracket_open);
    assert(n[n.size()-1]==bracket_close);
    if (n.size() == 3)
    {
      if (n[1]>1)
      {
        std::vector<int> next(n);
        --next[1];
        assert(IsNewick(next));
        newicks.push_back(next);
      }
      return newicks;
    }
    assert(n.size()==4);
    if (n[1] == 1)
    {
      std::vector<int> next
        = CreateVector(
            static_cast<int>(bracket_open),
            n[2]+1,
            static_cast<int>(bracket_close));
      assert(IsNewick(next));
      newicks.push_back(next);
    }
    else
    {
      std::vector<int> next(n);
      --next[1];
      assert(IsNewick(next));
      newicks.push_back(next);
    }
    if (n[2] == 1)
    {
      std::vector<int> next
        = CreateVector(
            static_cast<int>(bracket_open),
            n[1]+1,
            static_cast<int>(bracket_close));
      assert(IsNewick(next));
      newicks.push_back(next);
    }
    else
    {
      std::vector<int> next(n);
      --next[2];
      assert(IsNewick(next));
      newicks.push_back(next);
    }
    return newicks;
  }

  //newick is complex
  //Generate other Newicks and their coefficients
  const int sz = n.size();
  for (int i=0; i!=sz; ++i)
  {
    const int x = n[i];
    if (x < 0) //x is not a digit
    {
      continue;
    }
    if (x == 1)
    {
      //If x is not next to a digit, there is no simplification
      if (n[i-1]<0 && n[i+1]<1) continue;
      //If next to the x in a digit, merge these and remove their brackets
      //Is the 1 left of a value?
      if (n[i-1]==bracket_open)
      {
        assert(n[i+1] > 0);
        const int new_value = n[i+1] + 1;
        std::vector<int> next(n.begin(),n.begin() + i - 1);
        next.push_back(new_value);
        assert(n[i+2] < 0);
        std::copy(n.begin() + i + 3,n.end(),std::back_inserter(next));
        #ifndef NDEBUG
        if (!IsNewick(next))
        {
          std::cerr << DumbNewickToString(n)  << '\n';
          std::cerr << DumbNewickToString(next) << '\n';
          InspectInvalidNewick(std::cerr,next);
        }
        #endif
        assert(IsNewick(next));
        newicks.push_back(next);
      }
      else
      {
        //Is the 1 to the right of a value?
        assert(n[i-1] > 0); //< The other value
        const int new_value = n[i-1] + 1;
        std::vector<int> next(n.begin(),n.begin()+i-2);
        next.push_back(new_value);
        assert(n[i+1] < 0);
        std::copy(n.begin() + i + 2,n.end(),std::back_inserter(next));
        assert(IsNewick(next));
        newicks.push_back(next);
      }
    }
    else
    {
      std::vector<int> next = n;
      --next[i];
      assert(IsNewick(next));
      newicks.push_back(next);
    }
  }
  return newicks;
}
//---------------------------------------------------------------------------
///GetSimplerNewicks creates simpler, derived Newicks from a Newick.
///From http://www.richelbilderbeek.nl/CppGetSimplerNewicks.htm
const std::vector<std::vector<int> >
  GetSimplerNewicks(const std::vector<int>& n)
{
  //#define DEBUG_GETSIMPLERNEWICKS
  assert(IsNewick(n));

  #ifdef DEBUG_GETSIMPLERNEWICKS
  TRACE_FUNC();
  #endif


  const std::vector<int> depths = GetDepth(n);

  std::vector<std::vector<int> > newicks;

  const int size = boost::numeric_cast<int>(n.size());
  for (int i = 0; i!=size; ++i)
  {
    assert(i >= 0);
    assert(i < size);
    if (n[i] < 1) continue;
    assert(n[i] > 0);
    if (n[i] > 1)
    {
      std::vector<int> new_newick(n);
      --new_newick[i];
      #ifdef DEBUG_GETSIMPLERNEWICKS
      {
        const std::string stored = Newick::NewickToString(new_newick);
        TRACE(stored);
      }
      #endif
      newicks.push_back(new_newick);
      continue;
    }
    assert(n[i] == 1); //Most difficult...
    const int depth = depths[i];
    //j must first decrement, later increment with the same code
    int j_end  = -1;
    int j_step = -1;
    for (int j=i-1; ; j+=j_step)
    {
      //j must first decrement, later increment with the same code
      if (j == j_end
        //|| depths[j] < depth
        || (depths[j] == depth && n[j] < 0))
      {
        if (j_step == -1)
        {
          j = i + 1;
          j_end = size;
          j_step = 1;
        }
        else
        {
          break;
        }
      }
      assert(i!=j);
      assert(j >= 0);
      assert(j < size);
      //Only take frequencies of the same depth into account
      if (n[j] < 1 || depths[j] != depth) continue;
      std::vector<int> new_newick_with_zero(n);
      --new_newick_with_zero[i];
      assert(new_newick_with_zero[i] == 0);
      ++new_newick_with_zero[j];
      //Remove brackets after possibly lonely value
      #ifdef DEBUG_GETSIMPLERNEWICKS
      {
        const std::string newick_str_with_zeroes = Newick::DumbNewickToString(new_newick_with_zero);
        TRACE(newick_str_with_zeroes);
        const std::string dist_i_j = boost::lexical_cast<std::string>(std::abs(i - j));
        TRACE(dist_i_j)
      }
      #endif
      //If there is only one or two values between
      //the brackets, and one of these values was a
      //1 becoming added to the other, nullify the
      //1 and both brackets:
      //'((1,1),2)' -> '(00102)' -> '(1,2)'
      if (std::abs(i - j) == 1)
        //|| (std::abs(i - j) == 2 && n[i] == 1))
      {
        const int index_bracket_open  = std::min(i,j) - 1;
        const int index_bracket_close = std::max(i,j) + 1;
        if ( new_newick_with_zero[index_bracket_open]  == Newick::bracket_open
          && new_newick_with_zero[index_bracket_close] == Newick::bracket_close)
        {
          new_newick_with_zero[index_bracket_open]  = 0;
          new_newick_with_zero[index_bracket_close] = 0;
        }
      }
      #ifdef DEBUG_GETSIMPLERNEWICKS
      {
        const std::string newick_str_with_more_zeroes = Newick::DumbNewickToString(new_newick_with_zero);
        TRACE(newick_str_with_more_zeroes);
      }
      #endif
      //Remove decremented i and possibly nullified brackets
      std::vector<int> new_newick;
      std::remove_copy(
        new_newick_with_zero.begin(),
        new_newick_with_zero.end(),
        std::back_inserter(new_newick),
        0);
      //Add brackets if these are removed
      if (new_newick.front() != Newick::bracket_open
        || new_newick.back() != Newick::bracket_close)
      {
        new_newick = Surround(new_newick);
      }
      #ifdef DEBUG_GETSIMPLERNEWICKS
      {
        const std::string newick_str_done = Newick::DumbNewickToString(new_newick);
        TRACE(newick_str_done);
      }
      #endif
      #define DEBUG_2436964926435498753298216832187
      #ifdef  DEBUG_2436964926435498753298216832187
      if (!IsNewick(new_newick))
      {
        TRACE(Newick::DumbNewickToString(new_newick));
      }
      #endif

      assert(IsNewick(new_newick));
      newicks.push_back(new_newick);
      continue;
    }
  }

  /*
  const int size = boost::numeric_cast<int>(n.size());
  for (int bracket_open = 0; bracket_open!=size; ++bracket_open)
  {
    if (n[bracket_open] != Newick::bracket_open) continue;
    for (int bracket_close = bracket_open + 1; bracket_close!=size; ++bracket_close)
    {
      if (n[bracket_close] == Newick::bracket_open) break;
      if (n[bracket_close]  > 0) continue;
      assert(bracket_open < bracket_close);
      assert(bracket_open + 1 != bracket_close);
      assert(n[bracket_open] == Newick::bracket_open);
      assert(n[bracket_open +1] > 0);
      assert(n[bracket_close-1] > 0);
      assert(n[bracket_close] == Newick::bracket_close);
      #ifdef DEBUG_GETSIMPLERNEWICKS
      std::clog
        << "Bracket open: "
        << bracket_open
        << '\n'
        << "Bracket close: "
        << bracket_close
        << '\n';
      #endif
      for (int i=bracket_open+1; i!=bracket_close; ++i)
      {
        assert(n[i] > 0);
        if (n[i] > 1)
        {
          std::vector<int> new_newick(n);
          --new_newick[i];
          #ifdef DEBUG_GETSIMPLERNEWICKS
          std::clog << "Store: " << Newick::NewickToString(new_newick) << '\n';
          #endif
          newicks.push_back(new_newick);
          continue;
        }
        assert(n[i] == 1); //Most difficult...
        for (int j=bracket_open+1; j!=bracket_close; ++j)
        {
          if (i==j) continue;
          if (n[j] < 1) continue;
          //Decrement index i to zero
          //Increment index j
          std::vector<int> new_newick_with_zero(n);
          --new_newick_with_zero[i];
          assert(new_newick_with_zero[i] == 0);
          ++new_newick_with_zero[j];
          //Remove brackets after possibly lonely value
          #ifdef DEBUG_GETSIMPLERNEWICKS
          std::clog << "1: "
            << Newick::DumbNewickToString(new_newick_with_zero)
            << " -> ["
            << i
            << "]="
            << new_newick_with_zero[i]
            << " - ["
            << j
            << "]="
            << new_newick_with_zero[j]
            << " = "
            << std::abs(i-j)
            << '\n';
          #endif
          if (std::abs(i - j) == 1)
          {
            const int index_bracket_open  = std::min(i,j) - 1;
            const int index_bracket_close = std::max(i,j) + 1;
            #ifdef DEBUG_GETSIMPLERNEWICKS
            std::clog
              << "["
              << index_bracket_open
              << "]-["
              << index_bracket_close
              << "]\n";
            #endif

            if ( new_newick_with_zero[index_bracket_open]  == Newick::bracket_open
              && new_newick_with_zero[index_bracket_close] == Newick::bracket_close)
            {
              new_newick_with_zero[index_bracket_open]  = 0;
              new_newick_with_zero[index_bracket_close] = 0;
              #ifdef DEBUG_GETSIMPLERNEWICKS
              std::clog << "2.5: " << Newick::DumbNewickToString(new_newick_with_zero) << '\n';
              #endif
            }
          }
          #ifdef DEBUG_GETSIMPLERNEWICKS
          std::clog << "2: " << Newick::DumbNewickToString(new_newick_with_zero) << '\n';
          #endif
          //Remove decremented i and possibly nullified brackets
          std::vector<int> new_newick;
          std::remove_copy(
            new_newick_with_zero.begin(),
            new_newick_with_zero.end(),
            std::back_inserter(new_newick),
            0);
          //Add brackets if these are removed
          if (new_newick.front() != Newick::bracket_open
            || new_newick.back() != Newick::bracket_close)
          {
            new_newick = Surround(new_newick);
          }
          assert(IsNewick(new_newick));
          #ifdef DEBUG_GETSIMPLERNEWICKS
          std::clog << "Store: " << Newick::NewickToString(new_newick) << '\n';
          #endif
          newicks.push_back(new_newick);
          continue;
        }
      }
      //If a closing bracket is found,
      //stop looking for the following closing bracket
      break;
    }
  }
  */
  return newicks;
}
//---------------------------------------------------------------------------
///GetSimplerNewicksFrequencyPairs creates simpler, derived Newicks from a Newick.
///Its simpler Newicks are identical to those created by GetSimplerNewicks.
///From http://www.richelbilderbeek.nl/CppGetSimplerNewicksFrequencyPairs.htm
const std::vector<std::pair<std::vector<int>,int> >
  GetSimplerNewicksFrequencyPairs(const std::vector<int>& n)
{
  //#ifdef __GXX_EXPERIMENTAL_CXX0X__
  #ifdef DEBUG_TEMP_REMOVE_2738236826438

  //#define DEBUG_GETSIMPLERNEWICKSFREQUENCYPAIRS
  #ifdef DEBUG_GETSIMPLERNEWICKSFREQUENCYPAIRS
  TRACE_FUNC();
  #endif

  assert(IsNewick(n));

  std::vector<std::pair<std::vector<int>,int> > newicks;
  const std::vector<int> depths = GetDepth(n);


  const int size = boost::numeric_cast<int>(n.size());
  for (int i = 0; i!=size; ++i)
  {
    assert(i >= 0);
    assert(i < size);
    if (n[i] < 1) continue;
    assert(n[i] > 0);
    if (n[i] > 1)
    {
      std::vector<int> new_newick(n);
      --new_newick[i];
      #ifdef DEBUG_GETSIMPLERNEWICKS
      {
        const std::string stored = Newick::NewickToString(new_newick);
        TRACE(stored);
      }
      #endif
      newicks.push_back( { new_newick,n[i] } );
      continue;
    }
    assert(n[i] == 1); //Most difficult...
    const int depth = depths[i];
    //j must first decrement, later increment with the same code
    int j_end  = -1;
    int j_step = -1;
    for (int j=i-1; ; j+=j_step)
    {
      //j must first decrement, later increment with the same code
      if (j == j_end
        //|| depths[j] < depth
        || (depths[j] == depth && n[j] < 0))
      {
        if (j_step == -1)
        {
          j = i + 1;
          j_end = size;
          j_step = 1;
        }
        else
        {
          break;
        }
      }
      assert(i!=j);
      assert(j >= 0);
      assert(j < size);
      //Only take frequencies of the same depth into account
      if (n[j] < 1 || depths[j] != depth) continue;
      std::vector<int> new_newick_with_zero(n);
      --new_newick_with_zero[i];
      assert(new_newick_with_zero[i] == 0);
      ++new_newick_with_zero[j];
      //Remove brackets after possibly lonely value
      #ifdef DEBUG_GETSIMPLERNEWICKS
      {
        const std::string newick_str_with_zeroes = Newick::DumbNewickToString(new_newick_with_zero);
        TRACE(newick_str_with_zeroes);
        const std::string dist_i_j = boost::lexical_cast<std::string>(std::abs(i - j));
        TRACE(dist_i_j)
      }
      #endif
      //If there is only one or two values between
      //the brackets, and one of these values was a
      //1 becoming added to the other, nullify the
      //1 and both brackets:
      //'((1,1),2)' -> '(00102)' -> '(1,2)'
      if (std::abs(i - j) == 1)
        //|| (std::abs(i - j) == 2 && n[i] == 1))
      {
        const int index_bracket_open  = std::min(i,j) - 1;
        const int index_bracket_close = std::max(i,j) + 1;
        if ( new_newick_with_zero[index_bracket_open]  == Newick::bracket_open
          && new_newick_with_zero[index_bracket_close] == Newick::bracket_close)
        {
          new_newick_with_zero[index_bracket_open]  = 0;
          new_newick_with_zero[index_bracket_close] = 0;
        }
      }
      #ifdef DEBUG_GETSIMPLERNEWICKS
      {
        const std::string newick_str_with_more_zeroes = Newick::DumbNewickToString(new_newick_with_zero);
        TRACE(newick_str_with_more_zeroes);
      }
      #endif
      //Remove decremented i and possibly nullified brackets
      std::vector<int> new_newick;
      std::remove_copy(
        new_newick_with_zero.begin(),
        new_newick_with_zero.end(),
        std::back_inserter(new_newick),
        0);
      //Add brackets if these are removed
      if (new_newick.front() != Newick::bracket_open
        || new_newick.back() != Newick::bracket_close)
      {
        new_newick = Surround(new_newick);
      }
      #ifdef DEBUG_GETSIMPLERNEWICKS
      {
        const std::string newick_str_done = Newick::DumbNewickToString(new_newick);
        TRACE(newick_str_done);
      }
      #endif
      #define DEBUG_2436964926435498753298216832187
      #ifdef  DEBUG_2436964926435498753298216832187
      if (!IsNewick(new_newick))
      {
        TRACE(Newick::DumbNewickToString(new_newick));
      }
      #endif

      assert(IsNewick(new_newick));
      newicks.push_back( { new_newick, 1 } );
      continue;
    }
  }

  return newicks;

  /*
  const int size = boost::numeric_cast<int>(n.size());
  for (int i = 0; i!=size; ++i)
  {
    if (n[i] < 1) continue;
    assert(n[i] > 0);
    if (n[i] > 1)
    {
      std::vector<int> new_newick(n);
      --new_newick[i];
      #ifdef DEBUG_GETSIMPLERNEWICKSFREQUENCYPAIRS
      std::clog << "Store: " << Newick::NewickToString(new_newick) << '\n';
      #endif
      newicks.push_back(std::make_pair(new_newick,n[i]));
      continue;
    }
    assert(n[i] == 1); //Most difficult...
    const int depth = depths[i];
    const int j_start = FindPosBefore(n,Newick::bracket_open,i);
    const int j_end   = FindPosAfter( n,Newick::bracket_close,i);
    assert(j_start >= 0);
    assert(j_end >= 0);
    assert(j_start <= boost::numeric_cast<int>(n.size()));
    assert(j_end <= boost::numeric_cast<int>(n.size()));
    for (int j=j_start; j!=j_end; ++j)
    {
      if (i==j) continue;
      if (n[j] < 1) continue;
      if (depths[j] != depth) continue;
      //Decrement index i to zero
      //Increment index j
      std::vector<int> new_newick_with_zero(n);
      --new_newick_with_zero[i];
      assert(new_newick_with_zero[i] == 0);
      ++new_newick_with_zero[j];
      //Remove brackets after possibly lonely value
      #ifdef DEBUG_GETSIMPLERNEWICKSFREQUENCYPAIRS
      std::clog << "1: "
        << Newick::DumbNewickToString(new_newick_with_zero)
        << " -> ["
        << i
        << "]="
        << new_newick_with_zero[i]
        << " - ["
        << j
        << "]="
        << new_newick_with_zero[j]
        << " = "
        << std::abs(i-j)
        << '\n';
      #endif
      if (std::abs(i - j) == 1)
      {
        const int index_bracket_open  = std::min(i,j) - 1;
        const int index_bracket_close = std::max(i,j) + 1;
        #ifdef DEBUG_GETSIMPLERNEWICKSFREQUENCYPAIRS
        std::clog
          << "["
          << index_bracket_open
          << "]-["
          << index_bracket_close
          << "]\n";
        #endif
        if ( new_newick_with_zero[index_bracket_open]  == Newick::bracket_open
          && new_newick_with_zero[index_bracket_close] == Newick::bracket_close)
        {
          new_newick_with_zero[index_bracket_open]  = 0;
          new_newick_with_zero[index_bracket_close] = 0;
          #ifdef DEBUG_GETSIMPLERNEWICKSFREQUENCYPAIRS
          std::clog << "2.5: " << Newick::DumbNewickToString(new_newick_with_zero) << '\n';
          #endif
        }
      }
      #ifdef DEBUG_GETSIMPLERNEWICKS
      std::clog << "2: " << Newick::DumbNewickToString(new_newick_with_zero) << '\n';
      #endif
      //Remove decremented i and possibly nullified brackets
      std::vector<int> new_newick;
      std::remove_copy(
        new_newick_with_zero.begin(),
        new_newick_with_zero.end(),
        std::back_inserter(new_newick),
        0);
      //Add brackets if these are removed
      if (new_newick.front() != Newick::bracket_open
        || new_newick.back() != Newick::bracket_close)
      {
        new_newick = Surround(new_newick);
      }
      #ifdef DEBUG_GETSIMPLERNEWICKSFREQUENCYPAIRS
      std::clog << "Store: " << Newick::DumbNewickToString(new_newick) << '\n';
      #endif
      assert(IsNewick(new_newick));
      assert(n[i] == 1);
      newicks.push_back(std::make_pair(new_newick,1));
      //newicks.push_back(new_newick);
      continue;
    }
  }
  */
  #else
  return NewickCpp98::GetSimplerNewicksFrequencyPairs(n);
  #endif
}
//---------------------------------------------------------------------------

///GetSimplerBinaryNewicksFrequencyPairs creates simpler, derived Newicks from a
///binary Newick as well as the frequency that is simplified.
///From http://www.richelbilderbeek.nl/CppGetSimplerBinaryNewicksFrequencyPairs.htm
const std::vector<std::pair<std::vector<int>,int> >
  GetSimplerBinaryNewicksFrequencyPairs(
  const std::vector<int>& n)
{
  assert(IsNewick(n));
  assert(IsBinaryNewick(n));
  std::vector<std::pair<std::vector<int>,int> > v;

  assert(IsNewick(n));
  assert(IsBinaryNewick(n));

  //If newick is simple (by counting the number of opening brackets),
  //there are no simpler Newicks
  if (std::count( n.begin(),n.end(),
    static_cast<int>(bracket_open))==1)
  {
    //Simple newicks do not need to be simplified
    assert(n.size()==3 || n.size() == 4);
    assert(n[0]==bracket_open);
    assert(n[n.size()-1]==bracket_close);
    if (n.size() == 3)
    {
      if (n[1]>1)
      {
        std::vector<int> next(n);
        --next[1];
        assert(IsNewick(next));
        v.push_back(std::make_pair(next,n[1]));
      }
      return v;
    }
    assert(n.size()==4);
    if (n[1] == 1)
    {
      std::vector<int> next
        = CreateVector(
            static_cast<int>(bracket_open),
            n[2]+1,
            static_cast<int>(bracket_close));
      assert(IsNewick(next));
      v.push_back(std::make_pair(next,1));
    }
    else
    {
      std::vector<int> next(n);
      --next[1];
      assert(IsNewick(next));
      v.push_back(std::make_pair(next,n[1]));
    }
    if (n[2] == 1)
    {
      std::vector<int> next
        = CreateVector(
            static_cast<int>(bracket_open),
            n[1]+1,
            static_cast<int>(bracket_close));
      assert(IsNewick(next));
      v.push_back(std::make_pair(next,1));
    }
    else
    {
      std::vector<int> next(n);
      --next[2];
      assert(IsNewick(next));
      v.push_back(std::make_pair(next,n[2]));
    }
    return v;
  }

  //newick is complex
  //Generate other Newicks and their coefficients
  const int sz = n.size();
  for (int i=0; i!=sz; ++i)
  {
    const int x = n[i];
    if (x < 0) //x is not a digit
    {
      continue;
    }
    if (x == 1)
    {
      //If x is not next to a digit, there is no simplification
      if (n[i-1]<0 && n[i+1]<1) continue;
      //If next to the x in a digit, merge these and remove their brackets
      //Is the 1 left of a value?
      if (n[i-1]==bracket_open)
      {
        assert(n[i+1] > 0);
        const int new_value = n[i+1] + 1;
        std::vector<int> next(n.begin(),n.begin() + i - 1);
        next.push_back(new_value);
        assert(n[i+2] < 0);
        std::copy(n.begin() + i + 3,n.end(),std::back_inserter(next));
        #ifndef NDEBUG
        if (!IsNewick(next))
        {
          std::cerr << DumbNewickToString(n)  << '\n';
          std::cerr << DumbNewickToString(next) << '\n';
          InspectInvalidNewick(std::cerr,next);
        }
        #endif
        assert(IsNewick(next));
        v.push_back(std::make_pair(next,x));
      }
      else
      {
        //Is the 1 to the right of a value?
        assert(n[i-1] > 0); //< The other value
        const int new_value = n[i-1] + 1;
        std::vector<int> next(n.begin(),n.begin()+i-2);
        next.push_back(new_value);
        assert(n[i+1] < 0);
        std::copy(n.begin() + i + 2,n.end(),std::back_inserter(next));
        assert(IsNewick(next));
        v.push_back(std::make_pair(next,x));
      }
    }
    else
    {
      std::vector<int> next = n;
      --next[i];
      assert(IsNewick(next));
      v.push_back(std::make_pair(next,x));
    }
  }
  return v;
}
//---------------------------------------------------------------------------
///InspectInvalidNewick writes the cause of the Newick invalidity
///to the std::ostream.
///From http://www.richelbilderbeek.nl/CppInspectInvalidNewick.htm
void InspectInvalidNewick(std::ostream& os, const std::vector<int>& v)
{
  os << "InspectInvalidNewick on: "
    << DumbNewickToString(v) << '\n';
  try
  {
    CheckNewick(v);
  }
  catch (std::exception& e)
  {
    os << "Invalidity caused by: " << e.what() << '\n';
  }
}
//---------------------------------------------------------------------------
///IsNewick returns true if a std::string is a valid Newick
///and false otherwise.
///From http://www.richelbilderbeek.nl/CppIsNewick.htm
bool IsNewick(const std::string& s)
{
  try
  {
    CheckNewick(s);
  }
  catch (const std::exception& e)
  {
    return false;
  }
  return true;
}
//---------------------------------------------------------------------------
///IsSimple returns true if the Newick std::vector contains
///leaves only. For example, the Newick '(1,2,3)' is simple,
///the Newick '((1,2),3)' is not simple
///From http://www.richelbilderbeek.nl/CppIsNewick.htm
bool IsSimple(const std::vector<int>& v)
{
  assert(IsNewick(v));
  //A Newick is simple if it contains no '(' after the initial one
  return std::count(
    v.begin()+1,v.end(),
    static_cast<int>(bracket_open)
  ) == 0;
}
//---------------------------------------------------------------------------
///IsBinaryNewick checks if a Newick is a binary tree,
///that is: each node splits in two (not more) branches
///From http://www.richelbilderbeek.nl/CppIsBinaryNewick.htm
bool IsBinaryNewick(std::vector<int> v)
{
  assert(IsNewick(v));
  if (IsUnaryNewick(v)) return false;

  while (1)
  {
    const int sz = boost::numeric_cast<int>(v.size());
    //Only numbers?
    if (IsSimple(v))
    {
      //Binary Newick has size 4, for example '(1,2)'
      return sz == 4;
    }
    if (GetLeafMaxArity(v) > 2) return false;
    v = Newick::ReplaceLeave(v,42);
  }
}
//---------------------------------------------------------------------------
///IsNewick returns true if a std::vector<int> is a valid Newick
///and false otherwise.
///From http://www.richelbilderbeek.nl/CppIsNewick.htm
bool IsNewick(const std::vector<int>& v)
{
  try
  {
    CheckNewick(v);
  }
  catch (...)
  {
    return false;
  }
  return true;
}
//---------------------------------------------------------------------------
///IsTrinaryNewick checks if a Newick is a trinary tree,
///that is: each node splits in three or less branches
///From http://www.richelbilderbeek.nl/CppIsTrinaryNewick.htm
bool IsTrinaryNewick(std::vector<int> v)
{
  assert(IsNewick(v));
  if (IsUnaryNewick(v)) return false;
  if (IsBinaryNewick(v)) return false;

  bool trinarity_found = false;

  while (1)
  {
    const int sz = boost::numeric_cast<int>(v.size());
    //Only numbers?
    if (IsSimple(v))
    {
      //Ternary Newick has size 5, for example '(1,2,3)'
      return trinarity_found || sz == 5;
    }
    const int leaf_max_arity = GetLeafMaxArity(v);
    if (leaf_max_arity > 3) return false;
    if (leaf_max_arity == 3) trinarity_found = true;

    v = Newick::ReplaceLeave(v,42);
  }
}
//---------------------------------------------------------------------------
///IsUnaryNewick checks if a Newick is a unary tree,
///that is: there is only one node.
///From http://www.richelbilderbeek.nl/CppIsUnaryNewick.htm
bool IsUnaryNewick(const std::vector<int>& v)
{
  assert(IsNewick(v));
  return v.size() == 3
    && v[0] == Newick::bracket_open
    && v[1] >  0
    && v[2] == Newick::bracket_close;
}
//---------------------------------------------------------------------------
///NewickToString converts a Newick std::vector<int> to a
///standard-format std::string.
///From http://www.richelbilderbeek.nl/CppNewickToString.htm
const std::string NewickToString(const std::vector<int>& v)
{
  assert(v.size() > 2 && "A Newick must at least have one single value");
  assert(v[0] == bracket_open
    && "A std::vector<int> Newick must start with a bracket_open");
  assert(v[v.size() - 1] == bracket_close
    && "A std::vector<int> Newick must end with a bracket_close");
  std::string s;
  s.reserve(2 * v.size()); //Just a guess
  const int sz = v.size();
  for (int i=0; i!=sz; ++i)
  {
    const int x = v[i];
    if (x >= 0)
    {
      s+=boost::lexical_cast<std::string>(x);
      assert(i+1<sz && "Must not end with number");
      const int next = v[i+1];
      if (next > 0 || next == bracket_open)
      {
        s+=",";
      }
    }
    else if (x==bracket_open)
    {
      s+="(";
    }
    else if (x==bracket_close)
    {
      s+=")";
      //Final closing bracket?
      if (i+1==sz) break;
      const int next = v[i+1];
      if (next > 0 || next == bracket_open)
      {
        s+=",";
      }
    }
    else
    {
      assert(!"Should not get here"
        && "A std::vector<int> Newick must consist of brackets and values only");
    }
  }
  return s;
}
//---------------------------------------------------------------------------
///SortNewick orders a Newick is such a way
///that all opening brackets are at the left side.
///For example (1,(2,3)) becomes ((2,3),1)
/*
std::string SortNewick(const std::string& newick)
{
  assert(IsNewick(newick));
  //All leaves are 'cut' by replacing them with an x
  std::string s = newick;
  std::string n = "";
  //Find initial leaf and replace it with x
  {
    const boost::regex r("\\(\\d+,\\d+\\)");
    std::string::const_iterator start = s.begin();
    const std::string::const_iterator end = s.end();
    boost::match_results<std::string::const_iterator> what;
    boost::regex_search(start, end, what, r);
    n = what.str();
    s = boost::regex_replace(s,r,"x");
  }
  //When all leaves are cut, s == 'x'
  while (s!="x")
  {
    //Obtain leaf with x
    const boost::regex r("(\\(x,\\d+\\))|(\\(\\d+,x\\))");
    std::string::const_iterator start = s.begin();
    const std::string::const_iterator end = s.end();
    boost::match_results<std::string::const_iterator> what;
    //Search for inner leaf
    boost::regex_search(start, end, what, r);
    const std::string l = what.str();
    //Search leaf for digit
    boost::regex_search(l.begin(), l.end(), what,boost::regex("\\d+"));
    //Add digit to n
    n = "(" + n + "," + what.str() + ")";
    //Replace the leaf by an x
    s = boost::regex_replace(s,r,"x");
  }
  return n;
}
*/
//---------------------------------------------------------------------------
///ReplaceLeave replaces the first leaf that it finds by a value.
///For example, using ReplaceLeave on '((1,2),(3,4))' with a value
///of 42 results in '(42,(3,4))'.
const std::vector<int> ReplaceLeave(
  const std::vector<int>& newick,
  const int value)
{
  assert(IsNewick(newick) && "Only a valid Newick can have its leaves replaced");
  assert(!IsSimple(newick) && "There must a leaf to simplify");
  typedef std::vector<int>::const_iterator Iterator;
  const Iterator end = newick.end();
  for (Iterator from = newick.begin(); from!=end; ++from)
  {
    if (*from != Newick::bracket_open) continue;

    for (Iterator to = from + 1; to!=end; ++to)
    {
      if (*to > 0) continue;
      if (*to == Newick::bracket_open) break;
      if (*to == Newick::bracket_close)
      {
        //Found
        std::vector<int> new_newick(newick.begin(),from);
        new_newick.push_back(value);
        std::copy(to + 1,newick.end(),std::back_inserter(new_newick));
        assert(IsNewick(new_newick));
        return new_newick;
      }
    }
  }
  assert(!"Should not get here");
  throw std::logic_error("Should not get here");
}
//---------------------------------------------------------------------------
///StringToNewick converts a std::string to a Newick std::vector<int>
///StringToNewick assumes that the input is well-formed and
///has both trailing and tailing brackets.
///From http://www.richelbilderbeek.nl/CppNewickToVector.htm
const std::vector<int> StringToNewick(const std::string& newick)
{
  assert(IsNewick(newick));
  assert(!newick.empty()
    && "s must not be empty");
  assert(newick[              0]=='('
    && "s must begin with a '('");
  assert(newick[newick.size()-1]==')'
    && "s must end with a ')'");

  std::vector<int> v;
  int value = 0;

  BOOST_FOREACH(const char i,newick)
  {
    if (i == '(')
    {
      if (value!=0) v.push_back(value);
      value = 0;
      v.push_back(bracket_open);
      continue;
    }
    if (i == ')')
    {
      if (value!=0) v.push_back(value);
      value = 0;
      v.push_back(bracket_close);
      continue;
    }
    if (i == ',')
    {
      if (value!=0) v.push_back(value);
      value = 0;
      continue;
    }
    assert(i >= '0' && i <= '9'); //Should be a number
    value*=10;
    value+=boost::lexical_cast<int>(i);
  }
  assert(value == 0 && "Final bracket close must set value to zero");
  return v;
}
//---------------------------------------------------------------------------
///Surround surrounds the Newick with brackets
const std::vector<int> Surround(const std::vector<int>& newick)
{
  std::vector<int> new_newick;
  new_newick.push_back(Newick::bracket_open);
  std::copy(newick.begin(),newick.end(),std::back_inserter(new_newick));
  new_newick.push_back(Newick::bracket_close);
  return new_newick;
}
//---------------------------------------------------------------------------
///Test tests all Newick functions
void Test()
{
  //#ifdef __GXX_EXPERIMENTAL_CXX0X__
  ///\note
  ///The tests below must be put back in again once
  #ifdef DEBUG_TEMP_REMOVE_2738236826438


  TRACE("Testing basic Newick functionality");
  //Check difference between C++98 and C++0x
  assert(Newick::CreateValidTrinaryNewicks() == NewickCpp98::CreateValidTrinaryNewicks());
  assert(Newick::GetKnownProbabilities() == NewickCpp98::GetKnownProbabilities());



  //Check conversions from std::string to std::vector #1
  {
    const std::vector<int> v = Newick::StringToNewick("(11,(22,33))");
    assert(v.size() == 7);
    assert(v[0]==Newick::bracket_open);
    assert(v[1]==11);
    assert(v[2]==Newick::bracket_open);
    assert(v[3]==22);
    assert(v[4]==33);
    assert(v[5]==Newick::bracket_close);
    assert(v[6]==Newick::bracket_close);
  }
  //Check if well-formed Newicks are accepted
  {
    const std::vector<std::string> v = Newick::CreateValidNewicks();
    BOOST_FOREACH(const std::string& s,v)
    {
      #ifdef TRACE_REJECTED_NEWICKS
      const std::string debug = "I must be accepted: " + s;
      TRACE(debug);
      #endif
      assert(Newick::IsNewick(s));
      const std::vector<int> v = Newick::StringToNewick(s);
      assert(Newick::IsNewick(v));
    }
  }
  //Check if ill-formed Newicks are rejected
  {
    #ifndef NDEBUG
    const std::vector<std::string> v = Newick::CreateInvalidNewicks();
    BOOST_FOREACH(const std::string& s,v)
    {
      #ifdef TRACE_REJECTED_NEWICKS
      const std::string debug = "I must be rejected: " + s;
      TRACE(debug);
      #endif
      assert(!Newick::IsNewick(s));
      //Cannot test if std::vector<int> versions are rejected,
      //because Newick::StringToNewick assumes a valid Newick
      //const std::vector<int> v = Newick::StringToNewick(s);
      //assert(!Newick::IsNewick(v));
    }
    #endif
  }
  //Check conversions from std::string to std::vector #2
  {
    const std::vector<int> v = Newick::StringToNewick("((11,22),33)");
    assert(v.size() == 7);
    assert(v[0]==Newick::bracket_open);
    assert(v[1]==Newick::bracket_open);
    assert(v[2]==11);
    assert(v[3]==22);
    assert(v[4]==Newick::bracket_close);
    assert(v[5]==33);
    assert(v[6]==Newick::bracket_close);
  }
  assert(Newick::CalcNumOfSymmetriesBinary(StringToNewick("(1,(3,1))"))==0);
  assert(Newick::CalcNumOfSymmetriesBinary(StringToNewick("(3,(1,1))"))==1);
  assert(Newick::CalcNumOfSymmetriesBinary(StringToNewick("(1,((1,1),(1,1)))"))==3);
  assert(Newick::CalcNumOfSymmetriesBinary(StringToNewick("(1,((1,1),(2,2)))"))==2);
  assert(Newick::CalcNumOfSymmetriesBinary(StringToNewick("(1,(2,3))"))==0);
  assert(Newick::CalcNumOfSymmetriesBinary(StringToNewick("(99,99)"))==1);
  assert(Newick::CalcNumOfSymmetriesBinary(StringToNewick("(3,(2,2))"))==1);
  assert(Newick::CalcNumOfSymmetriesBinary(StringToNewick("(2,(2,2))"))==1);
  assert(Newick::CalcNumOfSymmetriesBinary(StringToNewick("((3,3),(2,2))"))==2);
  assert(Newick::CalcNumOfSymmetriesBinary(StringToNewick("((3,3),(3,3))"))==3);
  assert(Newick::CalcNumOfSymmetriesBinary(StringToNewick("((3,3),(3,4))"))==1);
  assert(Newick::CalcNumOfSymmetriesBinary(StringToNewick("(((3,3),(4,4)),5)"))==2);
  assert(Newick::CalcNumOfSymmetriesBinary(StringToNewick("(((3,3),(5,5)),5)"))==2);
  assert(Newick::CalcNumOfSymmetriesBinary(StringToNewick("(((5,5),(5,5)),5)"))==3);
  assert(Newick::CalcNumOfSymmetriesBinary(StringToNewick("(((5,5),(5,5)),(4,4))"))==4);
  assert(Newick::CalcNumOfSymmetriesBinary(StringToNewick("(((5,5),(4,4)),(4,4))"))==3);
  assert(Newick::CalcNumOfSymmetriesBinary(StringToNewick("(((4,4),(4,4)),(4,4))"))==4);
  assert(Newick::CalcNumOfCombinationsBinary(StringToNewick("(3,(1,1))"))==10);
  assert(Newick::CalcNumOfCombinationsBinary(StringToNewick("(1,(3,1))"))==20);
  assert(Newick::CalcNumOfCombinationsBinary(StringToNewick("(1,(1,(1,(1,1))))"))==60);
  assert(Newick::CalcNumOfCombinationsBinary(StringToNewick("(1,((1,1),(1,1)))"))==15);
  assert(bigIntegerToString(Newick::FactorialBigInt(1))=="1");
  assert(bigIntegerToString(Newick::FactorialBigInt(2))=="2");
  assert(bigIntegerToString(Newick::FactorialBigInt(3))=="6");
  assert(bigIntegerToString(Newick::FactorialBigInt(4))=="24");
  assert(bigIntegerToString(Newick::FactorialBigInt(5))=="120");
  assert(bigIntegerToString(Newick::FactorialBigInt(6))=="720");
  assert(Newick::GetLeafMaxArity(Newick::StringToNewick("(1)"))   == 1);
  assert(Newick::GetLeafMaxArity(Newick::StringToNewick("(12)"))  == 1);
  assert(Newick::GetLeafMaxArity(Newick::StringToNewick("(123)")) == 1);
  assert(Newick::GetLeafMaxArity(Newick::StringToNewick("(1,2)"))   == 2);
  assert(Newick::GetLeafMaxArity(Newick::StringToNewick("(12,2)"))  == 2);
  assert(Newick::GetLeafMaxArity(Newick::StringToNewick("(123,2)")) == 2);
  assert(Newick::GetLeafMaxArity(Newick::StringToNewick("(1,(1,2))"))   == 2);
  assert(Newick::GetLeafMaxArity(Newick::StringToNewick("(1,(12,2))"))  == 2);
  assert(Newick::GetLeafMaxArity(Newick::StringToNewick("(1,(123,2))")) == 2);
  assert(Newick::GetLeafMaxArity(Newick::StringToNewick("((1,2),3)"))   == 2);
  assert(Newick::GetLeafMaxArity(Newick::StringToNewick("((12,2),3)"))  == 2);
  assert(Newick::GetLeafMaxArity(Newick::StringToNewick("((123,2),3)")) == 2);
  assert(Newick::GetLeafMaxArity(Newick::StringToNewick("(1,2,3)"))   == 3);
  assert(Newick::GetLeafMaxArity(Newick::StringToNewick("(12,2,3)"))  == 3);
  assert(Newick::GetLeafMaxArity(Newick::StringToNewick("(123,2,3)")) == 3);
  assert(Newick::GetLeafMaxArity(Newick::StringToNewick("(1,(1,2,3))"))   == 3);
  assert(Newick::GetLeafMaxArity(Newick::StringToNewick("(1,(12,2,3))"))  == 3);
  assert(Newick::GetLeafMaxArity(Newick::StringToNewick("(1,(123,2,3))")) == 3);
  assert(Newick::GetLeafMaxArity(Newick::StringToNewick("((1,2,3),4)"))   == 3);
  assert(Newick::GetLeafMaxArity(Newick::StringToNewick("((12,2,3),4)"))  == 3);
  assert(Newick::GetLeafMaxArity(Newick::StringToNewick("((123,2,3),4)")) == 3);

  assert(fuzzy_equal_to()(  2.0,Newick::CalcDenominator(Newick::StringToNewick("(1,1)"),10.0)));
  assert(fuzzy_equal_to()(  6.0,Newick::CalcDenominator(Newick::StringToNewick("((1,1),1)"),10.0)));
  assert(fuzzy_equal_to()( 26.0,Newick::CalcDenominator(Newick::StringToNewick("(1,2)"),10.0)));
  assert(fuzzy_equal_to()( 32.0,Newick::CalcDenominator(Newick::StringToNewick("((1,1),2)"),10.0)));
  assert(fuzzy_equal_to()( 32.0,Newick::CalcDenominator(Newick::StringToNewick("(2,(1,1))"),10.0)));
  assert(fuzzy_equal_to()( 50.0,Newick::CalcDenominator(Newick::StringToNewick("((1,1),3)"),10.0)));
  assert(fuzzy_equal_to()( 80.0,Newick::CalcDenominator(Newick::StringToNewick("((1,2),3)"),10.0)));
  assert(fuzzy_equal_to()( 80.0,Newick::CalcDenominator(Newick::StringToNewick("((3,1),2)"),10.0)));
  assert(fuzzy_equal_to()( 80.0,Newick::CalcDenominator(Newick::StringToNewick("((2,3),1)"),10.0)));
  assert(fuzzy_equal_to()(102.0,Newick::CalcDenominator(Newick::StringToNewick("((2,1),4)"),10.0)));
  assert(fuzzy_equal_to()(152.0,Newick::CalcDenominator(Newick::StringToNewick("(2,(1,(3,3)))"),10.0)));
  assert(fuzzy_equal_to()(162.0,Newick::CalcDenominator(Newick::StringToNewick("((2,3),4)"),10.0)));
  assert(fuzzy_equal_to()(180.0,Newick::CalcDenominator(Newick::StringToNewick("((1,2),(3,4))"),10.0)));
  assert(fuzzy_equal_to()(180.0,Newick::CalcDenominator(Newick::StringToNewick("((4,1),(2,3))"),10.0)));
  assert(fuzzy_equal_to()(180.0,Newick::CalcDenominator(Newick::StringToNewick("((3,4),(1,2))"),10.0)));
  assert(fuzzy_equal_to()(180.0,Newick::CalcDenominator(Newick::StringToNewick("((2,3),(4,1))"),10.0)));
  {
    const std::vector<int> v = { 0,1,2,3,4,5,6 };
    assert(Newick::FindPosAfter(v,3,2)==3);
    assert(Newick::FindPosAfter(v,4,2)==4);
    assert(Newick::FindPosAfter(v,5,2)==5);
    assert(Newick::FindPosAfter(v,6,2)==6);
    assert(Newick::FindPosBefore(v,3,4)==3);
    assert(Newick::FindPosBefore(v,2,4)==2);
    assert(Newick::FindPosBefore(v,1,4)==1);
    assert(Newick::FindPosBefore(v,0,4)==0);
  }
  {
    const std::vector<int> v = Newick::GetDepth(Newick::StringToNewick("(1,(2,2))"));
    const std::vector<int> w = Newick::GetDepth(Newick::StringToNewick("(9,(9,9))"));
    const std::vector<int> x = { 0,0,1,1,1,1,0 };
    assert(v == x);
    assert(w == x);
  }
  {
    const std::vector<int> v = Newick::GetDepth(Newick::StringToNewick("((2,2),1)"));
    const std::vector<int> w = { 0,1,1,1,1,0,0 };
    assert(v == w);
  }
  {
    const std::vector<int> v = Newick::GetDepth(Newick::StringToNewick("(1,(2,2),1)"));
    const std::vector<int> w = { 0,0,1,1,1,1,0,0 };
    assert(v == w);
  }
  {
    const std::vector<int> v = Newick::GetDepth(Newick::StringToNewick("(1,(2,3),4,(5,6))"));
    const std::vector<int> w = { 0,0,1,1,1,1,0,1,1,1,1,0 };
    assert(v == w);
  }
  {
    const std::vector<int> v = Newick::GetDepth(Newick::StringToNewick("(1,(2,3),(5,6))"));
    const std::vector<int> w = { 0,0,1,1,1,1,1,1,1,1,0 };
    assert(v == w);
  }
  {
    const std::vector<int> v = Newick::GetDepth(Newick::StringToNewick("(1,(2,(3,4)),((5,6),7))"));
    const std::vector<int> w = { 0,0,1,1,2,2,2,2,1,1,2,2,2,2,1,1,0 };
    assert(v == w);
  }
  //Test GetRootBranches
  {
    const std::vector<std::vector<int> > v = Newick::GetRootBranches(Newick::StringToNewick("(1,2)"));
    assert(v.size() == 2);
    assert(std::find(v.begin(),v.end(),
      Newick::StringToNewick("(1)")) != v.end());
    assert(std::find(v.begin(),v.end(),
      Newick::StringToNewick("(2)")) != v.end());
  }
  {
    const std::vector<std::vector<int> > v = Newick::GetRootBranches(Newick::StringToNewick("(1,(2,3))"));
    assert(v.size() == 2);
    assert(std::find(v.begin(),v.end(),
      Newick::StringToNewick("(1)")) != v.end());
    assert(std::find(v.begin(),v.end(),
      Newick::StringToNewick("(2,3)")) != v.end());
  }
  {
    const std::vector<std::vector<int> > v = Newick::GetRootBranches(Newick::StringToNewick("(1,2,(3,4))"));
    assert(v.size() == 3);
    assert(std::find(v.begin(),v.end(),
      Newick::StringToNewick("(1)")) != v.end());
    assert(std::find(v.begin(),v.end(),
      Newick::StringToNewick("(2)")) != v.end());
    assert(std::find(v.begin(),v.end(),
      Newick::StringToNewick("(3,4)")) != v.end());
  }
  //Compare C++98 and C++0x version
  {
    const std::vector<std::string> v = Newick::CreateValidBinaryNewicks();
    BOOST_FOREACH(const std::string& s,v)
    {
      const std::vector<int> n = Newick::StringToNewick(s);
      assert(Newick::GetRootBranches(n) == NewickCpp98::GetRootBranches(n));
    }
  }

  //Check if binary and trinary Newicks are detected correctly
  {
    const std::vector<std::string> v = Newick::CreateValidBinaryNewicks();
    BOOST_FOREACH(const std::string& s,v)
    {
      const std::vector<int> n = Newick::StringToNewick(s);
      assert(Newick::IsBinaryNewick(n));
    }
  }
  //Check if unary Newicks are detected correctly
  {
    const std::vector<std::string> v = Newick::CreateValidUnaryNewicks();
    BOOST_FOREACH(const std::string& s,v)
    {
      const std::vector<int> n = Newick::StringToNewick(s);
      assert( Newick::GetLeafMaxArity(n)<=1);
      assert( Newick::IsUnaryNewick(n));
      assert(!Newick::IsBinaryNewick(n));
      assert(!Newick::IsTrinaryNewick(n));
    }
  }
  //Check if binary Newicks are detected correctly
  {
    const std::vector<std::string> v = Newick::CreateValidBinaryNewicks();
    BOOST_FOREACH(const std::string& s,v)
    {
      const std::vector<int> n = Newick::StringToNewick(s);
      assert( Newick::GetLeafMaxArity(n)<=2);
      assert(!Newick::IsUnaryNewick(n));
      assert( Newick::IsBinaryNewick(n));
      assert(!Newick::IsTrinaryNewick(n));
    }
  }
  //Check if trinary Newicks are detected correctly
  {
    const std::vector<std::string> v = Newick::CreateValidTrinaryNewicks();
    BOOST_FOREACH(const std::string& s,v)
    {
      //TRACE(s);
      const std::vector<int> n = Newick::StringToNewick(s);
      assert( Newick::GetLeafMaxArity(n)<=3);
      assert(!Newick::IsUnaryNewick(n));
      assert(!Newick::IsBinaryNewick(n));
      assert( Newick::IsTrinaryNewick(n));
    }
  }
  //Test binary Newick
  {
    const std::string s("(1,(2,3))");
    const std::vector<std::vector<int> > n = GetSimplerNewicks(Newick::StringToNewick(s));
    //#define DEBUG_1_BO_1_2_3_BC
    #ifdef  DEBUG_1_BO_1_2_3_BC
    BOOST_FOREACH(const auto& t,n)
    {
      TRACE(Newick::NewickToString(t));
    }
    #endif
    assert(n.size() == 2);
    assert(std::find(n.begin(),n.end(),Newick::StringToNewick("(1,(1,3))"))
      != n.end());
    assert(std::find(n.begin(),n.end(),Newick::StringToNewick("(1,(2,2))"))
      != n.end());
  }
  {
    const std::string s("(1,(2,3,4))");
    const std::vector<std::vector<int> > n = GetSimplerNewicks(Newick::StringToNewick(s));
    assert(n.size() == 3);
    assert(std::find(n.begin(),n.end(),Newick::StringToNewick("(1,(1,3,4))"))
      != n.end());
    assert(std::find(n.begin(),n.end(),Newick::StringToNewick("(1,(2,2,4))"))
      != n.end());
    assert(std::find(n.begin(),n.end(),Newick::StringToNewick("(1,(2,3,3))"))
      != n.end());
  }
  {
    const std::string s("(1,(1,3,4))");
    const std::vector<std::vector<int> > n = GetSimplerNewicks(Newick::StringToNewick(s));
    //#define DEBUG_1_BO_1_3_4_BC
    #ifdef  DEBUG_1_BO_1_3_4_BC
    TRACE(boost::lexical_cast<std::string>(n.size()));
    BOOST_FOREACH(const auto& t,n)
    {
      TRACE(Newick::NewickToString(t));
    }
    #endif
    assert(n.size() == 4);
    assert(std::find(n.begin(),n.end(),Newick::StringToNewick("(1,(4,4))"))
      != n.end());
    assert(std::find(n.begin(),n.end(),Newick::StringToNewick("(1,(3,5))"))
      != n.end());
    assert(std::find(n.begin(),n.end(),Newick::StringToNewick("(1,(1,2,4))"))
      != n.end());
    assert(std::find(n.begin(),n.end(),Newick::StringToNewick("(1,(1,3,3))"))
      != n.end());
  }
  {
    const std::string s("(1,(1,3,4))");
    const std::vector<std::pair<std::vector<int>,int> > n
      = GetSimplerNewicksFrequencyPairs(Newick::StringToNewick(s));
    typedef std::pair<std::vector<int>,int> Pair;
    #ifdef TRACE_GETSIMPLERNEWICKSFREQUENCYPAIRS_1_134
    BOOST_FOREACH(const Pair& p, n)
    {
      std::cout << Newick::NewickToString(p.first) << '\n';
    }
    #endif
    assert(n.size() == 4);
    assert(std::find(n.begin(),n.end(),
      std::make_pair(Newick::StringToNewick("(1,(4,4))"),1))
      != n.end());
    assert(std::find(n.begin(),n.end(),
      std::make_pair(Newick::StringToNewick("(1,(3,5))"),1))
      != n.end());
    assert(std::find(n.begin(),n.end(),
      std::make_pair(Newick::StringToNewick("(1,(1,2,4))"),3))
      != n.end());
    assert(std::find(n.begin(),n.end(),
      std::make_pair(Newick::StringToNewick("(1,(1,3,3))"),4))
      != n.end());
  }
  {
    const std::string s("((1,1),2)");
    const std::vector<std::vector<int> > n = Newick::GetSimplerNewicks(
      Newick::StringToNewick(s));
    //#define DEBUG_BO_1_1_BC_2
    #ifdef  DEBUG_BO_1_1_BC_2
    BOOST_FOREACH(const auto& t,n)
    {
      TRACE(Newick::NewickToString(t));
    }
    #endif
    assert(n.size() == 3);
    assert(std::find(n.begin(),n.end(),
      Newick::StringToNewick("(2,2)"))
      != n.end());
    assert(std::find(n.begin(),n.end(),
      Newick::StringToNewick("((1,1),1)"))
      != n.end());
  }
  {
    const std::string s("((1,1),2)");
    typedef std::pair<std::vector<int>,int> Pair;
    const std::vector<Pair> n
      = GetSimplerNewicksFrequencyPairs(Newick::StringToNewick(s));
    #ifdef TRACE_GETSIMPLERNEWICKSFREQUENCYPAIRS_11_2
    BOOST_FOREACH(const Pair& p, n)
    {
      std::clog << Newick::NewickToString(p.first) << '\n';
    }
    #endif
    assert(n.size() == 3);
    assert(std::find(n.begin(),n.end(),
      std::make_pair(Newick::StringToNewick("(2,2)"),1))
      != n.end());
    assert(std::find(n.begin(),n.end(),
      std::make_pair(Newick::StringToNewick("((1,1),1)"),2))
      != n.end());
  }
  {
    const std::string s("((2,1),4)");
    const std::vector<std::vector<int> > n = Newick::GetSimplerNewicks(
      Newick::StringToNewick(s));
    assert(n.size() == 3);
    assert(std::find(n.begin(),n.end(),
      Newick::StringToNewick("(3,4)"))
      != n.end());
    assert(std::find(n.begin(),n.end(),
      Newick::StringToNewick("((1,1),4)"))
      != n.end());
    assert(std::find(n.begin(),n.end(),
      Newick::StringToNewick("((2,1),3)"))
      != n.end());
  }
  {
    const std::string s("((2,1),4)");
    typedef std::pair<std::vector<int>,int> Pair;
    const std::vector<Pair> n
      = GetSimplerNewicksFrequencyPairs(Newick::StringToNewick(s));
    #ifdef TRACE_GETSIMPLERNEWICKSFREQUENCYPAIRS_21_2
    BOOST_FOREACH(const Pair& p, n)
    {
      TRACE(Newick::NewickToString(p.first));
    }
    #endif
    assert(n.size() == 3);
    assert(std::find(n.begin(),n.end(),
      std::make_pair(Newick::StringToNewick("(3,4)"),1))
      != n.end());
    assert(std::find(n.begin(),n.end(),
      std::make_pair(Newick::StringToNewick("((1,1),4)"),2))
      != n.end());
    assert(std::find(n.begin(),n.end(),
      std::make_pair(Newick::StringToNewick("((2,1),3)"),4))
      != n.end());
  }
  {
    const std::string s("((2,3),4)");
    const std::vector<std::vector<int> > n = Newick::GetSimplerNewicks(
      Newick::StringToNewick(s));
    assert(n.size() == 3);
    assert(std::find(n.begin(),n.end(),
      Newick::StringToNewick("((1,3),4)"))
      != n.end());
    assert(std::find(n.begin(),n.end(),
      Newick::StringToNewick("((2,2),4)"))
      != n.end());
    assert(std::find(n.begin(),n.end(),
      Newick::StringToNewick("((2,3),3)"))
      != n.end());
  }
  {
    const std::string s("((2,3),4)");
    typedef std::pair<std::vector<int>,int> Pair;
    const std::vector<Pair> n
      = GetSimplerNewicksFrequencyPairs(Newick::StringToNewick(s));
    #ifdef TRACE_GETSIMPLERNEWICKSFREQUENCYPAIRS_23_4
    BOOST_FOREACH(const Pair& p, n)
    {
      std::cout << Newick::NewickToString(p.first) << '\n';
    }
    #endif
    assert(n.size() == 3);
    assert(std::find(n.begin(),n.end(),
      std::make_pair(Newick::StringToNewick("((1,3),4)"),2))
      != n.end());
    assert(std::find(n.begin(),n.end(),
      std::make_pair(Newick::StringToNewick("((2,2),4)"),3))
      != n.end());
    assert(std::find(n.begin(),n.end(),
      std::make_pair(Newick::StringToNewick("((2,3),3)"),4))
      != n.end());
  }
  //Compare GetSimplerNewicks and
  //GetSimplerNewicksFrequencyPairs
  {
    const std::vector<std::string> newicks
      = Newick::CreateValidNewicks();
    BOOST_FOREACH(const std::string& newick_str,newicks)
    {
      const std::vector<int> newick
        = Newick::StringToNewick(newick_str);
      const std::vector<std::vector<int> > v1
        = Newick::GetSimplerNewicks(newick);
      const std::vector<std::pair<std::vector<int>,int> > v2
        = Newick::GetSimplerNewicksFrequencyPairs(newick);
      assert(v1.size() == v2.size());
      const int size = boost::numeric_cast<int>(v1.size());
      for (int i=0; i!=size; ++i)
      {
        #define DEBUG_COMPARE_GSN_VS_GSNFP
        #ifdef  DEBUG_COMPARE_GSN_VS_GSNFP
        if (v1[i] != v2[i].first)
        {
          TRACE("ERROR: DIFFERENT NEWICK SIMPLIFICATIONS");
          TRACE(Newick::NewickToString(newick));
          TRACE(Newick::NewickToString(v1[i]));
          TRACE(Newick::NewickToString(v2[i].first));
        }
        #endif
        assert(v1[i] == v2[i].first);
      }
      assert(Newick::GetSimplerNewicksFrequencyPairs(newick)
        == NewickCpp98::GetSimplerNewicksFrequencyPairs(newick));
    }
  }
  #else
  NewickCpp98::Test();
  #endif
}
//---------------------------------------------------------------------------
} //namespace Newick
//---------------------------------------------------------------------------

 

 

 

 

 

newickcpp98.h

 

//---------------------------------------------------------------------------
/*
NewickCpp98, C++98 Newick functions
Copyright (C) 2010-2011 Richel Bilderbeek

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program.If not, see <http://www.gnu.org/licenses/>.
*/
//---------------------------------------------------------------------------
//From http://www.richelbilderbeek.nl/CppNewick.htm
//---------------------------------------------------------------------------
#ifndef NEWICKCPP98_H
#define NEWICKCPP98_H
//---------------------------------------------------------------------------
#ifdef _WIN32
//See http://www.richelbilderbeek.nl/CppCompileErrorSwprintfHasNotBeenDeclared.htm
#undef __STRICT_ANSI__
#endif
//---------------------------------------------------------------------------
#include <string>
#include <vector>
//---------------------------------------------------------------------------
#include <boost/tuple/tuple.hpp>
//---------------------------------------------------------------------------
namespace NewickCpp98 {
//---------------------------------------------------------------------------
//Functions that do not use the C++0x standard
const std::vector<std::string> CreateValidTrinaryNewicks();
std::vector<boost::tuple<std::string,double,double> > GetKnownProbabilities();
const std::vector<std::pair<std::vector<int>,int> > GetSimplerNewicksFrequencyPairs(const std::vector<int>& n);
const std::vector<std::vector<int> > GetRootBranches(const std::vector<int>& n);
void Test();
//---------------------------------------------------------------------------
} //~namespace NewickCpp98
//---------------------------------------------------------------------------
#endif // NEWICKCPP98_H

 

 

 

 

 

newickcpp98.cpp

 

//---------------------------------------------------------------------------
#include <boost/numeric/conversion/cast.hpp>
#include <boost/tuple/tuple.hpp>
#include <boost/tuple/tuple_comparison.hpp>
//---------------------------------------------------------------------------
#include "fuzzy_equal_to.h"
#include "newick.h"
#include "newickcpp98.h"
//---------------------------------------------------------------------------
namespace NewickCpp98 {
//---------------------------------------------------------------------------
///CreateValidTrinaryNewicks creates std::strings
///that can be converted to a TrinaryNewickVector.
///From http://www.richelbilderbeek.nl/CppCreateValidTinaryNewicks.htm
const std::vector<std::string> CreateValidTrinaryNewicks()
{
  std::vector<std::string> v;
  v.push_back("(1,1,1)");
  v.push_back("(1,2,3)");
  v.push_back("((1,1),1,1)");
  v.push_back("(1,(1,1),1)");
  v.push_back("(1,1,(1,1))");
  v.push_back("(1,(2,3,4))");
  v.push_back("(1,2,(3,4))");
  v.push_back("(1,2,(3,4,5))");
  v.push_back("((1,2,3),4,5)");
  v.push_back("(11,22,33)");
  v.push_back("(11,(22,33,44))");
  v.push_back("(11,22,(33,44))");
  v.push_back("(11,22,(33,44,55))");
  v.push_back("((11,22,33),44,55)");
  v.push_back("((1,2),(3,4),(5,6))");
  v.push_back("((1,2,3),(4,5),(6,7))");
  v.push_back("((1,2),(3,4,5),(6,7))");
  v.push_back("((1,2),(3,4),(5,6,7))");
  v.push_back("((1,2,3),(4,5),(6,7))");
  v.push_back("((1,2),(3,4,5),(6,7))");
  v.push_back("((1,2),(3,4),(5,6,7))");
  v.push_back("((1,2,3),(4,5,6),(7,8))");
  v.push_back("((1,2),(3,4,5),(6,7,8))");
  v.push_back("((1,2,3),(4,5),(6,7,8))");
  v.push_back("((1,2,3),(4,5,6),(7,8,9))");
  v.push_back("((11,22,33),(44,55,66),(77,88,99))");
  return v;
}
//---------------------------------------------------------------------------
std::vector<boost::tuple<std::string,double,double> > GetKnownProbabilities()
{
  std::vector<boost::tuple<std::string,double,double> > v;

  //Sum equals 1
  v.push_back(boost::make_tuple("(1)"  , 10.0, 1.0000000));
  //Sum equals 2
  v.push_back(boost::make_tuple("(2)"  , 10.0, 0.0909091));
  v.push_back(boost::make_tuple("(1,1)", 10.0, 0.9090909));
  //Sum equals 3
  v.push_back(boost::make_tuple("(3)"  , 10.0, 0.0151515));
  v.push_back(boost::make_tuple("(1,2)", 10.0, 0.0757576));
  v.push_back(boost::make_tuple("(2,1)", 10.0, 0.0757576));
  v.push_back(boost::make_tuple("(1,(1,1))", 10.0, 0.2525253));
  v.push_back(boost::make_tuple("((1,1),1)", 10.0, 0.2525253));
  //Trinary
  v.push_back(boost::make_tuple("(1,1,1)"  , 10.0, 0.7575758));
  //Sum equals 4
  v.push_back(boost::make_tuple("(4)"  , 10.0, 0.0034965));
  v.push_back(boost::make_tuple("(1,3)", 10.0, 0.0116550));
  v.push_back(boost::make_tuple("(2,2)", 10.0, 0.0058275));
  v.push_back(boost::make_tuple("(3,1)", 10.0, 0.0116550));
  v.push_back(boost::make_tuple("(1,(1,2))", 10.0, 0.0194250));
  v.push_back(boost::make_tuple("(1,(2,1))", 10.0, 0.0194250));
  v.push_back(boost::make_tuple("(2,(1,1))", 10.0, 0.0194250));
  v.push_back(boost::make_tuple("((1,2),1)", 10.0, 0.0194250));
  v.push_back(boost::make_tuple("((2,1),1)", 10.0, 0.0194250));
  v.push_back(boost::make_tuple("((1,1),2)", 10.0, 0.0194250));
  //Trinary
  v.push_back(boost::make_tuple("(1,1,2)", 10.0, 0.0582751));
  v.push_back(boost::make_tuple("(1,2,1)", 10.0, 0.0582751));
  v.push_back(boost::make_tuple("(2,1,1)", 10.0, 0.0582751));
  v.push_back(boost::make_tuple("(1,1,(1,1))", 10.0, 0.1295001));   //(1)(confirmed)
  v.push_back(boost::make_tuple("(1,(1,1),1)", 10.0, 0.1295001));
  v.push_back(boost::make_tuple("((1,1),1,1)", 10.0, 0.1295001));
  v.push_back(boost::make_tuple("(1,(1,1,1))", 10.0, 0.0971251));   //(2)(confirmed)
  v.push_back(boost::make_tuple("((1,1,1),1)", 10.0, 0.0971251));
  //Quadrary
  v.push_back(boost::make_tuple("(1,1,1,1)", 10.0, 0.5827505));
  //Sum equals 5
  v.push_back(boost::make_tuple("(1,4)", 10.0, 0.0024975));
  v.push_back(boost::make_tuple("(2,3)", 10.0, 0.0008325));
  v.push_back(boost::make_tuple("(3,2)", 10.0, 0.0008325));
  v.push_back(boost::make_tuple("(4,1)", 10.0, 0.0024975));
  v.push_back(boost::make_tuple("(1,(1,3))", 10.0, 0.0028305));
  v.push_back(boost::make_tuple("(1,(2,2))", 10.0, 0.0012950));
  v.push_back(boost::make_tuple("(1,(3,1))", 10.0, 0.0028305));
  v.push_back(boost::make_tuple("(2,(1,2))", 10.0, 0.0014338));
  v.push_back(boost::make_tuple("(2,(2,1))", 10.0, 0.0014338));
  v.push_back(boost::make_tuple("(3,(1,1))", 10.0, 0.0026640));
  //Trinary
  v.push_back(boost::make_tuple("(1,1,(1,2))"  , 10.0, 0.0092731));   //(3)(confirmed)
  v.push_back(boost::make_tuple("(1,1,(2,1))"  , 10.0, 0.0092731));
  v.push_back(boost::make_tuple("(1,1,(1,1,1))", 10.0, 0.0348263));   //(4)(confirmed)
  v.push_back(boost::make_tuple("(1,(1,1,1),1)", 10.0, 0.0348263));
  v.push_back(boost::make_tuple("((1,1,1),1,1)", 10.0, 0.0348263));
  v.push_back(boost::make_tuple("(2,(1,1,1))"  , 10.0, 0.0070069));   //(5)(confirmed)
  v.push_back(boost::make_tuple("((1,1,1),2)"  , 10.0, 0.0070069));
  v.push_back(boost::make_tuple("(1,1,1,(1,1))", 10.0, 0.0692918));   //(6)(confirmed)
  v.push_back(boost::make_tuple("(1,2,(1,1))"  , 10.0, 0.0092223));   //(7)(confirmed)
  v.push_back(boost::make_tuple("(2,1,(1,1))"  , 10.0, 0.0092223));
  v.push_back(boost::make_tuple("(1,(1,1),2)"  , 10.0, 0.0092223));
  v.push_back(boost::make_tuple("(2,(1,1),1)"  , 10.0, 0.0092223));
  v.push_back(boost::make_tuple("((1,1),1,2)"  , 10.0, 0.0092223));
  v.push_back(boost::make_tuple("((1,1),2,1)"  , 10.0, 0.0092223));
  v.push_back(boost::make_tuple("(1,(1,1,2))"  , 10.0, 0.0069190));   //(9)(confirmed)
  v.push_back(boost::make_tuple("(1,(1,2,1))"  , 10.0, 0.0069190));
  v.push_back(boost::make_tuple("(1,(2,1,1))"  , 10.0, 0.0069190));
  v.push_back(boost::make_tuple("((1,1,2),1)"  , 10.0, 0.0069190));
  v.push_back(boost::make_tuple("((1,2,1),1)"  , 10.0, 0.0069190));
  v.push_back(boost::make_tuple("((2,1,1),1)"  , 10.0, 0.0069190));
  //Quadrary
  v.push_back(boost::make_tuple("(1,(1,1,1,1))", 10.0, 0.0415140));   //(8)(confirmed)
  //Pentary
  v.push_back(boost::make_tuple("(1,1,1,1,1)"  , 10.0, 0.4162504));
  //Sum equals 6
  v.push_back(boost::make_tuple("(1,5)", 10.0, 0.0006660));
  v.push_back(boost::make_tuple("(2,4)", 10.0, 0.0001665));
  v.push_back(boost::make_tuple("(3,3)", 10.0, 0.0001110));
  v.push_back(boost::make_tuple("(1,(1,4))", 10.0, 0.0005804));
  v.push_back(boost::make_tuple("(1,(2,3))", 10.0, 0.0001679));
  v.push_back(boost::make_tuple("(1,(3,2))", 10.0, 0.0001679));
  v.push_back(boost::make_tuple("(1,(4,1))", 10.0, 0.0005804));
  v.push_back(boost::make_tuple("(2,(1,3))", 10.0, 0.0001991));
  v.push_back(boost::make_tuple("(2,(2,2))", 10.0, 0.0000925));
  v.push_back(boost::make_tuple("(2,(3,1))", 10.0, 0.0001991));
  v.push_back(boost::make_tuple("(3,(1,2))", 10.0, 0.0001880));
  v.push_back(boost::make_tuple("(3,(2,1))", 10.0, 0.0001880));
  v.push_back(boost::make_tuple("(4,(1,1))", 10.0, 0.0005043));
  //Trinary
  v.push_back(boost::make_tuple("(1,1,(1,3))"  , 10.0, 0.0012712));
  v.push_back(boost::make_tuple("(1,1,(2,2))"  , 10.0, 0.0005563));
  v.push_back(boost::make_tuple("(1,1,(3,1))"  , 10.0, 0.0012712));
  v.push_back(boost::make_tuple("(1,(1,3),1)"  , 10.0, 0.0012712));
  v.push_back(boost::make_tuple("(1,(2,2),1)"  , 10.0, 0.0005563));
  v.push_back(boost::make_tuple("(1,(3,1),1)"  , 10.0, 0.0012712));
  v.push_back(boost::make_tuple("((1,3),1,1)"  , 10.0, 0.0012712));
  v.push_back(boost::make_tuple("((2,2),1,1)"  , 10.0, 0.0005563));
  v.push_back(boost::make_tuple("((3,1),1,1)"  , 10.0, 0.0012712));
  v.push_back(boost::make_tuple("(1,2,(1,2))"  , 10.0, 0.0006346));
  v.push_back(boost::make_tuple("(2,1,(1,2))"  , 10.0, 0.0006346));
  v.push_back(boost::make_tuple("(1,2,(2,1))"  , 10.0, 0.0006346));
  v.push_back(boost::make_tuple("(2,1,(2,1))"  , 10.0, 0.0006346));
  v.push_back(boost::make_tuple("(1,(2,1),2)"  , 10.0, 0.0006346));
  v.push_back(boost::make_tuple("(1,(1,2),2)"  , 10.0, 0.0006346));
  v.push_back(boost::make_tuple("(2,(2,1),1)"  , 10.0, 0.0006346));
  v.push_back(boost::make_tuple("(2,(1,2),1)"  , 10.0, 0.0006346));
  v.push_back(boost::make_tuple("(1,3,(1,1))"  , 10.0, 0.0011913));
  v.push_back(boost::make_tuple("(1,(1,1),3)"  , 10.0, 0.0011913));
  v.push_back(boost::make_tuple("((1,1),1,3)"  , 10.0, 0.0011913));
  v.push_back(boost::make_tuple("(3,(1,1),1)"  , 10.0, 0.0011913));
  v.push_back(boost::make_tuple("((1,1),3,1)"  , 10.0, 0.0011913));
  v.push_back(boost::make_tuple("(1,1,(1,1,2))", 10.0, 0.0023165));
  v.push_back(boost::make_tuple("(1,1,(1,2,1))", 10.0, 0.0023165));
  v.push_back(boost::make_tuple("(1,1,(2,1,1))", 10.0, 0.0023165));
  v.push_back(boost::make_tuple("(1,(1,1,2),1)", 10.0, 0.0023165));
  v.push_back(boost::make_tuple("(1,(1,2,1),1)", 10.0, 0.0023165));
  v.push_back(boost::make_tuple("(1,(2,1,1),1)", 10.0, 0.0023165));
  v.push_back(boost::make_tuple("((1,1,2),1,1)", 10.0, 0.0023165));
  v.push_back(boost::make_tuple("((1,2,1),1,1)", 10.0, 0.0023165));
  v.push_back(boost::make_tuple("((2,1,1),1,1)", 10.0, 0.0023165));
  v.push_back(boost::make_tuple("(1,2,(1,1,1))", 10.0, 0.0023323));
  v.push_back(boost::make_tuple("(2,1,(1,1,1))", 10.0, 0.0023323));
  v.push_back(boost::make_tuple("(1,(1,1,1),2)", 10.0, 0.0023323));
  v.push_back(boost::make_tuple("(2,(1,1,1),1)", 10.0, 0.0023323));
  v.push_back(boost::make_tuple("((1,1,1),1,2)", 10.0, 0.0023323));
  v.push_back(boost::make_tuple("((1,1,1),2,1)", 10.0, 0.0023323));
  //Quadrary
  v.push_back(boost::make_tuple("(1,(1,1,1,2))"  , 10.0, 0.0027574));
  v.push_back(boost::make_tuple("(1,(1,1,2,1))"  , 10.0, 0.0027574));
  v.push_back(boost::make_tuple("(1,(1,2,1,1))"  , 10.0, 0.0027574));
  v.push_back(boost::make_tuple("(1,(2,1,1,1))"  , 10.0, 0.0027574));
  v.push_back(boost::make_tuple("((1,1,1,2),1)"  , 10.0, 0.0027574));
  v.push_back(boost::make_tuple("((1,1,2,1),1)"  , 10.0, 0.0027574));
  v.push_back(boost::make_tuple("((1,2,1,1),1)"  , 10.0, 0.0027574));
  v.push_back(boost::make_tuple("((2,1,1,1),1)"  , 10.0, 0.0027574));
  v.push_back(boost::make_tuple("(2,(1,1,1,1))"  , 10.0, 0.0028154));
  v.push_back(boost::make_tuple("((1,1,1,1),2)"  , 10.0, 0.0028154));
  //Pentary
  v.push_back(boost::make_tuple("(1,(1,1,1,1,1))", 10.0, 0.0183824));   //(7)
  v.push_back(boost::make_tuple("((1,1,1,1,1),1)", 10.0, 0.0183824));
  //Hexary
  v.push_back(boost::make_tuple("(1,1,1,1,1,1)"  , 10.0, 0.2775003));
  return v;
}
//---------------------------------------------------------------------------
///GetRootBranches obtains the root branches from a non-unary Newick.
///Examples:
///(1,2)               -> { 1     , 2             }
///(1,2,3)             -> { 1     , 2     , 3     }
///((1,1),(2,2),(3,3)) -> { (1,1) , (2,2) , (3,3) }
///From http://www.richelbilderbeek.nl/CppGetRootBranchesBinary.htm
const std::vector<std::vector<int> >
  GetRootBranches(const std::vector<int>& n)
{
  //#define DEBUG_GETROOTBRANCHES
  #ifdef  DEBUG_GETROOTBRANCHES
  TRACE_FUNC();
  TRACE(Newick::NewickToString(n));
  #endif

  assert(Newick::IsNewick(n));
  assert(!Newick::IsUnaryNewick(n));

  const int size = boost::numeric_cast<int>(n.size());
  std::vector<std::vector<int> > v;

  if (Newick::IsSimple(n))
  {
    for (int i=1; i!=size-1; ++i) //Skip brackets
    {
      v.push_back(
        Newick::CreateVector(
          static_cast<int>(Newick::bracket_open),
          n[i],
          static_cast<int>(Newick::bracket_close)));
    }
    assert(Newick::IsNewick(v.back()));
    assert(v.size() > 1);
    return v;
  }
  //Complex newick
  assert(!Newick::IsSimple(n));
  const std::vector<int> depth = Newick::GetDepth(n);

  assert(depth.size() == n.size());
  //Search for open and closing brackets in depth 1
  for (int i=0; i!=size; ++i)
  {
    if (depth[i] == 0 && n[i] > 0)
    {
      //C++0x initialization list
      std::vector<int> tmp;
      tmp.push_back(static_cast<int>(Newick::bracket_open));
      tmp.push_back(n[i]);
      tmp.push_back(static_cast<int>(Newick::bracket_close));
      v.push_back(tmp);
      assert(Newick::IsNewick(v.back()));
      continue;
    }
    if (depth[i] != 1 || n[i]!=Newick::bracket_open) continue;
    for (int j=i+1; j!=size; ++j)
    {
      if (depth[j] != 1 || n[j]!=Newick::bracket_close) continue;
      std::vector<int> w;
      w.push_back(Newick::bracket_open);
      std::copy(n.begin() + i + 1,n.begin() + j,std::back_inserter(w));
      w.push_back(Newick::bracket_close);
      assert(Newick::IsNewick(w));
      v.push_back(w);
      //Set from index i after current end
      i = j;
      break;
    }
  }
  assert(v.size() > 1);
  return v;
}
//---------------------------------------------------------------------------
///GetSimplerNewicksFrequencyPairs creates simpler, derived Newicks from a Newick.
///Its simpler Newicks are identical to those created by GetSimplerNewicks.
///From http://www.richelbilderbeek.nl/CppGetSimplerNewicksFrequencyPairs.htm
const std::vector<std::pair<std::vector<int>,int> >
  GetSimplerNewicksFrequencyPairs(const std::vector<int>& n)
{
  //#define DEBUG_GETSIMPLERNEWICKSFREQUENCYPAIRS
  #ifdef DEBUG_GETSIMPLERNEWICKSFREQUENCYPAIRS
  TRACE_FUNC();
  #endif

  assert(Newick::IsNewick(n));

  std::vector<std::pair<std::vector<int>,int> > newicks;
  const std::vector<int> depths = Newick::GetDepth(n);


  const int size = boost::numeric_cast<int>(n.size());
  for (int i = 0; i!=size; ++i)
  {
    assert(i >= 0);
    assert(i < size);
    if (n[i] < 1) continue;
    assert(n[i] > 0);
    if (n[i] > 1)
    {
      std::vector<int> new_newick(n);
      --new_newick[i];
      #ifdef DEBUG_GETSIMPLERNEWICKS
      {
        const std::string stored = Newick::NewickToString(new_newick);
        TRACE(stored);
      }
      #endif

      newicks.push_back( std::make_pair(new_newick,n[i]) );
      continue;
    }
    assert(n[i] == 1); //Most difficult...
    const int depth = depths[i];
    //j must first decrement, later increment with the same code
    int j_end  = -1;
    int j_step = -1;
    for (int j=i-1; ; j+=j_step)
    {
      //j must first decrement, later increment with the same code
      if (j == j_end
        //|| depths[j] < depth
        || (depths[j] == depth && n[j] < 0))
      {
        if (j_step == -1)
        {
          j = i + 1;
          j_end = size;
          j_step = 1;
        }
        else
        {
          break;
        }
      }
      assert(i!=j);
      assert(j >= 0);
      assert(j < size);
      //Only take frequencies of the same depth into account
      if (n[j] < 1 || depths[j] != depth) continue;
      std::vector<int> new_newick_with_zero(n);
      --new_newick_with_zero[i];
      assert(new_newick_with_zero[i] == 0);
      ++new_newick_with_zero[j];
      //Remove brackets after possibly lonely value
      #ifdef DEBUG_GETSIMPLERNEWICKS
      {
        const std::string newick_str_with_zeroes = Newick::DumbNewickToString(new_newick_with_zero);
        TRACE(newick_str_with_zeroes);
        const std::string dist_i_j = boost::lexical_cast<std::string>(std::abs(i - j));
        TRACE(dist_i_j)
      }
      #endif
      //If there is only one or two values between
      //the brackets, and one of these values was a
      //1 becoming added to the other, nullify the
      //1 and both brackets:
      //'((1,1),2)' -> '(00102)' -> '(1,2)'
      if (std::abs(i - j) == 1)
        //|| (std::abs(i - j) == 2 && n[i] == 1))
      {
        const int index_bracket_open  = std::min(i,j) - 1;
        const int index_bracket_close = std::max(i,j) + 1;
        if ( new_newick_with_zero[index_bracket_open]  == Newick::bracket_open
          && new_newick_with_zero[index_bracket_close] == Newick::bracket_close)
        {
          new_newick_with_zero[index_bracket_open]  = 0;
          new_newick_with_zero[index_bracket_close] = 0;
        }
      }
      #ifdef DEBUG_GETSIMPLERNEWICKS
      {
        const std::string newick_str_with_more_zeroes = Newick::DumbNewickToString(new_newick_with_zero);
        TRACE(newick_str_with_more_zeroes);
      }
      #endif
      //Remove decremented i and possibly nullified brackets
      std::vector<int> new_newick;
      std::remove_copy(
        new_newick_with_zero.begin(),
        new_newick_with_zero.end(),
        std::back_inserter(new_newick),
        0);
      //Add brackets if these are removed
      if (new_newick.front() != Newick::bracket_open
        || new_newick.back() != Newick::bracket_close)
      {
        new_newick = Newick::Surround(new_newick);
      }
      #ifdef DEBUG_GETSIMPLERNEWICKS
      {
        const std::string newick_str_done = Newick::DumbNewickToString(new_newick);
        TRACE(newick_str_done);
      }
      #endif
      #define DEBUG_2436964926435498753298216832187
      #ifdef  DEBUG_2436964926435498753298216832187
      if (!Newick::IsNewick(new_newick))
      {
        TRACE(Newick::DumbNewickToString(new_newick));
      }
      #endif

      assert(Newick::IsNewick(new_newick));
      newicks.push_back(std::make_pair(new_newick, 1));
      continue;
    }
  }

  return newicks;

  /*
  const int size = boost::numeric_cast<int>(n.size());
  for (int i = 0; i!=size; ++i)
  {
    if (n[i] < 1) continue;
    assert(n[i] > 0);
    if (n[i] > 1)
    {
      std::vector<int> new_newick(n);
      --new_newick[i];
      #ifdef DEBUG_GETSIMPLERNEWICKSFREQUENCYPAIRS
      std::clog << "Store: " << Newick::NewickToString(new_newick) << '\n';
      #endif
      newicks.push_back(std::make_pair(new_newick,n[i]));
      continue;
    }
    assert(n[i] == 1); //Most difficult...
    const int depth = depths[i];
    const int j_start = FindPosBefore(n,Newick::bracket_open,i);
    const int j_end   = FindPosAfter( n,Newick::bracket_close,i);
    assert(j_start >= 0);
    assert(j_end >= 0);
    assert(j_start <= boost::numeric_cast<int>(n.size()));
    assert(j_end <= boost::numeric_cast<int>(n.size()));
    for (int j=j_start; j!=j_end; ++j)
    {
      if (i==j) continue;
      if (n[j] < 1) continue;
      if (depths[j] != depth) continue;
      //Decrement index i to zero
      //Increment index j
      std::vector<int> new_newick_with_zero(n);
      --new_newick_with_zero[i];
      assert(new_newick_with_zero[i] == 0);
      ++new_newick_with_zero[j];
      //Remove brackets after possibly lonely value
      #ifdef DEBUG_GETSIMPLERNEWICKSFREQUENCYPAIRS
      std::clog << "1: "
        << Newick::DumbNewickToString(new_newick_with_zero)
        << " -> ["
        << i
        << "]="
        << new_newick_with_zero[i]
        << " - ["
        << j
        << "]="
        << new_newick_with_zero[j]
        << " = "
        << std::abs(i-j)
        << '\n';
      #endif
      if (std::abs(i - j) == 1)
      {
        const int index_bracket_open  = std::min(i,j) - 1;
        const int index_bracket_close = std::max(i,j) + 1;
        #ifdef DEBUG_GETSIMPLERNEWICKSFREQUENCYPAIRS
        std::clog
          << "["
          << index_bracket_open
          << "]-["
          << index_bracket_close
          << "]\n";
        #endif
        if ( new_newick_with_zero[index_bracket_open]  == Newick::bracket_open
          && new_newick_with_zero[index_bracket_close] == Newick::bracket_close)
        {
          new_newick_with_zero[index_bracket_open]  = 0;
          new_newick_with_zero[index_bracket_close] = 0;
          #ifdef DEBUG_GETSIMPLERNEWICKSFREQUENCYPAIRS
          std::clog << "2.5: " << Newick::DumbNewickToString(new_newick_with_zero) << '\n';
          #endif
        }
      }
      #ifdef DEBUG_GETSIMPLERNEWICKS
      std::clog << "2: " << Newick::DumbNewickToString(new_newick_with_zero) << '\n';
      #endif
      //Remove decremented i and possibly nullified brackets
      std::vector<int> new_newick;
      std::remove_copy(
        new_newick_with_zero.begin(),
        new_newick_with_zero.end(),
        std::back_inserter(new_newick),
        0);
      //Add brackets if these are removed
      if (new_newick.front() != Newick::bracket_open
        || new_newick.back() != Newick::bracket_close)
      {
        new_newick = Surround(new_newick);
      }
      #ifdef DEBUG_GETSIMPLERNEWICKSFREQUENCYPAIRS
      std::clog << "Store: " << Newick::DumbNewickToString(new_newick) << '\n';
      #endif
      assert(IsNewick(new_newick));
      assert(n[i] == 1);
      newicks.push_back(std::make_pair(new_newick,1));
      //newicks.push_back(new_newick);
      continue;
    }
  }
  */
}
//---------------------------------------------------------------------------

///Test tests all Newick functions
void Test()
{
  TRACE("Testing basic Newick functionality");
  //Check difference between C++98 and C++0x
  assert(Newick::CreateValidTrinaryNewicks() == NewickCpp98::CreateValidTrinaryNewicks());
  assert(Newick::GetKnownProbabilities() == NewickCpp98::GetKnownProbabilities());

  //Check conversions from std::string to std::vector #1
  {
    const std::vector<int> v = Newick::StringToNewick("(11,(22,33))");
    assert(v.size() == 7);
    assert(v[0]==Newick::bracket_open);
    assert(v[1]==11);
    assert(v[2]==Newick::bracket_open);
    assert(v[3]==22);
    assert(v[4]==33);
    assert(v[5]==Newick::bracket_close);
    assert(v[6]==Newick::bracket_close);
  }
  //Check if well-formed Newicks are accepted
  {
    const std::vector<std::string> v = Newick::CreateValidNewicks();
    BOOST_FOREACH(const std::string& s,v)
    {
      #ifdef TRACE_REJECTED_NEWICKS
      const std::string debug = "I must be accepted: " + s;
      TRACE(debug);
      #endif
      assert(Newick::IsNewick(s));
      const std::vector<int> v = Newick::StringToNewick(s);
      assert(Newick::IsNewick(v));
    }
  }

  //Check if ill-formed Newicks are rejected
  {
    #ifndef NDEBUG
    const std::vector<std::string> v = Newick::CreateInvalidNewicks();
    BOOST_FOREACH(const std::string& s,v)
    {
      #ifdef TRACE_REJECTED_NEWICKS
      const std::string debug = "I must be rejected: " + s;
      TRACE(debug);
      #endif
      assert(!Newick::IsNewick(s));
      //Cannot test if std::vector<int> versions are rejected,
      //because Newick::StringToNewick assumes a valid Newick
      //const std::vector<int> v = Newick::StringToNewick(s);
      //assert(!Newick::IsNewick(v));
    }
    #endif
  }
  //Check conversions from std::string to std::vector #2
  {
    const std::vector<int> v = Newick::StringToNewick("((11,22),33)");
    assert(v.size() == 7);
    assert(v[0]==Newick::bracket_open);
    assert(v[1]==Newick::bracket_open);
    assert(v[2]==11);
    assert(v[3]==22);
    assert(v[4]==Newick::bracket_close);
    assert(v[5]==33);
    assert(v[6]==Newick::bracket_close);
  }
  assert(Newick::CalcNumOfSymmetriesBinary(Newick::StringToNewick("(1,(3,1))"))==0);
  assert(Newick::CalcNumOfSymmetriesBinary(Newick::StringToNewick("(3,(1,1))"))==1);
  assert(Newick::CalcNumOfSymmetriesBinary(Newick::StringToNewick("(1,((1,1),(1,1)))"))==3);
  assert(Newick::CalcNumOfSymmetriesBinary(Newick::StringToNewick("(1,((1,1),(2,2)))"))==2);
  assert(Newick::CalcNumOfSymmetriesBinary(Newick::StringToNewick("(1,(2,3))"))==0);
  assert(Newick::CalcNumOfSymmetriesBinary(Newick::StringToNewick("(99,99)"))==1);
  assert(Newick::CalcNumOfSymmetriesBinary(Newick::StringToNewick("(3,(2,2))"))==1);
  assert(Newick::CalcNumOfSymmetriesBinary(Newick::StringToNewick("(2,(2,2))"))==1);
  assert(Newick::CalcNumOfSymmetriesBinary(Newick::StringToNewick("((3,3),(2,2))"))==2);
  assert(Newick::CalcNumOfSymmetriesBinary(Newick::StringToNewick("((3,3),(3,3))"))==3);
  assert(Newick::CalcNumOfSymmetriesBinary(Newick::StringToNewick("((3,3),(3,4))"))==1);
  assert(Newick::CalcNumOfSymmetriesBinary(Newick::StringToNewick("(((3,3),(4,4)),5)"))==2);
  assert(Newick::CalcNumOfSymmetriesBinary(Newick::StringToNewick("(((3,3),(5,5)),5)"))==2);
  assert(Newick::CalcNumOfSymmetriesBinary(Newick::StringToNewick("(((5,5),(5,5)),5)"))==3);
  assert(Newick::CalcNumOfSymmetriesBinary(Newick::StringToNewick("(((5,5),(5,5)),(4,4))"))==4);
  assert(Newick::CalcNumOfSymmetriesBinary(Newick::StringToNewick("(((5,5),(4,4)),(4,4))"))==3);
  assert(Newick::CalcNumOfSymmetriesBinary(Newick::StringToNewick("(((4,4),(4,4)),(4,4))"))==4);
  assert(Newick::CalcNumOfCombinationsBinary(Newick::StringToNewick("(3,(1,1))"))==10);
  assert(Newick::CalcNumOfCombinationsBinary(Newick::StringToNewick("(1,(3,1))"))==20);
  assert(Newick::CalcNumOfCombinationsBinary(Newick::StringToNewick("(1,(1,(1,(1,1))))"))==60);
  assert(Newick::CalcNumOfCombinationsBinary(Newick::StringToNewick("(1,((1,1),(1,1)))"))==15);
  assert(bigIntegerToString(Newick::FactorialBigInt(1))=="1");
  assert(bigIntegerToString(Newick::FactorialBigInt(2))=="2");
  assert(bigIntegerToString(Newick::FactorialBigInt(3))=="6");
  assert(bigIntegerToString(Newick::FactorialBigInt(4))=="24");
  assert(bigIntegerToString(Newick::FactorialBigInt(5))=="120");
  assert(bigIntegerToString(Newick::FactorialBigInt(6))=="720");
  assert(Newick::GetLeafMaxArity(Newick::StringToNewick("(1)"))   == 1);
  assert(Newick::GetLeafMaxArity(Newick::StringToNewick("(12)"))  == 1);
  assert(Newick::GetLeafMaxArity(Newick::StringToNewick("(123)")) == 1);
  assert(Newick::GetLeafMaxArity(Newick::StringToNewick("(1,2)"))   == 2);
  assert(Newick::GetLeafMaxArity(Newick::StringToNewick("(12,2)"))  == 2);
  assert(Newick::GetLeafMaxArity(Newick::StringToNewick("(123,2)")) == 2);
  assert(Newick::GetLeafMaxArity(Newick::StringToNewick("(1,(1,2))"))   == 2);
  assert(Newick::GetLeafMaxArity(Newick::StringToNewick("(1,(12,2))"))  == 2);
  assert(Newick::GetLeafMaxArity(Newick::StringToNewick("(1,(123,2))")) == 2);
  assert(Newick::GetLeafMaxArity(Newick::StringToNewick("((1,2),3)"))   == 2);
  assert(Newick::GetLeafMaxArity(Newick::StringToNewick("((12,2),3)"))  == 2);
  assert(Newick::GetLeafMaxArity(Newick::StringToNewick("((123,2),3)")) == 2);
  assert(Newick::GetLeafMaxArity(Newick::StringToNewick("(1,2,3)"))   == 3);
  assert(Newick::GetLeafMaxArity(Newick::StringToNewick("(12,2,3)"))  == 3);
  assert(Newick::GetLeafMaxArity(Newick::StringToNewick("(123,2,3)")) == 3);
  assert(Newick::GetLeafMaxArity(Newick::StringToNewick("(1,(1,2,3))"))   == 3);
  assert(Newick::GetLeafMaxArity(Newick::StringToNewick("(1,(12,2,3))"))  == 3);
  assert(Newick::GetLeafMaxArity(Newick::StringToNewick("(1,(123,2,3))")) == 3);
  assert(Newick::GetLeafMaxArity(Newick::StringToNewick("((1,2,3),4)"))   == 3);
  assert(Newick::GetLeafMaxArity(Newick::StringToNewick("((12,2,3),4)"))  == 3);
  assert(Newick::GetLeafMaxArity(Newick::StringToNewick("((123,2,3),4)")) == 3);

  assert(fuzzy_equal_to()(  2.0,Newick::CalcDenominator(Newick::StringToNewick("(1,1)"),10.0)));
  assert(fuzzy_equal_to()(  6.0,Newick::CalcDenominator(Newick::StringToNewick("((1,1),1)"),10.0)));
  assert(fuzzy_equal_to()( 26.0,Newick::CalcDenominator(Newick::StringToNewick("(1,2)"),10.0)));
  assert(fuzzy_equal_to()( 32.0,Newick::CalcDenominator(Newick::StringToNewick("((1,1),2)"),10.0)));
  assert(fuzzy_equal_to()( 32.0,Newick::CalcDenominator(Newick::StringToNewick("(2,(1,1))"),10.0)));
  assert(fuzzy_equal_to()( 50.0,Newick::CalcDenominator(Newick::StringToNewick("((1,1),3)"),10.0)));
  assert(fuzzy_equal_to()( 80.0,Newick::CalcDenominator(Newick::StringToNewick("((1,2),3)"),10.0)));
  assert(fuzzy_equal_to()( 80.0,Newick::CalcDenominator(Newick::StringToNewick("((3,1),2)"),10.0)));
  assert(fuzzy_equal_to()( 80.0,Newick::CalcDenominator(Newick::StringToNewick("((2,3),1)"),10.0)));
  assert(fuzzy_equal_to()(102.0,Newick::CalcDenominator(Newick::StringToNewick("((2,1),4)"),10.0)));
  assert(fuzzy_equal_to()(152.0,Newick::CalcDenominator(Newick::StringToNewick("(2,(1,(3,3)))"),10.0)));
  assert(fuzzy_equal_to()(162.0,Newick::CalcDenominator(Newick::StringToNewick("((2,3),4)"),10.0)));
  assert(fuzzy_equal_to()(180.0,Newick::CalcDenominator(Newick::StringToNewick("((1,2),(3,4))"),10.0)));
  assert(fuzzy_equal_to()(180.0,Newick::CalcDenominator(Newick::StringToNewick("((4,1),(2,3))"),10.0)));
  assert(fuzzy_equal_to()(180.0,Newick::CalcDenominator(Newick::StringToNewick("((3,4),(1,2))"),10.0)));
  assert(fuzzy_equal_to()(180.0,Newick::CalcDenominator(Newick::StringToNewick("((2,3),(4,1))"),10.0)));
  {
    /*
    const std::vector<int> v = { 0,1,2,3,4,5,6 };
    assert(Newick::FindPosAfter(v,3,2)==3);
    assert(Newick::FindPosAfter(v,4,2)==4);
    assert(Newick::FindPosAfter(v,5,2)==5);
    assert(Newick::FindPosAfter(v,6,2)==6);
    assert(Newick::FindPosBefore(v,3,4)==3);
    assert(Newick::FindPosBefore(v,2,4)==2);
    assert(Newick::FindPosBefore(v,1,4)==1);
    assert(Newick::FindPosBefore(v,0,4)==0);
    */
  }
  /* C++98
  {
    const std::vector<int> v = Newick::GetDepth(Newick::StringToNewick("(1,(2,2))"));
    const std::vector<int> w = Newick::GetDepth(Newick::StringToNewick("(9,(9,9))"));
    const std::vector<int> x = { 0,0,1,1,1,1,0 };
    assert(v == x);
    assert(w == x);
  }
  {
    const std::vector<int> v = Newick::GetDepth(Newick::StringToNewick("((2,2),1)"));
    const std::vector<int> w = { 0,1,1,1,1,0,0 };
    assert(v == w);
  }
  {
    const std::vector<int> v = Newick::GetDepth(Newick::StringToNewick("(1,(2,2),1)"));
    const std::vector<int> w = { 0,0,1,1,1,1,0,0 };
    assert(v == w);
  }
  {
    const std::vector<int> v = Newick::GetDepth(Newick::StringToNewick("(1,(2,3),4,(5,6))"));
    const std::vector<int> w = { 0,0,1,1,1,1,0,1,1,1,1,0 };
    assert(v == w);
  }
  {
    const std::vector<int> v = Newick::GetDepth(Newick::StringToNewick("(1,(2,3),(5,6))"));
    const std::vector<int> w = { 0,0,1,1,1,1,1,1,1,1,0 };
    assert(v == w);
  }
  {
    const std::vector<int> v = Newick::GetDepth(Newick::StringToNewick("(1,(2,(3,4)),((5,6),7))"));
    const std::vector<int> w = { 0,0,1,1,2,2,2,2,1,1,2,2,2,2,1,1,0 };
    assert(v == w);
  }
  */
  //Test GetRootBranches
  {
    const std::vector<std::vector<int> > v = Newick::GetRootBranches(Newick::StringToNewick("(1,2)"));
    assert(v.size() == 2);
    assert(std::find(v.begin(),v.end(),
      Newick::StringToNewick("(1)")) != v.end());
    assert(std::find(v.begin(),v.end(),
      Newick::StringToNewick("(2)")) != v.end());
  }
  {
    const std::vector<std::vector<int> > v = Newick::GetRootBranches(Newick::StringToNewick("(1,(2,3))"));
    assert(v.size() == 2);
    assert(std::find(v.begin(),v.end(),
      Newick::StringToNewick("(1)")) != v.end());
    assert(std::find(v.begin(),v.end(),
      Newick::StringToNewick("(2,3)")) != v.end());
  }
  {
    const std::vector<std::vector<int> > v = Newick::GetRootBranches(Newick::StringToNewick("(1,2,(3,4))"));
    assert(v.size() == 3);
    assert(std::find(v.begin(),v.end(),
      Newick::StringToNewick("(1)")) != v.end());
    assert(std::find(v.begin(),v.end(),
      Newick::StringToNewick("(2)")) != v.end());
    assert(std::find(v.begin(),v.end(),
      Newick::StringToNewick("(3,4)")) != v.end());
  }
  //Compare C++98 and C++0x version
  {
    const std::vector<std::string> v = Newick::CreateValidBinaryNewicks();
    BOOST_FOREACH(const std::string& s,v)
    {
      const std::vector<int> n = Newick::StringToNewick(s);
      assert(Newick::GetRootBranches(n) == NewickCpp98::GetRootBranches(n));
    }
  }

  //Check if binary and trinary Newicks are detected correctly
  {
    const std::vector<std::string> v = Newick::CreateValidBinaryNewicks();
    BOOST_FOREACH(const std::string& s,v)
    {
      const std::vector<int> n = Newick::StringToNewick(s);
      assert(Newick::IsBinaryNewick(n));
    }
  }
  //Check if unary Newicks are detected correctly
  {
    const std::vector<std::string> v = Newick::CreateValidUnaryNewicks();
    BOOST_FOREACH(const std::string& s,v)
    {
      const std::vector<int> n = Newick::StringToNewick(s);
      assert( Newick::GetLeafMaxArity(n)<=1);
      assert( Newick::IsUnaryNewick(n));
      assert(!Newick::IsBinaryNewick(n));
      assert(!Newick::IsTrinaryNewick(n));
    }
  }
  //Check if binary Newicks are detected correctly
  {
    const std::vector<std::string> v = Newick::CreateValidBinaryNewicks();
    BOOST_FOREACH(const std::string& s,v)
    {
      const std::vector<int> n = Newick::StringToNewick(s);
      assert( Newick::GetLeafMaxArity(n)<=2);
      assert(!Newick::IsUnaryNewick(n));
      assert( Newick::IsBinaryNewick(n));
      assert(!Newick::IsTrinaryNewick(n));
    }
  }
  //Check if trinary Newicks are detected correctly
  {
    const std::vector<std::string> v = Newick::CreateValidTrinaryNewicks();
    BOOST_FOREACH(const std::string& s,v)
    {
      //TRACE(s);
      const std::vector<int> n = Newick::StringToNewick(s);
      assert( Newick::GetLeafMaxArity(n)<=3);
      assert(!Newick::IsUnaryNewick(n));
      assert(!Newick::IsBinaryNewick(n));
      assert( Newick::IsTrinaryNewick(n));
    }
  }
  //Test binary Newick
  {
    const std::string s("(1,(2,3))");
    const std::vector<std::vector<int> > n = Newick::GetSimplerNewicks(Newick::StringToNewick(s));
    //#define DEBUG_1_BO_1_2_3_BC
    #ifdef  DEBUG_1_BO_1_2_3_BC
    BOOST_FOREACH(const auto& t,n)
    {
      TRACE(Newick::NewickToString(t));
    }
    #endif
    assert(n.size() == 2);
    assert(std::find(n.begin(),n.end(),Newick::StringToNewick("(1,(1,3))"))
      != n.end());
    assert(std::find(n.begin(),n.end(),Newick::StringToNewick("(1,(2,2))"))
      != n.end());
  }
  {
    const std::string s("(1,(2,3,4))");
    const std::vector<std::vector<int> > n = Newick::GetSimplerNewicks(Newick::StringToNewick(s));
    assert(n.size() == 3);
    assert(std::find(n.begin(),n.end(),Newick::StringToNewick("(1,(1,3,4))"))
      != n.end());
    assert(std::find(n.begin(),n.end(),Newick::StringToNewick("(1,(2,2,4))"))
      != n.end());
    assert(std::find(n.begin(),n.end(),Newick::StringToNewick("(1,(2,3,3))"))
      != n.end());
  }
  {
    const std::string s("(1,(1,3,4))");
    const std::vector<std::vector<int> > n = Newick::GetSimplerNewicks(Newick::StringToNewick(s));
    //#define DEBUG_1_BO_1_3_4_BC
    #ifdef  DEBUG_1_BO_1_3_4_BC
    TRACE(boost::lexical_cast<std::string>(n.size()));
    BOOST_FOREACH(const auto& t,n)
    {
      TRACE(Newick::NewickToString(t));
    }
    #endif
    assert(n.size() == 4);
    assert(std::find(n.begin(),n.end(),Newick::StringToNewick("(1,(4,4))"))
      != n.end());
    assert(std::find(n.begin(),n.end(),Newick::StringToNewick("(1,(3,5))"))
      != n.end());
    assert(std::find(n.begin(),n.end(),Newick::StringToNewick("(1,(1,2,4))"))
      != n.end());
    assert(std::find(n.begin(),n.end(),Newick::StringToNewick("(1,(1,3,3))"))
      != n.end());
  }
  {
    const std::string s("(1,(1,3,4))");
    const std::vector<std::pair<std::vector<int>,int> > n
      = GetSimplerNewicksFrequencyPairs(Newick::StringToNewick(s));
    typedef std::pair<std::vector<int>,int> Pair;
    #ifdef TRACE_GETSIMPLERNEWICKSFREQUENCYPAIRS_1_134
    BOOST_FOREACH(const Pair& p, n)
    {
      std::cout << Newick::NewickToString(p.first) << '\n';
    }
    #endif
    assert(n.size() == 4);
    assert(std::find(n.begin(),n.end(),
      std::make_pair(Newick::StringToNewick("(1,(4,4))"),1))
      != n.end());
    assert(std::find(n.begin(),n.end(),
      std::make_pair(Newick::StringToNewick("(1,(3,5))"),1))
      != n.end());
    assert(std::find(n.begin(),n.end(),
      std::make_pair(Newick::StringToNewick("(1,(1,2,4))"),3))
      != n.end());
    assert(std::find(n.begin(),n.end(),
      std::make_pair(Newick::StringToNewick("(1,(1,3,3))"),4))
      != n.end());
  }
  {
    const std::string s("((1,1),2)");
    const std::vector<std::vector<int> > n = Newick::GetSimplerNewicks(
      Newick::StringToNewick(s));
    //#define DEBUG_BO_1_1_BC_2
    #ifdef  DEBUG_BO_1_1_BC_2
    BOOST_FOREACH(const auto& t,n)
    {
      TRACE(Newick::NewickToString(t));
    }
    #endif
    assert(n.size() == 3);
    assert(std::find(n.begin(),n.end(),
      Newick::StringToNewick("(2,2)"))
      != n.end());
    assert(std::find(n.begin(),n.end(),
      Newick::StringToNewick("((1,1),1)"))
      != n.end());
  }
  {
    const std::string s("((1,1),2)");
    typedef std::pair<std::vector<int>,int> Pair;
    const std::vector<Pair> n
      = GetSimplerNewicksFrequencyPairs(Newick::StringToNewick(s));
    #ifdef TRACE_GETSIMPLERNEWICKSFREQUENCYPAIRS_11_2
    BOOST_FOREACH(const Pair& p, n)
    {
      std::clog << Newick::NewickToString(p.first) << '\n';
    }
    #endif
    assert(n.size() == 3);
    assert(std::find(n.begin(),n.end(),
      std::make_pair(Newick::StringToNewick("(2,2)"),1))
      != n.end());
    assert(std::find(n.begin(),n.end(),
      std::make_pair(Newick::StringToNewick("((1,1),1)"),2))
      != n.end());
  }
  {
    const std::string s("((2,1),4)");
    const std::vector<std::vector<int> > n = Newick::GetSimplerNewicks(
      Newick::StringToNewick(s));
    assert(n.size() == 3);
    assert(std::find(n.begin(),n.end(),
      Newick::StringToNewick("(3,4)"))
      != n.end());
    assert(std::find(n.begin(),n.end(),
      Newick::StringToNewick("((1,1),4)"))
      != n.end());
    assert(std::find(n.begin(),n.end(),
      Newick::StringToNewick("((2,1),3)"))
      != n.end());
  }
  {
    const std::string s("((2,1),4)");
    typedef std::pair<std::vector<int>,int> Pair;
    const std::vector<Pair> n
      = GetSimplerNewicksFrequencyPairs(Newick::StringToNewick(s));
    #ifdef TRACE_GETSIMPLERNEWICKSFREQUENCYPAIRS_21_2
    BOOST_FOREACH(const Pair& p, n)
    {
      TRACE(Newick::NewickToString(p.first));
    }
    #endif
    assert(n.size() == 3);
    assert(std::find(n.begin(),n.end(),
      std::make_pair(Newick::StringToNewick("(3,4)"),1))
      != n.end());
    assert(std::find(n.begin(),n.end(),
      std::make_pair(Newick::StringToNewick("((1,1),4)"),2))
      != n.end());
    assert(std::find(n.begin(),n.end(),
      std::make_pair(Newick::StringToNewick("((2,1),3)"),4))
      != n.end());
  }
  {
    const std::string s("((2,3),4)");
    const std::vector<std::vector<int> > n = Newick::GetSimplerNewicks(
      Newick::StringToNewick(s));
    assert(n.size() == 3);
    assert(std::find(n.begin(),n.end(),
      Newick::StringToNewick("((1,3),4)"))
      != n.end());
    assert(std::find(n.begin(),n.end(),
      Newick::StringToNewick("((2,2),4)"))
      != n.end());
    assert(std::find(n.begin(),n.end(),
      Newick::StringToNewick("((2,3),3)"))
      != n.end());
  }
  {
    const std::string s("((2,3),4)");
    typedef std::pair<std::vector<int>,int> Pair;
    const std::vector<Pair> n
      = GetSimplerNewicksFrequencyPairs(Newick::StringToNewick(s));
    #ifdef TRACE_GETSIMPLERNEWICKSFREQUENCYPAIRS_23_4
    BOOST_FOREACH(const Pair& p, n)
    {
      std::cout << Newick::NewickToString(p.first) << '\n';
    }
    #endif
    assert(n.size() == 3);
    assert(std::find(n.begin(),n.end(),
      std::make_pair(Newick::StringToNewick("((1,3),4)"),2))
      != n.end());
    assert(std::find(n.begin(),n.end(),
      std::make_pair(Newick::StringToNewick("((2,2),4)"),3))
      != n.end());
    assert(std::find(n.begin(),n.end(),
      std::make_pair(Newick::StringToNewick("((2,3),3)"),4))
      != n.end());
  }
  //Compare GetSimplerNewicks and
  //GetSimplerNewicksFrequencyPairs
  {
    const std::vector<std::string> newicks
      = Newick::CreateValidNewicks();
    BOOST_FOREACH(const std::string& newick_str,newicks)
    {
      const std::vector<int> newick
        = Newick::StringToNewick(newick_str);
      const std::vector<std::vector<int> > v1
        = Newick::GetSimplerNewicks(newick);
      const std::vector<std::pair<std::vector<int>,int> > v2
        = Newick::GetSimplerNewicksFrequencyPairs(newick);
      assert(v1.size() == v2.size());
      const int size = boost::numeric_cast<int>(v1.size());
      for (int i=0; i!=size; ++i)
      {
        #define DEBUG_COMPARE_GSN_VS_GSNFP
        #ifdef  DEBUG_COMPARE_GSN_VS_GSNFP
        if (v1[i] != v2[i].first)
        {
          TRACE("ERROR: DIFFERENT NEWICK SIMPLIFICATIONS");
          TRACE(