@@ -111,7 +111,7 @@ Atom* AtomicDomain::front()

     return mAtoms.front();

 }

-Atom* AtomicDomain::randomAtom(GapsRng *rng)

+Atom* AtomicDomain::randomAtom(GapsRng *rng, const SmallPairedHashSetU64 &moves)

 {

     GAPS_ASSERT(size() > 0);

     GAPS_ASSERT(isSorted());

@@ -126,6 +126,12 @@ Atom* AtomicDomain::randomAtom(GapsRng *rng)

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

     unsigned leftIndex = index == 0 ? 0 : index - 1;

     index = (index == mAtoms.size()) ? index - 1 : index;

+

+    if (moves.contains(mAtoms[leftIndex]->pos) || moves.contains(mAtoms[index]->pos))

+    {

+        return NULL;

+    }

+

     if (std::abs(mAtoms[leftIndex]->pos - pos) < std::abs(mAtoms[index]->pos - pos))

     {

         index = leftIndex;

@@ -151,12 +157,14 @@ AtomNeighborhood AtomicDomain::randomAtomWithNeighbors(GapsRng *rng)

         index = leftIndex;

     }

+    //gaps_printf("index: %d, size: %d\n", index, mAtoms.size());

     Atom* left = (index == 0) ? NULL : mAtoms[index - 1];

     Atom* right = (index == mAtoms.size() - 1) ? NULL : mAtoms[index + 1];

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

 }

-AtomNeighborhood AtomicDomain::randomAtomWithRightNeighbor(GapsRng *rng)

+AtomNeighborhood AtomicDomain::randomAtomWithRightNeighbor(GapsRng *rng,

+const SmallPairedHashSetU64 &moves)

 {

     GAPS_ASSERT(size() > 0);

@@ -167,6 +175,12 @@ AtomNeighborhood AtomicDomain::randomAtomWithRightNeighbor(GapsRng *rng)

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

     unsigned leftIndex = index == 0 ? 0 : index - 1;

     index = (index == mAtoms.size()) ? index - 1 : index;

+

+    if (moves.contains(mAtoms[leftIndex]->pos) || moves.contains(mAtoms[index]->pos))

+    {

+        return AtomNeighborhood(NULL, NULL, NULL);

+    }

+

     if (std::abs(mAtoms[leftIndex]->pos - pos) < std::abs(mAtoms[index]->pos - pos))

     {

         index = leftIndex;

@@ -200,6 +214,22 @@ uint64_t AtomicDomain::size() const

     return mAtoms.size();

 }

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

+{

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

+    std::vector<Atom*>::iterator it;

+    #pragma omp critical(AtomicInsertOrErase)

+    {

+        GAPS_ASSERT(!vecContains(mAtoms, pos));

+

+        it = std::lower_bound(mAtoms.begin(), mAtoms.end(), pos, compareAtomLower);

+        it = mAtoms.insert(it, newAtom);

+

+        GAPS_ASSERT(vecContains(mAtoms, pos));

+    }

+    return *it;

+}

+

 void AtomicDomain::erase(uint64_t pos)

 {

     Atom *a = NULL;

@@ -210,7 +240,7 @@ void AtomicDomain::erase(uint64_t pos)

         std::vector<Atom*>::iterator it = std::lower_bound(mAtoms.begin(),

             mAtoms.end(), pos, compareAtomLower);

-        Atom *a = *it;

+        a = *it;

         mAtoms.erase(it);

         GAPS_ASSERT(!vecContains(mAtoms, pos));

@@ -218,7 +248,8 @@ void AtomicDomain::erase(uint64_t pos)

     delete a;

 }

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

+// for moving across a later birth (already present in domain)

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

 {

     #pragma omp critical(AtomicInsertOrErase)

     {

@@ -248,23 +279,7 @@ void AtomicDomain::move(uint64_t src, uint64_t dest)

         GAPS_ASSERT(!vecContains(mAtoms, src));

         GAPS_ASSERT(vecContains(mAtoms, dest));

     }

-}

-

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

-{

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

-    std::vector<Atom*>::iterator it;

-    #pragma omp critical(AtomicInsertOrErase)

-    {

-        GAPS_ASSERT(!vecContains(mAtoms, pos));

-

-        it = std::lower_bound(mAtoms.begin(), mAtoms.end(), pos, compareAtomLower);

-        it = mAtoms.insert(it, newAtom);

-

-        GAPS_ASSERT(vecContains(mAtoms, pos));

-    }

-    return *it;

-}

+}*/

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

 {

@@ -1,8 +1,9 @@

 #ifndef __COGAPS_ATOMIC_DOMAIN_H__

 #define __COGAPS_ATOMIC_DOMAIN_H__

-#include "utils/Archive.h"

+#include "data_structures/HashSets.h"

 #include "math/Random.h"

+#include "utils/Archive.h"

 #include <vector>

@@ -44,9 +45,9 @@ public:

     // access atoms

     Atom* front();

-    Atom* randomAtom(GapsRng *rng);

+    Atom* randomAtom(GapsRng *rng, const SmallPairedHashSetU64 &moves);

     AtomNeighborhood randomAtomWithNeighbors(GapsRng *rng);

-    AtomNeighborhood randomAtomWithRightNeighbor(GapsRng *rng);

+    AtomNeighborhood randomAtomWithRightNeighbor(GapsRng *rng, const SmallPairedHashSetU64 &moves);

     Atom* getLeftNeighbor(uint64_t pos);

     Atom* getRightNeighbor(uint64_t pos);

@@ -57,7 +58,7 @@ public:

     // these need to happen concurrently without invalidating pointers

     void erase(uint64_t pos);

-    void move(uint64_t src, uint64_t dest);

+    //void move(uint64_t src, uint64_t dest);

     // iterators

     std::vector<Atom*>::iterator begin() { return mAtoms.begin(); }

@@ -94,6 +94,10 @@ GapsResult GapsRunner::run(bool printThreads)

     }

     #endif

