@@ -46,13 +46,13 @@ uint64_t AtomicDomain::size() const

 }

 // O(1)

-Atom& AtomicDomain::left(const Atom &atom)

+Atom& AtomicDomain::_left(const Atom &atom)

 {

     return mAtoms[atom.leftNdx - 1];

 }

 // O(1)

-Atom& AtomicDomain::right(const Atom &atom)

+Atom& AtomicDomain::_right(const Atom &atom)

 {

     return mAtoms[atom.rightNdx - 1];

 }

@@ -84,8 +84,6 @@ bool AtomicDomain::hasRight(const Atom &atom) const

 // O(logN)

 Atom AtomicDomain::insert(uint64_t pos, float mass)

 {

-    //omp_set_lock(&lock);

-

     // insert position into map

     std::map<uint64_t, uint64_t>::iterator iter, iterLeft, iterRight;

     iter = mAtomPositions.insert(std::pair<uint64_t, uint64_t>(pos, size())).first;

@@ -98,21 +96,19 @@ Atom AtomicDomain::insert(uint64_t pos, float mass)

     {

         --iterLeft;

         atom.leftNdx = iterLeft->second + 1;

-        left(atom).rightNdx = size() + 1;

+        _left(atom).rightNdx = size() + 1;

     }

     if (++iter != mAtomPositions.end())

     {

         ++iterRight;

         atom.rightNdx = iterRight->second + 1;

-        right(atom).leftNdx = size() + 1;

+        _right(atom).leftNdx = size() + 1;

     } 

     // add atom to vector

     mAtoms.push_back(atom);

     mUsedPositions.insert(pos);

-    //omp_unset_lock(&lock);

-

     return atom;

 }

@@ -121,19 +117,17 @@ Atom AtomicDomain::insert(uint64_t pos, float mass)

 // swap with last atom in vector, pop off the back

 void AtomicDomain::erase(uint64_t pos)

 {

-    //omp_set_lock(&lock);

-

     // get vector index of this atom and erase it

     uint64_t index = mAtomPositions.at(pos);

     // connect neighbors of atom to be deleted

     if (hasLeft(mAtoms[index]))

     {

-        left(mAtoms[index]).rightNdx = mAtoms[index].rightNdx;

+        _left(mAtoms[index]).rightNdx = mAtoms[index].rightNdx;

     }

     if (hasRight(mAtoms[index]))

     {

-        right(mAtoms[index]).leftNdx = mAtoms[index].leftNdx;

+        _right(mAtoms[index]).leftNdx = mAtoms[index].leftNdx;

     }

     // replace with atom from back

@@ -149,11 +143,11 @@ void AtomicDomain::erase(uint64_t pos)

         // update moved atom's neighbors

         if (hasLeft(mAtoms[index]))

         {

-            left(mAtoms[index]).rightNdx = index + 1;

+            _left(mAtoms[index]).rightNdx = index + 1;

         }

         if (hasRight(mAtoms[index]))

         {

-            right(mAtoms[index]).leftNdx = index + 1;

+            _right(mAtoms[index]).leftNdx = index + 1;

         }

     }

@@ -161,14 +155,58 @@ void AtomicDomain::erase(uint64_t pos)

     mAtomPositions.erase(pos);

     mAtoms.pop_back();

     mUsedPositions.erase(pos);

+}

+

+void AtomicDomain::cacheInsert(uint64_t pos, float mass) const

+{

+    unsigned ndx = 0;

+

+    #pragma omp critical(atomicInsert)

+    {

+        ndx = mInsertCacheIndex++;

+    }

+

+    mInsertCache[ndx] = RawAtom(pos, mass);

+}

+

+void AtomicDomain::cacheErase(uint64_t pos) const

+{

+    unsigned ndx = 0;

+

+    #pragma omp critical(atomicErase)

+    {

+        ndx = mEraseCacheIndex++;

+    }

+

+    mEraseCache[ndx] = pos;

+}

+

+void AtomicDomain::resetCache(unsigned n)

