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

 # necessary yml

 language: r

+cache: packages

 sudo: required

 # use bioconductor

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

 # Generated by roxygen2: do not edit by hand

 export(CoGAPS)

+<<<<<<< HEAD

 export(CoGAPSFromFile)

+=======

+>>>>>>> develop

 export(CoGapsFromCheckpoint)

 export(GWCoGAPS)

 export(GWCoGapsFromCheckpoint)

@@ -20,7 +20,6 @@

 #' @param nEquil number of iterations for burn-in

 #' @param nSample number of iterations for sampling

 #' @param nOutputs how often to print status into R by iterations

-#' @param nSnapshots the number of individual samples to capture

 #' @param alphaA sparsity parameter for A domain

 #' @param alphaP sparsity parameter for P domain

 #' @param maxGibbmassA limit truncated normal to max size

@@ -116,21 +115,20 @@ checkpointFile="gaps_checkpoint.out", nCores=1, ...)

 #'  continues the run from that point

 #' @param D data matrix

 #' @param S uncertainty matrix

+#' @param nFactor number of patterns

+#' @param nIter number of iterations

+#' @param checkpointFile path to checkpoint file

 #' @param path path to checkpoint file

-#' @param checkpointFile name for future checkpooints made

 #' @return list with A and P matrix estimates

 #' @examples

 #' data(SimpSim)

 #' result <- CoGAPS(SimpSim.D, SimpSim.S, nFactor=3, nOutputs=250)

-CoGapsFromCheckpoint <- function(D, S, path, checkpointFile=NA)

+CoGapsFromCheckpoint <- function(D, S, nFactor, nIter, checkpointFile)

 {

-    if (is.na(checkpointFile))

-        checkpointFile <- path

-    cogapsFromCheckpoint_cpp(D, S, path, checkpointFile)

+    cogapsFromCheckpoint_cpp(D, S, nFactor, nIter, nIter, checkpointFile)

 }

 #' CoGAPS with file input for matrix

-#' @export

 #'

 #' @param D file path for data matrix

 #' @return list with A and P matrix estimates

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

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

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

+<<<<<<< HEAD

 cogapsFromFile_cpp <- function(D, nPatterns, maxIter, outputFrequency, seed, alphaA, alphaP, maxGibbsMassA, maxGibbsMassP, messages, singleCellRNASeq) {

     .Call('_CoGAPS_cogapsFromFile_cpp', PACKAGE = 'CoGAPS', D, nPatterns, maxIter, outputFrequency, seed, alphaA, alphaP, maxGibbsMassA, maxGibbsMassP, messages, singleCellRNASeq)

 }

 cogaps_cpp <- function(D, nPatterns, maxIter, outputFrequency, seed, alphaA, alphaP, maxGibbsMassA, maxGibbsMassP, messages, singleCellRNASeq, nCores) {

     .Call('_CoGAPS_cogaps_cpp', PACKAGE = 'CoGAPS', D, nPatterns, maxIter, outputFrequency, seed, alphaA, alphaP, maxGibbsMassA, maxGibbsMassP, messages, singleCellRNASeq, nCores)

+=======

+cogaps_cpp <- function(D, S, nFactor, nEquil, nEquilCool, nSample, nOutputs, 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, alphaA, alphaP, maxGibbmassA, maxGibbmassP, seed, messages, singleCellRNASeq, whichMatrixFixed, FP, checkpointInterval, cptFile, pumpThreshold, nPumpSamples, nCores)

+}

+

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

+    .Call('_CoGAPS_cogapsFromCheckpoint_cpp', PACKAGE = 'CoGAPS', D, S, nFactor, nEquil, nSample, cptFile)

+}

+

