@@ -45,7 +45,7 @@ CoGAPS <- function(D, S, nFactor=7, nEquil=1000, nSample=1000, nOutputs=1000,

 nSnapshots=0, alphaA=0.01, alphaP=0.01, maxGibbmassA=100, maxGibbmassP=100,

 seed=-1, messages=TRUE, singleCellRNASeq=FALSE, whichMatrixFixed='N',

 fixedPatterns=matrix(0), checkpointInterval=0, 

-checkpointFile="gaps_checkpoint.out", ...)

+checkpointFile="gaps_checkpoint.out", nCores=1, ...)

 {

     # get v2 arguments

     oldArgs <- list(...)

@@ -100,7 +100,8 @@ checkpointFile="gaps_checkpoint.out", ...)

     result <- cogaps_cpp(D, S, nFactor, nEquil, nEquil/10, nSample, nOutputs,

         nSnapshots, alphaA, alphaP, maxGibbmassA, maxGibbmassP, seed, messages,

         singleCellRNASeq, whichMatrixFixed, fixedPatterns, checkpointInterval,

-        checkpointFile, which(thresholdEnum==pumpThreshold), nPumpSamples)

+        checkpointFile, which(thresholdEnum==pumpThreshold), nPumpSamples,

+        nCores)

     # label matrices and return list

     patternNames <- paste('Patt', 1:nFactor, sep='')

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

 # Generated by using Rcpp::compileAttributes() -> do not edit by hand

 # Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393

-cogaps_cpp <- function(D, S, nFactor, nEquil, nEquilCool, nSample, nOutputs, nSnapshots, alphaA, alphaP, maxGibbmassA, maxGibbmassP, seed, messages, singleCellRNASeq, whichMatrixFixed, FP, checkpointInterval, cptFile, pumpThreshold, nPumpSamples) {

-    .Call('_CoGAPS_cogaps_cpp', PACKAGE = 'CoGAPS', D, S, nFactor, nEquil, nEquilCool, nSample, nOutputs, nSnapshots, alphaA, alphaP, maxGibbmassA, maxGibbmassP, seed, messages, singleCellRNASeq, whichMatrixFixed, FP, checkpointInterval, cptFile, pumpThreshold, nPumpSamples)

+cogaps_cpp <- function(D, S, nFactor, nEquil, nEquilCool, nSample, nOutputs, nSnapshots, alphaA, alphaP, maxGibbmassA, maxGibbmassP, seed, messages, singleCellRNASeq, whichMatrixFixed, FP, checkpointInterval, cptFile, pumpThreshold, nPumpSamples, nCores) {

+    .Call('_CoGAPS_cogaps_cpp', PACKAGE = 'CoGAPS', D, S, nFactor, nEquil, nEquilCool, nSample, nOutputs, nSnapshots, alphaA, alphaP, maxGibbmassA, maxGibbmassP, seed, messages, singleCellRNASeq, whichMatrixFixed, FP, checkpointInterval, cptFile, pumpThreshold, nPumpSamples, nCores)

 }

 cogapsFromCheckpoint_cpp <- function(D, S, nFactor, nEquil, nSample, fileName, cptFile) {

@@ -9,7 +9,7 @@ CoGAPS(D, S, nFactor = 7, nEquil = 1000, nSample = 1000,

   maxGibbmassA = 100, maxGibbmassP = 100, seed = -1, messages = TRUE,

   singleCellRNASeq = FALSE, whichMatrixFixed = "N",

   fixedPatterns = matrix(0), checkpointInterval = 0,

-  checkpointFile = "gaps_checkpoint.out", ...)

+  checkpointFile = "gaps_checkpoint.out", nCores = 1, ...)

 }

 \arguments{

 \item{D}{data matrix}

@@ -2,9 +2,17 @@

 #include "AtomicDomain.h"

 #include "math/Random.h"

+#include <omp.h>

 #include <stdint.h>

 #include <utility>

+//omp_lock_t lock;

+

+AtomicDomain::AtomicDomain()

+{

+    //omp_init_lock(&lock);

+}

+

 // O(1)

 Atom AtomicDomain::front() const

 {

@@ -76,6 +84,8 @@ bool AtomicDomain::hasRight(const Atom &atom) const

 // O(logN)

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

 {

+    //omp_set_lock(&lock);

+

     // insert position into map

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

     iter = mAtomPositions.insert(std::pair<uint64_t, uint64_t>(pos, size())).first;

@@ -100,6 +110,9 @@ Atom AtomicDomain::insert(uint64_t pos, float mass)

     // add atom to vector

     mAtoms.push_back(atom);

     mUsedPositions.insert(pos);

+

+    //omp_unset_lock(&lock);

+

     return atom;

 }

@@ -108,6 +121,8 @@ Atom AtomicDomain::insert(uint64_t pos, float mass)

 // swap with last atom in vector, pop off the back

 void AtomicDomain::erase(uint64_t pos)

 {

+    //omp_set_lock(&lock);

+

     // get vector index of this atom and erase it

     uint64_t index = mAtomPositions.at(pos);

@@ -146,10 +161,14 @@ void AtomicDomain::erase(uint64_t pos)

     mAtomPositions.erase(pos);

     mAtoms.pop_back();

     mUsedPositions.erase(pos);

+

+    //omp_unset_lock(&lock);

 }

 // O(logN)

 void AtomicDomain::updateMass(uint64_t pos, float newMass)

 {

+    //omp_set_lock(&lock);

     mAtoms[mAtomPositions.at(pos)].mass = newMass;

+    //omp_unset_lock(&lock);

 }

\ No newline at end of file

@@ -60,6 +60,8 @@ private:

 public:

+    AtomicDomain();

+

     void setDomainSize(uint64_t size) { mDomainSize = size; }

     // access atoms

@@ -10,13 +10,14 @@ unsigned nEquilCool, unsigned nSample, unsigned nOutputs, unsigned nSnapshots,

 float alphaA, float alphaP, float maxGibbmassA, float maxGibbmassP,

 unsigned seed, bool messages, bool singleCellRNASeq, char whichMatrixFixed,

 const Rcpp::NumericMatrix &FP, unsigned checkpointInterval,

-const std::string &cptFile, unsigned pumpThreshold, unsigned nPumpSamples)

+const std::string &cptFile, unsigned pumpThreshold, unsigned nPumpSamples,

+unsigned nCores)

 {

     // create internal state from parameters and run from there

     GapsRunner runner(D, S, nFactor, nEquil, nEquilCool, nSample,

         nOutputs, nSnapshots, alphaA, alphaP, maxGibbmassA, maxGibbmassP, seed,

         messages, singleCellRNASeq, checkpointInterval, cptFile,

-        whichMatrixFixed, FP);

+        whichMatrixFixed, FP, nCores);

     return runner.run();

 }

