BipartitionList.cpp

Go to the documentation of this file.

// SPDX-FileCopyrightText: The Bio++ Development Group
//
// SPDX-License-Identifier: CECILL-2.1
 
#include <Bpp/Exceptions.h>
#include <Bpp/Io/FileTools.h>
#include <Bpp/Numeric/Random/RandomTools.h>
#include <Bpp/Text/TextTools.h>
 
#include "BipartitionList.h"
#include "BipartitionTools.h"
#include "TreeTemplate.h"
 
// From SeqLib:
#include <Bpp/Seq/Alphabet/Alphabet.h>
#include <Bpp/Seq/Alphabet/AlphabetTools.h>
 
 
using namespace bpp;
 
// From the STL:
#include <iostream>
#include <climits> // defines CHAR_BIT
 
using namespace std;
 
/****************************************************************/
/* utilitary classes required for sorting elements/bipartitions */
/****************************************************************/
 
class StringAndInt
{
public:
  int ind;
  string str;
 
public:
  StringAndInt() : ind(0),
    str() {}
};
 
bool operator<(StringAndInt sai1, StringAndInt sai2)
{
  if (sai1.str < sai2.str)
    return true;
  return false;
}
 
/******************************************************************************/
 
class IntAndInt
{
public:
  size_t ind;
  int val;
};
 
bool operator<(IntAndInt iai1, IntAndInt iai2)
{
  if (iai1.val < iai2.val)
    return true;
  return false;
}
 
 
/******************************************************************************/
 
BipartitionList::BipartitionList(const Tree& tr, bool sorted, std::vector<int>* index) :
  bitBipartitionList_(),
  elements_(),
  sorted_(sorted)
{
  size_t nbbip;
 
  elements_ = tr.getLeavesNames();
 
  if (tr.isRooted())
    nbbip = tr.getNumberOfNodes() - 2;
  else
    nbbip = tr.getNumberOfNodes() - 1;
 
  if (sorted)
    std::sort(elements_.begin(), elements_.end());
 
  size_t lword  = static_cast<size_t>(BipartitionTools::LWORD);
  size_t nbword = (elements_.size() + lword - 1) / lword;
  size_t nbint  = nbword * lword / (CHAR_BIT * sizeof(int));
 
  for (size_t i = 0; i < nbbip; i++)
  {
    bitBipartitionList_.push_back(new int[nbint]);
    for (size_t j = 0; j < nbint; j++)
    {
      bitBipartitionList_[i][j] = 0;
    }
  }
 
  size_t cpt = 0;
  vector<string> underlyingNames;
  const Tree* tree = &tr;
  const TreeTemplate<Node>* ttree = dynamic_cast<const TreeTemplate<Node>*>(tree);
  if (ttree)
  {
    // Gain some time...
    buildBitBipartitions(ttree->getRootNode(), bitBipartitionList_, elements_, &cpt, index);
  }
  else
  {
    TreeTemplate<Node> tmp(tr);
    buildBitBipartitions(tmp.getRootNode(), bitBipartitionList_, elements_, &cpt, index);
  }
}
 
/******************************************************************************/
 
BipartitionList::BipartitionList(
    const std::vector<std::string>& elements,
    const std::vector<int*>& bitBipL) :
  bitBipartitionList_(),
  elements_(elements),
  sorted_()
{
  size_t lword  = static_cast<size_t>(BipartitionTools::LWORD);
  size_t nbword = (elements.size() + lword - 1) / lword;
  size_t nbint  = nbword * lword / (CHAR_BIT * sizeof(int));
 
  for (size_t i = 0; i < bitBipL.size(); i++)
  {
    bitBipartitionList_.push_back(new int[nbint]);
    for (size_t j = 0; j < nbint; j++)
    {
      bitBipartitionList_[i][j] = bitBipL[i][j];
    }
  }
 
  vector<string> cpelements_ = elements;
  std::sort(cpelements_.begin(), cpelements_.end());
  if (cpelements_ == elements)
    sorted_ = true;
  else
    sorted_ = false;
}
 