+cogapsFromFile_cpp <- function(D) {

+    invisible(.Call('_CoGAPS_cogapsFromFile_cpp', PACKAGE = 'CoGAPS', D))

+>>>>>>> develop

 }

 getBuildReport_cpp <- function() {

@@ -639,6 +639,7 @@ ac_subst_files=''

 ac_user_opts='

 enable_option_checking

 enable_debug

+enable_debug_tsan

 enable_warnings

 enable_simd

 '

@@ -1273,6 +1274,7 @@ Optional Features:

   --disable-FEATURE       do not include FEATURE (same as --enable-FEATURE=no)

   --enable-FEATURE[=ARG]  include FEATURE [ARG=yes]

   --enable-debug          build debug version of CoGAPS

+  --enable-debug-tsan     build debug version of CoGAPS

   --enable-warnings       compile CoGAPS with warning messages

   --enable-simd           specify simd instruction set (sse, avx)

@@ -2984,6 +2986,15 @@ else

 fi

+# Check if compiling debug version with thread sanitizer

+# Check whether --enable-debug_tsan was given.

+if test "${enable_debug_tsan+set}" = set; then :

+  enableval=$enable_debug_tsan; build_debug_tsan=yes

+else

+  build_debug_tsan=no

+fi

+

+

 # Check if compiler warnings should be turned on

 # Check whether --enable-warnings was given.

 if test "${enable_warnings+set}" = set; then :

@@ -4177,7 +4188,21 @@ echo "building on $ax_cv_cxx_compiler_vendor compiler version $ax_cv_cxx_compile

 # set compile flags for debug build

 if test "x$build_debug" = "xyes" ; then

     echo "Building Debug Version of CoGAPS"

-    GAPS_CPP_FLAGS+=" -DGAPS_DEBUG "

+    GAPS_CPP_FLAGS+=" -DGAPS_DEBUG -fno-omit-frame-pointer "

+    if test "x$ax_cv_cxx_compiler_vendor" = "xclang" ; then

+        GAPS_CPP_FLAGS+=" -fsanitize=address -fsanitize=undefined "

+        GAPS_LIBS+=" -fsanitize=address -fsanitize=undefined "

+    elif test "x$ax_cv_cxx_compiler_vendor" = "xgnu" ; then

+        GAPS_CPP_FLAGS+=" -fsanitize=address -fsanitize=undefined "

+        GAPS_LIBS+=" -fsanitize=address -fsanitize=undefined "

+    fi

+elif test "x$build_debug_tsan" = "xyes" ; then

+    echo "Building Debug Version of CoGAPS with thread sanitizer"

+    GAPS_CPP_FLAGS+=" -DGAPS_DEBUG -fno-omit-frame-pointer "

+    if test "x$ax_cv_cxx_compiler_vendor" = "xclang" ; then

+        GAPS_CPP_FLAGS+=" -fsanitize=thread "

+        GAPS_LIBS+=" -fsanitize=thread "

+    fi

 fi

 # set compile flags if warnings enabled

@@ -19,6 +19,10 @@ AC_PROG_CXX

 AC_ARG_ENABLE(debug, [AC_HELP_STRING([--enable-debug],

     [build debug version of CoGAPS])], [build_debug=yes], [build_debug=no])

+# Check if compiling debug version with thread sanitizer

+AC_ARG_ENABLE(debug_tsan, [AC_HELP_STRING([--enable-debug-tsan],

+    [build debug version of CoGAPS])], [build_debug_tsan=yes], [build_debug_tsan=no])

+

 # Check if compiler warnings should be turned on

 AC_ARG_ENABLE(warnings, [AC_HELP_STRING([--enable-warnings],

     [compile CoGAPS with warning messages])], [warnings=yes], [warnings=no])

@@ -47,7 +51,21 @@ echo "building on $ax_cv_cxx_compiler_vendor compiler version $ax_cv_cxx_compile

 # set compile flags for debug build

 if test "x$build_debug" = "xyes" ; then

     echo "Building Debug Version of CoGAPS"

-    GAPS_CPP_FLAGS+=" -DGAPS_DEBUG "

+    GAPS_CPP_FLAGS+=" -DGAPS_DEBUG -fno-omit-frame-pointer "

+    if test "x$ax_cv_cxx_compiler_vendor" = "xclang" ; then

+        GAPS_CPP_FLAGS+=" -fsanitize=address -fsanitize=undefined "

+        GAPS_LIBS+=" -fsanitize=address -fsanitize=undefined "

+    elif test "x$ax_cv_cxx_compiler_vendor" = "xgnu" ; then

+        GAPS_CPP_FLAGS+=" -fsanitize=address -fsanitize=undefined "

+        GAPS_LIBS+=" -fsanitize=address -fsanitize=undefined "

+    fi

+elif test "x$build_debug_tsan" = "xyes" ; then

+    echo "Building Debug Version of CoGAPS with thread sanitizer"

+    GAPS_CPP_FLAGS+=" -DGAPS_DEBUG -fno-omit-frame-pointer "

+    if test "x$ax_cv_cxx_compiler_vendor" = "xclang" ; then

+        GAPS_CPP_FLAGS+=" -fsanitize=thread "

+        GAPS_LIBS+=" -fsanitize=thread "

+    fi

 fi

 # set compile flags if warnings enabled

@@ -5,11 +5,19 @@

 \title{CoGAPS Matrix Factorization Algorithm}

 \usage{

 CoGAPS(D, S, nFactor = 7, nEquil = 250, nSample = 250, nOutputs = 1000,

+<<<<<<< HEAD

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

   maxGibbmassP = 100, seed = NA, messages = TRUE,

   singleCellRNASeq = FALSE, whichMatrixFixed = "N",

   fixedPatterns = matrix(0), checkpointInterval = 0,

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

+=======

+  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",

+  nCores = 1, ...)

+>>>>>>> develop

 }

 \arguments{

 \item{D}{data matrix}

@@ -25,8 +33,6 @@ greater than or equal to the number of rows of FP}

 \item{nOutputs}{how often to print status into R by iterations}

-\item{nSnapshots}{the number of individual samples to capture}

-

 \item{alphaA}{sparsity parameter for A domain}

 \item{alphaP}{sparsity parameter for P domain}

@@ -4,16 +4,20 @@

 \alias{CoGapsFromCheckpoint}

 \title{Restart CoGAPS from Checkpoint File}

 \usage{

-CoGapsFromCheckpoint(D, S, path, checkpointFile = NA)

+CoGapsFromCheckpoint(D, S, nFactor, nIter, checkpointFile)

 }

 \arguments{

 \item{D}{data matrix}

 \item{S}{uncertainty matrix}

-\item{path}{path to checkpoint file}

+\item{nFactor}{number of patterns}

+

+\item{nIter}{number of iterations}

-\item{checkpointFile}{name for future checkpooints made}

+\item{checkpointFile}{path to checkpoint file}

+

+\item{path}{path to checkpoint file}

 }

 \value{

 list with A and P matrix estimates

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

 #include <fstream>

+#include <boost/random/mersenne_twister.hpp>

+

 // flags for opening an archive

 #define ARCHIVE_READ  std::ios::in

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

@@ -39,6 +41,48 @@ public:

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

+    template<typename T>

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

+    {

+        ar.mStream.write(reinterpret_cast<char*>(&val), sizeof(T)); // NOLINT

+        return ar;

+    }

+

+    template<typename T>

+    friend Archive& readFromArchive(Archive &ar, T &val)

+    {

+        ar.mStream.read(reinterpret_cast<char*>(&val), sizeof(T)); // NOLINT

+        return ar;

+    }    

+

+    // explicitly define which types can be automatically written/read

+    // don't have C++11 and don't want to add another dependency on boost,

+    // so no template tricks

+

+    friend Archive& operator<<(Archive &ar, bool val)     { return writeToArchive(ar, val); }

+    friend Archive& operator<<(Archive &ar, int val)      { return writeToArchive(ar, val); }

+    friend Archive& operator<<(Archive &ar, unsigned val) { return writeToArchive(ar, val); }

+    friend Archive& operator<<(Archive &ar, uint64_t val) { return writeToArchive(ar, val); }

+    friend Archive& operator<<(Archive &ar, int64_t val)  { return writeToArchive(ar, val); }

+    friend Archive& operator<<(Archive &ar, float val)    { return writeToArchive(ar, val); }

+    friend Archive& operator<<(Archive &ar, double val)   { return writeToArchive(ar, val); }

+    friend Archive& operator<<(Archive &ar, boost::random::mt11213b val) { return writeToArchive(ar, val); }

+

+    friend Archive& operator>>(Archive &ar, bool &val)     { return readFromArchive(ar, val); }

+    friend Archive& operator>>(Archive &ar, int &val)      { return readFromArchive(ar, val); }

+    friend Archive& operator>>(Archive &ar, unsigned &val) { return readFromArchive(ar, val); }

+    friend Archive& operator>>(Archive &ar, uint64_t &val) { return readFromArchive(ar, val); }

+    friend Archive& operator>>(Archive &ar, int64_t &val)  { return readFromArchive(ar, val); }

+    friend Archive& operator>>(Archive &ar, float &val)    { return readFromArchive(ar, val); }

+    friend Archive& operator>>(Archive &ar, double &val)   { return readFromArchive(ar, val); }

+    friend Archive& operator>>(Archive &ar, boost::random::mt11213b &val) { return readFromArchive(ar, val); }

+

+/*

+    friend Archive& operator>>(Archive &ar, uint64_t &val)

+    {

+        return readPrimitiveFromArchive(ar, val);

+    }

+

     template<typename T>

     friend Archive& operator<<(Archive &ar, T val)

     {

@@ -46,12 +90,14 @@ public:

         return ar;

     }

+

     template<typename T>

     friend Archive& operator>>(Archive &ar, T &val)

     {

         ar.mStream.read(reinterpret_cast<char*>(&val), sizeof(T)); // NOLINT

         return ar;

     }

+*/

 };

 #endif

@@ -43,8 +43,8 @@ public:

     }

     // serialization

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

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

