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

 Package: CoGAPS

-Version: 3.3.24

+Version: 3.3.31

 Date: 2018-04-24

 Title: Coordinated Gene Activity in Pattern Sets

 Author: Thomas Sherman, Wai-shing Lee, Conor Kelton, Ondrej Maxian, Jacob Carey,

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

 #include "AtomicDomain.h"

-#include "GapsAssert.h"

+#include "utils/GapsAssert.h"

 #include <vector>

@@ -136,34 +136,6 @@ uint64_t AtomicDomain::size() const

     return mAtoms.size();

 }

-void AtomicDomain::cacheInsert(uint64_t pos, float mass) const

-{

-    unsigned ndx = 0;

-    #pragma omp critical(atomicInsert)

-    {

-        ndx = mInsertCacheIndex++;

-    }

-    mInsertCache[ndx] = Atom(pos, mass);

-}

-

-void AtomicDomain::cacheErase(uint64_t pos) const

-{

-    unsigned ndx = 0;

-    #pragma omp critical(atomicErase)

-    {

-        ndx = mEraseCacheIndex++;

-    }

-    mEraseCache[ndx] = pos;

-}

-

-void AtomicDomain::resetCache(unsigned n)

-{

-    mInsertCacheIndex = 0;

-    mEraseCacheIndex = 0;

-    mInsertCache.resize(n);

-    mEraseCache.resize(n);

-}

-

 void AtomicDomain::erase(uint64_t pos)

 {

     GAPS_ASSERT(size() > 0);

@@ -181,22 +153,6 @@ void AtomicDomain::insert(uint64_t pos, float mass)

     mAtoms.insert(it, Atom(pos, mass));

 }

-void AtomicDomain::flushCache()

-{

-    for (unsigned i = 0; i < mEraseCacheIndex; ++i)

-    {

-        erase(mEraseCache[i]);

-    }

-

-    for (unsigned i = 0; i < mInsertCacheIndex; ++i)

-    {

-        insert(mInsertCache[i].pos, mInsertCache[i].mass);

-    }

-

-    mInsertCache.clear();

-    mEraseCache.clear();

-}

-

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

 {

     ar << domain.mDomainLength << domain.mRng << domain.mAtoms.size();

@@ -1,7 +1,7 @@

 #ifndef __COGAPS_ATOMIC_DOMAIN_H__

 #define __COGAPS_ATOMIC_DOMAIN_H__

-#include "Archive.h"

+#include "utils/Archive.h"

 #include "math/Random.h"

 #include <vector>

@@ -47,11 +47,8 @@ public:

     uint64_t randomFreePosition() const;

     uint64_t size() const;

-    // modify domain

-    void cacheInsert(uint64_t pos, float mass) const;

-    void cacheErase(uint64_t pos) const;

-    void flushCache();

-    void resetCache(unsigned n);

+    void erase(uint64_t pos);

+    void insert(uint64_t pos, float mass);

     // serialization

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

@@ -67,23 +64,8 @@ private:

     // domain storage, sorted vector

     std::vector<Atom> mAtoms;

-    // holds cache of operations

-    mutable std::vector<Atom> mInsertCache;

-    mutable std::vector<uint64_t> mEraseCache;

-

-    // current index in the operation cache

-    mutable unsigned mInsertCacheIndex;

-    mutable unsigned mEraseCacheIndex;

-

     // random number generator

     mutable GapsRng mRng;

-

-    void erase(uint64_t pos);

-    void insert(uint64_t pos, float mass);

-

-    // serialization

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

-    friend Archive& operator>>(Archive &ar, AtomicDomain &domain);

 };

 #endif

\ No newline at end of file

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

 #include "GapsRunner.h"

 #include "math/SIMD.h"

+#include "utils/GlobalConfig.h"

 #include <Rcpp.h>

 #include <string>

new file mode 100644

@@ -0,0 +1,54 @@

+#include "GapsResult.h"

