@@ -34,8 +34,8 @@

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

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

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

-fixedPatterns=matrix(0), checkpointInterval=0, checkpointFile="gaps_checkpoint.out",

-...)

+fixedPatterns=matrix(0), checkpointInterval=0, 

+checkpointFile="gaps_checkpoint.out", pumpThreshold="unique", ...)

 {

     # get v2 arguments

     oldArgs <- list(...)

@@ -69,12 +69,13 @@ fixedPatterns=matrix(0), checkpointInterval=0, checkpointFile="gaps_checkpoint.o

         stop('invalid number of rows for fixedPatterns')

     if (whichMatrixFixed == 'P' & ncol(fixedPatterns) != ncol(D))

         stop('invalid number of columns for fixedPatterns')

+    thresholdEnum <- c("unique", "cut")

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

         singleCellRNASeq, whichMatrixFixed, fixedPatterns, checkpointInterval,

-        checkpointFile)

+        checkpointFile, which(thresholdEnum==pumpThreshold))

     # label matrices and return list

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

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

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

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

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

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

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

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

 }

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

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

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

   singleCellRNASeq = FALSE, whichMatrixFixed = "N",

   fixedPatterns = matrix(0), checkpointInterval = 0,

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

+  checkpointFile = "gaps_checkpoint.out", pumpThreshold = "unique", ...)

 }

 \arguments{

 \item{D}{data matrix}

@@ -14,44 +14,75 @@ float gaps::algo::sum(const Vector &vec)

     return sum;

 }

-Vector gaps::algo::scalarMultiple(const Vector &vec, float n)

+float gaps::algo::min(const Vector &vec)

 {

-    Vector retVec(vec.size());

+    float min = vec[0];

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

     {

-        retVec[i] = vec[i] * n;

+        min = std::min(min, vec[i]);

     }

-    return retVec;

+    return min;

 }

-Vector gaps::algo::squaredScalarMultiple(const Vector &vec, float n)

+float gaps::algo::max(const Vector &vec)

 {

-    Vector retVec(vec.size());

+    float max = vec[0];

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

     {

-        retVec[i] = GAPS_SQ(vec[i] * n);

+        max = std::max(max, vec[i]);

     }

-    return retVec;

+    return max;

 }

-Vector gaps::algo::scalarDivision(const Vector &vec, float n)

+float gaps::algo::dot(const Vector &A, const Vector &B)

 {

-    Vector temp(vec.size());

+    float dotProd = 0.f;

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

+    {

+        dotProd += A[i] * B[i];

+    }

+    return dotProd;

+}

+

+unsigned gaps::algo::whichMin(const Vector &vec)

+{

+    float min = vec[0];

+    unsigned minNdx = 0;

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

+    {

+        if (vec[i] < min)

+        {

+            min = vec[i];

+            minNdx = i;

+        }

+    }

+    return minNdx;

+}

+

+Vector gaps::algo::rank(Vector vec)

+{

+    std::vector< std::pair<float, float> > sortVec(vec.size());

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

+    {

+        sortVec[i] = std::pair<float, float>(vec[i], i);

+    }

+    

+    std::sort(sortVec.begin(), sortVec.end());

+    Vector ranks(vec.size());

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

     {

-        temp[i] = vec[i] / n;

+        ranks[i] = sortVec[i].second;

     }

-    return temp;

+    return ranks;

 }

-Vector gaps::algo::squaredScalarDivision(const Vector &vec, float n)

+Vector gaps::algo::elementSq(Vector vec)

 {

-    Vector retVec(vec.size());

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

     {

-        retVec[i] = GAPS_SQ(vec[i] / n);

+        vec[i] *= vec[i];

     }

-    return retVec;

+    return vec;

 }

 bool gaps::algo::isRowZero(const RowMatrix &mat, unsigned row)