/******************************************************************************/
 
BipartitionList::BipartitionList(const BipartitionList& bipL) :
  bitBipartitionList_(),
  elements_(bipL.elements_),
  sorted_(bipL.sorted_)
{
  size_t lword  = static_cast<size_t>(BipartitionTools::LWORD);
  size_t nbword = (bipL.getNumberOfElements() + lword - 1) / lword;
  size_t nbint  = nbword * lword / (CHAR_BIT * sizeof(int));
 
  bitBipartitionList_.resize(bipL.getNumberOfBipartitions());
  vector<int*> bitBipL = bipL.getBitBipartitionList();
  for (size_t i = 0; i < bipL.getNumberOfBipartitions(); i++)
  {
    bitBipartitionList_[i] = new int[nbint];
    for (size_t j = 0; j < nbint; j++)
    {
      bitBipartitionList_[i][j] = bitBipL[i][j];
    }
  }
}
 
/******************************************************************************/
 
BipartitionList& BipartitionList::operator=(const BipartitionList& bipL)
{
  size_t lword  = static_cast<size_t>(BipartitionTools::LWORD);
  size_t nbword = (bipL.getNumberOfElements() + lword - 1) / lword;
  size_t nbint  = nbword * lword / (CHAR_BIT * sizeof(int));
 
  for (size_t i = 0; i < bitBipartitionList_.size(); i++)
  {
    delete[] bitBipartitionList_[i];
  }
  bitBipartitionList_.resize(bipL.getNumberOfBipartitions());
  vector<int*> bitBipL = bipL.getBitBipartitionList();
  for (size_t i = 0; i < bipL.getNumberOfBipartitions(); i++)
  {
    bitBipartitionList_[i] = new int[nbint];
    for (size_t j = 0; j < nbint; j++)
    {
      bitBipartitionList_[i][j] = bitBipL[i][j];
    }
  }
 
  elements_ = bipL.elements_;
  sorted_   = bipL.sorted_;
  return *this;
}
 
/******************************************************************************/
 
BipartitionList::~BipartitionList()
{
  for (size_t i = 0; i < bitBipartitionList_.size(); i++)
  {
    delete[] bitBipartitionList_[i];
  }
}
 
/******************************************************************************/
 
map<string, bool> BipartitionList::getBipartition(size_t i) const
{
  map<string, bool> bip;
 
  if (i >= bitBipartitionList_.size())
    throw Exception("Bipartition index exceeds BipartitionList size");
 
  for (size_t j = 0; j < elements_.size(); j++)
  {
    if (BipartitionTools::testBit(bitBipartitionList_[i], static_cast<int>(j)))
      bip[elements_[j]] = true;
    else
      bip[elements_[j]] = false;
  }
  return bip;
}
 
/******************************************************************************/
 
int* BipartitionList::getBitBipartition(size_t i)
{
  if (i >= bitBipartitionList_.size())
    throw Exception("Bipartition index exceeds BipartitionList size");
 
  return bitBipartitionList_[i];
}
 
/******************************************************************************/
 
bool BipartitionList::haveSameElementsThan(map<string, bool>& bipart) const
{
  vector<string> elements = elements_;
  vector<string> keys;
 
  map<string, bool>::iterator it;
 
  for (it = bipart.begin(); it != bipart.end(); it++)
  {
    keys.push_back(it->first);
  }
 
  std::sort(elements.begin(), elements.end());
  std::sort(keys.begin(), keys.end());
 
  if (elements == keys)
    return true;
  return false;
}
 
/******************************************************************************/
 
