@@ -79,6 +79,9 @@ private:

     Atom& _left(const Atom &atom);

     Atom& _right(const Atom &atom);

+    void erase(uint64_t pos);

+    Atom insert(uint64_t pos, float mass);

+

 public:

     void setDomainSize(uint64_t size) { mDomainSize = size; }

@@ -95,8 +98,6 @@ public:

     bool hasRight(const Atom &atom) const;

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

@@ -41,6 +41,7 @@ GapsResult GapsDispatcher::run()

     if (mPrintMessages)

     {

         unsigned availableCores = omp_get_max_threads();

+        mNumCoresPerSet = std::min(availableCores, mNumCoresPerSet);

         gaps_printf("Running on %d out of %d available cores\n",

             mNumCoresPerSet, availableCores);

     }

@@ -300,14 +300,11 @@ unsigned col)

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

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

     {

-        mDomain.updateMass(pos, mass);

-        mMatrix(row, col) += mass;

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

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

         mQueue.acceptBirth();

     }

     else

     {

-        mDomain.cacheErase(pos);

         mQueue.rejectBirth();

     }

 }

@@ -136,7 +136,6 @@ bool ProposalQueue::birth(AtomicDomain &domain)

     mQueue.push_back(AtomicProposal('B', pos));

     mUsedIndices.insert(pos / mDimensionSize);

     mUsedPositions.insert(pos);

-    domain.insert(pos, 0.f);

     ++mMaxAtoms;

     return true;

 }