+{

+    mInsertCacheIndex = 0;

+    mEraseCacheIndex = 0;

+    mInsertCache.resize(n);

+    mEraseCache.resize(n);

+}

+

+void AtomicDomain::flushCache()

+{

+    for (unsigned i = 0; i < mEraseCacheIndex; ++i)

+    {

+        erase(mEraseCache[i]);

+    }

+

+    for (unsigned i = 0; i < mInsertCacheIndex; ++i)

+    {

+        insert(mInsertCache[i].pos, mInsertCache[i].mass);

+    }

-    //omp_unset_lock(&lock);

+    mInsertCache.clear();

+    mEraseCache.clear();

 }

 // O(logN)

 void AtomicDomain::updateMass(uint64_t pos, float newMass)

 {

-    //omp_set_lock(&lock);

     mAtoms[mAtomPositions.at(pos)].mass = newMass;

-    //omp_unset_lock(&lock);

 }

\ No newline at end of file

@@ -43,6 +43,15 @@ public:

     }

 };

+struct RawAtom

+{

+    uint64_t pos;

+    float mass;

+

+    RawAtom() : pos(0), mass(0.f) {}

+    RawAtom(uint64_t p, float m) : pos(p), mass(m) {}

+};

+

 // data structure that holds atoms

 class AtomicDomain

 {

@@ -58,6 +67,15 @@ private:

     // TODO google_dense_set - first profile and benchmark

     boost::unordered_set<uint64_t> mUsedPositions;

+    mutable std::vector<RawAtom> mInsertCache;

+    mutable std::vector<uint64_t> mEraseCache;

+

+    mutable unsigned mInsertCacheIndex;

+    mutable unsigned mEraseCacheIndex;

+

+    Atom& _left(const Atom &atom);

+    Atom& _right(const Atom &atom);

+

 public:

     AtomicDomain();

@@ -70,8 +88,6 @@ public:

     uint64_t randomFreePosition() const;

     uint64_t size() const;

-    Atom& left(const Atom &atom);

-    Atom& right(const Atom &atom);

     const Atom& left(const Atom &atom) const;

     const Atom& right(const Atom &atom) const;

     bool hasLeft(const Atom &atom) const;

@@ -80,11 +96,15 @@ public:

     // modify domain

     Atom insert(uint64_t pos, float mass);

     void erase(uint64_t pos);

+    void cacheInsert(uint64_t pos, float mass) const;

+    void cacheErase(uint64_t pos) const;

     void updateMass(uint64_t pos, float newMass);

+    void flushCache();

+    void resetCache(unsigned n);

     // serialization

     friend Archive& operator<<(Archive &ar, AtomicDomain &domain);

     friend Archive& operator>>(Archive &ar, AtomicDomain &domain);

 };

-#endif

+#endif

\ No newline at end of file

@@ -185,6 +185,7 @@ void GibbsSampler<T, MatA, MatB>::update(unsigned nSteps, unsigned nCores)

         mNumQueues += 1.f;

         mAvgQueue = mQueue.size() / mNumQueues + mAvgQueue * (mNumQueues - 1.f) / mNumQueues;

         n += mQueue.size();

+        mDomain.resetCache(mQueue.size());

         //Rprintf("round: %d\n", count++);

         #pragma omp parallel for num_threads(nCores)

@@ -192,6 +193,7 @@ void GibbsSampler<T, MatA, MatB>::update(unsigned nSteps, unsigned nCores)

         {

             processProposal(mQueue[i]);

         }

+        mDomain.flushCache();

         mQueue.clear(1);

         GAPS_ASSERT(n <= nSteps);

     }

@@ -229,23 +231,17 @@ void GibbsSampler<T, MatA, MatB>::processProposal(const AtomicProposal &prop)

 template <class T, class MatA, class MatB>

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

 {

-    #pragma omp critical(gibbs)

-    {

-        mDomain.insert(pos, mass);

-        mMatrix(row, col) += mass;

-        impl()->updateAPMatrix(row, col, mass);

-    }

+    mDomain.cacheInsert(pos, mass);

+    mMatrix(row, col) += mass;

+    impl()->updateAPMatrix(row, col, mass);

 }

 template <class T, class MatA, class MatB>

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

 {

-    #pragma omp critical(gibbs)

-    {

-        mDomain.erase(pos);

-        mMatrix(row, col) += -mass;

-        impl()->updateAPMatrix(row, col, -mass);

-    }

+    mDomain.cacheErase(pos);

+    mMatrix(row, col) += -mass;

+    impl()->updateAPMatrix(row, col, -mass);

 }

 // add an atom at pos, calculate mass either with an exponential distribution

