RandomTools.h

Go to the documentation of this file.

// SPDX-FileCopyrightText: The Bio++ Development Group
//
// SPDX-License-Identifier: CECILL-2.1

#ifndef BPP_NUMERIC_RANDOM_RANDOMTOOLS_H
#define BPP_NUMERIC_RANDOM_RANDOMTOOLS_H


#include "../../Exceptions.h"
#include "../VectorExceptions.h"
#include "../VectorTools.h"

// From the STL:
#include <cmath>
#include <cassert>
#include <ctime>
#include <vector>
#include <algorithm>
#include <random>

namespace bpp
{
class RandomTools
{
public:
  RandomTools() {}
  virtual ~RandomTools() {}

private:
  static double incompletebetafe(double a,
      double b,
      double x,
      double big,
      double biginv);
  static double incompletebetafe2(double a,
      double b,
      double x,
      double big,
      double biginv);
  static double incompletebetaps(double a,
      double b,
      double x,
      double maxgam);

public:
  static std::random_device RANDOM_DEVICE;
  static std::mt19937 DEFAULT_GENERATOR;

  static void setSeed(std::mt19937::result_type seed)
  {
    DEFAULT_GENERATOR.seed(seed);
  }

  static double giveRandomNumberBetweenZeroAndEntry(double entry)
  {
    std::uniform_real_distribution<double> dis(0, entry);
    return dis(DEFAULT_GENERATOR);
  }

  static bool flipCoin(double prob = 0.5)
  {
    std::bernoulli_distribution d(prob);
    return d(DEFAULT_GENERATOR);
  }

  template<class intType>
  static intType giveIntRandomNumberBetweenZeroAndEntry(intType entry)
  {
    if (entry == 0)
      throw Exception("RandomTools::giveIntRandomNumberBetweenZeroAndEntry. Entry must be at least 1.");
    std::uniform_int_distribution<intType> dis(0, entry - 1);
    return dis(DEFAULT_GENERATOR);
  }

  static double randGaussian(double mean, double variance)
  {
    std::normal_distribution<double> dis(mean, sqrt(variance));
    return dis(DEFAULT_GENERATOR);
  }

  static double randGamma(double alpha)
  {
    std::gamma_distribution<double> dis(alpha, 1.);
    return dis(DEFAULT_GENERATOR);
  }

  static double randGamma(double alpha, double beta)
  {
    std::gamma_distribution<double> dis(alpha, beta);
    return dis(DEFAULT_GENERATOR);
  }

  static double randBeta(double alpha, double beta);

  static double randExponential(double mean)
  {
    std::exponential_distribution<double> dis(mean);
    return dis(DEFAULT_GENERATOR);
  }

  template<class T>
  static T pickOne(std::vector<T>& v, bool replace = false)
  {
    if (v.empty())
      throw EmptyVectorException<T>("RandomTools::pickOne: input vector is empty", &v);
    size_t pos = RandomTools::giveIntRandomNumberBetweenZeroAndEntry<size_t>(v.size());
    if (replace)
      return v[pos];
    else
    {
      T e = v[pos];
      v[pos] = v.back();
      v.pop_back();
      return e;
    }
  }

  template<class T>
  static T pickOne(const std::vector<T>& v)
  {
    if (v.empty())
      throw EmptyVectorException<T>("RandomTools::pickOne: input vector is empty", &v);
    size_t pos = RandomTools::giveIntRandomNumberBetweenZeroAndEntry<size_t>(v.size());
    return v[pos];
  }

  template<class T>
  static void getSample(const std::vector<T>& vin, std::vector<T>& vout, bool replace = false)
  {
    if (vout.size() > vin.size() && !replace)
      throw IndexOutOfBoundsException("RandomTools::getSample: size exceeded v.size.", vout.size(), 0, vin.size());
    if (replace)
    {
      for (size_t i = 0; i < vout.size(); ++i)
      {
        vout[i] = pickOne(vin);
      }
    }
    else
    {
      std::vector<size_t> hat(vin.size());
      std::iota(hat.begin(), hat.end(), 0);
      std::shuffle(hat.begin(), hat.end(), DEFAULT_GENERATOR);
      for (size_t i = 0; i < vout.size(); i++)
      {
        vout[i] = vin[hat[i]];
      }
    }
  }

