@@ -14,7 +14,7 @@ Atom AtomicDomain::front() const

 // O(1)

 Atom AtomicDomain::randomAtom() const

 {

-    uint64_t ndx = gaps::random::uniform64(0, mAtoms.size() - 1);

+    uint64_t ndx = mRng.uniform64(0, mAtoms.size() - 1);

     return mAtoms[ndx];

 }

@@ -24,7 +24,7 @@ uint64_t AtomicDomain::randomFreePosition() const

     uint64_t pos = 0;

     do

     {

-        pos = gaps::random::uniform64(0, mDomainSize);

+        pos = mRng.uniform64(0, mDomainSize);

     } while (mPositionLookup.count(pos) > 0);

     return pos;

 }

@@ -2,6 +2,7 @@

 #define __COGAPS_ATOMIC_DOMAIN_H__

 #include "Archive.h"

+#include "math/Random.h"

 #include <boost/unordered_map.hpp>

@@ -15,6 +16,7 @@ class AtomicDomain;

 // Want the neighbors to be pointers, but can't use pointers since vector

 // is resized, shifted integers allow for 0 to be "null" and 1 to be the

 // first index. AtomicDomain is responsible for managing this correctly.

+// TODO use get/set neighbor

 struct Atom

 {

 private:    

@@ -75,6 +77,8 @@ private:

     mutable unsigned mInsertCacheIndex;

     mutable unsigned mEraseCacheIndex;

+    mutable GapsRng mRng;

+

     Atom& _left(const Atom &atom);

     Atom& _right(const Atom &atom);

@@ -61,7 +61,7 @@ unsigned getNumPatterns(const Rcpp::List &allParams)

     {

         std::string file(Rcpp::as<std::string>(allParams["checkpointInFile"]));

         Archive ar(file, ARCHIVE_READ);

-        gaps::random::load(ar);

+        GapsRng::load(ar);

         ar >> nPatterns;

         ar.close();

     }

@@ -119,7 +119,7 @@ const Rcpp::Nullable<Rcpp::NumericMatrix> &fixedMatrix, bool isMaster)

     {

         std::string file(Rcpp::as<std::string>(allParams["checkpointInFile"]));

         Archive ar(file, ARCHIVE_READ);

-        gaps::random::load(ar);

+        GapsRng::load(ar);

         ar >> runner;

         ar.close();

     }

@@ -135,7 +135,7 @@ const Rcpp::Nullable<Rcpp::NumericMatrix> &fixedMatrix, bool isMaster)

         // set parameters that would be saved in the checkpoint

         const Rcpp::S4 &gapsParams(allParams["gaps"]);

-        gaps::random::setSeed(gapsParams.slot("seed"));

+        GapsRng::setSeed(gapsParams.slot("seed"));

         runner.recordSeed(gapsParams.slot("seed"));

         runner.setMaxIterations(gapsParams.slot("nIterations"));

         runner.setSparsity(gapsParams.slot("alphaA"),

@@ -137,8 +137,8 @@ void GapsRunner::runOnePhase()

         }

         // number of updates per iteration is poisson 

-        unsigned nA = gaps::random::poisson(gaps::max(mASampler.nAtoms(), 10));

-        unsigned nP = gaps::random::poisson(gaps::max(mPSampler.nAtoms(), 10));

+        unsigned nA = mRng.poisson(gaps::max(mASampler.nAtoms(), 10));

+        unsigned nP = mRng.poisson(gaps::max(mPSampler.nAtoms(), 10));

         updateSampler(nA, nP);

         if (mPhase == 'S')

@@ -245,7 +245,7 @@ void GapsRunner::createCheckpoint()

         // create checkpoint file

         Archive ar(mCheckpointOutFile, ARCHIVE_WRITE);

-        gaps::random::save(ar);

+        GapsRng::save(ar);

         ar << mNumPatterns << mSeed << mASampler << mPSampler << mStatistics

             << mFixedMatrix << mMaxIterations << mPhase << mCurrentIteration

             << mNumUpdatesA << mNumUpdatesP;