void BipartitionList::addBipartition(map<string, bool>& bipart, bool checkElements)
{
  if (checkElements && !BipartitionList::haveSameElementsThan(bipart))
    throw Exception("Distinct bipartition element sets");
 
  size_t lword  = static_cast<size_t>(BipartitionTools::LWORD);
  size_t nbword = (elements_.size() + lword - 1) / lword;
  size_t nbint  = nbword * lword / (CHAR_BIT * sizeof(int));
  bitBipartitionList_.push_back(new int[nbint]);
  size_t ind    = bitBipartitionList_.size() - 1;
  for (size_t j = 0; j < nbint; j++)
  {
    bitBipartitionList_[ind][j] = 0;
  }
 
  for (size_t i = 0; i < elements_.size(); i++)
  {
    if (bipart[elements_[i]] == true)
      BipartitionTools::bit1(bitBipartitionList_[ind], static_cast<int>(i));
    else
      BipartitionTools::bit0(bitBipartitionList_[ind], static_cast<int>(i));
  }
}
 
/******************************************************************************/
 
void BipartitionList::deleteBipartition(size_t i)
{
  if (i >= bitBipartitionList_.size())
    throw Exception("Bipartition index exceeds BipartitionList size");
 
  delete[] bitBipartitionList_[i];
  bitBipartitionList_.erase(bitBipartitionList_.begin() + static_cast<ptrdiff_t>(i));
}
 
/******************************************************************************/
 
bool BipartitionList::containsBipartition(map<string, bool>& bipart, bool checkElements) const
{
  size_t i, j;
  bool dac, padac;
 
  if (checkElements && !BipartitionList::haveSameElementsThan(bipart))
    throw Exception("Distinct bipartition element sets");
 
  for (i = 0; i < bitBipartitionList_.size(); i++)
  {
    dac = padac = false;
    for (j = 0; j < elements_.size(); j++)
    {
      if (BipartitionTools::testBit(bitBipartitionList_[i], static_cast<int>(j)))
      {
        if (bipart[elements_[j]])
          dac = true;
        else
          padac = true;
      }
      else
      {
        if (bipart[elements_[j]])
          padac = true;
        else
          dac = true;
      }
      if (dac && padac)
        break;
    }
    if (j == elements_.size())
      return true;
  }
  return false;
}
 
/******************************************************************************/
 
bool BipartitionList::areIdentical(size_t k1, size_t k2) const
{
  bool dac, padac;
 
  if (k1 >= bitBipartitionList_.size())
    throw Exception("Bipartition index exceeds BipartitionList size");
  if (k2 >= bitBipartitionList_.size())
    throw Exception("Bipartition index exceeds BipartitionList size");
 
  dac = padac = false;
  for (size_t j = 0; j < elements_.size(); j++)
  {
    if (BipartitionTools::testBit(bitBipartitionList_[k1], static_cast<int>(j)))
    {
      if (BipartitionTools::testBit(bitBipartitionList_[k2], static_cast<int>(j)))
        dac = true;
      else
        padac = true;
    }
    else
    {
      if (BipartitionTools::testBit(bitBipartitionList_[k2], static_cast<int>(j)))
        padac = true;
      else
        dac = true;
    }
    if (dac && padac)
      return false;
  }
  return true;
}
 
/******************************************************************************/
 
bool BipartitionList::areCompatible(size_t k1, size_t k2) const
{
  bool uu, uz, zu, zz;
 
  if (k1 >= bitBipartitionList_.size())
    throw Exception("Bipartition index exceeds BipartitionList size");
  if (k2 >= bitBipartitionList_.size())
    throw Exception("Bipartition index exceeds BipartitionList size");
 
  uu = uz = zu = zz = false;
 
  for (size_t j = 0; j < elements_.size(); j++)
  {
    if (BipartitionTools::testBit(bitBipartitionList_[k1], static_cast<int>(j)))
    {
      if (BipartitionTools::testBit(bitBipartitionList_[k2], static_cast<int>(j)))
        uu = true;
      else
        uz = true;
    }
    else
    {
      if (BipartitionTools::testBit(bitBipartitionList_[k2], static_cast<int>(j)))
        zu = true;
      else
        zz = true;
    }
    if (uu && uz && zu && zz)
      return false;
  }
 
  return true;
}
 
/******************************************************************************/
 
bool BipartitionList::areAllCompatible() const
{
  for (size_t i = 0; i < bitBipartitionList_.size(); i++)
  {
    for (size_t j = i + 1; j < bitBipartitionList_.size(); j++)
    {
      if (!BipartitionList::areCompatible(i, j))
        return false;
    }
  }
  return true;
}
 