+    friend Archive& operator<<(Archive &ar, Atom &a);

+    friend Archive& operator>>(Archive &ar, Atom &a);

 };

 struct RawAtom

@@ -85,7 +85,7 @@ bool singleCellRNASeq, unsigned nCores)

 /*

 Rcpp::List cogapsFromCheckpoint_cpp(const Rcpp::NumericMatrix &D,

 const Rcpp::NumericMatrix &S, unsigned nFactor, unsigned nEquil,

-unsigned nSample, const std::string &fileName, const std::string &cptFile)

+unsigned nSample, const std::string &cptFile)

 {   

     GapsRunner runner(D, S, nFactor, nEquil, nSample, cptFile);

     return runner.run();

@@ -19,6 +19,15 @@ mPSampler(data, data.pmax(0.1f), nPatterns, alphaP, maxGibbsMassP),

 mStatistics(data.nRow(), data.nCol(), nPatterns),

 mSamplePhase(false), mNumUpdatesA(0), mNumUpdatesP(0)

 {

+    if (mFixedMatrix == 'A')

+    {

+        mASampler.setMatrix(ColMatrix(FP));

+    }

+    else if (mFixedMatrix == 'P')

+    {

+        mPSampler.setMatrix(RowMatrix(FP));

+    }

+

     mASampler.sync(mPSampler);

     mPSampler.sync(mASampler);

 }

@@ -100,4 +109,4 @@ void GapsRunner::setUncertainty(const RowMatrix &S)

 {

     mASampler.setUncertainty(S);

     mPSampler.setUncertainty(S);

-}

\ No newline at end of file

+}

@@ -1,11 +1,11 @@

 #include "GapsStatistics.h"

 #include "math/Algorithms.h"

-GapsStatistics::GapsStatistics(unsigned nRow, unsigned nCol, unsigned nFactor)

+GapsStatistics::GapsStatistics(unsigned nRow, unsigned nCol, unsigned nFactor, PumpThreshold t)

     : mAMeanMatrix(nRow, nFactor), mAStdMatrix(nRow, nFactor),

         mPMeanMatrix(nFactor, nCol), mPStdMatrix(nFactor, nCol),

         mStatUpdates(0), mNumPatterns(nFactor), mPumpMatrix(nRow, nCol),

-        mPumpThreshold(PUMP_CUT), mPumpStatUpdates(0)

+        mPumpThreshold(t), mPumpStatUpdates(0)

 {}

 void GapsStatistics::update(const AmplitudeGibbsSampler &ASampler,

@@ -45,15 +45,9 @@ static unsigned geneThreshold(const ColMatrix &rankMatrix, unsigned pat)

     return static_cast<unsigned>(cutRank);

 }

-void GapsStatistics::updatePump(const AmplitudeGibbsSampler &ASampler,

-const PatternGibbsSampler &PSampler)

+void GapsStatistics::patternMarkers(ColMatrix normedA, RowMatrix normedP,

+ColMatrix &statMatrix)

 {

-    mPumpStatUpdates++;

-

-    // copy matrices so we can modify locally

-    ColMatrix normedA(ASampler.mMatrix);

-    RowMatrix normedP(PSampler.mMatrix);

-

     // helpful notation

     unsigned nGenes = normedA.nRow();

     unsigned nPatterns = normedA.nCol();

@@ -101,7 +95,7 @@ const PatternGibbsSampler &PSampler)

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

         {

             unsigned minNdx = gaps::algo::whichMin(rsStat.getRow(i));

-            mPumpMatrix(i,minNdx)++;

+            statMatrix(i,minNdx)++;

         }

     }

     else if (mPumpThreshold == PUMP_CUT)

@@ -119,13 +113,20 @@ const PatternGibbsSampler &PSampler)

             {

                 if (rankMatrix(i,j) <= cutRank)

                 {

-                    mPumpMatrix(i,j)++;

+                    statMatrix(i,j)++;

                 }

             }

         }

     }

 }

+void GapsStatistics::updatePump(const AmplitudeGibbsSampler &ASampler,

+const PatternGibbsSampler &PSampler)

+{

+    mPumpStatUpdates++;

+    patternMarkers(ASampler.mMatrix, PSampler.mMatrix, mPumpMatrix);

+}

+

 ColMatrix GapsStatistics::AMean() const

 {

     return mAMeanMatrix / mStatUpdates;

@@ -157,6 +158,21 @@ float GapsStatistics::meanChiSq(const AmplitudeGibbsSampler &ASampler) const

         M);

 }

+Rcpp::NumericMatrix GapsStatistics::pumpMatrix() const

+{

+    unsigned denom = mPumpStatUpdates != 0 ? mPumpStatUpdates : 1.f;

+    return (mPumpMatrix / denom).rMatrix();

+}

+

+Rcpp::NumericMatrix GapsStatistics::meanPattern()

+{

+    ColMatrix Amean(mAMeanMatrix / static_cast<float>(mStatUpdates));

+    RowMatrix Pmean(mPMeanMatrix / static_cast<float>(mStatUpdates));

+    ColMatrix mat(mAMeanMatrix.nRow(), mAMeanMatrix.nCol());

+    patternMarkers(Amean, Pmean, mat);

+    return mat.rMatrix();

+}

+

 Archive& operator<<(Archive &ar, GapsStatistics &stat)

 {

     ar << stat.mAMeanMatrix << stat.mAStdMatrix << stat.mPMeanMatrix

@@ -28,7 +28,7 @@ private:

 public:

-    GapsStatistics(unsigned nRow, unsigned nCol, unsigned nFactor);

+    GapsStatistics(unsigned nRow, unsigned nCol, unsigned nFactor, PumpThreshold t=PUMP_CUT);

     ColMatrix AMean() const;

     RowMatrix PMean() const;

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

     void updatePump(const AmplitudeGibbsSampler &ASampler,

         const PatternGibbsSampler &PSampler);

+    void patternMarkers(ColMatrix normedA, RowMatrix normedP, ColMatrix &statMatrix);

+

     // serialization

     friend Archive& operator<<(Archive &ar, GapsStatistics &stat);

     friend Archive& operator>>(Archive &ar, GapsStatistics &stat);

@@ -5,9 +5,6 @@ AmplitudeGibbsSampler::AmplitudeGibbsSampler(const RowMatrix &D,

 const RowMatrix &S, unsigned nFactor, float alpha, float maxGibbsMass)

     : GibbsSampler(D, S, D.nRow(), nFactor, alpha)

 {

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

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

-    mMaxGibbsMass = maxGibbsMass / mLambda;

     mQueue.setDimensionSize(mBinSize, mNumCols);

 }

@@ -102,9 +99,6 @@ PatternGibbsSampler::PatternGibbsSampler(const RowMatrix &D,

 const RowMatrix &S, unsigned nFactor, float alpha, float maxGibbsMass)

     : GibbsSampler(D, S, nFactor, D.nCol(), alpha)

 {

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

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

-    mMaxGibbsMass = maxGibbsMass / mLambda;

     mQueue.setDimensionSize(mBinSize , mNumRows);

 }

@@ -18,6 +18,15 @@ class GapsStatistics;

 /************************** GIBBS SAMPLER INTERFACE **************************/

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