@@ -27,11 +27,13 @@ namespace gaps

 namespace algo

 {

     // vector algorithms    

+    unsigned whichMin(const Vector &vec);

     float sum(const Vector &vec);

-    Vector scalarMultiple(const Vector &vec, float n);

-    Vector squaredScalarMultiple(const Vector &vec, float n);

-    Vector scalarDivision(const Vector &vec, float n);

-    Vector squaredScalarDivision(const Vector &vec, float n);

+    float min(const Vector &vec);

+    float max(const Vector &vec);

+    float dot(const Vector &A, const Vector &B);

+    Vector rank(Vector vec);

+    Vector elementSq(Vector vec);

     // generic matrix algorithms

     template<class GenericMatrix>

@@ -43,9 +45,6 @@ namespace algo

     template<class GenericMatrix>

     float nonZeroMean(const GenericMatrix &mat);

-    template<class GenericMatrix>

-    GenericMatrix scalarDivision(const GenericMatrix &mat, float n);

-

     template<class GenericMatrix>

     GenericMatrix computeStdDev(const GenericMatrix &stdMat,

         const GenericMatrix &meanMat, unsigned nUpdates);

@@ -56,7 +55,7 @@ namespace algo

     void matrixMultiplication(TwoWayMatrix &C, const ColMatrix &A,

         const RowMatrix &B);

-    // chi2

+    // chiSq / 2

     float loglikelihood(const TwoWayMatrix &D, const TwoWayMatrix &S,

         const TwoWayMatrix &AP);

@@ -111,26 +110,12 @@ float gaps::algo::nonZeroMean(const GenericMatrix &mat)

     return sum / (float)nNonZeros;

 }

-template<class GenericMatrix>

-GenericMatrix gaps::algo::scalarDivision(const GenericMatrix &mat, float n)

-{

-    GenericMatrix retMat(mat.nRow(), mat.nCol());

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

-    {

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

-        {

-            retMat(r,c) = mat(r,c) / n;

-        }

-    }

-    return retMat;

-}

-

 template<class GenericMatrix>

 GenericMatrix gaps::algo::computeStdDev(const GenericMatrix &stdMat,

 const GenericMatrix &meanMat, unsigned nUpdates)

 {

     GenericMatrix retMat(stdMat.nRow(), stdMat.nCol());

-    float meanTerm = 0.0;

+    float meanTerm = 0.f;

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

     {

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

@@ -129,14 +129,26 @@ AtomicProposal AtomicSupport::proposeExchange() const

     float mass2 = at(pos2);

     // calculate new mass

-    float pupper = gaps::random::p_gamma(mass1 + mass2, 2.0, 1.0 / mLambda);

-    float newMass = gaps::random::q_gamma(gaps::random::uniform(0.0, pupper),

-        2.0, 1.0 / mLambda);

+    float pupper = gaps::random::p_gamma(mass1 + mass2, 2.f, 1.f / mLambda);

+    float newMass = gaps::random::q_gamma(gaps::random::uniform(0.f, pupper),

+        2.f, 1.f / mLambda);

     // calculate mass changes

     float delta1 = mass1 > mass2 ? newMass - mass1 : mass2 - newMass;

     float delta2 = mass2 > mass1 ? newMass - mass2 : mass1 - newMass;

+    // set first position to be postive change

+    float delta1_cached = delta1;

+    uint64_t pos1_cached = pos1;

+    pos1 = delta1_cached > 0 ? pos1 : pos2;

+    pos2 = delta1_cached > 0 ? pos2 : pos1_cached;

+    delta1 = delta1_cached > 0 ? delta1 : delta2;

+    delta2 = delta1_cached > 0 ? delta2 : delta1_cached;

+    //delta1 = at(pos1) + (delta1 > 0) ? delta1 : delta2;

+    //delta2 = at(pos2) + (delta1 > 0) ? delta2 : delta1_temp;

+    //delta1 -= at(pos1);

+    //delta2 -= at(pos2);

+

     // preserve total mass

     return AtomicProposal(mLabel, 'E', pos1, delta1, pos2, delta2);

 }

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

     void removeAtom(uint64_t pos); // O(logN)

     AtomNeighbors getNeighbors(uint64_t pos) const; // O(logN)

-    // convert atomic position to row/col of the matrix

-    uint64_t getRow(uint64_t pos) const;

-    uint64_t getCol(uint64_t pos) const;

-

     // get atomic positions at random

     uint64_t randomFreePosition() const;

     uint64_t randomAtomPosition() const;

@@ -115,6 +111,10 @@ public:

     // convert an AtomicProposal to an ElementChange

     MatrixChange getMatrixChange(const AtomicProposal &prop) const;

+    // convert atomic position to row/col of the matrix

+    uint64_t getRow(uint64_t pos) const;

+    uint64_t getCol(uint64_t pos) const;

+

     // getters

     float alpha() const {return mAlpha;}

     float lambda() const {return mLambda;}

@@ -100,8 +100,8 @@ static void takeSnapshots(GapsInternalState &state)

 {

     if (state.nSnapshots && !((state.iter+1)%(state.nSample/state.nSnapshots)))

     {

-        state.snapshotsA.push_back(state.sampler.getNormedMatrix('A'));

-        state.snapshotsP.push_back(state.sampler.getNormedMatrix('P'));

+        state.snapshotsA.push_back(state.sampler.normedAMatrix().rMatrix());

+        state.snapshotsP.push_back(state.sampler.normedPMatrix().rMatrix());

     }    

 }

@@ -135,9 +135,9 @@ static void runSampPhase(GapsInternalState &state)

         makeCheckpointIfNeeded(state);

         updateSampler(state);

         state.sampler.updateStatistics();

+        state.sampler.updatePumpStatistics();

         takeSnapshots(state);

         displayStatus(state, "Samp: ", state.nSample);

-        state.sampler.updatePumpStatistics();

         storeSamplerInfo(state, state.nAtomsASample, state.nAtomsPSample,

             state.chi2VecSample);

     }