@@ -5,7 +5,7 @@ unsigned nFactor, unsigned nEquil, unsigned nCool, unsigned nSample,

 unsigned nOutputs, unsigned nSnapshots, float alphaA, float alphaP,

 float maxGibbsMassA, float maxGibbsMassP, uint32_t seed, bool messages,

 bool singleCellRNASeq, unsigned cptInterval, const std::string &cptFile,

-char whichMatrixFixed, const Rcpp::NumericMatrix &FP)

+char whichMatrixFixed, const Rcpp::NumericMatrix &FP, unsigned nCores)

     :

 mChiSqEquil(nEquil), mNumAAtomsEquil(nEquil), mNumPAtomsEquil(nEquil),

 mChiSqSample(nSample), mNumAAtomsSample(nSample), mNumPAtomsSample(nSample),

@@ -16,7 +16,8 @@ mCheckpointInterval(cptInterval), mCheckpointFile(cptFile),

 mNumUpdatesA(0), mNumUpdatesP(0),

 mASampler(D, S, nFactor, alphaA, maxGibbsMassA),

 mPSampler(D, S, nFactor, alphaP, maxGibbsMassP),

-mStatistics(D.nrow(), D.ncol(), nFactor)

+mStatistics(D.nrow(), D.ncol(), nFactor),

+mNumCores(nCores)

 {

     mASampler.sync(mPSampler);

     mPSampler.sync(mASampler);

@@ -137,11 +138,11 @@ void GapsRunner::runSampPhase()

 void GapsRunner::updateSampler()

 {

     mNumUpdatesA += mIterA;

-    mASampler.update(mIterA);

+    mASampler.update(mIterA, mNumCores);

     mPSampler.sync(mASampler);

     mNumUpdatesP += mIterP;

-    mPSampler.update(mIterP);

+    mPSampler.update(mIterP, mNumCores);

     mASampler.sync(mPSampler);

 }

@@ -63,6 +63,8 @@ public:

     PatternGibbsSampler mPSampler;

     GapsStatistics mStatistics;

+    unsigned mNumCores;

+

     void createCheckpoint();

     void makeCheckpointIfNeeded();

     void displayStatus(const std::string &type, unsigned nIterTotal);

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

         unsigned nOutputs, unsigned nSnapshots, float alphaA, float alphaP,

         float maxGibbsMassA, float maxGibbsMassP, uint32_t seed, bool messages,

         bool singleCellRNASeq, unsigned cptInterval, const std::string &cptFile,

-        char whichMatrixFixed, const Rcpp::NumericMatrix &FP);