+class GibbsSampler;

+

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

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

+

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

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

+

 template <class T, class MatA, class MatB>

 class GibbsSampler

 {

@@ -99,11 +108,11 @@ public:

     uint64_t nAtoms() const;

     void setUncertainty(const RowMatrix &S) { mSMatrix = MatB(S); }

+    void setUncertainty(const std::string &path) { mSMatrix = MatB(path); }

-    void setUncertainty(FileParser &p)

-    {

-        mSMatrix = MatB(p);

-    }

+    // serialization

+    friend Archive& operator<< <T, MatA, MatB> (Archive &ar, GibbsSampler &samp);

+    friend Archive& operator>> <T, MatA, MatB> (Archive &ar, GibbsSampler &samp);

 };

 class AmplitudeGibbsSampler : public GibbsSampler<AmplitudeGibbsSampler, ColMatrix, RowMatrix>

@@ -203,6 +212,11 @@ mAvgQueue(0.f), mNumQueues(0.f)

         % static_cast<uint64_t>(mNumRows * mNumCols);

     mQueue.setDomainSize(std::numeric_limits<uint64_t>::max() - remain);

     mDomain.setDomainSize(std::numeric_limits<uint64_t>::max() - remain);

+

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

+        gaps::algo::mean(mDMatrix);

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

