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

 # Generated by roxygen2: do not edit by hand

 export(CoGAPS)

-export(CoGAPSFromFile)

+export(CoGapsFromCheckpoint)

 export(GWCoGAPS)

 export(GWCoGapsFromCheckpoint)

 export(binaryA)

@@ -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

@@ -43,7 +42,7 @@

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

 #' @export

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

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

+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, ...)

@@ -99,7 +98,7 @@ checkpointFile="gaps_checkpoint.out", nCores=1, ...)

     # run algorithm with call to C++ code

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

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

+        alphaA, alphaP, maxGibbmassA, maxGibbmassP, seed, messages,

         singleCellRNASeq, whichMatrixFixed, fixedPatterns, checkpointInterval,

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

         nCores)

@@ -120,21 +119,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,12 @@

 # 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, 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)

+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, fileName, cptFile) {

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

+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) {

@@ -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

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

 # include macros

-m4_builtin([include], ax_check_compile_flag.m4)

-m4_builtin([include], ax_compiler_vendor.m4)

-m4_builtin([include], ax_compiler_version.m4)

-m4_builtin([include], ax_openmp.m4)

+m4_builtin([include], m4/ax_check_compile_flag.m4)

+m4_builtin([include], m4/ax_compiler_vendor.m4)

+m4_builtin([include], m4/ax_compiler_version.m4)

+m4_builtin([include], m4/ax_openmp.m4)

 # get version of CoGAPS from DESCRIPTION file

 AC_INIT(CoGAPS, m4_esyscmd_s([awk -e '/^Version:/ {print $2}' Repo/DESCRIPTION]))

@@ -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

@@ -4,12 +4,12 @@

 \alias{CoGAPS}

 \title{CoGAPS Matrix Factorization Algorithm}

 \usage{

-CoGAPS(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", nCores = 1, ...)

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

+  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, ...)

 }

 \arguments{

 \item{D}{data matrix}

@@ -25,8 +25,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}{path to checkpoint file}

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

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

 }

 \value{

 list with A and P matrix estimates

@@ -25,4 +29,8 @@ Restart CoGAPS from Checkpoint File

 loads the state of a previous CoGAPS run from a file and

  continues the run from that point

 }

+\examples{

+data(SimpSim)

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

+}

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

 #include <Rcpp.h>

 #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)

@@ -40,6 +42,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)

     {

@@ -47,12 +91,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

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

 // [[Rcpp::export]]

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

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

-unsigned nEquilCool, unsigned nSample, unsigned nOutputs, unsigned nSnapshots,

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

@@ -16,23 +16,24 @@ 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,

+        nOutputs, alphaA, alphaP, maxGibbmassA, maxGibbmassP, seed,

         messages, singleCellRNASeq, checkpointInterval, cptFile,

-        whichMatrixFixed, FP, nCores);

+        whichMatrixFixed, FP, nCores, static_cast<PumpThreshold>(pumpThreshold),

+        nPumpSamples);

     return runner.run();

 }

 // [[Rcpp::export]]

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

 }

 // [[Rcpp::export]]

-void cogapsFromFile_cpp(const std::string &D)

+void cogapsFromFile_cpp(const std::string &D) // NOLINT

 {

     // TODO implement

 }

@@ -27,8 +27,7 @@

 #define GAPS_ERROR(msg)                             \

     do {                                            \

-        Rcpp::Rcout << "error: " << msg << '\n';    \

-        Rcpp::stop("CoGAPS aborted");                              \

+        Rcpp::stop(msg);                            \

     } while(0)

 #endif

\ No newline at end of file

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

+#include "GapsAssert.h"

 #include "GapsRunner.h"

-#include "math/SIMD.h"

 #include "math/Random.h"

-#include "GapsAssert.h"

+#include "math/SIMD.h"

 #include <algorithm>

 #ifdef __GAPS_OPENMP__

