Browse code

Snapshots working for both phases

Tom Sherman authored on 05/09/2019 00:20:30
Showing1 changed files
... ...
@@ -9,10 +9,10 @@
9 9
 #include "data_structures/Matrix.h"
10 10
 #include "utils/GapsAssert.h"
11 11
 
12
-class Archive;
13
-
14 12
 #define GAPS_SQ(x) ((x) * (x))
15 13
 
14
+class Archive;
15
+
16 16
 class GapsStatistics
17 17
 {
18 18
 public:
... ...
@@ -21,6 +21,8 @@ public:
21 21
     void update(const DataModel &AModel, const DataModel &PModel);
22 22
     template <class DataModel>
23 23
     void updatePump(const DataModel &AModel);
24
+    template <class DataModel>
25
+    void takeSnapshot(GapsAlgorithmPhase whichPhase, const DataModel &AModel, const DataModel &PModel);
24 26
     Matrix Amean() const;
25 27
     Matrix Pmean() const;
26 28
     Matrix Asd() const;
... ...
@@ -33,8 +35,8 @@ public:
33 35
     std::vector<unsigned> atomHistory(char m) const;
34 36
     float meanChiSq(const DenseNormalModel &model) const;
35 37
     float meanChiSq(const SparseNormalModel &model) const;
36
-    void takeSnapshot();
37
-    const std::vector<Matrix>& getSnapshots(char whichMatrix) const;
38
+    const std::vector<Matrix>& getEquilibrationSnapshots(char whichMatrix) const;
39
+    const std::vector<Matrix>& getSamplingSnapshots(char whichMatrix) const;
38 40
     friend Archive& operator<<(Archive &ar, const GapsStatistics &stat);
39 41
     friend Archive& operator>>(Archive &ar, GapsStatistics &stat);
40 42
 private:
... ...
@@ -43,8 +45,10 @@ private:
43 45
     Matrix mPMeanMatrix;
44 46
     Matrix mPStdMatrix;
45 47
     Matrix mPumpMatrix;
46
-    std::vector<Matrix> mSnapshotsA;
47
-    std::vector<Matrix> mSnapshotsP;
48
+    std::vector<Matrix> mEquilibrationSnapshotsA;
49
+    std::vector<Matrix> mEquilibrationSnapshotsP;
50
+    std::vector<Matrix> mSamplingSnapshotsA;
51
+    std::vector<Matrix> mSamplingSnapshotsP;
48 52
     std::vector<float> mChisqHistory;
49 53
     std::vector<unsigned> mAtomHistoryA;
50 54
     std::vector<unsigned> mAtomHistoryP;
... ...
@@ -140,4 +144,20 @@ void GapsStatistics::update(const DataModel &AModel, const DataModel &PModel)
140 144
     }
141 145
 }
142 146
 
143
-#endif
144 147
\ No newline at end of file
148
+template <class DataModel>
149
+void GapsStatistics::takeSnapshot(GapsAlgorithmPhase whichPhase, const DataModel &AModel,
150
+const DataModel &PModel)
151
+{
152
+    if (whichPhase == GAPS_EQUILIBRATION_PHASE)
153
+    {
154
+        mEquilibrationSnapshotsA.push_back(AModel.mMatrix.getMatrix());
155
+        mEquilibrationSnapshotsP.push_back(PModel.mMatrix.getMatrix());
156
+    }
157
+    else if (whichPhase == GAPS_SAMPLING_PHASE)
158
+    {
159
+        mSamplingSnapshotsA.push_back(AModel.mMatrix.getMatrix());
160
+        mSamplingSnapshotsP.push_back(PModel.mMatrix.getMatrix());
161
+    }
162
+}
163
+
164
+#endif // __COGAPS_GAPS_STATISTICS_H__
145 165
\ No newline at end of file
Browse code

Added option to take snapshots during sampling

The argument 'nSnapshots' specifies how many samples of the A and P matrix should be saved. The snapshots are equally spaced out during the sampling phase. This is useful for various post-run analysis techniques but should primarily be used to test ideas, not as part of an official analysis.

Tom Sherman authored on 03/09/2019 15:23:59
Showing1 changed files
... ...
@@ -33,6 +33,8 @@ public:
33 33
     std::vector<unsigned> atomHistory(char m) const;
34 34
     float meanChiSq(const DenseNormalModel &model) const;
35 35
     float meanChiSq(const SparseNormalModel &model) const;
36
+    void takeSnapshot();
37
+    const std::vector<Matrix>& getSnapshots(char whichMatrix) const;
36 38
     friend Archive& operator<<(Archive &ar, const GapsStatistics &stat);
37 39
     friend Archive& operator>>(Archive &ar, GapsStatistics &stat);
38 40
 private:
... ...
@@ -41,6 +43,8 @@ private:
41 43
     Matrix mPMeanMatrix;
42 44
     Matrix mPStdMatrix;
43 45
     Matrix mPumpMatrix;
46
+    std::vector<Matrix> mSnapshotsA;
47
+    std::vector<Matrix> mSnapshotsP;
44 48
     std::vector<float> mChisqHistory;
45 49
     std::vector<unsigned> mAtomHistoryA;
46 50
     std::vector<unsigned> mAtomHistoryP;
Browse code

clean up linter warnings

Tom Sherman authored on 24/06/2019 19:44:02
Showing1 changed files
... ...
@@ -1,63 +1,50 @@
1 1
 #ifndef __COGAPS_GAPS_STATISTICS_H__
2 2
 #define __COGAPS_GAPS_STATISTICS_H__
3 3
 
4
-#include "GapsParameters.h"
5 4
 #include "gibbs_sampler/DenseNormalModel.h"
6 5
 #include "gibbs_sampler/SparseNormalModel.h"
7 6
 #include "math/Math.h"
8 7
 #include "math/MatrixMath.h"
9 8
 #include "math/VectorMath.h"
10 9
 #include "data_structures/Matrix.h"
10
+#include "utils/GapsAssert.h"
11
+
12
+class Archive;
11 13
 
12 14
 #define GAPS_SQ(x) ((x) * (x))
13 15
 
14 16
 class GapsStatistics
