new file mode 100644

@@ -0,0 +1,25 @@

+# All conflicts needed to be considered for a truly equal algo

+

+## Birth

+- (ignore for now) track all atoms that are selected in death/move

+    /exchange, if a position is selected that conflicts with one 

+    of those atoms, we can't proceed since it is possible that 

+    position will be free when the birth happens

+- track previous births to make sure no overwriting happens

+

+## Death

+- (ignore for now) if any atom is birthed before this step then

+    it should be possible to select it

+- track all atoms selected for death/exchange if any atom is

+    destroyed before this step we can't destroy it again

+

+## Move

+- Same selection problem as death

+- can't have any atoms birth/death/move/exchange within the

+    selected atom's boundaries

+

+## Exchange

+- Same selection problem as death

+- Similiar problem as move with selecting exchanged atom

+

+can't ignore birthed atoms, death/move/exchange fail if randomAtom returns a newly birthed atom

\ No newline at end of file

@@ -14,7 +14,10 @@ Atom AtomicDomain::front() const

 // O(1)

 Atom AtomicDomain::randomAtom() const

 {

-    return mAtoms[gaps::random::uniform64(0, mAtoms.size() - 1)];

+    uint64_t ndx = gaps::random::uniform64(0, mAtoms.size() - 1);

+    GAPS_ASSERT(mAtoms.size() >= 1);

+    GAPS_ASSERT(ndx < mAtoms.size());

+    return mAtoms[ndx];

 }

 // Average Case O(1)

@@ -71,7 +74,7 @@ bool AtomicDomain::hasRight(const Atom &atom) const

 }

 // O(logN)

-void AtomicDomain::insert(uint64_t pos, float mass)

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

 {

     // insert position into map

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

@@ -97,6 +100,7 @@ void AtomicDomain::insert(uint64_t pos, float mass)

     // add atom to vector

     mAtoms.push_back(atom);

     mUsedPositions.insert(pos);

+    return atom;

 }

 // O(logN)

@@ -76,7 +76,7 @@ public:

     bool hasRight(const Atom &atom) const;

     // modify domain

-    void insert(uint64_t pos, float mass);

+    Atom insert(uint64_t pos, float mass);

     void erase(uint64_t pos);

     void updateMass(uint64_t pos, float newMass);

@@ -90,7 +90,9 @@ Rcpp::List GapsRunner::run()

         Rcpp::Named("chiSqValues") = chi2Vec.rVec(),

         Rcpp::Named("randSeed") = mSeed,

         Rcpp::Named("numUpdates") = mNumUpdatesA + mNumUpdatesP,

-        Rcpp::Named("meanChi2") = meanChiSq

+        Rcpp::Named("meanChi2") = meanChiSq,

+        Rcpp::Named("AAvgQueue") = mASampler.getAvgQueue(),

+        Rcpp::Named("PAvgQueue") = mPSampler.getAvgQueue()

     );

 }

@@ -8,6 +8,7 @@ const Rcpp::NumericMatrix &S, unsigned nFactor, float alpha, float maxGibbsMass)

     float meanD = gaps::algo::mean(mDMatrix);

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

     mMaxGibbsMass = maxGibbsMass / mLambda;

+    mQueue.setDimensionSize(D.nrow());

 }

 unsigned AmplitudeGibbsSampler::getRow(uint64_t pos) const

@@ -110,6 +111,7 @@ const Rcpp::NumericMatrix &S, unsigned nFactor, float alpha, float maxGibbsMass)

     float meanD = gaps::algo::mean(mDMatrix);

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

     mMaxGibbsMass = maxGibbsMass / mLambda;

+    mQueue.setDimensionSize(D.ncol());

 }

 unsigned PatternGibbsSampler::getRow(uint64_t pos) const

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

 #include <Rcpp.h>

 #include <algorithm>

+#include <omp.h>

 // forward declarations needed for friend classes

 class AmplitudeGibbsSampler;

