random.hpp

Go to the documentation of this file.

#ifndef VMCMC_RANDOM_H_
#define VMCMC_RANDOM_H_

#include <vmcmc/logger.hpp>

#include <random>
#include <type_traits>
#include <cmath>
#include <atomic>
#include <thread>

#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/triangular.hpp>

namespace vmcmc
{

namespace ublas = boost::numeric::ublas;

template <typename EngineT>
class RandomPrototype
{
public:
    typedef EngineT engine_type;
    typedef typename engine_type::result_type result_type;

    static void Seed(result_type seed);

    static RandomPrototype& Instance();

protected:
    static std::atomic<result_type> sSeed;

    RandomPrototype();
    virtual ~RandomPrototype();

    RandomPrototype(const RandomPrototype& other) = delete;
    RandomPrototype(RandomPrototype&& other) = delete;
    void operator=(const RandomPrototype& other) = delete;
    void operator=(RandomPrototype&& other) = delete;

public:
    engine_type& GetEngine() { return fEngine; }

    template <typename FloatT>
    typename std::enable_if<std::is_floating_point<FloatT>::value, FloatT>::type
    Uniform(FloatT min, FloatT max, bool minIncluded, bool maxIncluded);

    template <typename FloatT = double>
    typename std::enable_if<std::is_floating_point<FloatT>::value, FloatT>::type
    Uniform(FloatT min = 0.0, FloatT max = 1.0);

    template <typename IntegerT>
    typename std::enable_if<std::is_integral<IntegerT>::value, IntegerT>::type
    Uniform(IntegerT inclMin, IntegerT inclMax);

    template <typename FloatT = double>
    bool Bool(FloatT probability = 0.5);

    template <typename FloatT = double>
    FloatT Normal(FloatT mean = 0.0, FloatT sigma = 1.0);

    template <typename FloatT = double>
    FloatT StudentT(FloatT n = 1.0, FloatT mean = 0.0, FloatT sigma = 1.0);

    template <typename FloatT>
    inline FloatT Exponential(FloatT tau);

    template <typename IntegerT = uint32_t>
    typename std::enable_if<std::is_integral<IntegerT>::value, IntegerT>::type
    Poisson(double mean);

    template <typename FloatT>
    typename std::enable_if<std::is_floating_point<FloatT>::value, FloatT>::type
    Poisson(FloatT mean);

    template <typename ProbRangeT, typename IndexType = uint32_t>
    IndexType Discrete(const ProbRangeT& probabilities);

    template <typename DistributionT>
    typename DistributionT::result_type FromDistribution(DistributionT& dist = DistributionT());

    template <typename DistributionT, typename VectorT, typename MatrixT>
    VectorT FromMultiVariateDistribution(DistributionT& dist, const VectorT& mean, const MatrixT& cholesky);

    template <typename DistributionT, typename VectorT>
    VectorT FromMultiVariateDistribution(DistributionT& dist, const VectorT& mean, const VectorT& sigma);

public:
    static constexpr result_type min() { return engine_type::min(); }
    static constexpr result_type max() { return engine_type::max(); }

