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

 #include "data_structures/Matrix.h"

 #include "utils/GapsAssert.h"

-class Archive;

-

 #define GAPS_SQ(x) ((x) * (x))

+class Archive;

+

 class GapsStatistics

 {

 public:

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

     void update(const DataModel &AModel, const DataModel &PModel);

     template <class DataModel>

     void updatePump(const DataModel &AModel);

+    template <class DataModel>

+    void takeSnapshot(GapsAlgorithmPhase whichPhase, const DataModel &AModel, const DataModel &PModel);

     Matrix Amean() const;

     Matrix Pmean() const;

     Matrix Asd() const;

@@ -33,8 +35,8 @@ public:

     std::vector<unsigned> atomHistory(char m) const;

     float meanChiSq(const DenseNormalModel &model) const;

     float meanChiSq(const SparseNormalModel &model) const;

-    void takeSnapshot();

-    const std::vector<Matrix>& getSnapshots(char whichMatrix) const;

+    const std::vector<Matrix>& getEquilibrationSnapshots(char whichMatrix) const;

+    const std::vector<Matrix>& getSamplingSnapshots(char whichMatrix) const;

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

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

 private:

@@ -43,8 +45,10 @@ private:

     Matrix mPMeanMatrix;

     Matrix mPStdMatrix;

     Matrix mPumpMatrix;

-    std::vector<Matrix> mSnapshotsA;

-    std::vector<Matrix> mSnapshotsP;

+    std::vector<Matrix> mEquilibrationSnapshotsA;

+    std::vector<Matrix> mEquilibrationSnapshotsP;

+    std::vector<Matrix> mSamplingSnapshotsA;

+    std::vector<Matrix> mSamplingSnapshotsP;

     std::vector<float> mChisqHistory;

     std::vector<unsigned> mAtomHistoryA;

     std::vector<unsigned> mAtomHistoryP;

@@ -140,4 +144,20 @@ void GapsStatistics::update(const DataModel &AModel, const DataModel &PModel)

     }

 }

-#endif

\ No newline at end of file

+template <class DataModel>

+void GapsStatistics::takeSnapshot(GapsAlgorithmPhase whichPhase, const DataModel &AModel,

+const DataModel &PModel)

+{

+    if (whichPhase == GAPS_EQUILIBRATION_PHASE)

+    {

+        mEquilibrationSnapshotsA.push_back(AModel.mMatrix.getMatrix());

+        mEquilibrationSnapshotsP.push_back(PModel.mMatrix.getMatrix());

+    }

+    else if (whichPhase == GAPS_SAMPLING_PHASE)

+    {

+        mSamplingSnapshotsA.push_back(AModel.mMatrix.getMatrix());

+        mSamplingSnapshotsP.push_back(PModel.mMatrix.getMatrix());

+    }

+}

+

+#endif // __COGAPS_GAPS_STATISTICS_H__

\ No newline at end of file

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

     std::vector<unsigned> atomHistory(char m) const;

     float meanChiSq(const DenseNormalModel &model) const;

     float meanChiSq(const SparseNormalModel &model) const;

+    void takeSnapshot();

+    const std::vector<Matrix>& getSnapshots(char whichMatrix) const;

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

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

 private:

@@ -41,6 +43,8 @@ private:

     Matrix mPMeanMatrix;

     Matrix mPStdMatrix;

     Matrix mPumpMatrix;

+    std::vector<Matrix> mSnapshotsA;

+    std::vector<Matrix> mSnapshotsP;

     std::vector<float> mChisqHistory;

     std::vector<unsigned> mAtomHistoryA;

     std::vector<unsigned> mAtomHistoryP;

@@ -1,63 +1,50 @@

 #ifndef __COGAPS_GAPS_STATISTICS_H__

 #define __COGAPS_GAPS_STATISTICS_H__

-#include "GapsParameters.h"

 #include "gibbs_sampler/DenseNormalModel.h"

 #include "gibbs_sampler/SparseNormalModel.h"

 #include "math/Math.h"

 #include "math/MatrixMath.h"

 #include "math/VectorMath.h"

 #include "data_structures/Matrix.h"

+#include "utils/GapsAssert.h"

+