@@ -193,8 +193,8 @@ static Rcpp::List runCogaps(GapsInternalState &state)

         Rcpp::Named("randSeed") = state.seed,

         Rcpp::Named("numUpdates") = state.nUpdatesA + state.nUpdatesP,

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

-        Rcpp::Named("PumpStats") = state.sampler.getPumpMatrix(),

-        Rcpp::Named("meanPatternAssignment") = state.sampler.getMeanPattern()

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

+        Rcpp::Named("meanPatternAssignment") = state.sampler.meanPattern()

     );

 }

@@ -205,7 +205,7 @@ unsigned nEquilCool, unsigned nSample, unsigned nOutputs, unsigned nSnapshots,

 float alphaA, float alphaP, float maxGibbmassA, float maxGibbmassP, int seed,

 bool messages, bool singleCellRNASeq, char whichMatrixFixed,

 const Rcpp::NumericMatrix &FP, unsigned checkpointInterval,

-const std::string &cptFile)

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

 {

     // set seed

     uint32_t seedUsed = static_cast<uint32_t>(seed);

@@ -220,7 +220,8 @@ const std::string &cptFile)

     // create internal state from parameters and run from there

     GapsInternalState state(D, S, nFactor, nEquil, nEquilCool, nSample,

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

-        messages, singleCellRNASeq, whichMatrixFixed, FP, checkpointInterval);

+        messages, singleCellRNASeq, whichMatrixFixed, FP, checkpointInterval,

+        static_cast<PumpThreshold>(pumpThreshold));

     checkpointFile = cptFile;

     return runCogaps(state);

 }

@@ -12,23 +12,24 @@ mDMatrix(D), mSMatrix(S), mAPMatrix(D.nrow(), D.ncol()),

 mAMatrix(D.nrow(), nFactor), mPMatrix(nFactor, D.ncol()),

 mADomain('A', D.nrow(), nFactor), mPDomain('P', nFactor, D.ncol()),

 mAMeanMatrix(D.nrow(), nFactor), mAStdMatrix(D.nrow(), nFactor),

-mPMeanMatrix(nFactor, D.ncol()), mPStdMatrix(nFactor, D.ncol())

+mPMeanMatrix(nFactor, D.ncol()), mPStdMatrix(nFactor, D.ncol()),

+mPumpMatrix(D.nrow(), nFactor)

 {}

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

 const Rcpp::NumericMatrix &S, unsigned nFactor, float alphaA, float alphaP,

 float maxGibbmassA, float maxGibbmassP, bool singleCellRNASeq,

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