@@ -34,10 +34,11 @@ static std::vector< std::vector<unsigned> > sampleIndices(unsigned n, unsigned n

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

 unsigned nFactor, unsigned nEquil, unsigned nCool, unsigned nSample,

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

+unsigned nOutputs, 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, unsigned nCores)

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

+PumpThreshold pumpThreshold, unsigned nPumpSamples)

     :

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

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

@@ -46,11 +47,21 @@ mSampleIter(nSample), mNumOutputs(nOutputs),

 mPrintMessages(messages), mCurrentIter(0), mPhase(GAPS_BURN), mSeed(seed),

 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),

-mNumCores(nCores)

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

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

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

+mNumCores(nCores), mFixedMatrix(whichMatrixFixed),

+mNumPumpSamples(nPumpSamples)

 {

+    if (mFixedMatrix == 'A')

+    {

+        mASampler.setMatrix(ColMatrix(FP));

+    }

+    else if (mFixedMatrix == 'P')

+    {

+        mPSampler.setMatrix(RowMatrix(FP));

+    }

+

     mASampler.sync(mPSampler);

     mPSampler.sync(mASampler);

     gaps::random::setSeed(seed);

@@ -61,20 +72,21 @@ unsigned nFactor, unsigned nEquil, unsigned nSample, const std::string &cptFile)

     :

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

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

-mASampler(D, S, nFactor), mPSampler(D, S, nFactor),

+mCheckpointFile(cptFile), mASampler(D, S, nFactor), mPSampler(D, S, nFactor),

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

 {

     Archive ar(cptFile, ARCHIVE_READ);

     gaps::random::load(ar);

+    unsigned storedPhase = 0;

     ar >> mChiSqEquil >> mNumAAtomsEquil >> mNumPAtomsEquil >> mChiSqSample

         >> mNumAAtomsSample >> mNumPAtomsSample >> mIterA >> mIterP

         >> mEquilIter >> mCoolIter >> mSampleIter >> mNumOutputs

-        >> mPrintMessages >> mCurrentIter >> mPhase >> mSeed >> mLastCheckpoint

-        >> mCheckpointInterval >> mCheckpointFile >> mNumUpdatesA

-        >> mNumUpdatesP >> mASampler >> mPSampler >> mStatistics >> mNumCores

-        >> mStartTime;

+        >> mPrintMessages >> mCurrentIter >> storedPhase >> mSeed

+        >> mCheckpointInterval >> mNumUpdatesA

+        >> mNumUpdatesP >> mASampler >> mPSampler >> mStatistics >> mNumCores;

+    mPhase = static_cast<GapsPhase>(storedPhase);

     mASampler.sync(mPSampler);

     mPSampler.sync(mASampler);

 }

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

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

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

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

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

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

+        Rcpp::Named("pumpStats") = mStatistics.pumpMatrix(),

+        Rcpp::Named("meanPatternAssignment") = mStatistics.meanPattern()

     );

 }

@@ -173,6 +187,10 @@ void GapsRunner::runSampPhase()

         makeCheckpointIfNeeded();

         updateSampler();

         mStatistics.update(mASampler, mPSampler);

+        if (mNumPumpSamples != 0 && ((mCurrentIter + 1) % (mSampleIter / mNumPumpSamples)) == 0)

+        {

+            mStatistics.updatePump(mASampler, mPSampler);

+        }

         storeSamplerInfo(mNumAAtomsSample, mNumPAtomsSample, mChiSqSample);

         displayStatus("Samp: ", mSampleIter);

     }

@@ -215,13 +233,19 @@ double GapsRunner::estPercentComplete()

 void GapsRunner::updateSampler()

 {

-    mNumUpdatesA += mIterA;

-    mASampler.update(mIterA, mNumCores);

-    mPSampler.sync(mASampler);

+    if (mFixedMatrix != 'A')

+    {

+        mNumUpdatesA += mIterA;

+        mASampler.update(mIterA, mNumCores);

+        mPSampler.sync(mASampler);

+    }

-    mNumUpdatesP += mIterP;

-    mPSampler.update(mIterP, mNumCores);

-    mASampler.sync(mPSampler);

+    if (mFixedMatrix != 'P')

+    {

+        mNumUpdatesP += mIterP;

+        mPSampler.update(mIterP, mNumCores);

+        mASampler.sync(mPSampler);

+    }

 }

 void GapsRunner::storeSamplerInfo(Vector &atomsA, Vector &atomsP, Vector &chi2)

