SubstitutionRegister.h

Go to the documentation of this file.

// SPDX-FileCopyrightText: The Bio++ Development Group
//
// SPDX-License-Identifier: CECILL-2.1
 
#ifndef BPP_PHYL_MAPPING_SUBSTITUTIONREGISTER_H
#define BPP_PHYL_MAPPING_SUBSTITUTIONREGISTER_H
 
 
#include "../Model/StateMap.h"
 
// From bpp-core:
#include <Bpp/Clonable.h>
#include <Bpp/Numeric/Matrix/Matrix.h>
#include <Bpp/Text/StringTokenizer.h>
 
// From bpp-seq:
#include <Bpp/Seq/Alphabet/NucleicAlphabet.h>
#include <Bpp/Seq/Alphabet/ProteicAlphabet.h>
#include <Bpp/Seq/Alphabet/CodonAlphabet.h>
#include <Bpp/Seq/GeneticCode/GeneticCode.h>
 
// From the STL:
#include <vector>
#include <string>
#include <algorithm>
 
namespace bpp
{
class SubstitutionRegisterInterface :
  public virtual Clonable
{
public:
  SubstitutionRegisterInterface() {}
  virtual ~SubstitutionRegisterInterface() {}
 
  virtual SubstitutionRegisterInterface* clone() const = 0;
 
public:
  virtual const Alphabet& alphabet() const = 0;
 
  virtual std::shared_ptr<const Alphabet> getAlphabet() const = 0;
 
  virtual const StateMapInterface& stateMap() const = 0;
 
  virtual std::shared_ptr<const StateMapInterface> getStateMap() const = 0;
 
  virtual size_t getNumberOfSubstitutionTypes() const = 0;
 
  virtual const std::string& getName() const = 0;
 
  virtual size_t getType(size_t fromState, size_t toState) const = 0;
 
  virtual std::string getTypeName(size_t type) const = 0;
};
 
class AbstractSubstitutionRegister :
  public virtual SubstitutionRegisterInterface
{
protected:
  std::shared_ptr<const StateMapInterface> stateMap_;
  std::string name_;
 
public:
  AbstractSubstitutionRegister(std::shared_ptr<const StateMapInterface> stateMap, const std::string& name) :
    stateMap_(stateMap), name_(name)
  {}
 
  AbstractSubstitutionRegister(const AbstractSubstitutionRegister& asr) :
    stateMap_(asr.stateMap_), name_(asr.name_)
  {}
 
  AbstractSubstitutionRegister& operator=(const AbstractSubstitutionRegister& asr)
  {
    stateMap_ = asr.stateMap_;
    name_ = asr.name_;
    return *this;
  }
 
  virtual ~AbstractSubstitutionRegister() {}
 
public:
  const StateMapInterface& stateMap() const override { return *stateMap_; }
 
  std::shared_ptr<const StateMapInterface> getStateMap() const override { return stateMap_; }
 
  const Alphabet& alphabet() const override { return *stateMap_->getAlphabet(); }
 
  std::shared_ptr<const Alphabet> getAlphabet() const override { return stateMap_->getAlphabet(); }
 
  const std::string& getName() const override
  {
    return name_;
  }
};
 
class TotalSubstitutionRegister :
  public AbstractSubstitutionRegister
{
public:
  TotalSubstitutionRegister(std::shared_ptr<const StateMapInterface> stateMap) :
    AbstractSubstitutionRegister(stateMap, "Total")
  {}
 
  TotalSubstitutionRegister* clone() const override { return new TotalSubstitutionRegister(*this); }
 
public:
  size_t getNumberOfSubstitutionTypes() const override { return 1; }
 
  size_t getType(size_t fromState, size_t toState) const override
  {
    return fromState == toState ? 0 : 1;
  }
 
  std::string getTypeName(size_t type) const override
  {
    if (type == 0)
    {
      return "nosub";
    }
    else if (type == 1)
    {
      return "substitution";
    }
    else
    {
      throw Exception("TotalSubstitutionRegister::getTypeName. Bad substitution type.");
    }
  }
};
 
class CompleteSubstitutionRegister :
  public AbstractSubstitutionRegister
{
private:
  std::shared_ptr<const SubstitutionRegisterInterface> preg_;
 
  bool isRegComplete_;
 
public:
  CompleteSubstitutionRegister(std::shared_ptr<const SubstitutionRegisterInterface> reg) :
    AbstractSubstitutionRegister(reg->getStateMap(), reg->getName()),
    preg_(reg),
    isRegComplete_(true)
  {
    size_t size = reg->alphabet().getSize();
    for (size_t i = 0; i < size; ++i)
    {
      for (size_t j = 0; j < size; ++j)
      {
        if ((i != j) && reg->getType(i, j) == 0)
        {
          isRegComplete_ = false;
          return;
        }
      }
    }
  }
 
  CompleteSubstitutionRegister* clone() const override { return new CompleteSubstitutionRegister(*this); }
 
  CompleteSubstitutionRegister(const CompleteSubstitutionRegister& csr) :
    AbstractSubstitutionRegister(csr),
    preg_(csr.preg_),
    isRegComplete_(csr.isRegComplete_)
  {}
 
  CompleteSubstitutionRegister& operator=(const CompleteSubstitutionRegister& csr)
  {
    AbstractSubstitutionRegister::operator=(csr);
    preg_ = csr.preg_;
    isRegComplete_ = csr.isRegComplete_;
    return *this;
  }
 
  virtual ~CompleteSubstitutionRegister() {}
 
public:
  size_t getNumberOfSubstitutionTypes() const override
  {
    return preg_->getNumberOfSubstitutionTypes() + (isRegComplete_ ? 0 : 1);
  }
 
  size_t getType(size_t fromState, size_t toState) const override
  {
    if (fromState == toState)
      return 0;
 
    size_t t = preg_->getType(fromState, toState);
    if (t == 0)
      return getNumberOfSubstitutionTypes();
    else
      return t;
  }
 
  std::string getTypeName(size_t type) const override
  {
    try
    {
      return preg_->getTypeName(type);
    }
    catch (Exception& e)
    {
      if (type == getNumberOfSubstitutionTypes())
        return "Completion_substitution";
      else
        throw Exception("CompleteSubstitutionRegister::getTypeName. Bad substitution type.");
    }
  }
};
 
 
class VectorOfSubstitutionRegisters :
  public AbstractSubstitutionRegister
{
private:
  std::vector< std::shared_ptr<SubstitutionRegisterInterface>> vSubReg_;
 
public:
  VectorOfSubstitutionRegisters(std::shared_ptr<const StateMapInterface> stateMap) :
    AbstractSubstitutionRegister(stateMap, "Combination"),
    vSubReg_()
  {}
 
  VectorOfSubstitutionRegisters(const VectorOfSubstitutionRegisters& vosr) :
    AbstractSubstitutionRegister(vosr),
    vSubReg_(vosr.vSubReg_)
  {}
 
  VectorOfSubstitutionRegisters* clone() const override
  {
    return new VectorOfSubstitutionRegisters(*this);
  }
 
  virtual ~VectorOfSubstitutionRegisters() {}
 
  void addRegister(std::shared_ptr<SubstitutionRegisterInterface> reg)
  {
    if (reg)
    {
      if (reg->stateMap() != stateMap())
        throw Exception("VectorOfSubstitutionRegisters::addRegister : mismatch between state maps");
 
      vSubReg_.push_back(reg);
    }
  }
 
  size_t getType(size_t i, size_t j) const override
  {
    if (i == j)
      return 0;
 
    size_t x = 0;
 
    for (size_t p = 0; p < vSubReg_.size(); ++p)
    {
      x *= vSubReg_[p]->getNumberOfSubstitutionTypes();
      size_t z = vSubReg_[p]->getType(i, j);
      if (z == 0)
        return 0;
      x += z - 1;
    }
 
    return x + 1;
  }
 
  size_t getNumberOfSubstitutionTypes() const override
  {
    if (vSubReg_.size() == 0)
      return 0;
 
    size_t n = 1;
 
    for (size_t i = 0; i < vSubReg_.size(); ++i)
    {
      n *= vSubReg_[i]->getNumberOfSubstitutionTypes();
    }
 
    return n;
  }
 
  std::string getTypeName(size_t type) const override
  {
    if (type == 0)
      return "nosub";
 
    std::string res = "";
    size_t ty = type - 1;
 
    for (size_t p = vSubReg_.size(); p > 0; p--)
    {
      if (p != vSubReg_.size())
        res = "_X_" + res;
      res = vSubReg_[p - 1]->getTypeName((ty % vSubReg_[p - 1]->getNumberOfSubstitutionTypes()) + 1) + res;
 
      ty /= vSubReg_[p - 1]->getNumberOfSubstitutionTypes();
    }
 
    return res;
  }
};
 
 
class GeneralSubstitutionRegister :
  public AbstractSubstitutionRegister
{
protected:
  size_t size_;
 
  RowMatrix<size_t> matrix_;
 
  std::map<size_t, std::map<size_t, std::vector<size_t>>> types_;
 
public:
  GeneralSubstitutionRegister(std::shared_ptr<const StateMapInterface> stateMap) :
    AbstractSubstitutionRegister(stateMap, "General"),
    size_(stateMap->getNumberOfModelStates()),
    matrix_(size_, size_),
    types_()
  {}
 
  GeneralSubstitutionRegister(std::shared_ptr<const StateMapInterface> stateMap, const RowMatrix<size_t>& matrix) :
    AbstractSubstitutionRegister(stateMap, "General"),
    size_(stateMap->getNumberOfModelStates()),
    matrix_(matrix),
    types_()
  {
    if (matrix_.getNumberOfRows() != size_)
      throw DimensionException("GeneralSubstitutionRegister", size_, matrix_.getNumberOfRows());
    if (matrix_.getNumberOfColumns() != size_)
      throw DimensionException("GeneralSubstitutionRegister", size_, matrix_.getNumberOfColumns());
    updateTypes_();
  }
 
  GeneralSubstitutionRegister(const GeneralSubstitutionRegister& gsr) :
    AbstractSubstitutionRegister(gsr),
    size_(gsr.size_),
    matrix_(gsr.matrix_),
    types_(gsr.types_)
  {}
 
  GeneralSubstitutionRegister& operator=(const GeneralSubstitutionRegister& gsr)
  {
    AbstractSubstitutionRegister::operator=(gsr);
    size_   = gsr.size_;
    matrix_ = gsr.matrix_;
    types_  = gsr.types_;
    return *this;
  }
 
  GeneralSubstitutionRegister* clone() const override
  {
    return new GeneralSubstitutionRegister(*this);
  }
 
  virtual ~GeneralSubstitutionRegister() {}
 
  size_t getType(size_t i, size_t j) const override
  {
    return matrix_(i, j);
  }
 
  size_t getNumberOfSubstitutionTypes() const override
  {
    return types_.find(0) == types_.end() ? types_.size() : types_.size() - 1;
  }
 
  std::string getTypeName(size_t type) const override
  {
    if (types_.find(type) != types_.end())
      return TextTools::toString(type);
 
    throw Exception("Bad type number " + TextTools::toString(type) + " in GeneralSubstitutionRegister::getTypeName.");
  }
 
protected:
  void updateTypes_();
};
 
 
class SelectedSubstitutionRegister :
  public GeneralSubstitutionRegister
{
  std::map<size_t, std::string> categoryNames_;
 
public:
  SelectedSubstitutionRegister (std::shared_ptr<const StateMapInterface> stateMap, std::string selectedSubstitutions) :
    GeneralSubstitutionRegister(stateMap),
    categoryNames_()
  {
    selectedSubstitutions.erase(std::remove(selectedSubstitutions.begin(), selectedSubstitutions.end(), ' '), selectedSubstitutions.end());
    size_t typeSubs = 0;
    size_t coord1 = 0;
    size_t coord2 = 0;
    std::string codon1 = "";
    std::string codon2 = "";
    StringTokenizer subsGroup(selectedSubstitutions, "()");
    while (subsGroup.hasMoreToken())
    {
      typeSubs++;
      StringTokenizer namesSubs(subsGroup.nextToken(), ":");
      if (namesSubs.numberOfRemainingTokens() == 2)
      {
        categoryNames_[typeSubs] = namesSubs.nextToken();
      }
      else if (namesSubs.numberOfRemainingTokens() == 1)
      {
        categoryNames_[typeSubs] = TextTools::toString(typeSubs);
      }
      StringTokenizer substitutions(namesSubs.nextToken(), ",");
      while (substitutions.hasMoreToken())
      {
        StringTokenizer coordinates(substitutions.nextToken(), "->");
        codon1 = coordinates.nextToken();
        codon2 = coordinates.nextToken();
        coord1 = stateMap->alphabet().getStateIndex(codon1);
        coord2 = stateMap->alphabet().getStateIndex(codon2);
        this->matrix_(coord1, coord2) = typeSubs;
      }
    }
    updateTypes_();
  }
 
  SelectedSubstitutionRegister* clone() const { return new SelectedSubstitutionRegister(*this); }
 
  ~SelectedSubstitutionRegister() {}
 
  std::string getTypeName(size_t type) const
  {
    if (types_.find(type) != types_.end())
      return TextTools::toString(categoryNames_.find(type)->second);
 
    throw Exception("Bad type number " + TextTools::toString(type) + " in GeneralSubstitutionRegister::getTypeName.");
  }
};
 
class AAInteriorSubstitutionRegister :
  public GeneralSubstitutionRegister
{
  std::map<std::string, size_t> categoryCorrespondance_;
  std::shared_ptr<const GeneticCode> genCode_;
 
public:
  AAInteriorSubstitutionRegister(
      std::shared_ptr<const StateMapInterface> stateMap,
      std::shared_ptr<const GeneticCode> genCode) :
    GeneralSubstitutionRegister(stateMap),
    categoryCorrespondance_(),
    genCode_(genCode)
  {
    updateMatrix_();
    updateTypes_();
  }
 
  AAInteriorSubstitutionRegister(const AAInteriorSubstitutionRegister& ais) :
    GeneralSubstitutionRegister(ais),
    categoryCorrespondance_(ais.categoryCorrespondance_),
    genCode_(ais.genCode_)
  {}
 
  AAInteriorSubstitutionRegister& operator=(const AAInteriorSubstitutionRegister& ais)
  {
    GeneralSubstitutionRegister::operator=(ais);
    genCode_ = ais.genCode_;
    categoryCorrespondance_ = ais.categoryCorrespondance_;
 
    return *this;
  }
 
  AAInteriorSubstitutionRegister* clone() const { return new AAInteriorSubstitutionRegister(*this); }
 
  ~AAInteriorSubstitutionRegister() {}
 
 
  std::string getTypeName(size_t type) const
  {
    if (types_.find(type) != types_.end())
    {
      for (std::map<std::string, size_t>::const_iterator it = categoryCorrespondance_.begin(); it != categoryCorrespondance_.end(); it++)
      {
        if (it->second == type)
          return TextTools::toString(it->first);
      }
    }
    throw Exception("Bad type number " + TextTools::toString(type) + " in GeneralSubstitutionRegister::getTypeName.");
  }
 
protected:
  void updateMatrix_()
  {
    size_t categoryIndex = 1;
    for (size_t i = 1; i <= stateMap().getNumberOfModelStates(); ++i)
    {
      int state1 = stateMap().alphabet().getStateAt(i).getNum();
      for (size_t j = i + 1; j <= stateMap().getNumberOfModelStates(); ++j)
      {
        int state2 = stateMap().alphabet().getStateAt(j).getNum();
        if (!(genCode_->isStop(state1)) && !(genCode_->isStop(state2)))
        {
          if (genCode_->translate(state1) == genCode_->translate(state2))
          {
            std::string aminoAcid = genCode_->getTargetAlphabet()->intToChar(genCode_->translate(state1));
            if (categoryCorrespondance_.find(aminoAcid) == categoryCorrespondance_.end())
            {
              categoryCorrespondance_[aminoAcid] = categoryIndex;
              categoryIndex++;
            }
            matrix_(i, j) = categoryCorrespondance_[aminoAcid];
            matrix_(j, i) = categoryCorrespondance_[aminoAcid];
          }
        }
      }
    }
  }
};
 
class AAExteriorSubstitutionRegister :
  public GeneralSubstitutionRegister
{
  std::map<std::string, size_t> categoryCorrespondance_;
  std::shared_ptr<const GeneticCode> genCode_;
 
public:
  AAExteriorSubstitutionRegister(
      std::shared_ptr<const StateMapInterface> stateMap,
      std::shared_ptr<const GeneticCode> genCode) :
    GeneralSubstitutionRegister(stateMap),
    categoryCorrespondance_(),
    genCode_(genCode)
  {
    updateMatrix_();
    updateTypes_();
  }
 
  AAExteriorSubstitutionRegister(const AAExteriorSubstitutionRegister& ais) :
    GeneralSubstitutionRegister(ais),
    categoryCorrespondance_(ais.categoryCorrespondance_),
    genCode_(ais.genCode_)
  {}
 
  AAExteriorSubstitutionRegister& operator=(const AAExteriorSubstitutionRegister& ais)
  {
    GeneralSubstitutionRegister::operator=(ais);
    categoryCorrespondance_ = ais.categoryCorrespondance_;
    genCode_ = ais.genCode_;
 
    return *this;
  }
 
  AAExteriorSubstitutionRegister* clone() const { return new AAExteriorSubstitutionRegister(*this); }
 
  ~AAExteriorSubstitutionRegister() {}
 
  std::string getTypeName(size_t type) const
  {
    if (types_.find(type) != types_.end())
    {
      for (std::map<std::string, size_t>::const_iterator it = categoryCorrespondance_.begin(); it != categoryCorrespondance_.end(); it++)
      {
        if (it->second == type)
          return TextTools::toString(it->first);
      }
    }
    throw Exception("Bad type number " + TextTools::toString(type) + " in GeneralSubstitutionRegister::getTypeName.");
  }
 
protected:
  void updateMatrix_()
  {
    size_t categoryIndex = 1;
    for (size_t i = 1; i <= stateMap().getNumberOfModelStates(); ++i)
    {
      int state1 = stateMap().alphabet().getStateAt(i).getNum();
      for (size_t j = i + 1; j <= stateMap().getNumberOfModelStates(); ++j)
      {
        int state2 = stateMap().alphabet().getStateAt(j).getNum();
        if (!(genCode_->isStop(state1)) && !(genCode_->isStop(state2)))
        {
          if (genCode_->translate(state1) != genCode_->translate(state2))
          {
            std::string aminoAcid1 = genCode_->getTargetAlphabet()->intToChar(genCode_->translate(state1));
            std::string aminoAcid2 = genCode_->getTargetAlphabet()->intToChar(genCode_->translate(state2));
            bool AA1IsNotInGroup = ((categoryCorrespondance_.find(aminoAcid1 + "->" + aminoAcid2) == categoryCorrespondance_.end()));
            bool AA2IsNotInGroup = ((categoryCorrespondance_.find(aminoAcid2 + "->" + aminoAcid1) == categoryCorrespondance_.end()));
            if (AA1IsNotInGroup)
            {
              categoryCorrespondance_[aminoAcid1 + "->" + aminoAcid2] = categoryIndex;
              categoryIndex++;
            }
            this->matrix_(i, j) = categoryCorrespondance_[aminoAcid1 + "->" + aminoAcid2];
            if (AA2IsNotInGroup)
            {
              categoryCorrespondance_[aminoAcid2 + "->" + aminoAcid1] = categoryIndex;
              categoryIndex++;
            }
            this->matrix_(j, i) = categoryCorrespondance_[aminoAcid2 + "->" + aminoAcid1];
          }
        }
      }
    }
  }
};
 
class TsTvSubstitutionRegister :
  public AbstractSubstitutionRegister
{
private:
  std::shared_ptr<const GeneticCode> genCode_;
 
public:
  TsTvSubstitutionRegister(std::shared_ptr<const StateMapInterface> stateMap) :
    AbstractSubstitutionRegister(stateMap, "TsTv"),
    genCode_()
  {}
 
  TsTvSubstitutionRegister(std::shared_ptr<const StateMapInterface> stateMap, std::shared_ptr<const GeneticCode> genCode) :
    AbstractSubstitutionRegister(stateMap, "TsTv"),
    genCode_(genCode)
  {}
 
  TsTvSubstitutionRegister(const TsTvSubstitutionRegister& reg) :
    AbstractSubstitutionRegister(reg),
    genCode_(reg.genCode_)
  {}
 
  TsTvSubstitutionRegister& operator=(const TsTvSubstitutionRegister& reg)
  {
    AbstractSubstitutionRegister::operator=(reg);
    genCode_ = reg.genCode_;
    return *this;
  }
 
  TsTvSubstitutionRegister* clone() const { return new TsTvSubstitutionRegister(*this); }
 
public:
  size_t getNumberOfSubstitutionTypes() const { return 2; }
 
  size_t getType(size_t fromState, size_t toState) const
  {
    int x = stateMap().getAlphabetStateAsInt(fromState);
    int y = stateMap().getAlphabetStateAsInt(toState);
    if (x == y)
      return 0; // nothing happens
 
 
    if (genCode_)
    {
      auto cAlpha = genCode_->getCodonAlphabet();
      if (cAlpha->isGap(x)
          || cAlpha->isGap(y)
          || genCode_->isStop(x)
          || genCode_->isStop(y))
        return 0;
 
      size_t countPos = 0;
      if (cAlpha->getFirstPosition(x) != cAlpha->getFirstPosition(y))
        countPos++;
      if (cAlpha->getSecondPosition(x) != cAlpha->getSecondPosition(y))
        countPos++;
      if (cAlpha->getThirdPosition(x) != cAlpha->getThirdPosition(y))
        countPos++;
      if (countPos > 1)
        return 0;
    }
 
    int d = std::abs(x - y);
 
    while (d != 0)
    {
      if (d == 2)
        return 1; // This is a transition
      else
        d /= 4;
    }
 
    return 2; // This is a transversion
  }
 
  std::string getTypeName(size_t type) const
  {
    if (type == 0)
    {
      return "nosub";
    }
    else if (type == 1)
    {
      return "Ts";
    }
    else if (type == 2)
    {
      return "Tv";
    }
    else
    {
      throw Exception("TsTvSubstitutionRegister::getTypeName. Bad substitution type.");
    }
  }
};
 
class SWSubstitutionRegister :
  public AbstractSubstitutionRegister
{
private:
  std::shared_ptr<const GeneticCode> genCode_;
 
public:
  SWSubstitutionRegister(std::shared_ptr<const StateMapInterface> stateMap) :
    AbstractSubstitutionRegister(stateMap, "SW"),
    genCode_()
  {}
 
  SWSubstitutionRegister(
      std::shared_ptr<const StateMapInterface> stateMap,
      std::shared_ptr<const GeneticCode> genCode) :
    AbstractSubstitutionRegister(stateMap, "SW"),
    genCode_(genCode)
  {}
 
  SWSubstitutionRegister(const SWSubstitutionRegister& reg) :
    AbstractSubstitutionRegister(reg),
    genCode_(reg.genCode_)
  {}
 
  SWSubstitutionRegister& operator=(const SWSubstitutionRegister& reg)
  {
    AbstractSubstitutionRegister::operator=(reg);
    genCode_ = reg.genCode_;
    return *this;
  }
 
  SWSubstitutionRegister* clone() const { return new SWSubstitutionRegister(*this); }
 
public:
  size_t getNumberOfSubstitutionTypes() const { return 4; }
 
  size_t getType(size_t fromState, size_t toState) const
  {
    int x = stateMap().getAlphabetStateAsInt(fromState);
    int y = stateMap().getAlphabetStateAsInt(toState);
    if (x == y)
      return 0; // nothing happens
 
    int nd, na;
    if (genCode_)
    {
      auto cAlpha = genCode_->getCodonAlphabet();
      if (cAlpha->isGap(x)
          || cAlpha->isGap(y)
          || genCode_->isStop(x)
          || genCode_->isStop(y))
        return 0;
 
      nd = cAlpha->getFirstPosition(x);
      na = cAlpha->getFirstPosition(y);
      if (na != nd)
      {
        if (cAlpha->getSecondPosition(x) != cAlpha->getSecondPosition(y)
            || (cAlpha->getThirdPosition(x) != cAlpha->getThirdPosition(y)))
          return 0;
      }
      else
      {
        nd = cAlpha->getSecondPosition(x);
        na = cAlpha->getSecondPosition(y);
        if (na != nd)
        {
          if (cAlpha->getThirdPosition(x) != cAlpha->getThirdPosition(y))
            return 0;
        }
        else
        {
          nd = cAlpha->getThirdPosition(x);
          na = cAlpha->getThirdPosition(y);
        }
      }
    }
    else
    {
      nd = x;
      na = y;
    }
 
    switch (nd)
    {
    case 0:
      return (na == 3) ? 4 : 3;
    case 1:
      return (na == 2) ? 1 : 2;
    case 2:
      return (na == 1) ? 1 : 2;
    case 3:
      return (na == 0) ? 4 : 3;
    default:
      return 0;
    }
  }
 
  std::string getTypeName(size_t type) const
  {
    switch (type)
    {
    case 0:
      return "nosub";
    case 1:
      return "S->S";
    case 2:
      return "S->W";
    case 3:
      return "W->S";
    case 4:
      return "W->W";
    default:
      throw Exception("SWSubstitutionRegister::getTypeName. Bad substitution type.");
    }
  }
};
 
 
class DnDsSubstitutionRegister :
  public AbstractSubstitutionRegister
{
private:
  std::shared_ptr<const GeneticCode> genCode_;
  bool countMultiple_;
 
public:
  DnDsSubstitutionRegister(
      std::shared_ptr<const StateMapInterface> stateMap,
      std::shared_ptr<const GeneticCode> genCode,
      bool countMultiple = false) :
    AbstractSubstitutionRegister(stateMap, "DnDs"),
    genCode_(genCode),
    countMultiple_(countMultiple)
  {}
 
  DnDsSubstitutionRegister(const DnDsSubstitutionRegister& reg) :
    AbstractSubstitutionRegister(reg),
    genCode_(reg.genCode_),
    countMultiple_(reg.countMultiple_)
  {}
 
  DnDsSubstitutionRegister& operator=(const DnDsSubstitutionRegister& reg)
  {
    AbstractSubstitutionRegister::operator=(reg);
    genCode_ = reg.genCode_;
    countMultiple_ = reg.countMultiple_;
    return *this;
  }
 
  DnDsSubstitutionRegister* clone() const { return new DnDsSubstitutionRegister(*this); }
 
public:
  size_t getNumberOfSubstitutionTypes() const { return 2; }
 
  size_t getType(size_t fromState, size_t toState) const
  {
    int x = stateMap().getAlphabetStateAsInt(fromState);
    int y = stateMap().getAlphabetStateAsInt(toState);
 
    auto cAlpha = genCode_->getCodonAlphabet();
    if (cAlpha->isGap(x)
        || cAlpha->isGap(y)
        || genCode_->isStop(x)
        || genCode_->isStop(y))
      return 0;
    if (x == y)
      return 0; // nothing happens
    if (!countMultiple_)
    {
      size_t countPos = 0;
      if (cAlpha->getFirstPosition(x) != cAlpha->getFirstPosition(y))
        countPos++;
      if (cAlpha->getSecondPosition(x) != cAlpha->getSecondPosition(y))
        countPos++;
      if (cAlpha->getThirdPosition(x) != cAlpha->getThirdPosition(y))
        countPos++;
      if (countPos > 1)
        return 0;
    }
    return genCode_->areSynonymous(x, y) ? 1 : 2;
  }
 
  std::string getTypeName (size_t type) const
  {
    if (type == 0)
    {
      return "nosub";
    }
    else if (type == 1)
    {
      return "dS";
    }
    else if (type == 2)
    {
      return "dN";
    }
    else
    {
      throw Exception("DnDsSubstitutionRegister::getTypeName. Bad substitution type.");
    }
  }
};
 
class KrKcSubstitutionRegister :
  public AbstractSubstitutionRegister
{
private:
  std::vector< std::vector<bool>> types_;
 
  // In case of codon model
  std::shared_ptr<const GeneticCode> genCode_;
 
public:
  KrKcSubstitutionRegister(std::shared_ptr<const StateMapInterface> stateMap) :
    AbstractSubstitutionRegister(stateMap, "KrKc"),
    types_(20),
    genCode_()
  {
    init();
  }
 
  KrKcSubstitutionRegister(
      std::shared_ptr<const StateMapInterface> stateMap,
      std::shared_ptr<const GeneticCode> genCode) :
    AbstractSubstitutionRegister(stateMap, "KrKc"),
    types_(20),
    genCode_(genCode)
  {
    init();
  }
 
  KrKcSubstitutionRegister(const KrKcSubstitutionRegister& kreg) :
    AbstractSubstitutionRegister(kreg),
    types_(kreg.types_),
    genCode_(kreg.genCode_)
  {}
 
  KrKcSubstitutionRegister& operator=(const KrKcSubstitutionRegister& kreg)
  {
    AbstractSubstitutionRegister::operator=(kreg);
    types_ = kreg.types_;
    genCode_ = kreg.genCode_;
 
    return *this;
  }
 
  void init()
  {
    for (size_t i = 0; i < 20; ++i)
    {
      types_[i].resize(20);
      for (size_t j = 0; j < 20; ++j)
      {
        types_[i][j] = false;
      }
    }
 
    types_[0][7] = true;
    types_[0][14] = true;
    types_[0][19] = true;
    types_[7][0] = true;
    types_[7][14] = true;
    types_[7][19] = true;
    types_[14][0] = true;
    types_[14][7] = true;
    types_[14][19] = true;
    types_[19][0] = true;
    types_[19][7] = true;
    types_[19][14] = true;
 
    types_[1][5] = true;
    types_[1][6] = true;
    types_[1][8] = true;
    types_[1][11] = true;
    types_[1][17] = true;
    types_[1][18] = true;
    types_[5][6] = true;
    types_[5][8] = true;
    types_[5][11] = true;
    types_[5][1] = true;
    types_[5][17] = true;
    types_[5][18] = true;
    types_[6][8] = true;
    types_[6][11] = true;
    types_[6][5] = true;
    types_[6][1] = true;
    types_[6][17] = true;
    types_[6][18] = true;
    types_[8][6] = true;
    types_[8][11] = true;
    types_[8][5] = true;
    types_[8][1] = true;
    types_[8][17] = true;
    types_[8][18] = true;
    types_[11][1] = true;
    types_[11][5] = true;
    types_[11][6] = true;
    types_[11][8] = true;
    types_[11][17] = true;
    types_[11][18] = true;
    types_[17][1] = true;
    types_[17][5] = true;
    types_[17][6] = true;
    types_[17][8] = true;
    types_[17][11] = true;
    types_[17][18] = true;
    types_[18][1] = true;
    types_[18][6] = true;
    types_[18][8] = true;
    types_[18][11] = true;
    types_[18][5] = true;
    types_[18][17] = true;
 
    types_[2][3] = true;
    types_[2][4] = true;
    types_[2][15] = true;
    types_[2][16] = true;
    types_[3][2] = true;
    types_[3][4] = true;
    types_[3][15] = true;
    types_[3][16] = true;
    types_[4][2] = true;
    types_[4][3] = true;
    types_[4][15] = true;
    types_[4][16] = true;
    types_[15][2] = true;
    types_[15][3] = true;
    types_[15][4] = true;
    types_[15][16] = true;
    types_[16][2] = true;
    types_[16][3] = true;
    types_[16][4] = true;
    types_[16][15] = true;
 
    types_[9][10] = true;
    types_[9][12] = true;
    types_[9][13] = true;
    types_[10][9] = true;
    types_[10][12] = true;
    types_[10][13] = true;
    types_[12][9] = true;
    types_[12][10] = true;
    types_[12][13] = true;
    types_[13][9] = true;
    types_[13][10] = true;
    types_[13][12] = true;
  }
 
  KrKcSubstitutionRegister* clone() const { return new KrKcSubstitutionRegister(*this); }
 
public:
  size_t getNumberOfSubstitutionTypes() const { return 2; }
 
  size_t getType(size_t fromState, size_t toState) const
  {
    if (fromState == toState)
      return 0; // nothing happens
    int x = stateMap().getAlphabetStateAsInt(fromState);
    int y = stateMap().getAlphabetStateAsInt(toState);
 
    if (genCode_) // Codon model
    {
      if (genCode_->getSourceAlphabet()->isGap(x)
          || genCode_->getSourceAlphabet()->isGap(y)
          || genCode_->isStop(x)
          || genCode_->isStop(y))
        return 0;
      if (genCode_->areSynonymous(x, y))
        return 0;
 
      // avoid multiple substitutions
      auto cAlpha = genCode_->getCodonAlphabet();
 
      size_t countPos = 0;
      if (cAlpha->getFirstPosition(x) != cAlpha->getFirstPosition(y))
        countPos++;
      if (cAlpha->getSecondPosition(x) != cAlpha->getSecondPosition(y))
        countPos++;
      if (cAlpha->getThirdPosition(x) != cAlpha->getThirdPosition(y))
        countPos++;
      if (countPos > 1)
        return 0;
 
      x = genCode_->translate(x);
      y = genCode_->translate(y);
    }
 
    return types_[static_cast<size_t>(x)][static_cast<size_t>(y)] ? 1 : 2;
  }
 
  std::string getTypeName(size_t type) const
  {
    if (type == 0)
    {
      return "nosub";
    }
    else if (type == 1)
    {
      return "Kc";
    }
    else if (type == 2)
    {
      return "Kr";
    }
    else
    {
      throw Exception("KrKcSubstitutionRegister::getTypeName. Bad substitution type.");
    }
  }
};
} // end of namespace bpp.
#endif // BPP_PHYL_MAPPING_SUBSTITUTIONREGISTER_H