+

+#include "GapsStatistics.h"

+#include "file_parser/FileParser.h"

+

+GapsResult::GapsResult(const GapsStatistics &stat)

+    :

+Amean(stat.Amean()), Asd(stat.Asd()), Pmean(stat.Pmean()),

+Psd(stat.Psd()), meanChiSq(0.f), seed(0)

+{}

+

+void GapsResult::writeToFile(const std::string &fullPath)

+{

+    std::size_t pos = fullPath.find_last_of('.');

+    std::string base = fullPath.substr(0, pos);

+

+    switch (FileParser::fileType(fullPath))

+    {

+        case GAPS_CSV: return writeCsv(base);

+        case GAPS_TSV: return writeTsv(base);

+        case GAPS_GCT: return writeGct(base);

+        default : GAPS_ERROR("Invalid File Type");

+    }

+}

+

+void GapsResult::writeCsv(const std::string &path)

+{

+    unsigned nPatterns = Amean.nCol();

+    std::string label("_" + gaps::to_string(nPatterns) + "_");

+    FileParser::writeToCsv(path + label + "Amean.csv", Amean);

+    FileParser::writeToCsv(path + label + "Pmean.csv", Pmean);

+    FileParser::writeToCsv(path + label + "Asd.csv", Asd);

+    FileParser::writeToCsv(path + label + "Psd.csv", Psd);

+}

+

+void GapsResult::writeTsv(const std::string &path)

+{

+    unsigned nPatterns = Amean.nCol();

+    std::string label("_" + gaps::to_string(nPatterns) + "_");

+    FileParser::writeToCsv(path + label + "Amean.tsv", Amean);

+    FileParser::writeToCsv(path + label + "Pmean.tsv", Pmean);

+    FileParser::writeToCsv(path + label + "Asd.tsv", Asd);

+    FileParser::writeToCsv(path + label + "Psd.tsv", Psd);

+}

+

+void GapsResult::writeGct(const std::string &path)

+{

+    unsigned nPatterns = Amean.nCol();

+    std::string label("_" + gaps::to_string(nPatterns) + "_");

+    FileParser::writeToCsv(path + label + "Amean.gct", Amean);

+    FileParser::writeToCsv(path + label + "Pmean.gct", Pmean);

+    FileParser::writeToCsv(path + label + "Asd.gct", Asd);

+    FileParser::writeToCsv(path + label + "Psd.gct", Psd);

+}

\ No newline at end of file

new file mode 100644

@@ -0,0 +1,28 @@

+#ifndef __COGAPS_GAPS_RESULT__

+#define __COGAPS_GAPS_RESULT__

+

+#include "GapsStatistics.h"

+#include "data_structures/Matrix.h"

+

+#include <stdint.h>

+#include <string>

+

+struct GapsResult

+{

+    ColMatrix Amean;

+    ColMatrix Asd;

+    ColMatrix Pmean;

+    ColMatrix Psd;

+    

+    float meanChiSq;

+    uint32_t seed;

+

+    GapsResult(const GapsStatistics &stat);

+

+    void writeToFile(const std::string &fullPath);

+    void writeCsv(const std::string &path);

+    void writeTsv(const std::string &path);

+    void writeGct(const std::string &path);

+};

+

+#endif // __COGAPS_GAPS_RESULT__

\ No newline at end of file

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

-#include "GapsAssert.h"

-#include "GapsPrint.h"

 #include "GapsRunner.h"

 #include "math/Random.h"

 #include "math/SIMD.h"

+#include "utils/GapsAssert.h"

+#include "utils/GapsPrint.h"

 #ifdef __GAPS_R_BUILD__

 #include <Rcpp.h>

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

 #ifndef __COGAPS_GAPS_RUNNER_H__

 #define __COGAPS_GAPS_RUNNER_H__

+#include "GapsResult.h"

 #include "GapsStatistics.h"

 #include "GibbsSampler.h"