@@ -49,6 +50,9 @@ protected:

     unsigned mNumCols;

     uint64_t mBinSize;

+    float mAvgQueue;

+    float mNumQueues;

+

     T* impl();

     void processProposal(const AtomicProposal &prop);

@@ -75,6 +79,7 @@ public:

     void update(unsigned nSteps);

     void setAnnealingTemp(float temp);

+    float getAvgQueue() const { return mAvgQueue; }

     float chi2() const;

     uint64_t nAtoms() const;

@@ -148,7 +153,8 @@ const Rcpp::NumericMatrix &S, unsigned nrow, unsigned ncol, float alpha)

     :

 mMatrix(nrow, ncol), mOtherMatrix(NULL), mDMatrix(D), mSMatrix(S),

 mAPMatrix(D.nrow(), D.ncol()), mQueue(nrow * ncol, alpha),

-mAnnealingTemp(0.f), mNumRows(nrow), mNumCols(ncol)

+mAnnealingTemp(0.f), mNumRows(nrow), mNumCols(ncol),

+mAvgQueue(0.f), mNumQueues(0.f)

 {

     mBinSize = std::numeric_limits<uint64_t>::max()

         / static_cast<uint64_t>(mNumRows * mNumCols);

@@ -167,9 +173,23 @@ T* GibbsSampler<T, MatA, MatB>::impl()

 template <class T, class MatA, class MatB>

 void GibbsSampler<T, MatA, MatB>::update(unsigned nSteps)

 {

-    for (unsigned n = 0; n < nSteps; ++n)

+    unsigned n = 0;

+    while (n < nSteps)

     {

-        processProposal(mQueue.makeProposal(mDomain));

+        GAPS_ASSERT(nSteps - (mQueue.size() + n) >= 0);

+        mQueue.populate(mDomain, nSteps - (mQueue.size() + n));

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

+        GAPS_ASSERT(mQueue.size() > 0);

+        mNumQueues += 1.f;

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

+        n += mQueue.size();

+        //#pragma omp parallel for

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

+        {

+            processProposal(mQueue[i]);

+        }

+        mQueue.clear(1);

+        GAPS_ASSERT(n <= nSteps);

     }

 }

@@ -226,7 +246,14 @@ unsigned col)

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

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

     {

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

+        mDomain.updateMass(pos, mass);

+        mMatrix(row, col) += mass;

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

+    }

+    else

+    {

+        mDomain.erase(pos);

+        mQueue.rejectBirth();

     }

 }

@@ -243,6 +270,11 @@ unsigned col)

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

     {

         addMass(pos, mass, row, col);

+        mQueue.rejectDeath();

+    }

+    else

+    {

+        mQueue.acceptDeath();

     }

 }

@@ -310,10 +342,11 @@ float m2, unsigned r1, unsigned c1, unsigned r2, unsigned c2)

             if (u < deltaLL * mAnnealingTemp)

             {

                 acceptExchange(p1, m1, delta, p2, m2, -delta, r1, c1, r2, c2);

+                return;

             }

-            return;

         }

     }

+    mQueue.rejectDeath();

 }

 // helper function for acceptExchange, used to conditionally update the mass

@@ -325,12 +358,13 @@ float delta)

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

     {

         mDomain.erase(pos);

-        return -mass;

+        mQueue.acceptDeath();

+        return false;

     }

     else

     {

         mDomain.updateMass(pos, mass + delta);

-        return delta;

+        return true;

     }

 }

