Go back to Richel Bilderbeek's homepage.

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

 

 

 

 

 

(C++) TwoDigitNewick

 

STLQt CreatorLubuntu

 

TwoDigitNewick is a Newick class.

Technical facts

 

 

 

 

 

 

./CppTwoDigitNewick/CppTwoDigitNewick.pri

 

INCLUDEPATH += \
    ../../Classes/CppTwoDigitNewick

SOURCES += \
    ../../Classes/CppTwoDigitNewick/twodigitnewick.cpp \
    ../../Classes/CppTwoDigitNewick/twodigitnewickderivative.cpp \
    ../../Classes/CppTwoDigitNewick/twodigitnewickindexer.cpp \
    ../../Classes/CppTwoDigitNewick/twodigitnewickindextable.cpp \
    ../../Classes/CppTwoDigitNewick/twodigitnewicks.cpp

HEADERS  += \
    ../../Classes/CppTwoDigitNewick/twodigitnewick.h \
    ../../Classes/CppTwoDigitNewick/twodigitnewickderivative.h \
    ../../Classes/CppTwoDigitNewick/twodigitnewickindexer.h \
    ../../Classes/CppTwoDigitNewick/twodigitnewickindextable.h \
    ../../Classes/CppTwoDigitNewick/twodigitnewicks.h

OTHER_FILES += \
    ../../Classes/CppTwoDigitNewick/Licence.txt

 

 

 

 

 

./CppTwoDigitNewick/twodigitnewick.h

 