+#ifdef __GAPS_DEBUG_NO_QUEUE__

+    gaps_printf("Running without a queue\n");

+#endif

+

     // cascade through phases, allows algorithm to be resumed in either phase

     GAPS_ASSERT(mPhase == 'C' || mPhase == 'S');

     switch (mPhase)

@@ -23,8 +23,9 @@ void GibbsSampler::setSparsity(float alpha, bool singleCell)

         gaps::algo::mean(mDMatrix);

     mAlpha = alpha;

-    mQueue.setAlpha(alpha);

     mLambda = alpha * std::sqrt(mNumPatterns / meanD);

+    mQueue.setAlpha(alpha);

+    mQueue.setLambda(mLambda);

 }

 void GibbsSampler::setMaxGibbsMass(float max)

@@ -72,8 +73,11 @@ void GibbsSampler::update(unsigned nSteps, unsigned nCores)

     while (n < nSteps)

     {

         // populate queue, prepare domain for this queue

+        #ifdef __GAPS_DEBUG_NO_QUEUE__

+        mQueue.populate(mDomain, 1);

+        #else

         mQueue.populate(mDomain, nSteps - n);

-        //mQueue.populate(mDomain, 1);

+        #endif

         n += mQueue.size();

         // process all proposed updates

@@ -117,7 +121,6 @@ void GibbsSampler::birth(const AtomicProposal &prop)

         : prop.rng.exponential(mLambda);

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

-    //gaps_printf("mass: %.16f, accept: %d\n", mass, mass >= gaps::epsilon);

     if (mass >= gaps::epsilon)

     {

         mQueue.acceptBirth();

@@ -152,7 +155,6 @@ void GibbsSampler::death(const AtomicProposal &prop)

     // accept/reject rebirth

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

-    //gaps_printf("deltaLL: %.16f\n", deltaLL);

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

     {

         mQueue.rejectDeath();

@@ -176,7 +178,7 @@ void GibbsSampler::move(const AtomicProposal &prop)

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

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

     {

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

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

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

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

     }

@@ -1,4 +1,4 @@

-PKG_CPPFLAGS = @GAPS_CPP_FLAGS@ 

+PKG_CPPFLAGS = @GAPS_CPP_FLAGS@ -D__GAPS_DEBUG_NO_QUEUE__

 PKG_CXXFLAGS = @GAPS_CXX_FLAGS@

 PKG_LIBS = @GAPS_LIBS@

@@ -1,4 +1,5 @@

 #include "ProposalQueue.h"

+#include "math/Algorithms.h"

 #include "math/Random.h"

 #include "utils/GapsAssert.h"

@@ -30,11 +31,17 @@ void ProposalQueue::setAlpha(float alpha)

     mAlpha = static_cast<double>(alpha);

 }

+void ProposalQueue::setLambda(float lambda)

+{

+    mLambda = lambda;

+}

+

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

 {

     GAPS_ASSERT(mQueue.empty());

     GAPS_ASSERT(mUsedAtoms.isEmpty());

     GAPS_ASSERT(mUsedMatrixIndices.isEmpty());

+    GAPS_ASSERT(mProposedMoves.isEmpty());

     GAPS_ASSERT(mMinAtoms == mMaxAtoms);

     GAPS_ASSERT(mMaxAtoms == domain.size());

@@ -47,10 +54,6 @@ void ProposalQueue::populate(AtomicDomain &domain, unsigned limit)

             success = false;

             mUseCachedRng = true;

         }

-        else

-        {

-            //gaps_printf("%c %llu\n", mQueue.back().type, mQueue.back().atom1->pos);

-        }

     }

 }