15 17
 {
16 18
 public:
17
-
18 19
     GapsStatistics(unsigned nGenes, unsigned nSamples, unsigned nPatterns);
19
-
20 20
     template <class DataModel>
21 21
     void update(const DataModel &AModel, const DataModel &PModel);
22
-
23 22
     template <class DataModel>
24 23
     void updatePump(const DataModel &AModel);
25
-
26 24
     Matrix Amean() const;
27 25
     Matrix Pmean() const;
28 26
     Matrix Asd() const;
29 27
     Matrix Psd() const;
30
-
31 28
     Matrix pumpMatrix() const;
32 29
     Matrix meanPattern() const;
33
-
34 30
     void addChiSq(float chisq);
35 31
     void addAtomCount(unsigned atomA, unsigned atomP);
36
-
37 32
     std::vector<float> chisqHistory() const;
38 33
     std::vector<unsigned> atomHistory(char m) const;
39
-    
40 34
     float meanChiSq(const DenseNormalModel &model) const;
41 35
     float meanChiSq(const SparseNormalModel &model) const;
42
-
43
-    // serialization
44 36
     friend Archive& operator<<(Archive &ar, const GapsStatistics &stat);
45 37
     friend Archive& operator>>(Archive &ar, GapsStatistics &stat);
46
-
47 38
 private:
48
-
49 39
     Matrix mAMeanMatrix;
50 40
     Matrix mAStdMatrix;
51 41
     Matrix mPMeanMatrix;
52 42
     Matrix mPStdMatrix;
53 43
     Matrix mPumpMatrix;
54
-   
55 44
     std::vector<float> mChisqHistory;
56 45
     std::vector<unsigned> mAtomHistoryA;
57 46
     std::vector<unsigned> mAtomHistoryP;
58
-
59 47
     PumpThreshold mPumpThreshold;
60
-
61 48
     unsigned mStatUpdates;
62 49
     unsigned mNumPatterns;
63 50
     unsigned mPumpUpdates;
... ...
@@ -68,7 +55,6 @@ inline void pumpMatrixUniqueThreshold(const FactorMatrix &AMatrix, Matrix *statM
68 55
 {
69 56
     GAPS_ASSERT(statMatrix->nRow() == AMatrix.nRow());
70 57
     GAPS_ASSERT(statMatrix->nCol() == AMatrix.nCol());
71
-
72 58
     std::vector<float> maxValues(AMatrix.nRow(), 0.f);
73 59
     std::vector<unsigned> maxIndices(AMatrix.nRow(), 0);
74 60
     for (unsigned j = 0; j < AMatrix.nCol(); ++j)
... ...
@@ -82,10 +68,9 @@ inline void pumpMatrixUniqueThreshold(const FactorMatrix &AMatrix, Matrix *statM
82 68
             }
83 69
         }
84 70
     }
85
-
86 71
     for (unsigned i = 0; i < AMatrix.nRow(); ++i)
87 72
     {
88
-        statMatrix->operator()(i,maxIndices[i])++;
73
+        statMatrix->operator()(i, maxIndices[i])++;
89 74
     }
90 75
 }
91 76
 