@@ -341,8 +375,18 @@ void GibbsSampler<T, MatA, MatB>::acceptExchange(uint64_t p1, float m1,

 float d1, uint64_t p2, float m2, float d2, unsigned r1, unsigned c1,

 unsigned r2, unsigned c2)

 {

-    d1 = updateAtomMass(p1, m1, d1);

-    d2 = updateAtomMass(p2, m2, d2);

+    bool b1 = updateAtomMass(p1, m1, d1);

+    bool b2 = updateAtomMass(p2, m2, d2);

+    GAPS_ASSERT(b1 || b2);

+    

+    if (!b1) d1 = -m1;

+    if (!b2) d2 = -m2;

+

+    if (b1 && b2)

+    {

+        mQueue.rejectDeath();

+    }

+

     mMatrix(r1, c1) += d1;

     mMatrix(r2, c2) += d2;

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

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

-PKG_CPPFLAGS =  -DSIMD -msse4.2 -DBOOST_MATH_PROMOTE_DOUBLE_POLICY=0

+PKG_CPPFLAGS =  -DSIMD -DGAPS_DEBUG -msse4.2 -DBOOST_MATH_PROMOTE_DOUBLE_POLICY=0 -fopenmp

 OBJECTS =   AtomicDomain.o \

             Cogaps.o \

             GapsRunner.o \

@@ -12,11 +12,60 @@ void ProposalQueue::setDomainSize(uint64_t size)

     mDomainSize = size;

 }

+void ProposalQueue::setDimensionSize(unsigned size)

+{

+    mDimensionSize = size;

+}

+

 void ProposalQueue::setAlpha(float alpha)

 {

     mAlpha = alpha;

 }

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

+{

+    unsigned nIter = 0;

+    while (nIter++ < limit && makeProposal(domain));

+}

+

+void ProposalQueue::clear(unsigned n)

+{

+    //mQueue.erase(mQueue.end() - n, mQueue.end());

+    //mUsedIndices.erase(mUsedIndices.end() - n, mUsedIndices.end());

+    //mUsedPositions.erase(mUsedPositions.end() - n, mUsedPositions.end());

+    mQueue.clear();

+    mUsedPositions.clear();

+    mUsedIndices.clear();

+    GAPS_ASSERT(mMaxAtoms - mMinAtoms <= mQueue.size());

+}

+

+unsigned ProposalQueue::size() const

+{

+    return mQueue.size();

+}

+

+const AtomicProposal& ProposalQueue::operator[](int n) const

+{

+    //return mQueue[mQueue.size() - 1 - n];

+    return mQueue[n];

+}

+

+void ProposalQueue::acceptDeath()

+{

+    mMaxAtoms--;

+}

+

+void ProposalQueue::rejectDeath()

+{

+    mMinAtoms++;

+}

+

+void ProposalQueue::rejectBirth()

+{

+    mMinAtoms--;

+    mMaxAtoms--;

+}

+

 float ProposalQueue::deathProb(uint64_t nAtoms) const

 {

     double size = static_cast<double>(mDomainSize);

@@ -25,55 +74,162 @@ float ProposalQueue::deathProb(uint64_t nAtoms) const

     return static_cast<double>(nAtoms) / (static_cast<double>(nAtoms) + term2);

 }

-AtomicProposal ProposalQueue::makeProposal(const AtomicDomain &domain)

+bool ProposalQueue::makeProposal(AtomicDomain &domain)

 {

-    // always birth when no atoms exist

-    if (domain.size() == 0)

+    // special indeterminate cases

+    if ((mMinAtoms == 0 && mMaxAtoms > 0) || (mMinAtoms < 2 && mMaxAtoms >= 2))

     {

-        return birth(domain);

+        return false;

     }

-    float bdProb = domain.size() < 2 ? 0.6667f : 0.5f;

-    float u = gaps::random::uniform();

-    if (u <= bdProb)

+    // always birth when no atoms exist

+    if (mMaxAtoms == 0)

     {

-        return gaps::random::uniform() < deathProb(domain.size()) ?

-            death(domain) : birth(domain);

+        return birth(domain);

     }

-    else if (u < 0.75f || domain.size() < 2)

+

+    float bdProb = mMaxAtoms < 2 ? 0.6667f : 0.5f;

+    float u1 = gaps::random::uniform(); // cache these values if a failure

+    float lowerBound = deathProb(mMinAtoms);

+    float upperBound = deathProb(mMaxAtoms);

+

+    if (u1 <= bdProb)

     {

-        return move(domain);

+        float u2 = gaps::random::uniform(); // happens, needed to prevent bias

+        if (u2 >= upperBound)

+        {

+            return birth(domain);

+        }

+        if (u2 < lowerBound)

+        {

+            return death(domain);

+        }

+        return false;

     }

-    else

+    else if (u1 >= bdProb)

     {

-        return exchange(domain);

+        return (u1 < 0.75f || mMaxAtoms < 2) ? move(domain) : exchange(domain);

     }

+    return false;

 }