+char whichMatrixFixed, const Rcpp::NumericMatrix &FP, PumpThreshold pumpThreshold)

     :

 mDMatrix(D), mSMatrix(S), mAPMatrix(D.nrow(), D.ncol()),

 mAMatrix(D.nrow(), nFactor), mPMatrix(nFactor, D.ncol()),

 mADomain('A', D.nrow(), nFactor), mPDomain('P', nFactor, D.ncol()),

 mAMeanMatrix(D.nrow(), nFactor), mAStdMatrix(D.nrow(), nFactor),

 mPMeanMatrix(nFactor, D.ncol()), mPStdMatrix(nFactor, D.ncol()),

-mPumpMatrix(D.nrow(), nFactor),

-mStatUpdates(0), mMaxGibbsMassA(maxGibbmassA), mMaxGibbsMassP(maxGibbmassP),

-mAnnealingTemp(1.0), mSingleCellRNASeq(singleCellRNASeq),

-mNumFixedPatterns(0), mFixedMat(whichMatrixFixed)

+mPumpMatrix(D.nrow(), nFactor), mPumpThreshold(pumpThreshold), mStatUpdates(0),

+mPumpStatUpdates(0), mMaxGibbsMassA(maxGibbmassA), mMaxGibbsMassP(maxGibbmassP),

+mAnnealingTemp(1.0), mSingleCellRNASeq(singleCellRNASeq), mNumFixedPatterns(0),

+mFixedMat(whichMatrixFixed)

 {

     float meanD = mSingleCellRNASeq ? gaps::algo::nonZeroMean(mDMatrix)

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

@@ -47,8 +48,7 @@ mNumFixedPatterns(0), mFixedMat(whichMatrixFixed)

         ColMatrix temp(FP);

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

         {

-            mAMatrix.getCol(j) = gaps::algo::scalarDivision(temp.getCol(j),

-                gaps::algo::sum(temp.getCol(j)));

+            mAMatrix.getCol(j) = temp.getCol(j) / gaps::algo::sum(temp.getCol(j));

         }

     }

     else if (mFixedMat == 'P')

@@ -57,44 +57,12 @@ mNumFixedPatterns(0), mFixedMat(whichMatrixFixed)

         RowMatrix temp(FP);

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

         {

-            mPMatrix.getRow(i) = gaps::algo::scalarDivision(temp.getRow(i),

-                gaps::algo::sum(temp.getRow(i)));

+            mPMatrix.getRow(i) = temp.getRow(i) / gaps::algo::sum(temp.getRow(i));

         }

     }

     gaps::algo::matrixMultiplication(mAPMatrix, mAMatrix, mPMatrix);

 }

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

-        normVec[r] = (normVec[r] == 0) ? 1.f : normVec[r];

-    }

-

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

-    }

-}

-

 float GibbsSampler::getGibbsMass(const MatrixChange &change)

 {

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

@@ -309,27 +277,10 @@ bool GibbsSampler::exchange(AtomicSupport &domain, AtomicProposal &proposal)

         return false;

     }

-    uint64_t pos1 = proposal.delta1 > 0 ? proposal.pos1 : proposal.pos2;

-    uint64_t pos2 = proposal.delta1 > 0 ? proposal.pos2 : proposal.pos1;

-

-    float mass1 = domain.at(pos1);

-    float mass2 = domain.at(pos2);

-

-    float newMass1 = mass1 + std::max(proposal.delta1, proposal.delta2);

-    float newMass2 = mass2 + std::min(proposal.delta1, proposal.delta2);

-

-    unsigned row1 = change.delta1 > 0 ? change.row1 : change.row2;

-    unsigned row2 = change.delta1 > 0 ? change.row2 : change.row1;

-

-    unsigned col1 = change.delta1 > 0 ? change.col1 : change.col2;

-    unsigned col2 = change.delta1 > 0 ? change.col2 : change.col1;

-

-    change.row1 = row1;

-    change.col1 = col1;

-    change.row2 = row2;

-    change.col2 = col2;

-    change.delta1 = newMass1 - mass1;

-    change.delta2 = newMass2 - mass2;

+    float mass1 = domain.at(proposal.pos1);

+    float mass2 = domain.at(proposal.pos2);

+    float newMass1 = mass1 + proposal.delta1;

