@@ -32,4 +32,5 @@

 ^src/cpp_tests/testGibbsSampler.o

 ^src/cpp_tests/testMatrix.o

 ^src/cpp_tests/testProposalQueue.o

-^src/cpp_tests/testRandom.o

\ No newline at end of file

+^src/cpp_tests/testRandom.o

+^src/cpp_tests/testSerialization.o

\ No newline at end of file

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

 #' data(SimpSim)

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

 #' @export

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

+CoGAPS <- function(D, S, nFactor=7, nEquil=250, nSample=250, 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, 

@@ -114,6 +114,7 @@ checkpointFile="gaps_checkpoint.out", nCores=1, ...)

 }

 #' Restart CoGAPS from Checkpoint File

+#' @export

 #'

 #' @details loads the state of a previous CoGAPS run from a file and

 #'  continues the run from that point

@@ -122,6 +123,9 @@ checkpointFile="gaps_checkpoint.out", nCores=1, ...)

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

 {

     if (is.na(checkpointFile))

@@ -103,4 +103,4 @@ NULL

 #' @docType data

 #' @name tf2ugFC

 #' @usage tf2ugFC

-NULL

\ No newline at end of file

+NULL

@@ -46,4 +46,4 @@ plotGAPS <- function(A, P, outputPDF="")

     {

         dev.off()

     }

-}

\ No newline at end of file

+}

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

 ac_user_opts='

 enable_option_checking

 enable_debug

+enable_warnings

 enable_simd

 '

       ac_precious_vars='build_alias

@@ -1272,6 +1273,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-warnings       compile CoGAPS with warning messages

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

 Some influential environment variables:

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

 fi

+# Check if compiler warnings should be turned on

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

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

+  enableval=$enable_warnings; warnings=yes

+else

+  warnings=no

+fi

+

+

 # Check if specific version of SIMD instructions requested

 # Check whether --enable-simd was given.

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

@@ -4167,6 +4178,10 @@ echo "building on $ax_cv_cxx_compiler_vendor compiler version $ax_cv_cxx_compile

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

     echo "Building Debug Version of CoGAPS"

     GAPS_CPP_FLAGS+=" -DGAPS_DEBUG "

+fi

+

+# set compile flags if warnings enabled

+if test "x$warnings" = "xyes" ; then

     { $as_echo "$as_me:${as_lineno-$LINENO}: checking whether C++ compiler accepts -Wall" >&5

 $as_echo_n "checking whether C++ compiler accepts -Wall... " >&6; }

 if ${ax_cv_check_cxxflags___Wall+:} false; then :

@@ -4272,6 +4287,76 @@ else

   :

 fi

+    { $as_echo "$as_me:${as_lineno-$LINENO}: checking whether C++ compiler accepts -Wno-unused-parameter" >&5

+$as_echo_n "checking whether C++ compiler accepts -Wno-unused-parameter... " >&6; }

+if ${ax_cv_check_cxxflags___Wno_unused_parameter+:} false; then :

+  $as_echo_n "(cached) " >&6

+else

+

+  ax_check_save_flags=$CXXFLAGS

+  CXXFLAGS="$CXXFLAGS  -Wno-unused-parameter"

+  cat confdefs.h - <<_ACEOF >conftest.$ac_ext

+/* end confdefs.h.  */

+

+int

+main ()

+{

+

+  ;

+  return 0;

+}

+_ACEOF

+if ac_fn_cxx_try_compile "$LINENO"; then :

+  ax_cv_check_cxxflags___Wno_unused_parameter=yes

+else

+  ax_cv_check_cxxflags___Wno_unused_parameter=no

+fi

+rm -f core conftest.err conftest.$ac_objext conftest.$ac_ext

+  CXXFLAGS=$ax_check_save_flags

+fi

+{ $as_echo "$as_me:${as_lineno-$LINENO}: result: $ax_cv_check_cxxflags___Wno_unused_parameter" >&5

+$as_echo "$ax_cv_check_cxxflags___Wno_unused_parameter" >&6; }

+if test "x$ax_cv_check_cxxflags___Wno_unused_parameter" = xyes; then :

+  GAPS_CXX_FLAGS+=" -Wno-unused-parameter "

+else

+  :

+fi

+