@@ -11,66 +12,6 @@

 namespace bpt = boost::posix_time;

 #define bpt_now() bpt::microsec_clock::local_time()

-struct GapsResult

-{

-    ColMatrix Amean;

-    ColMatrix Asd;

-    RowMatrix Pmean;

-    RowMatrix Psd;

-    

-    float meanChiSq;

-    uint32_t seed;

-

-    GapsResult(const GapsStatistics &stat) :

-        Amean(stat.Amean()), Asd(stat.Asd()), Pmean(stat.Pmean()),

-        Psd(stat.Psd()), meanChiSq(0.f), seed(0)

-    {}

-

-    void writeToFile(const std::string &fullPath)

-    {

-        std::size_t pos = fullPath.find_last_of('.');

-        std::string base = fullPath.substr(0, pos);

-

-        switch (FileParser::fileType(fullPath))

-        {

-            case GAPS_CSV: return writeCsv(base);

-            case GAPS_TSV: return writeTsv(base);

-            case GAPS_GCT: return writeGct(base);

-            default: GAPS_ERROR("Invalid file type\n");

-        }

-    }

-

-    void writeCsv(const std::string &path)

-    {

-        unsigned nPatterns = Amean.nCol();

-        std::string label("_" + gaps::to_string(nPatterns) + "_");

-        FileParser::writeToCsv(path + label + "Amean.csv", Amean);

-        FileParser::writeToCsv(path + label + "Pmean.csv", Pmean);

-        FileParser::writeToCsv(path + label + "Asd.csv", Asd);

-        FileParser::writeToCsv(path + label + "Psd.csv", Psd);

-    }

-

-    void writeTsv(const std::string &path)

-    {

-        unsigned nPatterns = Amean.nCol();

-        std::string label("_" + gaps::to_string(nPatterns) + "_");

-        FileParser::writeToCsv(path + label + "Amean.tsv", Amean);

-        FileParser::writeToCsv(path + label + "Pmean.tsv", Pmean);

-        FileParser::writeToCsv(path + label + "Asd.tsv", Asd);

-        FileParser::writeToCsv(path + label + "Psd.tsv", Psd);

-    }

-

-    void writeGct(const std::string &path)

-    {

-        unsigned nPatterns = Amean.nCol();

-        std::string label("_" + gaps::to_string(nPatterns) + "_");

-        FileParser::writeToCsv(path + label + "Amean.gct", Amean);

-        FileParser::writeToCsv(path + label + "Pmean.gct", Pmean);

-        FileParser::writeToCsv(path + label + "Asd.gct", Asd);

-        FileParser::writeToCsv(path + label + "Psd.gct", Psd);

-    }

-};

