gsl_rng.h

Go to the documentation of this file.

/*
 * $Id: gsl_rng.h 17393 2014-08-23 20:19:14Z predoehl $
 */

#ifndef GSL_RNG_H_INCLUDED_LIBKJB_UOFARIZONAVISION
#define GSL_RNG_H_INCLUDED_LIBKJB_UOFARIZONAVISION

#include <l_cpp/l_exception.h>
#include <gsl_cpp/gsl_util.h>
#include <vector>
#include <string>

#ifdef KJB_HAVE_GSL
#include "gsl/gsl_rng.h"    /* this is a GSL header, not our wrapper */
#else
#warning "Compiling GNU GSL wrapper without GNU GSL; it will not run properly"
typedef void gsl_rng;
typedef void gsl_rng_type;
gsl_rng_type    *gsl_rng_mt19937,
                *gsl_rng_ranlxs0,
                *gsl_rng_ranlxs1,
                *gsl_rng_ranlxs2,
                *gsl_rng_ranlxd1,
                *gsl_rng_ranlxd2,
                *gsl_rng_cmrg,
                *gsl_rng_mrg,
                *gsl_rng_taus2,
                *gsl_rng_gfsr4;
#endif


namespace kjb {


enum {
    GSL_RNG_MT19937,        
    GSL_RNG_RANLXS0,        
    GSL_RNG_RANLXS1,        
    GSL_RNG_RANLXS2,        
    GSL_RNG_RANLXD1,        
    GSL_RNG_RANLXD2,        
    GSL_RNG_CMRG,           
    GSL_RNG_MRG,            
    GSL_RNG_TAUS2,          
    GSL_RNG_GFSR4,          
};


#ifdef KJB_HAVE_GSL
std::string gsl_rng_serialize_implementation( const gsl_rng* );

void gsl_rng_deserialize_implementation( gsl_rng*, const std::string& );
#endif


template< unsigned KIND >
class Gsl_rng_basic {

    gsl_rng *m_rng;

    // I assume there will be fewer than, say, 100 kinds of RNG.
    static char kind()
    {
        return char( '0' + KIND );
    }

public:

    Gsl_rng_basic( const gsl_rng_type* type )
#ifdef KJB_HAVE_GSL
    :   m_rng( gsl_rng_alloc( type ) )
    {
        ETX_2( 00 == m_rng, "Gsl_rng_basic ctor: bad_alloc" );
    }
#else
    {
        KJB_THROW_2( Missing_dependency, "GNU GSL" );
    }
#endif

    Gsl_rng_basic< KIND >( const Gsl_rng_basic< KIND >& rng )
#ifdef KJB_HAVE_GSL
    :   m_rng( gsl_rng_clone( rng.m_rng ) )
    {
        ETX_2( 00 == m_rng, "Gsl_rng_basic copy ctor: bad_alloc" );
    }
#else
    {
        KJB_THROW_2( Missing_dependency, "GNU GSL" );
    }
#endif

    Gsl_rng_basic< KIND >& operator=( const Gsl_rng_basic< KIND >& that )
    {
        if ( this != &that )
        {
#ifdef KJB_HAVE_GSL
            gsl_rng_memcpy( m_rng, that.m_rng );
#endif
        }
        return *this;
    }

    void seed( unsigned long seed_val ) const
    {
#ifdef KJB_HAVE_GSL
        gsl_rng_set( m_rng, seed_val );
#endif
    }

    std::string serialize() const
    {
#ifdef KJB_HAVE_GSL
        return gsl_rng_serialize_implementation( m_rng ) + kind();
#endif
    }

    void deserialize( const std::string& state ) const
    {
#ifdef KJB_HAVE_GSL
        if ( state.size() < 2 || state[ state.size() - 1 ] != kind() )
        {
            GSL_ETX( GSL_EFAILED );
        }

        gsl_rng_deserialize_implementation( m_rng,
                            std::string( state.begin(), state.end() - 1 ) );
#endif
    }

    unsigned long get() const
    {
#ifdef KJB_HAVE_GSL
        return gsl_rng_get( m_rng );
#endif
    }

    const char* name() const
    {
#ifdef KJB_HAVE_GSL
        return gsl_rng_name( m_rng );
#endif
    }

    unsigned long max() const
    {
#ifdef KJB_HAVE_GSL
        return gsl_rng_max( m_rng );
#endif
    }

    unsigned long min() const
    {
#ifdef KJB_HAVE_GSL
        return gsl_rng_min( m_rng );
#endif
    }

    double uniform() const
    {
#ifdef KJB_HAVE_GSL
        return gsl_rng_uniform( m_rng );
#endif
    }

    double uniform_pos() const
    {
#ifdef KJB_HAVE_GSL
        return gsl_rng_uniform_pos( m_rng );
#endif
    }

    unsigned long uniform_int( unsigned long end_value ) const
    {
#ifdef KJB_HAVE_GSL
        unsigned long dev = gsl_rng_uniform_int( m_rng, end_value );

        if ( 0 == dev && end_value > max() )
        {
            GSL_ETX( GSL_EINVAL );
        }

        return dev;
#endif
    }

    operator const gsl_rng*() const
    {
        return m_rng;
    }

#ifdef KJB_HAVE_GSL
    ~Gsl_rng_basic()
    {
        gsl_rng_free( m_rng );
    }
#endif

};


#define Gsl_rng_template( Foo, foo, FOO )       \
    struct Foo : public Gsl_rng_basic< FOO >    \
    {                                           \
        Foo() : Gsl_rng_basic< FOO >( foo ) {}  \
    }

Gsl_rng_template(   Gsl_rng_mt19937,    gsl_rng_mt19937,    GSL_RNG_MT19937 );

Gsl_rng_template(   Gsl_rng_ranlxs0,    gsl_rng_ranlxs0,    GSL_RNG_RANLXS0 );

Gsl_rng_template(   Gsl_rng_ranlxs1,    gsl_rng_ranlxs1,    GSL_RNG_RANLXS1 );

Gsl_rng_template(   Gsl_rng_ranlxs2,    gsl_rng_ranlxs2,    GSL_RNG_RANLXS2 );

Gsl_rng_template(   Gsl_rng_ranlxd1,    gsl_rng_ranlxd1,    GSL_RNG_RANLXD1 );

Gsl_rng_template(   Gsl_rng_ranlxd2,    gsl_rng_ranlxd2,    GSL_RNG_RANLXD2 );

Gsl_rng_template(   Gsl_rng_cmrg,       gsl_rng_cmrg,       GSL_RNG_CMRG    );

Gsl_rng_template(   Gsl_rng_mrg,        gsl_rng_mrg,        GSL_RNG_MRG     );

Gsl_rng_template(   Gsl_rng_taus2,      gsl_rng_taus2,      GSL_RNG_TAUS2   );

Gsl_rng_template(   Gsl_rng_gfsr4,      gsl_rng_gfsr4,      GSL_RNG_GFSR4   );


typedef Gsl_rng_mt19937 Gsl_rng_default;

} // end ns kjb

#endif /* GSL_RNG_H_INCLUDED_LIBKJB_UOFARIZONAVISION */