+    mMaxGibbsMass = maxGibbsMass / mLambda;

 }

 template <class T, class MatA, class MatB>

@@ -503,4 +517,28 @@ uint64_t GibbsSampler<T, MatA, MatB>::nAtoms() const

     return mDomain.size();

 }

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

+void GibbsSampler<T, MatA, MatB>::setMatrix(const MatA &mat)

+{   

+    mMatrix = mat;

+}

+

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

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

+{

+    ar << samp.mMatrix << samp.mAPMatrix << samp.mQueue << samp.mDomain << samp.mLambda << samp.mMaxGibbsMass

+        << samp.mAnnealingTemp << samp.mNumRows << samp.mNumCols << samp.mBinSize << samp.mAvgQueue

+        << samp.mNumQueues;

+    return ar;

+}

+

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

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

+{

+    ar >> samp.mMatrix >> samp.mAPMatrix >> samp.mQueue >> samp.mDomain >> samp.mLambda >> samp.mMaxGibbsMass

+        >> samp.mAnnealingTemp >> samp.mNumRows >> samp.mNumCols >> samp.mBinSize >> samp.mAvgQueue

+        >> samp.mNumQueues;

+    return ar;

+}

+

 #endif

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

     mQueue.clear();

     mUsedPositions.clear();

     mUsedIndices.clear();

