src/math/Random.cpp
1c81468c
 // [[Rcpp::depends(BH)]]
53d667e0
 
b0eb9022
 #include "Random.h"
1d70f59e
 #include "../GapsAssert.h"
351b251a
 
0f279c94
 // TODO remove boost dependency
561abd98
 // TODO make concurrent random generation safe
 // spawn a random generator for each thread
 // - no way to acheive consistent results across different nThreads
0f279c94
 
2630fb70
 
351b251a
 #include <boost/math/distributions/exponential.hpp>
2630fb70
 #include <boost/math/distributions/gamma.hpp>
0733d393
 #include <boost/math/distributions/normal.hpp>
 
 #include <boost/random/exponential_distribution.hpp>
 #include <boost/random/mersenne_twister.hpp>
 #include <boost/random/normal_distribution.hpp>
 #include <boost/random/poisson_distribution.hpp>
 #include <boost/random/uniform_01.hpp>
 #include <boost/random/uniform_real_distribution.hpp>
3b232069
 
70ce1bd7
 #include <algorithm>
dd3a9441
 #include <stdint.h>
03f50947
 
d465f5e7
 #ifdef __GAPS_OPENMP__
03f50947
     #include <omp.h>
 #else
     typedef int omp_int_t;
     inline omp_int_t omp_get_thread_num() { return 0; }
     inline omp_int_t omp_get_max_threads() { return 1; }
 #endif
53d667e0
 
2cfdfcd8
 #define Q_GAMMA_THRESHOLD 0.000001f
c37dd577
 #define Q_GAMMA_MIN_VALUE 0.0
 
0f279c94
 //typedef boost::random::mt19937 RNGType;
 typedef boost::random::mt11213b RNGType; // should be faster
3b232069
 
0733d393
 static std::vector<RNGType>& rng()
 {
     static std::vector<RNGType> allRngs(omp_get_max_threads());
     return allRngs;
 }
1c81468c
 
3b232069
 void gaps::random::save(Archive &ar)
ae66804f
 {
8269c03c
     for (unsigned i = 0; i < rng().size(); ++i)
     {
         ar << rng().at(i);
     }
ae66804f
 }
 
3b232069
 void gaps::random::load(Archive &ar)
ae66804f
 {
8269c03c
     for (unsigned i = 0; i < rng().size(); ++i)
     {
         ar >> rng().at(i);
     }
ae66804f
 }
 
075fd5dd
 void gaps::random::setSeed(uint32_t seed)
1c81468c
 {
e5b7277b
     unsigned n = omp_get_max_threads();
 
0733d393
     rng().at(0).seed(seed);
e5b7277b
 
     boost::random::uniform_int_distribution<uint32_t> seedDist(0,
         std::numeric_limits<uint32_t>::max());
 
     for (unsigned i = 1; i < n; ++i)
     {
0733d393
         uint32_t newSeed = seedDist(rng().at(i-1));
         rng().at(i).seed(newSeed);
e5b7277b
     }
1c81468c
 }
 
7b8ff090
 float gaps::random::normal(float mean, float var)
1c81468c
 {
7b8ff090
     boost::random::normal_distribution<float> dist(mean, var);
0733d393
     return dist(rng().at(omp_get_thread_num()));
1c81468c
 }
 
7b8ff090
 int gaps::random::poisson(float lambda)
1c81468c
 {
4968e606
     boost::random::poisson_distribution<> dist(lambda);
0733d393
     return dist(rng().at(omp_get_thread_num()));
1c81468c
 }
 
7b8ff090
 float gaps::random::exponential(float lambda)
1c81468c
 {
4968e606
     boost::random::exponential_distribution<> dist(lambda);
0733d393
     return dist(rng().at(omp_get_thread_num()));
1c81468c
 }
 
7b8ff090
 float gaps::random::uniform()
1c81468c
 {
0733d393
     boost::random::uniform_01<RNGType&> u01_dist(rng().at(omp_get_thread_num()));
e5b7277b
     return u01_dist();
1c81468c
 }
090d993c
 
7b8ff090
 float gaps::random::uniform(float a, float b)