+class Archive;

 #define GAPS_SQ(x) ((x) * (x))

 class GapsStatistics

 {

 public:

-

     GapsStatistics(unsigned nGenes, unsigned nSamples, unsigned nPatterns);

-

     template <class DataModel>

     void update(const DataModel &AModel, const DataModel &PModel);

-

     template <class DataModel>

     void updatePump(const DataModel &AModel);

-

     Matrix Amean() const;

     Matrix Pmean() const;

     Matrix Asd() const;

     Matrix Psd() const;

-

     Matrix pumpMatrix() const;

     Matrix meanPattern() const;

-

     void addChiSq(float chisq);

     void addAtomCount(unsigned atomA, unsigned atomP);

-

     std::vector<float> chisqHistory() const;

     std::vector<unsigned> atomHistory(char m) const;

-    

     float meanChiSq(const DenseNormalModel &model) const;

     float meanChiSq(const SparseNormalModel &model) const;

-

-    // serialization

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

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

-

 private:

-

     Matrix mAMeanMatrix;

     Matrix mAStdMatrix;

     Matrix mPMeanMatrix;

     Matrix mPStdMatrix;

     Matrix mPumpMatrix;

-   

     std::vector<float> mChisqHistory;

     std::vector<unsigned> mAtomHistoryA;

     std::vector<unsigned> mAtomHistoryP;

-

     PumpThreshold mPumpThreshold;

-

     unsigned mStatUpdates;

     unsigned mNumPatterns;

     unsigned mPumpUpdates;

@@ -68,7 +55,6 @@ inline void pumpMatrixUniqueThreshold(const FactorMatrix &AMatrix, Matrix *statM

 {

     GAPS_ASSERT(statMatrix->nRow() == AMatrix.nRow());

     GAPS_ASSERT(statMatrix->nCol() == AMatrix.nCol());

-

     std::vector<float> maxValues(AMatrix.nRow(), 0.f);

     std::vector<unsigned> maxIndices(AMatrix.nRow(), 0);

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

@@ -82,10 +68,9 @@ inline void pumpMatrixUniqueThreshold(const FactorMatrix &AMatrix, Matrix *statM

             }

         }

     }

-

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

     {

-        statMatrix->operator()(i,maxIndices[i])++;

+        statMatrix->operator()(i, maxIndices[i])++;

     }

 }