-

 class GapsRunner

 {

 private:

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

 #ifndef __COGAPS_GIBBS_SAMPLER_H__

 #define __COGAPS_GIBBS_SAMPLER_H__

-#include "Archive.h"

+#include "utils/Archive.h"

 #include "AtomicDomain.h"

-#include "GapsAssert.h"

+#include "utils/GapsAssert.h"

 #include "ProposalQueue.h"

 #include "data_structures/Matrix.h"

 #include "math/Algorithms.h"

@@ -43,28 +43,33 @@ protected:

     MatA mMatrix;

     MatB* mOtherMatrix;

-    ProposalQueue mQueue;

+    GapsRng mPropRng;

     AtomicDomain mDomain;

+    float mAlpha;

     float mLambda;

     float mMaxGibbsMass;

     float mAnnealingTemp;

     unsigned mNumRows;

     unsigned mNumCols;

+

+    uint64_t mNumBins;

     uint64_t mBinSize;

+    uint64_t mDomainLength;

     float mAvgQueue;

     float mNumQueues;

     T* impl();

-    void processProposal(AtomicProposal *prop);

+    float deathProb(uint64_t nAtoms) const;

+    void makeAndProcessProposal();

-    void birth(AtomicProposal *prop);

-    void death(AtomicProposal *prop);

-    void move(AtomicProposal *prop);

-    void exchange(AtomicProposal *prop);

+    void birth(GapsRng *rng);

+    void death(GapsRng *rng);

+    void move(GapsRng *rng);

+    void exchange(GapsRng *rng);

     bool updateAtomMass(Atom *atom, float delta);

     void acceptExchange(AtomicProposal *prop, float d1, unsigned r1,

@@ -199,15 +204,19 @@ mDMatrix(data, transposeData, partitionRows, indices),

 mSMatrix(mDMatrix.pmax(0.1f, 0.1f)),

 mAPMatrix(mDMatrix.nRow(), mDMatrix.nCol()),

 mMatrix(amp ? mDMatrix.nRow() : nPatterns, amp ? nPatterns : mDMatrix.nCol()),

-mOtherMatrix(NULL), mQueue(amp ? mDMatrix.nRow() : mDMatrix.nCol(), nPatterns),

-mDomain(mMatrix.nRow() * mMatrix.nCol()), mLambda(0.f),

-mMaxGibbsMass(100.f), mAnnealingTemp(1.f), mNumRows(mMatrix.nRow()),

-mNumCols(mMatrix.nCol()), mAvgQueue(0.f), mNumQueues(0.f)

+mOtherMatrix(NULL),

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

+mLambda(0.f),

+mMaxGibbsMass(100.f),

+mAnnealingTemp(1.f),

+mNumRows(mMatrix.nRow()),

+mNumCols(mMatrix.nCol()),

+mAvgQueue(0.f),

+mNumQueues(0.f),

+mNumBins(mMatrix.nRow() * mMatrix.nCol()),

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

+mDomainLength(mBinSize * mNumBins)

 {

-    // calculate atomic domain size

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

-        / static_cast<uint64_t>(mNumRows * mNumCols);

-

     // default sparsity parameters

     setSparsity(0.01, false);

 }

@@ -223,14 +232,13 @@ bool transpose, bool partitionRows, const std::vector<unsigned> &indices)

 template <class T, class MatA, class MatB>

 void GibbsSampler<T, MatA, MatB>::setSparsity(float alpha, bool singleCell)

 {

-    mQueue.setAlpha(alpha);

-

     float meanD = singleCell ? gaps::algo::nonZeroMean(mDMatrix) :

         gaps::algo::mean(mDMatrix);

     unsigned nPatterns = mDMatrix.nRow() == mMatrix.nRow() ? mMatrix.nCol() :

         mMatrix.nRow();

+    mAlpha = alpha;

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

 }

@@ -288,104 +296,93 @@ 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);

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

+        makeAndProcessProposal();

+        ++n;

     }

 }

 template <class T, class MatA, class MatB>

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

+float GibbsSampler<T, MatA, MatB>::deathProb(uint64_t nAtoms) const

 {

-    switch (prop->type)

+    double size = static_cast<double>(mDomainLength);

+    double term1 = (size - static_cast<double>(nAtoms)) / size;

+    double term2 = mAlpha * static_cast<double>(mNumBins) * term1;

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

+}

+

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

+void GibbsSampler<T, MatA, MatB>::makeAndProcessProposal()

+{

+    GapsRng rng; // need to make this independent of nThreads (use streams?)

+    if (mDomain.size() < 2) // always birth with 0 or 1 atoms

     {

-        case 'B':

-            birth(prop);

-            break;

-        case 'D':

-            death(prop);

-            break;

-        case 'M':

-            move(prop);

-            break;

-        case 'E':

-            exchange(prop);

-            break;

+        return birth(&rng);

     }

+

+    //float u1 = rng.uniform();

+    float u1 = 0.f;

+    if (u1 <= 0.5f)

+    {

+        return rng.uniform() < deathProb(mDomain.size()) ? death(&rng) : birth(&rng);

+    }

+    return (u1 < 0.75f) ? move(&rng) : exchange(&rng);

 }

 // add an atom at pos, calculate mass either with an exponential distribution

 // or with the gibbs mass calculation

 template <class T, class MatA, class MatB>

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