@@ -280,10 +304,10 @@ void GapsRunner::createCheckpoint()

     ar << mChiSqEquil << mNumAAtomsEquil << mNumPAtomsEquil << mChiSqSample

         << mNumAAtomsSample << mNumPAtomsSample << mIterA << mIterP

         << mEquilIter << mCoolIter << mSampleIter << mNumOutputs

-        << mPrintMessages << mCurrentIter << mPhase << mSeed << mLastCheckpoint

-        << mCheckpointInterval << mCheckpointFile << mNumUpdatesA

-        << mNumUpdatesP << mASampler << mPSampler << mStatistics << mNumCores

-        << mStartTime;

+        << mPrintMessages << mCurrentIter << static_cast<unsigned>(mPhase)

+        << mSeed 

+        << mCheckpointInterval << mNumUpdatesA

+        << mNumUpdatesP << mASampler << mPSampler << mStatistics << mNumCores;

     ar.close();

     // display time it took to create checkpoint

@@ -17,9 +17,9 @@ namespace bpt = boost::posix_time;

 enum GapsPhase

 {

-    GAPS_BURN,

-    GAPS_COOL,

-    GAPS_SAMP

+    GAPS_BURN=0,

+    GAPS_COOL=1,

+    GAPS_SAMP=2

 };

 // Manages all data and parameters for running CoGAPS; contains the high-level

@@ -69,6 +69,10 @@ public:

     bpt::ptime mStartTime;

+    char mFixedMatrix;

+

+    unsigned mNumPumpSamples;

+

     void createCheckpoint();

     void makeCheckpointIfNeeded();

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

@@ -84,10 +88,11 @@ public:

     // construct from parameters

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

         unsigned nFactor, unsigned nEquil, unsigned nCool, unsigned nSample,

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

+        unsigned nOutputs, 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, unsigned nCores);

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

+        PumpThreshold pumpThreshold, unsigned nPumpSamples);

     // construct from checkpoint file

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

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

+}

+

 Rcpp::NumericMatrix GapsStatistics::AMean() const

 {

     return (mAMeanMatrix / mStatUpdates).rMatrix();

@@ -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,12 +28,14 @@ private:

 public:

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

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

     Rcpp::NumericMatrix AMean() const;

     Rcpp::NumericMatrix AStd() const;

     Rcpp::NumericMatrix PMean() const;

     Rcpp::NumericMatrix PStd() const;

+    Rcpp::NumericMatrix pumpMatrix() const;

+    Rcpp::NumericMatrix meanPattern();

     float meanChiSq(const AmplitudeGibbsSampler &ASampler) const;

@@ -43,6 +45,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);

@@ -2,12 +2,10 @@

 #include "math/SIMD.h"

 AmplitudeGibbsSampler::AmplitudeGibbsSampler(const Rcpp::NumericMatrix &D,

-const Rcpp::NumericMatrix &S, unsigned nFactor, float alpha, float maxGibbsMass)

-    : GibbsSampler(D, S, D.nrow(), nFactor, alpha)

+const Rcpp::NumericMatrix &S, unsigned nFactor, float alpha, float maxGibbsMass,

+bool singleCell)

+    : GibbsSampler(D, S, D.nrow(), nFactor, alpha, maxGibbsMass, singleCell, nFactor)

 {

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

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

-    mMaxGibbsMass = maxGibbsMass / mLambda;

     mQueue.setDimensionSize(mBinSize, mNumCols);

 }