@@ -97,6 +97,8 @@ private:

     unsigned mNumUpdatesA;

     unsigned mNumUpdatesP;

+

+    mutable GapsRng mRng;

     void runOnePhase();

     void updateSampler(unsigned nA, unsigned nP);

@@ -65,19 +65,20 @@ protected:

     void addMass(uint64_t pos, float mass, unsigned row, unsigned col);

     void removeMass(uint64_t pos, float mass, unsigned row, unsigned col);

-    void birth(uint64_t pos, unsigned row, unsigned col);

-    void death(uint64_t pos, float mass, unsigned row, unsigned col);

+    void birth(uint64_t pos, unsigned row, unsigned col, GapsRng *rng);

+    void death(uint64_t pos, float mass, unsigned row, unsigned col, GapsRng *rng);

     void move(uint64_t src, float mass, uint64_t dest, unsigned r1, unsigned c1,

-        unsigned r2, unsigned c2);

+        unsigned r2, unsigned c2, GapsRng *rng);

     void exchange(uint64_t p1, float m1, uint64_t p2, float m2, unsigned r1,

-        unsigned c1, unsigned r2, unsigned c2);

+        unsigned c1, unsigned r2, unsigned c2, GapsRng *rng);

     bool updateAtomMass(uint64_t pos, float mass, float delta);

     void acceptExchange(uint64_t p1, float m1, float d1, uint64_t p2, float m2,

         float d2, unsigned r1, unsigned c1, unsigned r2, unsigned c2);

-    std::pair<float, bool> gibbsMass(AlphaParameters alpha);

-    std::pair<float, bool> gibbsMass(AlphaParameters alpha, float m1, float m2);

+    std::pair<float, bool> gibbsMass(AlphaParameters alpha, GapsRng *rng);

+    std::pair<float, bool> gibbsMass(AlphaParameters alpha, float m1, float m2,

+        GapsRng *rng);

 public:

@@ -399,16 +400,16 @@ void GibbsSampler<T, MatA, MatB>::processProposal(const AtomicProposal &prop)

     switch (prop.type)

     {

         case 'B':

-            birth(prop.pos1, r1, c1);

+            birth(prop.pos1, r1, c1, &prop.rng);

             break;

         case 'D':

-            death(prop.pos1, prop.mass1, r1, c1);

+            death(prop.pos1, prop.mass1, r1, c1, &prop.rng);

             break;

         case 'M':

-            move(prop.pos1, prop.mass1, prop.pos2, r1, c1, r2, c2);

+            move(prop.pos1, prop.mass1, prop.pos2, r1, c1, r2, c2, &prop.rng);

             break;

         case 'E':

-            exchange(prop.pos1, prop.mass1, prop.pos2, prop.mass2, r1, c1, r2, c2);

+            exchange(prop.pos1, prop.mass1, prop.pos2, prop.mass2, r1, c1, r2, c2, &prop.rng);

             break;

     }

 }

@@ -433,18 +434,18 @@ void GibbsSampler<T, MatA, MatB>::removeMass(uint64_t pos, float mass, unsigned

 // or with the gibbs mass calculation

 template <class T, class MatA, class MatB>

 void GibbsSampler<T, MatA, MatB>::birth(uint64_t pos, unsigned row,

-unsigned col)

+unsigned col, GapsRng *rng)

 {

     // calculate proposed mass

     float mass = 0.f;

     if (impl()->canUseGibbs(row, col))

     {

         AlphaParameters alpha = impl()->alphaParameters(row, col);

-        mass = gibbsMass(alpha).first;

+        mass = gibbsMass(alpha, rng).first;

     }

     else

     {

-        mass = gaps::random::exponential(mLambda);

+        mass = rng->exponential(mLambda);

     }

     // accept mass as long as it's non-zero

@@ -463,7 +464,7 @@ unsigned col)

 // the original mass or the gibbs mass calculation

 template <class T, class MatA, class MatB>

 void GibbsSampler<T, MatA, MatB>::death(uint64_t pos, float mass, unsigned row,

-unsigned col)