+    float newMass2 = mass2 + proposal.delta2;

     if (canUseGibbs(change))

     {

@@ -352,10 +303,7 @@ bool GibbsSampler::exchange(AtomicSupport &domain, AtomicProposal &proposal)

                 proposal.delta1 = gaps::random::q_norm(u, mean, sd);

                 proposal.delta1 = std::max(proposal.delta1, -mass1);

                 proposal.delta1 = std::min(proposal.delta1, mass2);

-                proposal.pos1 = pos1;

-                proposal.pos2 = pos2;

                 proposal.delta2 = -proposal.delta1;

-

                 return evaluateChange(domain, proposal, 0.0, true);

             }

         }

@@ -367,26 +315,21 @@ bool GibbsSampler::exchange(AtomicSupport &domain, AtomicProposal &proposal)

     float pnew = gaps::random::d_gamma(pnewMass, 2.0, 1.f / domain.lambda());

     float pold = gaps::random::d_gamma(poldMass, 2.0, 1.f / domain.lambda());

-    proposal.pos1 = pos1;

-    proposal.delta1 = newMass1 - mass1;

-    proposal.pos2 = pos2;

-    proposal.delta2 = newMass2 - mass2;

-

-    if (pold == 0.0 && pnew != 0.0)

+    if (pold == 0.f && pnew != 0.f)

     {

-        return evaluateChange(domain, proposal, 0.0, true);

+        return evaluateChange(domain, proposal, 0.f, true);

     }

     else

     {

-        float priorLL = pold == 0.0 ? 0.0 : log(pnew / pold);

-        return evaluateChange(domain, proposal, std::log(gaps::random::uniform())

-            - priorLL);

+        float priorLL = (pold == 0.f) ? 0.f : log(pnew / pold);

+        float rejectProb = std::log(gaps::random::uniform()) - priorLL;

+        return evaluateChange(domain, proposal, rejectProb);

     }

 }

 Rcpp::NumericMatrix GibbsSampler::AMeanRMatrix() const

 {

-    return gaps::algo::scalarDivision(mAMeanMatrix, mStatUpdates).rMatrix();

+    return (mAMeanMatrix / mStatUpdates).rMatrix();

 }

 Rcpp::NumericMatrix GibbsSampler::AStdRMatrix() const

@@ -397,7 +340,7 @@ Rcpp::NumericMatrix GibbsSampler::AStdRMatrix() const

 Rcpp::NumericMatrix GibbsSampler::PMeanRMatrix() const

 {

-    return gaps::algo::scalarDivision(mPMeanMatrix, mStatUpdates).rMatrix();

+    return (mPMeanMatrix / mStatUpdates).rMatrix();

 }

 Rcpp::NumericMatrix GibbsSampler::PStdRMatrix() const

@@ -406,10 +349,24 @@ Rcpp::NumericMatrix GibbsSampler::PStdRMatrix() const

         mStatUpdates).rMatrix();

 }

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

+{

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

+}

+

+Rcpp::NumericMatrix GibbsSampler::meanPattern()

+{

+    ColMatrix Amean(mAMeanMatrix / (float)mStatUpdates);

+    RowMatrix Pmean(mPMeanMatrix / (float)mStatUpdates);

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

+    patternMarkers(Amean, Pmean, mat);

+    return mat.rMatrix();

+}

