RandomTools.h

Go to the documentation of this file.

//
// File: RandomTools.h
// Authors:
//   Julien Dutheil
//   Sylvain Gaillard
// Last modified: 2008-11-06 00:00:00
//
 
/*
  Copyright or © or Copr. Bio++ Development Team, (November 17, 2004)
  
  This software is a computer program whose purpose is to provide classes
  for numerical calculus.
  
  This software is governed by the CeCILL license under French law and
  abiding by the rules of distribution of free software. You can use,
  modify and/ or redistribute the software under the terms of the CeCILL
  license as circulated by CEA, CNRS and INRIA at the following URL
  "http://www.cecill.info".
  
  As a counterpart to the access to the source code and rights to copy,
  modify and redistribute granted by the license, users are provided only
  with a limited warranty and the software's author, the holder of the
  economic rights, and the successive licensors have only limited
  liability.
  
  In this respect, the user's attention is drawn to the risks associated
  with loading, using, modifying and/or developing or reproducing the
  software by the user in light of its specific status of free software,
  that may mean that it is complicated to manipulate, and that also
  therefore means that it is reserved for developers and experienced
  professionals having in-depth computer knowledge. Users are therefore
  encouraged to load and test the software's suitability as regards their
  requirements in conditions enabling the security of their systems and/or
  data to be ensured and, more generally, to use and operate it in the
  same conditions as regards security.
  
  The fact that you are presently reading this means that you have had
  knowledge of the CeCILL license and that you accept its terms.
*/
 
#ifndef BPP_NUMERIC_RANDOM_RANDOMTOOLS_H
#define BPP_NUMERIC_RANDOM_RANDOMTOOLS_H
 
 
#include "../../Exceptions.h"
#include "../VectorExceptions.h"
#include "../VectorTools.h"
#include "RandomFactory.h"
 
// From the STL:
#include <cmath>
#include <cassert>
#include <ctime>
#include <vector>
#include <algorithm>
 
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 RandomFactory* DEFAULT_GENERATOR;
 
  static double giveRandomNumberBetweenZeroAndEntry(double entry, const RandomFactory& generator = * DEFAULT_GENERATOR);
 
  static bool flipCoin(const RandomFactory& generator = * DEFAULT_GENERATOR);
 
  template<class intType>
  static intType giveIntRandomNumberBetweenZeroAndEntry(intType entry, const RandomFactory& generator = * DEFAULT_GENERATOR)
  {
    return static_cast<intType>(giveRandomNumberBetweenZeroAndEntry(static_cast<double>(entry), generator));
  }
 
  static void setSeed(long seed);
 
  static double randGaussian(double mean, double variance, const RandomFactory& generator = * DEFAULT_GENERATOR);
 
  static double randGamma(double dblAlpha, const RandomFactory& generator = * DEFAULT_GENERATOR);
 
  static double randGamma(double alpha, double beta, const RandomFactory& generator = * DEFAULT_GENERATOR);
 
  static double randBeta(double alpha, double beta, const RandomFactory& generator = * DEFAULT_GENERATOR);
 
  static double randExponential(double mean, const RandomFactory& generator = * 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::random_shuffle(hat.begin(), hat.end());
      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, const RandomFactory& generator);
  static double DblGammaLessThanOne(double dblAlpha, const RandomFactory& generator);
};
} // end of namespace bpp.
#endif // BPP_NUMERIC_RANDOM_RANDOMTOOLS_H