  template<class T>
  static T pickOne(std::vector<T>& v, std::vector<double>& w, bool replace = false)
  {
    if (v.empty())
      throw EmptyVectorException<T>("RandomTools::pickOne (with weight): input vector is empty", &v);
    // Compute cumulative sum of weights:
    std::vector<double> sumw = VectorTools::cumSum(w);
    // Convert to cumulative distribution:
    sumw /= sumw.back();
    // Get random positions:
    double prob = RandomTools::giveRandomNumberBetweenZeroAndEntry(1.0);
    size_t pos = v.size() - 1;
    for (size_t i = 0; i < v.size(); ++i)
    {
      if (prob < sumw[i])
      {
        pos = i;
        break;
      }
    }
    if (replace)
      return v[pos];
    else
    {
      T e = v[pos];
      v[pos] = v.back();
      v.pop_back();
      w[pos] = w.back();
      w.pop_back();
      return e;
    }
  }

  template<class T>
  static T pickOne(const std::vector<T>& v, const std::vector<double>& w)
  {
    if (v.empty())
      throw EmptyVectorException<T>("RandomTools::pickOne (with weight): input vector is empty", &v);
    // Compute cumulative sum of weights:
    std::vector<double> sumw = VectorTools::cumSum(w);
    // Convert to cumulative distribution:
    sumw /= sumw.back();
    // Get random positions:
    double prob = RandomTools::giveRandomNumberBetweenZeroAndEntry(1.0);
    size_t pos = v.size() - 1;
    for (size_t i = 0; i < v.size(); ++i)
    {
      if (prob < sumw[i])
      {
        pos = i;
        break;
      }
    }
    return v[pos];
  }

  static size_t pickFromCumSum(const std::vector<double>& w)
  {
    double prob = RandomTools::giveRandomNumberBetweenZeroAndEntry(1.0);
    size_t pos = 0;
    while (pos < w.size() - 1)
    {
      if (prob <= w[pos])
        return pos;
      pos += 1;
    }
    return w.size() - 1;
  }

  template<class T>
  static void getSample(const std::vector<T>& vin, const std::vector<double>& w, std::vector<T>& vout, bool replace = false)
  {
    if (vout.size() > vin.size() && !replace)
      throw IndexOutOfBoundsException("RandomTools::getSample (with weights): size exceeded v.size.", vout.size(), 0, vin.size());
    std::vector<size_t> hat(vin.size());
    std::iota(hat.begin(), hat.end(), 0);
    if (replace)
    {
      for (size_t i = 0; i < vout.size(); i++)
      {
        vout[i] = vin[pickOne(hat, w)];
      }
    }
    else
    {
      std::vector<double> w2(w); // non const copy
      for (size_t i = 0; i < vout.size(); i++)
      {
        vout[i] = vin[pickOne(hat, w2, false)];
      }
    }
  }

  static std::vector<size_t> randMultinomial(size_t n, const std::vector<double>& probs);

  static double qNorm(double prob);

  static double qNorm(double prob, double mu, double sigma);

  static double lnGamma(double alpha) { return std::lgamma(alpha); }

  static double incompleteGamma(double x, double alpha, double ln_gamma_alpha);


  static double qChisq(double prob, double v);

  static double pChisq(double x, double v)
  {
    if (x < 0) return 0;
    return pGamma(x, v / 2, 0.5);
  }

  static double qGamma(double prob, double alpha, double beta)
  {
    return qChisq(prob, 2.0 * (alpha)) / (2.0 * (beta));
  }

  static double pGamma(double x, double alpha, double beta)
  {
    if (alpha < 0) throw Exception("RandomTools::pGamma. Negative alpha is not allowed.");
    if (beta < 0) throw Exception("RandomTools::pGamma. Negative beta is not allowed.");
    if (alpha == 0.) return 1.;
    return incompleteGamma(beta * x, alpha, lnGamma(alpha));
  }

  static double pNorm(double z);

  /* @brief Normal cumulative function.
   *
   * Returns Prob{x<=z} where x ~ N(mu,sigma^2)

   * @param z the value.
   * @param mu The mean of the distribution
   * @param sigma The standard deviation of the distribution
   * @return The corresponding probability.
   */

  static double pNorm(double z, double mu, double sigma);

  static double lnBeta(double alpha, double beta);

  static double incompleteBeta(double x, double alpha, double beta);
  static double pBeta(double x, double alpha, double beta)
  {
    return incompleteBeta(x, alpha, beta);
  }

  static double qBeta(double prob, double alpha, double beta);

private:
  static double DblGammaGreaterThanOne(double dblAlpha);
  static double DblGammaLessThanOne(double dblAlpha);
};
} // end of namespace bpp.
#endif // BPP_NUMERIC_RANDOM_RANDOMTOOLS_H