/******************************************************************************/
 
bool BipartitionList::areAllCompatibleWith(map<string, bool>& bipart, bool checkElements) const
{
  if (checkElements && !haveSameElementsThan(bipart))
    throw Exception("Distinct bipartition element sets");
  size_t nbBip = bitBipartitionList_.size();
  const_cast<BipartitionList*>(this)->addBipartition(bipart, false);
 
  for (size_t i = 0; i < nbBip; i++)
  {
    if (!areCompatible(i, nbBip))
    {
      const_cast<BipartitionList*>(this)->deleteBipartition(nbBip);
      return false;
    }
  }
  const_cast<BipartitionList*>(this)->deleteBipartition(nbBip);
  return true;
}
 
/******************************************************************************/
 
void BipartitionList::sortElements()
{
  vector<StringAndInt> relements_;
  StringAndInt sai;
  size_t nbbip;
 
  for (size_t i = 0; i < elements_.size(); i++)
  {
    sai.str = elements_[i];
    sai.ind = static_cast<int>(i);
    relements_.push_back(sai);
  }
 
  std::sort(relements_.begin(), relements_.end());
 
  for (size_t i = 0; i < elements_.size(); i++)
  {
    elements_[i] = relements_[i].str;
  }
 
  nbbip = bitBipartitionList_.size();
  bitBipartitionList_.resize(2 * nbbip);
  size_t lword  = static_cast<size_t>(BipartitionTools::LWORD);
  size_t nbword = (elements_.size() + lword - 1) / lword;
  size_t nbint  = nbword * lword / (CHAR_BIT * sizeof(int));
 
  for (size_t j = nbbip; j < 2 * nbbip; j++)
  {
    bitBipartitionList_[j] = new int[nbint];
    for (size_t k = 0; k < nbint; k++)
    {
      bitBipartitionList_[j][k] = 0;
    }
    for (size_t i = 0; i < elements_.size(); i++)
    {
      if (BipartitionTools::testBit(bitBipartitionList_[j - nbbip], relements_[i].ind))
        BipartitionTools::bit1(bitBipartitionList_[j], static_cast<int>(i));
      else
        BipartitionTools::bit0(bitBipartitionList_[j], static_cast<int>(i));
    }
  }
 
  for (size_t j = 0; j < nbbip; j++)
  {
    delete[] bitBipartitionList_[j];
  }
 
  bitBipartitionList_.erase(bitBipartitionList_.begin(), bitBipartitionList_.begin() + static_cast<ptrdiff_t>(nbbip));
  sorted_ = true;
}
 
/******************************************************************************/
 
size_t BipartitionList::getPartitionSize(size_t k) const
{
  size_t size = 0;
  if (k >= bitBipartitionList_.size())
    throw Exception("Bipartition index exceeds BipartitionList size");
 
  for (size_t i = 0; i < elements_.size(); i++)
  {
    if (BipartitionTools::testBit(bitBipartitionList_[k], static_cast<int>(i)))
      size++;
  }
 
  if (size <= elements_.size() / 2)
    return size;
  else
    return elements_.size() - size;
}
 
/******************************************************************************/
 
void BipartitionList::removeTrivialBipartitions()
{
  size_t size = bitBipartitionList_.size();
  for (size_t i = size; i > 0; i--)
  {
    if (BipartitionList::getPartitionSize(i - 1) < 2)
      BipartitionList::deleteBipartition(i - 1);
  }
}
 
/******************************************************************************/
 
void BipartitionList::addTrivialBipartitions(bool checkExisting)
{
  map<string, bool> bip;
 
  for (size_t i = 0; i < elements_.size(); i++)
  {
    bip[elements_[i]] = false;
  }
  for (size_t i = 0; i < elements_.size(); i++)
  {
    bip[elements_[i]] = true;
    if (checkExisting && BipartitionList::containsBipartition(bip, false))
      continue;
    BipartitionList::addBipartition(bip, false);
    bip[elements_[i]] = false;
  }
}
 