+void GibbsSampler<T, MatA, MatB>::birth(GapsRng *rng)

 {

-    unsigned row = impl()->getRow(prop->birthPos);

-    unsigned col = impl()->getCol(prop->birthPos);

+    uint64_t pos = mDomain.randomFreePosition();

+    AtomicProposal prop = AtomicProposal('B', pos);

+

+    unsigned row = impl()->getRow(prop.birthPos);

+    unsigned col = impl()->getCol(prop.birthPos);

     // calculate proposed mass

     float mass = 0.f;

     if (impl()->canUseGibbs(row, col))

     {

         AlphaParameters alpha = impl()->alphaParameters(row, col);

-        mass = gibbsMass(alpha, &(prop->rng)).value; // 0 if it fails

+        mass = gibbsMass(alpha, &(prop.rng)).value; // 0 if it fails

     }

     else

     {

-        mass = prop->rng.exponential(mLambda);

+        mass = prop.rng.exponential(mLambda);

     }

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

     if (mass >= gaps::epsilon)

     {

-        mDomain.cacheInsert(prop->birthPos, mass);

+        mDomain.insert(prop.birthPos, mass);

         mMatrix(row, col) += mass;

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

-        mQueue.acceptBirth();

-    }

-    else

-    {

-        mQueue.rejectBirth();

     }

 }

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

 // the original mass or the gibbs mass calculation

 template <class T, class MatA, class MatB>

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