@@ -61,6 +64,7 @@ void ProposalQueue::clear()

     mQueue.clear();

     mUsedMatrixIndices.clear();

     mUsedAtoms.clear();

+    mProposedMoves.clear();

 }

 unsigned ProposalQueue::size() const

@@ -110,7 +114,6 @@ bool ProposalQueue::makeProposal(AtomicDomain &domain)

 {

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

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

-    //gaps_printf("u1=%f, u2=%f\n", mU1, mU2);

     mUseCachedRng = false;

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

@@ -138,15 +141,14 @@ bool ProposalQueue::makeProposal(AtomicDomain &domain)

         }

         return false; // can't determine B/D since range is too wide

     }

-    //return false;

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

-    return exchange(domain);

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

 }

 bool ProposalQueue::birth(AtomicDomain &domain)

 {

     AtomicProposal prop('B');

     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

@@ -154,6 +156,12 @@ bool ProposalQueue::birth(AtomicDomain &domain)

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

     }

+    if (mProposedMoves.overlap(pos))

+    {

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

+        return false; // this birth would break assumption moves doesn't re-order domain

+    }

+

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

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

     if (mUsedMatrixIndices.contains(prop.r1))

@@ -173,7 +181,14 @@ bool ProposalQueue::birth(AtomicDomain &domain)

 bool ProposalQueue::death(AtomicDomain &domain)

 {

     AtomicProposal prop('D');

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

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

+

+    if (prop.atom1 == NULL)

+    {

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

+        return false; // can't select atom

+    }

+        

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

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

@@ -196,7 +211,8 @@ bool ProposalQueue::move(AtomicDomain &domain)

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

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

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

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

+

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

     {

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

@@ -218,25 +234,31 @@ bool ProposalQueue::move(AtomicDomain &domain)

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

     {

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

-        return true;

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

+        return true; // automatically accept moves in same bin

     }

     mQueue.push_back(prop);

     mUsedMatrixIndices.insert(prop.r1);

     mUsedMatrixIndices.insert(prop.r2);

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

-    mUsedAtoms.insert(prop.pos);

+    mProposedMoves.insert(prop.atom1->pos, prop.pos);

     return true;

 }

 bool ProposalQueue::exchange(AtomicDomain &domain)

 {

     AtomicProposal prop('E');

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

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

+

+    if (hood.center == NULL)

+    {

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

+        return false; // can't select atom

+    }

+

     prop.atom1 = hood.center;

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

-

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

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

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

@@ -250,14 +272,14 @@ bool ProposalQueue::exchange(AtomicDomain &domain)

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

     {

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

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

+        float delta = (prop.atom1->mass > prop.atom2->mass) ? newMass - prop.atom1->mass : prop.atom2->mass - newMass;

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

+        {

+            prop.atom1->mass += delta;

+            prop.atom2->mass -= delta;

+        }        

+        return true;

     }

     mQueue.push_back(prop);

@@ -35,6 +35,7 @@ public:

     // initialize

     ProposalQueue(unsigned nrow, unsigned ncol);

     void setAlpha(float alpha);

+    void setLambda(float lambda);

     // modify/access queue

     void populate(AtomicDomain &domain, unsigned limit);

@@ -58,6 +59,7 @@ private:

     FixedHashSetU32 mUsedMatrixIndices;

     SmallHashSetU64 mUsedAtoms;

+    SmallPairedHashSetU64 mProposedMoves;

     mutable GapsRng mRng;

@@ -70,6 +72,7 @@ private:

     double mDomainLength; // length of entire atomic domain

     double mNumBins; // number of matrix elements

+    float mLambda;

     float mU1;

     float mU2;

@@ -70,3 +70,49 @@ const std::vector<uint64_t>& SmallHashSetU64::vec()

 {

     return mSet;

 }

+

+///////////////////////////// SmallPairedHashSetU64 ////////////////////////////

+

+SmallPairedHashSetU64::SmallPairedHashSetU64() {}

+

+

+void SmallPairedHashSetU64::insert(uint64_t a, uint64_t b)

+{

+    mSet.push_back(a < b ? PositionPair(a, b) : PositionPair(b, a));

+}

+

+void SmallPairedHashSetU64::clear()

+{

+    mSet.clear();

+}

+

+bool SmallPairedHashSetU64::overlap(uint64_t pos)

+{

+    unsigned sz = mSet.size();

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

+    {

+        if (mSet[i].a < pos && pos < mSet[i].b)

+        {

+            return true;

+        }

+    }

+    return false;

+}

+

+bool SmallPairedHashSetU64::contains(uint64_t pos) const

+{

+    unsigned sz = mSet.size();

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

+    {

+        if (mSet[i].a == pos || mSet[i].b == pos)

+        {

+            return true;

+        }

+    }

+    return false;

+}

+

+bool SmallPairedHashSetU64::isEmpty()

+{

+    return mSet.empty();

+}

\ No newline at end of file

@@ -8,7 +8,7 @@

 #include <boost/unordered_set.hpp>

-#ifdef __GAPS_OPENMP__

+#ifdef __GAPS_OPENMP__ // defined in global config

 #include <omp.h>

 #endif

@@ -46,4 +46,29 @@ private:

     std::vector<uint64_t> mSet;

 };

+struct PositionPair

+{

+    uint64_t a;

+    uint64_t b;

+

+    PositionPair(uint64_t inA, uint64_t inB) : a(inA), b(inB) {}

+};

+

+class SmallPairedHashSetU64

+{

+public:

+

+    SmallPairedHashSetU64();

+

+    void insert(uint64_t a, uint64_t b);

+    void clear();

+    bool contains(uint64_t pos) const;

+    bool overlap(uint64_t pos);

+    bool isEmpty();

+

+private:

+

+    std::vector<PositionPair> mSet;

+};

+

 #endif // __COGAPS_HASH_SETS_H__

@@ -22,9 +22,9 @@ static Xoroshiro128plus seeder;

 /////////////////////////////// OptionalFloat //////////////////////////////////

-OptionalFloat::OptionalFloat() : mHasValue(false), mValue(0.f) {}

+OptionalFloat::OptionalFloat() : mValue(0.f), mHasValue(false) {}

-OptionalFloat::OptionalFloat(float f) : mHasValue(true), mValue(f) {}

+OptionalFloat::OptionalFloat(float f) : mValue(f), mHasValue(true) {}

 float OptionalFloat::value()

 {