+    { $as_echo "$as_me:${as_lineno-$LINENO}: checking whether C++ compiler accepts -Wno-unused-function" >&5

+$as_echo_n "checking whether C++ compiler accepts -Wno-unused-function... " >&6; }

+if ${ax_cv_check_cxxflags___Wno_unused_function+:} false; then :

+  $as_echo_n "(cached) " >&6

+else

+

+  ax_check_save_flags=$CXXFLAGS

+  CXXFLAGS="$CXXFLAGS  -Wno-unused-function"

+  cat confdefs.h - <<_ACEOF >conftest.$ac_ext

+/* end confdefs.h.  */

+

+int

+main ()

+{

+

+  ;

+  return 0;

+}

+_ACEOF

+if ac_fn_cxx_try_compile "$LINENO"; then :

+  ax_cv_check_cxxflags___Wno_unused_function=yes

+else

+  ax_cv_check_cxxflags___Wno_unused_function=no

+fi

+rm -f core conftest.err conftest.$ac_objext conftest.$ac_ext

+  CXXFLAGS=$ax_check_save_flags

+fi

+{ $as_echo "$as_me:${as_lineno-$LINENO}: result: $ax_cv_check_cxxflags___Wno_unused_function" >&5

+$as_echo "$ax_cv_check_cxxflags___Wno_unused_function" >&6; }

+if test "x$ax_cv_check_cxxflags___Wno_unused_function" = xyes; then :

+  GAPS_CXX_FLAGS+=" -Wno-unused-function "

+else

+  :

+fi

+

 fi

 # set compile flags for SIMD

@@ -1,77 +1,87 @@

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

-

-# get version of CoGAPS from DESCRIPTION file

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

-

-# get C++ compiler from R configuration

-CXX=`"${R_HOME}/bin/R" CMD config CXX`

-

-# Switch to a C++ compiler, and check if it works.

-AC_LANG(C++)

-AC_REQUIRE_CPP

-AC_PROG_CXX

-

-# Check if compiling debug version

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

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

-

-# Check if specific version of SIMD instructions requested

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

-    [specify simd instruction set (sse, avx)])],

-    [simd_version=$enableval], [simd_version=sse])

-

-# default CoGAPS specific flags

-GAPS_CPP_FLAGS=" -DBOOST_MATH_PROMOTE_DOUBLE_POLICY=0 "

-GAPS_CXX_FLAGS=

-GAPS_LIBS=

-

-# get compiler info

-AX_COMPILER_VENDOR

-AX_COMPILER_VERSION

-AX_OPENMP

-

-# set openmp flags

-GAPS_CXX_FLAGS+=" $OPENMP_CXXFLAGS "

-GAPS_LIBS+=" $OPENMP_CXXFLAGS "

-

-echo "building on $ax_cv_cxx_compiler_vendor compiler version $ax_cv_cxx_compiler_version"

-

-# set compile flags for debug build

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

-    echo "Building Debug Version of CoGAPS"

-    GAPS_CPP_FLAGS+=" -DGAPS_DEBUG "

-    AX_CHECK_COMPILE_FLAG([-Wall], [GAPS_CXX_FLAGS+=" -Wall "])

-    AX_CHECK_COMPILE_FLAG([-Wextra], [GAPS_CXX_FLAGS+=" -Wextra "])

-    AX_CHECK_COMPILE_FLAG([-Werror], [GAPS_CXX_FLAGS+=" -Werror "])

-fi

-

-# set compile flags for SIMD

-if test "x$simd_version" = "xsse" ; then

-    echo "Compiling with SSE instructions"

-    AX_CHECK_COMPILE_FLAG([-msse4.2], [GAPS_CXX_FLAGS+=" -msse4.2 "])

-elif test "x$simd_version" = "xyes" ; then

-    echo "Compiling with SSE instructions"

-    AX_CHECK_COMPILE_FLAG([-msse4.2], [GAPS_CXX_FLAGS+=" -msse4.2 "])

-elif test "x$simd_version" = "xavx" ; then

-    echo "Compiling with AVX instructions"

-    AX_CHECK_COMPILE_FLAG([-mavx], [GAPS_CXX_FLAGS+=" -mavx "])    

-elif test "x$simd_version" = "xno" ; then

-    echo "Compiling without SIMD instructions"

-else

-    echo "Error: Invalid SIMD type"