@@ -99,12 +97,10 @@ unsigned r2, unsigned c2, float m2)

 }

 PatternGibbsSampler::PatternGibbsSampler(const Rcpp::NumericMatrix &D,

-const Rcpp::NumericMatrix &S, unsigned nFactor, float alpha, float maxGibbsMass)

-    : GibbsSampler(D, S, nFactor, D.ncol(), alpha)

+const Rcpp::NumericMatrix &S, unsigned nFactor, float alpha, float maxGibbsMass,

+bool singleCell)

+    : GibbsSampler(D, S, nFactor, D.ncol(), alpha, maxGibbsMass, singleCell, nFactor)

 {

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

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

-    mMaxGibbsMass = maxGibbsMass / mLambda;

     mQueue.setDimensionSize(mBinSize , mNumRows);

 }

@@ -20,6 +20,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

 {

@@ -74,7 +83,8 @@ protected:

 public:

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

-        unsigned nrow, unsigned ncol, float alpha);

+        unsigned nrow, unsigned ncol, float alpha, float maxGibbsMass,

+        bool singleCell, unsigned nFactor);

     void update(unsigned nSteps, unsigned nCores);

     void setAnnealingTemp(float temp);

@@ -82,6 +92,12 @@ public:

     float chi2() const;

     uint64_t nAtoms() const;

+

+    void setMatrix(const MatA &mat);

+

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

@@ -101,7 +117,7 @@ public:

     AmplitudeGibbsSampler(const Rcpp::NumericMatrix &D,

         const Rcpp::NumericMatrix &S, unsigned nFactor, float alpha=0.f,

-        float maxGibbsMass=0.f);

+        float maxGibbsMass=0.f, bool singleCell=false);

     void sync(PatternGibbsSampler &sampler);

@@ -131,7 +147,7 @@ public:

     PatternGibbsSampler(const Rcpp::NumericMatrix &D,

         const Rcpp::NumericMatrix &S, unsigned nFactor, float alpha=0.f,

-        float maxGibbsMass=0.f);

+        float maxGibbsMass=0.f, bool singleCell=false);

     void sync(AmplitudeGibbsSampler &sampler);

@@ -148,7 +164,8 @@ public:

 template <class T, class MatA, class MatB>

 GibbsSampler<T, MatA, MatB>::GibbsSampler(const Rcpp::NumericMatrix &D,

-const Rcpp::NumericMatrix &S, unsigned nrow, unsigned ncol, float alpha)

+const Rcpp::NumericMatrix &S, unsigned nrow, unsigned ncol, float alpha,

+float maxGibbsMass, bool singleCell, unsigned nFactor)

     :

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

 mAPMatrix(D.nrow(), D.ncol()), mQueue(nrow * ncol, alpha), mLambda(0.f),

@@ -161,6 +178,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>

@@ -461,4 +483,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;

 }

@@ -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, 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) {

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

@@ -18,7 +18,6 @@ BEGIN_RCPP

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

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

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

-    Rcpp::traits::input_parameter< unsigned >::type nSnapshots(nSnapshotsSEXP);

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

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

     Rcpp::traits::input_parameter< float >::type maxGibbmassA(maxGibbmassASEXP);

@@ -33,13 +32,13 @@ BEGIN_RCPP

     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, nSnapshots, alphaA, alphaP, maxGibbmassA, maxGibbmassP, seed, messages, singleCellRNASeq, whichMatrixFixed, FP, checkpointInterval, cptFile, pumpThreshold, nPumpSamples, nCores));

+    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& fileName, const std::string& cptFile);

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

+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) {

 BEGIN_RCPP

     Rcpp::RObject rcpp_result_gen;

     Rcpp::RNGScope rcpp_rngScope_gen;

@@ -48,9 +47,8 @@ BEGIN_RCPP

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

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

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

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

     return rcpp_result_gen;

 END_RCPP

 }

@@ -86,8 +84,8 @@ END_RCPP

 }

 static const R_CallMethodDef CallEntries[] = {

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

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

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

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