@@ -11,6 +11,7 @@

 ^.git

+^src/Makevars((?!\.in))

 ^src/AtomicDomain.o

 ^src/Cogaps.o

 ^src/GapsRunner.o

@@ -24,6 +24,12 @@ static bool vecContains(const std::vector<Atom*> &vec, uint64_t pos)

     return std::binary_search(vec.begin(), vec.end(), &temp, compareAtom);

 }

+// check if a position in contained in a vector of positions

+static bool vecContains(const std::vector<uint64_t> &vec, uint64_t pos)

+{

+    return std::binary_search(vec.begin(), vec.end(), pos);

+}

+

 // used in debug mode to check if vector is always sorted

 static bool isSorted(const std::vector<Atom*> &vec)

 {

@@ -132,11 +138,16 @@ AtomNeighborhood AtomicDomain::randomAtomWithRightNeighbor(GapsRng *rng)

     return AtomNeighborhood(NULL, mAtoms[index], right);

 }

-uint64_t AtomicDomain::randomFreePosition(GapsRng *rng) const

+uint64_t AtomicDomain::randomFreePosition(GapsRng *rng,

+const std::vector<uint64_t> &possibleDeaths) const

 {

     uint64_t pos = rng->uniform64(1, mDomainLength);

     while (vecContains(mAtoms, pos))

     {

+        if (vecContains(possibleDeaths, pos))

+        {

+            return 0; // might actually be a free position

+        }

         pos = rng->uniform64(1, mDomainLength);

     }

     return pos;

@@ -163,6 +174,25 @@ void AtomicDomain::erase(uint64_t pos)

     delete a;

 }

+void AtomicDomain::move(uint64_t src, uint64_t dest)

+{

+    GAPS_ASSERT(size() > 0);

+    GAPS_ASSERT(vecContains(mAtoms, src));

+

+    #pragma omp critical(AtomicInsertOrErase)

+    {

+        std::vector<Atom*>::iterator it = std::lower_bound(mAtoms.begin(), mAtoms.end(), src, compareAtomLower);

+        unsigned ndx = std::distance(mAtoms.begin(), it);

+        while (ndx + 1 < mAtoms.size() && dest > mAtoms[ndx + 1]->pos)

+        {

+            Atom* temp = mAtoms[ndx];

+            mAtoms[ndx] = mAtoms[ndx + 1];

+            mAtoms[ndx + 1] = temp;

+            ++ndx;

+        }

+    }

+}

+

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

 {

     Atom *newAtom = new Atom(pos, mass);

@@ -33,6 +33,8 @@ struct AtomNeighborhood

     bool hasRight();

 };

+class ProposalQueue; // needed for friend

+

 class AtomicDomain

 {

 public:

@@ -49,12 +51,13 @@ public:

     Atom* getLeftNeighbor(uint64_t pos);

     Atom* getRightNeighbor(uint64_t pos);

-    uint64_t randomFreePosition(GapsRng *rng) const;

+    uint64_t randomFreePosition(GapsRng *rng,

+        const std::vector<uint64_t> &possibleDeaths) const;

     uint64_t size() const;

     // these need to happen concurrently without invalidating pointers

     void erase(uint64_t pos);

-    Atom* insert(uint64_t pos, float mass);

+    void move(uint64_t src, uint64_t dest);

     // serialization

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

@@ -64,6 +67,10 @@ public:

 private:

 #endif

+    // only the proposal queue can insert

+    friend class ProposalQueue;

+    Atom* insert(uint64_t pos, float mass);

+

     // size of atomic domain to ensure all bins are equal length

     uint64_t mDomainLength;

@@ -62,6 +62,7 @@ unsigned getNumPatterns(const Rcpp::List &allParams)

     {

         std::string file(Rcpp::as<std::string>(allParams["checkpointInFile"]));

         Archive ar(file, ARCHIVE_READ);

+        GapsRng::load(ar);

         ar >> nPatterns;

         ar.close();

     }

@@ -99,6 +100,7 @@ const Rcpp::Nullable<Rcpp::NumericMatrix> &fixedMatrix, bool isMaster)

 {

     // calculate essential parameters needed for constructing GapsRunner

     const Rcpp::S4 &gapsParams(allParams["gaps"]);

+    GapsRng::setSeed(gapsParams.slot("seed"));

     unsigned nPatterns = getNumPatterns(allParams); // TODO clarify this sets the checkpoint seed as well

     bool printThreads = !processDistributedParameters(allParams).first;

     bool partitionRows = processDistributedParameters(allParams).second;

@@ -106,7 +108,7 @@ const Rcpp::Nullable<Rcpp::NumericMatrix> &fixedMatrix, bool isMaster)

     // construct GapsRunner

     GapsRunner runner(data, allParams["transposeData"], nPatterns,

-        partitionRows, cIndices, gapsParams.slot("seed"));