+

 float GibbsSampler::meanChiSq() const

 {

-    ColMatrix Amean = gaps::algo::scalarDivision(mAMeanMatrix, mStatUpdates);

-    RowMatrix Pmean = gaps::algo::scalarDivision(mPMeanMatrix, mStatUpdates);

+    ColMatrix Amean = mAMeanMatrix / (float)mStatUpdates;

+    RowMatrix Pmean = mPMeanMatrix / (float)mStatUpdates;

     TwoWayMatrix Mmean(Amean.nRow(), Pmean.nCol());

     gaps::algo::matrixMultiplication(Mmean, Amean, Pmean);

     return 2.f * gaps::algo::loglikelihood(mDMatrix, mSMatrix, Mmean);

@@ -418,155 +375,127 @@ float GibbsSampler::meanChiSq() const

 void GibbsSampler::updateStatistics()

 {

     mStatUpdates++;

+    unsigned nPatterns = mAMatrix.nCol();

-    Vector normVec(mPMatrix.nRow());

-        

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

+    Vector normVec(nPatterns);

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

     {

-        normVec[r] = gaps::algo::sum(mPMatrix.getRow(r));

-        normVec[r] = normVec[r] == 0 ? 1.f : normVec[r];

-    }

+        normVec[j] = gaps::algo::sum(mPMatrix.getRow(j));

+        normVec[j] = normVec[j] == 0 ? 1.f : normVec[j];

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

-    {

-        mAMeanMatrix.getCol(c) += gaps::algo::scalarMultiple(mAMatrix.getCol(c),

-            normVec[c]);

-        mAStdMatrix.getCol(c) += gaps::algo::squaredScalarMultiple(mAMatrix.getCol(c),

-            normVec[c]);

-    }

+        Vector quot(mPMatrix.getRow(j) / normVec[j]);

+        mPMeanMatrix.getRow(j) += quot;

+        mPStdMatrix.getRow(j) += gaps::algo::elementSq(quot);

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

-    {

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

-            normVec[r]);

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

-            normVec[r]);

+        Vector prod(mAMatrix.getCol(j) * normVec[j]);

+        mAMeanMatrix.getCol(j) += prod;

+        mAStdMatrix.getCol(j) += gaps::algo::elementSq(prod); 

     }

 }

-ColMatrix normedAMatrix()

+void GibbsSampler::updatePumpStatistics()

+{

+    mPumpStatUpdates++;

+    patternMarkers(normedAMatrix(), normedPMatrix(), mPumpMatrix);

+}

+

+ColMatrix GibbsSampler::normedAMatrix() const

 {

     ColMatrix normedA(mAMatrix);

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

     {

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

-        normedA.getCol(j) += gaps::algo::scalarMultiple(normedA.getCol(j), factor);

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

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

     }

     return normedA;

 }

-RowMatrix normedPMatrix()

+RowMatrix GibbsSampler::normedPMatrix() const

 {

     RowMatrix normedP(mPMatrix);

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

     {

         float factor = gaps::algo::sum(mPMatrix.getRow(i));

-        normedP.getRow(i) += gaps::algo::scalarMultiple(normedP.getRow(i), factor);

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

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

     }

     return normedP;

 }

-void GibbsSampler::updatePumpStatistics(bool uniqueThreshold)

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

 {

-    ColMatrix normedA(normedAMatrix());

-    RowMatrix normedP(normedPMatrix());

-

-    // scale A matrix, assuming P matrix is not scaled

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

-    {

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

-        normedA.getCol(j) = gaps::algo::scalarMultiple(normedA.getCol(j), scale);

-    }

-    

-    // get scaled A matrix

-    RowMatrix diff(normedA.nRow(), normedA.nCol());

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

+    float cutRank = rankMatrix.nRow();

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

     {

-        float rowMax = 0.f;

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

-        {

-            rowMax = std::max(normedA(i,j), rowMax);

-        }

-        

-        if (rowMax > 0.f)

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

         {

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

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

             {

-                diff(i,j) = normedA(i,j) / rowMax - mLP[j];

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

             }

         }

-        else

-        {

-            diff.getRow(i) = gaps::algo::scalarMultiple(mLP, -1.f);

-        }

     }

+    return static_cast<unsigned>(cutRank);

+}

-    // compute sstat

-    RowMatrix sStat(normedA.nRow(), normedA.nCol());

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

-    {

-        sStat.getRow(i) = std::sqrt(gaps::algo::dotProduct(diff.getRow(i), diff.getRow(i)));

-    }

+void GibbsSampler::patternMarkers(RowMatrix normedA, RowMatrix normedP,

+ColMatrix &statMatrix)