@@ -95,11 +80,9 @@ inline void pumpMatrixCutThreshold(const FactorMatrix &AMatrix, Matrix *statMatr

 {

     GAPS_ASSERT(statMatrix->nRow() == AMatrix.nRow());

     GAPS_ASSERT(statMatrix->nCol() == AMatrix.nCol());

-

-    // we need to access data in columns due to matrix data layout

     std::vector<float> maxValues(AMatrix.nRow(), 0.f);

     std::vector<unsigned> maxIndices(AMatrix.nRow(), 0);

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

+    for (unsigned j = 0; j < AMatrix.nCol(); ++j) // col-major ordering

     {

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

         {

@@ -110,7 +93,6 @@ inline void pumpMatrixCutThreshold(const FactorMatrix &AMatrix, Matrix *statMatr

             }

         }

     }

-

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

     {

         statMatrix->operator()(i,maxIndices[i])++;

@@ -131,29 +113,26 @@ void GapsStatistics::updatePump(const DataModel &AModel)

     }

 }

+// TODO is there precision loss? use double?

 template <class DataModel>

 void GapsStatistics::update(const DataModel &AModel, const DataModel &PModel)

 {

     GAPS_ASSERT(mNumPatterns == AModel.mMatrix.nCol());

     GAPS_ASSERT(mNumPatterns == PModel.mMatrix.nCol());

-

-    // precision loss? use double?

     ++mStatUpdates;

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

     {

         float norm = gaps::max(PModel.mMatrix.getCol(j));

-        norm = norm == 0.f ? 1.f : norm;

-        GAPS_ASSERT(norm > 0.f);

-

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

         Vector quot(PModel.mMatrix.getCol(j) / norm);

-        GAPS_ASSERT(gaps::min(quot) >= 0.f);

         mPMeanMatrix.getCol(j) += quot;

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

-

         Vector prod(AModel.mMatrix.getCol(j) * norm);

-        GAPS_ASSERT(gaps::min(prod) >= 0.f);

         mAMeanMatrix.getCol(j) += prod;

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

+        GAPS_ASSERT(norm > 0.f);

+        GAPS_ASSERT(gaps::min(quot) >= 0.f);

+        GAPS_ASSERT(gaps::min(prod) >= 0.f);

     }

 }

@@ -44,9 +44,7 @@ public:

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

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

-#ifndef GAPS_INTERNAL_TESTS

 private:

-#endif

     Matrix mAMeanMatrix;

     Matrix mAStdMatrix;

@@ -136,12 +136,11 @@ void GapsStatistics::updatePump(const DataModel &AModel)

 template <class DataModel>

 void GapsStatistics::update(const DataModel &AModel, const DataModel &PModel)

 {

-    ++mStatUpdates;

-

     GAPS_ASSERT(mNumPatterns == AModel.mMatrix.nCol());

     GAPS_ASSERT(mNumPatterns == PModel.mMatrix.nCol());

     // precision loss? use double?

+    ++mStatUpdates;

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

     {

         float norm = gaps::max(PModel.mMatrix.getCol(j));

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

 #include "GapsParameters.h"

 #include "gibbs_sampler/DenseNormalModel.h"

-#include "gibbs_sampler/DenseNormalWithUncertaintyModel.h"

 #include "gibbs_sampler/SparseNormalModel.h"

 #include "math/Math.h"

 #include "math/MatrixMath.h"

@@ -40,7 +39,6 @@ public:

     float meanChiSq(const DenseNormalModel &model) const;

     float meanChiSq(const SparseNormalModel &model) const;

-    float meanChiSq(const DenseNormalWithUncertaintyModel &model) const;

     // serialization

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

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

 #include "GapsParameters.h"

 #include "gibbs_sampler/DenseNormalModel.h"

+#include "gibbs_sampler/DenseNormalWithUncertaintyModel.h"

 #include "gibbs_sampler/SparseNormalModel.h"

 #include "math/Math.h"

 #include "math/MatrixMath.h"

@@ -39,6 +40,7 @@ public:

     float meanChiSq(const DenseNormalModel &model) const;

     float meanChiSq(const SparseNormalModel &model) const;

+    float meanChiSq(const DenseNormalWithUncertaintyModel &model) const;

     // serialization

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

@@ -2,10 +2,8 @@

 #define __COGAPS_GAPS_STATISTICS_H__

 #include "GapsParameters.h"

-#include "gibbs_sampler/GibbsSampler.h"

-#include "gibbs_sampler/SingleThreadedGibbsSampler.h"

-#include "gibbs_sampler/DenseStoragePolicy.h"

-#include "gibbs_sampler/SparseStoragePolicy.h"

+#include "gibbs_sampler/DenseNormalModel.h"

+#include "gibbs_sampler/SparseNormalModel.h"

 #include "math/Math.h"

 #include "math/MatrixMath.h"

 #include "math/VectorMath.h"

@@ -19,11 +17,11 @@ public:

     GapsStatistics(unsigned nGenes, unsigned nSamples, unsigned nPatterns);

-    template <class Sampler>

-    void update(const Sampler &ASampler, const Sampler &PSampler);

+    template <class DataModel>

+    void update(const DataModel &AModel, const DataModel &PModel);

-    template <class Sampler>

-    void updatePump(const Sampler &ASampler);

+    template <class DataModel>

+    void updatePump(const DataModel &AModel);

     Matrix Amean() const;

     Matrix Pmean() const;

@@ -39,14 +37,8 @@ public:

     std::vector<float> chisqHistory() const;

     std::vector<unsigned> atomHistory(char m) const;

-    // TODO covert to these

-    float meanChiSq(const DenseStorage &storage) const;

-    float meanChiSq(const SparseStorage &storage) const;

-

-    float meanChiSq(const GibbsSampler<DenseStorage> &PSampler) const;

-    float meanChiSq(const GibbsSampler<SparseStorage> &PSampler) const;

-    float meanChiSq(const SingleThreadedGibbsSampler<DenseStorage> &PSampler) const;

-    float meanChiSq(const SingleThreadedGibbsSampler<SparseStorage> &PSampler) const;

+    float meanChiSq(const DenseNormalModel &model) const;

+    float meanChiSq(const SparseNormalModel &model) const;

     // serialization

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

@@ -127,41 +119,41 @@ inline void pumpMatrixCutThreshold(const FactorMatrix &AMatrix, Matrix *statMatr

     }

 }

-template <class Sampler>

-void GapsStatistics::updatePump(const Sampler &ASampler)

+template <class DataModel>

+void GapsStatistics::updatePump(const DataModel &AModel)

 {

     ++mPumpUpdates;

     if (mPumpThreshold == PUMP_UNIQUE)

     {

-        pumpMatrixCutThreshold(ASampler.mMatrix, &mPumpMatrix);

+        pumpMatrixCutThreshold(AModel.mMatrix, &mPumpMatrix);

     }

     else

     {

-        pumpMatrixUniqueThreshold(ASampler.mMatrix, &mPumpMatrix);

+        pumpMatrixUniqueThreshold(AModel.mMatrix, &mPumpMatrix);

     }

 }

-template <class Sampler>

-void GapsStatistics::update(const Sampler &ASampler, const Sampler &PSampler)

+template <class DataModel>

+void GapsStatistics::update(const DataModel &AModel, const DataModel &PModel)

 {

-    mStatUpdates++;

+    ++mStatUpdates;

-    GAPS_ASSERT(mNumPatterns == ASampler.mMatrix.nCol());

-    GAPS_ASSERT(mNumPatterns == PSampler.mMatrix.nCol());

+    GAPS_ASSERT(mNumPatterns == AModel.mMatrix.nCol());

+    GAPS_ASSERT(mNumPatterns == PModel.mMatrix.nCol());

     // precision loss? use double?

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

     {

-        float norm = gaps::max(PSampler.mMatrix.getCol(j));

+        float norm = gaps::max(PModel.mMatrix.getCol(j));

         norm = norm == 0.f ? 1.f : norm;

         GAPS_ASSERT(norm > 0.f);

-        Vector quot(PSampler.mMatrix.getCol(j) / norm);

+        Vector quot(PModel.mMatrix.getCol(j) / norm);

         GAPS_ASSERT(gaps::min(quot) >= 0.f);

         mPMeanMatrix.getCol(j) += quot;

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

-        Vector prod(ASampler.mMatrix.getCol(j) * norm);

+        Vector prod(AModel.mMatrix.getCol(j) * norm);

         GAPS_ASSERT(gaps::min(prod) >= 0.f);

         mAMeanMatrix.getCol(j) += prod;

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

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

 #include "GapsParameters.h"

 #include "gibbs_sampler/GibbsSampler.h"

+#include "gibbs_sampler/SingleThreadedGibbsSampler.h"

 #include "gibbs_sampler/DenseStoragePolicy.h"

 #include "gibbs_sampler/SparseStoragePolicy.h"

 #include "math/Math.h"

@@ -38,8 +39,14 @@ public:

     std::vector<float> chisqHistory() const;

     std::vector<unsigned> atomHistory(char m) const;

+    // TODO covert to these

+    float meanChiSq(const DenseStorage &storage) const;

+    float meanChiSq(const SparseStorage &storage) const;

+

     float meanChiSq(const GibbsSampler<DenseStorage> &PSampler) const;

     float meanChiSq(const GibbsSampler<SparseStorage> &PSampler) const;

+    float meanChiSq(const SingleThreadedGibbsSampler<DenseStorage> &PSampler) const;

+    float meanChiSq(const SingleThreadedGibbsSampler<SparseStorage> &PSampler) const;

     // serialization

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

@@ -1,12 +1,14 @@

 #ifndef __COGAPS_GAPS_STATISTICS_H__

 #define __COGAPS_GAPS_STATISTICS_H__

+#include "GapsParameters.h"

+#include "gibbs_sampler/GibbsSampler.h"

+#include "gibbs_sampler/DenseStoragePolicy.h"

+#include "gibbs_sampler/SparseStoragePolicy.h"

 #include "math/Math.h"

 #include "math/MatrixMath.h"

 #include "math/VectorMath.h"

 #include "data_structures/Matrix.h"

-#include "gibbs_sampler/DenseGibbsSampler.h"

-#include "gibbs_sampler/SparseGibbsSampler.h"

 #define GAPS_SQ(x) ((x) * (x))

@@ -36,8 +38,8 @@ public:

     std::vector<float> chisqHistory() const;

     std::vector<unsigned> atomHistory(char m) const;

-    float meanChiSq(const DenseGibbsSampler &PSampler) const;

-    float meanChiSq(const SparseGibbsSampler &PSampler) const;

+    float meanChiSq(const GibbsSampler<DenseStorage> &PSampler) const;

+    float meanChiSq(const GibbsSampler<SparseStorage> &PSampler) const;

     // serialization

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

@@ -19,11 +19,17 @@ public:

     template <class Sampler>

     void update(const Sampler &ASampler, const Sampler &PSampler);

+    template <class Sampler>

+    void updatePump(const Sampler &ASampler);

+

     Matrix Amean() const;

     Matrix Pmean() const;

     Matrix Asd() const;

     Matrix Psd() const;

+    Matrix pumpMatrix() const;

+    Matrix meanPattern() const;

+

     void addChiSq(float chisq);

     void addAtomCount(unsigned atomA, unsigned atomP);

@@ -45,15 +51,87 @@ private:

     Matrix mAStdMatrix;

     Matrix mPMeanMatrix;

     Matrix mPStdMatrix;

-    

+    Matrix mPumpMatrix;

+   

     std::vector<float> mChisqHistory;

     std::vector<unsigned> mAtomHistoryA;

     std::vector<unsigned> mAtomHistoryP;

+    PumpThreshold mPumpThreshold;

+

     unsigned mStatUpdates;

     unsigned mNumPatterns;

+    unsigned mPumpUpdates;

 };

+template <class FactorMatrix>

+inline void pumpMatrixUniqueThreshold(const FactorMatrix &AMatrix, Matrix *statMatrix)

+{

+    GAPS_ASSERT(statMatrix->nRow() == AMatrix.nRow());

+    GAPS_ASSERT(statMatrix->nCol() == AMatrix.nCol());

+

+    std::vector<float> maxValues(AMatrix.nRow(), 0.f);

+    std::vector<unsigned> maxIndices(AMatrix.nRow(), 0);

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

+    {

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

+        {

+            if (maxValues[i] < AMatrix(i,j))

+            {

+                maxValues[i] = AMatrix(i,j);

+                maxIndices[i] = j;

+            }

+        }

+    }

+

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

+    {

+        statMatrix->operator()(i,maxIndices[i])++;

+    }

+}

+

+// implemented the same as unique for now - TODO

+template <class FactorMatrix>

+inline void pumpMatrixCutThreshold(const FactorMatrix &AMatrix, Matrix *statMatrix)

+{

+    GAPS_ASSERT(statMatrix->nRow() == AMatrix.nRow());

+    GAPS_ASSERT(statMatrix->nCol() == AMatrix.nCol());

+

+    // we need to access data in columns due to matrix data layout

+    std::vector<float> maxValues(AMatrix.nRow(), 0.f);

+    std::vector<unsigned> maxIndices(AMatrix.nRow(), 0);

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

+    {

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

+        {

+            if (maxValues[i] < AMatrix(i,j))

+            {

+                maxValues[i] = AMatrix(i,j);

+                maxIndices[i] = j;

+            }

+        }

+    }

+

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

+    {

+        statMatrix->operator()(i,maxIndices[i])++;

+    }

+}

+

+template <class Sampler>

+void GapsStatistics::updatePump(const Sampler &ASampler)

+{

+    ++mPumpUpdates;

+    if (mPumpThreshold == PUMP_UNIQUE)

+    {

+        pumpMatrixCutThreshold(ASampler.mMatrix, &mPumpMatrix);

+    }

+    else

+    {

+        pumpMatrixUniqueThreshold(ASampler.mMatrix, &mPumpMatrix);

+    }

+}

+

 template <class Sampler>

 void GapsStatistics::update(const Sampler &ASampler, const Sampler &PSampler)

 {

@@ -23,6 +23,12 @@ public:

     Matrix Pmean() const;

     Matrix Asd() const;

     Matrix Psd() const;

+

+    void addChiSq(float chisq);

+    void addAtomCount(unsigned atomA, unsigned atomP);

+

+    std::vector<float> chisqHistory() const;

+    std::vector<unsigned> atomHistory(char m) const;

     float meanChiSq(const DenseGibbsSampler &PSampler) const;

     float meanChiSq(const SparseGibbsSampler &PSampler) const;

@@ -40,6 +46,10 @@ private:

     Matrix mPMeanMatrix;

     Matrix mPStdMatrix;

+    std::vector<float> mChisqHistory;

+    std::vector<unsigned> mAtomHistoryA;

+    std::vector<unsigned> mAtomHistoryP;

+

     unsigned mStatUpdates;

     unsigned mNumPatterns;

 };

old mode 100644

new mode 100755

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

 class GapsStatistics

 {

-private:

-

-    Matrix mAMeanMatrix;

-    Matrix mAStdMatrix;

-    Matrix mPMeanMatrix;

-    Matrix mPStdMatrix;

-    

-    unsigned mStatUpdates;

-    unsigned mNumPatterns;

-

 public:

     GapsStatistics(unsigned nGenes, unsigned nSamples, unsigned nPatterns);

@@ -40,6 +30,18 @@ public:

     // serialization

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

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

+

+#ifndef GAPS_INTERNAL_TESTS

+private:

+#endif

+

+    Matrix mAMeanMatrix;

+    Matrix mAStdMatrix;

+    Matrix mPMeanMatrix;

+    Matrix mPStdMatrix;

+    

+    unsigned mStatUpdates;

+    unsigned mNumPatterns;

 };

 template <class Sampler>

@@ -47,31 +47,26 @@ void GapsStatistics::update(const Sampler &ASampler, const Sampler &PSampler)

 {

     mStatUpdates++;

-    // update     

-    // precision loss? use double?

-    DEBUG_PING

     GAPS_ASSERT(mNumPatterns == ASampler.mMatrix.nCol());

     GAPS_ASSERT(mNumPatterns == PSampler.mMatrix.nCol());

+    // precision loss? use double?

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

     {

         float norm = gaps::max(PSampler.mMatrix.getCol(j));

         norm = norm == 0.f ? 1.f : norm;

         GAPS_ASSERT(norm > 0.f);

-        DEBUG_PING

         Vector quot(PSampler.mMatrix.getCol(j) / norm);

         GAPS_ASSERT(gaps::min(quot) >= 0.f);

         mPMeanMatrix.getCol(j) += quot;

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

-        DEBUG_PING

         Vector prod(ASampler.mMatrix.getCol(j) * norm);

         GAPS_ASSERT(gaps::min(prod) >= 0.f);

         mAMeanMatrix.getCol(j) += prod;

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

     }

-    DEBUG_PING

 }

 #endif

\ No newline at end of file

@@ -38,7 +38,7 @@ public:

     float meanChiSq(const SparseGibbsSampler &PSampler) const;

     // serialization

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

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

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

 };

@@ -49,22 +49,29 @@ void GapsStatistics::update(const Sampler &ASampler, const Sampler &PSampler)

     // update     

     // precision loss? use double?

+    DEBUG_PING

+    GAPS_ASSERT(mNumPatterns == ASampler.mMatrix.nCol());

+    GAPS_ASSERT(mNumPatterns == PSampler.mMatrix.nCol());

+

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

     {

         float norm = gaps::max(PSampler.mMatrix.getCol(j));

         norm = norm == 0.f ? 1.f : norm;

         GAPS_ASSERT(norm > 0.f);

+        DEBUG_PING

         Vector quot(PSampler.mMatrix.getCol(j) / norm);

         GAPS_ASSERT(gaps::min(quot) >= 0.f);

         mPMeanMatrix.getCol(j) += quot;

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

+        DEBUG_PING

         Vector prod(ASampler.mMatrix.getCol(j) * norm);

         GAPS_ASSERT(gaps::min(prod) >= 0.f);

         mAMeanMatrix.getCol(j) += prod;

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

     }

+    DEBUG_PING

 }

 #endif

\ No newline at end of file

@@ -48,18 +48,22 @@ void GapsStatistics::update(const Sampler &ASampler, const Sampler &PSampler)

     mStatUpdates++;

     // update     

+    // precision loss? use double?

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

     {

         float norm = gaps::max(PSampler.mMatrix.getCol(j));

         norm = norm == 0.f ? 1.f : norm;

+        GAPS_ASSERT(norm > 0.f);

         Vector quot(PSampler.mMatrix.getCol(j) / norm);

+        GAPS_ASSERT(gaps::min(quot) >= 0.f);

         mPMeanMatrix.getCol(j) += quot;

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

         Vector prod(ASampler.mMatrix.getCol(j) * norm);

+        GAPS_ASSERT(gaps::min(prod) >= 0.f);

         mAMeanMatrix.getCol(j) += prod;

-        mAStdMatrix.getCol(j) += gaps::elementSq(prod); // precision loss? use double?

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

     }

 }

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

 #ifndef __COGAPS_GAPS_STATISTICS_H__

 #define __COGAPS_GAPS_STATISTICS_H__

