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

 #ifndef __COGAPS_ARCHIVE_H__

 #define __COGAPS_ARCHIVE_H__

+#include "GapsAssert.h"

+

 #include <fstream>

 #include <stdint.h>

@@ -9,8 +11,8 @@

 #define ARCHIVE_WRITE (std::ios::out | std::ios::trunc)

 // magic number written to beginning of archive files

-// needs to be updated everytime to method of checkpointing changes

-#define ARCHIVE_MAGIC_NUM 0xCE45D32A

+// needs to be updated everytime the method of checkpointing changes

+#define ARCHIVE_MAGIC_NUM 0xCE45D32B // v3.3.22

 class Archive

 {

@@ -21,10 +23,28 @@ private:

 public:

     Archive(const std::string &path, std::ios_base::openmode flags)

-        : mStream(path.c_str(), std::ios::binary | flags)

-    {}

+        :

+    mStream(path.c_str(), std::ios::binary | flags)

+    {

+        if (flags == ARCHIVE_WRITE)

+        {

+            *this << static_cast<uint32_t>(ARCHIVE_MAGIC_NUM);

+        }

+        else // read

+        {

+            uint32_t magic = 0;

+            *this >> magic;

+            if (magic != ARCHIVE_MAGIC_NUM)

+            {

+                GAPS_ERROR("incompatible checkpoint file\n");

+            }            

+        }

+    }

-    void close() {mStream.close();}

+    void close()

+    {

+        mStream.close();

+    }

     template<typename T>

     friend Archive& writeToArchive(Archive &ar, T val)

@@ -3,6 +3,40 @@

 #include <vector>

+////////////////////////////////// HELPER //////////////////////////////////////

+

+// used with std::lower_bound

+static bool compareAtomLower(const Atom &lhs, uint64_t pos)

+{

+    return lhs.pos < pos;

+}

+

+// used with std::binary_search

+static bool compareAtom(const Atom &lhs, const Atom &rhs)

+{

+    return lhs.pos < rhs.pos;

+}

+

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

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

+{

+    Atom temp(pos, 0.f);

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

+}

+

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

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

+{

+    for (unsigned i = 1; i < vec.size(); ++i)

+    {

+        if (vec[i].pos <= vec[i-1].pos)

+        {

+            return false;

+        }

+    }

+    return true;

+}

+

 ////////////////////////////////// ATOM ////////////////////////////////////////

 Atom::Atom() : pos(0), mass(0.f) {}

@@ -46,7 +80,7 @@ bool AtomNeighborhood::hasRight()

 ////////////////////////////// ATOMIC DOMAIN ///////////////////////////////////

-AtomicDomain::AtomicDomain(unsigned nBins)

+AtomicDomain::AtomicDomain(uint64_t nBins)

 {

     uint64_t binLength = std::numeric_limits<uint64_t>::max() / nBins;

     mDomainLength = binLength * nBins;

@@ -55,12 +89,15 @@ AtomicDomain::AtomicDomain(unsigned nBins)

 Atom* AtomicDomain::front()

 {

     GAPS_ASSERT(size() > 0);

+

     return &(mAtoms.front());

 }

 Atom* AtomicDomain::randomAtom()

 {

     GAPS_ASSERT(size() > 0);

+    GAPS_ASSERT(isSorted(mAtoms));

+

     unsigned index = mRng.uniform32(0, mAtoms.size() - 1);

     return &(mAtoms[index]);

 }

@@ -68,6 +105,7 @@ Atom* AtomicDomain::randomAtom()

 AtomNeighborhood AtomicDomain::randomAtomWithNeighbors()

 {

     GAPS_ASSERT(size() > 0);

+

     unsigned index = mRng.uniform32(0, mAtoms.size() - 1);

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

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

@@ -77,32 +115,12 @@ AtomNeighborhood AtomicDomain::randomAtomWithNeighbors()

 AtomNeighborhood AtomicDomain::randomAtomWithRightNeighbor()

 {

     GAPS_ASSERT(size() > 0);

+

     unsigned index = mRng.uniform32(0, mAtoms.size() - 1);

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

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

 }

-static bool compareAtomLower(const Atom &lhs, uint64_t pos)

-{

-    return lhs.pos < pos;

-}

-

-static bool compareAtomUpper(uint64_t pos, const Atom &lhs)

-{

-    return lhs.pos < pos;

-}

-

-static bool compareAtom(const Atom &lhs, const Atom &rhs)

-{

-    return lhs.pos < rhs.pos;

-}

-

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

-{

-    Atom temp(pos, 0.f);

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

-}

-

 uint64_t AtomicDomain::randomFreePosition() const

 {

     uint64_t pos = mRng.uniform64(0, mDomainLength);

@@ -149,6 +167,8 @@ void AtomicDomain::resetCache(unsigned n)

 void AtomicDomain::erase(uint64_t pos)

 {

     GAPS_ASSERT(size() > 0);

+    GAPS_ASSERT(vecContains(mAtoms, pos));

+

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

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

     mAtoms.erase(it);

@@ -191,7 +211,7 @@ Archive& operator<<(Archive &ar, AtomicDomain &domain)

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

 {

     Atom temp;

-    unsigned size = 0;

+    uint64_t size = 0;

     ar >> domain.mDomainLength >> domain.mRng >> size;

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

     {

@@ -36,7 +36,7 @@ class AtomicDomain

 {

 public:

-    AtomicDomain(unsigned nBins);

+    AtomicDomain(uint64_t nBins);

     // access atoms

     Atom* front();

@@ -64,7 +64,7 @@ private:

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

     uint64_t mDomainLength;

-    // domain storage, specialized hash map

+    // domain storage, sorted vector

     std::vector<Atom> mAtoms;

     // holds cache of operations

@@ -48,4 +48,4 @@

     #define DEBUG_PING   

 #endif

-#endif

\ No newline at end of file

+#endif // __COGAPS_GAPS_ASSERT_H__

\ No newline at end of file

@@ -253,37 +253,6 @@ void GibbsSampler<T, MatA, MatB>::setMatrix(const Matrix &mat)

     mMatrix = mat;

 }

-template <class T, class MatA, class MatB>

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

-{

-    unsigned n = 0;

-    while (n < nSteps)

-    {

-        // populate queue, prepare domain for this queue

-        //mQueue.populate(mDomain, nSteps - n);

-        mQueue.populate(mDomain, 1);

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

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

-        }

-

-        mDomain.flushCache();

-        mQueue.clear();

-    }

-}

-

 template <class T, class MatA, class MatB>

 unsigned GibbsSampler<T, MatA, MatB>::dataRows() const

 {

@@ -308,50 +277,45 @@ uint64_t GibbsSampler<T, MatA, MatB>::nAtoms() const

     return mDomain.size();

 }

-#ifdef GAPS_DEBUG

 template <class T, class MatA, class MatB>

-float GibbsSampler<T, MatA, MatB>::getAvgQueue() const

-{

-    return mAvgQueue;

-}

-

-template <class T, class MatA, class MatB>

-bool GibbsSampler<T, MatA, MatB>::internallyConsistent()

+T* GibbsSampler<T, MatA, MatB>::impl()

 {

-    return true;

+    return static_cast<T*>(this);

 }

-#endif

 template <class T, class MatA, class MatB>

-Archive& operator<<(Archive &ar, GibbsSampler<T, MatA, MatB> &sampler)

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

 {

-    ar << sampler.mMatrix << sampler.mAPMatrix << sampler.mQueue <<

-        sampler.mDomain << sampler.mLambda << sampler.mMaxGibbsMass <<

-        sampler.mAnnealingTemp << sampler.mNumRows << sampler.mNumCols <<

-        sampler.mBinSize << sampler.mAvgQueue << sampler.mNumQueues;

-    return ar;

-}

+    unsigned n = 0;

+    while (n < nSteps)

+    {

+        // populate queue, prepare domain for this queue

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

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

+        n += mQueue.size();

+        

+        // update average queue count

+        #ifdef GAPS_DEBUG

+        mNumQueues += 1.f;

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

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

+        #endif

-template <class T, class MatA, class MatB>

-Archive& operator>>(Archive &ar, GibbsSampler<T, MatA, MatB> &sampler)

-{

-    ar >> sampler.mMatrix >> sampler.mAPMatrix >> sampler.mQueue >>

-        sampler.mDomain >> sampler.mLambda >> sampler.mMaxGibbsMass >>

-        sampler.mAnnealingTemp >> sampler.mNumRows >> sampler.mNumCols >>

-        sampler.mBinSize >> sampler.mAvgQueue >> sampler.mNumQueues;

-    return ar;

-}

+        // process all proposed updates

+        #pragma omp parallel for num_threads(nCores)

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

+        {

+            processProposal(&mQueue[i]);

+        }

-template <class T, class MatA, class MatB>

-T* GibbsSampler<T, MatA, MatB>::impl()

-{

-    return static_cast<T*>(this);

+        mDomain.flushCache();

+        mQueue.clear();

+    }

 }

 template <class T, class MatA, class MatB>

 void GibbsSampler<T, MatA, MatB>::processProposal(AtomicProposal *prop)

 {

-    unsigned r1 = 0, c1 = 0, r2 = 0, c2 = 0;

     switch (prop->type)

     {

         case 'B':

@@ -413,8 +377,9 @@ void GibbsSampler<T, MatA, MatB>::death(AtomicProposal *prop)

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

     // kill off atom

-    mMatrix(row, col) = gaps::max(mMatrix(row, col) - prop->atom1->mass, 0.f);

-    impl()->updateAPMatrix(row, col, -1.f * prop->atom1->mass);

+    float newVal = gaps::max(mMatrix(row, col) - prop->atom1->mass, 0.f);

+    impl()->updateAPMatrix(row, col, newVal - mMatrix(row, col));

+    mMatrix(row, col) = newVal;

     // calculate rebirth mass

     float rebirthMass = prop->atom1->mass;

@@ -459,13 +424,13 @@ void GibbsSampler<T, MatA, MatB>::move(AtomicProposal *prop)

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

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

     {

-        //mDomain.cacheMove(prop->atom1->pos, prop->moveDest, prop->atom1->mass);

-        prop->atom1->pos = prop->moveDest; // temporary for vector domain

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

-        mMatrix(r1, c1) = gaps::max(mMatrix(r1, c1) - prop->atom1->mass, 0.f);

-        mMatrix(r2, c2) += prop->atom1->mass;

+        float newVal = gaps::max(mMatrix(r1, c1) - prop->atom1->mass, 0.f);

+        impl()->updateAPMatrix(r1, c1, newVal - mMatrix(r1, c1));

+        mMatrix(r1, c1) = newVal;

-        impl()->updateAPMatrix(r1, c1, -1.f * prop->atom1->mass);

+        mMatrix(r2, c2) += prop->atom1->mass;

         impl()->updateAPMatrix(r2, c2, prop->atom1->mass);

     }

 }

@@ -565,10 +530,15 @@ float d1, unsigned r1, unsigned c1, unsigned r2, unsigned c2)

         mQueue.rejectDeath();

     }

-    mMatrix(r1, c1) += d1;

-    mMatrix(r2, c2) += d2;

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

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

+    // ensure matrix values don't go negative (truncation error at fault)

+    float v1 = gaps::max(mMatrix(r1, c1) + d1, 0.f);

+    impl()->updateAPMatrix(r1, c1, v1 - mMatrix(r1, c1));

+    mMatrix(r1, c1) = v1;

+

+

+    float v2 = gaps::max(mMatrix(r2, c2) + d2, 0.f);

+    impl()->updateAPMatrix(r2, c2, v2 - mMatrix(r2, c2));

+    mMatrix(r2, c2) = v2;

 }

 template <class T, class MatA, class MatB>

@@ -618,4 +588,36 @@ float m1, float m2, GapsRng *rng)

     return std::pair<float, bool>(0.f, false);

 }

-#endif

+#ifdef GAPS_DEBUG

+template <class T, class MatA, class MatB>

+float GibbsSampler<T, MatA, MatB>::getAvgQueue() const

+{

+    return mAvgQueue;

+}

+

+template <class T, class MatA, class MatB>

+bool GibbsSampler<T, MatA, MatB>::internallyConsistent()

+{

+    return true;

+}

+#endif // GAPS_DEBUG

+

+template <class T, class MatA, class MatB>

+Archive& operator<<(Archive &ar, GibbsSampler<T, MatA, MatB> &s)

+{

+    ar << s.mMatrix << s.mAPMatrix << s.mQueue << s.mDomain << s.mLambda

+        << s.mMaxGibbsMass << s.mAnnealingTemp << s.mNumRows << s.mNumCols

+        << s.mBinSize << s.mAvgQueue << s.mNumQueues;

+    return ar;

+}

+

+template <class T, class MatA, class MatB>

+Archive& operator>>(Archive &ar, GibbsSampler<T, MatA, MatB> &s)

+{

+    ar >> s.mMatrix >> s.mAPMatrix >> s.mQueue >> s.mDomain >> s.mLambda

+        >> s.mMaxGibbsMass >> s.mAnnealingTemp >> s.mNumRows >> s.mNumCols

+        >> s.mBinSize >> s.mAvgQueue >> s.mNumQueues;

+    return ar;

+}

+

+#endif // __COGAPS_GIBBS_SAMPLER_H__

@@ -2,12 +2,38 @@

 #include "ProposalQueue.h"

 #include "math/Random.h"

+//////////////////////////////// AtomicProposal ////////////////////////////////

+

+// birth

+AtomicProposal::AtomicProposal(char t, uint64_t pos)

+    : type(t), birthPos(pos), moveDest(0), atom1(NULL), atom2(NULL)

+{}

+    

+// death

+AtomicProposal::AtomicProposal(char t, Atom *atom)

+    : type(t), birthPos(0), moveDest(0), atom1(atom), atom2(NULL)

+{}

+

+// move

+AtomicProposal::AtomicProposal(char t, Atom *atom, uint64_t dest)

+    : type(t), birthPos(0), moveDest(dest), atom1(atom), atom2(NULL)

+{}

+

+// exchange

+AtomicProposal::AtomicProposal(char t, Atom *a1, Atom *a2)

+    : type(t), birthPos(0), moveDest(0), atom1(a1), atom2(a2)

+{}

+

+

+//////////////////////////////// ProposalQueue /////////////////////////////////

+

 ProposalQueue::ProposalQueue(unsigned primaryDimSize, unsigned secondaryDimSize)

     :

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

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

 mSecondaryDimLength(mBinLength * secondaryDimSize),

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

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

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

 {

     mUsedIndices.setDimensionSize(primaryDimSize);

 }

@@ -19,6 +45,9 @@ void ProposalQueue::setAlpha(float alpha)

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

 {

+    GAPS_ASSERT(mMinAtoms == mMaxAtoms);

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

+

     unsigned nIter = 0;

     bool success = true;

     while (nIter++ < limit && success)

@@ -33,10 +62,11 @@ void ProposalQueue::populate(AtomicDomain &domain, unsigned limit)

 void ProposalQueue::clear()

 {

+    GAPS_ASSERT(mMinAtoms == mMaxAtoms);

+

     mQueue.clear();

     mUsedPositions.clear();

     mUsedIndices.clear();

-    GAPS_ASSERT(mMinAtoms == mMaxAtoms);

 }

 unsigned ProposalQueue::size() const

@@ -46,6 +76,9 @@ unsigned ProposalQueue::size() const

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

 {

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

+    GAPS_ASSERT(n < mQueue.size());

+

     return mQueue[n];

 }

@@ -130,7 +163,8 @@ unsigned ProposalQueue::secondaryIndex(uint64_t pos) const

 }

 // TODO add atoms with empty mass? fill in mass in gibbssampler?

-// inserting invalidates previous pointers

+// inserting invalidates previous pointers, but not inserting

+// prevents them from being selected for death

 bool ProposalQueue::birth(AtomicDomain &domain)

 {

     uint64_t pos = domain.randomFreePosition();

@@ -198,7 +232,7 @@ bool ProposalQueue::move(AtomicDomain &domain)

 bool ProposalQueue::exchange(AtomicDomain &domain)

 {

     AtomNeighborhood hood = domain.randomAtomWithRightNeighbor();

-    Atom* a1 = hood.center;    

+    Atom* a1 = hood.center;

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

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

@@ -242,20 +276,18 @@ bool ProposalQueue::exchange(AtomicDomain &domain)

     return true;

 }

-Archive& operator<<(Archive &ar, ProposalQueue &queue)

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

 {

-    ar << queue.mMinAtoms << queue.mMaxAtoms << queue.mNumBins

-        << queue.mBinLength << queue.mSecondaryDimLength << queue.mDomainLength

-        << queue.mSecondaryDimSize << queue.mAlpha << queue.mRng

-        << queue.mU1 << queue.mU2 << queue.mUseCachedRng;

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

+        << q.mSecondaryDimLength << q.mDomainLength << q.mSecondaryDimSize

+        << q.mAlpha << q.mRng;

     return ar;

 }

-Archive& operator>>(Archive &ar, ProposalQueue &queue)

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

 {

-    ar >> queue.mMinAtoms >> queue.mMaxAtoms >> queue.mNumBins

-        >> queue.mBinLength >> queue.mSecondaryDimLength >> queue.mDomainLength

-        >> queue.mSecondaryDimSize >> queue.mAlpha >> queue.mRng

-        >> queue.mU1 >> queue.mU2 >> queue.mUseCachedRng;

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

+        >> q.mSecondaryDimLength >> q.mDomainLength >> q.mSecondaryDimSize

+        >> q.mAlpha >> q.mRng;

     return ar;

 }

@@ -21,25 +21,10 @@ struct AtomicProposal

     mutable GapsRng rng;

-    // birth

-    AtomicProposal(char t, uint64_t pos)

-        : type(t), birthPos(pos), moveDest(0), atom1(NULL), atom2(NULL)

-    {}

-        

-    // death

-    AtomicProposal(char t, Atom *atom)

-        : type(t), birthPos(0), moveDest(0), atom1(atom), atom2(NULL)

-    {}

-

-    // move

-    AtomicProposal(char t, Atom *atom, uint64_t dest)

-        : type(t), birthPos(0), moveDest(dest), atom1(atom), atom2(NULL)

-    {}

-

-    // exchange

-    AtomicProposal(char t, Atom *a1, Atom *a2)

-        : type(t), birthPos(0), moveDest(0), atom1(a1), atom2(a2)

-    {}

+    AtomicProposal(char t, uint64_t pos); // birth

+    AtomicProposal(char t, Atom *atom); // death

+    AtomicProposal(char t, Atom *atom, uint64_t dest); // move

+    AtomicProposal(char t, Atom *a1, Atom *a2); // exchange

 };

 class ProposalQueue

@@ -226,8 +226,8 @@ things like hardware and data size. The best approach is to play around with

 different values and see how it effects the estimated time.

 ```{r}

-CoGAPS(SimpSim.data, nIterations=10000, outputFrequency=5000, nThreads=1)

-CoGAPS(SimpSim.data, nIterations=10000, outputFrequency=5000, nThreads=4)

+CoGAPS(SimpSim.data, nIterations=10000, outputFrequency=5000, nThreads=1, seed=5)

+CoGAPS(SimpSim.data, nIterations=10000, outputFrequency=5000, nThreads=4, seed=5)

 ```

 Note this method relies on CoGAPS being compiled with OpenMP support, use