@@ -70,9 +70,8 @@ protected:

     void acceptExchange(AtomicProposal *prop, float d1, unsigned r1,

         unsigned c1, unsigned r2, unsigned c2);

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

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

-        GapsRng *rng);

+    OptionalFloat gibbsMass(AlphaParameters alpha, GapsRng *rng);

+    OptionalFloat gibbsMass(AlphaParameters alpha, float m1, float m2, GapsRng *rng);

 public:

@@ -346,7 +345,7 @@ void GibbsSampler<T, MatA, MatB>::birth(AtomicProposal *prop)

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

     {

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

-        mass = gibbsMass(alpha, &(prop->rng)).first;

+        mass = gibbsMass(alpha, &(prop->rng)).value; // 0 if it fails

     }

     else

     {

@@ -386,10 +385,10 @@ void GibbsSampler<T, MatA, MatB>::death(AtomicProposal *prop)

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

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

     {

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

-        if (gMass.second)

+        OptionalFloat gMass = gibbsMass(alpha, &(prop->rng));

+        if (gMass.hasValue)

         {

-            rebirthMass = gMass.first;

+            rebirthMass = gMass.value;

         }

     }

@@ -454,16 +453,15 @@ void GibbsSampler<T, MatA, MatB>::exchange(AtomicProposal *prop)

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

     {

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

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

-        if (gMass.second)

+        OptionalFloat gMass = gibbsMass(alpha, m1, m2, &(prop->rng));

+        if (gMass.hasValue)

         {

-            acceptExchange(prop, gMass.first, r1, c1, r2, c2);

+            acceptExchange(prop, gMass.value, r1, c1, r2, c2);

             return;

         }

     }

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

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

+    float newMass = prop->rng.truncGammaUpper(m1 + m2, 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);

@@ -537,7 +535,7 @@ float d1, unsigned r1, unsigned c1, unsigned r2, unsigned c2)

 }

 template <class T, class MatA, class MatB>

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

+OptionalFloat GibbsSampler<T, MatA, MatB>::gibbsMass(AlphaParameters alpha,

 GapsRng *rng)

 {        

     alpha.s *= mAnnealingTemp;

@@ -553,14 +551,17 @@ GapsRng *rng)

         {

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

+            if (gMass >= gaps::epsilon)

+            {

+                return OptionalFloat(gMass);

+            }

         }

     }

-    return std::pair<float, bool>(0.f, false);

+    return OptionalFloat();

 }

 template <class T, class MatA, class MatB>

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

+OptionalFloat GibbsSampler<T, MatA, MatB>::gibbsMass(AlphaParameters alpha,

 float m1, float m2, GapsRng *rng)

 {

     alpha.s *= mAnnealingTemp;

@@ -576,11 +577,11 @@ float m1, float m2, GapsRng *rng)

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

         {

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

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

+            return OptionalFloat(gMass);

         }

     }

-    return std::pair<float, bool>(0.f, false);

+    return OptionalFloat();

 }

 #ifdef GAPS_DEBUG

@@ -8,8 +8,8 @@

 namespace gaps

 {

     const float epsilon = 1.0e-10f;

-    const float pi = 3.14159265358979323846264f;

-    const float pi_double = 3.14159265358979323846264;

+    const float pi = 3.1415926535897932384626433832795f;

+    const float pi_double = 3.1415926535897932384626433832795;

     float min(float a, float b);

     unsigned min(unsigned a, unsigned b);

@@ -237,14 +237,16 @@ float GapsRng::inverseNormSample(float a, float b, float mean, float sd)

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

 }

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

+float GapsRng::truncGammaUpper(float b, float shape, float scale)

 {

-    float u = uniform(a, b);

+    float upper = gaps::p_gamma(b, shape, scale);

+

+    float u = uniform(0.f, upper);

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

     {

-        u = uniform(a, b);

+        u = uniform(0.f, upper);

     }

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

+    return gaps::q_gamma(u, shape, scale);

 }

 Archive& operator<<(Archive &ar, GapsRng &gen)

@@ -303,3 +305,9 @@ double gaps::lgamma(double x)

 {

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

 }

+

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

+{

+    return std::exp((d - mean) * (d - mean) / (-2.f * sd * sd))

+        / std::sqrt(2.f * gaps::pi * sd * sd);

+}

@@ -9,8 +9,8 @@

 struct OptionalFloat

 {

-    bool hasValue;

     float value;

+    bool hasValue;

     OptionalFloat() : hasValue(false), value(0.f) {}

     OptionalFloat(float f) : hasValue(true), value(f) {}

@@ -26,6 +26,13 @@ namespace gaps

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

+

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

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

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

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

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

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

 }

 // used for seeding individual rngs

@@ -66,7 +73,7 @@ public:

     float exponential(float lambda);

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

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

+    float truncGammaUpper(float b, float shape, float scale);

     static void setSeed(uint64_t seed);

     static void save(Archive &ar);