-    exit -1

-fi

-

-# export variables containing flags

-AC_SUBST(GAPS_CPP_FLAGS)

-AC_SUBST(GAPS_CXX_FLAGS)

-AC_SUBST(GAPS_LIBS)

-

-# create makefile, output configure file

-AC_CONFIG_FILES([src/Makevars])

-AC_OUTPUT

-

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

+

+# get version of CoGAPS from DESCRIPTION file

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

+

+# get C++ compiler from R configuration

+CXX=`"${R_HOME}/bin/R" CMD config CXX`

+

+# Switch to a C++ compiler, and check if it works.

+AC_LANG(C++)

+AC_REQUIRE_CPP

+AC_PROG_CXX

+

+# Check if compiling debug version

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

+    [build debug version of CoGAPS])], [build_debug=yes], [build_debug=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])

+

+# Check if specific version of SIMD instructions requested

+AC_ARG_ENABLE(simd, [AC_HELP_STRING([--enable-simd],

+    [specify simd instruction set (sse, avx)])],

+    [simd_version=$enableval], [simd_version=sse])

+

+# default CoGAPS specific flags

+GAPS_CPP_FLAGS=" -DBOOST_MATH_PROMOTE_DOUBLE_POLICY=0 "

+GAPS_CXX_FLAGS=

+GAPS_LIBS=

+

+# get compiler info

+AX_COMPILER_VENDOR

+AX_COMPILER_VERSION

+AX_OPENMP

+

+# set openmp flags

+GAPS_CXX_FLAGS+=" $OPENMP_CXXFLAGS "

+GAPS_LIBS+=" $OPENMP_CXXFLAGS "

+

+echo "building on $ax_cv_cxx_compiler_vendor compiler version $ax_cv_cxx_compiler_version"

+

+# set compile flags for debug build

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

+    echo "Building Debug Version of CoGAPS"

+    GAPS_CPP_FLAGS+=" -DGAPS_DEBUG "

+fi

+

+# set compile flags if warnings enabled

+if test "x$warnings" = "xyes" ; then

+    AX_CHECK_COMPILE_FLAG([-Wall], [GAPS_CXX_FLAGS+=" -Wall "])

+    AX_CHECK_COMPILE_FLAG([-Wextra], [GAPS_CXX_FLAGS+=" -Wextra "])

+    AX_CHECK_COMPILE_FLAG([-Werror], [GAPS_CXX_FLAGS+=" -Werror "])

+    AX_CHECK_COMPILE_FLAG([-Wno-unused-parameter], [GAPS_CXX_FLAGS+=" -Wno-unused-parameter "])

+    AX_CHECK_COMPILE_FLAG([-Wno-unused-function], [GAPS_CXX_FLAGS+=" -Wno-unused-function "])

+fi

+

+# set compile flags for SIMD

+if test "x$simd_version" = "xsse" ; then

+    echo "Compiling with SSE instructions"

+    AX_CHECK_COMPILE_FLAG([-msse4.2], [GAPS_CXX_FLAGS+=" -msse4.2 "])

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

+    echo "Compiling with SSE instructions"

+    AX_CHECK_COMPILE_FLAG([-msse4.2], [GAPS_CXX_FLAGS+=" -msse4.2 "])

+elif test "x$simd_version" = "xavx" ; then

+    echo "Compiling with AVX instructions"

+    AX_CHECK_COMPILE_FLAG([-mavx], [GAPS_CXX_FLAGS+=" -mavx "])    

+elif test "x$simd_version" = "xno" ; then

+    echo "Compiling without SIMD instructions"

+else

+    echo "Error: Invalid SIMD type"

+    exit -1

+fi

+

+# export variables containing flags

+AC_SUBST(GAPS_CPP_FLAGS)

+AC_SUBST(GAPS_CXX_FLAGS)

+AC_SUBST(GAPS_LIBS)

+

+# create makefile, output configure file

+AC_CONFIG_FILES([src/Makevars])

+AC_OUTPUT

+

@@ -198,3 +198,41 @@ void AtomicDomain::updateMass(uint64_t pos, float newMass)

 {

     mAtoms[mAtomPositions.at(pos)].mass = newMass; // TODO at is C++11

 }

+

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

+{

+    ar << a.leftNdx << a.rightNdx << a.pos << a.mass;

+    return ar;

+}