-#include "GibbsSampler.h"

+#include "math/Math.h"

+#include "math/MatrixMath.h"

+#include "math/VectorMath.h"

 #include "data_structures/Matrix.h"

+#include "gibbs_sampler/DenseGibbsSampler.h"

+#include "gibbs_sampler/SparseGibbsSampler.h"

+

+#define GAPS_SQ(x) ((x) * (x))

 class GapsStatistics

 {

@@ -18,20 +24,43 @@ private:

 public:

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

+    GapsStatistics(unsigned nGenes, unsigned nSamples, unsigned nPatterns);

-    void update(const GibbsSampler *ASampler, const GibbsSampler *PSampler);

+    template <class Sampler>

+    void update(const Sampler &ASampler, const Sampler &PSampler);

     Matrix Amean() const;

     Matrix Pmean() const;

     Matrix Asd() const;

     Matrix Psd() const;

-

-    float meanChiSq(const GibbsSampler *PSampler) const;

+    

+    float meanChiSq(const DenseGibbsSampler &PSampler) const;

+    float meanChiSq(const SparseGibbsSampler &PSampler) const;

     // serialization

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

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

 };

+template <class Sampler>

+void GapsStatistics::update(const Sampler &ASampler, const Sampler &PSampler)

+{

+    mStatUpdates++;

+

+    // update     

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

+    {

+        float norm = gaps::max(PSampler.mMatrix.getCol(j));

+        norm = norm == 0.f ? 1.f : norm;

+

+        Vector quot(PSampler.mMatrix.getCol(j) / norm);

+        mPMeanMatrix.getCol(j) += quot;

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

+

+        Vector prod(ASampler.mMatrix.getCol(j) * norm);

+        mAMeanMatrix.getCol(j) += prod;

+        mAStdMatrix.getCol(j) += gaps::elementSq(prod); // precision loss? use double?

+    }

+}