/******************************************************************************/
 
void BipartitionList::sortByPartitionSize()
{
  vector<int*> sortedBitBipL;
  vector<IntAndInt> iaiVec;
  IntAndInt iai;
 
  for (size_t i = 0; i < bitBipartitionList_.size(); i++)
  {
    iai.ind = i;
    iai.val = static_cast<int>(BipartitionList::getPartitionSize(i));
    iaiVec.push_back(iai);
  }
 
  std::sort(iaiVec.begin(), iaiVec.end());
 
  for (size_t i = 0; i < bitBipartitionList_.size(); i++)
  {
    sortedBitBipL.push_back(bitBipartitionList_[iaiVec[i].ind]);
  }
 
  bitBipartitionList_ = sortedBitBipL;
}
 
/******************************************************************************/
 
void BipartitionList::flip(size_t k)
{
  if (k >= bitBipartitionList_.size())
    throw Exception("Bipartition index exceeds BipartitionList size");
  size_t lword = static_cast<size_t>(BipartitionTools::LWORD);
  size_t nbword = (elements_.size() + lword - 1) / lword;
  size_t nbint = nbword * lword / (CHAR_BIT * sizeof(int));
  int* flipbip = new int[nbint];
  for (size_t i = 0; i < nbint; i++)
  {
    flipbip[i] = 0;
  }
  BipartitionTools::bitNot(flipbip, bitBipartitionList_[k], nbint);
  delete[] bitBipartitionList_[k];
  bitBipartitionList_[k] = flipbip;
}
 
/******************************************************************************/
 
void BipartitionList::removeRedundantBipartitions()
{
  bool deletion = true;
 
  while (deletion)
  {
    deletion = false;
    for (size_t i = 0; i < bitBipartitionList_.size(); i++)
    {
      for (size_t j = i + 1; j < bitBipartitionList_.size(); j++)
      {
        if (BipartitionList::areIdentical(i, j))
        {
          BipartitionList::deleteBipartition(j);
          deletion = true;
          break;
        }
      }
      if (deletion)
        break;
    }
  }
}
 
/******************************************************************************/
 
unique_ptr<TreeTemplate<Node>> BipartitionList::toTree() const
{
  unique_ptr<BipartitionList> sortedBipL;
  vector<int*> sortedBitBipL;
  int* bip;
  vector<Node*> vecNd, sonNd;
  vector<bool> alive;
  size_t lword, nbword, nbint, ii;
 
  /* check, copy and prepare bipartition list */
 
  if (!BipartitionList::areAllCompatible())
    throw Exception("Trying to build a tree from incompatible bipartitions");
 
  sortedBipL.reset(dynamic_cast<BipartitionList*>(clone()));
  for (size_t i = 0; i < sortedBipL->getNumberOfBipartitions(); i++)
  {
    if (sortedBipL->getPartitionSize(i) > sortedBipL->getNumberOfElements() / 2)
      sortedBipL->flip(i);
  }
  sortedBipL->sortByPartitionSize();
  sortedBipL->removeRedundantBipartitions();
  sortedBitBipL = sortedBipL->getBitBipartitionList();
 
  for (size_t i = 0; i < sortedBipL->getNumberOfBipartitions(); i++)
  {
    alive.push_back(true);
  }
  vecNd.resize(sortedBipL->getNumberOfBipartitions() + 1);
  lword  = static_cast<size_t>(BipartitionTools::LWORD);
  nbword = (elements_.size() + lword - 1) / lword;
  nbint  = nbword * lword / (CHAR_BIT * sizeof(int));
  bip    = new int[1]; bip[0] = 0;
 
  /* main loop: create one node per bipartition */
  for (size_t i = 0; i < sortedBipL->getNumberOfBipartitions(); i++)
  {
    if (sortedBipL->getPartitionSize(i) == 1)
    { // terminal
      for (size_t j = 0; j < sortedBipL->getNumberOfElements(); j++)
      {
        if (BipartitionTools::testBit(sortedBitBipL[i], static_cast<int>(j)))
        {
          vecNd[i] = new Node(elements_[j]);
          break;
        }
      }
    }
    else
    { // internal
      sonNd.clear();
      for (size_t j = 0; j < i; j++)
      {
        if (alive[j])
        {
          for (ii = 0; ii < nbint; ii++)
          {
            BipartitionTools::bitOr(bip, sortedBitBipL[j] + ii, sortedBitBipL[i] + ii, 1);
            if (bip[0] != sortedBitBipL[i][ii])
              break;
          }
          if (ii == nbint)
          {
            sonNd.push_back(vecNd[j]);
            alive[j] = false;
          }
        }
      }
      vecNd[i] = new Node();
      for (size_t k = 0; k < sonNd.size(); k++)
      {
        vecNd[i]->addSon(sonNd[k]);
      }
    }
  }
 
  /* create last node, which joins alive bipartitions = fatherless nodes */
  Node* rootNd = new Node();
  for (size_t i = 0; i < sortedBipL->getNumberOfBipartitions(); i++)
  {
    if (alive[i])
      rootNd->addSon(vecNd[i]);
  }
 
  /* construct tree and return */
  auto tr = make_unique<TreeTemplate<Node>>(rootNd);
  tr->resetNodesId();
 
  return tr;
}
 