+

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

+{

+    ar >> a.leftNdx >> a.rightNdx >> a.pos >> a.mass;

+    return ar;

+}

+

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

+{

+    unsigned nAtoms = domain.mAtoms.size();

+    ar << nAtoms << domain.mDomainSize;

+

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

+    {

+        ar << domain.mAtoms[i];

+    }

+    return ar;

+}

+

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

+{

+    unsigned nAtoms = 0;

+    ar >> nAtoms >> domain.mDomainSize;

+

+    Atom a(0, 0.f);

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

+    {

+        ar >> a;

+        domain.insert(a.pos, a.mass);

+    }

+    return ar;

+}

@@ -17,8 +17,7 @@ class AtomicDomain;

 // first index. AtomicDomain is responsible for managing this correctly.

 struct Atom

 {

-//private:    

-public: // TODO

+private:    

     friend class AtomicDomain;

     uint64_t leftNdx; // shift these up 1, 0 == no neighbor

@@ -42,6 +41,10 @@ public:

     {

         return !(pos == other.pos);

     }

+

+    // serialization

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

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

 };

 struct RawAtom

@@ -10,11 +10,9 @@

 // create "nSets" vectors where each vector contains a vector of indices in the

 // range [0,n)

-// see createGWCoGAPSSets.R - here we just sample indices, that function

-// samples gene names as well

 static std::vector< std::vector<unsigned> > sampleIndices(unsigned n, unsigned nSets)

 {

-    unsigned setSize = (int) n / nSets;

+    unsigned setSize = n / nSets;

     std::vector< std::vector<unsigned> > sampleIndices;

     std::vector<unsigned> toBeSampled;

     for (unsigned i = 1; i < n; ++i)

@@ -32,51 +30,6 @@ static std::vector< std::vector<unsigned> > sampleIndices(unsigned n, unsigned n

     sampleIndices.push_back(toBeSampled);

     return sampleIndices;

-

-    /*

-    std::vector< std::vector<unsigned> > sampleIndices;

-    std::vector<unsigned> sampled;

-

-    for (unsigned i = 0; i < (n - 1) % nSets; ++i)

-    {

-        std::vector<unsigned> set;

-        for (unsigned i = 0; i < ((unsigned) (n - 1) / nSets) + 1; ++i)

-        {

-            while(true)

-            {

-                unsigned sample = gaps::random::uniform64(1, n);

-                if (find(sampled.begin(), sampled.end(), sample) == sampled.end())

-                {

-                    set.push_back(sample);

-                    sampled.push_back(sample);

-                    break;

-                }

-            }

-        }

-        sampleIndices.push_back(set);

-    }

-

-    for (unsigned i = (n - 1) % nSets; i < nSets; ++i)

-    {

-        std::vector<unsigned> set;

-        for (unsigned i = 0; i < (unsigned) (n - 1) / nSets; ++i)

-        {

-            while(true)

-            {

-                unsigned sample = gaps::random::uniform64(1, n);

-                if (find(sampled.begin(), sampled.end(), sample) == sampled.end())

-                {

-                    set.push_back(sample);

-                    sampled.push_back(sample);

-                    break;

-                }

-            }

-        }

-        sampleIndices.push_back(set);

-    }

-

-    return sampleIndices;

-    */

 }

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

@@ -89,7 +42,7 @@ char whichMatrixFixed, const Rcpp::NumericMatrix &FP, unsigned nCores)

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

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

 mIterA(10), mIterP(10), mEquilIter(nEquil), mCoolIter(nCool),

-mSampleIter(nSample), mNumPatterns(nFactor), mNumOutputs(nOutputs),

+mSampleIter(nSample), mNumOutputs(nOutputs),

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

 mCheckpointInterval(cptInterval), mCheckpointFile(cptFile),

 mNumUpdatesA(0), mNumUpdatesP(0),

@@ -111,18 +64,19 @@ mChiSqSample(nSample), mNumAAtomsSample(nSample), mNumPAtomsSample(nSample),

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

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

 {

-    //Archive ar(cptFile, ARCHIVE_READ);

-    //gaps::random::load(ar);

-

-   //ar >> mChiSqEquil >> mNumAAtomsEquil >> mNumPAtomsEquil >> mChiSqSample

-   //    >> mNumAAtomsSample >> mNumPAtomsSample >> mIterA >> mIterP

-   //    >> mEquilIter >> mCoolIter >> mSampleIter >> mNumPatterns >> mNumOutputs

-   //    >> mPrintMessages >> mCurrentIter >> mPhase >> mSeed

-   //    >> mCheckpointInterval >> mCheckpointFile >> mNumUpdatesA

-   //    >> mNumUpdatesP >> mASampler >> mPSampler >> mStatistics;

-

-    //mASampler.sync(mPSampler);

-    //mPSampler.sync(mASampler);

+    Archive ar(cptFile, ARCHIVE_READ);

+    gaps::random::load(ar);

+

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

+

+    mASampler.sync(mPSampler);

+    mPSampler.sync(mASampler);

 }

 // execute the steps of the algorithm, return list to R

@@ -224,7 +178,6 @@ void GapsRunner::runSampPhase()

     }

 }

-

 // sum coef * log(i) for i = 1 to total, fit coef from number of atoms

 // approximates sum of number of atoms

 // this should be proportional to total number of updates

@@ -316,29 +269,30 @@ void GapsRunner::displayStatus(const std::string &type, unsigned nIterTotal)

 void GapsRunner::createCheckpoint()

 {

     // create backup file

-    //std::rename(mCheckpointFile.c_str(), (mCheckpointFile + ".backup").c_str());

+    std::rename(mCheckpointFile.c_str(), (mCheckpointFile + ".backup").c_str());

     // record starting time

-    //bpt::ptime start = bpt_now();

+    bpt::ptime start = bpt_now();

     // save state to file, write magic number at beginning

-    //Archive ar(mCheckpointFile, ARCHIVE_WRITE);

-    //gaps::random::save(ar);

-    //ar << mChiSqEquil << mNumAAtomsEquil << mNumPAtomsEquil << mChiSqSample

-    //    << mNumAAtomsSample << mNumPAtomsSample << mIterA << mIterP

-    //    << mEquilIter << mCoolIter << mSampleIter << mNumPatterns << mNumOutputs

-    //    << mPrintMessages << mCurrentIter << mPhase << mSeed

-    //    << mCheckpointInterval << mNumUpdatesA << mNumUpdatesP << mASampler

-    //    << mPSampler << mStatistics;

-    //ar.close();

+    Archive ar(mCheckpointFile, ARCHIVE_WRITE);

+    gaps::random::save(ar);

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

+    ar.close();

     // display time it took to create checkpoint

-    //bpt::time_duration diff = bpt_now() - start;

-    //double elapsed = diff.total_milliseconds() / 1000.;

-    //Rprintf("created checkpoint in %.3f seconds\n", elapsed);

+    bpt::time_duration diff = bpt_now() - start;

+    double elapsed = diff.total_milliseconds() / 1000.;

+    Rprintf("created checkpoint in %.3f seconds\n", elapsed);

     // delete backup file

-    //std::remove((mCheckpointFile + ".backup").c_str());

+    std::remove((mCheckpointFile + ".backup").c_str());

 }

 void GapsRunner::makeCheckpointIfNeeded()

@@ -46,7 +46,7 @@ public:

     unsigned mCoolIter;

     unsigned mSampleIter;

-    unsigned mNumPatterns;

+    //unsigned mNumPatterns;

     unsigned mNumOutputs;

     bool mPrintMessages;

@@ -4,7 +4,8 @@

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

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

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

-        mNumPatterns(nFactor), mStatUpdates(0)

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

+        mPumpThreshold(PUMP_CUT), mPumpStatUpdates(0)

 {}

 void GapsStatistics::update(const AmplitudeGibbsSampler &ASampler,

@@ -12,6 +13,7 @@ const PatternGibbsSampler &PSampler)

 {

     mStatUpdates++;

+    // update     

     for (unsigned j = 0; j < mNumPatterns; ++j)

     {

         float norm = gaps::algo::sum(PSampler.mMatrix.getRow(j));

@@ -27,6 +29,103 @@ const PatternGibbsSampler &PSampler)

     }

 }

+static unsigned geneThreshold(const ColMatrix &rankMatrix, unsigned pat)

+{

+    float cutRank = rankMatrix.nRow();

+    for (unsigned i = 0; i < rankMatrix.nRow(); ++i)

+    {

+        for (unsigned j = 0; j < rankMatrix.nCol(); ++j)

+        {

+            if (j != pat && rankMatrix(i,j) <= rankMatrix(i,pat))

+            {

+                cutRank = std::min(cutRank, std::max(0.f, rankMatrix(i,pat)-1));

+            }

+        }

+    }

+    return static_cast<unsigned>(cutRank);

+}

+

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

+const PatternGibbsSampler &PSampler)

+{

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

+

+    // norm matrices

+    for (unsigned j = 0; j < nPatterns; ++j)

+    {

+        float factor = gaps::algo::sum(normedP.getRow(j));

+        factor = (factor == 0) ? 1.f : factor;

+

+        normedP.getRow(j) = normedP.getRow(j) / factor;

+        float scale = gaps::algo::max(normedP.getRow(j));

+        normedA.getCol(j) = normedA.getCol(j) * factor * scale;

+    }

+

+    RowMatrix scaledA(nGenes, nPatterns);

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

+    {

+        for (unsigned j = 0; j < nPatterns; ++j)

+        {

+            scaledA(i,j) = normedA(i,j);

+        }

+    }

+   

+    // compute sstat

+    RowMatrix rsStat(nGenes, nPatterns);

+    ColMatrix csStat(nGenes, nPatterns);

+    Vector lp(nPatterns), diff(nPatterns);

+    for (unsigned j = 0; j < nPatterns; ++j)

+    {

+        lp[j] = 1.f;

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

+        {

+            float geneMax = gaps::algo::max(scaledA.getRow(i));

+            diff = geneMax > 0.f ? scaledA.getRow(i) / geneMax - lp : lp * -1.f;

+            rsStat(i,j) = std::sqrt(gaps::algo::dot(diff, diff));

+            csStat(i,j) = rsStat(i,j);

+        }

+        lp[j] = 0.f;

+    }

+

+    // update PUMP matrix

+    if (mPumpThreshold == PUMP_UNIQUE)

+    {

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

+        {

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

+            mPumpMatrix(i,minNdx)++;

+        }

+    }

+    else if (mPumpThreshold == PUMP_CUT)

+    {

+        ColMatrix rankMatrix(nGenes, nPatterns);

+        for (unsigned j = 0; j < nPatterns; ++j)

+        {

+            rankMatrix.getCol(j) = gaps::algo::rank(csStat.getCol(j));

+        }

+        

+        for (unsigned j = 0; j < nPatterns; ++j)

+        {

+            unsigned cutRank = geneThreshold(rankMatrix, j);

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

+            {

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

+                {

+                    mPumpMatrix(i,j)++;

+                }

+            }

+        }

+    }

+}

+

 Rcpp::NumericMatrix GapsStatistics::AMean() const

 {

     return (mAMeanMatrix / mStatUpdates).rMatrix();

@@ -58,3 +157,17 @@ float GapsStatistics::meanChiSq(const AmplitudeGibbsSampler &ASampler) const

         M);

 }

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

+{

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

+        << stat.mPStdMatrix << stat.mStatUpdates << stat.mNumPatterns;

+    return ar;

+}

+

+Archive& operator>>(Archive &ar, GapsStatistics &stat)

+{

+    ar >> stat.mAMeanMatrix >> stat.mAStdMatrix >> stat.mPMeanMatrix

+        >> stat.mPStdMatrix >> stat.mStatUpdates >> stat.mNumPatterns;

+    return ar;

+}

+

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

 #include "GibbsSampler.h"

 #include "data_structures/Matrix.h"

+enum PumpThreshold

+{

+    PUMP_UNIQUE=1,

+    PUMP_CUT=2

+};

+

 class GapsStatistics

 {

 private:

@@ -16,6 +22,10 @@ private:

     unsigned mStatUpdates;

     unsigned mNumPatterns;

+    ColMatrix mPumpMatrix;

+    PumpThreshold mPumpThreshold;

+    unsigned mPumpStatUpdates;

+

 public:

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

@@ -29,6 +39,13 @@ public:

     void update(const AmplitudeGibbsSampler &ASampler,

         const PatternGibbsSampler &PSampler);

+

+    void updatePump(const AmplitudeGibbsSampler &ASampler,

+        const PatternGibbsSampler &PSampler);

+

+    // serialization

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

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

 };

 #endif

\ No newline at end of file

@@ -172,7 +172,6 @@ T* GibbsSampler<T, MatA, MatB>::impl()

 template <class T, class MatA, class MatB>

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

 {

-    static unsigned count = 0;

     unsigned n = 0;

     while (n < nSteps)

     {

new file mode 100644

@@ -0,0 +1,27 @@

+PKG_CPPFLAGS =  -DBOOST_MATH_PROMOTE_DOUBLE_POLICY=0  

+PKG_CXXFLAGS =  -fopenmp  -msse4.2 

+PKG_LIBS =  -fopenmp 

+

+OBJECTS =   AtomicDomain.o \

+            Cogaps.o \

+            GapsRunner.o \

+            GapsStatistics.o \

+            GibbsSampler.o \

+            ProposalQueue.o \

+            RcppExports.o \

+            test-runner.o \

+            data_structures/Matrix.o \

+            data_structures/Vector.o \

+            file_parser/CsvParser.o \

+            file_parser/TsvParser.o \

+            file_parser/MtxParser.o \

+            math/Algorithms.o \

+            math/Random.o \

+            cpp_tests/testAlgorithms.o \

+            cpp_tests/testEfficientSets.o \

+            cpp_tests/testFileParsers.o \

+            cpp_tests/testGapsRunner.o \

+            cpp_tests/testMatrix.o \

+            cpp_tests/testRandom.o \

+            cpp_tests/testSerialization.o

+

@@ -18,12 +18,10 @@ OBJECTS =   AtomicDomain.o \

             math/Algorithms.o \

             math/Random.o \

             cpp_tests/testAlgorithms.o \

-            cpp_tests/testAtomicDomain.o \

             cpp_tests/testEfficientSets.o \

             cpp_tests/testFileParsers.o \

             cpp_tests/testGapsRunner.o \

-            cpp_tests/testGibbsSampler.o \

             cpp_tests/testMatrix.o \

-            cpp_tests/testProposalQueue.o \

-            cpp_tests/testRandom.o

+            cpp_tests/testRandom.o \

+            cpp_tests/testSerialization.o

@@ -37,7 +37,6 @@ void ProposalQueue::populate(AtomicDomain &domain, unsigned limit)

     }

 }

-// TODO efficiently allow clearing a range of proposals

 void ProposalQueue::clear()

 {

     mQueue.clear();

@@ -220,3 +219,17 @@ bool ProposalQueue::exchange(AtomicDomain &domain)

     --mMinAtoms;

     return true;

 }

+

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

+{

+    ar << queue.mMinAtoms << queue.mMaxAtoms << queue.mNumBins

+        << queue.mDimensionSize << queue.mDomainSize << queue.mAlpha;

+    return ar;

+}

+

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

+{

+    ar >> queue.mMinAtoms >> queue.mMaxAtoms >> queue.mNumBins

+        >> queue.mDimensionSize >> queue.mDomainSize >> queue.mAlpha;

+    return ar;

+}

@@ -23,16 +23,6 @@ struct AtomicProposal

     {}

 };

-// Note that mUsedIndices can only contain 1 entry per row/col - otherwise

-// there will be a conflict - one option is to make a static allocation at

-// the start the size of the total rows/cols, make it a vector of bool or char

-// and store 0/1 - clear takes a long time though. could store a vector of used

-// indices as well, use the vector to clear quickly

-

-// could use unique indentifier integer - increment with each clear, only return

-// true if equal to indentifier

-

-// generate single atomic proposal for now

 class ProposalQueue

 {

 private:

@@ -4,25 +4,18 @@

 #define MAT_SUM(nR, nC) ((nR + nC - 2) * nR * nC / 2.f)

-#if 0

-

 TEST_CASE("Test Algorithms.h")

 {

     unsigned nrow = 25;

     unsigned ncol = 20;

     unsigned nfactor = 7;

-    Vector v(nrow);

-    TwoWayMatrix D(nrow, ncol), S(nrow, ncol), AP(nrow, ncol);

     ColMatrix A(nrow, nfactor);

     RowMatrix P(nfactor, ncol);

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

     {

-        v[i] = i;

         for (unsigned j = 0; j < ncol; ++j)

         {

-            D.set(i, j, (i + j) / 10.f);