+        partitionRows, cIndices);

     // set uncertainty

     if (!isNull(uncertainty))

@@ -120,6 +122,7 @@ const Rcpp::Nullable<Rcpp::NumericMatrix> &fixedMatrix, bool isMaster)

     {

         std::string file(Rcpp::as<std::string>(allParams["checkpointInFile"]));

         Archive ar(file, ARCHIVE_READ);

+        GapsRng::load(ar);

         ar >> runner;

         ar.close();

     }

@@ -134,6 +137,7 @@ const Rcpp::Nullable<Rcpp::NumericMatrix> &fixedMatrix, bool isMaster)

         }

         // set parameters that would be saved in the checkpoint

+        runner.recordSeed(gapsParams.slot("seed"));

         runner.setMaxIterations(gapsParams.slot("nIterations"));

         runner.setSparsity(gapsParams.slot("alphaA"),

             gapsParams.slot("alphaP"), gapsParams.slot("singleCell"));

@@ -27,6 +27,11 @@ void GapsRunner::setFixedMatrix(char which, const Matrix &mat)

     }

 }

+void GapsRunner::recordSeed(uint32_t seed)

+{

+    mSeed = seed;

+}

+

 uint32_t GapsRunner::getSeed() const

 {

     return mSeed;

@@ -52,8 +52,7 @@ public:

     template <class DataType>

     GapsRunner(const DataType &data, bool transposeData, unsigned nPatterns,

-        bool partitionRows, const std::vector<unsigned> &indices,

-        uint32_t seed);

+        bool partitionRows, const std::vector<unsigned> &indices);

     template <class DataType>

     void setUncertainty(const DataType &unc, bool transposeData,

@@ -61,6 +60,7 @@ public:

     void setFixedMatrix(char which, const Matrix &mat);

+    void recordSeed(uint32_t seed);

     uint32_t getSeed() const;

     void setMaxIterations(unsigned nIterations);

@@ -82,8 +82,7 @@ public:

 // problem with passing file parser - need to read it twice

 template <class DataType>

 GapsRunner::GapsRunner(const DataType &data, bool transposeData,

-unsigned nPatterns, bool partitionRows, const std::vector<unsigned> &indices,

-uint32_t seed)

+unsigned nPatterns, bool partitionRows, const std::vector<unsigned> &indices)

     :

 mASampler(data, !transposeData, nPatterns,!partitionRows, indices),

 mPSampler(data, transposeData, nPatterns, partitionRows, indices),

@@ -91,14 +90,10 @@ mStatistics(mPSampler.dataRows(), mPSampler.dataCols(), nPatterns),

 mFixedMatrix('N'), mMaxIterations(1000), mMaxThreads(1), mPrintMessages(true),

 mOutputFrequency(500), mCheckpointOutFile("gaps_checkpoint.out"),

 mCheckpointInterval(0), mPhase('C'), mCurrentIteration(0),

-mNumPatterns(nPatterns), mSeed(seed), mNumUpdatesA(0), mNumUpdatesP(0),

-mRng(seed)

+mNumPatterns(nPatterns), mSeed(0), mNumUpdatesA(0), mNumUpdatesP(0)

 {

     mASampler.sync(mPSampler);

     mPSampler.sync(mASampler);

-

-    mASampler.setSeed(mRng.uniform64());

-    mPSampler.setSeed(mRng.uniform64());

 }

 template <class DataType>

@@ -44,11 +44,6 @@ void GibbsSampler::setMatrix(const Matrix &mat)

     mMatrix = mat;

 }

-void GibbsSampler::setSeed(uint64_t seed)

-{

-    mSeeder.seed(seed);

-}

-

 float GibbsSampler::chi2() const

 {   

     return 2.f * gaps::algo::loglikelihood(mDMatrix, mSMatrix, mAPMatrix);

@@ -79,26 +74,19 @@ void GibbsSampler::update(unsigned nSteps, unsigned nCores)

         mQueue.populate(mDomain, nSteps - n);

         n += mQueue.size();

-        // update average queue count

-        #ifdef GAPS_DEBUG

-        mNumQueues += 1.f;

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

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

-        #endif

-

         // process all proposed updates

         #pragma omp parallel for num_threads(nCores)

         for (unsigned i = 0; i < mQueue.size(); ++i)

         {

-            processProposal(&mQueue[i]);

+            processProposal(mQueue[i]);

         }

         mQueue.clear();

     }

 }