... ...
@@ -95,11 +80,9 @@ inline void pumpMatrixCutThreshold(const FactorMatrix &AMatrix, Matrix *statMatr
95 80
 {
96 81
     GAPS_ASSERT(statMatrix->nRow() == AMatrix.nRow());
97 82
     GAPS_ASSERT(statMatrix->nCol() == AMatrix.nCol());
98
-
99
-    // we need to access data in columns due to matrix data layout
100 83
     std::vector<float> maxValues(AMatrix.nRow(), 0.f);
101 84
     std::vector<unsigned> maxIndices(AMatrix.nRow(), 0);
102
-    for (unsigned j = 0; j < AMatrix.nCol(); ++j)
85
+    for (unsigned j = 0; j < AMatrix.nCol(); ++j) // col-major ordering
103 86
     {
104 87
         for (unsigned i = 0; i < AMatrix.nRow(); ++i)
105 88
         {
... ...
@@ -110,7 +93,6 @@ inline void pumpMatrixCutThreshold(const FactorMatrix &AMatrix, Matrix *statMatr
110 93
             }
111 94
         }
112 95
     }
113
-
114 96
     for (unsigned i = 0; i < AMatrix.nRow(); ++i)
115 97
     {
116 98
         statMatrix->operator()(i,maxIndices[i])++;
... ...
@@ -131,29 +113,26 @@ void GapsStatistics::updatePump(const DataModel &AModel)
131 113
     }
132 114
 }
133 115
 
116
+// TODO is there precision loss? use double?
134 117
 template <class DataModel>
135 118
 void GapsStatistics::update(const DataModel &AModel, const DataModel &PModel)
136 119
 {
137 120
     GAPS_ASSERT(mNumPatterns == AModel.mMatrix.nCol());
138 121
     GAPS_ASSERT(mNumPatterns == PModel.mMatrix.nCol());
139
-
140
-    // precision loss? use double?
141 122
     ++mStatUpdates;
142 123
     for (unsigned j = 0; j < mNumPatterns; ++j)
143 124
     {
144 125
         float norm = gaps::max(PModel.mMatrix.getCol(j));
145
-        norm = norm == 0.f ? 1.f : norm;
146
-        GAPS_ASSERT(norm > 0.f);
147
-
126
+        norm = (norm == 0.f) ? 1.f : norm;
148 127
         Vector quot(PModel.mMatrix.getCol(j) / norm);
149
-        GAPS_ASSERT(gaps::min(quot) >= 0.f);
150 128
         mPMeanMatrix.getCol(j) += quot;
151 129
         mPStdMatrix.getCol(j) += gaps::elementSq(quot);
152
-
153 130
         Vector prod(AModel.mMatrix.getCol(j) * norm);
154
-        GAPS_ASSERT(gaps::min(prod) >= 0.f);
155 131
         mAMeanMatrix.getCol(j) += prod;
156 132
         mAStdMatrix.getCol(j) += gaps::elementSq(prod);
133
+        GAPS_ASSERT(norm > 0.f);
134
+        GAPS_ASSERT(gaps::min(quot) >= 0.f);
135
+        GAPS_ASSERT(gaps::min(prod) >= 0.f);
157 136
     }
158 137
 }
159 138
 
Browse code

removed internal tests option - this behavior should be tested in debug mode

Tom Sherman authored on 24/06/2019 14:14:03
Showing1 changed files
... ...
@@ -44,9 +44,7 @@ public:
44 44
     friend Archive& operator<<(Archive &ar, const GapsStatistics &stat);
45 45
     friend Archive& operator>>(Archive &ar, GapsStatistics &stat);
46 46
 
47
-#ifndef GAPS_INTERNAL_TESTS
48 47
 private:
49
-#endif
50 48
 
51 49
     Matrix mAMeanMatrix;
52 50
     Matrix mAStdMatrix;
Browse code

tests mostly passing, except when running on many threads and using the sparseOptimization flag there is not consistent results for the same seed

Tom Sherman authored on 21/06/2019 19:44:45
Showing1 changed files
... ...
@@ -136,12 +136,11 @@ void GapsStatistics::updatePump(const DataModel &AModel)
136 136
 template <class DataModel>
137 137
 void GapsStatistics::update(const DataModel &AModel, const DataModel &PModel)
138 138
 {
139
-    ++mStatUpdates;
140
-
141 139
     GAPS_ASSERT(mNumPatterns == AModel.mMatrix.nCol());
142 140
     GAPS_ASSERT(mNumPatterns == PModel.mMatrix.nCol());
143 141
 
144 142
     // precision loss? use double?
143
+    ++mStatUpdates;
145 144
     for (unsigned j = 0; j < mNumPatterns; ++j)
146 145
     {
147 146
         float norm = gaps::max(PModel.mMatrix.getCol(j));
Browse code

all sampler/model combinations seem to be working

Tom Sherman authored on 03/06/2019 17:59:25
Showing1 changed files
... ...
@@ -3,7 +3,6 @@
3 3
 
4 4
 #include "GapsParameters.h"
5 5
 #include "gibbs_sampler/DenseNormalModel.h"
6
-#include "gibbs_sampler/DenseNormalWithUncertaintyModel.h"
7 6
 #include "gibbs_sampler/SparseNormalModel.h"
8 7
 #include "math/Math.h"
9 8
 #include "math/MatrixMath.h"
... ...
@@ -40,7 +39,6 @@ public:
40 39
     
41 40
     float meanChiSq(const DenseNormalModel &model) const;
42 41
     float meanChiSq(const SparseNormalModel &model) const;
43
-    float meanChiSq(const DenseNormalWithUncertaintyModel &model) const;
44 42
 
45 43
     // serialization
46 44
     friend Archive& operator<<(Archive &ar, const GapsStatistics &stat);
Browse code

added separate data model for when uncertainty is passed - not implemented

Tom Sherman authored on 23/05/2019 18:55:10
Showing1 changed files
... ...
@@ -3,6 +3,7 @@
3 3
 
4 4
 #include "GapsParameters.h"
5 5
 #include "gibbs_sampler/DenseNormalModel.h"
6
+#include "gibbs_sampler/DenseNormalWithUncertaintyModel.h"
6 7
 #include "gibbs_sampler/SparseNormalModel.h"
7 8
 #include "math/Math.h"
8 9
 #include "math/MatrixMath.h"
... ...
@@ -39,6 +40,7 @@ public:
39 40
     
40 41
     float meanChiSq(const DenseNormalModel &model) const;
41 42
     float meanChiSq(const SparseNormalModel &model) const;
43
+    float meanChiSq(const DenseNormalWithUncertaintyModel &model) const;
42 44
 
43 45
     // serialization
44 46
     friend Archive& operator<<(Archive &ar, const GapsStatistics &stat);
Browse code

infrastructure in place for extending samplers and models

Tom Sherman authored on 22/05/2019 21:27:14
Showing1 changed files
... ...
@@ -2,10 +2,8 @@
2 2
 #define __COGAPS_GAPS_STATISTICS_H__
3 3
 
4 4
 #include "GapsParameters.h"
5
-#include "gibbs_sampler/GibbsSampler.h"
6
-#include "gibbs_sampler/SingleThreadedGibbsSampler.h"
7
-#include "gibbs_sampler/DenseStoragePolicy.h"
8
-#include "gibbs_sampler/SparseStoragePolicy.h"
5
+#include "gibbs_sampler/DenseNormalModel.h"
6
+#include "gibbs_sampler/SparseNormalModel.h"
9 7
 #include "math/Math.h"
10 8
 #include "math/MatrixMath.h"
11 9
 #include "math/VectorMath.h"
... ...
@@ -19,11 +17,11 @@ public:
19 17
 
20 18
     GapsStatistics(unsigned nGenes, unsigned nSamples, unsigned nPatterns);
21 19
 
22
-    template <class Sampler>
23
-    void update(const Sampler &ASampler, const Sampler &PSampler);
20
+    template <class DataModel>
21
+    void update(const DataModel &AModel, const DataModel &PModel);
24 22
 
25
-    template <class Sampler>
26
-    void updatePump(const Sampler &ASampler);
23
+    template <class DataModel>
24
+    void updatePump(const DataModel &AModel);
27 25
 
28 26
     Matrix Amean() const;
29 27
     Matrix Pmean() const;
... ...
@@ -39,14 +37,8 @@ public:
39 37
     std::vector<float> chisqHistory() const;
40 38
     std::vector<unsigned> atomHistory(char m) const;
41 39
     
42
-    // TODO covert to these
43
-    float meanChiSq(const DenseStorage &storage) const;
44
-    float meanChiSq(const SparseStorage &storage) const;
45
-
46
-    float meanChiSq(const GibbsSampler<DenseStorage> &PSampler) const;
47
-    float meanChiSq(const GibbsSampler<SparseStorage> &PSampler) const;
48
-    float meanChiSq(const SingleThreadedGibbsSampler<DenseStorage> &PSampler) const;
49
-    float meanChiSq(const SingleThreadedGibbsSampler<SparseStorage> &PSampler) const;
40
+    float meanChiSq(const DenseNormalModel &model) const;
41
+    float meanChiSq(const SparseNormalModel &model) const;
50 42
 
51 43
     // serialization
52 44
     friend Archive& operator<<(Archive &ar, const GapsStatistics &stat);
... ...
@@ -127,41 +119,41 @@ inline void pumpMatrixCutThreshold(const FactorMatrix &AMatrix, Matrix *statMatr
127 119
     }
128 120
 }
129 121
 
130
-template <class Sampler>
131
-void GapsStatistics::updatePump(const Sampler &ASampler)
122
+template <class DataModel>
123
+void GapsStatistics::updatePump(const DataModel &AModel)
132 124
 {
133 125
     ++mPumpUpdates;
134 126
     if (mPumpThreshold == PUMP_UNIQUE)
135 127
     {
136
-        pumpMatrixCutThreshold(ASampler.mMatrix, &mPumpMatrix);
128
+        pumpMatrixCutThreshold(AModel.mMatrix, &mPumpMatrix);
137 129
     }
138 130
     else
139 131
     {
140
-        pumpMatrixUniqueThreshold(ASampler.mMatrix, &mPumpMatrix);
132
+        pumpMatrixUniqueThreshold(AModel.mMatrix, &mPumpMatrix);
141 133
     }
142 134
 }
143 135
 
144
-template <class Sampler>
145
-void GapsStatistics::update(const Sampler &ASampler, const Sampler &PSampler)
136
+template <class DataModel>
137
+void GapsStatistics::update(const DataModel &AModel, const DataModel &PModel)
146 138
 {
147
-    mStatUpdates++;
139
+    ++mStatUpdates;
148 140
 
149
-    GAPS_ASSERT(mNumPatterns == ASampler.mMatrix.nCol());
150
-    GAPS_ASSERT(mNumPatterns == PSampler.mMatrix.nCol());
141
+    GAPS_ASSERT(mNumPatterns == AModel.mMatrix.nCol());
142
+    GAPS_ASSERT(mNumPatterns == PModel.mMatrix.nCol());
151 143
 
152 144
     // precision loss? use double?
153 145
     for (unsigned j = 0; j < mNumPatterns; ++j)
154 146
     {
155
-        float norm = gaps::max(PSampler.mMatrix.getCol(j));
147
+        float norm = gaps::max(PModel.mMatrix.getCol(j));
156 148
         norm = norm == 0.f ? 1.f : norm;
157 149
         GAPS_ASSERT(norm > 0.f);
158 150
 
159
-        Vector quot(PSampler.mMatrix.getCol(j) / norm);
151
+        Vector quot(PModel.mMatrix.getCol(j) / norm);
160 152
         GAPS_ASSERT(gaps::min(quot) >= 0.f);
161 153
         mPMeanMatrix.getCol(j) += quot;
162 154
         mPStdMatrix.getCol(j) += gaps::elementSq(quot);
163 155
 
164
-        Vector prod(ASampler.mMatrix.getCol(j) * norm);
156
+        Vector prod(AModel.mMatrix.getCol(j) * norm);
165 157
         GAPS_ASSERT(gaps::min(prod) >= 0.f);
166 158
         mAMeanMatrix.getCol(j) += prod;
167 159
         mAStdMatrix.getCol(j) += gaps::elementSq(prod);
Browse code

single-threaded sampler working

Tom Sherman authored on 22/05/2019 17:38:53
Showing1 changed files
... ...
@@ -3,6 +3,7 @@
3 3
 
4 4
 #include "GapsParameters.h"
5 5
 #include "gibbs_sampler/GibbsSampler.h"
6
+#include "gibbs_sampler/SingleThreadedGibbsSampler.h"
6 7
 #include "gibbs_sampler/DenseStoragePolicy.h"
7 8
 #include "gibbs_sampler/SparseStoragePolicy.h"
8 9
 #include "math/Math.h"
... ...
@@ -38,8 +39,14 @@ public:
38 39
     std::vector<float> chisqHistory() const;
39 40
     std::vector<unsigned> atomHistory(char m) const;
40 41
     
42
+    // TODO covert to these
43
+    float meanChiSq(const DenseStorage &storage) const;
44
+    float meanChiSq(const SparseStorage &storage) const;
45
+
41 46
     float meanChiSq(const GibbsSampler<DenseStorage> &PSampler) const;
42 47
     float meanChiSq(const GibbsSampler<SparseStorage> &PSampler) const;
48
+    float meanChiSq(const SingleThreadedGibbsSampler<DenseStorage> &PSampler) const;
49
+    float meanChiSq(const SingleThreadedGibbsSampler<SparseStorage> &PSampler) const;
43 50
 
44 51
     // serialization
45 52
     friend Archive& operator<<(Archive &ar, const GapsStatistics &stat);
Browse code

implemented policy based design for sparse/dense storage in the GibbsSampler

Tom Sherman authored on 08/01/2019 22:16:52
Showing1 changed files
... ...
@@ -1,12 +1,14 @@
1 1
 #ifndef __COGAPS_GAPS_STATISTICS_H__
2 2
 #define __COGAPS_GAPS_STATISTICS_H__
3 3
 
4
+#include "GapsParameters.h"
5
+#include "gibbs_sampler/GibbsSampler.h"
6
+#include "gibbs_sampler/DenseStoragePolicy.h"
7
+#include "gibbs_sampler/SparseStoragePolicy.h"
4 8
 #include "math/Math.h"
5 9
 #include "math/MatrixMath.h"
6 10
 #include "math/VectorMath.h"
7 11
 #include "data_structures/Matrix.h"
8
-#include "gibbs_sampler/DenseGibbsSampler.h"
9
-#include "gibbs_sampler/SparseGibbsSampler.h"
10 12
 
11 13
 #define GAPS_SQ(x) ((x) * (x))
12 14
 
... ...
@@ -36,8 +38,8 @@ public:
36 38
     std::vector<float> chisqHistory() const;
37 39
     std::vector<unsigned> atomHistory(char m) const;
38 40
     
39
-    float meanChiSq(const DenseGibbsSampler &PSampler) const;
40
-    float meanChiSq(const SparseGibbsSampler &PSampler) const;
41
+    float meanChiSq(const GibbsSampler<DenseStorage> &PSampler) const;
42
+    float meanChiSq(const GibbsSampler<SparseStorage> &PSampler) const;
41 43
 
42 44
     // serialization
43 45
     friend Archive& operator<<(Archive &ar, const GapsStatistics &stat);
Browse code

added PUMP back to COGAPS

Tom Sherman authored on 21/12/2018 00:52:37
Showing1 changed files
... ...
@@ -19,11 +19,17 @@ public:
19 19
     template <class Sampler>
20 20
     void update(const Sampler &ASampler, const Sampler &PSampler);
21 21
 
22
+    template <class Sampler>
23
+    void updatePump(const Sampler &ASampler);
24
+
22 25
     Matrix Amean() const;
23 26
     Matrix Pmean() const;
24 27
     Matrix Asd() const;
25 28
     Matrix Psd() const;
26 29
 
30
+    Matrix pumpMatrix() const;
31
+    Matrix meanPattern() const;
32
+
27 33
     void addChiSq(float chisq);
28 34
     void addAtomCount(unsigned atomA, unsigned atomP);
29 35
 
... ...
@@ -45,15 +51,87 @@ private:
45 51
     Matrix mAStdMatrix;
46 52
     Matrix mPMeanMatrix;
47 53
     Matrix mPStdMatrix;
48
-    
54
+    Matrix mPumpMatrix;
55
+   
49 56
     std::vector<float> mChisqHistory;
50 57
     std::vector<unsigned> mAtomHistoryA;
51 58
     std::vector<unsigned> mAtomHistoryP;
52 59
 
60
+    PumpThreshold mPumpThreshold;
61
+
53 62
     unsigned mStatUpdates;
54 63
     unsigned mNumPatterns;
64
+    unsigned mPumpUpdates;
55 65
 };
56 66
 
67
+template <class FactorMatrix>
68
+inline void pumpMatrixUniqueThreshold(const FactorMatrix &AMatrix, Matrix *statMatrix)
69
+{
70
+    GAPS_ASSERT(statMatrix->nRow() == AMatrix.nRow());
71
+    GAPS_ASSERT(statMatrix->nCol() == AMatrix.nCol());
72
+
73
+    std::vector<float> maxValues(AMatrix.nRow(), 0.f);
74
+    std::vector<unsigned> maxIndices(AMatrix.nRow(), 0);
75
+    for (unsigned j = 0; j < AMatrix.nCol(); ++j)
76
+    {
77
+        for (unsigned i = 0; i < AMatrix.nRow(); ++i)
78
+        {
79
+            if (maxValues[i] < AMatrix(i,j))
80
+            {
81
+                maxValues[i] = AMatrix(i,j);
82
+                maxIndices[i] = j;
83
+            }
84
+        }
85
+    }
86
+
87
+    for (unsigned i = 0; i < AMatrix.nRow(); ++i)
88
+    {
89
+        statMatrix->operator()(i,maxIndices[i])++;
90
+    }
91
+}
92
+
93
+// implemented the same as unique for now - TODO
94
+template <class FactorMatrix>
95
+inline void pumpMatrixCutThreshold(const FactorMatrix &AMatrix, Matrix *statMatrix)
96
+{
97
+    GAPS_ASSERT(statMatrix->nRow() == AMatrix.nRow());
98
+    GAPS_ASSERT(statMatrix->nCol() == AMatrix.nCol());
99
+
100
+    // we need to access data in columns due to matrix data layout
101
+    std::vector<float> maxValues(AMatrix.nRow(), 0.f);
102
+    std::vector<unsigned> maxIndices(AMatrix.nRow(), 0);
103
+    for (unsigned j = 0; j < AMatrix.nCol(); ++j)
104
+    {
105
+        for (unsigned i = 0; i < AMatrix.nRow(); ++i)
106
+        {
107
+            if (maxValues[i] < AMatrix(i,j))
108
+            {
109
+                maxValues[i] = AMatrix(i,j);
110
+                maxIndices[i] = j;
111
+            }
112
+        }
113
+    }
114
+
115
+    for (unsigned i = 0; i < AMatrix.nRow(); ++i)
116
+    {
117
+        statMatrix->operator()(i,maxIndices[i])++;
118
+    }
119
+}
120
+
121
+template <class Sampler>
122
+void GapsStatistics::updatePump(const Sampler &ASampler)
123
+{
124
+    ++mPumpUpdates;
125
+    if (mPumpThreshold == PUMP_UNIQUE)
126
+    {
127
+        pumpMatrixCutThreshold(ASampler.mMatrix, &mPumpMatrix);
128
+    }
129
+    else
130
+    {
131
+        pumpMatrixUniqueThreshold(ASampler.mMatrix, &mPumpMatrix);
132
+    }
133
+}
134
+
57 135
 template <class Sampler>
58 136
 void GapsStatistics::update(const Sampler &ASampler, const Sampler &PSampler)
59 137
 {
Browse code

store atom and chisq history

Tom Sherman authored on 15/11/2018 21:58:54
Showing1 changed files
... ...
@@ -23,6 +23,12 @@ public:
23 23
     Matrix Pmean() const;
24 24
     Matrix Asd() const;
25 25
     Matrix Psd() const;
26
+
27
+    void addChiSq(float chisq);
28
+    void addAtomCount(unsigned atomA, unsigned atomP);
29
+
30
+    std::vector<float> chisqHistory() const;
31
+    std::vector<unsigned> atomHistory(char m) const;
26 32
     
27 33
     float meanChiSq(const DenseGibbsSampler &PSampler) const;
28 34
     float meanChiSq(const SparseGibbsSampler &PSampler) const;
... ...
@@ -40,6 +46,10 @@ private:
40 46
     Matrix mPMeanMatrix;
41 47
     Matrix mPStdMatrix;
42 48
     
49
+    std::vector<float> mChisqHistory;
50
+    std::vector<unsigned> mAtomHistoryA;
51
+    std::vector<unsigned> mAtomHistoryP;
52
+
43 53
     unsigned mStatUpdates;
44 54
     unsigned mNumPatterns;
45 55
 };
Browse code

updated config to commit file permissions

Tom Sherman authored on 29/10/2018 19:56:14
Showing1 changed files
1 1
old mode 100644
2 2
new mode 100755
Browse code

tests passing

Tom Sherman authored on 25/10/2018 22:31:44
Showing1 changed files
... ...
@@ -12,16 +12,6 @@
12 12
 
13 13
 class GapsStatistics
14 14
 {
15
-private:
16
-
17
-    Matrix mAMeanMatrix;
18
-    Matrix mAStdMatrix;
19
-    Matrix mPMeanMatrix;
20
-    Matrix mPStdMatrix;
21
-    
22
-    unsigned mStatUpdates;
23
-    unsigned mNumPatterns;
24
-
25 15
 public:
26 16
 
27 17
     GapsStatistics(unsigned nGenes, unsigned nSamples, unsigned nPatterns);
... ...
@@ -40,6 +30,18 @@ public:
40 30
     // serialization
41 31
     friend Archive& operator<<(Archive &ar, const GapsStatistics &stat);
42 32
     friend Archive& operator>>(Archive &ar, GapsStatistics &stat);
33
+
34
+#ifndef GAPS_INTERNAL_TESTS
35
+private:
36
+#endif
37
+
38
+    Matrix mAMeanMatrix;
39
+    Matrix mAStdMatrix;
40
+    Matrix mPMeanMatrix;
41
+    Matrix mPStdMatrix;
42
+    
43
+    unsigned mStatUpdates;
44
+    unsigned mNumPatterns;
43 45
 };
44 46
 
45 47
 template <class Sampler>
Browse code

all tests passing besides checkpoints

Tom Sherman authored on 23/10/2018 23:49:14
Showing1 changed files
... ...
@@ -47,31 +47,26 @@ void GapsStatistics::update(const Sampler &ASampler, const Sampler &PSampler)
47 47
 {
48 48
     mStatUpdates++;
49 49
 
50
-    // update     
51
-    // precision loss? use double?
52
-    DEBUG_PING
53 50
     GAPS_ASSERT(mNumPatterns == ASampler.mMatrix.nCol());
54 51
     GAPS_ASSERT(mNumPatterns == PSampler.mMatrix.nCol());
55 52
 
53
+    // precision loss? use double?
56 54
     for (unsigned j = 0; j < mNumPatterns; ++j)
57 55
     {
58 56
         float norm = gaps::max(PSampler.mMatrix.getCol(j));
59 57
         norm = norm == 0.f ? 1.f : norm;
60 58
         GAPS_ASSERT(norm > 0.f);
61 59
 
62
-        DEBUG_PING
63 60
         Vector quot(PSampler.mMatrix.getCol(j) / norm);
64 61
         GAPS_ASSERT(gaps::min(quot) >= 0.f);
65 62
         mPMeanMatrix.getCol(j) += quot;
66 63
         mPStdMatrix.getCol(j) += gaps::elementSq(quot);
67 64
 
68
-        DEBUG_PING
69 65
         Vector prod(ASampler.mMatrix.getCol(j) * norm);
70 66
         GAPS_ASSERT(gaps::min(prod) >= 0.f);
71 67
         mAMeanMatrix.getCol(j) += prod;
72 68
         mAStdMatrix.getCol(j) += gaps::elementSq(prod);
73 69
     }
74
-    DEBUG_PING
75 70
 }
76 71
 
77 72
 #endif
78 73
\ No newline at end of file
Browse code

use free functions instead of class for gaps runner

Tom Sherman authored on 23/10/2018 21:37:09
Showing1 changed files
... ...
@@ -38,7 +38,7 @@ public:
38 38
     float meanChiSq(const SparseGibbsSampler &PSampler) const;
39 39
 
40 40
     // serialization
41
-    friend Archive& operator<<(Archive &ar, GapsStatistics &stat);
41
+    friend Archive& operator<<(Archive &ar, const GapsStatistics &stat);
42 42
     friend Archive& operator>>(Archive &ar, GapsStatistics &stat);
43 43
 };
44 44
 
... ...
@@ -49,22 +49,29 @@ void GapsStatistics::update(const Sampler &ASampler, const Sampler &PSampler)
49 49
 
50 50
     // update     
51 51
     // precision loss? use double?
52
+    DEBUG_PING
53
+    GAPS_ASSERT(mNumPatterns == ASampler.mMatrix.nCol());
54
+    GAPS_ASSERT(mNumPatterns == PSampler.mMatrix.nCol());
55
+
52 56
     for (unsigned j = 0; j < mNumPatterns; ++j)
53 57
     {
54 58
         float norm = gaps::max(PSampler.mMatrix.getCol(j));
55 59
         norm = norm == 0.f ? 1.f : norm;
56 60
         GAPS_ASSERT(norm > 0.f);
57 61
 
62
+        DEBUG_PING
58 63
         Vector quot(PSampler.mMatrix.getCol(j) / norm);
59 64
         GAPS_ASSERT(gaps::min(quot) >= 0.f);
60 65
         mPMeanMatrix.getCol(j) += quot;
61 66
         mPStdMatrix.getCol(j) += gaps::elementSq(quot);
62 67
 
68
+        DEBUG_PING
63 69
         Vector prod(ASampler.mMatrix.getCol(j) * norm);
64 70
         GAPS_ASSERT(gaps::min(prod) >= 0.f);
65 71
         mAMeanMatrix.getCol(j) += prod;
66 72
         mAStdMatrix.getCol(j) += gaps::elementSq(prod);
67 73
     }
74
+    DEBUG_PING
68 75
 }
69 76
 
70 77
 #endif
71 78
\ No newline at end of file
Browse code

restoring some tests

Tom Sherman authored on 22/10/2018 19:25:14
Showing1 changed files
... ...
@@ -48,18 +48,22 @@ void GapsStatistics::update(const Sampler &ASampler, const Sampler &PSampler)
48 48
     mStatUpdates++;
49 49
 
50 50
     // update     
51
+    // precision loss? use double?
51 52
     for (unsigned j = 0; j < mNumPatterns; ++j)
52 53
     {
53 54
         float norm = gaps::max(PSampler.mMatrix.getCol(j));
54 55
         norm = norm == 0.f ? 1.f : norm;
56
+        GAPS_ASSERT(norm > 0.f);
55 57
 
56 58
         Vector quot(PSampler.mMatrix.getCol(j) / norm);
59
+        GAPS_ASSERT(gaps::min(quot) >= 0.f);
57 60
         mPMeanMatrix.getCol(j) += quot;
58 61
         mPStdMatrix.getCol(j) += gaps::elementSq(quot);
59 62
 
60 63
         Vector prod(ASampler.mMatrix.getCol(j) * norm);
64
+        GAPS_ASSERT(gaps::min(prod) >= 0.f);
61 65
         mAMeanMatrix.getCol(j) += prod;
62
-        mAStdMatrix.getCol(j) += gaps::elementSq(prod); // precision loss? use double?
66
+        mAStdMatrix.getCol(j) += gaps::elementSq(prod);
63 67
     }
64 68
 }