+void GibbsSampler<T, MatA, MatB>::death(GapsRng *rng)

 {

+    Atom* a = mDomain.randomAtom();

+    AtomicProposal prop = AtomicProposal('D', a);

+

     // calculate bin for this atom

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

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

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

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

     // kill off atom

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

+    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;

+    float rebirthMass = prop.atom1->mass;

     AlphaParameters alpha = impl()->alphaParameters(row, col);

     if (impl()->canUseGibbs(row, col))

     {

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

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

         if (gMass.hasValue)

         {

             rebirthMass = gMass.value;

@@ -394,43 +391,48 @@ void GibbsSampler<T, MatA, MatB>::death(AtomicProposal *prop)

     // accept/reject rebirth

     float deltaLL = rebirthMass * (alpha.su - alpha.s * rebirthMass / 2.f);

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

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

     {

-        prop->atom1->mass = rebirthMass;

+        prop.atom1->mass = rebirthMass;

         mMatrix(row, col) += rebirthMass;

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

-        mQueue.rejectDeath();

     }

     else

     {

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

-        mQueue.acceptDeath();

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

     }

 }

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

 template <class T, class MatA, class MatB>

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

+void GibbsSampler<T, MatA, MatB>::move(GapsRng *rng)

 {

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

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

-    unsigned r2 = impl()->getRow(prop->moveDest);

-    unsigned c2 = impl()->getCol(prop->moveDest);

-

-    GAPS_ASSERT(r1 != r2 || c1 != c2); // same bin handled during proposal

-

-    float deltaLL = impl()->computeDeltaLL(r1, c1, -1.f * prop->atom1->mass,

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

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

+    AtomNeighborhood hood = mDomain.randomAtomWithNeighbors();

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

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

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

+    AtomicProposal prop = AtomicProposal('M', hood.center, newLocation);

+

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

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

+    unsigned r2 = impl()->getRow(prop.moveDest);

+    unsigned c2 = impl()->getCol(prop.moveDest);

+

+    if (r1 != r2 || c1 != c2)

     {

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

+        float deltaLL = impl()->computeDeltaLL(r1, c1, -1.f * prop.atom1->mass,

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

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

+        {

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

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

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

-        mMatrix(r1, c1) = newVal;

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

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

+            mMatrix(r1, c1) = newVal;

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

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

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

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

+        }

     }

 }

@@ -438,53 +440,57 @@ void GibbsSampler<T, MatA, MatB>::move(AtomicProposal *prop)

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

 // the exchange

 template <class T, class MatA, class MatB>

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

+void GibbsSampler<T, MatA, MatB>::exchange(GapsRng *rng)

 {

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

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

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

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

-

-    GAPS_ASSERT(r1 != r2 || c1 != c2); // same bin handled during proposal

+    AtomNeighborhood hood = mDomain.randomAtomWithRightNeighbor();

+    Atom* a1 = hood.center;

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

+    AtomicProposal prop = AtomicProposal('E', a1, a2);

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

-    float m2 = prop->atom2->mass;

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

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

+    unsigned r2 = impl()->getRow(prop.atom2->pos);

+    unsigned c2 = impl()->getCol(prop.atom2->pos);

-    if (impl()->canUseGibbs(r1, c1, r2, c2))

+    if (r1 != r2 || c1 != c2)

     {

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

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

-        if (gMass.hasValue)

+        float m1 = prop.atom1->mass;

+        float m2 = prop.atom2->mass;

+

+        if (impl()->canUseGibbs(r1, c1, r2, c2))

         {

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

-            return;

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

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

+            if (gMass.hasValue)

+            {

+                acceptExchange(&prop, gMass.value, r1, c1, r2, c2);

+                return;

+            }

         }

-    }

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

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

-    float delta = m1 > m2 ? newMass - m1 : m2 - newMass; // change larger mass

-    float pOldMass = 2.f * newMass > m1 + m2 ? gaps::max(m1, m2) : gaps::min(m1, m2);

+        float delta = m1 > m2 ? newMass - m1 : m2 - newMass; // change larger mass

+        float pOldMass = 2.f * newMass > m1 + m2 ? gaps::max(m1, m2) : gaps::min(m1, m2);

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

-    float pOld = gaps::d_gamma(pOldMass, 2.f, 1.f / mLambda);

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

+        float pOld = gaps::d_gamma(pOldMass, 2.f, 1.f / mLambda);

-    if (pOld == 0.f && pNew != 0.f) // special case

-    {

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

-        return;

-    }

+        if (pOld == 0.f && pNew != 0.f) // special case

+        {

+            acceptExchange(&prop, delta, r1, c1, r2, c2);

+            return;

+        }

-    float deltaLL = impl()->computeDeltaLL(r1, c1, delta, r2, c2, -delta);

-    float priorLL = (pOld == 0.f) ? 1.f : pOld / pNew;

-    float u = std::log(prop->rng.uniform() * priorLL);

-    if (u < deltaLL * mAnnealingTemp)

-    {

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

-        return;

+        float deltaLL = impl()->computeDeltaLL(r1, c1, delta, r2, c2, -delta);

+        float priorLL = (pOld == 0.f) ? 1.f : pOld / pNew;

+        float u = std::log(prop.rng.uniform() * priorLL);

+        if (u < deltaLL * mAnnealingTemp)

+        {

+            acceptExchange(&prop, delta, r1, c1, r2, c2);

+            return;

+        }

     }

-

-    mQueue.rejectDeath();

 }

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

@@ -495,8 +501,7 @@ bool GibbsSampler<T, MatA, MatB>::updateAtomMass(Atom *atom, float delta)

     if (atom->mass + delta < gaps::epsilon)

     {

         DEBUG_PING // want to know if this ever happens

-        mDomain.cacheErase(atom->pos);

-        mQueue.acceptDeath();

+        mDomain.erase(atom->pos);

         return false;

     }

     atom->mass += delta;

@@ -518,11 +523,6 @@ float d1, unsigned r1, unsigned c1, unsigned r2, unsigned c2)

     if (!b1) { d1 = -1.f * prop->atom1->mass; }

     if (!b2) { d2 = -1.f * prop->atom2->mass; }

-    if (b1 && b2)

-    {

-        mQueue.rejectDeath();

-    }

-

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

@@ -601,7 +601,7 @@ bool GibbsSampler<T, MatA, MatB>::internallyConsistent()

 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

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

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

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

     return ar;

@@ -610,7 +610,7 @@ Archive& operator<<(Archive &ar, GibbsSampler<T, MatA, MatB> &s)

 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

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

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

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

     return ar;

@@ -4,6 +4,7 @@ PKG_LIBS = @GAPS_LIBS@

 OBJECTS =   AtomicDomain.o \

             Cogaps.o \

+            GapsResult.o \

             GapsRunner.o \

             GapsStatistics.o \

             GibbsSampler.o \

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

-#include "GapsAssert.h"

+#include "utils/GapsAssert.h"

 #include "ProposalQueue.h"

 #include "math/Random.h"

@@ -1,7 +1,7 @@

 #ifndef __GAPS_PROPOSAL_QUEUE_H__

 #define __GAPS_PROPOSAL_QUEUE_H__

-#include "Archive.h"

+#include "utils/Archive.h"

 #include "AtomicDomain.h"

 #include "data_structures/EfficientSets.h"

 #include "math/Random.h"

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

 #include "catch.h"

 #include "../AtomicDomain.h"

-#include "../GapsPrint.h"

+#include "../utils/GapsPrint.h"

 TEST_CASE("AtomicDomain")

 {

@@ -1,12 +1,12 @@

 #include "catch.h"

-#include "../Archive.h"

+#include "../utils/Archive.h"

 #include "../data_structures/Matrix.h"

 #include "../GibbsSampler.h"

 #include "../math/Random.h"

 #if 0

-TEST_CASE("Test Archive.h")

+TEST_CASE("Test utils/Archive.h")

 {

     SECTION("Reading/Writing to an Archive")

     {

@@ -1,7 +1,7 @@

 #include "Matrix.h"

 #include "../file_parser/CsvParser.h"

 #include "../file_parser/MatrixElement.h"

-#include "../GapsAssert.h"

+#include "../utils/GapsAssert.h"

 template <class MatA, class MatB>

 inline void copyMatrix(MatA &dest, const MatB &source)

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

 #ifndef __COGAPS_MATRIX_H__

 #define __COGAPS_MATRIX_H__

-#include "../Archive.h"

-#include "../GapsAssert.h"

+#include "../utils/Archive.h"

+#include "../utils/GapsAssert.h"

 #include "../file_parser/FileParser.h"

 #include "../math/Math.h"

 #include "Vector.h"

@@ -1,7 +1,7 @@

 #ifndef __COGAPS_VECTOR_H__

 #define __COGAPS_VECTOR_H__

-#include "../Archive.h"

+#include "../utils/Archive.h"

 #include <boost/align/aligned_allocator.hpp>

 #include <vector>

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

 #include "CsvParser.h"

-#include "../GapsAssert.h"

+#include "../utils/GapsAssert.h"

 #include <fstream>

 #include <iostream>

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

-#include "../GapsAssert.h"

+#include "../utils/GapsAssert.h"

 #include "CsvParser.h"

 #include "FileParser.h"

 #include "MtxParser.h"

@@ -1,7 +1,7 @@

 #ifndef __COGAPS_FILE_PARSER_H__

 #define __COGAPS_FILE_PARSER_H__

-#include "../GapsAssert.h"

+#include "../utils/GapsAssert.h"

 #include "MatrixElement.h"

 #include <fstream>

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

 #include "GctParser.h"

-#include "../GapsAssert.h"

+#include "../utils/GapsAssert.h"

 #include "../math/Algorithms.h"

 #include <fstream>

@@ -1,7 +1,7 @@

 #include "MatrixElement.h"

 #include "MtxParser.h"

-#include "../GapsAssert.h"

+#include "../utils/GapsAssert.h"