-AtomicProposal ProposalQueue::birth(const AtomicDomain &domain)

+static bool isInVector(const std::vector<uint64_t> &vec, uint64_t n)

+{

+    return std::find(vec.begin(), vec.end(), n) != vec.end();

+}

+

+bool ProposalQueue::birth(AtomicDomain &domain)

 {

     uint64_t pos = domain.randomFreePosition();

-    return AtomicProposal('B', pos);

+    if (mUsedIndices.count(pos / mDimensionSize))

+    {

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

+    }

+

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

+    mUsedIndices.insert(pos / mDimensionSize);

+    mUsedPositions.insert(pos);

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

+    ++mMinAtoms;

+    ++mMaxAtoms;

+    return true;

 }

-AtomicProposal ProposalQueue::death(const AtomicDomain &domain)

+bool ProposalQueue::death(AtomicDomain &domain)

 {

     Atom a = domain.randomAtom();

-    return AtomicProposal('D', a.pos, a.mass);

+    if (mUsedIndices.count(a.pos / mDimensionSize))

+    {

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

+    }

+

+    mQueue.push_back(AtomicProposal('D', a.pos, a.mass));

+    mUsedIndices.insert(a.pos / mDimensionSize);

+    mUsedPositions.insert(a.pos);

+    --mMinAtoms;

+    return true;

 }

-AtomicProposal ProposalQueue::move(const AtomicDomain &domain)

+bool ProposalQueue::move(AtomicDomain &domain)

 {

     Atom a = domain.randomAtom();

     uint64_t lbound = domain.hasLeft(a) ? domain.left(a).pos : 0;

     uint64_t rbound = domain.hasRight(a) ? domain.right(a).pos : mDomainSize;

+

+    uint64_t low = *(mUsedPositions.lower_bound(lbound));

+    if (low != *(mUsedPositions.end()) && low <= rbound)

+    {

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

+    }

+

     uint64_t newLocation = gaps::random::uniform64(lbound, rbound - 1);

-    return AtomicProposal('M', a.pos, a.mass, newLocation);

+    if (mUsedIndices.count(a.pos / mDimensionSize) || mUsedIndices.count(newLocation / mDimensionSize))

+    {

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

+    }

+

+    mQueue.push_back(AtomicProposal('M', a.pos, a.mass, newLocation));

+    mUsedIndices.insert(a.pos / mDimensionSize);

+    mUsedIndices.insert(newLocation / mDimensionSize);

+    mUsedPositions.insert(a.pos);

+    mUsedPositions.insert(newLocation);

+    return true;

 }

-AtomicProposal ProposalQueue::exchange(const AtomicDomain &domain)

