@@ -147,23 +147,36 @@ void AsynchronousGibbsSampler<DataModel>::birth(const AtomicProposal &prop)

 template <class DataModel>

 void AsynchronousGibbsSampler<DataModel>::death(const AtomicProposal &prop)

 {

-    // try to do a rebirth in the place of this atom

+    // determine mass to attempt rebirth with

+    float rebirthMass = prop.atom1->mass(); // default rebirth mass == no change to atom

+    AlphaParameters alpha = DataModel::alphaParametersWithChange(prop.r1, prop.c1,

+        -1.f * prop.atom1->mass()) * DataModel::annealingTemp();

     if (DataModel::canUseGibbs(prop.c1))

     {

-        OptionalFloat mass = DataModel::sampleDeathAndRebirth(prop.r1, prop.c1,

-            -1.f * prop.atom1->mass(), &(prop.rng));

-        if (mass.hasValue() && mass.value() >= gaps::epsilon)

+        OptionalFloat gMass = gibbsMass(alpha, 0.f, DataModel::maxGibbsMass(), &(prop.rng),

+            DataModel::lambda());

+        if (gMass.hasValue())

         {

-            mQueue.rejectDeath();

-            DataModel::safelyChangeMatrix(prop.r1, prop.c1, mass.value() - prop.atom1->mass());

-            prop.atom1->updateMass(mass.value());

-            return;

+            rebirthMass = gMass.value();

         }

     }

-    // if rebirth fails, then kill off atom

-    mQueue.acceptDeath();

-    DataModel::safelyChangeMatrix(prop.r1, prop.c1, -1.f * prop.atom1->mass());

-    mDomain.cacheErase(prop.atom1);

+    // handle accept/reject of the rebirth

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

+    if (std::log(prop.rng.uniform()) < deltaLL) // accept

+    {

+        mQueue.rejectDeath();

+        if (rebirthMass != prop.atom1->mass())

+        {

+            DataModel::safelyChangeMatrix(prop.r1, prop.c1, rebirthMass - prop.atom1->mass());

+            prop.atom1->updateMass(rebirthMass);

+        }

+    }

+    else // reject

+    {

+        mQueue.acceptDeath();

+        DataModel::safelyChangeMatrix(prop.r1, prop.c1, -1.f * prop.atom1->mass());

+        mDomain.cacheErase(prop.atom1);

+    }

 }

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

@@ -82,11 +82,21 @@ uint64_t DenseNormalModel::nPatterns() const

     return mMatrix.nCol();

 }

+float DenseNormalModel::annealingTemp() const

+{

+    return mAnnealingTemp;

+}

+

 float DenseNormalModel::lambda() const

 {

     return mLambda;

 }

+float DenseNormalModel::maxGibbsMass() const

+{

+    return mMaxGibbsMass;

+}

+

 bool DenseNormalModel::canUseGibbs(unsigned col) const

 {

     return !gaps::isVectorZero(mOtherMatrix->getCol(col));

@@ -33,7 +33,9 @@ protected:

     friend class GapsStatistics;

     uint64_t nElements() const;

     uint64_t nPatterns() const;

+    float annealingTemp() const;

     float lambda() const;

+    float maxGibbsMass() const;

     bool canUseGibbs(unsigned col) const;

     bool canUseGibbs(unsigned c1, unsigned c2) const;

     void changeMatrix(unsigned row, unsigned col, float delta);

@@ -157,23 +157,34 @@ void SingleThreadedGibbsSampler<DataModel>::death()

     AtomType *atom = mDomain.randomAtom(&mRng);

     unsigned row = (atom->pos() / mBinLength) / mNumPatterns;

     unsigned col = (atom->pos() / mBinLength) % mNumPatterns;

-

-    // try to do a rebirth in the place of this atom

+    // determine mass to attempt rebirth with

+    float rebirthMass = atom->mass(); // default rebirth mass == no change to atom

+    AlphaParameters alpha = DataModel::alphaParametersWithChange(row, col, -1.f * atom->mass())

+        * DataModel::annealingTemp();

     if (DataModel::canUseGibbs(col))

     {

-        OptionalFloat mass = DataModel::sampleDeathAndRebirth(row, col,

-            -1.f * atom->mass(), &mRng);

-        if (mass.hasValue() && mass.value() >= gaps::epsilon)

+        OptionalFloat gMass = gibbsMass(alpha, 0.f, DataModel::maxGibbsMass(), &mRng,

+            DataModel::lambda());

+        if (gMass.hasValue())

         {

-            DataModel::safelyChangeMatrix(row, col, mass.value() - atom->mass());

-            atom->updateMass(mass.value());

-            return;

+            rebirthMass = gMass.value();

         }

     }

-

-    // if rebirth fails, then kill off atom

-    DataModel::safelyChangeMatrix(row, col, -1.f * atom->mass());

-    mDomain.erase(atom);

+    // handle accept/reject of the rebirth

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

+    if (std::log(mRng.uniform()) < deltaLL) // accept

+    {

+        if (rebirthMass != atom->mass())

+        {

+            DataModel::safelyChangeMatrix(row, col, rebirthMass - atom->mass());

+            atom->updateMass(rebirthMass);

+        }

+    }

+    else // reject

+    {

+        DataModel::safelyChangeMatrix(row, col, -1.f * atom->mass());

+        mDomain.erase(atom);

+    }

 }

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

@@ -74,11 +74,21 @@ uint64_t SparseNormalModel::nPatterns() const

     return mMatrix.nCol();

 }

+float SparseNormalModel::annealingTemp() const

+{

+    return mAnnealingTemp;

+}

+

 float SparseNormalModel::lambda() const

 {

     return mLambda;

 }

+float SparseNormalModel::maxGibbsMass() const

+{

+    return mMaxGibbsMass;

+}

+

 bool SparseNormalModel::canUseGibbs(unsigned col) const

 {

     return !gaps::isVectorZero(mOtherMatrix->getCol(col));

@@ -33,7 +33,9 @@ public:

 protected:

     uint64_t nElements() const;

     uint64_t nPatterns() const;

+    float annealingTemp() const;

     float lambda() const;

+    float maxGibbsMass() const;

     bool canUseGibbs(unsigned col) const;

     bool canUseGibbs(unsigned c1, unsigned c2) const;

     void changeMatrix(unsigned row, unsigned col, float delta);