-void GibbsSampler::processProposal(AtomicProposal *prop)

+void GibbsSampler::processProposal(const AtomicProposal &prop)

 {

-    switch (prop->type)

+    switch (prop.type)

     {

         case 'B':

             birth(prop);

@@ -117,45 +105,40 @@ void GibbsSampler::processProposal(AtomicProposal *prop)

 // add an atom at a random position, calculate mass either with an

 // exponential distribution or with the gibbs mass distribution

-void GibbsSampler::birth(AtomicProposal *prop)

+void GibbsSampler::birth(const AtomicProposal &prop)

 {

-    unsigned row = getRow(prop->atom1->pos);

-    unsigned col = getCol(prop->atom1->pos);

-

     // calculate proposed mass

-    float mass = canUseGibbs(col)

-        ? gibbsMass(alphaParameters(row, col), &(prop->rng)).value()

-        : prop->rng.exponential(mLambda);

+    float mass = canUseGibbs(prop.c1)

+        ? gibbsMass(alphaParameters(prop.r1, prop.c1), &(prop.rng)).value()

+        : prop.rng.exponential(mLambda);

     // accept mass as long as it's non-zero

     if (mass >= gaps::epsilon)

     {

         mQueue.acceptBirth();

-        prop->atom1->mass = mass;

-        changeMatrix(row, col, mass);

+        prop.atom1->mass = mass;

+        changeMatrix(prop.r1, prop.c1, mass);

     }

     else

     {

         mQueue.rejectBirth();

-        mDomain.erase(prop->atom1->pos);

+        mDomain.erase(prop.atom1->pos);

     }

 }

 // automatically accept death, attempt a rebirth at the same position, using

 // the original mass or the gibbs mass distribution

-void GibbsSampler::death(AtomicProposal *prop)

+void GibbsSampler::death(const AtomicProposal &prop)

 {

-    unsigned row = getRow(prop->atom1->pos);

-    unsigned col = getCol(prop->atom1->pos);

-

     // calculate alpha parameters assuming atom dies

-    AlphaParameters alpha = alphaParametersWithChange(row, col, -prop->atom1->mass);

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

+        -prop.atom1->mass);

     // try to calculate rebirth mass using gibbs distribution, otherwise exponetial

-    float rebirthMass = prop->atom1->mass;

-    if (canUseGibbs(col))

+    float rebirthMass = prop.atom1->mass;

+    if (canUseGibbs(prop.c1))

     {

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

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

         if (gMass.hasValue())

         {

             rebirthMass = gMass.value();

@@ -164,84 +147,73 @@ void GibbsSampler::death(AtomicProposal *prop)

     // accept/reject rebirth

     float deltaLL = getDeltaLL(alpha, rebirthMass) * mAnnealingTemp;

-    if (std::log(prop->rng.uniform()) < deltaLL)

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

     {

         mQueue.rejectDeath();

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

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

         {

-            safelyChangeMatrix(row, col, rebirthMass - prop->atom1->mass);

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

         }

-        prop->atom1->mass = rebirthMass;

+        prop.atom1->mass = rebirthMass;

     }

     else

     {

         mQueue.acceptDeath();

-        safelyChangeMatrix(row, col, -prop->atom1->mass);

-        mDomain.erase(prop->atom1->pos);

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

+        mDomain.erase(prop.atom1->pos);

     }

 }

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

-void GibbsSampler::move(AtomicProposal *prop)

+void GibbsSampler::move(const AtomicProposal &prop)

 {

-    unsigned r1 = getRow(prop->atom1->pos);

-    unsigned c1 = getCol(prop->atom1->pos);

-    unsigned r2 = getRow(prop->pos);

-    unsigned c2 = getCol(prop->pos);

-    GAPS_ASSERT(r1 != r2 || c1 != c2);

-

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

-    if (std::log(prop->rng.uniform()) < getDeltaLL(alpha, -prop->atom1->mass) * mAnnealingTemp)

+    AlphaParameters alpha = alphaParameters(prop.r1, prop.c1, prop.r2, prop.c2);

+    if (std::log(prop.rng.uniform()) < getDeltaLL(alpha, -prop.atom1->mass) * mAnnealingTemp)

     {

-        prop->atom1->pos = prop->pos;

-        safelyChangeMatrix(r1, c1, -prop->atom1->mass);

-        changeMatrix(r2, c2, prop->atom1->mass);

+        //prop.atom1->pos = prop.pos;

+        mDomain.move(prop.atom1->pos, prop.pos);

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

+        changeMatrix(prop.r2, prop.c2, prop.atom1->mass);

     }

 }

 // exchange some amount of mass between two positions, note it is possible

 // for one of the atoms to be deleted if it's mass becomes too small

-void GibbsSampler::exchange(AtomicProposal *prop)

+void GibbsSampler::exchange(const AtomicProposal &prop)

 {

-    unsigned r1 = getRow(a1->pos);

-    unsigned c1 = getCol(a1->pos);

-    unsigned r2 = getRow(a2->pos);

-    unsigned c2 = getCol(a2->pos);

-    GAPS_ASSERT(r1 != r2 || c1 != c2);

-

     // attempt gibbs distribution exchange

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

-    if (canUseGibbs(c1, c2))

+    AlphaParameters alpha = alphaParameters(prop.r1, prop.c1, prop.r2, prop.c2);

+    if (canUseGibbs(prop.c1, prop.c2))

     {

-        OptionalFloat gMass = gibbsMass(alpha, prop->atom1->mass,

-            prop->atom2->mass, &(prop->rng));

+        OptionalFloat gMass = gibbsMass(alpha, prop.atom1->mass,

+            prop.atom2->mass, &(prop.rng));

         if (gMass.hasValue())

         {

-            acceptExchange(prop->atom1, prop->atom2, gMass.value(), r1, c1, r2, c2);

+            acceptExchange(prop, gMass.value());

             return;

         }

     }

     // resort to metropolis-hastings if gibbs fails

-    exchangeUsingMetropolisHastings(prop, alpha, r1, c1, r2, c2);

+    exchangeUsingMetropolisHastings(prop, alpha);

 }

-void GibbsSampler::exchangeUsingMetropolisHastings(AtomicProposal *prop,

-AlphaParameters alpha, unsigned r1, unsigned c1, unsigned r2, unsigned c2)

+void GibbsSampler::exchangeUsingMetropolisHastings(const AtomicProposal &prop,

+AlphaParameters alpha)

 {

     // compute amount of mass to be exchanged

-    float totalMass = prop->atom1->mass + prop->atom2->mass;

-    float newMass = prop->rng.truncGammaUpper(totalMass, 2.f, 1.f / mLambda);

+    float totalMass = prop.atom1->mass + prop.atom2->mass;

+    float newMass = prop.rng.truncGammaUpper(totalMass, 2.f, 1.f / mLambda);

     // compute amount to change atom1 by - always change larger mass to newMass

-    float delta = (prop->atom1->mass > prop->atom2->mass)

-        ? newMass - prop->atom1->mass

-        : prop->atom2->mass - newMass;

+    float delta = (prop.atom1->mass > prop.atom2->mass)

+        ? newMass - prop.atom1->mass

+        : prop.atom2->mass - newMass;

     // choose mass for priorLL calculation

     float oldMass = (2.f * newMass > totalMass)

-        ? gaps::max(prop->atom1->mass, prop->atom2->mass)

-        : gaps::min(prop->atom1->mass, prop->atom2->mass);

+        ? gaps::max(prop.atom1->mass, prop.atom2->mass)

+        : gaps::min(prop.atom1->mass, prop.atom2->mass);

     // calculate priorLL

     float pNew = gaps::d_gamma(newMass, 2.f, 1.f / mLambda);

@@ -250,24 +222,23 @@ AlphaParameters alpha, unsigned r1, unsigned c1, unsigned r2, unsigned c2)

     // accept/reject

     float deltaLL = getDeltaLL(alpha, delta) * mAnnealingTemp;

-    if (priorLL == 0.f || std::log(prop->rng.uniform() * priorLL) < deltaLL)

+    if (priorLL == 0.f || std::log(prop.rng.uniform() * priorLL) < deltaLL)

     {

-        acceptExchange(prop->atom1, prop->atom2, delta, r1, c1, r2, c2);

+        acceptExchange(prop, delta);

         return;

     }

 }

 // helper function for exchange step

-void GibbsSampler::acceptExchange(Atom *a1, Atom *a2, float delta,

-unsigned r1, unsigned c1, unsigned r2, unsigned c2)

+void GibbsSampler::acceptExchange(const AtomicProposal &prop, float delta)

 {

-    if (a1->mass + delta > gaps::epsilon && a2->mass - delta > gaps::epsilon)

+    if (prop.atom1->mass + delta > gaps::epsilon && prop.atom2->mass - delta > gaps::epsilon)

     {

-        a1->mass += delta;

-        a2->mass -= delta;

+        prop.atom1->mass += delta;

+        prop.atom2->mass -= delta;

-        changeMatrix(r1, c1, delta);

-        changeMatrix(r2, c2, -delta);

+        changeMatrix(prop.r1, prop.c1, delta);

+        changeMatrix(prop.r2, prop.c2, -delta);

     }

 }

@@ -319,7 +290,15 @@ float m1, float m2, GapsRng *rng)

     return OptionalFloat();

 }