65 69
 
Browse code

dense sampler appears to be working

Tom Sherman authored on 15/10/2018 20:52:28
Showing1 changed files
... ...
@@ -1,8 +1,14 @@
1 1
 #ifndef __COGAPS_GAPS_STATISTICS_H__
2 2
 #define __COGAPS_GAPS_STATISTICS_H__
3 3
 
4
-#include "GibbsSampler.h"
4
+#include "math/Math.h"
5
+#include "math/MatrixMath.h"
6
+#include "math/VectorMath.h"
5 7
 #include "data_structures/Matrix.h"
8
+#include "gibbs_sampler/DenseGibbsSampler.h"
9
+#include "gibbs_sampler/SparseGibbsSampler.h"
10
+
11
+#define GAPS_SQ(x) ((x) * (x))
6 12
 
7 13
 class GapsStatistics
8 14
 {
... ...
@@ -18,20 +24,43 @@ private:
18 24
 
19 25
 public:
20 26
 
21
-    GapsStatistics(unsigned nRow, unsigned nCol, unsigned nPatterns);
27
+    GapsStatistics(unsigned nGenes, unsigned nSamples, unsigned nPatterns);
22 28
 
23
-    void update(const GibbsSampler *ASampler, const GibbsSampler *PSampler);
29
+    template <class Sampler>
30
+    void update(const Sampler &ASampler, const Sampler &PSampler);
24 31
 
25 32
     Matrix Amean() const;
26 33
     Matrix Pmean() const;
27 34
     Matrix Asd() const;
28 35
     Matrix Psd() const;
29
-
30
-    float meanChiSq(const GibbsSampler *PSampler) const;
36
+    
37
+    float meanChiSq(const DenseGibbsSampler &PSampler) const;
38
+    float meanChiSq(const SparseGibbsSampler &PSampler) const;
31 39
 
32 40
     // serialization
33 41
     friend Archive& operator<<(Archive &ar, GapsStatistics &stat);
34 42
     friend Archive& operator>>(Archive &ar, GapsStatistics &stat);
35 43
 };