+

 #endif

\ No newline at end of file

@@ -8,10 +8,10 @@ class GapsStatistics

 {

 private:

-    ColMatrix mAMeanMatrix;

-    ColMatrix mAStdMatrix;

-    ColMatrix mPMeanMatrix;

-    ColMatrix mPStdMatrix;

+    Matrix mAMeanMatrix;

+    Matrix mAStdMatrix;

+    Matrix mPMeanMatrix;

+    Matrix mPStdMatrix;

     unsigned mStatUpdates;

     unsigned mNumPatterns;

@@ -20,14 +20,14 @@ public:

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

-    void update(const GibbsSampler &ASampler, const GibbsSampler &PSampler);

+    void update(const GibbsSampler *ASampler, const GibbsSampler *PSampler);

-    ColMatrix Amean() const;

-    ColMatrix Pmean() const;

-    ColMatrix Asd() const;

-    ColMatrix Psd() const;

+    Matrix Amean() const;

+    Matrix Pmean() const;

+    Matrix Asd() const;

+    Matrix Psd() const;

-    float meanChiSq(const GibbsSampler &PSampler) const;

+    float meanChiSq(const GibbsSampler *PSampler) const;

     // serialization

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

@@ -20,15 +20,14 @@ public:

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

-    void update(const AmplitudeGibbsSampler &ASampler,

-        const PatternGibbsSampler &PSampler);

+    void update(const GibbsSampler &ASampler, const GibbsSampler &PSampler);

     ColMatrix Amean() const;

     ColMatrix Pmean() const;

     ColMatrix Asd() const;

     ColMatrix Psd() const;

-    float meanChiSq(const PatternGibbsSampler &ASampler) const;

+    float meanChiSq(const GibbsSampler &PSampler) const;

     // serialization

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

@@ -4,50 +4,31 @@

 #include "GibbsSampler.h"

 #include "data_structures/Matrix.h"

-enum PumpThreshold

-{

-    PUMP_UNIQUE=1,

-    PUMP_CUT=2

-};

-

 class GapsStatistics

 {

 private:

     ColMatrix mAMeanMatrix;

     ColMatrix mAStdMatrix;

-    RowMatrix mPMeanMatrix;

-    RowMatrix mPStdMatrix;

+    ColMatrix mPMeanMatrix;

+    ColMatrix mPStdMatrix;

     unsigned mStatUpdates;