gsl.hpp

Go to the documentation of this file.

//============================================================================
// vSMC/include/vsmc/rng/gsl.hpp
//----------------------------------------------------------------------------
//                         vSMC: Scalable Monte Carlo
//----------------------------------------------------------------------------
// Copyright (c) 2013,2014, Yan Zhou
// All rights reserved.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are met:
//
//   Redistributions of source code must retain the above copyright notice,
//   this list of conditions and the following disclaimer.
//
//   Redistributions in binary form must reproduce the above copyright notice,
//   this list of conditions and the following disclaimer in the documentation
//   and/or other materials provided with the distribution.
//
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS AS IS
// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
// POSSIBILITY OF SUCH DAMAGE.
//============================================================================

#ifndef VSMC_RNG_GSL_HPP
#define VSMC_RNG_GSL_HPP

#include <vsmc/rng/internal/common.hpp>
#include <vsmc/rng/generator_wrapper.hpp>
#include <gsl/gsl_rng.h>

#define VSMC_DEFINE_RNG_GSL_RNG_TYPE_POINTER(Generator, pointer) \
    template <> struct GSLRngTypePointer< GSL_RNG_TYPE_##Generator >         \
    {static const ::gsl_rng_type *get () {return ::gsl_rng_##pointer ;}};

#define VSMC_DEFINE_RNG_GSL_RNG_MIN_MAX(Generator, Min, Max) \
template <> struct GSLRngMinMax< GSL_RNG_TYPE_##Generator >                  \
{                                                                            \
    static VSMC_CONSTEXPR const uint32_t _Min =                              \
    static_cast<uint32_t>(Min##UL);                                          \
    static VSMC_CONSTEXPR const uint32_t _Max =                              \
    static_cast<uint32_t>(Max##UL);                                          \
    static VSMC_CONSTEXPR uint32_t min VSMC_MNE () {return _Min;}            \
    static VSMC_CONSTEXPR uint32_t max VSMC_MNE () {return _Max;}            \
}; // VSMC_DEFINE_RNG_GSL_RNG_MIN_MAX

namespace vsmc {

enum GSLRngType {
    GSL_RNG_TYPE_MT19937,   
    GSL_RNG_TYPE_RANLXS0,   
    GSL_RNG_TYPE_RANLXS1,   
    GSL_RNG_TYPE_RANLXS2,   
    GSL_RNG_TYPE_RANLXD1,   
    GSL_RNG_TYPE_RANLXD2,   
    GSL_RNG_TYPE_RANLUX,    
    GSL_RNG_TYPE_RANLUX389, 
    GSL_RNG_TYPE_CMRG,      
    GSL_RNG_TYPE_MRG,       
    GSL_RNG_TYPE_TAUS,      
    GSL_RNG_TYPE_TAUS2,     
    GSL_RNG_TYPE_GFSR4      
}; // enum GSLRngType

namespace internal {

template <GSLRngType> struct GSLRngTypePointer;

VSMC_DEFINE_RNG_GSL_RNG_TYPE_POINTER(MT19937,   mt19937)
VSMC_DEFINE_RNG_GSL_RNG_TYPE_POINTER(RANLXS0,   ranlxs0)
VSMC_DEFINE_RNG_GSL_RNG_TYPE_POINTER(RANLXS1,   ranlxs1)
VSMC_DEFINE_RNG_GSL_RNG_TYPE_POINTER(RANLXS2,   ranlxs2)
VSMC_DEFINE_RNG_GSL_RNG_TYPE_POINTER(RANLXD1,   ranlxd1)
VSMC_DEFINE_RNG_GSL_RNG_TYPE_POINTER(RANLXD2,   ranlxd2)
VSMC_DEFINE_RNG_GSL_RNG_TYPE_POINTER(RANLUX,    ranlux)
VSMC_DEFINE_RNG_GSL_RNG_TYPE_POINTER(RANLUX389, ranlux389)
VSMC_DEFINE_RNG_GSL_RNG_TYPE_POINTER(CMRG,      cmrg)
VSMC_DEFINE_RNG_GSL_RNG_TYPE_POINTER(MRG,       mrg)
VSMC_DEFINE_RNG_GSL_RNG_TYPE_POINTER(TAUS,      taus)
VSMC_DEFINE_RNG_GSL_RNG_TYPE_POINTER(TAUS2,     taus2)
VSMC_DEFINE_RNG_GSL_RNG_TYPE_POINTER(GFSR4,     gfsr4)

template <GSLRngType RngType> struct GSLRngMinMax;

VSMC_DEFINE_RNG_GSL_RNG_MIN_MAX(MT19937,   0, 0xFFFFFFFF)
VSMC_DEFINE_RNG_GSL_RNG_MIN_MAX(RANLXS0,   0, 0x00FFFFFF)
VSMC_DEFINE_RNG_GSL_RNG_MIN_MAX(RANLXS1,   0, 0x00FFFFFF)
VSMC_DEFINE_RNG_GSL_RNG_MIN_MAX(RANLXS2,   0, 0x00FFFFFF)
VSMC_DEFINE_RNG_GSL_RNG_MIN_MAX(RANLXD1,   0, 0xFFFFFFFF)
VSMC_DEFINE_RNG_GSL_RNG_MIN_MAX(RANLXD2,   0, 0xFFFFFFFF)
VSMC_DEFINE_RNG_GSL_RNG_MIN_MAX(RANLUX,    0, 0x00FFFFFF)
VSMC_DEFINE_RNG_GSL_RNG_MIN_MAX(RANLUX389, 0, 0x00FFFFFF)
VSMC_DEFINE_RNG_GSL_RNG_MIN_MAX(CMRG,      0, 0x7FFFFFFE)
VSMC_DEFINE_RNG_GSL_RNG_MIN_MAX(MRG,       0, 0x7FFFFFFE)
VSMC_DEFINE_RNG_GSL_RNG_MIN_MAX(TAUS,      0, 0xFFFFFFFF)
VSMC_DEFINE_RNG_GSL_RNG_MIN_MAX(TAUS2,     0, 0xFFFFFFFF)
VSMC_DEFINE_RNG_GSL_RNG_MIN_MAX(GFSR4,     0, 0xFFFFFFFF)
} // namespace vsmc::internal

template <GSLRngType RngType>
class GSLGenerator
{
    public :

    GSLGenerator () :
        rng_(::gsl_rng_alloc(internal::GSLRngTypePointer<RngType>::get())) {}

    GSLGenerator (const GSLGenerator<RngType> &other) :
        rng_(::gsl_rng_clone(other.rng_)) {}

    GSLGenerator<RngType> &operator= (const GSLGenerator<RngType> &other)
    {
        if (this != &other)
            ::gsl_rng_memcpy(rng_, other.rng_);

        return *this;
    }

#if VSMC_HAS_CXX11_RVALUE_REFERENCES
    GSLGenerator (GSLGenerator<RngType> &&other) :
        rng_(other.rng_) {other.rng_ = VSMC_NULLPTR;}

    GSLGenerator<RngType> &operator= (GSLGenerator<RngType> &&other)
    {
        using std::swap;

        if (this != &other)
            swap(rng_, other.rng_);

        return *this;
    }
#endif

    ~GSLGenerator ()
    {
        if (rng_ != VSMC_NULLPTR)
            ::gsl_rng_free(rng_);
    }

    void seed (unsigned long s)
    {
        if (rng_ != VSMC_NULLPTR)
            ::gsl_rng_set(rng_, s);
    }

    unsigned long generate ()
    {
        if (rng_ != VSMC_NULLPTR)
            return ::gsl_rng_get(rng_);

        return 0;
    }

    unsigned long min VSMC_MNE () const {return ::gsl_rng_min(rng_);}
    unsigned long max VSMC_MNE () const {return ::gsl_rng_max(rng_);}

    private :

    ::gsl_rng *rng_;
}; // class GSLGenerator

template <GSLRngType RngType>
class GSLEngine :
    public GeneratorWrapper<uint32_t, GSLGenerator<RngType>,
    internal::GSLRngMinMax<RngType> >
{
    typedef GeneratorWrapper<uint32_t, GSLGenerator<RngType>,
    internal::GSLRngMinMax<RngType> > base;

    public :

    typedef uint32_t result_type;

    GSLEngine (result_type s = 0) : base(s) {seed(s);}

    template <typename SeedSeq>
    explicit GSLEngine (SeedSeq &seq,
            typename cxx11::enable_if<internal::is_seed_seq<SeedSeq,
            result_type, GSLEngine<RngType>
            >::value>::type * = VSMC_NULLPTR) : base(seq) {seed(seq);}

    void seed (result_type s)
    {this->generator().seed(static_cast<unsigned long>(s));}

    template <typename SeedSeq>
    void seed (SeedSeq &seq,
            typename cxx11::enable_if<internal::is_seed_seq<SeedSeq,
            result_type, GSLEngine<RngType>
            >::value>::type * = VSMC_NULLPTR)
    {
        unsigned long s;
        seq.generate(&s, &s + 1);
        this->generator().seed(s);
    }
}; // class GSLEngine

typedef GSLEngine<GSL_RNG_TYPE_MT19937> GSL_MT19937;

typedef GSLEngine<GSL_RNG_TYPE_RANLXS0> GSL_RANLXS0;

typedef GSLEngine<GSL_RNG_TYPE_RANLXS1> GSL_RANLXS1;

typedef GSLEngine<GSL_RNG_TYPE_RANLXS2> GSL_RANLXS2;

typedef GSLEngine<GSL_RNG_TYPE_RANLXD1> GSL_RANLXD1;

typedef GSLEngine<GSL_RNG_TYPE_RANLXD2> GSL_RANLXD2;

typedef GSLEngine<GSL_RNG_TYPE_RANLUX> GSL_RANLUX;

typedef GSLEngine<GSL_RNG_TYPE_RANLUX389> GSL_RANLUX389;

typedef GSLEngine<GSL_RNG_TYPE_CMRG> GSL_CMRG;

typedef GSLEngine<GSL_RNG_TYPE_MRG> GSL_MRG;

typedef GSLEngine<GSL_RNG_TYPE_TAUS> GSL_TAUS;

typedef GSLEngine<GSL_RNG_TYPE_TAUS2> GSL_TAUS2;

typedef GSLEngine<GSL_RNG_TYPE_GFSR4> GSL_GFSR4;

} // namespace vsmc

#endif // VSMC_RNG_GSL_HPP