/******************************************************************************/
 
vector<string> BipartitionList::buildBitBipartitions(const Node* nd, vector<int*>& bitbip, const vector<string>& elements, size_t* cpt, vector<int>* index) const
{
  vector<string> underelements_, retelements_;
 
  if (nd->getNumberOfSons() == 0)
    underelements_.push_back(nd->getName());
 
  for (size_t i = 0; i < nd->getNumberOfSons(); i++)
  {
    retelements_ = BipartitionList::buildBitBipartitions(nd->getSon(i), bitbip, elements, cpt, index);
    for (size_t j = 0; j < retelements_.size(); j++)
    {
      underelements_.push_back(retelements_[j]);
    }
  }
 
  if (!nd->hasFather())
    return underelements_; // root node
 
  if (!nd->getFather()->hasFather())
  {
    size_t nbrootson = nd->getFather()->getNumberOfSons();
    if (nbrootson == 2 && nd == nd->getFather()->getSon(1))
      return underelements_; // son 2 of root node when root node has 2 sons
  }
 
  bool ones;
  if (underelements_.size() <= elements.size() / 2)
    ones = true;
  else
    ones = false;
 
  for (size_t i = 0; i < elements.size(); i++)
  {
    if (ones)
      BipartitionTools::bit0(bitbip[*cpt], static_cast<int>(i));
    else
      BipartitionTools::bit1(bitbip[*cpt], static_cast<int>(i));
  }
 
  for (size_t i = 0; i < underelements_.size(); i++)
  {
    size_t taxa_ind = 0;
    while (underelements_[i] != elements[taxa_ind])
      taxa_ind++;
    if (ones)
      BipartitionTools::bit1(bitbip[*cpt], static_cast<int>(taxa_ind));
    else
      BipartitionTools::bit0(bitbip[*cpt], static_cast<int>(taxa_ind));
  }
 
  (*cpt)++;
 
  if (index)
    index->push_back(nd->getId());
 
  return underelements_;
}
 
/******************************************************************************/
 
RowMatrix<int> BipartitionList::toMatrix() const
{
  vector< map<string, bool>> bipl;
  for (size_t i = 0; i < getNumberOfBipartitions(); i++)
  {
    bipl.push_back(getBipartition(i));
  }
 
  vector<string> el = getElementNames();
 
  RowMatrix<int> mat(el.size(), getNumberOfBipartitions());
 
  for (size_t j = 0; j < el.size(); j++)
  {
    for (size_t i = 0; i < getNumberOfBipartitions(); i++)
    {
      mat(j, i) = bipl[i][el[j]];
    }
  }
  return mat;
}
 
/******************************************************************************/