+    GAPS_ASSERT(mMinAtoms == mMaxAtoms);

 }

 unsigned ProposalQueue::size() const

@@ -136,7 +137,6 @@ bool ProposalQueue::birth(AtomicDomain &domain)

     mUsedIndices.insert(pos / mDimensionSize);

     mUsedPositions.insert(pos);

     domain.insert(pos, 0.f);

-    ++mMinAtoms; // TODO could be rejected

     ++mMaxAtoms;

     return true;

 }

@@ -5,6 +5,7 @@

 using namespace Rcpp;

+<<<<<<< HEAD

 // cogapsFromFile_cpp

 Rcpp::List cogapsFromFile_cpp(const std::string& D, unsigned nPatterns, unsigned maxIter, unsigned outputFrequency, unsigned seed, float alphaA, float alphaP, float maxGibbsMassA, float maxGibbsMassP, bool messages, bool singleCellRNASeq);

 RcppExport SEXP _CoGAPS_cogapsFromFile_cpp(SEXP DSEXP, SEXP nPatternsSEXP, SEXP maxIterSEXP, SEXP outputFrequencySEXP, SEXP seedSEXP, SEXP alphaASEXP, SEXP alphaPSEXP, SEXP maxGibbsMassASEXP, SEXP maxGibbsMassPSEXP, SEXP messagesSEXP, SEXP singleCellRNASeqSEXP) {

@@ -16,12 +17,28 @@ BEGIN_RCPP

     Rcpp::traits::input_parameter< unsigned >::type maxIter(maxIterSEXP);

     Rcpp::traits::input_parameter< unsigned >::type outputFrequency(outputFrequencySEXP);

     Rcpp::traits::input_parameter< unsigned >::type seed(seedSEXP);

+=======

+// cogaps_cpp

+Rcpp::List cogaps_cpp(const Rcpp::NumericMatrix& D, const Rcpp::NumericMatrix& S, unsigned nFactor, unsigned nEquil, unsigned nEquilCool, unsigned nSample, unsigned nOutputs, 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 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;

+    Rcpp::traits::input_parameter< const Rcpp::NumericMatrix& >::type D(DSEXP);

+    Rcpp::traits::input_parameter< const Rcpp::NumericMatrix& >::type S(SSEXP);

+    Rcpp::traits::input_parameter< unsigned >::type nFactor(nFactorSEXP);

+    Rcpp::traits::input_parameter< unsigned >::type nEquil(nEquilSEXP);

+    Rcpp::traits::input_parameter< unsigned >::type nEquilCool(nEquilCoolSEXP);

+    Rcpp::traits::input_parameter< unsigned >::type nSample(nSampleSEXP);

+    Rcpp::traits::input_parameter< unsigned >::type nOutputs(nOutputsSEXP);

+>>>>>>> develop

     Rcpp::traits::input_parameter< float >::type alphaA(alphaASEXP);

     Rcpp::traits::input_parameter< float >::type alphaP(alphaPSEXP);

     Rcpp::traits::input_parameter< float >::type maxGibbsMassA(maxGibbsMassASEXP);

     Rcpp::traits::input_parameter< float >::type maxGibbsMassP(maxGibbsMassPSEXP);

     Rcpp::traits::input_parameter< bool >::type messages(messagesSEXP);

     Rcpp::traits::input_parameter< bool >::type singleCellRNASeq(singleCellRNASeqSEXP);

+<<<<<<< HEAD

     rcpp_result_gen = Rcpp::wrap(cogapsFromFile_cpp(D, nPatterns, maxIter, outputFrequency, seed, alphaA, alphaP, maxGibbsMassA, maxGibbsMassP, messages, singleCellRNASeq));

     return rcpp_result_gen;

 END_RCPP

@@ -29,10 +46,27 @@ END_RCPP

 // cogaps_cpp

 Rcpp::List cogaps_cpp(const Rcpp::NumericMatrix& D, unsigned nPatterns, unsigned maxIter, unsigned outputFrequency, unsigned seed, float alphaA, float alphaP, float maxGibbsMassA, float maxGibbsMassP, bool messages, bool singleCellRNASeq, unsigned nCores);

 RcppExport SEXP _CoGAPS_cogaps_cpp(SEXP DSEXP, SEXP nPatternsSEXP, SEXP maxIterSEXP, SEXP outputFrequencySEXP, SEXP seedSEXP, SEXP alphaASEXP, SEXP alphaPSEXP, SEXP maxGibbsMassASEXP, SEXP maxGibbsMassPSEXP, SEXP messagesSEXP, SEXP singleCellRNASeqSEXP, SEXP nCoresSEXP) {

+=======

+    Rcpp::traits::input_parameter< char >::type whichMatrixFixed(whichMatrixFixedSEXP);

+    Rcpp::traits::input_parameter< const Rcpp::NumericMatrix& >::type FP(FPSEXP);

+    Rcpp::traits::input_parameter< unsigned >::type checkpointInterval(checkpointIntervalSEXP);

+    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::traits::input_parameter< unsigned >::type nCores(nCoresSEXP);

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

+    return rcpp_result_gen;

+END_RCPP

+}

+// cogapsFromCheckpoint_cpp

+Rcpp::List cogapsFromCheckpoint_cpp(const Rcpp::NumericMatrix& D, const Rcpp::NumericMatrix& S, unsigned nFactor, unsigned nEquil, unsigned nSample, const std::string& cptFile);

+RcppExport SEXP _CoGAPS_cogapsFromCheckpoint_cpp(SEXP DSEXP, SEXP SSEXP, SEXP nFactorSEXP, SEXP nEquilSEXP, SEXP nSampleSEXP, SEXP cptFileSEXP) {

+>>>>>>> develop

 BEGIN_RCPP

     Rcpp::RObject rcpp_result_gen;

     Rcpp::RNGScope rcpp_rngScope_gen;

     Rcpp::traits::input_parameter< const Rcpp::NumericMatrix& >::type D(DSEXP);

+<<<<<<< HEAD

     Rcpp::traits::input_parameter< unsigned >::type nPatterns(nPatternsSEXP);

     Rcpp::traits::input_parameter< unsigned >::type maxIter(maxIterSEXP);

     Rcpp::traits::input_parameter< unsigned >::type outputFrequency(outputFrequencySEXP);

@@ -45,6 +79,14 @@ BEGIN_RCPP

     Rcpp::traits::input_parameter< bool >::type singleCellRNASeq(singleCellRNASeqSEXP);

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

     rcpp_result_gen = Rcpp::wrap(cogaps_cpp(D, nPatterns, maxIter, outputFrequency, seed, alphaA, alphaP, maxGibbsMassA, maxGibbsMassP, messages, singleCellRNASeq, nCores));

+=======

+    Rcpp::traits::input_parameter< const Rcpp::NumericMatrix& >::type S(SSEXP);

+    Rcpp::traits::input_parameter< unsigned >::type nFactor(nFactorSEXP);

+    Rcpp::traits::input_parameter< unsigned >::type nEquil(nEquilSEXP);

+    Rcpp::traits::input_parameter< unsigned >::type nSample(nSampleSEXP);

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

+    rcpp_result_gen = Rcpp::wrap(cogapsFromCheckpoint_cpp(D, S, nFactor, nEquil, nSample, cptFile));

+>>>>>>> develop

     return rcpp_result_gen;

 END_RCPP

 }

@@ -70,8 +112,14 @@ END_RCPP

 }

 static const R_CallMethodDef CallEntries[] = {

+<<<<<<< HEAD

     {"_CoGAPS_cogapsFromFile_cpp", (DL_FUNC) &_CoGAPS_cogapsFromFile_cpp, 11},

     {"_CoGAPS_cogaps_cpp", (DL_FUNC) &_CoGAPS_cogaps_cpp, 12},

+=======

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

+    {"_CoGAPS_cogapsFromCheckpoint_cpp", (DL_FUNC) &_CoGAPS_cogapsFromCheckpoint_cpp, 6},

+    {"_CoGAPS_cogapsFromFile_cpp", (DL_FUNC) &_CoGAPS_cogapsFromFile_cpp, 1},

+>>>>>>> develop

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

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

     {NULL, NULL, 0}

@@ -1,2 +1,3 @@

 #include "catch.h"

 #include "../GapsRunner.h"

+

@@ -5,8 +5,8 @@

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

 #include "Vector.h"

-#include <vector>

 #include <algorithm>

+#include <vector>

 // forward declarations

 class RowMatrix;

@@ -9,14 +9,14 @@ struct MatrixElement

     unsigned dim[2];

     float val;

-    MatrixElement(unsigned r, unsigned c, float v)

+    MatrixElement(unsigned r, unsigned c, float v) // NOLINT

         : val(v)

     {

         dim[0] = r;

         dim[1] = c;

     }

-    MatrixElement(unsigned r, unsigned c, const std::string &s)

+    MatrixElement(unsigned r, unsigned c, const std::string &s) // NOLINT

         :  val(0.f)

     {

         dim[0] = r;

@@ -27,17 +27,17 @@ struct MatrixElement

     unsigned operator[](unsigned i)

     {

-        return dim[i];

+        return dim[i]; // NOLINT

     }

     unsigned row() const

     {

-        return dim[0];

+        return dim[0]; // NOLINT

     }

     unsigned col() const

     {

-        return dim[1];

+        return dim[1]; // NOLINT

     }

     float value() const

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

-#include "MtxParser.h"

 #include "MatrixElement.h"

+#include "MtxParser.h"

+

 #include "../GapsAssert.h"

 #include <sstream>

@@ -20,8 +20,8 @@

 #include <boost/random/uniform_01.hpp>

 #include <boost/random/uniform_real_distribution.hpp>

-#include <stdint.h>

 #include <algorithm>

+#include <stdint.h>

 #ifdef __GAPS_OPENMP__

     #include <omp.h>

@@ -108,20 +108,20 @@ public:

 #if defined( __GAPS_AVX512__ )

     float scalar()

     {

-        float* ra = reinterpret_cast<float*>(&mData);

+        float* ra = reinterpret_cast<float*>(&mData); // NOLINT

         return ra[0] + ra[1] + ra[2] + ra[3] + ra[4] + ra[5] + ra[6] + ra[7] +

             ra[8] + ra[9] + ra[10] + ra[11] + ra[12] + ra[13] + ra[14] + ra[15];

     }

 #elif defined( __GAPS_AVX__ )

     float scalar()

     {

-        float* ra = reinterpret_cast<float*>(&mData);

+        float* ra = reinterpret_cast<float*>(&mData); // NOLINT

         return ra[0] + ra[1] + ra[2] + ra[3] + ra[4] + ra[5] + ra[6] + ra[7];

     }

 #elif defined( __GAPS_SSE__ )

     float scalar()

     {

-        float* ra = reinterpret_cast<float*>(&mData);

+        float* ra = reinterpret_cast<float*>(&mData); // NOLINT

         return ra[0] + ra[1] + ra[2] + ra[3];

     }

 #else