@@ -538,24 +538,12 @@ template <class T, class MatA, class MatB>

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

 GapsRng *rng)

 {        

-    alpha.s *= mAnnealingTemp;

-    alpha.su *= mAnnealingTemp;

-

+    alpha *= mAnnealingTemp;

     if (alpha.s > gaps::epsilon)

     {

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

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

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

-

-        if (pLower < 1.f)

-        {

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

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

-            if (gMass >= gaps::epsilon)

-            {

-                return OptionalFloat(gMass);

-            }

-        }

+        return rng->truncNormal(gaps::epsilon, mMaxGibbsMass / mLambda, mean, sd);

     }

     return OptionalFloat();

 }

@@ -564,22 +552,12 @@ template <class T, class MatA, class MatB>

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

 float m1, float m2, GapsRng *rng)

 {

-    alpha.s *= mAnnealingTemp;

-    alpha.su *= mAnnealingTemp;

-

+    alpha *= mAnnealingTemp;

     if (alpha.s > gaps::epsilon)

     {

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

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

-        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 = rng->inverseNormSample(pLower, pUpper, mean, sd);

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

-            return OptionalFloat(gMass);

-        }

+        return rng->truncNormal(-m1, m2, mean, sd);

     }

     return OptionalFloat();

 }

@@ -19,9 +19,15 @@ struct AlphaParameters

     AlphaParameters operator+(const AlphaParameters &other) const

     {

         float rs = s + other.s;

-        float rsu = su - other.su; // weird

+        float rsu = su - other.su; // not a typo

         return AlphaParameters(rs, rsu);

     }

+

+    void operator*=(float v)

+    {

+        s *= v;

+        su *= v;

+    }

 };

 namespace gaps

@@ -227,14 +227,24 @@ float GapsRng::exponential(float lambda)

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

 }

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

+OptionalFloat GapsRng::truncNormal(float a, float b, float mean, float sd)

 {

-    float u = uniform(a, b);

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

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

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

+

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

     {

-        u = uniform(a, b);

+        float u = uniform(pLower, pUpper);

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

+        {

+            u = uniform(pLower, pUpper);

+        }

+        float ret = gaps::q_norm(u, mean, sd);

+        GAPS_ASSERT(ret >= a);

+        GAPS_ASSERT(ret <= b);

+        return OptionalFloat(ret);

     }

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

+    OptionalFloat();

 }

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

@@ -72,7 +72,7 @@ public:

     int poisson(double lambda);

     float exponential(float lambda);

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

+    OptionalFloat truncNormal(float a, float b, float mean, float sd);

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

     static void setSeed(uint64_t seed);