075fd5dd
 {
59870cd4
     if (a == b)
     {
         return a;
     }
0733d393
     boost::random::uniform_real_distribution<> dist(a,b);
     return dist(rng().at(omp_get_thread_num()));
59870cd4
 }
 
7475fc95
 uint64_t gaps::random::uniform64()
 {
     boost::random::uniform_int_distribution<uint64_t> dist(0,
         std::numeric_limits<uint64_t>::max());
0733d393
     return dist(rng().at(omp_get_thread_num()));
7475fc95
 }
 
 uint64_t gaps::random::uniform64(uint64_t a, uint64_t b)
59870cd4
 {
     if (a == b)
     {
         return a;
     }
0733d393
     boost::random::uniform_int_distribution<uint64_t> dist(a,b);
     return dist(rng().at(omp_get_thread_num()));
075fd5dd
 }
 
7b8ff090
 float gaps::random::d_gamma(float d, float shape, float scale)
c37dd577
 {
     boost::math::gamma_distribution<> gam(shape, scale);
     return pdf(gam, d);
 }
 
7b8ff090
 float gaps::random::p_gamma(float p, float shape, float scale)
c37dd577
 {
     boost::math::gamma_distribution<> gam(shape, scale);
     return cdf(gam, p);
 }
 
7b8ff090
 float gaps::random::q_gamma(float q, float shape, float scale)
c37dd577
 {
     if (q < Q_GAMMA_THRESHOLD)
     {
         return Q_GAMMA_MIN_VALUE;
     }
0733d393
     boost::math::gamma_distribution<> gam(shape, scale);
     return quantile(gam, q);
c37dd577
 }
 
7b8ff090
 float gaps::random::d_norm(float d, float mean, float sd)
c37dd577
 {
     boost::math::normal_distribution<> norm(mean, sd);
     return pdf(norm, d);
 }
 
7b8ff090
 float gaps::random::q_norm(float q, float mean, float sd)
c37dd577
 {
     boost::math::normal_distribution<> norm(mean, sd);
     return quantile(norm, q);
 }
 
7b8ff090
 float gaps::random::p_norm(float p, float mean, float sd)
c37dd577
 {
     boost::math::normal_distribution<> norm(mean, sd);
     return cdf(norm, p);
 }
49a5b154
 
 float gaps::random::inverseNormSample(float a, float b, float mean, float sd)
 {
     float u = gaps::random::uniform(a, b);
     while (u == 0.f || u == 1.f)
     {
         u = gaps::random::uniform(a, b);
     }
     return gaps::random::q_norm(u, mean, sd);
 }
 
 float gaps::random::inverseGammaSample(float a, float b, float mean, float sd)
 {
     float u = gaps::random::uniform(a, b);
     while (u == 0.f || u == 1.f)
     {
         u = gaps::random::uniform(a, b);
     }
     return gaps::random::q_gamma(u, mean, sd);
4e746870
 }
 
70ce1bd7
 std::vector<unsigned> sample(std::vector<unsigned> &elements, unsigned n) {
     std::vector<unsigned> sampleVect;
     for (unsigned i = 0; i < n; ++i)
     {
f082ea7a
         GAPS_ASSERT(n <= elements.size());
70ce1bd7
         unsigned sampleIndex = gaps::random::uniform64(0, elements.size());
         sampleVect.push_back(elements.at(sampleIndex));
f082ea7a
 
         elements[sampleIndex] = elements.back();
         elements.pop_back();
70ce1bd7
     }
     return sampleVect;
 
     /*
4e746870
     std::vector<unsigned> sampleVect;
     std::vector<unsigned> sampledIndices;
     for (unsigned i = 0; i < n; ++i)
     {
         while(true)
         {
             unsigned sampleIndex = gaps::random::uniform64(0, elements.size());
             if (find(sampledIndices.begin(), sampledIndices.end(), sampleIndex) == sampledIndices.end())
             {
                 sampleVect.push_back(elements.at(sampleIndex));
                 sampledIndices.push_back(sampleIndex);
                 break;
             }
         }
     }
     return sampleVect;
70ce1bd7
     */
4e746870
 }