randist.hpp

Go to the documentation of this file.

/*
 * $Id: randist.hpp 293 2012-12-17 20:27:36Z jdl3 $
 * Copyright (C) 2010, 2020 John D Lamb
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 * 
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */
 
#ifndef CCGSL_RANDIST_HPP
#define CCGSL_RANDIST_HPP
 
#include<new>
#ifdef __GXX_EXPERIMENTAL_CXX0X__
#include<stdexcept>
#include<type_traits>
#endif
#include<gsl/gsl_randist.h>
#include"rng.hpp"
 
namespace gsl {
  namespace ran {
    class discrete_t {
    public:
      discrete_t(){
    ccgsl_pointer = 0;
    // just plausibly we could fail to allocate count
    try { count = new size_t; } catch( std::bad_alloc& e ){
      throw e;
    }
    *count = 1; // initially there is just one reference to ccgsl_pointer
      }
      // Refines random access container
      // Refines assignable
#ifndef DOXYGEN_SKIP
      explicit discrete_t( size_t K, double const* P  ){
    ccgsl_pointer = gsl_ran_discrete_preproc( K, P );
    // just plausibly we could allocate gsl_ran_discrete_t but not count
    try { count = new size_t; } catch( std::bad_alloc& e ){
      // try to tidy up before rethrowing
      gsl_ran_discrete_free( ccgsl_pointer );
      throw e;
    }
    *count = 1; // initially there is just one reference to ccgsl_pointer
      }
#endif // DOXYGEN_SKIP
      template<typename ARRAY>
      explicit discrete_t( ARRAY const& P ){
    ccgsl_pointer = gsl_ran_discrete_preproc( P.size(), P.data() );
    // just plausibly we could allocate gsl_ran_discrete_t but not count
    try { count = new size_t; } catch( std::bad_alloc& e ){
      // try to tidy up before rethrowing
      gsl_ran_discrete_free( ccgsl_pointer );
      throw e;
    }
    *count = 1; // initially there is just one reference to ccgsl_pointer
      }
      explicit discrete_t( gsl_ran_discrete_t* v ){
    ccgsl_pointer = v;
    // just plausibly we could fail to allocate count: no further action needed.
    count = new size_t;
    *count = 1; // initially there is just one reference to ccgsl_pointer
      }
      // copy constructor
      discrete_t( discrete_t const& v ){ ccgsl_pointer = v.ccgsl_pointer;
    count = v.count; if( count != 0 ) ++*count; }
      // assignment operator
      discrete_t& operator=( discrete_t const& v ){
    // first, possibly delete anything pointed to by this
    if( --*count == 0 ){
      if( ccgsl_pointer != 0 ) gsl_ran_discrete_free( ccgsl_pointer );
      delete count;
    } // Then copy
    ccgsl_pointer = v.ccgsl_pointer; count = v.count; if( count != 0 ) ++*count; return *this;
      }
      // destructor
      ~discrete_t(){
    if( --*count == 0 ){
      // could have allocated null pointer
      if( ccgsl_pointer != 0 ) gsl_ran_discrete_free( ccgsl_pointer );
      delete count;
    }
      }
      // Refines equality comparable
      // == operator
      bool operator==( discrete_t const& v ) const { return ccgsl_pointer == v.ccgsl_pointer; }
      // != operator
      bool operator!=( discrete_t const& v ) const { return not operator==( v ); }
      // Refines forward container
      // Refines less than comparable
      // operator<
      bool operator<( discrete_t const& v ) const { return ccgsl_pointer < v.ccgsl_pointer; }
      // operator>
      bool operator>( discrete_t const& v ) const { return ccgsl_pointer > v.ccgsl_pointer; }
      // operator<=
      bool operator<=( discrete_t const& v ) const { return ccgsl_pointer <= v.ccgsl_pointer; }
      // operator>=
      bool operator>=( discrete_t const& v ) const { return ccgsl_pointer >= v.ccgsl_pointer; }
      bool empty() const { return ccgsl_pointer == 0; }
      // swap() --- should work even if sizes don't match
      void swap( discrete_t& v ){
    gsl_ran_discrete_t* tmp = ccgsl_pointer; ccgsl_pointer = v.ccgsl_pointer;
    v.ccgsl_pointer = tmp; size_t* tmp2 = count; count = v.count; v.count = tmp2; }
    private:
      gsl_ran_discrete_t* ccgsl_pointer;
      size_t* count;
    public:
      // shared reference functions
      gsl_ran_discrete_t* get() const { return ccgsl_pointer; }
      bool unique() const { return count != 0 and *count == 1; }
      size_t use_count() const { return count == 0 ? 0 : *count; }
#ifdef __GXX_EXPERIMENTAL_CXX0X__
      explicit
#endif
      operator bool() const { return ccgsl_pointer != 0; }
    };
    // Functions
    inline unsigned int bernoulli( rng const& r, double p ){
      return gsl_ran_bernoulli( r.get(), p ); } 
    inline double bernoulli_pdf( unsigned int const k, double p ){
      return gsl_ran_bernoulli_pdf( k, p ); } 
    inline double beta( rng const& r, double const a, double const b ){
      return gsl_ran_beta( r.get(), a, b ); } 
    inline double beta_pdf( double const x, double const a, double const b ){
      return gsl_ran_beta_pdf( x, a, b ); } 
    inline unsigned int binomial( rng const& r, double p, unsigned int n ){
      return gsl_ran_binomial( r.get(), p, n ); } 
    inline unsigned int binomial_knuth( rng const& r, double p, unsigned int n ){
      return gsl_ran_binomial_knuth( r.get(), p, n ); } 
    inline unsigned int binomial_tpe( rng const& r, double p, unsigned int n ){
      return gsl_ran_binomial_tpe( r.get(), p, n ); } 
    inline double binomial_pdf( unsigned int const k, double const p, unsigned int const n ){
      return gsl_ran_binomial_pdf( k, p, n ); } 
    inline double exponential( rng const& r, double const mu ){
      return gsl_ran_exponential( r.get(), mu ); } 
    inline double exponential_pdf( double const x, double const mu ){
      return gsl_ran_exponential_pdf( x, mu ); } 
    inline double exppow( rng const& r, double const a, double const b ){
      return gsl_ran_exppow( r.get(), a, b ); } 
    inline double exppow_pdf( double const x, double const a, double const b ){
      return gsl_ran_exppow_pdf( x, a, b ); } 
    inline double cauchy( rng const& r, double const a ){
      return gsl_ran_cauchy( r.get(), a ); } 
    inline double cauchy_pdf( double const x, double const a ){
      return gsl_ran_cauchy_pdf( x, a ); } 
    inline double chisq( rng const& r, double const nu ){
      return gsl_ran_chisq( r.get(), nu ); } 
    inline double chisq_pdf( double const x, double const nu ){
      return gsl_ran_chisq_pdf( x, nu ); } 
    inline void dirichlet( rng const& r, size_t const K, double const alpha[], double theta[] ){
      gsl_ran_dirichlet( r.get(), K, alpha, theta ); } 
    inline double dirichlet_pdf( size_t const K, double const alpha[], double const theta[] ){
      return gsl_ran_dirichlet_pdf( K, alpha, theta ); } 
    inline double dirichlet_lnpdf( size_t const K, double const alpha[], double const theta[] ){
      return gsl_ran_dirichlet_lnpdf( K, alpha, theta ); } 
    inline double erlang( rng const& r, double const a, double const n ){
      return gsl_ran_erlang( r.get(), a, n ); } 
    inline double erlang_pdf( double const x, double const a, double const n ){
      return gsl_ran_erlang_pdf( x, a, n ); } 
    inline double fdist( rng const& r, double const nu1, double const nu2 ){
      return gsl_ran_fdist( r.get(), nu1, nu2 ); } 
    inline double fdist_pdf( double const x, double const nu1, double const nu2 ){
      return gsl_ran_fdist_pdf( x, nu1, nu2 ); } 
    inline double flat( rng const& r, double const a, double const b ){
      return gsl_ran_flat( r.get(), a, b ); } 
    inline double flat_pdf( double x, double const a, double const b ){
      return gsl_ran_flat_pdf( x, a, b ); } 
    inline double gamma( rng const& r, double const a, double const b ){
      return gsl_ran_gamma( r.get(), a, b ); } 
    inline double gamma_int( rng const& r, unsigned int const a ){
      return gsl_ran_gamma_int( r.get(), a ); } 
    inline double gamma_pdf( double const x, double const a, double const b ){
      return gsl_ran_gamma_pdf( x, a, b ); } 
    inline double gamma_mt( rng const& r, double const a, double const b ){
      return gsl_ran_gamma_mt( r.get(), a, b ); } 
    inline double gamma_knuth( rng const& r, double const a, double const b ){
      return gsl_ran_gamma_knuth( r.get(), a, b ); } 
    inline double gaussian( rng const& r, double const sigma ){
      return gsl_ran_gaussian( r.get(), sigma ); } 
    inline double gaussian_ratio_method( rng const& r, double const sigma ){
      return gsl_ran_gaussian_ratio_method( r.get(), sigma ); } 
    inline double gaussian_ziggurat( rng const& r, double const sigma ){
      return gsl_ran_gaussian_ziggurat( r.get(), sigma ); } 
    inline double gaussian_pdf( double const x, double const sigma ){
      return gsl_ran_gaussian_pdf( x, sigma ); } 
    inline double ugaussian( rng const& r ){
      return gsl_ran_ugaussian( r.get() ); } 
    inline double ugaussian_ratio_method( rng const& r ){
      return gsl_ran_ugaussian_ratio_method( r.get() ); } 
    inline double ugaussian_pdf( double const x ){
      return gsl_ran_ugaussian_pdf( x ); } 
    inline double gaussian_tail( rng const& r, double const a, double const sigma ){
      return gsl_ran_gaussian_tail( r.get(), a, sigma ); } 
    inline double gaussian_tail_pdf( double const x, double const a, double const sigma ){
      return gsl_ran_gaussian_tail_pdf( x, a, sigma ); } 
    inline double ugaussian_tail( rng const& r, double const a ){
      return gsl_ran_ugaussian_tail( r.get(), a ); } 
    inline double ugaussian_tail_pdf( double const x, double const a ){
      return gsl_ran_ugaussian_tail_pdf( x, a ); } 
#ifndef DOXYGEN_SKIP
    inline void bivariate_gaussian( rng const& r, double sigma_x, double sigma_y,
                    double rho, double* x, double* y ){
      gsl_ran_bivariate_gaussian( r.get(), sigma_x, sigma_y, rho, x, y ); } 
#endif // DOXYGEN_SKIP
    inline void bivariate_gaussian( rng const& r, double sigma_x, double sigma_y,
                    double rho, double& x, double& y ){
      gsl_ran_bivariate_gaussian( r.get(), sigma_x, sigma_y, rho, &x, &y ); } 
    inline double bivariate_gaussian_pdf( double const x, double const y, double const sigma_x,
                      double const sigma_y, double const rho ){
      return gsl_ran_bivariate_gaussian_pdf( x, y, sigma_x, sigma_y, rho ); } 
    inline double landau( rng const& r ){
      return gsl_ran_landau( r.get() ); } 
    inline double landau_pdf( double const x ){ return gsl_ran_landau_pdf( x ); } 
    inline unsigned int geometric( rng const& r, double const p ){
      return gsl_ran_geometric( r.get(), p ); } 
    inline double geometric_pdf( unsigned int const k, double const p ){
      return gsl_ran_geometric_pdf( k, p ); } 
    inline unsigned int hypergeometric( rng const& r, unsigned int n1, unsigned int n2, unsigned int t ){
      return gsl_ran_hypergeometric( r.get(), n1, n2, t ); } 
    inline double hypergeometric_pdf( unsigned int const k, unsigned int const n1,
                      unsigned int const n2, unsigned int t ){
      return gsl_ran_hypergeometric_pdf( k, n1, n2, t ); } 
    inline double gumbel1( rng const& r, double const a, double const b ){
      return gsl_ran_gumbel1( r.get(), a, b ); } 
    inline double gumbel1_pdf( double const x, double const a, double const b ){
      return gsl_ran_gumbel1_pdf( x, a, b ); } 
    inline double gumbel2( rng const& r, double const a, double const b ){
      return gsl_ran_gumbel2( r.get(), a, b ); } 
    inline double gumbel2_pdf( double const x, double const a, double const b ){
      return gsl_ran_gumbel2_pdf( x, a, b ); } 
    inline double logistic( rng const& r, double const a ){
      return gsl_ran_logistic( r.get(), a ); } 
    inline double logistic_pdf( double const x, double const a ){
      return gsl_ran_logistic_pdf( x, a ); } 
    inline double lognormal( rng const& r, double const zeta, double const sigma ){
      return gsl_ran_lognormal( r.get(), zeta, sigma ); } 
    inline double lognormal_pdf( double const x, double const zeta, double const sigma ){
      return gsl_ran_lognormal_pdf( x, zeta, sigma ); } 
    inline unsigned int logarithmic( rng const& r, double const p ){
      return gsl_ran_logarithmic( r.get(), p ); } 
    inline double logarithmic_pdf( unsigned int const k, double const p ){
      return gsl_ran_logarithmic_pdf( k, p ); } 
    inline void multinomial( rng const& r, size_t const K, unsigned int const N,
                 double const p[], unsigned int n[] ){
      gsl_ran_multinomial( r.get(), K, N, p, n ); } 
    inline double multinomial_pdf( size_t const K, double const p[], unsigned int const n[] ){
      return gsl_ran_multinomial_pdf( K, p, n ); } 
    inline double multinomial_lnpdf( size_t const K, double const p[], unsigned int const n[] ){
      return gsl_ran_multinomial_lnpdf( K, p, n ); } 
    inline unsigned int negative_binomial( rng const& r, double p, double n ){
      return gsl_ran_negative_binomial( r.get(), p, n ); } 
    inline double negative_binomial_pdf( unsigned int const k, double const p, double n ){
      return gsl_ran_negative_binomial_pdf( k, p, n ); } 
    inline unsigned int pascal( rng const& r, double p, unsigned int n ){
      return gsl_ran_pascal( r.get(), p, n ); } 
    inline double pascal_pdf( unsigned int const k, double const p, unsigned int n ){
      return gsl_ran_pascal_pdf( k, p, n ); } 
    inline double pareto( rng const& r, double a, double const b ){
      return gsl_ran_pareto( r.get(), a, b ); } 
    inline double pareto_pdf( double const x, double const a, double const b ){
      return gsl_ran_pareto_pdf( x, a, b ); } 
    inline unsigned int poisson( rng const& r, double mu ){
      return gsl_ran_poisson( r.get(), mu ); } 
    inline void poisson_array( rng const& r, size_t n, unsigned int array[], double mu ){
      gsl_ran_poisson_array( r.get(), n, array, mu ); } 
    inline double poisson_pdf( unsigned int const k, double const mu ){
      return gsl_ran_poisson_pdf( k, mu ); } 
    inline double rayleigh( rng const& r, double const sigma ){
      return gsl_ran_rayleigh( r.get(), sigma ); } 
    inline double rayleigh_pdf( double const x, double const sigma ){
      return gsl_ran_rayleigh_pdf( x, sigma ); } 
    inline double rayleigh_tail( rng const& r, double const a, double const sigma ){
      return gsl_ran_rayleigh_tail( r.get(), a, sigma ); } 
    inline double rayleigh_tail_pdf( double const x, double const a, double const sigma ){
      return gsl_ran_rayleigh_tail_pdf( x, a, sigma ); } 
    inline double tdist( rng const& r, double const nu ){
      return gsl_ran_tdist( r.get(), nu ); } 
    inline double tdist_pdf( double const x, double const nu ){
      return gsl_ran_tdist_pdf( x, nu ); } 
    inline double laplace( rng const& r, double const a ){
      return gsl_ran_laplace( r.get(), a ); } 
    inline double laplace_pdf( double const x, double const a ){
      return gsl_ran_laplace_pdf( x, a ); } 
    inline double levy( rng const& r, double const c, double const alpha ){
      return gsl_ran_levy( r.get(), c, alpha ); } 
    inline double levy_skew( rng const& r, double const c, double const alpha, double const beta ){
      return gsl_ran_levy_skew( r.get(), c, alpha, beta ); } 
    inline double weibull( rng const& r, double const a, double const b ){
      return gsl_ran_weibull( r.get(), a, b ); } 
    inline double weibull_pdf( double const x, double const a, double const b ){
      return gsl_ran_weibull_pdf( x, a, b ); } 
#ifndef DOXYGEN_SKIP
    inline void dir_2d( rng const& r, double* x, double* y ){
      gsl_ran_dir_2d( r.get(), x, y ); } 
#endif // DOXYGEN_SKIP
    inline void dir_2d( rng const& r, double& x, double& y ){
      gsl_ran_dir_2d( r.get(), &x, &y ); } 
#ifndef DOXYGEN_SKIP
    inline void dir_2d_trig_method( rng const& r, double* x, double* y ){
      gsl_ran_dir_2d_trig_method( r.get(), x, y ); } 
#endif // DOXYGEN_SKIP
    inline void dir_2d_trig_method( rng const& r, double& x, double& y ){
      gsl_ran_dir_2d_trig_method( r.get(), &x, &y ); } 
#ifndef DOXYGEN_SKIP
    inline void dir_3d( rng const& r, double* x, double* y, double* z ){
      gsl_ran_dir_3d( r.get(), x, y, z ); } 
#endif // DOXYGEN_SKIP
    inline void dir_3d( rng const& r, double& x, double& y, double& z ){
      gsl_ran_dir_3d( r.get(), &x, &y, &z ); } 
#ifndef DOXYGEN_SKIP
    inline void dir_nd( rng const& r, size_t n, double* x ){
      gsl_ran_dir_nd( r.get(), n, x ); } 
#endif // DOXYGEN_SKIP
    template<typename DATA>
    inline void dir_nd( rng const& r, DATA& x ){
      gsl_ran_dir_nd( r.get(), x.size(), x.data() ); } 
#ifndef DOXYGEN_SKIP
    inline void shuffle( rng const& r, void* base, size_t nmembm, size_t size ){
      gsl_ran_shuffle( r.get(), base, nmembm, size ); } 
#endif // DOXYGEN_SKIP
    template<typename ARRAY>
    inline void shuffle( rng const& r, ARRAY& base ){
      gsl_ran_shuffle( r.get(), base.data(), base.size(),
               sizeof( typename ARRAY::value_type ) ); } 
#ifndef DOXYGEN_SKIP
    inline int choose( rng const& r, void* dest, size_t k, void* src, size_t n, size_t size ){
      return gsl_ran_choose( r.get(), dest, k, src, n, size ); } 
#endif // DOXYGEN_SKIP
    template<typename ARRAY1, typename ARRAY2>
    inline int choose( rng const& r, ARRAY1& dest, ARRAY2& src ){
#ifdef __GXX_EXPERIMENTAL_CXX0X__
      if( not std::is_same<typename ARRAY1::value_type,typename ARRAY2::value_type>::value )
    return GSL_EFAILED;
#endif
      return gsl_ran_choose( r.get(), dest.data(), dest.size(), src.data(), src.size(),
                 sizeof( typename ARRAY1::value_type ) ); } 
#ifndef DOXYGEN_SKIP
    inline void sample( rng const& r, void* dest, size_t k, void* src, size_t n, size_t size ){
      gsl_ran_sample( r.get(), dest, k, src, n, size ); } 
#endif // DOXYGEN_SKIP
    template<typename ARRAY1, typename ARRAY2>
    inline void sample( rng const& r, ARRAY1& dest, ARRAY2& src ){
#ifdef __GXX_EXPERIMENTAL_CXX0X__
      if( not std::is_same<typename ARRAY1::value_type,typename ARRAY2::value_type>::value )
    throw std::invalid_argument( "ran::sample: arrays hold different types of data." );
#endif
      gsl_ran_sample( r.get(), dest.data(), dest.size(), src.data(), src.size(),
              sizeof( typename ARRAY1::value_type ) ); } 
    inline size_t discrete( rng const& r, discrete_t const& g ){
      return gsl_ran_discrete( r.get(), g.get() ); } 
    inline double discrete_pdf( size_t k, discrete_t const& g ){
      return gsl_ran_discrete_pdf( k, g.get() ); } 
  }
}
 
#endif