BppORateDistributionFormat.cpp

Go to the documentation of this file.

// SPDX-FileCopyrightText: The Bio++ Development Group
//
// SPDX-License-Identifier: CECILL-2.1
 
#include "../Model/RateDistribution/ConstantRateDistribution.h"
#include "../Model/RateDistribution/ExponentialDiscreteRateDistribution.h"
#include "../Model/RateDistribution/GammaDiscreteRateDistribution.h"
#include "../Model/RateDistribution/GaussianDiscreteRateDistribution.h"
#include "BppORateDistributionFormat.h"
 
// From bpp-core:
#include <Bpp/Text/KeyvalTools.h>
#include <Bpp/Numeric/Prob/InvariantMixedDiscreteDistribution.h>
#include <Bpp/Numeric/Prob/SimpleDiscreteDistribution.h>
#include <Bpp/Numeric/Prob/MixtureOfDiscreteDistributions.h>
#include <Bpp/Io/BppOParametrizableFormat.h>
 
using namespace bpp;
 
// From the STL:
#include <iomanip>
 
using namespace std;
 
 
unique_ptr<DiscreteDistributionInterface> BppORateDistributionFormat::readDiscreteDistribution(
    const std::string& distDescription,
    bool parseArguments)
{
  unparsedArguments_.clear();
  string distName;
  map<string, string> args;
  KeyvalTools::parseProcedure(distDescription, distName, args);
  unique_ptr<DiscreteDistributionInterface> rDist;
 
  if (distName == "Uniform")
    throw Exception("BppO Warning: Uniform distribution is deprecated, use Constant instead.");
 
  if ((distName == "InvariantMixed") || (distName == "Invariant"))
  {
    // We have to parse the nested distribution first:
    string nestedDistDescription = args["dist"];
    if (TextTools::isEmpty(nestedDistDescription))
      throw Exception("BppORateDistributionFormat::read. Missing argument 'dist' for distribution 'Invariant'.");
    if (verbose_)
      ApplicationTools::displayResult("Invariant Mixed distribution", distName );
    BppORateDistributionFormat nestedReader(false);
    auto nestedDistribution = nestedReader.readDiscreteDistribution(nestedDistDescription, false);
    map<string, string> unparsedArgumentsNested(nestedReader.getUnparsedArguments());
 
    // Now we create the Invariant rate distribution:
    rDist = make_unique<InvariantMixedDiscreteDistribution>(std::move(nestedDistribution), 0.1, 0.000001);
 
    // Then we update the parameter set:
    for (auto it : unparsedArgumentsNested)
    {
      unparsedArguments_["Invariant." + it.first] = it.second;
    }
 
    if (args.find("p") != args.end())
      unparsedArguments_["Invariant.p"] = args["p"];
  }
  else if (distName == "Constant")
  {
    if (!allowConstant_)
      throw Exception("BppORateDistributionFormat::read(). Constant distribution not allowed.");
 
    if (args.find("value") != args.end())
      ApplicationTools::displayMessage("Found argument 'value' in Constant distribution. Constant distribution is defined to have an average of 1.");
    rDist = make_unique<ConstantRateDistribution>();
  }
  else if (distName == "Simple")
  {
    if (args.find("values") == args.end())
      throw Exception("Missing argument 'values' in Simple distribution");
    if (args.find("probas") == args.end())
      throw Exception("Missing argument 'probas' in Simple distribution");
    vector<double> probas, values;
 
    string rf = args["values"];
    StringTokenizer strtok(rf.substr(1, rf.length() - 2), ",");
    while (strtok.hasMoreToken())
      values.push_back(TextTools::toDouble(strtok.nextToken()));
 
    rf = args["probas"];
    StringTokenizer strtok2(rf.substr(1, rf.length() - 2), ",");
    while (strtok2.hasMoreToken())
      probas.push_back(TextTools::toDouble(strtok2.nextToken()));
 
    std::map<size_t, std::vector<double>> ranges;
 
    if (args.find("ranges") != args.end())
    {
      string rr = args["ranges"];
      StringTokenizer strtok3(rr.substr(1, rr.length() - 2), ",");
      string desc;
      double deb, fin;
      size_t num;
      size_t po, pf, ppv;
      while (strtok3.hasMoreToken())
      {
        desc = strtok3.nextToken();
        po = desc.find("[");
        ppv = desc.find(";");
        pf = desc.find("]");
        num = (size_t)(TextTools::toInt(desc.substr(1, po - 1)));
        deb = TextTools::toDouble(desc.substr(po + 1, ppv - po - 1));
        fin = TextTools::toDouble(desc.substr(ppv + 1, pf - ppv - 1));
        vector<double> vd;
        vd.push_back(deb);
        vd.push_back(fin);
        ranges[num] = vd;
      }
    }
    if (ranges.size() == 0)
      rDist = make_unique<SimpleDiscreteDistribution>(values, probas);
    else
      rDist = make_unique<SimpleDiscreteDistribution>(values, ranges, probas);
 
    vector<string> v = rDist->getParameters().getParameterNames();
 
    for (size_t i = 0; i < v.size(); ++i)
    {
      unparsedArguments_[v[i]] = TextTools::toString(rDist->getParameterValue(rDist->getParameterNameWithoutNamespace(v[i])));
    }
  }
  else if (distName == "Mixture")
  {
    if (args.find("probas") == args.end())
      throw Exception("Missing argument 'probas' in Mixture distribution");
    vector<double> probas;
    vector<unique_ptr<DiscreteDistributionInterface>> v_pdd;
    string rf = args["probas"];
    StringTokenizer strtok2(rf.substr(1, rf.length() - 2), ",");
    while (strtok2.hasMoreToken())
      probas.push_back(TextTools::toDouble(strtok2.nextToken()));
 
    vector<string> v_nestedDistrDescr;
 
    size_t nbd = 0;
    while (args.find("dist" + TextTools::toString(++nbd)) != args.end())
      v_nestedDistrDescr.push_back(args["dist" + TextTools::toString(nbd)]);
 
    if (v_nestedDistrDescr.size() != probas.size())
      throw Exception("Number of distributions (keyword 'dist" + TextTools::toString(probas.size()) + "') do not fit the number of probabilities");
 
    for (unsigned i = 0; i < v_nestedDistrDescr.size(); i++)
    {
      BppORateDistributionFormat nestedReader(false);
      auto pdd = nestedReader.readDiscreteDistribution(v_nestedDistrDescr[i], false);
      map<string, string> unparsedArgumentsNested(nestedReader.getUnparsedArguments());
 
      for (auto& it : unparsedArgumentsNested)
      {
        unparsedArguments_[distName + "." + TextTools::toString(i + 1) + "_" + it.first] = it.second;
      }
      v_pdd.push_back(std::move(pdd));
    }
    rDist = make_unique<MixtureOfDiscreteDistributions>(v_pdd, probas);
  }
  else
  {
    if (args.find("n") == args.end())
      throw Exception("Missing argument 'n' (number of classes) in " + distName
            + " distribution");
    size_t nbClasses = TextTools::to<size_t>(args["n"]);
 
    if (distName == "Gamma")
    {
      rDist = make_unique<GammaDiscreteRateDistribution>(nbClasses, 1.);
 
      if (args.find("alpha") != args.end())
        unparsedArguments_["Gamma.alpha"] = args["alpha"];
      if (args.find("beta") != args.end())
        throw Exception("Found argument 'beta' in Gamma distribution. Gamma distribution is defined to have an average of 1, with beta=alpha.");
    }
    else if (distName == "Gaussian")
    {
      rDist = make_unique<GaussianDiscreteRateDistribution>(nbClasses, 1);
      if (args.find("mu") != args.end())
        throw Exception("Found argument 'mu' in Gaussian distribution. Gaussian distribution is defined to have an average of 1, with mu=1.");
      if (args.find("sigma") != args.end())
        unparsedArguments_["Gaussian.sigma"] = args["sigma"];
    }
    else if (distName == "Exponential")
    {
      rDist = make_unique<ExponentialDiscreteRateDistribution>(nbClasses);
      if (args.find("lambda") != args.end())
        throw Exception("Found argument 'lambda' in Exponential distribution. Exponential distribution is defined to have an average of 1, with lambda=1.");
    }
    else
    {
      throw Exception("Unsupported rate distribution: " + distName + ".");
    }
  }
  if (verbose_)
  {
    ApplicationTools::displayResult("Distribution", distName);
    ApplicationTools::displayResult("Number of classes", TextTools::toString((int)rDist->getNumberOfCategories()));
  }
 
  if (parseArguments)
    initialize_(*rDist);
 
  return rDist;
}
 
 
void BppORateDistributionFormat::writeDiscreteDistribution(
    const DiscreteDistributionInterface& dist,
    OutputStream& out,
    std::map<std::string, std::string>& globalAliases,
    std::vector<std::string>& writtenNames) const
{
  bool comma = false;
 
 
  out << dist.getName() + "(";
 
  auto* invar = dynamic_cast<const InvariantMixedDiscreteDistribution*>(&dist);
  if (invar)
  {
    auto& pd = invar->variableSubDistribution();
    out << "dist=";
    writeDiscreteDistribution(pd, out, globalAliases, writtenNames);
    comma = true;
  }
  else
  {
    try
    {
      auto& mix = dynamic_cast<const MixtureOfDiscreteDistributions&>(dist);
      size_t nd = mix.getNumberOfDistributions();
      for (size_t i = 0; i < nd; ++i)
      {
        if (comma)
          out << ",";
        out << "dist" + TextTools::toString(i + 1) + "=";
        writeDiscreteDistribution(mix.nDistribution(i), out, globalAliases, writtenNames);
        comma = true;
      }
      out << ",probas=(";
      for (size_t i = 0; i < nd; ++i)
      {
        out << mix.getNProbability(i);
        if (i != nd - 1)
          out << ",";
      }
      out << ")";
      for (size_t i = 1; i < nd; ++i)
      {
        writtenNames.push_back(mix.getNamespace() + "theta" + TextTools::toString(i));
      }
    }
    catch (bad_cast&)
    {}
  }
  if (dynamic_cast<const ExponentialDiscreteRateDistribution*>(&dist) ||
      dynamic_cast<const GammaDiscreteRateDistribution*>(&dist) ||
      dynamic_cast<const GaussianDiscreteRateDistribution*>(&dist))
  {
    if (comma)
      out << ",";
    out << "n="  << dist.getNumberOfCategories();
    comma = true;
  }
 
  try
  {
    auto& ps = dynamic_cast<const SimpleDiscreteDistribution&>(dist);
    size_t nd = ps.getNumberOfCategories();
    if (comma)
      out << ",";
    out << "values=(";
    for (size_t i = 0; i < nd; ++i)
    {
      out << ps.getCategory(i);
      if (i != nd - 1)
        out << ",";
    }
    out << "),probas=(";
    for (size_t i = 0; i < nd; ++i)
    {
      out << ps.getProbability(i);
      if (i != nd - 1)
        out << ",";
    }
    out << ")";
 
    auto range = ps.getRanges();
    if (range.size() != 0)
    {
      out << ", ranges=(";
      auto it = range.begin();
      while (it != range.end())
      {
        out << "V" << TextTools::toString(it->first);
        out << "[" << TextTools::toString(it->second[0]) << ";" << TextTools::toString(it->second[1]) << "]";
        it++;
        if (it != range.end())
          out << ",";
      }
    }
 
    out << ")";
 
    for (size_t i = 1; i < nd; ++i)
    {
      writtenNames.push_back(ps.getNamespace() + "theta" + TextTools::toString(i));
    }
    for (size_t i = 1; i < nd + 1; ++i)
    {
      writtenNames.push_back(ps.getNamespace() + "V" + TextTools::toString(i));
    }
 
    comma = true;
  }
  catch (bad_cast&)
  {}
 
 
  // Writing the parameters
  BppOParametrizableFormat bOP;
  bOP.write(dynamic_cast<const ParameterAliasable&>(dist), out, globalAliases, dist.getIndependentParameters().getParameterNames(), writtenNames, true, comma);
  out << ")";
}