-// needed to prevent negative values in matrix

+// here mass + delta is guaranteed to be positive

+void GibbsSampler::changeMatrix(unsigned row, unsigned col, float delta)

+{

+    mMatrix(row, col) += delta;

+    GAPS_ASSERT(mMatrix(row, col) >= 0.f);

+    updateAPMatrix(row, col, delta);

+}

+

+// delta could be negative, this is needed to prevent negative values in matrix

 void GibbsSampler::safelyChangeMatrix(unsigned row, unsigned col, float delta)

 {

     float newVal = gaps::max(mMatrix(row, col) + delta, 0.f);

@@ -350,16 +329,6 @@ void GibbsSampler::updateAPMatrix(unsigned row, unsigned col, float delta)

     }

 }

-unsigned GibbsSampler::getRow(uint64_t pos) const

-{

-    return pos / (mBinSize * mNumPatterns); // nCol == nPatterns

-}

-

-unsigned GibbsSampler::getCol(uint64_t pos) const

-{

-    return (pos / mBinSize) % mNumPatterns; // nCol == nPatterns

-}

-

 bool GibbsSampler::canUseGibbs(unsigned col) const

 {

     return !gaps::algo::isVectorZero(mOtherMatrix->colPtr(col),

@@ -2,6 +2,7 @@

 #define __COGAPS_GIBBS_SAMPLER_H__

 #include "AtomicDomain.h"

+#include "ProposalQueue.h"

 #include "data_structures/Matrix.h"

 #include "math/Algorithms.h"

 #include "math/Random.h"

@@ -27,7 +28,6 @@ public:

     void setMaxGibbsMass(float max);

     void setAnnealingTemp(float temp);

     void setMatrix(const Matrix &mat);

-    void setSeed(uint64_t seed);

     float chi2() const;

     uint64_t nAtoms() const;

@@ -68,21 +68,22 @@ private:

     uint64_t mBinSize;

     uint64_t mDomainLength;

-    void makeAndProcessProposal(GapsRng *rng);

+    void processProposal(const AtomicProposal &prop);

     float deathProb(uint64_t nAtoms) const;

-    void birth(GapsRng *rng);

-    void death(GapsRng *rng);

-    void move(GapsRng *rng);

-    void exchange(GapsRng *rng);

-

-    void acceptExchange(Atom *a1, Atom *a2, float d1, unsigned r1,

-        unsigned c1, unsigned r2, unsigned c2);

+    void birth(const AtomicProposal &prop);

+    void death(const AtomicProposal &prop);

+    void move(const AtomicProposal &prop);

+    void exchange(const AtomicProposal &prop);

+    void exchangeUsingMetropolisHastings(const AtomicProposal &prop,

+        AlphaParameters alpha);

+    void acceptExchange(const AtomicProposal &prop, float delta);

     bool updateAtomMass(Atom *atom, float delta);

     OptionalFloat gibbsMass(AlphaParameters alpha, GapsRng *rng);

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

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

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

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

@@ -105,6 +106,7 @@ mAPMatrix(mDMatrix.nRow(), mDMatrix.nCol()),

 mMatrix(mDMatrix.nCol(), nPatterns),

 mOtherMatrix(NULL),

 mDomain(mMatrix.nRow() * mMatrix.nCol()),

+mQueue(mMatrix.nRow(), mMatrix.nCol()),

 mLambda(0.f),

 mMaxGibbsMass(100.f),

 mAnnealingTemp(1.f),

@@ -8,8 +8,10 @@ OBJECTS =   AtomicDomain.o \

             GapsRunner.o \

             GapsStatistics.o \

             GibbsSampler.o \

+            ProposalQueue.o \

             RcppExports.o \

             test-runner.o \

+            data_structures/HashSets.o \

             data_structures/Matrix.o \

             data_structures/Vector.o \

             file_parser/CsvParser.o \

@@ -21,7 +23,6 @@ OBJECTS =   AtomicDomain.o \

             math/Random.o \

             cpp_tests/testAlgorithms.o \

             cpp_tests/testAtomicDomain.o \

-            cpp_tests/testEfficientSets.o \

             cpp_tests/testFileParsers.o \

             cpp_tests/testMatrix.o \

             cpp_tests/testRandom.o \

@@ -1,6 +1,6 @@

-#include "GapsAssert.h"

 #include "ProposalQueue.h"

 #include "math/Random.h"

+#include "utils/GapsAssert.h"

 //////////////////////////////// AtomicProposal ////////////////////////////////

@@ -11,20 +11,24 @@ AtomicProposal::AtomicProposal(char t)

 //////////////////////////////// ProposalQueue /////////////////////////////////

-ProposalQueue::ProposalQueue(unsigned primaryDimSize, unsigned secondaryDimSize)

+ProposalQueue::ProposalQueue(unsigned nrow, unsigned ncol)

     :

-mMinAtoms(0), mMaxAtoms(0), mNumBins(primaryDimSize * secondaryDimSize),

-mBinLength(std::numeric_limits<uint64_t>::max() / mNumBins),

-mSecondaryDimLength(mBinLength * secondaryDimSize),

-mDomainLength(mBinLength * mNumBins), mSecondaryDimSize(secondaryDimSize),

-mAlpha(0.f), mU1(0.f), mU2(0.f), mUseCachedRng(false)

-{

-    mUsedIndices.setDimensionSize(primaryDimSize);

-}

+mUsedMatrixIndices(nrow),

+mMinAtoms(0),

+mMaxAtoms(0),

+mBinLength(std::numeric_limits<uint64_t>::max() / static_cast<uint64_t>(nrow * ncol)),

+mNumCols(ncol),

+mAlpha(0.0),

+mDomainLength(static_cast<double>(mBinLength * static_cast<uint64_t>(nrow * ncol))),

+mNumBins(static_cast<double>(nrow * ncol)),

+mU1(0.f),

+mU2(0.f),

+mUseCachedRng(false)

+{}

 void ProposalQueue::setAlpha(float alpha)

 {

-    mAlpha = alpha;

+    mAlpha = static_cast<double>(alpha);

 }

 void ProposalQueue::populate(AtomicDomain &domain, unsigned limit)

@@ -36,9 +40,9 @@ void ProposalQueue::populate(AtomicDomain &domain, unsigned limit)

     bool success = true;

     while (nIter++ < limit && success)

     {

-        success = makeProposal(domain);

-        if (!success)

+        if (!makeProposal(domain))

         {

+            success = false;

             mUseCachedRng = true;

         }

     }

@@ -49,8 +53,8 @@ void ProposalQueue::clear()

     GAPS_ASSERT(mMinAtoms == mMaxAtoms);

     mQueue.clear();

-    mUsedPositions.clear();

-    mUsedIndices.clear();

+    mUsedMatrixIndices.clear();

+    mUsedAtoms.clear();

 }

 unsigned ProposalQueue::size() const

@@ -69,25 +73,25 @@ AtomicProposal& ProposalQueue::operator[](int n)

 void ProposalQueue::acceptDeath()

 {

     #pragma omp atomic

-    mMaxAtoms--;

+    --mMaxAtoms;

 }

 void ProposalQueue::rejectDeath()

 {

     #pragma omp atomic

-    mMinAtoms++;

+    ++mMinAtoms;

 }

 void ProposalQueue::acceptBirth()

 {

     #pragma omp atomic

-    mMinAtoms++;

+    ++mMinAtoms;

 }

 void ProposalQueue::rejectBirth()

 {

     #pragma omp atomic

-    mMaxAtoms--;

+    --mMaxAtoms;

 }

 float ProposalQueue::deathProb(double nAtoms) const

@@ -98,16 +102,14 @@ float ProposalQueue::deathProb(double nAtoms) const

 bool ProposalQueue::makeProposal(AtomicDomain &domain)

 {

-    // special indeterminate case

     if (mMinAtoms < 2 && mMaxAtoms >= 2)

     {

-        return false;

+        return false; // special indeterminate case

     }

-    // always birth when no atoms exist

     if (mMaxAtoms < 2)

     {

-        return birth(domain);

+        return birth(domain); // always birth when no atoms exist

     }

     mU1 = mUseCachedRng ? mU1 : mRng.uniform();

@@ -131,24 +133,29 @@ bool ProposalQueue::makeProposal(AtomicDomain &domain)

     }

     return (mU1 < 0.75f) ? move(domain) : exchange(domain);

 }