+        char whichMatrixFixed, const Rcpp::NumericMatrix &FP, unsigned nCores);

     // construct from checkpoint file

     GapsRunner(const Rcpp::NumericMatrix &D, const Rcpp::NumericMatrix &S,

@@ -77,7 +77,7 @@ protected:

 public:

-    void update(unsigned nSteps);

+    void update(unsigned nSteps, unsigned nCores);

     void setAnnealingTemp(float temp);

     float getAvgQueue() const { return mAvgQueue; }

@@ -171,19 +171,23 @@ T* GibbsSampler<T, MatA, MatB>::impl()

 }

 template <class T, class MatA, class MatB>

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

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

 {

+    static unsigned count = 0;

     unsigned n = 0;

     while (n < nSteps)

     {

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

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

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

+        

+        //mQueue.populate(mDomain, std::min(2u, nSteps - (mQueue.size() + n)));

         GAPS_ASSERT(mQueue.size() > 0);

         mNumQueues += 1.f;

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

         n += mQueue.size();

-        //#pragma omp parallel for

+        //Rprintf("round: %d\n", count++);

+

+        #pragma omp parallel for num_threads(nCores)

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

         {

             processProposal(mQueue[i]);

@@ -199,6 +203,7 @@ void GibbsSampler<T, MatA, MatB>::processProposal(const AtomicProposal &prop)

     unsigned r1 = impl()->getRow(prop.pos1);

     unsigned c1 = impl()->getCol(prop.pos1);

     unsigned r2 = 0, c2 = 0;

+    //Rprintf("type: %c pos1: %lu\n", prop.type, prop.pos1);

     switch (prop.type)

     {

         case 'B':

@@ -213,6 +218,7 @@ void GibbsSampler<T, MatA, MatB>::processProposal(const AtomicProposal &prop)

             move(prop.pos1, prop.mass1, prop.pos2, r1, c1, r2, c2);

             break;

         case 'E':

+            //Rprintf("pos2: %lu\n", prop.pos2);

             r2 = impl()->getRow(prop.pos2);

             c2 = impl()->getCol(prop.pos2);

             exchange(prop.pos1, prop.mass1, prop.pos2, prop.mass2, r1, c1, r2, c2);

@@ -223,17 +229,23 @@ void GibbsSampler<T, MatA, MatB>::processProposal(const AtomicProposal &prop)

 template <class T, class MatA, class MatB>

 void GibbsSampler<T, MatA, MatB>::addMass(uint64_t pos, float mass, unsigned row, unsigned col)

 {

-    mDomain.insert(pos, mass);

-    mMatrix(row, col) += mass;

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

+    #pragma omp critical(gibbs)

+    {

+        mDomain.insert(pos, mass);

+        mMatrix(row, col) += mass;

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

+    }

 }

 template <class T, class MatA, class MatB>

 void GibbsSampler<T, MatA, MatB>::removeMass(uint64_t pos, float mass, unsigned row, unsigned col)

 {

-    mDomain.erase(pos);

-    mMatrix(row, col) += -mass;

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

+    #pragma omp critical(gibbs)

+    {

+        mDomain.erase(pos);

+        mMatrix(row, col) += -mass;

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

+    }

 }

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

@@ -244,16 +256,19 @@ unsigned col)

 {

     float mass = impl()->canUseGibbs(row, col) ? gibbsMass(row, col, mass)

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

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

+    #pragma omp critical(gibbs)

     {

-        mDomain.updateMass(pos, mass);

-        mMatrix(row, col) += mass;

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

-    }

-    else

-    {

-        mDomain.erase(pos);

-        mQueue.rejectBirth();

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

+        {

+            mDomain.updateMass(pos, mass);

+            mMatrix(row, col) += mass;

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

+        }

+        else

+        {

+            mDomain.erase(pos);

+            mQueue.rejectBirth();

+        }

     }

 }

@@ -263,6 +278,7 @@ template <class T, class MatA, class MatB>

 void GibbsSampler<T, MatA, MatB>::death(uint64_t pos, float mass, unsigned row,

 unsigned col)

 {

+    GAPS_ASSERT(mass > 0.f);

     removeMass(pos, mass, row, col);

     float newMass = impl()->canUseGibbs(row, col) ? gibbsMass(row, col, -mass) : 0.f;

     mass = newMass < gaps::algo::epsilon ? mass : newMass;

@@ -283,6 +299,7 @@ template <class T, class MatA, class MatB>

 void GibbsSampler<T, MatA, MatB>::move(uint64_t src, float mass, uint64_t dest,

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

 {

+    GAPS_ASSERT(mass > 0.f);

     if (r1 != r2 || c1 != c2) // automatically reject if change in same bin

     {

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

@@ -301,6 +318,8 @@ template <class T, class MatA, class MatB>

 void GibbsSampler<T, MatA, MatB>::exchange(uint64_t p1, float m1, uint64_t p2,

 float m2, unsigned r1, unsigned c1, unsigned r2, unsigned c2)

 {

+    GAPS_ASSERT(m1 > 0.f);

+    GAPS_ASSERT(m2 > 0.f);

     float pUpper = gaps::random::p_gamma(m1 + m2, 2.f, 1.f / mLambda);

     float newMass = gaps::random::inverseGammaSample(0.f, pUpper, 2.f, 1.f / mLambda);

     if (r1 != r2 || c1 != c2) // automatically reject if change in same bin

@@ -355,17 +374,22 @@ template <class T, class MatA, class MatB>

 float GibbsSampler<T, MatA, MatB>::updateAtomMass(uint64_t pos, float mass,

 float delta)

 {

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

-    {

-        mDomain.erase(pos);

-        mQueue.acceptDeath();

-        return false;

-    }

-    else

+    bool ret_val = false;

+    #pragma omp critical(gibbs)

     {

-        mDomain.updateMass(pos, mass + delta);

-        return true;

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

+        {

+            mDomain.erase(pos);

+            mQueue.acceptDeath();

+            ret_val = false;

+        }

+        else

+        {

+            mDomain.updateMass(pos, mass + delta);

+            ret_val = true;

+        }

     }

+    return ret_val;

 }

 // helper function for exchange step, updates the atomic domain, matrix, and

@@ -375,6 +399,7 @@ 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)

 {

+    //Rprintf("p1: %lu p2: %lu\n", p1, p2);

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

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

     GAPS_ASSERT(b1 || b2);

@@ -387,10 +412,13 @@ unsigned r2, unsigned c2)

         mQueue.rejectDeath();

     }

-    mMatrix(r1, c1) += d1;

-    mMatrix(r2, c2) += d2;

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

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

+    #pragma omp critical(gibbs)

+    {

+        mMatrix(r1, c1) += d1;

+        mMatrix(r2, c2) += d2;

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

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

+    }

 }

 template <class T, class MatA, class MatB>

@@ -52,17 +52,22 @@ const AtomicProposal& ProposalQueue::operator[](int n) const

 void ProposalQueue::acceptDeath()

 {

+    #pragma omp atomic

     mMaxAtoms--;

 }

 void ProposalQueue::rejectDeath()

 {

+    #pragma omp atomic

     mMinAtoms++;

 }

 void ProposalQueue::rejectBirth()

 {

+    #pragma omp atomic

     mMinAtoms--;

+

+    #pragma omp atomic

     mMaxAtoms--;

 }

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

 using namespace Rcpp;

 // cogaps_cpp

-Rcpp::List cogaps_cpp(const Rcpp::NumericMatrix& D, const Rcpp::NumericMatrix& S, unsigned nFactor, unsigned nEquil, unsigned nEquilCool, unsigned nSample, unsigned nOutputs, unsigned nSnapshots, float alphaA, float alphaP, float maxGibbmassA, float maxGibbmassP, unsigned seed, bool messages, bool singleCellRNASeq, char whichMatrixFixed, const Rcpp::NumericMatrix& FP, unsigned checkpointInterval, const std::string& cptFile, unsigned pumpThreshold, unsigned nPumpSamples);

-RcppExport SEXP _CoGAPS_cogaps_cpp(SEXP DSEXP, SEXP SSEXP, SEXP nFactorSEXP, SEXP nEquilSEXP, SEXP nEquilCoolSEXP, SEXP nSampleSEXP, SEXP nOutputsSEXP, SEXP nSnapshotsSEXP, SEXP alphaASEXP, SEXP alphaPSEXP, SEXP maxGibbmassASEXP, SEXP maxGibbmassPSEXP, SEXP seedSEXP, SEXP messagesSEXP, SEXP singleCellRNASeqSEXP, SEXP whichMatrixFixedSEXP, SEXP FPSEXP, SEXP checkpointIntervalSEXP, SEXP cptFileSEXP, SEXP pumpThresholdSEXP, SEXP nPumpSamplesSEXP) {

+Rcpp::List cogaps_cpp(const Rcpp::NumericMatrix& D, const Rcpp::NumericMatrix& S, unsigned nFactor, unsigned nEquil, unsigned nEquilCool, unsigned nSample, unsigned nOutputs, unsigned nSnapshots, float alphaA, float alphaP, float maxGibbmassA, float maxGibbmassP, unsigned seed, bool messages, bool singleCellRNASeq, char whichMatrixFixed, const Rcpp::NumericMatrix& FP, unsigned checkpointInterval, const std::string& cptFile, unsigned pumpThreshold, unsigned nPumpSamples, unsigned nCores);

+RcppExport SEXP _CoGAPS_cogaps_cpp(SEXP DSEXP, SEXP SSEXP, SEXP nFactorSEXP, SEXP nEquilSEXP, SEXP nEquilCoolSEXP, SEXP nSampleSEXP, SEXP nOutputsSEXP, SEXP nSnapshotsSEXP, SEXP alphaASEXP, SEXP alphaPSEXP, SEXP maxGibbmassASEXP, SEXP maxGibbmassPSEXP, SEXP seedSEXP, SEXP messagesSEXP, SEXP singleCellRNASeqSEXP, SEXP whichMatrixFixedSEXP, SEXP FPSEXP, SEXP checkpointIntervalSEXP, SEXP cptFileSEXP, SEXP pumpThresholdSEXP, SEXP nPumpSamplesSEXP, SEXP nCoresSEXP) {

 BEGIN_RCPP

     Rcpp::RObject rcpp_result_gen;

     Rcpp::RNGScope rcpp_rngScope_gen;

@@ -32,7 +32,8 @@ BEGIN_RCPP

     Rcpp::traits::input_parameter< const std::string& >::type cptFile(cptFileSEXP);

     Rcpp::traits::input_parameter< unsigned >::type pumpThreshold(pumpThresholdSEXP);

     Rcpp::traits::input_parameter< unsigned >::type nPumpSamples(nPumpSamplesSEXP);

-    rcpp_result_gen = Rcpp::wrap(cogaps_cpp(D, S, nFactor, nEquil, nEquilCool, nSample, nOutputs, nSnapshots, alphaA, alphaP, maxGibbmassA, maxGibbmassP, seed, messages, singleCellRNASeq, whichMatrixFixed, FP, checkpointInterval, cptFile, pumpThreshold, nPumpSamples));

+    Rcpp::traits::input_parameter< unsigned >::type nCores(nCoresSEXP);

+    rcpp_result_gen = Rcpp::wrap(cogaps_cpp(D, S, nFactor, nEquil, nEquilCool, nSample, nOutputs, nSnapshots, alphaA, alphaP, maxGibbmassA, maxGibbmassP, seed, messages, singleCellRNASeq, whichMatrixFixed, FP, checkpointInterval, cptFile, pumpThreshold, nPumpSamples, nCores));

     return rcpp_result_gen;

 END_RCPP

 }

@@ -75,7 +76,7 @@ END_RCPP

 }

 static const R_CallMethodDef CallEntries[] = {

-    {"_CoGAPS_cogaps_cpp", (DL_FUNC) &_CoGAPS_cogaps_cpp, 21},

+    {"_CoGAPS_cogaps_cpp", (DL_FUNC) &_CoGAPS_cogaps_cpp, 22},

     {"_CoGAPS_cogapsFromCheckpoint_cpp", (DL_FUNC) &_CoGAPS_cogapsFromCheckpoint_cpp, 7},

     {"_CoGAPS_getBuildReport_cpp", (DL_FUNC) &_CoGAPS_getBuildReport_cpp, 0},

     {"_CoGAPS_run_catch_unit_tests", (DL_FUNC) &_CoGAPS_run_catch_unit_tests, 0},

@@ -10,7 +10,7 @@ TEST_CASE("Test Top Level CoGAPS Call")

     GapsRunner runner(D, S, 3, 500, 100, 500, 250, 0, 0.01f, 0.01f, 100.f,

         100.f, 123, false, false, 0, "gaps_checkpoint.out", 'N',

-        Rcpp::NumericMatrix(1,1));

+        Rcpp::NumericMatrix(1,1), 1);

     REQUIRE_NOTHROW(runner.run());

 }