36 44
 
45
+template <class Sampler>
46
+void GapsStatistics::update(const Sampler &ASampler, const Sampler &PSampler)
47
+{
48
+    mStatUpdates++;
49
+
50
+    // update     
51
+    for (unsigned j = 0; j < mNumPatterns; ++j)
52
+    {
53
+        float norm = gaps::max(PSampler.mMatrix.getCol(j));
54
+        norm = norm == 0.f ? 1.f : norm;
55
+
56
+        Vector quot(PSampler.mMatrix.getCol(j) / norm);
57
+        mPMeanMatrix.getCol(j) += quot;
58
+        mPStdMatrix.getCol(j) += gaps::elementSq(quot);
59
+
60
+        Vector prod(ASampler.mMatrix.getCol(j) * norm);
61
+        mAMeanMatrix.getCol(j) += prod;
62
+        mAStdMatrix.getCol(j) += gaps::elementSq(prod); // precision loss? use double?
63
+    }
64
+}
65
+
37 66
 #endif
38 67
\ No newline at end of file
Browse code

polymorphic structure

Tom Sherman authored on 02/10/2018 20:54:21
Showing1 changed files
... ...
@@ -8,10 +8,10 @@ class GapsStatistics
8 8
 {
9 9
 private:
10 10
 
11
-    ColMatrix mAMeanMatrix;
12
-    ColMatrix mAStdMatrix;
13
-    ColMatrix mPMeanMatrix;
14
-    ColMatrix mPStdMatrix;
11
+    Matrix mAMeanMatrix;
12
+    Matrix mAStdMatrix;
13
+    Matrix mPMeanMatrix;
14
+    Matrix mPStdMatrix;
15 15
     
16 16
     unsigned mStatUpdates;
17 17
     unsigned mNumPatterns;
... ...
@@ -20,14 +20,14 @@ public:
20 20
 
21 21
     GapsStatistics(unsigned nRow, unsigned nCol, unsigned nPatterns);
22 22
 
23
-    void update(const GibbsSampler &ASampler, const GibbsSampler &PSampler);
23
+    void update(const GibbsSampler *ASampler, const GibbsSampler *PSampler);
24 24
 
25
-    ColMatrix Amean() const;
26
-    ColMatrix Pmean() const;
27
-    ColMatrix Asd() const;
28
-    ColMatrix Psd() const;
25
+    Matrix Amean() const;
26
+    Matrix Pmean() const;
27
+    Matrix Asd() const;
28
+    Matrix Psd() const;
29 29
 
30
-    float meanChiSq(const GibbsSampler &PSampler) const;
30
+    float meanChiSq(const GibbsSampler *PSampler) const;
31 31
 
32 32
     // serialization
33 33
     friend Archive& operator<<(Archive &ar, GapsStatistics &stat);
Browse code

Single Gibbs Sampler - consistent with old version

Tom Sherman authored on 29/08/2018 15:12:10
Showing1 changed files
... ...
@@ -20,15 +20,14 @@ public:
20 20
 
21 21
     GapsStatistics(unsigned nRow, unsigned nCol, unsigned nPatterns);
22 22
 
23
-    void update(const AmplitudeGibbsSampler &ASampler,
24
-        const PatternGibbsSampler &PSampler);
23
+    void update(const GibbsSampler &ASampler, const GibbsSampler &PSampler);
25 24
 
26 25
     ColMatrix Amean() const;
27 26
     ColMatrix Pmean() const;
28 27
     ColMatrix Asd() const;
29 28
     ColMatrix Psd() const;
30 29
 
31
-    float meanChiSq(const PatternGibbsSampler &ASampler) const;
30
+    float meanChiSq(const GibbsSampler &PSampler) const;
32 31
 
33 32
     // serialization
34 33
     friend Archive& operator<<(Archive &ar, GapsStatistics &stat);
Browse code

nearly consistent with ColMatrix

Tom Sherman authored on 28/08/2018 21:41:41
Showing1 changed files
... ...
@@ -4,50 +4,31 @@
4 4
 #include "GibbsSampler.h"
5 5
 #include "data_structures/Matrix.h"
6 6
 
7
-enum PumpThreshold
8
-{
9
-    PUMP_UNIQUE=1,
10
-    PUMP_CUT=2
11
-};
12
-
13 7
 class GapsStatistics
14 8
 {
15 9
 private:
16 10
 
17 11
     ColMatrix mAMeanMatrix;
18 12
     ColMatrix mAStdMatrix;
19
-    RowMatrix mPMeanMatrix;
20
-    RowMatrix mPStdMatrix;
13
+    ColMatrix mPMeanMatrix;
14
+    ColMatrix mPStdMatrix;
21 15
     
22 16
     unsigned mStatUpdates;
23 17
     unsigned mNumPatterns;
24 18
 
25
-    ColMatrix mPumpMatrix;
26
-    PumpThreshold mPumpThreshold;
27
-    unsigned mPumpStatUpdates;
28
-
29 19
 public:
30 20
 
31
-    GapsStatistics(unsigned nRow, unsigned nCol, unsigned nPatterns,
32
-        PumpThreshold t=PUMP_CUT);
21
+    GapsStatistics(unsigned nRow, unsigned nCol, unsigned nPatterns);
33 22
 
34 23
     void update(const AmplitudeGibbsSampler &ASampler,
35 24
         const PatternGibbsSampler &PSampler);
36 25
 
37 26
     ColMatrix Amean() const;
38
-    RowMatrix Pmean() const;
27
+    ColMatrix Pmean() const;
39 28
     ColMatrix Asd() const;
40
-    RowMatrix Psd() const;
41
-
42
-    float meanChiSq(const AmplitudeGibbsSampler &ASampler) const;
43
-
44
-    // PUMP statistics
45
-    void updatePump(const AmplitudeGibbsSampler &ASampler,
46
-        const PatternGibbsSampler &PSampler);
29
+    ColMatrix Psd() const;
47 30
 
48
-    RowMatrix pumpMatrix() const;
49
-    RowMatrix meanPattern();
50
-    void patternMarkers(ColMatrix normedA, RowMatrix normedP, ColMatrix &statMatrix);
31
+    float meanChiSq(const PatternGibbsSampler &ASampler) const;
51 32
 
52 33
     // serialization
53 34
     friend Archive& operator<<(Archive &ar, GapsStatistics &stat);
Browse code

new R Interface working

Tom Sherman authored on 02/07/2018 00:58:55
Showing1 changed files
... ...
@@ -28,24 +28,25 @@ private:
28 28
 
29 29
 public:
30 30
 
31
-    GapsStatistics(unsigned nRow, unsigned nCol, unsigned nFactor, PumpThreshold t=PUMP_CUT);
32
-
33
-    ColMatrix AMean() const;
34
-    RowMatrix PMean() const;
35
-    ColMatrix AStd() const;
36
-    RowMatrix PStd() const;
37
-
38
-    float meanChiSq(const AmplitudeGibbsSampler &ASampler) const;
31
+    GapsStatistics(unsigned nRow, unsigned nCol, unsigned nPatterns,
32
+        PumpThreshold t=PUMP_CUT);
39 33
 
40 34
     void update(const AmplitudeGibbsSampler &ASampler,
41 35
         const PatternGibbsSampler &PSampler);
42 36
 
37
+    ColMatrix Amean() const;
38
+    RowMatrix Pmean() const;
39
+    ColMatrix Asd() const;
40
+    RowMatrix Psd() const;
41
+
42
+    float meanChiSq(const AmplitudeGibbsSampler &ASampler) const;
43
+
44
+    // PUMP statistics
43 45
     void updatePump(const AmplitudeGibbsSampler &ASampler,
44 46
         const PatternGibbsSampler &PSampler);
45 47
 
46 48
     RowMatrix pumpMatrix() const;
47 49
     RowMatrix meanPattern();
48
-
49 50
     void patternMarkers(ColMatrix normedA, RowMatrix normedP, ColMatrix &statMatrix);
50 51
 
51 52
     // serialization
Browse code

compiling and running after merging with develop

Tom Sherman authored on 18/06/2018 22:20:06
Showing1 changed files
... ...
@@ -43,6 +43,9 @@ public:
43 43
     void updatePump(const AmplitudeGibbsSampler &ASampler,
44 44
         const PatternGibbsSampler &PSampler);
45 45
 
46
+    RowMatrix pumpMatrix() const;
47
+    RowMatrix meanPattern();
48
+
46 49
     void patternMarkers(ColMatrix normedA, RowMatrix normedP, ColMatrix &statMatrix);
47 50
 
48 51
     // serialization
Tom Sherman authored on 18/06/2018 21:11:30
Showing0 changed files
Browse code

framework with dispatcher compiling and running

Tom Sherman authored on 18/06/2018 13:49:06
Showing1 changed files
... ...
@@ -30,10 +30,10 @@ public:
30 30
 
31 31
     GapsStatistics(unsigned nRow, unsigned nCol, unsigned nFactor);
32 32
 
33
-    Rcpp::NumericMatrix AMean() const;
34
-    Rcpp::NumericMatrix AStd() const;
35
-    Rcpp::NumericMatrix PMean() const;
36
-    Rcpp::NumericMatrix PStd() const;
33
+    ColMatrix AMean() const;
34
+    RowMatrix PMean() const;
35
+    ColMatrix AStd() const;
36
+    RowMatrix PStd() const;
37 37
 
38 38
     float meanChiSq(const AmplitudeGibbsSampler &ASampler) const;
39 39
 
Browse code

added back PUMP and FixedPatterns option

Tom Sherman authored on 13/06/2018 21:35:14
Showing1 changed files
... ...
@@ -28,12 +28,14 @@ private:
28 28
 
29 29
 public:
30 30
 
31
-    GapsStatistics(unsigned nRow, unsigned nCol, unsigned nFactor);
31
+    GapsStatistics(unsigned nRow, unsigned nCol, unsigned nFactor, PumpThreshold t=PUMP_CUT);
32 32
 
33 33
     Rcpp::NumericMatrix AMean() const;
34 34
     Rcpp::NumericMatrix AStd() const;
35 35
     Rcpp::NumericMatrix PMean() const;
36 36
     Rcpp::NumericMatrix PStd() const;
37
+    Rcpp::NumericMatrix pumpMatrix() const;
38
+    Rcpp::NumericMatrix meanPattern();
37 39
 
38 40
     float meanChiSq(const AmplitudeGibbsSampler &ASampler) const;
39 41
 
... ...
@@ -43,6 +45,8 @@ public:
43 45
     void updatePump(const AmplitudeGibbsSampler &ASampler,
44 46
         const PatternGibbsSampler &PSampler);
45 47
 
48
+    void patternMarkers(ColMatrix normedA, RowMatrix normedP, ColMatrix &statMatrix);
49
+
46 50
     // serialization
47 51
     friend Archive& operator<<(Archive &ar, GapsStatistics &stat);
48 52
     friend Archive& operator>>(Archive &ar, GapsStatistics &stat);
Browse code

restore serialization

Tom Sherman authored on 11/06/2018 23:46:44
Showing1 changed files
... ...
@@ -4,6 +4,12 @@
4 4
 #include "GibbsSampler.h"
5 5
 #include "data_structures/Matrix.h"
6 6
 
7
+enum PumpThreshold
8
+{
9
+    PUMP_UNIQUE=1,
10
+    PUMP_CUT=2
11
+};
12
+
7 13
 class GapsStatistics
8 14
 {
9 15
 private:
... ...
@@ -16,6 +22,10 @@ private:
16 22
     unsigned mStatUpdates;
17 23
     unsigned mNumPatterns;
18 24
 
25
+    ColMatrix mPumpMatrix;
26
+    PumpThreshold mPumpThreshold;
27
+    unsigned mPumpStatUpdates;
28
+
19 29
 public:
20 30
 
21 31
     GapsStatistics(unsigned nRow, unsigned nCol, unsigned nFactor);
... ...
@@ -29,6 +39,13 @@ public:
29 39
 
30 40
     void update(const AmplitudeGibbsSampler &ASampler,
31 41
         const PatternGibbsSampler &PSampler);
42
+
43
+    void updatePump(const AmplitudeGibbsSampler &ASampler,
44
+        const PatternGibbsSampler &PSampler);
45
+
46
+    // serialization
47
+    friend Archive& operator<<(Archive &ar, GapsStatistics &stat);
48
+    friend Archive& operator>>(Archive &ar, GapsStatistics &stat);
32 49
 };
33 50
 
34 51
 #endif
35 52
\ No newline at end of file
Browse code

added crude time estimate

Tom Sherman authored on 22/05/2018 15:05:30
Showing1 changed files
... ...
@@ -2,7 +2,7 @@
2 2
 #define __COGAPS_GAPS_STATISTICS_H__
3 3
 
4 4
 #include "GibbsSampler.h"
5
-#include "math/Matrix.h"
5
+#include "data_structures/Matrix.h"
6 6
 
7 7
 class GapsStatistics
8 8
 {
Browse code

separated math into own folder

Tom Sherman authored on 03/05/2018 21:05:24
Showing1 changed files
... ...
@@ -2,7 +2,7 @@
2 2
 #define __COGAPS_GAPS_STATISTICS_H__
3 3
 
4 4
 #include "GibbsSampler.h"
5
-#include "Matrix.h"
5
+#include "math/Matrix.h"
6 6
 
7 7
 class GapsStatistics
8 8
 {
Browse code

still inconsistent

sherman5 authored on 17/04/2018 22:28:29
Showing1 changed files
... ...
@@ -25,7 +25,7 @@ public:
25 25
     Rcpp::NumericMatrix PMean() const;
26 26
     Rcpp::NumericMatrix PStd() const;
27 27
 
28
-    float meanChiSq() const;
28
+    float meanChiSq(const AmplitudeGibbsSampler &ASampler) const;
29 29
 
30 30
     void update(const AmplitudeGibbsSampler &ASampler,
31 31
         const PatternGibbsSampler &PSampler);
Browse code

running cleanly, algo still off

sherman5 authored on 17/04/2018 21:37:06
Showing1 changed files
... ...
@@ -1,5 +1,5 @@
1
-#ifndef __GAPS_GAPS_STATISTICS_H__
2
-#define __GAPS_GAPS_STATISTICS_H__
1
+#ifndef __COGAPS_GAPS_STATISTICS_H__
2
+#define __COGAPS_GAPS_STATISTICS_H__
3 3
 
4 4
 #include "GibbsSampler.h"
5 5
 #include "Matrix.h"
Browse code

gibbs sampler interface

sherman5 authored on 19/03/2018 23:23:12
Showing1 changed files
... ...
@@ -13,16 +13,19 @@ private:
13 13
     RowMatrix mPMeanMatrix;
14 14
     RowMatrix mPStdMatrix;
15 15
     
16
-    unsigned mStatUpdates;    
16
+    unsigned mStatUpdates;
17
+    unsigned mNumPatterns;
17 18
 
18 19
 public:
19 20
 
20 21
     GapsStatistics(unsigned nRow, unsigned nCol, unsigned nFactor);
21 22
 
22
-    Rcpp::NumericMatrix AMean();
23
-    Rcpp::NumericMatrix AStd();
24
-    Rcpp::NumericMatrix PMean();
25
-    Rcpp::NumericMatrix PStd();
23
+    Rcpp::NumericMatrix AMean() const;
24
+    Rcpp::NumericMatrix AStd() const;
25
+    Rcpp::NumericMatrix PMean() const;
26
+    Rcpp::NumericMatrix PStd() const;
27
+
28
+    float meanChiSq() const;
26 29
 
27 30
     void update(const AmplitudeGibbsSampler &ASampler,
28 31
         const PatternGibbsSampler &PSampler);
Browse code

separated gibbs samplers

sherman5 authored on 07/03/2018 22:59:17
Showing1 changed files
1 1
new file mode 100644
... ...
@@ -0,0 +1,31 @@
1
+#ifndef __GAPS_GAPS_STATISTICS_H__
2
+#define __GAPS_GAPS_STATISTICS_H__
3
+
4
+#include "GibbsSampler.h"
5
+#include "Matrix.h"
6
+
7
+class GapsStatistics
8
+{
9
+private:
10
+
11
+    ColMatrix mAMeanMatrix;
12
+    ColMatrix mAStdMatrix;
13
+    RowMatrix mPMeanMatrix;
14
+    RowMatrix mPStdMatrix;
15
+    
16
+    unsigned mStatUpdates;    
17
+
18
+public:
19
+
20
+    GapsStatistics(unsigned nRow, unsigned nCol, unsigned nFactor);
21
+
22
+    Rcpp::NumericMatrix AMean();
23
+    Rcpp::NumericMatrix AStd();
24
+    Rcpp::NumericMatrix PMean();
25
+    Rcpp::NumericMatrix PStd();
26
+
27
+    void update(const AmplitudeGibbsSampler &ASampler,
28
+        const PatternGibbsSampler &PSampler);
29
+};
30
+
31
+#endif
0 32
\ No newline at end of file