    result_type operator()();

private:
    engine_type fEngine;
};

template <typename EngineT>
void RandomPrototype<EngineT>::Seed(result_type seed)
{
    sSeed = (seed == 0) ? std::random_device()() : seed;
}

template <typename EngineT>
std::atomic<typename EngineT::result_type> RandomPrototype<EngineT>::sSeed( 0 );



// TODO: Alright, Xcode as of version 7.3 does not support thread_local.
// At least, non-POD types are not supported.
// Possible workarounds: use mutexes or boost::thread_specific_ptr.
// For now I use a mutex lock, if thread_local is not supported by the compiler.

#ifdef NO_TLS

#include <mutex>

template <typename EngineT>
inline RandomPrototype<EngineT>& RandomPrototype<EngineT>::Instance()
{
    static RandomPrototype sInstance;
    return sInstance;
}

template <typename EngineT>
inline auto RandomPrototype<EngineT>::operator()() -> result_type
{
    static std::mutex sMtx;
    std::lock_guard<std::mutex> lock(sMtx);

    return fEngine();
}

#else

template <typename EngineT>
inline RandomPrototype<EngineT>& RandomPrototype<EngineT>::Instance()
{
    static thread_local RandomPrototype sInstance;
    return sInstance;
}

template <typename EngineT>
inline auto RandomPrototype<EngineT>::operator()() -> result_type
{
    return fEngine();
}

#endif // NO_TLS


template <typename EngineT>
inline RandomPrototype<EngineT>::RandomPrototype() :
    fEngine()
{
    fEngine.seed( sSeed++ );
}

template <typename EngineT>
inline RandomPrototype<EngineT>::~RandomPrototype()
{ }

template <typename EngineT>
template <typename FloatT>
typename std::enable_if<std::is_floating_point<FloatT>::value, FloatT>::type
inline RandomPrototype<EngineT>::Uniform(FloatT min, FloatT max, bool minIncluded, bool maxIncluded)
{
    if (minIncluded) {
        if (!maxIncluded) {
            // pass
        }
        else {
            max = std::nextafter(max, (max > min) ? std::numeric_limits<FloatT>::max() : -std::numeric_limits<FloatT>::max());
        }
    }
    else {
        if (!maxIncluded) {
            min = std::nextafter(min, (max > min) ? std::numeric_limits<FloatT>::max() : -std::numeric_limits<FloatT>::max());
        }
        else {
            std::swap(min, max);
        }
    }
    return std::uniform_real_distribution<FloatT>(min, max)(*this);
}

template <typename EngineT>
template <typename FloatT>
typename std::enable_if<std::is_floating_point<FloatT>::value, FloatT>::type
inline RandomPrototype<EngineT>::Uniform(FloatT min, FloatT max)
{
    return std::uniform_real_distribution<FloatT>(min, max)(*this);
}

template <typename EngineT>
template <typename IntegerT>
typename std::enable_if<std::is_integral<IntegerT>::value, IntegerT>::type
inline RandomPrototype<EngineT>::Uniform(IntegerT inclMin, IntegerT inclMax)
{
    return std::uniform_int_distribution<IntegerT>(inclMin, inclMax)(*this);
}

template <typename EngineT>
template <typename FloatT>
inline bool RandomPrototype<EngineT>::Bool(FloatT probability)
{
    return std::uniform_real_distribution<FloatT>(0.0, 1.0)(*this) < probability;
}

template <typename EngineT>
template <typename FloatT>
inline FloatT RandomPrototype<EngineT>::Normal(FloatT mean, FloatT sigma)
{
    return std::normal_distribution<FloatT>(mean, sigma)(*this);
}

template <typename EngineT>
template <typename DistributionT>
inline typename DistributionT::result_type RandomPrototype<EngineT>::FromDistribution(DistributionT& dist)
{
    return dist(*this);
}

template <typename EngineT>
template <typename DistributionT, typename VectorT, typename MatrixT>
inline VectorT RandomPrototype<EngineT>::FromMultiVariateDistribution(DistributionT& dist, const VectorT& mean, const MatrixT& cholesky)
{
    LOG_DEFINE("vmcmc.random");
    LOG_ASSERT( mean.size() == cholesky.size1() );

    VectorT noise( mean.size() );

    for (size_t i = 0; i < noise.size(); ++i)
        noise[i] = dist(*this);

    return mean + ublas::prod(noise, cholesky);
}

template <typename EngineT>
template <typename DistributionT, typename VectorT>
inline VectorT RandomPrototype<EngineT>::FromMultiVariateDistribution(DistributionT& dist, const VectorT& mean, const VectorT& sigma)
{
    LOG_DEFINE("vmcmc.random");
    LOG_ASSERT( mean.size() == sigma.size1() );

    VectorT noise( mean.size() );

    for (size_t i = 0; i < noise.size(); ++i)
        noise[i] = dist(*this);

    return mean + ublas::element_prod(noise, sigma);
}


template <typename EngineT>
template <typename IntegerT>
typename std::enable_if<std::is_integral<IntegerT>::value, IntegerT>::type
inline RandomPrototype<EngineT>::Poisson(double mean)
{
    return std::poisson_distribution<IntegerT>(mean)(*this);
}

template <typename EngineT>
template <typename FloatT>
typename std::enable_if<std::is_floating_point<FloatT>::value, FloatT>::type
inline RandomPrototype<EngineT>::Poisson(FloatT mean)
{
    if (mean > std::numeric_limits<uint64_t>::max() / 2.0)
        return Normal<FloatT>(mean, sqrt(mean));
    else
        return (FloatT) std::poisson_distribution<uint64_t>(mean)(*this);
}

template <typename EngineT>
template <typename FloatT>
inline FloatT RandomPrototype<EngineT>::Exponential(FloatT tau)
{
    return std::exponential_distribution<FloatT>(1.0/tau)(*this);
}

template <typename EngineT>
template <typename ProbRangeT, typename IndexType>
inline IndexType RandomPrototype<EngineT>::Discrete(const ProbRangeT& probabilities)
{
    return std::discrete_distribution<IndexType>(probabilities.begin(), probabilities.end())(*this);
}


using Random = RandomPrototype<std::mt19937>;

} /* namespace vmcmc */

#endif /* VMCMC_RANDOM_H_ */