@@ -256,19 +252,16 @@ unsigned col)

 {

     float mass = impl()->canUseGibbs(row, col) ? gibbsMass(row, col, mass)

         : gaps::random::exponential(mLambda);

-    #pragma omp critical(gibbs)

+    if (mass >= gaps::algo::epsilon)

     {

-        if (mass >= gaps::algo::epsilon)

-        {

-            mDomain.updateMass(pos, mass);

-            mMatrix(row, col) += mass;

-            impl()->updateAPMatrix(row, col, mass);

-        }

-        else

-        {

-            mDomain.erase(pos);

-            mQueue.rejectBirth();

-        }

+        mDomain.updateMass(pos, mass);

+        mMatrix(row, col) += mass;

+        impl()->updateAPMatrix(row, col, mass);

+    }

+    else

+    {

+        mDomain.cacheErase(pos);

+        mQueue.rejectBirth();

     }

 }

@@ -279,17 +272,24 @@ void GibbsSampler<T, MatA, MatB>::death(uint64_t pos, float mass, unsigned row,

 unsigned col)

 {

     GAPS_ASSERT(mass > 0.f);

-    removeMass(pos, mass, row, col);

+

+    //removeMass(pos, mass, row, col);

+    mMatrix(row, col) += -mass;

+    impl()->updateAPMatrix(row, col, -mass);

+

     float newMass = impl()->canUseGibbs(row, col) ? gibbsMass(row, col, -mass) : 0.f;

     mass = newMass < gaps::algo::epsilon ? mass : newMass;

     float deltaLL = impl()->computeDeltaLL(row, col, mass);

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

     {

-        addMass(pos, mass, row, col);

+        mDomain.updateMass(pos, mass);

+        mMatrix(row, col) += mass;

+        impl()->updateAPMatrix(row, col, mass);

         mQueue.rejectDeath();

     }

     else

     {

+        mDomain.cacheErase(pos);

         mQueue.acceptDeath();

     }

 }

@@ -374,22 +374,17 @@ template <class T, class MatA, class MatB>

 float GibbsSampler<T, MatA, MatB>::updateAtomMass(uint64_t pos, float mass,

 float delta)

 {

-    bool ret_val = false;

-    #pragma omp critical(gibbs)

+    if (mass + delta < gaps::algo::epsilon)

     {

-        if (mass + delta < gaps::algo::epsilon)

-        {

-            mDomain.erase(pos);

-            mQueue.acceptDeath();

-            ret_val = false;

-        }

-        else

-        {

-            mDomain.updateMass(pos, mass + delta);

-            ret_val = true;

-        }

+        mDomain.cacheErase(pos);

+        mQueue.acceptDeath();

+        return false;

+    }

+    else

+    {

+        mDomain.updateMass(pos, mass + delta);

+        return true;

     }

-    return ret_val;

 }

 // helper function for exchange step, updates the atomic domain, matrix, and

@@ -412,13 +407,10 @@ unsigned r2, unsigned c2)

         mQueue.rejectDeath();

     }

-    #pragma omp critical(gibbs)

-    {

-        mMatrix(r1, c1) += d1;

-        mMatrix(r2, c2) += d2;

-        impl()->updateAPMatrix(r1, c1, d1);

-        impl()->updateAPMatrix(r2, c2, d2);

-    }

+    mMatrix(r1, c1) += d1;

+    mMatrix(r2, c2) += d2;

+    impl()->updateAPMatrix(r1, c1, d1);

+    impl()->updateAPMatrix(r2, c2, d2);

 }

 template <class T, class MatA, class MatB>