-    

+

 bool ProposalQueue::birth(AtomicDomain &domain)

 {

     AtomicProposal prop('B');

-    prop.r1 = (prop.atom1->pos / mBinSize) / mNumCols;

-    prop.c1 = (prop.atom1->pos / mBinSize) % mNumCols;

-

-    prop.init(); // initialize rng

+    uint64_t pos = domain.randomFreePosition(&(prop.rng), mUsedAtoms.vec());

+    if (pos == 0)

+    {

+        DEBUG_PING // want to notify since this event should have near 0 probability

+        GapsRng::rollBackOnce(); // ensure same proposal next time

+        return false; // atom conflict, might have open spot if atom moves/dies

+    }

+    prop.r1 = (pos / mBinLength) / mNumCols;

+    prop.c1 = (pos / mBinLength) % mNumCols;

     if (mUsedMatrixIndices.contains(prop.r1))

     {

-        return false; // matrix conflict - can't compute alpha parameters

+        GapsRng::rollBackOnce(); // ensure same proposal next time

+        return false; // matrix conflict - can't compute gibbs mass

     }

-    prop.atom1 = mDomain.insert(domain.randomFreePosition(), 0.f);

+    prop.atom1 = domain.insert(pos, 0.f);

     mUsedMatrixIndices.insert(prop.r1);

     mUsedAtoms.insert(prop.atom1->pos);

-

     mQueue.push_back(prop);

     ++mMaxAtoms;

     return true;

@@ -157,18 +164,18 @@ bool ProposalQueue::birth(AtomicDomain &domain)

 bool ProposalQueue::death(AtomicDomain &domain)

 {

     AtomicProposal prop('D');

-    prop->atom1 = domain.randomAtom();

-    prop.r1 = (prop.atom1->pos / mBinSize) / mNumCols;

-    prop.c1 = (prop.atom1->pos / mBinSize) % mNumCols;

+    prop.atom1 = domain.randomAtom(&(prop.rng));

+    prop.r1 = (prop.atom1->pos / mBinLength) / mNumCols;

+    prop.c1 = (prop.atom1->pos / mBinLength) % mNumCols;

     if (mUsedMatrixIndices.contains(prop.r1))

     {

-        return false; // matrix conflict - can't compute alpha parameters

+        GapsRng::rollBackOnce(); // ensure same proposal next time

+        return false; // matrix conflict - can't compute gibbs mass or deltaLL

     }

+

     mUsedMatrixIndices.insert(prop.r1);

     mUsedAtoms.insert(prop.atom1->pos);

-

-    prop.init(); // initialize rng

     mQueue.push_back(prop);

     --mMinAtoms;

     return true;

@@ -177,107 +184,89 @@ bool ProposalQueue::death(AtomicDomain &domain)

 bool ProposalQueue::move(AtomicDomain &domain)

 {

     AtomicProposal prop('M');

+    AtomNeighborhood hood = domain.randomAtomWithNeighbors(&(prop.rng));

-    AtomNeighborhood hood = domain.randomAtomWithNeighbors();

-    

     uint64_t lbound = hood.hasLeft() ? hood.left->pos : 0;

     uint64_t rbound = hood.hasRight() ? hood.right->pos : mDomainLength;

-

-    uint64_t newLocation = mRng.uniform64(lbound + 1, rbound - 1);

-

     if (mUsedAtoms.contains(lbound) || mUsedAtoms.contains(rbound))

     {

+        GapsRng::rollBackOnce();

         return false; // atomic conflict - don't know neighbors

     }

-    if (primaryIndex(hood.center->pos) == primaryIndex(newLocation)

-    && secondaryIndex(hood.center->pos) == secondaryIndex(newLocation))

+    prop.pos = prop.rng.uniform64(lbound + 1, rbound - 1);

+    prop.atom1 = hood.center;

+    prop.r1 = (prop.atom1->pos / mBinLength) / mNumCols;

+    prop.c1 = (prop.atom1->pos / mBinLength) % mNumCols;

+    prop.r2 = (prop.pos / mBinLength) / mNumCols;

+    prop.c2 = (prop.pos / mBinLength) % mNumCols;

+

+    if (mUsedMatrixIndices.contains(prop.r1) || mUsedMatrixIndices.contains(prop.r2))

     {

-        hood.center->pos = newLocation; // automatically accept moves in same bin

-        return true;

+        GapsRng::rollBackOnce(); // ensure same proposal next time

+        return false; // matrix conflict - can't compute deltaLL

     }

-    if (mUsedIndices.contains(primaryIndex(hood.center->pos))

-    || mUsedIndices.contains(primaryIndex(newLocation)))

+    if (prop.r1 == prop.r2 && prop.c1 == prop.c2)

     {

-        return false; // matrix conflict - can't compute deltaLL

+        prop.atom1->pos = prop.pos; // automatically accept moves in same bin

+        return true;

     }

-    mQueue.push_back(AtomicProposal('M', hood.center, newLocation));

-    mUsedIndices.insert(primaryIndex(hood.center->pos));

-    mUsedIndices.insert(primaryIndex(newLocation));

-    mUsedPositions.insert(hood.center->pos);

-    mUsedPositions.insert(newLocation);

+    mQueue.push_back(prop);

+    mUsedMatrixIndices.insert(prop.r1);

+    mUsedMatrixIndices.insert(prop.r2);

+    mUsedAtoms.insert(prop.atom1->pos);

+    mUsedAtoms.insert(prop.pos);

     return true;

 }

 bool ProposalQueue::exchange(AtomicDomain &domain)

 {

-    AtomNeighborhood hood = domain.randomAtomWithRightNeighbor();

-    Atom* a1 = hood.center;

-    Atom* a2 = hood.hasRight() ? hood.right : domain.front();

+    AtomicProposal prop('E');

+    AtomNeighborhood hood = domain.randomAtomWithRightNeighbor(&(prop.rng));

+    prop.atom1 = hood.center;

+    prop.atom2 = hood.hasRight() ? hood.right : domain.front();

-    if (hood.hasRight()) // has neighbor

-    {

-        if (!mUsedPositions.isEmptyInterval(a1->pos, a2->pos))

-        {

-            return false; // atomic conflict - don't know right neighbor

-        }

-    }

-    else // exchange with first atom

-    {

-        if (!mUsedPositions.isEmptyInterval(a1->pos, mDomainLength))

-        {

-            return false; // atomic conflict - don't know right neighbor

-        }

-        

-        if (!mUsedPositions.isEmptyInterval(0, domain.front()->pos))

-        {

-            return false; // atomic conflict - don't know right neighbor

-        }

-    }

+    prop.r1 = (prop.atom1->pos / mBinLength) / mNumCols;

+    prop.c1 = (prop.atom1->pos / mBinLength) % mNumCols;

+    prop.r2 = (prop.atom2->pos / mBinLength) / mNumCols;

+    prop.c2 = (prop.atom2->pos / mBinLength) % mNumCols;

-    if (primaryIndex(a1->pos) == primaryIndex(a2->pos)

-    && secondaryIndex(a1->pos) == secondaryIndex(a2->pos))

+    if (mUsedMatrixIndices.contains(prop.r1) || mUsedMatrixIndices.contains(prop.r2))

     {

-        GapsRng rng;

-        float newMass = rng.truncGammaUpper(a1->mass + a2->mass, 2.f, 1.f / mLambda);

-        float delta = (a1->mass > a2->mass) ? newMass - a1->mass : a2->mass - newMass;

-        if (a1->mass + delta > gaps::epsilon && a2->mass - delta > gaps::epsilon)

-        {

-            a1->mass += delta;

-            a2->mass -= delta;

-        }        

-        return true;

+        GapsRng::rollBackOnce(); // ensure same proposal next time

+        return false; // matrix conflict - can't compute deltaLL or gibbs mass

     }

-    if (mUsedIndices.contains(primaryIndex(a1->pos))

-    || mUsedIndices.contains(primaryIndex(a2->pos)))

+    if (prop.r1 == prop.r2 && prop.c1 == prop.c2)

     {

-        return false; // matrix conflict - can't compute gibbs mass or deltaLL

+        //float newMass = prop.rng.truncGammaUpper(a1->mass + a2->mass, 2.f, 1.f / mLambda);

+        //float delta = (a1->mass > a2->mass) ? newMass - a1->mass : a2->mass - newMass;

+        //if (a1->mass + delta > gaps::epsilon && a2->mass - delta > gaps::epsilon)

+        //{

+        //    a1->mass += delta;

+        //    a2->mass -= delta;

+        //}        

+        return true; // TODO automatically accept exchanges in same bin

     }

-    mQueue.push_back(AtomicProposal('E', a1, a2));

-    mUsedIndices.insert(primaryIndex(a1->pos));

-    mUsedIndices.insert(primaryIndex(a2->pos));

-    mUsedPositions.insert(a1->pos);

-    mUsedPositions.insert(a2->pos);

-    --mMinAtoms;

+    mQueue.push_back(prop);

+    mUsedMatrixIndices.insert(prop.r1);

+    mUsedMatrixIndices.insert(prop.r2);

     return true;

 }

 Archive& operator<<(Archive &ar, ProposalQueue &q)

 {

-    ar << q.mMinAtoms << q.mMaxAtoms << q.mNumBins << q.mBinLength

-        << q.mSecondaryDimLength << q.mDomainLength << q.mSecondaryDimSize

-        << q.mAlpha << q.mRng;

+    ar << q.mRng << q.mMinAtoms << q.mMaxAtoms << q.mBinLength << q.mNumCols

+        << q.mAlpha << q.mDomainLength << q.mNumBins;

     return ar;

 }

 Archive& operator>>(Archive &ar, ProposalQueue &q)

 {

-    ar >> q.mMinAtoms >> q.mMaxAtoms >> q.mNumBins >> q.mBinLength