+bool ProposalQueue::exchange(AtomicDomain &domain)

 {

     Atom a1 = domain.randomAtom();

     Atom a2 = domain.hasRight(a1) ? domain.right(a1) : domain.front();

-    return AtomicProposal('E', a1.pos, a1.mass, a2.pos, a2.mass);

+

+    //if (domain.hasRight(a1))

+    //{

+    //    if (*(mUsedPositions.lower_bound(a1.pos)) <= a2.pos)

+    //    {

+    //        return false;

+    //    }

+    //}

+    //else

+    //{

+    //    if (*(mUsedPositions.end()) >= a1.pos || *(mUsedPositions.begin()) <= a2.pos)

+    //    {

+    //        return false;

+    //    }

+    //}

+

+    if (domain.hasRight(a1)) // has neighbor

+    {

+        std::set<uint64_t>::iterator low = mUsedPositions.lower_bound(a1.pos);

+        if (low != mUsedPositions.end() && *low <= a2.pos)

+        {

+            return false;

+        }

+    }

+    else // exchange with first atom

+    {

+        for (std::set<uint64_t>::iterator it = mUsedPositions.begin(); it != mUsedPositions.end(); ++it)

+        {

+            if (*it >= a1.pos || *it <= domain.front().pos)

+            {

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

+            }

+        }

+        //std::set<uint64_t>::iterator it = mUsedPositions.upper_bound(a1.pos);

+        //if (it != mUsedPositions.end() || *(mUsedPositions.begin()) <= a2.pos)

+        //{

+        //    return false;

+        //}

+    }

+

+    if (mUsedIndices.count(a1.pos / mDimensionSize) || mUsedIndices.count(a2.pos / mDimensionSize))

+    {

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

+    }

+

+    mQueue.push_back(AtomicProposal('E', a1.pos, a1.mass, a2.pos, a2.mass));

+    mUsedIndices.insert(a1.pos / mDimensionSize);

+    mUsedIndices.insert(a2.pos / mDimensionSize);

+    mUsedPositions.insert(a1.pos);

+    mUsedPositions.insert(a2.pos);

+    --mMinAtoms;

+    return true;

 }

\ No newline at end of file

@@ -4,6 +4,7 @@

 #include "Archive.h"

 #include "AtomicDomain.h"

+#include <boost/unordered_set.hpp>

 #include <stdint.h>

 #include <cstddef>

@@ -25,29 +26,50 @@ class ProposalQueue

 {

 private:

+    std::vector<AtomicProposal> mQueue; // not really a queue for now

+    

+    boost::unordered_set<uint64_t> mUsedIndices; // used rows/cols for A/P matrix

+    std::set<uint64_t> mUsedPositions; // used positions in atomic domain

+

+    uint64_t mMinAtoms;

+    uint64_t mMaxAtoms;

+

     uint64_t mNumBins;

+    uint64_t mDimensionSize; // rows of A, cols of P

     uint64_t mDomainSize;

+

     float mAlpha;

     float deathProb(uint64_t nAtoms) const;

-    AtomicProposal birth(const AtomicDomain &domain);

-    AtomicProposal death(const AtomicDomain &domain);

-    AtomicProposal move(const AtomicDomain &domain);

-    AtomicProposal exchange(const AtomicDomain &domain);

+    bool birth(AtomicDomain &domain);

+    bool death(AtomicDomain &domain);

+    bool move(AtomicDomain &domain);

+    bool exchange(AtomicDomain &domain);

+

+    bool makeProposal(AtomicDomain &domain);

 public:

     ProposalQueue(unsigned nBins, float alpha)

-        : mNumBins(nBins), mAlpha(alpha)

+        : mMinAtoms(0), mMaxAtoms(0), mNumBins(nBins), mAlpha(alpha)

     {}

     // set parameters

     void setNumBins(unsigned nBins);

     void setDomainSize(uint64_t size);

     void setAlpha(float alpha);

+    void setDimensionSize(unsigned nIndices);

     // modify/access queue

-    AtomicProposal makeProposal(const AtomicDomain &domain);

+    void populate(AtomicDomain &domain, unsigned limit);

+    void clear(unsigned n);

+    unsigned size() const;

+    const AtomicProposal& operator[](int n) const;

+

+    // update min/max atoms

+    void acceptDeath();

+    void rejectDeath();

+    void rejectBirth();

     // serialization

     friend Archive& operator<<(Archive &ar, ProposalQueue &queue);

@@ -5,6 +5,10 @@

     #define __x86_64__ 1

 #endif

+#if !defined(_OPENMP)

+    #warning "Compiler does not support OpenMP"

+#endif

+

 #ifndef SSE_INSTR_SET

     #ifndef SIMD

         #define SSE_INSTR_SET 0