+unsigned col, GapsRng *rng)

 {

     // kill off atom

     mMatrix(row, col) = gaps::max(mMatrix(row, col) - mass, 0.f);

@@ -474,7 +475,7 @@ unsigned col)

     AlphaParameters alpha = impl()->alphaParameters(row, col);

     if (impl()->canUseGibbs(row, col))

     {

-        std::pair<float, bool> gMass = gibbsMass(alpha);

+        std::pair<float, bool> gMass = gibbsMass(alpha, rng);

         if (gMass.second)

         {

             rebirthMass = gMass.first;

@@ -483,7 +484,7 @@ unsigned col)

     // accept/reject rebirth

     float deltaLL = rebirthMass * (alpha.su - alpha.s * rebirthMass / 2.f);

-    if (deltaLL * mAnnealingTemp >= std::log(gaps::random::uniform()))

+    if (deltaLL * mAnnealingTemp >= std::log(rng->uniform()))

     {

         mDomain.updateMass(pos, rebirthMass);

         mMatrix(row, col) += rebirthMass;

@@ -500,12 +501,12 @@ unsigned col)

 // move mass from src to dest in the atomic domain

 template <class T, class MatA, class MatB>

 void GibbsSampler<T, MatA, MatB>::move(uint64_t src, float mass, uint64_t dest,

-unsigned r1, unsigned c1, unsigned r2, unsigned c2)

+unsigned r1, unsigned c1, unsigned r2, unsigned c2, GapsRng *rng)

 {

     if (r1 != r2 || c1 != c2) // automatically reject if change in same bin

     {

         float deltaLL = impl()->computeDeltaLL(r1, c1, -mass, r2, c2, mass);

-        if (deltaLL * mAnnealingTemp > std::log(gaps::random::uniform()))

+        if (deltaLL * mAnnealingTemp > std::log(rng->uniform()))

         {

             removeMass(src, mass, r1, c1);

             addMass(dest, mass, r2, c2);

@@ -518,14 +519,14 @@ unsigned r1, unsigned c1, unsigned r2, unsigned c2)

 // the exchange

 template <class T, class MatA, class MatB>

 void GibbsSampler<T, MatA, MatB>::exchange(uint64_t p1, float m1, uint64_t p2,

-float m2, unsigned r1, unsigned c1, unsigned r2, unsigned c2)

+float m2, unsigned r1, unsigned c1, unsigned r2, unsigned c2, GapsRng *rng)

 {

     if (r1 != r2 || c1 != c2) // automatically reject if change in same bin

     {

         if (impl()->canUseGibbs(r1, c1, r2, c2))

         {

             AlphaParameters alpha = impl()->alphaParameters(r1, c1, r2, c2);

-            std::pair<float, bool> gMass = gibbsMass(alpha, m1, m2);

+            std::pair<float, bool> gMass = gibbsMass(alpha, m1, m2, rng);

             if (gMass.second)

             {

                 acceptExchange(p1, m1, gMass.first, p2, m2, -gMass.first, r1,

@@ -534,14 +535,14 @@ float m2, unsigned r1, unsigned c1, unsigned r2, unsigned c2)

             }

         }

-        float pUpper = gaps::random::p_gamma(m1 + m2, 2.f, 1.f / mLambda);

-        float newMass = gaps::random::inverseGammaSample(0.f, pUpper, 2.f, 1.f / mLambda);

+        float pUpper = gaps::p_gamma(m1 + m2, 2.f, 1.f / mLambda);

+        float newMass = rng->inverseGammaSample(0.f, pUpper, 2.f, 1.f / mLambda);

         float delta = m1 > m2 ? newMass - m1 : m2 - newMass; // change larger mass

         float pOldMass = 2.f * newMass > m1 + m2 ? gaps::max(m1, m2) : gaps::min(m1, m2);

-        float pNew = gaps::random::d_gamma(newMass, 2.f, 1.f / mLambda);

-        float pOld = gaps::random::d_gamma(pOldMass, 2.f, 1.f / mLambda);

+        float pNew = gaps::d_gamma(newMass, 2.f, 1.f / mLambda);

+        float pOld = gaps::d_gamma(pOldMass, 2.f, 1.f / mLambda);

         if (pOld == 0.f && pNew != 0.f) // special case

         {

@@ -550,7 +551,7 @@ float m2, unsigned r1, unsigned c1, unsigned r2, unsigned c2)

         }

         float deltaLL = impl()->computeDeltaLL(r1, c1, delta, r2, c2, -delta);

         float priorLL = (pOld == 0.f) ? 1.f : pOld / pNew;

-        float u = std::log(gaps::random::uniform() * priorLL);

+        float u = std::log(rng->uniform() * priorLL);

         if (u < deltaLL * mAnnealingTemp)

         {

             acceptExchange(p1, m1, delta, p2, m2, -delta, r1, c1, r2, c2);

@@ -605,7 +606,8 @@ unsigned r2, unsigned c2)

 }

 template <class T, class MatA, class MatB>

-std::pair<float, bool> GibbsSampler<T, MatA, MatB>::gibbsMass(AlphaParameters alpha)

+std::pair<float, bool> GibbsSampler<T, MatA, MatB>::gibbsMass(AlphaParameters alpha,

+GapsRng *rng)

 {        

     alpha.s *= mAnnealingTemp;

     alpha.su *= mAnnealingTemp;

@@ -614,11 +616,11 @@ std::pair<float, bool> GibbsSampler<T, MatA, MatB>::gibbsMass(AlphaParameters al

     {

         float mean = (alpha.su - mLambda) / alpha.s;

         float sd = 1.f / std::sqrt(alpha.s);

-        float pLower = gaps::random::p_norm(0.f, mean, sd);

+        float pLower = gaps::p_norm(0.f, mean, sd);

         if (pLower < 1.f)

         {

-            float m = gaps::random::inverseNormSample(pLower, 1.f, mean, sd);

+            float m = rng->inverseNormSample(pLower, 1.f, mean, sd);

             float gMass = gaps::min(m, mMaxGibbsMass / mLambda);

             return std::pair<float, bool>(gMass, gMass >= gaps::epsilon);

         }

@@ -628,7 +630,7 @@ std::pair<float, bool> GibbsSampler<T, MatA, MatB>::gibbsMass(AlphaParameters al

 template <class T, class MatA, class MatB>

 std::pair<float, bool> GibbsSampler<T, MatA, MatB>::gibbsMass(AlphaParameters alpha,

-float m1, float m2)

+float m1, float m2, GapsRng *rng)

 {

     alpha.s *= mAnnealingTemp;

     alpha.su *= mAnnealingTemp;

@@ -637,12 +639,12 @@ float m1, float m2)

     {

         float mean = alpha.su / alpha.s; // lambda cancels out

         float sd = 1.f / std::sqrt(alpha.s);

-        float pLower = gaps::random::p_norm(-m1, mean, sd);

-        float pUpper = gaps::random::p_norm(m2, mean, sd);

+        float pLower = gaps::p_norm(-m1, mean, sd);

+        float pUpper = gaps::p_norm(m2, mean, sd);

         if (!(pLower >  0.95f || pUpper < 0.05f))

         {

-            float delta = gaps::random::inverseNormSample(pLower, pUpper, mean, sd);

+            float delta = rng->inverseNormSample(pLower, pUpper, mean, sd);

             float gibbsMass = gaps::min(gaps::max(-m1, delta), m2); // conserve mass

             return std::pair<float, bool>(gibbsMass, true);

         }

@@ -103,8 +103,8 @@ bool ProposalQueue::makeProposal(AtomicDomain &domain)

     float bdProb = mMaxAtoms < 2 ? 0.6667f : 0.5f;

-    mU1 = mUseCachedRng ? mU1 : gaps::random::uniform();

-    mU2 = mUseCachedRng ? mU2: gaps::random::uniform();

+    mU1 = mUseCachedRng ? mU1 : mRng.uniform();

+    mU2 = mUseCachedRng ? mU2: mRng.uniform();

     mUseCachedRng = false;

     float lowerBound = deathProb(mMinAtoms);

@@ -166,7 +166,7 @@ bool ProposalQueue::move(AtomicDomain &domain)

         return false;

     }

-    uint64_t newLocation = gaps::random::uniform64(lbound, rbound - 1);

+    uint64_t newLocation = mRng.uniform64(lbound, rbound - 1);

     if (mUsedIndices.count(a.pos / mDimensionSize) || mUsedIndices.count(newLocation / mDimensionSize))

     {

         return false; // matrix conflict - can't compute deltaLL

@@ -4,6 +4,7 @@

 #include "Archive.h"

 #include "AtomicDomain.h"

 #include "data_structures/EfficientSets.h"

+#include "math/Random.h"

 #include <boost/unordered_set.hpp>

@@ -17,6 +18,7 @@ struct AtomicProposal

     float mass1;

     uint64_t pos2;

     float mass2;

+    mutable GapsRng rng;

     AtomicProposal(char t, uint64_t p1, float m1=0.f, uint64_t p2=0, float m2=0.f)

         : type(t), pos1(p1), mass1(m1), pos2(p2), mass2(m2)

@@ -45,6 +47,8 @@ private:

     float mU1;

     float mU2;

+    mutable GapsRng mRng;

+

     float deathProb(uint64_t nAtoms) const;

     bool birth(AtomicDomain &domain);

     bool death(AtomicDomain &domain);

@@ -4,6 +4,8 @@

 #define TEST_APPROX(x) Approx(x).epsilon(0.001)

+#if 0

+

 TEST_CASE("Test Random.h - Random Number Generation")

 {

     gaps::random::setSeed(0);

@@ -84,10 +86,12 @@ TEST_CASE("Test Random.h - Random Number Generation")

 TEST_CASE("Test Random.h - Distribution Calculations")

 {

-    REQUIRE(gaps::random::d_gamma(0.5f, 1.f, 1.f) == TEST_APPROX(0.607f));

-    REQUIRE(gaps::random::p_gamma(0.5f, 1.f, 1.f) == TEST_APPROX(0.394f));

-    REQUIRE(gaps::random::q_gamma(0.5f, 1.f, 1.f) == TEST_APPROX(0.693f));

-    REQUIRE(gaps::random::d_norm(0.5f, 0.f, 1.f) == TEST_APPROX(0.352f));

-    REQUIRE(gaps::random::q_norm(0.5f, 0.f, 1.f) == TEST_APPROX(0.000f));

-    REQUIRE(gaps::random::p_norm(0.5f, 0.f, 1.f) == TEST_APPROX(0.692f));

+    REQUIRE(gaps::d_gamma(0.5f, 1.f, 1.f) == TEST_APPROX(0.607f));

+    REQUIRE(gaps::p_gamma(0.5f, 1.f, 1.f) == TEST_APPROX(0.394f));

+    REQUIRE(gaps::q_gamma(0.5f, 1.f, 1.f) == TEST_APPROX(0.693f));

+    REQUIRE(gaps::d_norm(0.5f, 0.f, 1.f) == TEST_APPROX(0.352f));

+    REQUIRE(gaps::q_norm(0.5f, 0.f, 1.f) == TEST_APPROX(0.000f));

+    REQUIRE(gaps::p_norm(0.5f, 0.f, 1.f) == TEST_APPROX(0.692f));

 }

+

+#endif

\ No newline at end of file

@@ -4,6 +4,8 @@

 #include "../GibbsSampler.h"

 #include "../math/Random.h"

+#if 0

+

 TEST_CASE("Test Archive.h")

 {

     SECTION("Reading/Writing to an Archive")

@@ -164,4 +166,6 @@ TEST_CASE("Test Archive.h")

     // cleanup directory

     std::remove("test_ar.temp");

-}

\ No newline at end of file

+}

+

+#endif

\ No newline at end of file

@@ -20,149 +20,29 @@

 #include <algorithm>

 #include <stdint.h>

-#ifdef __GAPS_OPENMP__

-    #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

-

 #define Q_GAMMA_THRESHOLD 0.000001f

 #define Q_GAMMA_MIN_VALUE 0.0

-//typedef boost::random::mt19937 RNGType;

-//typedef boost::random::mt11213b RNGType; // should be faster

+const float maxU32AsFloat = static_cast<float>(std::numeric_limits<uint32_t>::max());

+const double maxU32AsDouble = static_cast<double>(std::numeric_limits<uint32_t>::max());

 static Xoroshiro128plus seeder;

-static std::vector<GapsRng>& rng()

+void GapsRng::save(Archive &ar)

 {

-    static std::vector<GapsRng> allRngs(omp_get_max_threads());

-    return allRngs;

+    ar << seeder;

 }

-void gaps::random::save(Archive &ar)

+void GapsRng::load(Archive &ar)

 {

-    for (unsigned i = 0; i < rng().size(); ++i)

-    {

-        ar << seeder;

-    }

+    ar >> seeder;

 }

-void gaps::random::load(Archive &ar)

-{

-    for (unsigned i = 0; i < rng().size(); ++i)

-    {

-        ar >> seeder;

-    }

-}

-

-void gaps::random::setSeed(uint32_t seed)

+void GapsRng::setSeed(uint64_t seed)

 {

     seeder.seed(seed);

 }

-int gaps::random::poisson(float lambda)

-{   

-    return rng().at(omp_get_thread_num()).poisson(lambda);

-}

-

-float gaps::random::exponential(float lambda)

-{

-    return rng().at(omp_get_thread_num()).exponential(lambda);

-}

-

-float gaps::random::uniform()

-{

-    return rng().at(omp_get_thread_num()).uniform();

-}

-

-float gaps::random::uniform(float a, float b)

-{

-    return rng().at(omp_get_thread_num()).uniform(a,b);

-}

-

-uint64_t gaps::random::uniform64()

-{

-    return rng().at(omp_get_thread_num()).uniform64();

-}

-

-uint64_t gaps::random::uniform64(uint64_t a, uint64_t b)

-{

-    return rng().at(omp_get_thread_num()).uniform64(a,b);

-}

-

-float gaps::random::d_gamma(float d, float shape, float scale)

-{

-    boost::math::gamma_distribution<> gam(shape, scale);

-    return pdf(gam, d);

-}

-

-float gaps::random::p_gamma(float p, float shape, float scale)

-{

-    boost::math::gamma_distribution<> gam(shape, scale);

-    return cdf(gam, p);

-}

-

-float gaps::random::q_gamma(float q, float shape, float scale)

-{

-    if (q < Q_GAMMA_THRESHOLD)

-    {

-        return Q_GAMMA_MIN_VALUE;

-    }

-    boost::math::gamma_distribution<> gam(shape, scale);

-    return quantile(gam, q);

-}

-

-float gaps::random::d_norm(float d, float mean, float sd)

-{

-    boost::math::normal_distribution<> norm(mean, sd);

-    return pdf(norm, d);

-}

-

-float gaps::random::q_norm(float q, float mean, float sd)

-{

-    boost::math::normal_distribution<> norm(mean, sd);

-    return quantile(norm, q);

-}

-

-float gaps::random::p_norm(float p, float mean, float sd)

-{

-    boost::math::normal_distribution<> norm(mean, sd);

-    return cdf(norm, p);

-}

-

-double gaps::lgamma(double x)

-{

-    return boost::math::lgamma(x);

-}

-

-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);

-}

-

-///////////////////////// NEW /////////////////////////////////////////////////

-

-const float maxU32AsFloat = static_cast<float>(std::numeric_limits<uint32_t>::max());

-const double maxU32AsDouble = static_cast<double>(std::numeric_limits<uint32_t>::max());

-

 static uint64_t rotl(const uint64_t x, int k)

 {

     return (x << k) | (x >> (64 - k));

@@ -352,3 +232,67 @@ float GapsRng::exponential(float lambda)

     return -1.f * std::log(uniform()) / lambda;

 }

+float GapsRng::inverseNormSample(float a, float b, float mean, float sd)

+{

+    float u = uniform(a, b);

+    while (u == 0.f || u == 1.f)

+    {

+        u = uniform(a, b);

+    }

+    return gaps::q_norm(u, mean, sd);

+}

+

+float GapsRng::inverseGammaSample(float a, float b, float mean, float sd)

+{

+    float u = uniform(a, b);

+    while (u == 0.f || u == 1.f)

+    {

+        u = uniform(a, b);

+    }

+    return gaps::q_gamma(u, mean, sd);

+}

+

+float gaps::d_gamma(float d, float shape, float scale)

+{

+    boost::math::gamma_distribution<> gam(shape, scale);

+    return pdf(gam, d);

+}

+

+float gaps::p_gamma(float p, float shape, float scale)

+{

+    boost::math::gamma_distribution<> gam(shape, scale);

+    return cdf(gam, p);

+}

+

+float gaps::q_gamma(float q, float shape, float scale)

+{

+    if (q < Q_GAMMA_THRESHOLD)

+    {

+        return Q_GAMMA_MIN_VALUE;

+    }

+    boost::math::gamma_distribution<> gam(shape, scale);

+    return quantile(gam, q);

+}

+

+float gaps::d_norm(float d, float mean, float sd)

+{

+    boost::math::normal_distribution<> norm(mean, sd);

+    return pdf(norm, d);

+}

+

+float gaps::q_norm(float q, float mean, float sd)

+{

+    boost::math::normal_distribution<> norm(mean, sd);

+    return quantile(norm, q);

+}

+

+float gaps::p_norm(float p, float mean, float sd)

+{

+    boost::math::normal_distribution<> norm(mean, sd);

+    return cdf(norm, p);

+}

+

+double gaps::lgamma(double x)

+{

+    return boost::math::lgamma(x);

+}

@@ -11,32 +11,13 @@ namespace gaps

 {

     double lgamma(double x);

-    namespace random

-    {

-        void setSeed(uint32_t seed);

-

-        float uniform();

-        float uniform(float a, float b);

-        uint64_t uniform64();

-        uint64_t uniform64(uint64_t a, uint64_t b);

-

-        float exponential(float lambda);

-        int poisson(float lambda);

-

-        float inverseNormSample(float a, float b, float mean, float sd);

-        float inverseGammaSample(float a, float b, float mean, float sd);

-

-        float d_gamma(float d, float shape, float scale);

-        float p_gamma(float p, float shape, float scale);

-        float q_gamma(float q, float shape, float scale);

-        float d_norm(float d, float mean, float sd);

-        float q_norm(float q, float mean, float sd);

-        float p_norm(float p, float mean, float sd);

-

-        void save(Archive &ar);

-        void load(Archive &ar);

-    } // namespace random

-} // namespace gaps

+    float d_gamma(float d, float shape, float scale);

+    float p_gamma(float p, float shape, float scale);

+    float q_gamma(float q, float shape, float scale);

+    float d_norm(float d, float mean, float sd);

+    float q_norm(float q, float mean, float sd);

+    float p_norm(float p, float mean, float sd);

+}

 // used for seeding individual rngs

 class Xoroshiro128plus

@@ -74,6 +55,13 @@ public:

     int poisson(double lambda);

     float exponential(float lambda);

+

+    float inverseNormSample(float a, float b, float mean, float sd);

+    float inverseGammaSample(float a, float b, float mean, float sd);

+

+    static void setSeed(uint64_t seed);

+    static void save(Archive &ar);

+    static void load(Archive &ar);

 private: