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

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

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

-cogaps <- function(DMatrix, SMatrix, nFactor, alphaA, alphaP, nEquil, nEquilCool, nSample, maxGibbsMassA, maxGibbsMassP, fixedPatterns, whichMatrixFixed, seed, messages, singleCellRNASeq) {

-    .Call('_CoGAPS_cogaps', PACKAGE = 'CoGAPS', DMatrix, SMatrix, nFactor, alphaA, alphaP, nEquil, nEquilCool, nSample, maxGibbsMassA, maxGibbsMassP, fixedPatterns, whichMatrixFixed, seed, messages, singleCellRNASeq)

+cogaps <- function(DMatrix, SMatrix, nFactor, alphaA, alphaP, nEquil, nEquilCool, nSample, maxGibbsMassA, maxGibbsMassP, fixedPatterns, whichMatrixFixed, seed, messages, singleCellRNASeq, numOutputs, numSnapshots) {

+    .Call('_CoGAPS_cogaps', PACKAGE = 'CoGAPS', DMatrix, SMatrix, nFactor, alphaA, alphaP, nEquil, nEquilCool, nSample, maxGibbsMassA, maxGibbsMassP, fixedPatterns, whichMatrixFixed, seed, messages, singleCellRNASeq, numOutputs, numSnapshots)

 }

 run_catch_unit_tests <- function() {

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

 #' given in Abins and Pbins

 #'@param fixedDomain character to indicate whether A or P is

 #' domain for relative probabilities

-#'@param sampleSnapshots Boolean to indicate whether to capture

-#' individual samples from Markov chain during sampling

 #'@param numSnapshots the number of individual samples to capture

 #'@param alphaA sparsity parameter for A domain

 #'@param nMaxA PRESENTLY UNUSED, future = limit number of atoms

@@ -61,8 +59,7 @@ gapsRun <- function(D, S, ABins = data.frame(), PBins = data.frame(),

                     nFactor = 7, simulation_id = "simulation",

                     nEquil = 1000, nSample = 1000, nOutR = 1000,

                     output_atomic = FALSE, fixedBinProbs = FALSE,

-                    fixedDomain = "N", sampleSnapshots = TRUE,

-                    numSnapshots = 100, alphaA = 0.01,

+                    fixedDomain = "N", numSnapshots = 0, alphaA = 0.01,

                     nMaxA = 100000, max_gibbmass_paraA = 100.0,

                     alphaP = 0.01, nMaxP = 100000, max_gibbmass_paraP = 100.0,

                     seed=-1, messages=TRUE, singleCellRNASeq=FALSE,

@@ -102,7 +99,8 @@ gapsRun <- function(D, S, ABins = data.frame(), PBins = data.frame(),

     # call to C++ Rcpp code

     cogapResult <- cogaps(as.matrix(D), as.matrix(S), nFactor, alphaA, alphaP,

         nEquil, floor(nEquil/10), nSample, max_gibbmass_paraA, max_gibbmass_paraP,

-        fixedPatterns, whichMatrixFixed, seed, messages, singleCellRNASeq)

+        fixedPatterns, whichMatrixFixed, seed, messages, singleCellRNASeq,

+        nOutR, numSnapshots)

     # convert returned files to matrices to simplify visualization and processing

     cogapResult$Amean <- as.matrix(cogapResult$Amean);

@@ -9,9 +9,9 @@ the data matrix}

 gapsRun(D, S, ABins = data.frame(), PBins = data.frame(), nFactor = 7,

   simulation_id = "simulation", nEquil = 1000, nSample = 1000,

   nOutR = 1000, output_atomic = FALSE, fixedBinProbs = FALSE,

-  fixedDomain = "N", sampleSnapshots = TRUE, numSnapshots = 100,

-  alphaA = 0.01, nMaxA = 1e+05, max_gibbmass_paraA = 100, alphaP = 0.01,

-  nMaxP = 1e+05, max_gibbmass_paraP = 100, seed = -1, messages = TRUE,

+  fixedDomain = "N", numSnapshots = 0, alphaA = 0.01, nMaxA = 1e+05,

+  max_gibbmass_paraA = 100, alphaP = 0.01, nMaxP = 1e+05,

+  max_gibbmass_paraP = 100, seed = -1, messages = TRUE,

   singleCellRNASeq = FALSE, fixedPatterns = matrix(0),

   whichMatrixFixed = "N")

 }

@@ -45,9 +45,6 @@ given in Abins and Pbins}

 \item{fixedDomain}{character to indicate whether A or P is

 domain for relative probabilities}

-\item{sampleSnapshots}{Boolean to indicate whether to capture

-individual samples from Markov chain during sampling}

-

 \item{numSnapshots}{the number of individual samples to capture}

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

@@ -138,6 +138,16 @@ Vector gaps::algo::squaredScalarMultiple(const Vector &vec, double n)

     return retVec;

 }

+Vector gaps::algo::squaredScalarDivision(const Vector &vec, double n)

+{

+    Vector retVec(vec.size());

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

+    {

+        retVec(i) = GAPS_SQ(vec(i) / n);

+    }

+    return retVec;

+}

+

 ColMatrix gaps::algo::scalarDivision(const ColMatrix &mat, double n)

 {

     ColMatrix retMat(mat.nRow(), mat.nCol());

@@ -11,16 +11,21 @@ enum GapsPhase

     GAPS_SAMPLING

 };

+static std::vector<Rcpp::NumericMatrix> snapshotsA;

+static std::vector<Rcpp::NumericMatrix> snapshotsP;

+

 static void runGibbsSampler(GibbsSampler &sampler, unsigned nIterTotal,

 unsigned& nIterA, unsigned& nIterP, Vector& chi2Vec, Vector& aAtomVec,

-Vector& pAtomVec, GapsPhase phase);

+Vector& pAtomVec, GapsPhase phase, unsigned numOutputs, bool messages,

+unsigned numSnapshots);

 // [[Rcpp::export]]

 Rcpp::List cogaps(Rcpp::NumericMatrix DMatrix, Rcpp::NumericMatrix SMatrix,

 unsigned nFactor, double alphaA, double alphaP, unsigned nEquil,

 unsigned nEquilCool, unsigned nSample, double maxGibbsMassA,

 double maxGibbsMassP, Rcpp::NumericMatrix fixedPatterns,

-char whichMatrixFixed, int seed, bool messages, bool singleCellRNASeq)

+char whichMatrixFixed, int seed, bool messages, bool singleCellRNASeq,

+unsigned numOutputs, unsigned numSnapshots)

 {

     // set seed

     uint32_t seedUsed = static_cast<uint32_t>(seed);

@@ -47,33 +52,32 @@ char whichMatrixFixed, int seed, bool messages, bool singleCellRNASeq)

     Vector nAtomsAEquil(nEquil);

     Vector nAtomsPEquil(nEquil);

     runGibbsSampler(sampler, nEquil, nIterA, nIterP, chi2Vec,

-        nAtomsAEquil, nAtomsPEquil, GAPS_CALIBRATION);

+        nAtomsAEquil, nAtomsPEquil, GAPS_CALIBRATION, numOutputs, messages,

+        numSnapshots);

     // cooling phase

     Vector trash(nEquilCool);

     runGibbsSampler(sampler, nEquilCool, nIterA, nIterP, trash,

-        trash, trash, GAPS_COOLING);

+        trash, trash, GAPS_COOLING, numOutputs, messages, numSnapshots);

     // sampling phase

     Vector chi2VecSample(nSample);

     Vector nAtomsASample(nSample);

     Vector nAtomsPSample(nSample);

     runGibbsSampler(sampler, nSample, nIterA, nIterP, chi2VecSample,

-        nAtomsASample, nAtomsPSample, GAPS_SAMPLING);

+        nAtomsASample, nAtomsPSample, GAPS_SAMPLING, numOutputs, messages,

+        numSnapshots);

     // combine chi2 vectors

     chi2Vec.concat(chi2VecSample);

-    //Just leave the snapshots as empty lists

-    Rcpp::List ASnapR = Rcpp::List::create();

-    Rcpp::List PSnapR = Rcpp::List::create();

     return Rcpp::List::create(

         Rcpp::Named("Amean") = sampler.AMeanRMatrix(),

         Rcpp::Named("Asd") = sampler.AStdRMatrix(),

         Rcpp::Named("Pmean") = sampler.PMeanRMatrix(),

         Rcpp::Named("Psd") = sampler.PStdRMatrix(),

-        Rcpp::Named("ASnapshots") = ASnapR,

-        Rcpp::Named("PSnapshots") = PSnapR,

+        Rcpp::Named("ASnapshots") = Rcpp::wrap(snapshotsA),

+        Rcpp::Named("PSnapshots") = Rcpp::wrap(snapshotsP),

         Rcpp::Named("atomsAEquil") = nAtomsAEquil.rVec(),

         Rcpp::Named("atomsASamp") = nAtomsASample.rVec(),

         Rcpp::Named("atomsPEquil") = nAtomsPEquil.rVec(),

@@ -85,7 +89,8 @@ char whichMatrixFixed, int seed, bool messages, bool singleCellRNASeq)

 static void runGibbsSampler(GibbsSampler &sampler, unsigned nIterTotal,

 unsigned& nIterA, unsigned& nIterP, Vector& chi2Vec, Vector& aAtomVec,

-Vector& pAtomVec, GapsPhase phase)

+Vector& pAtomVec, GapsPhase phase, unsigned numOutputs, bool messages,

+unsigned numSnapshots)

 {

     double tempDenom = (double)nIterTotal / 2.0;

@@ -109,16 +114,29 @@ Vector& pAtomVec, GapsPhase phase)

         if (phase == GAPS_SAMPLING)

         {

             sampler.updateStatistics();

+            if (numSnapshots > 0 && (i + 1) % (nIterTotal / numSnapshots) == 0)

+            {

+                snapshotsA.push_back(sampler.getNormedMatrix('A'));

+                snapshotsP.push_back(sampler.getNormedMatrix('P'));

+            }

         }

-        chi2Vec(i) = sampler.chi2();

         double numAtomsA = sampler.totalNumAtoms('A');

         double numAtomsP = sampler.totalNumAtoms('P');

         aAtomVec(i) = numAtomsA;

         pAtomVec(i) = numAtomsP;

+        chi2Vec(i) = sampler.chi2();

         if (phase != GAPS_COOLING)

         {

+            if ((i + 1) % numOutputs == 0 && messages)

+            {

+                std::string temp = phase == GAPS_CALIBRATION ? "Equil: " : "Samp: ";

+                std::cout << temp << i + 1 << " of " << nIterTotal << ", Atoms:"

+                    << numAtomsA << "(" << numAtomsP << ") Chi2 = "

+                    << sampler.chi2() << '\n';

+            }

+

             nIterA = gaps::random::poisson(std::max(numAtomsA, 10.0));

             nIterP = gaps::random::poisson(std::max(numAtomsP, 10.0));

         }

@@ -54,6 +54,40 @@ mStatUpdates(0), mFixedMat(whichMat)

     mChi2 = 2.0 * gaps::algo::loglikelihood(mDMatrix, mSMatrix, mAPMatrix);

 }

+Rcpp::NumericMatrix GibbsSampler::getNormedMatrix(char mat)

+{

+    Vector normVec(mPMatrix.nRow());

+    for (unsigned r = 0; r < mPMatrix.nRow(); ++r)

+    {

+        normVec(r) = gaps::algo::sum(mPMatrix.getRow(r));

+        if (normVec(r) == 0)

+        {

+            normVec(r) = 1.0;

+        }

+    }

+

+    if (mat == 'A')

+    {

+        ColMatrix normedA(mAMatrix);    

+        for (unsigned c = 0; c < normedA.nCol(); ++c)

+        {

+            normedA.getCol(c) += gaps::algo::scalarMultiple(normedA.getCol(c),

+                normVec(c));

+        }

+        return normedA.rMatrix();

+    }

+    else

+    {

+        RowMatrix normedP(mPMatrix);

+        for (unsigned r = 0; r < normedP.nRow(); ++r)

+        {

+            normedP.getRow(r) += gaps::algo::scalarDivision(normedP.getRow(r),

+                normVec(r));

+        }

+        return normedP.rMatrix();

+    }

+}

+

 double GibbsSampler::getGibbsMass(const MatrixChange &change)

 {

     // check if this change is death (only called in birth/death)

@@ -396,10 +430,10 @@ void GibbsSampler::updateStatistics()

     for (unsigned r = 0; r < mPMeanMatrix.nRow(); ++r)

     {

-        mPMeanMatrix.getRow(r) += gaps::algo::scalarMultiple(mPMatrix.getRow(r),

-            1.0 / normVec(r));

+        mPMeanMatrix.getRow(r) += gaps::algo::scalarDivision(mPMatrix.getRow(r),

+            normVec(r));

-        mPStdMatrix.getRow(r) += gaps::algo::squaredScalarMultiple(mPMatrix.getRow(r),

-            1.0 / normVec(r));

+        mPStdMatrix.getRow(r) += gaps::algo::squaredScalarDivision(mPMatrix.getRow(r),

+            normVec(r));

     }

 }

\ No newline at end of file

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

     Rcpp::NumericMatrix PStdRMatrix() const;

     void updateStatistics();

+

+    Rcpp::NumericMatrix getNormedMatrix(char mat);

 };

 #endif

@@ -51,6 +51,15 @@ RowMatrix::RowMatrix(const Rcpp::NumericMatrix &rmat)

     }

 }

+RowMatrix::RowMatrix(const RowMatrix &mat)

+: mNumRows(mat.nRow()), mNumCols(mat.nCol())

+{

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

+    {

+        mRows.push_back(mat.getRow(i));

+    }

+}

+

 Vector& RowMatrix::getRow(unsigned row)

 {

     return mRows[row];

@@ -116,6 +125,15 @@ ColMatrix::ColMatrix(const Rcpp::NumericMatrix &rmat)

     }

 }

+ColMatrix::ColMatrix(const ColMatrix &mat)

+: mNumRows(mat.nRow()), mNumCols(mat.nCol())

+{

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

+    {

+        mCols.push_back(mat.getCol(j));

+    }

+}

+

 Vector& ColMatrix::getCol(unsigned col)

 {

     return mCols[col];

@@ -74,6 +74,7 @@ public:

     RowMatrix(unsigned nrow, unsigned ncol);

     RowMatrix(const Rcpp::NumericMatrix& rmat);

+    RowMatrix(const RowMatrix &mat);

     unsigned nRow() const {return mNumRows;}

     unsigned nCol() const {return mNumCols;}

@@ -99,6 +100,7 @@ public:

     ColMatrix(unsigned nrow, unsigned ncol);

     ColMatrix(const Rcpp::NumericMatrix& rmat);

+    ColMatrix(const ColMatrix &mat);

     unsigned nRow() const {return mNumRows;}

     unsigned nCol() const {return mNumCols;}

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

 using namespace Rcpp;

 // cogaps

-Rcpp::List cogaps(Rcpp::NumericMatrix DMatrix, Rcpp::NumericMatrix SMatrix, unsigned nFactor, double alphaA, double alphaP, unsigned nEquil, unsigned nEquilCool, unsigned nSample, double maxGibbsMassA, double maxGibbsMassP, Rcpp::NumericMatrix fixedPatterns, char whichMatrixFixed, int seed, bool messages, bool singleCellRNASeq);

-RcppExport SEXP _CoGAPS_cogaps(SEXP DMatrixSEXP, SEXP SMatrixSEXP, SEXP nFactorSEXP, SEXP alphaASEXP, SEXP alphaPSEXP, SEXP nEquilSEXP, SEXP nEquilCoolSEXP, SEXP nSampleSEXP, SEXP maxGibbsMassASEXP, SEXP maxGibbsMassPSEXP, SEXP fixedPatternsSEXP, SEXP whichMatrixFixedSEXP, SEXP seedSEXP, SEXP messagesSEXP, SEXP singleCellRNASeqSEXP) {

+Rcpp::List cogaps(Rcpp::NumericMatrix DMatrix, Rcpp::NumericMatrix SMatrix, unsigned nFactor, double alphaA, double alphaP, unsigned nEquil, unsigned nEquilCool, unsigned nSample, double maxGibbsMassA, double maxGibbsMassP, Rcpp::NumericMatrix fixedPatterns, char whichMatrixFixed, int seed, bool messages, bool singleCellRNASeq, unsigned numOutputs, unsigned numSnapshots);

+RcppExport SEXP _CoGAPS_cogaps(SEXP DMatrixSEXP, SEXP SMatrixSEXP, SEXP nFactorSEXP, SEXP alphaASEXP, SEXP alphaPSEXP, SEXP nEquilSEXP, SEXP nEquilCoolSEXP, SEXP nSampleSEXP, SEXP maxGibbsMassASEXP, SEXP maxGibbsMassPSEXP, SEXP fixedPatternsSEXP, SEXP whichMatrixFixedSEXP, SEXP seedSEXP, SEXP messagesSEXP, SEXP singleCellRNASeqSEXP, SEXP numOutputsSEXP, SEXP numSnapshotsSEXP) {

 BEGIN_RCPP

     Rcpp::RObject rcpp_result_gen;

     Rcpp::RNGScope rcpp_rngScope_gen;

@@ -26,7 +26,9 @@ BEGIN_RCPP

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

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

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

-    rcpp_result_gen = Rcpp::wrap(cogaps(DMatrix, SMatrix, nFactor, alphaA, alphaP, nEquil, nEquilCool, nSample, maxGibbsMassA, maxGibbsMassP, fixedPatterns, whichMatrixFixed, seed, messages, singleCellRNASeq));

+    Rcpp::traits::input_parameter< unsigned >::type numOutputs(numOutputsSEXP);

+    Rcpp::traits::input_parameter< unsigned >::type numSnapshots(numSnapshotsSEXP);

+    rcpp_result_gen = Rcpp::wrap(cogaps(DMatrix, SMatrix, nFactor, alphaA, alphaP, nEquil, nEquilCool, nSample, maxGibbsMassA, maxGibbsMassP, fixedPatterns, whichMatrixFixed, seed, messages, singleCellRNASeq, numOutputs, numSnapshots));

     return rcpp_result_gen;

 END_RCPP

 }

@@ -42,7 +44,7 @@ END_RCPP

 }

 static const R_CallMethodDef CallEntries[] = {

-    {"_CoGAPS_cogaps", (DL_FUNC) &_CoGAPS_cogaps, 15},

+    {"_CoGAPS_cogaps", (DL_FUNC) &_CoGAPS_cogaps, 17},

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

     {NULL, NULL, 0}

 };