//---------------------------------------------------------------------------
/*
TestTwoDigitNewick, tool to test the two-digit-Newick architecture
Copyright (C) 2010-2015 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/ToolTestTwoDigitNewick.htm
//---------------------------------------------------------------------------
#ifndef TWODIGITNEWICK_H
#define TWODIGITNEWICK_H

#include <string>
#include <vector>

#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Weffc++"
#include "twodigitnewickderivative.h"
#pragma GCC diagnostic pop

namespace ribi {

///TwoDigitNewick contains all
///TwoDigitNewickDerivative that can be
///constructed from a phylogeny. For example,
///if from a certain phylogeny three derived
///phylogenies can be constructed, TwoDigitNewickDerivativesData
///will hold three elements
struct TwoDigitNewick
{
  //An empty TwoDigitNewick
  TwoDigitNewick();

  ///A TwoDogitNewick cannot be created without
  ///its derivatives: a TwoDigitNewick IS its
  ///derivatives in a way.
  ///sum_above_zero and sum_above_one are needed to
  ///calculate its denominator
  TwoDigitNewick(
    const std::vector<TwoDigitNewickDerivative>& derivatives,
    const int sum_above_zero,
    const int sum_above_one);

  ///Empty returns !IsComplete
  bool Empty() const;

  const std::vector<TwoDigitNewickDerivative>& GetDerivatives() const;
  ///IsComplete determines if the TwoDigitNewick is
  ///initialized completely
  bool IsComplete() const;
  bool IsProbabilityKnown() const;
  double GetDenominator() const;
  double GetProbability() const;
  int GetSumTermsAboveOne() const;
  int GetSumTermsAboveZero() const;
  void SetProbability(const double p);
  static void SetTheta(const double theta);

  static std::string GetVersion() noexcept;
  static std::vector<std::string> GetVersionHistory() noexcept;

  private:
  ///m_derivatives contains all the information
  ///needed to get to this Newick's derivatives
  std::vector<TwoDigitNewickDerivative> m_derivatives;

  ///m_probability denotes the probability
  ///a Newick exists.
  ///A negative value of m_probability denotes
  ///that it is not yet calculated.
  double m_probability;

  ///m_denominator constant for a Newick.
  double m_denominator;

  int m_sum_terms_above_zero;
  int m_sum_terms_above_one;

  static double sm_theta;

  double CalculateDenominator(
    const int sum_above_zero,
    const int sum_above_one) const;

  public:

  static double CalculateProbability(
    const std::string& newick,
    const double theta);

  private:
  #ifndef NDEBUG
  static void Test() noexcept;
  #endif

};

} //~namespace ribi

#endif // TWODIGITNEWICK_H

 

 

 

 

 

./CppTwoDigitNewick/twodigitnewick.cpp

 

//---------------------------------------------------------------------------
/*
TestTwoDigitNewick, tool to test the two-digit-Newick architecture
Copyright (C) 2010-2015 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/ToolTestTwoDigitNewick.htm
//---------------------------------------------------------------------------
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Weffc++"
#include "twodigitnewick.h"

#include <cassert>
#include <iostream>


#include <boost/numeric/conversion/cast.hpp>

#include "binarynewickvector.h"
#include "newick.h"
#include "twodigitnewickindexer.h"

#pragma GCC diagnostic pop

double ribi::TwoDigitNewick::sm_theta = -1.0;

ribi::TwoDigitNewick::TwoDigitNewick()
  : m_derivatives{},
    m_probability(-1.0),
    m_denominator(-1.0),
    m_sum_terms_above_zero(-1),
    m_sum_terms_above_one(-1)
{
  assert(this->Empty());
}

ribi::TwoDigitNewick::TwoDigitNewick(
  const std::vector<TwoDigitNewickDerivative>& derivatives,
  const int sum_above_zero,
  const int sum_above_one)
  : m_derivatives(derivatives),
    m_probability(-1.0),
    m_denominator(
      CalculateDenominator(sum_above_zero,sum_above_one)),
    m_sum_terms_above_zero(sum_above_zero),
    m_sum_terms_above_one(sum_above_one)
{

}

double ribi::TwoDigitNewick::CalculateDenominator(
  const int sum_above_zero,
  const int sum_above_one) const
{
  assert(sm_theta >= 0.0);
  assert(sum_above_zero >= 0);
  #ifndef NDEBUG
  if (sum_above_one < 0)
  {
    std::cerr << "Invalid sum_above_one: " << sum_above_one << '\n';
  }
  #endif
  assert(sum_above_one >= 0);
  const double d
    = boost::numeric_cast<double>(
      sum_above_zero * (sum_above_zero - 1))
    + (boost::numeric_cast<double>(sum_above_one)
       * sm_theta);

  return d;
}

double ribi::TwoDigitNewick::CalculateProbability(
  const std::string& newick_str,
  const double theta)
{
  assert(Newick().IsUnaryNewick(Newick().StringToNewick(newick_str))
      || Newick().IsBinaryNewick(Newick().StringToNewick(newick_str)));
  ribi::TwoDigitNewick::SetTheta(theta);
  const BinaryNewickVector n(newick_str);
  const TwoDigitNewickIndexer i(n,theta);

  return i.GetProbability();
}

bool ribi::TwoDigitNewick::Empty() const
{
  return m_derivatives.empty();
}

double ribi::TwoDigitNewick::GetDenominator() const
{
  assert(IsComplete());
  return m_denominator;
}

const std::vector<ribi::TwoDigitNewickDerivative>& ribi::TwoDigitNewick::GetDerivatives() const
{
  return m_derivatives;
}

double ribi::TwoDigitNewick::GetProbability() const
{
  assert(IsProbabilityKnown());
  return m_probability;
}

int ribi::TwoDigitNewick::GetSumTermsAboveOne() const
{
  assert(m_sum_terms_above_one >= 0);
  return m_sum_terms_above_one;
}

int ribi::TwoDigitNewick::GetSumTermsAboveZero() const
{
  assert(m_sum_terms_above_zero >= 0);
  return m_sum_terms_above_zero;
}

std::string ribi::TwoDigitNewick::GetVersion() noexcept
{
  return "1.1";
}

std::vector<std::string> ribi::TwoDigitNewick::GetVersionHistory() noexcept
{
  std::vector<std::string> v;
  v.push_back("2010-08-22: version 1.0: initial version");
  v.push_back("2011-02-20: version 1.1: added version history");
  return v;
}

bool ribi::TwoDigitNewick::IsComplete() const
{
  return (!m_derivatives.empty()
    && m_sum_terms_above_zero >= 0
    && m_sum_terms_above_one  >= 0);
}

bool ribi::TwoDigitNewick::IsProbabilityKnown() const
{
  return m_probability >= 0.0;
}

void ribi::TwoDigitNewick::SetProbability(const double p)
{
  //if (p < 0.0 || p > 1.0)
  //{
  //  TRACE(p);
  //}
  assert(p >= 0.0);
  assert(p <= 1.00001);
  m_probability = p;
}

void ribi::TwoDigitNewick::SetTheta(const double theta)
{
  assert(theta >= 0.0);
  sm_theta = theta;
}

#ifndef NDEBUG
void ribi::TwoDigitNewick::Test() noexcept
{
  ribi::TwoDigitNewick::SetTheta(10.0);
  const std::vector<std::string> v = Newick().CreateValidNewicks();
  for(const std::string& s: v)
  {
    if ( Newick().CalcComplexity(Newick().StringToNewick(s))
      >  BigInteger(1000000) )
    {
      continue;
    }
    if (Newick().IsBinaryNewick(Newick().StringToNewick(s)))
    {
      BinaryNewickVector n(s);
      TwoDigitNewickIndexer(n,10.0);
    }
  }
}
#endif // NDEBUG

 

 

 

 

 

./CppTwoDigitNewick/twodigitnewickderivative.h

 

//---------------------------------------------------------------------------
/*
TestTwoDigitNewick, tool to test the two-digit-Newick architecture
Copyright (C) 2010-2015 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/ToolTestTwoDigitNewick.htm
//---------------------------------------------------------------------------
#ifndef TWODIGITDERIVATIVE_H
#define TWODIGITDERIVATIVE_H

namespace ribi {

///TwoDigitNewickDerivative contains
///the index of the derived phylogeny
///and the value that must be changed
///to get there. For example, if for a
///complex phylogeny (of unknown index)
///a '3' must be changed to a '2' to get
///to the phylogeny with index 42,
///TwoDigitNewickDerivatives has an
///m_derived_index of 42 and an
///m_value_changed of 3.
///For calculations it is important to know which value
///has changed, sometimes with another.
///If no other value changes, set m_other_value_changed to zero.
///Some examples of using m_value_changed and m_other_value_changed:\n
///(1,1) ->   (2), then m_value_changed == 1 && m_other_value_changed == 1\n
///(1,2) ->   (3), then m_value_changed == 1 && m_other_value_changed == 2\n
///(1,2) -> (1,1), then m_value_changed == 2 && m_other_value_changed == 0\n
///(1,3) ->   (4), then m_value_changed == 1 && m_other_value_changed == 3\n
///(1,3) -> (1,2), then m_value_changed == 3 && m_other_value_changed == 0\n
///(2,3) -> (1,3), then m_value_changed == 2 && m_other_value_changed == 0\n
///(2,3) -> (2,2), then m_value_changed == 3 && m_other_value_changed == 0\n
///A valid TwoDigitNewickDerivative has either:\n
///- m_value_changed == 1 && m_other_value_changed  > 0\n
///- m_value_changed  > 1 && m_other_value_changed == 0\n
///This is checked at the constructor
struct TwoDigitNewickDerivative
{
  TwoDigitNewickDerivative(
    const int derived_index,
    const int value_changed,
    const int other_value);

  int m_derived_index;
  int m_value_changed;
  int m_other_value_changed;
};

} //~namespace ribi

#endif // TWODIGITDERIVATIVE_H

 

 

 

 

 

./CppTwoDigitNewick/twodigitnewickderivative.cpp

 

//---------------------------------------------------------------------------
/*
TestTwoDigitNewick, tool to test the two-digit-Newick architecture
Copyright (C) 2010-2015 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/ToolTestTwoDigitNewick.htm
//---------------------------------------------------------------------------
#include "twodigitnewickderivative.h"

#include <cassert>

ribi::TwoDigitNewickDerivative::TwoDigitNewickDerivative(
  const int derived_index,
  const int value_changed,
  const int other_value_changed)
  : m_derived_index(derived_index),
    m_value_changed(value_changed),
    m_other_value_changed(other_value_changed)
{
  assert(m_derived_index > 0);
  assert(m_value_changed > 0);
  assert( (m_value_changed == 1 && m_other_value_changed  > 0)
       || (m_value_changed  > 1 && m_other_value_changed == 0) );
}

 

 

 

 

 

./CppTwoDigitNewick/twodigitnewickindexer.h

 

//---------------------------------------------------------------------------
/*
TestTwoDigitNewick, tool to test the two-digit-Newick architecture
Copyright (C) 2010-2015 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/ToolTestTwoDigitNewick.htm
//---------------------------------------------------------------------------
#ifndef TWODIGITNEWICKINDEXER_H
#define TWODIGITNEWICKINDEXER_H

#include <cassert>
#include <vector>

#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Weffc++"
#include "twodigitnewickindextable.h"
#include "twodigitnewicks.h"
#include "twodigitnewick.h"
#include "binarynewickvector.h"
#pragma GCC diagnostic pop

namespace ribi {

///NewickIndex converts any newick to (X,Y)
struct TwoDigitNewickIndexer
{
  //TwoDigitNewickIndexer constructor does all the work
  TwoDigitNewickIndexer(const BinaryNewickVector& n, const double theta);

  ///ConstructNewick constructs a full BinaryNewickVector from
  ///the TwoDigitNewick at index i.
  ///ConstructNewick is for debugging purposes only,
  ///as it's the idea to work with non-full determined
  ///(that is: two-digit) Newicks
  const BinaryNewickVector ConstructNewick(const int i) const;

  int GetCalculatedToIndex() const { return m_calculated_to_index; }
  ///GetData allows a peek at the x-y ordered data
  const TwoDigitNewickIndexTable& GetIndexTable() const
  {
    return m_index_table;
  }
  const TwoDigitNewick& GetNewick(const int i) const
  {
    return m_newicks.GetNewick(i);
  }
  const TwoDigitNewick& GetNewick(const int x, const int y) const
  {
    assert(m_index_table.CanGetIndex(x,y));
    const int i = m_index_table.GetIndex(x,y);
    return GetNewick(i);
  }
  const TwoDigitNewicks& GetNewicks() const
  {
    return m_newicks;
  }
  ///GetProbability returns the probability of the BinaryNewickVector
  ///given at the constructor
  double GetProbability() const;

  int GetReserved() const { return m_reserved; }


  ///After a leaf has been cut, then
  void TryToCalculateNewNewicks();
  void UpdateCalculatedFromIndex();


  private:

  int m_calculated_to_index;
  int m_current_index;
  TwoDigitNewickIndexTable m_index_table;
  TwoDigitNewicks m_newicks;

  ///m_probability is the probability of the given Newick.
  ///m_probability is calculated in the constructor of TwoDigitNewickIndexer
  ///and can be obtained by GetProbability.
  double m_probability;

  const int m_reserved;
  const double m_theta;

  ///Calculate the Newick probability of Newick (a,b).
  ///Both a and b are simple, that is: no index
  ///for a complex Newick
  //double CalculateEwensProbability(
  //  const int a,const int b) const;

  ///CalculateReserved calculates the index
  ///that must be reserver
  int CalculateReserved(const BinaryNewickVector& n) const;

  const TwoDigitNewick CreateTwoDigitNewickDerivatives(const int x, const int y,
    const int sum_above_zero,
    const int sum_above_one);

  const TwoDigitNewick CreateTwoDigitNewickDerivativesSimpleSimple(
    const int x,
    const int y);

  const TwoDigitNewick CreateTwoDigitNewickDerivativesSimpleComplex(
    const int x,
    const int y,
    const int sum_above_zero,
    const int sum_above_one);

  const TwoDigitNewick CreateTwoDigitNewickDerivativesComplexComplex(
    const int x,
    const int y,
    const int sum_above_zero,
    const int sum_above_one);

  ///GetDeltaSumAboveZero calculates the delta in the
  ///TwoDigitNewick::m_sum_above_zero of a new Newick
  ///when an old_value is changed.
  int GetDeltaSumAboveZero(const int old_value) const;

  ///GetDeltaSumAboveOne calculates the delta in the
  ///TwoDigitNewick::m_sum_above_one of a new Newick
  ///when an old_value is changed.
  int GetDeltaSumAboveOne(const TwoDigitNewickDerivative& d) const;
  //int GetDeltaSumAboveOne(const int old_value) const;


  ///IsSimple determines if a Newick at index i is simple.
  ///A Newick is simple if F=(X,Y) where X and Y are
  ///less than the reserved index
  bool IsSimple(const int i) const;

  ///Feed the (x,y) indices of a Newick and obtain the
  ///summarized index.
  ///Feed also obtains the derived Newicks.
  ///Note that Feed must be fed simple Newicks first
  ///int Feed(const int x, const int y)
  int SummarizeNewick(
    const int x, const int y,
    const int sum_above_zero,
    const int sum_above_one);

  ///TryToCalculateNewNewick tries to calculate the probability
  ///of Newick with index i
  void TryToCalculateNewNewick(const int i);

};

///IsSimple determines if a binary Newick is simple,
//that is in the form '(X,Y)', where both X and Y are values
//bool IsSimple(const BinaryNewickVector& n);

} //namespace ribi

#endif // TWODIGITNEWICKINDEXER_H

 

 

 

 

 

./CppTwoDigitNewick/twodigitnewickindexer.cpp

 

//---------------------------------------------------------------------------
/*
TestTwoDigitNewick, tool to test the two-digit-Newick architecture
Copyright (C) 2010-2015 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/ToolTestTwoDigitNewick.htm
//---------------------------------------------------------------------------
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Weffc++"
#include "twodigitnewickindexer.h"

#include <algorithm>
#include <cassert>
#include <iostream>
#include <stack>


#include <boost/numeric/conversion/cast.hpp>

#include "newick.h"
#include "binarynewickvector.h"

#pragma GCC diagnostic pop

//TwoDigitNewickIndexer constructor does all the work
ribi::TwoDigitNewickIndexer::TwoDigitNewickIndexer(
  const BinaryNewickVector& n,
  const double theta)
  : m_calculated_to_index{CalculateReserved(n)},
    m_current_index{CalculateReserved(n)},
    m_index_table{},
    m_newicks{CalculateReserved(n),theta},
    m_probability{-1.0},
    m_reserved{CalculateReserved(n)},
    m_theta{theta}
{
  assert(m_reserved == m_newicks.Size());
  assert(m_current_index == m_newicks.Size());

  //If the Newick is simple
  if (Newick().IsSimple(n.Peek()))
  {
    //Calculate the Ewens probability only
    m_probability = Newick().CalcProbabilitySimpleNewick(n.Peek(),theta);
    return;
  }

  //Assume all reserved Newicks from index 2 are complete
  #ifndef NDEBUG
  {
    const int sz = m_newicks.Size();
    for (int i=2; i!=sz; ++i)
    {
      assert(GetNewick(i).IsComplete());
    }
  }
  #endif

  //Create all Newicks and derivatives, but do not calculate their
  //probabilities yet
  std::vector<int> v = n.Peek();
  while(v.size()>2) //Find a leaf and cut it until the string is empty
  {
    //Find a leaf
    //Find index i (starting opening bracket) and j (closing bracket)
    const std::size_t sz = v.size();
    std::size_t i = 0;
    std::size_t j = 0;
    for (i=0 ; i!=sz; ++i) //Index of opening bracket
    {
      if (v[i]!=Newick::bracket_open) continue;
      for (j=i+1; j!=sz; ++j)
      {
        if (v[j]==Newick::bracket_open) { j = 0; break; }
        if (v[j]!=Newick::bracket_close) continue;
        break;
      }
      if (j ==  0) continue; //j cannot be 0 after previous for loop
      break;
    }
    //Find simplest leaf
    assert(v[i]==Newick::bracket_open);
    assert(v[j]==Newick::bracket_close);
    std::vector<int> v_new(v.begin(),v.begin() + i);
    const int x = v[i+1];
    const int y = v[i+2];
    //Feed it and obtain simpler index
    //Indices 0 and 1 cannot be complete, because
    //they don't have sensible derivates to point to
    assert(x < 2 || this->GetNewick(x).IsComplete());
    assert(y < 2 || this->GetNewick(y).IsComplete());
    const int saz
      = this->GetNewick(x).GetSumTermsAboveZero()
      + this->GetNewick(y).GetSumTermsAboveZero();
    const int sao
      = this->GetNewick(x).GetSumTermsAboveOne()
      + this->GetNewick(y).GetSumTermsAboveOne();

    assert(saz >  0);
    assert(sao >= 0);

    v_new.push_back(SummarizeNewick(x,y,saz,sao));
    //Replace leaf with simpler index
    std::copy(v.begin() + j + 1, v.end(),std::back_inserter(v_new));
    v = v_new;
  }

  //Now all Newicks are created, but do not have their probabilities calculated
  const int sz = m_newicks.Size();
  //m_calculated_to_index == sz denotes that all Newicks'
  //probabilities are calculated
  while (m_calculated_to_index != sz)
  {
    //Try to calculate if new Newicks can be solved
    this->TryToCalculateNewNewicks();
    //If all is well, there will be new Newicks known
    this->UpdateCalculatedFromIndex();
    //When no new probabilities are calculated,
    //all Newicks are solved
  }
}

int ribi::TwoDigitNewickIndexer::CalculateReserved(const BinaryNewickVector& n) const
{
  const std::vector<int>& v = n.Peek();
  //Count the number of elements
  const int n_elements
    = std::count_if(v.begin(),v.end(),
      std::bind2nd(std::greater<int>(),0));
  const int max_element = *std::max_element(std::begin(v),std::end(v));
  //\todo: +1 needed?
  return n_elements + max_element + 1;
}

///ConstructNewick constructs a full BinaryNewickVector from
///the TwoDigitNewick at index i.
///ConstructNewick is for debugging purposes only,
///as it's the idea to work with non-full determined
///(that is: two-digit) Newicks
const ribi::BinaryNewickVector ribi::TwoDigitNewickIndexer::ConstructNewick(const int i) const
{
  std::vector<int> v;

  if (i < m_reserved)
  {
    v.push_back(Newick::bracket_open);
    //Newick '(0)' is not valid, so fake it as '(1)'
    v.push_back(i == 0 ? 1 : i);
    v.push_back(Newick::bracket_close);
    BinaryNewickVector n(v);
    return n;
  }

  //Search for index i in Indextable to get two digits
  {
    const std::pair<int,int> p = m_index_table.Find(i);
    v.push_back(Newick::bracket_open);
    v.push_back(p.first);
    v.push_back(p.second);
    v.push_back(Newick::bracket_close);
  }
  assert(Newick().IsNewick(v));

  //As long as there are not only reserved (that is: simple)
  //values in v, replace those by their simplers
  while (1)
  {
    const std::vector<int>::iterator i
      = std::find_if(v.begin(),v.end(),
          std::bind2nd(std::greater_equal<int>(),m_reserved));
    if (i== v.end()) break;
    assert(*i >= m_reserved);
    //Create a new std::vector from the v's begin to i
    std::vector<int> v_new( v.begin(),i);
    v_new.push_back(Newick::bracket_open);
    const std::pair<int,int> p = m_index_table.Find(*i);
    v_new.push_back(p.first);
    v_new.push_back(p.second);
    v_new.push_back(Newick::bracket_close);
    //Copy the remainder of v (from after i) to v_new
    std::copy(i+1,v.end(),std::back_inserter(v_new));
    std::swap(v,v_new);
    assert(Newick().IsNewick(v));
  }
  BinaryNewickVector n(v);
  return n;
}

const ribi::TwoDigitNewick ribi::TwoDigitNewickIndexer::CreateTwoDigitNewickDerivatives(
  const int x, const int y,
  const int sum_above_zero,
  const int sum_above_one)
{
  //std::clog << "Sum_above_one: " << sum_above_one << '\n';
  assert(sum_above_zero >= 0);
  assert(sum_above_one >= 0);
  if (IsSimple(x))
  {
    if (IsSimple(y))
    {
      return
        CreateTwoDigitNewickDerivativesSimpleSimple(x,y);
    }
    return
      CreateTwoDigitNewickDerivativesSimpleComplex(
        x,y,sum_above_zero,sum_above_one);
  }
  if (IsSimple(y))
  {
    return
      CreateTwoDigitNewickDerivativesSimpleComplex(
        y,x,sum_above_zero,sum_above_one);
  }
  return
    CreateTwoDigitNewickDerivativesComplexComplex(
      x,y,sum_above_zero,sum_above_one);
}

const ribi::TwoDigitNewick
  ribi::TwoDigitNewickIndexer::CreateTwoDigitNewickDerivativesSimpleSimple(
  const int x, const int y)
{
  assert(x <= y);
  assert(x <= m_reserved && y <= m_reserved);
  assert(IsSimple(x) && IsSimple(y));


  if (x==1 && y==1)
  {
    //(1,1) -> 2
    //'2' has reserved index 2
    std::vector<TwoDigitNewickDerivative> v;
    const int derived_index = 2;
    const int value_changed = 1;
    const int other_value_changed = 1;
    //Add this derivative twice, because there are two ways
    //to change (1,1) -> 2
    v.push_back(
      TwoDigitNewickDerivative(
        derived_index,value_changed,other_value_changed));
    v.push_back(
      TwoDigitNewickDerivative(
        derived_index,value_changed,other_value_changed));
    //saz = sum above zero
    const int saz = 2;
    //sao = sum above one
    const int sao = 0;
    TwoDigitNewick n(v,saz,sao);
    return n;
  }
  else if (x==1)
  {
    //(1,y) -> { (1,y-1), (y+1) }
    //'1','y' and 'y-1' are reserved indices
    assert(x == 1);
    assert(y > 1 && y < m_reserved);
    std::vector<TwoDigitNewickDerivative> v;
    //Create (1,y-1)
    {
      //saz = sum above zero
      const int saz = x + y - 1;
      //sao = sum above one
      const int sao = (y - 1 == 1 ? 0 : y - 1);
      assert(saz >  0);
      assert(sao >= 0);
      const int d_i = SummarizeNewick(x,y-1,saz,sao);
      const int value_changed = y;
      const int other_value_changed = 0; //<Only y changed
      v.push_back(TwoDigitNewickDerivative(d_i,value_changed,other_value_changed));
    }
    //Create (y+1)
    {
      ///\todo: this
      const int d_i = y + 1;
      const int value_changed = 1;       //Note the reversal of 1 and y
      const int other_value_changed = y; //Note the reversal of 1 and y
      v.push_back(TwoDigitNewickDerivative(d_i,value_changed,other_value_changed));
    }
    TwoDigitNewick n( v,x+y,y);
    return n;
  }
  else
  {
    //(x,y) -> { (x-1,y), (x,y-1) }
    //'y','y-1','x' and 'x-1' are reserved indices
    assert(x > 1 && x < m_reserved);
    assert(y > 1 && y < m_reserved);

    std::vector<TwoDigitNewickDerivative> v;

    //saz = sum above zero
    const int saz = x + y - 1;
    //sao = sum above one
    const int sao_left  = (x - 1 == 1 ? x + y - 2 : x + y - 1);
    const int sao_right = (y - 1 == 1 ? x + y - 2 : x + y - 1);

    assert(saz       >  0);
    assert(sao_left  >= 0);
    assert(sao_right >= 0);
    //Derive (x-1,y)
    {
      const int d_i_left  = SummarizeNewick(x-1,y,saz,sao_left );
      const int value_changed_left = x;
      const int other_value_changed_left = 0; //<Only x changed
      v.push_back(
        TwoDigitNewickDerivative(
          d_i_left,value_changed_left,other_value_changed_left));
    }
    //Derive (x,y-1)
    {
      const int d_i_right = SummarizeNewick(x,y-1,saz,sao_right);
      const int value_changed_right = y;
      const int other_value_changed_right = 0; //<Only y changed
      v.push_back(
        TwoDigitNewickDerivative(
          d_i_right,value_changed_right,other_value_changed_right));
      TwoDigitNewick n(v,x+y,x+y);
      return n;
    }
  }
}

const ribi::TwoDigitNewick ribi::TwoDigitNewickIndexer::CreateTwoDigitNewickDerivativesSimpleComplex(
  const int x, const int y,
  const int sum_above_zero,
  const int sum_above_one)
{
  assert(IsSimple(x) && !IsSimple(y));

  std::vector<TwoDigitNewickDerivative> v;
  assert(x>=0);
  if (x>1)
  {
    assert(x > 1);
    assert(y >= m_reserved && "So cannot work with y-1");
    //saz = sum above zero
    const int saz = sum_above_zero - 1;
    assert(saz >= 0);
    //sao = sum above one
    const int sao = (x - 1 == 1 ? sum_above_one - 2 : sum_above_one - 1);

    assert(saz >  0);
    assert(sao >= 0);

    const int i = SummarizeNewick(x-1,y,saz,sao);
    assert(i < m_newicks.Size());

    const int value_changed = x;
    const int other_value_changed = 0; //<Only x changed
    v.push_back(TwoDigitNewickDerivative(i,value_changed,other_value_changed));
  }
  assert(y < m_newicks.Size());
  assert(!m_newicks.Empty(y));
  //\todo: Find out why 'const TwoDigitNewickDerivativesData&' does not work
  const TwoDigitNewick v_derived = m_newicks.GetNewick(y);
  for(const TwoDigitNewickDerivative& i: v_derived.GetDerivatives())
  {
    assert(i.m_derived_index < m_newicks.Size()
      && "TwoDigitNewickDerivative index must be smaller than the number of derivatives");

    //dsaz = delta sum above zero
    const int dsaz = GetDeltaSumAboveZero(i.m_value_changed);
    //dsao = delta sum above one
    const int dsao = GetDeltaSumAboveOne(i);
    //saz = sum above zero
    const int saz = sum_above_zero + dsaz;
    //sao = sum above one
    const int sao = sum_above_one + dsao;

    assert(saz >  0);
    assert(sao >= 0);

    const int d_i = SummarizeNewick( x,i.m_derived_index,saz,sao);

    const int value_changed = i.m_value_changed;
    ///\todo: check that guess is right, that using
    ///i.m_other_value_changed is better than '0'
    const int other_value_changed = i.m_other_value_changed;
    v.push_back(TwoDigitNewickDerivative(d_i,value_changed,other_value_changed));
  }
  TwoDigitNewick n( v,sum_above_zero,sum_above_one);
  return n;
}

const ribi::TwoDigitNewick
  ribi::TwoDigitNewickIndexer::CreateTwoDigitNewickDerivativesComplexComplex(
  const int x, const int y,
  const int sum_above_zero,
  const int sum_above_one)
{
  assert(!IsSimple(x) && !IsSimple(y));
  std::vector<TwoDigitNewickDerivative> v;

  //Get (X',Y)
  {
    assert(x < boost::numeric_cast<int>(m_newicks.Size()));
    assert(!m_newicks.Empty(x));
    //\todo: Find out why 'const TwoDigitNewickDerivativesData&' does not work
    const TwoDigitNewick v_derived = m_newicks.GetNewick(x);
    for(const TwoDigitNewickDerivative& i: v_derived.GetDerivatives())
    {
      //dsaz = delta sum above zero
      const int dsaz = GetDeltaSumAboveZero(i.m_value_changed);
      //dsao = delta sum above one
      const int dsao = GetDeltaSumAboveOne(i);
      //saz = sum above zero
      const int saz = sum_above_zero + dsaz;
      //sao = sum above one
      const int sao = sum_above_one + dsao;

      assert(saz >  0);
      assert(sao >= 0);

      const int d_i = SummarizeNewick(y,i.m_derived_index,saz,sao);

      const int value_changed = i.m_value_changed;
      ///\todo: check that guess is right, that using
      ///i.m_other_value_changed is better than '0'
      const int other_value_changed = i.m_other_value_changed;

      v.push_back(TwoDigitNewickDerivative(d_i,value_changed,other_value_changed));
    }
  }
  //Get (X,Y')
  {
    assert(y < boost::numeric_cast<int>(m_newicks.Size()));
    assert(!m_newicks.Empty(y));
    //\todo: Find out why 'const TwoDigitNewickDerivativesData&' does not work
    const TwoDigitNewick v_derived = m_newicks.GetNewick(y);
    for(const TwoDigitNewickDerivative& i: v_derived.GetDerivatives())
    {
      //dsaz = delta sum above zero
      const int dsaz = GetDeltaSumAboveZero(i.m_value_changed);
      //dsao = delta sum above one
      const int dsao = GetDeltaSumAboveOne(i);
      //saz = sum above zero
      const int saz = sum_above_zero + dsaz;
      //sao = sum above one
      const int sao = sum_above_one + dsao;

      assert(saz >  0);
      assert(sao >= 0);

      const int d_i = SummarizeNewick(x,i.m_derived_index,saz,sao);
      const int value_changed = i.m_value_changed;
      ///\todo: check that guess is right, that using
      ///i.m_other_value_changed is better than '0'
      const int other_value_changed = i.m_other_value_changed;

      v.push_back(TwoDigitNewickDerivative(d_i,value_changed,other_value_changed));
    }
  }
  TwoDigitNewick n(v,sum_above_zero,sum_above_one);
  return n;
}

///GetDeltaSumAboveZero calculates the delta in the
///TwoDigitNewick::m_sum_above_zero of a new Newick
///when an old_value is changed.
int ribi::TwoDigitNewickIndexer::GetDeltaSumAboveZero(const int old_value) const
{
  assert(old_value > 0);
  return (old_value == 1 ? 0 : -1);
}

///GetDeltaSumAboveOne calculates the delta in the
///TwoDigitNewick::m_sum_above_one of a new Newick
///when an old_value is changed.
int ribi::TwoDigitNewickIndexer::GetDeltaSumAboveOne(const TwoDigitNewickDerivative& d) const
{
  const int x = d.m_value_changed;
  //std::clog << "GetDeltaSumAboveOne for x " << x << '\n';
  assert(x > 0);
  if (x  > 2) return -1;
  if (x == 2) return -2;
  assert(x == 1);
  return (d.m_other_value_changed == 1 ? 2 : 1);
}

///GetProbability returns the probability of the BinaryNewickVector
///given at the constructor
double ribi::TwoDigitNewickIndexer::GetProbability() const
{
  assert(m_probability >= 0.0);
  assert(m_probability <= 1.00001);
  assert(m_calculated_to_index == m_newicks.Size()
    && "Assume calculation is completed");

  return m_probability;
}

///IsSimple determines if an index is the index of
///a simple Newick
bool ribi::TwoDigitNewickIndexer::IsSimple(const int i) const
{
  return i < m_reserved;
}

///Allow for recursion
int ribi::TwoDigitNewickIndexer::SummarizeNewick(
  const int x, const int y,
  const int sum_above_zero,
  const int sum_above_one)
{
  assert(sum_above_zero >  0);
  assert(sum_above_one  >= 0);

  assert(x > 0 && y > 0);
  //Ensure proper ordering
  if (x > y)
  {
    return SummarizeNewick(
      y,x,sum_above_zero,sum_above_one);
  }
  assert(x <= y);

  //If index is known, return the index
  if (m_index_table.CanGetIndex(x,y))
  {
    const int i = m_index_table.GetIndex(x,y);
    if (i) return i;
  }

  //A new (x,y) pair is fed
  assert(x <= y);
  const int i = m_current_index;
  m_index_table.SetIndex(x,y,i);

  ++m_current_index;

  assert(sum_above_zero >= 0);
  assert(sum_above_one  >= 0);

  TwoDigitNewick n = CreateTwoDigitNewickDerivatives(
    x,y,sum_above_zero,sum_above_one);

  //Check if the Newick is simple
  //and its probability can be calculated with
  //the Ewens formula
  if (x < m_reserved && y < m_reserved)
  {
    std::vector<int> v;
    v.push_back(Newick::bracket_open);
    v.push_back(x);
    v.push_back(y);
    v.push_back(Newick::bracket_close);
    n.SetProbability(Newick().CalcProbabilitySimpleNewick(v,m_theta));
    //If the user requests simple Newicks to be solved,
    //perhaps the requested probability has just been calculated
    m_probability = n.GetProbability();
  }
  m_newicks.SetNewick(i,n);
  return i;
}

///TryToCalculateNewNewick tries to calculate the probability
///of Newick with index i
void ribi::TwoDigitNewickIndexer::TryToCalculateNewNewick(const int i)
{
  assert(i >= m_calculated_to_index);
  assert(i  < m_newicks.Size());
  ///\todo: why cannot use 'const TwoDigitNewick&'?
  ///\bug: this seems to give memory problems (and funny output)
  const TwoDigitNewick n = this->GetNewick(i);
  //Remember: m_calculated_to_index is not increased
  if (n.IsProbabilityKnown())
  {
    //Already know the probability
    return;
  }
  const std::vector<TwoDigitNewickDerivative> derivatives = n.GetDerivatives();
  //Check if of all derivates the probability is known
  for(const TwoDigitNewickDerivative& derivative: derivatives)
  {
    if (!GetNewick(derivative.m_derived_index).IsProbabilityKnown())
    {
      //Too bad, derived Newick's probability is unknown
      return;
    }
  }
  ///Calculate constants and probability
  double p = 0.0;
  //The denominator is that of the focal Newick,
  //and not its derivative(s)
  const double denominator = GetNewick(i).GetDenominator();
  for(const TwoDigitNewickDerivative& derivative: derivatives)
  {
    assert(GetNewick(derivative.m_derived_index).IsProbabilityKnown());

    const double c
      = (derivative.m_value_changed == 1
      ? m_theta
      : boost::numeric_cast<double>(derivative.m_value_changed) * boost::numeric_cast<double>(derivative.m_value_changed - 1))
      / denominator;
    const double p_this
      = c * GetNewick(derivative.m_derived_index).GetProbability();
    p+=p_this;
  }
  //std::clog << "Index " << i
  //  << " has been calculated to have the probability of "
  //  << p << '\n';
  m_probability = p;
  this->m_newicks.SetNewickProbability(i,p);
}

///TryToCalculateNewNewicks tries to calculate new Newick probabilities
///\warning: m_calculated_to_index is not increased
void ribi::TwoDigitNewickIndexer::TryToCalculateNewNewicks()
{
  const int sz = m_newicks.Size();
  //m_calculated_to_index == sz denotes that all Newicks'
  //probabilities are calculated
  assert(m_calculated_to_index <= sz);
  for (int i = m_calculated_to_index; i!=sz; ++i)
  {
    TryToCalculateNewNewick(i);
  }
}

void ribi::TwoDigitNewickIndexer::UpdateCalculatedFromIndex()
{
  const int sz = m_newicks.Size();
  while (m_calculated_to_index < sz)
  {
    assert(m_calculated_to_index < m_newicks.Size());
    ///\todo: find out why 'const TwoDigitNewick&' does not work
    const TwoDigitNewick n
      = m_newicks.GetNewick(m_calculated_to_index);
    if (n.IsProbabilityKnown())
    {
      //std::clog << "At index " << m_calculated_to_index
      //  << " the probability is known\n";
      ++m_calculated_to_index;
    }
    else
    {
      break;
    }
  }
}

 

 

 

 

 

./CppTwoDigitNewick/twodigitnewickindextable.h

 

//---------------------------------------------------------------------------
/*
TestTwoDigitNewick, tool to test the two-digit-Newick architecture
Copyright (C) 2010-2015 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/ToolTestTwoDigitNewick.htm
//---------------------------------------------------------------------------
#ifndef TWODIGITNEWICKINDEXERTABLE_H
#define TWODIGITNEWICKINDEXERTABLE_H

#include <vector>

namespace ribi {

///TwoDigitNewickIndexTable manages (X,Y) -> index
struct TwoDigitNewickIndexTable
{
  ///CanGetIndex returns if GetIndex can be called with success
  bool CanGetIndex(const int x, const int y) const;
  ///Find returns the x-y-coordinats of the Newick with index i
  const std::pair<int,int> Find(const int i) const;
  ///Get returns the internals of TwoDigitNewickIndexTable
  const std::vector<std::vector<int> >& Get() const { return m_v; }
  ///GetIndex returns m_index_table[x][y]
  int GetIndex(const int x, const int y) const;
  ///GetNumAllocated calculates the number of indices allocated
  int GetNumAllocated() const;
  ///GetNumUsed calculates the number of indices used
  int GetNumUsed() const;
  ///Set (x,y) to index z
  void SetIndex(const int x, const int y, const int z);
  private:
  ///m_index_table is the index table that maps (x,y) to a value,
  ///so that m_index_table[x][y] equals that value
  std::vector<std::vector<int> > m_v;
};

} //~namespace ribi

#endif // TWODIGITNEWICKINDEXERTABLE_H

 

 

 

 

 

./CppTwoDigitNewick/twodigitnewickindextable.cpp

 

//---------------------------------------------------------------------------
/*
TestTwoDigitNewick, tool to test the two-digit-Newick architecture
Copyright (C) 2010-2015 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/ToolTestTwoDigitNewick.htm
//---------------------------------------------------------------------------
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Weffc++"
#include "twodigitnewickindextable.h"


#include <algorithm>
#include <cassert>
#include <stdexcept>
#include <iostream>
#include <sstream>


#include <boost/numeric/conversion/cast.hpp>

#pragma GCC diagnostic pop

///CanGetData returns is GetData can be called with success
bool ribi::TwoDigitNewickIndexTable::CanGetIndex(const int x, const int y) const
{
  return x < boost::numeric_cast<int>(m_v.size())
      && y < boost::numeric_cast<int>(m_v[x].size());
}

///Find returns the x-y-coordinats of the Newick with index i
const std::pair<int,int> ribi::TwoDigitNewickIndexTable::Find(const int i) const
{
  const int maxx = boost::numeric_cast<int>(m_v.size());
  for (int x=0; x!=maxx; ++x)
  {
    const int maxy = boost::numeric_cast<int>(m_v[x].size());
    for (int y=0; y!=maxy; ++y)
    {
      if (m_v[x][y]==i) return std::make_pair(x,y);
    }
  }
  std::cerr << "Index " << i << " not found in ribi::TwoDigitNewickIndexTable::Find\n";
  assert(!"Should not get here");
  std::stringstream s;
  s << "Index " << i << " not found in ribi::TwoDigitNewickIndexTable::Find";
  throw std::logic_error(s.str());
}

///GetIndex returns m_index_table[x][y]
int ribi::TwoDigitNewickIndexTable::GetIndex(const int x, const int y) const
{
  assert(CanGetIndex(x,y));
  return m_v[x][y];
}

///GetNumAllocated calculates the number of indices allocated
int ribi::TwoDigitNewickIndexTable::GetNumAllocated() const
{
  int n_allocated = 0;
  for(const std::vector<int>& v: m_v)
  {
    n_allocated+=boost::numeric_cast<int>(v.size());
  }
  return n_allocated;
}

///GetNumUsed calculates the number of indices used
int ribi::TwoDigitNewickIndexTable::GetNumUsed() const
{
  int n_non_zero = 0;
  for(const std::vector<int>& v: m_v)
  {
    n_non_zero +=std::count_if(v.begin(),v.end(),
      std::bind2nd(std::greater<int>(),0));
  }
  return n_non_zero;
}

///SetIndex sets m_index_table[x][y] == z
///and resized the m_index_table is necessary
void ribi::TwoDigitNewickIndexTable::SetIndex(
  const int x,
  const int y,
  const int z)
{
  //Does x have a proper std::vector?
  if (x >= boost::numeric_cast<int>(m_v.size()))
  {
    m_v.resize(x + 1,std::vector<int>());
  }
  if (y >= boost::numeric_cast<int>(m_v[x].size()))
  {
    m_v[x].resize(y + 1,0);
  }
  assert(m_v[x][y] == 0);
  assert(z > 0);
  m_v[x][y] = z;
  assert(m_v[x][y] > 0);
}

 

 

 

 

 

./CppTwoDigitNewick/twodigitnewicks.h

 

//---------------------------------------------------------------------------
/*
TestTwoDigitNewick, tool to test the two-digit-Newick architecture
Copyright (C) 2010-2015 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/ToolTestTwoDigitNewick.htm
//---------------------------------------------------------------------------
#ifndef TWODIGITNEWICKS_H
#define TWODIGITNEWICKS_H

#include <vector>

#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Weffc++"
#include "twodigitnewick.h"
#pragma GCC diagnostic pop

namespace ribi {

///TwoDigitNewicks manages TwoDigitTwoDigitNewicks.
///TwoDigitTwoDigitNewicks manage F'(index) -> { indices }
struct TwoDigitNewicks
{
  TwoDigitNewicks(const int n_reserved, const double theta);
  ///Empty returns if m_v is empty
  ///Set a Newick at index i
  ///\warning: do not use push_back, because of recursive calls
  void SetNewick(const int i, const TwoDigitNewick& v);
  bool Empty() const { return m_v.empty(); }
  ///Empty returns if an index is empty
  bool Empty(const int i) const;
  ///Get returns the internals of TwoDigitNewicks
  const std::vector<TwoDigitNewick>& Get() const { return m_v; }
  ///GetNewick returns the Newick at index i.
  ///i is checked for its range.
  const TwoDigitNewick& GetNewick(const int i) const;
  ///Sets the derivatives of F(i) to { v }
  //void SetDerivatives(const int i,const TwoDigitNewick& v);
  ///Obtain the number of Derivatives

  ///SetNewickProbability sets the probability of the Newick
  ///at index i to p
  void SetNewickProbability(const int i,const double p);

  int Size() const;
private:
  ///m_derivatives contains for every indexed Newick
  ///its derived indexed TwoDigitNewicks
  std::vector<TwoDigitNewick> m_v;
};

} //~namespace ribi

#endif // TWODIGITNEWICKS_H

 

 

 

 

 

./CppTwoDigitNewick/twodigitnewicks.cpp

 

//---------------------------------------------------------------------------
/*
TestTwoDigitNewick, tool to test the two-digit-Newick architecture
Copyright (C) 2010-2015 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/ToolTestTwoDigitNewick.htm
//---------------------------------------------------------------------------
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Weffc++"
#include "twodigitnewicks.h"

#include <cassert>


#include <boost/numeric/conversion/cast.hpp>

#include "newick.h"
#include "binarynewickvector.h"

#pragma GCC diagnostic pop

ribi::TwoDigitNewicks::TwoDigitNewicks(const int n_reserved, const double theta)
  : m_v{}
{
  //Create derivatives of simplest and reserved TwoDigitNewicks

  //Create the Newick at index 0
  //this->SetNewick(0,TwoDigitNewick(std::vector<TwoDigitNewickDerivative>(TwoDigitNewickDerivative(),0,0));

  //Let i start at 1, because '(0)' is no valid Newick
  for (int i=1; i!=n_reserved; ++i)
  {
    std::vector<TwoDigitNewickDerivative> v;
    if (i > 1)
    {
      const int value_changed = i;
      const int other_value_changed = 0; //<There is no other value
      v.push_back(TwoDigitNewickDerivative(i-1,value_changed,other_value_changed));
    }
    //saz = sum above zero
    const int saz = i;
    //sao = sum above one
    const int sao = (i > 1 ? i : 0);

    TwoDigitNewick n(v,saz,sao);

    assert(i!=0 && "'(0)' is no valid Newick");

    n.SetProbability(Newick().CalcProbabilitySimpleNewick(
      Newick().CreateVector(
        static_cast<int>(Newick::bracket_open),
        i,
        static_cast<int>(Newick::bracket_close)),
      theta)
    );
    this->SetNewick(i,n);

    assert( (i < 2 || this->GetNewick(i).IsComplete())
     && "All newick with index >= 2 must be complete");
  }
}

///Empty returns if an index is empty
bool ribi::TwoDigitNewicks::Empty(const int i) const
{
  assert(i >= 0);
  assert(i < this->Size());
  return m_v[i].Empty();
}

const ribi::TwoDigitNewick& ribi::TwoDigitNewicks::GetNewick(
  const int i) const
{
  //Check if i is in range
  assert(i>=0);
  assert(i < this->Size());
  //Check if returned indices are okay
  #ifndef NDEBUG
  const TwoDigitNewick& v = m_v[i];
  for(const TwoDigitNewickDerivative& j: v.GetDerivatives())
  {
    assert(j.m_derived_index >= 0);
    assert(j.m_derived_index < this->Size() );
  }
  #endif
  return m_v[i];
}

void ribi::TwoDigitNewicks::SetNewick(const int i, const TwoDigitNewick& v)
{
  //std::clog << __LINE__ << " - " << i << '\n';
  //Allocate storage
  //TODO: replace by push_back
  if (i >= boost::numeric_cast<int>(m_v.size()))
  {
    m_v.resize(i + 1);
    assert(m_v[m_v.size()-1].Empty());
  }
  //
  //std::clog << "Adding a Newick at index " << i << '\n';
  assert(m_v[i].Empty());

  #ifndef NDEBUG
  for(const TwoDigitNewickDerivative& j: v.GetDerivatives())
  {
    assert(j.m_derived_index > 0);
    assert(j.m_derived_index < boost::numeric_cast<int>(m_v.size())
      && "Cannot set a derivative index "
         "bigger than the number of derivatives");
  }
  #endif

  m_v[i] = v;
}

void ribi::TwoDigitNewicks::SetNewickProbability(
  const int i,const double p)
{
  assert(i >= 0);
  assert(i  < Size());
  assert(p >= 0.0);
  assert(p <= 1.0);
  m_v[i].SetProbability(p);
}

int ribi::TwoDigitNewicks::Size() const
{
  return boost::numeric_cast<int>(m_v.size());
}

 

 

 

 

 

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

Go back to Richel Bilderbeek's homepage.

 

Valid XHTML 1.0 Strict

This page has been created by the tool CodeToHtml