+{

+    // helpful notation

+    unsigned nGenes = normedA.nRow();

+    unsigned nPatterns = normedA.nCol();

-    if (uniqueThreshold)

+    // scale A matrix

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

     {

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

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

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

         {

-            unsigned whichMin = 0;

-            float minVal = sStat(i,0);

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

-            {

-                if (sStat(i,j) < minVal)

-                {

-                    minVal = sStat(i,j);

-                    whichMin = j;

-                }

-            }

-            mPumpMatrix(i,whichMin)++;

+            normedA(i,j) *= scale;

         }

     }

-    else

+    

+    // compute sstat

+    TwoWayMatrix sStat(nGenes, nPatterns);

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

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

     {

-        ColMatrix rankMatrix(sStat.nRow(), sStat.nCol());

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

+        lp[j] = 1.f;

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

         {

-            rankMatrix.getCol(j) = gaps::algo::rank(sStat.getCol(j))

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

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

+            sStat.set(i, j, std::sqrt(gaps::algo::dot(diff, diff)));

         }

+        lp[j] = 0.f;

+    }

-        Vector geneThreshold(rankMatrix.nCol());

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

+    // update PUMP matrix

+    if (mPumpThreshold == PUMP_UNIQUE)

+    {

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

         {

-            unsigned maxRank = 0;

-            Vector rankNdx = gaps::algo::order(rankMatrix.getCol(j));

-            for (unsigned i = 0; i < rankNdx; ++i)

-            {

-                unsigned rank = rankMatrix(i,j);

-

-

-

-

-

-

-

-

-        ColMatrix sortedRanks(rankMatrix.nRow(), rankMatrix.nCol());

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

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

+            statMatrix(i,minNdx)++;

+        }

+    }

+    else if (mPumpThreshold == PUMP_CUT)

+    {

+        ColMatrix rankMatrix(nGenes, nPatterns);

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

         {

-            unsigned rank = rankMatrix

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

-            {

-

-            }

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

         }

-

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

+        

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

         {

-            float threshold = geneThreshold(rankMatrix, j);

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

+            unsigned cutRank = geneThreshold(rankMatrix, j);

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

             {

-                if (rankMatrix(i,j) < threshold)

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

                 {

-                    mPumpMatrix(i,j)++;

+                    statMatrix(i,j)++;

                 }

             }

         }

@@ -578,6 +507,7 @@ Archive& operator<<(Archive &ar, GibbsSampler &sampler)

     ar << sampler.mAMatrix << sampler.mPMatrix << sampler.mADomain

         << sampler.mPDomain << sampler.mAMeanMatrix << sampler.mAStdMatrix

         << sampler.mPMeanMatrix << sampler.mPStdMatrix << sampler.mStatUpdates

+        << sampler.mPumpMatrix << sampler.mPumpThreshold << sampler.mPumpStatUpdates

         << sampler.mMaxGibbsMassA << sampler.mMaxGibbsMassP

         << sampler.mAnnealingTemp << sampler.mSingleCellRNASeq

         << sampler.mNumFixedPatterns << sampler.mFixedMat;

@@ -590,6 +520,7 @@ Archive& operator>>(Archive &ar, GibbsSampler &sampler)

     ar >> sampler.mAMatrix >> sampler.mPMatrix >> sampler.mADomain

         >> sampler.mPDomain >> sampler.mAMeanMatrix >> sampler.mAStdMatrix

         >> sampler.mPMeanMatrix >> sampler.mPStdMatrix >> sampler.mStatUpdates

+        >> sampler.mPumpMatrix >> sampler.mPumpThreshold >> sampler.mPumpStatUpdates

         >> sampler.mMaxGibbsMassA >> sampler.mMaxGibbsMassP

         >> sampler.mAnnealingTemp >> sampler.mSingleCellRNASeq

         >> sampler.mNumFixedPatterns >> sampler.mFixedMat;

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

 #include <vector>

+enum PumpThreshold

+{

+    PUMP_UNIQUE=1,

+    PUMP_CUT=2

+};

+

 class GibbsSampler

 {

 private:

@@ -22,11 +28,11 @@ public:

     ColMatrix mAMeanMatrix, mAStdMatrix;

     RowMatrix mPMeanMatrix, mPStdMatrix;

+    unsigned mStatUpdates;

     ColMatrix mPumpMatrix;

-    Vector mLP;

-

-    unsigned mStatUpdates;

+    PumpThreshold mPumpThreshold;

+    unsigned mPumpStatUpdates;

     float mMaxGibbsMassA;

     float mMaxGibbsMassP;

@@ -60,11 +66,10 @@ public:

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

         unsigned nFactor);

-    GibbsSampler(const Rcpp::NumericMatrix &D,

-        const Rcpp::NumericMatrix &S, unsigned nFactor, float alphaA,

-        float alphaP, float maxGibbmassA, float maxGibbmassP,

-        bool singleCellRNASeq, char whichMatrixFixed,

-        const Rcpp::NumericMatrix &FP);

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

+        unsigned nFactor, float alphaA, float alphaP, float maxGibbmassA,

+        float maxGibbmassP, bool singleCellRNASeq, char whichMatrixFixed,

+        const Rcpp::NumericMatrix &FP, PumpThreshold pumpThreshold);

     void update(char matrixLabel);

@@ -72,21 +77,26 @@ public:

     void setAnnealingTemp(float temp);

     float chi2() const;

+    ColMatrix normedAMatrix() const;

+    RowMatrix normedPMatrix() const;

+

+    unsigned nRow() const {return mDMatrix.nRow();}

+    unsigned nCol() const {return mDMatrix.nCol();}

+    unsigned nFactor() const {return mAMatrix.nCol();}

+

+    // statistics

+    void updateStatistics();

+    void updatePumpStatistics();

     Rcpp::NumericMatrix AMeanRMatrix() const;

     Rcpp::NumericMatrix AStdRMatrix() const;

     Rcpp::NumericMatrix PMeanRMatrix() const;

     Rcpp::NumericMatrix PStdRMatrix() const;

+    Rcpp::NumericMatrix pumpMatrix() const;

+    Rcpp::NumericMatrix meanPattern();

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

     float meanChiSq() const;

-    void updateStatistics();

-    void updatePumpStatistics();

-

-    Rcpp::NumericMatrix getNormedMatrix(char mat);

-

-    unsigned nRow() const {return mDMatrix.nRow();}

-    unsigned nCol() const {return mDMatrix.nCol();}

-    unsigned nFactor() const {return mAMatrix.nCol();}

-

+    // serialization

     friend Archive& operator<<(Archive &ar, GibbsSampler &sampler);

     friend Archive& operator>>(Archive &ar, GibbsSampler &sampler);

 };

@@ -3,6 +3,7 @@

 #include "Archive.h"

 #include "Matrix.h"

+#include "GibbsSampler.h"

 #include <Rcpp.h>

@@ -46,7 +47,7 @@ struct GapsInternalState

     unsigned nUpdatesP;

     GibbsSampler sampler;

-

+    

     SnapshotList snapshotsA;

     SnapshotList snapshotsP;

@@ -55,7 +56,8 @@ struct GapsInternalState

         unsigned nS, unsigned nOut, unsigned nSnap, float alphaA, float alphaP,

         float maxGibbmassA, float maxGibbmassP, int sd, bool msgs,

         bool singleCellRNASeq, char whichMatrixFixed,

-        const Rcpp::NumericMatrix &FP, unsigned cptInterval)

+        const Rcpp::NumericMatrix &FP, unsigned cptInterval,

+        PumpThreshold pumpThreshold)

             :

         chi2VecEquil(nE), nAtomsAEquil(nE), nAtomsPEquil(nE),

         chi2VecSample(nS), nAtomsASample(nS), nAtomsPSample(nS),

@@ -64,7 +66,7 @@ struct GapsInternalState

         phase(GAPS_BURN), seed(sd), checkpointInterval(cptInterval),

         nUpdatesA(0), nUpdatesP(0),

         sampler(D, S, nF, alphaA, alphaP, maxGibbmassA, maxGibbmassP,

-            singleCellRNASeq, whichMatrixFixed, FP)

+            singleCellRNASeq, whichMatrixFixed, FP, pumpThreshold)

     {}

     GapsInternalState(const Rcpp::NumericMatrix &D,

@@ -50,6 +50,87 @@ void Vector::operator+=(const Vector &vec)

     }

 }

+void Vector::operator=(const Vector &vec)

+{

+    mValues = vec.mValues;

+}

+

+Vector Vector::operator+(Vector vec) const

+{