Browse code

clean up linter warnings

Tom Sherman authored on 24/06/2019 19:44:02
Showing56 changed files

... ...
@@ -2,6 +2,10 @@ language: r
2 2
 cache: packages
3 3
 
4 4
 matrix:
5
+  fast_finish: true
6
+  allow_failures:
7
+    - stage: Linting
8
+    - stage: Code Formatting
5 9
   include:
6 10
     - stage: Bioconductor Release Testing
7 11
       r: bioc-release
... ...
@@ -10,6 +14,12 @@ matrix:
10 14
     - stage: Test Running in Debug Mode
11 15
       install: R -e 'devtools::install_deps(dep = T, configure.args = c(debugMode = " --enable-debug " ))'
12 16
       script: R -f inst/scripts/debugRuns.R
17
+    - stage: Valgrind Analysis
18
+      script: R -d valgrind -e 'library(CoGAPS); data(GIST); CoGAPS(GIST.matrix, nIterations=1000, outputFrequency=250)'
19
+    - stage: Linting
20
+      script: echo "linting not implemented"
21
+    - stage: Code Formatting
22
+      script: echo "formatting not implemented"
13 23
 
14 24
 notifications:
15 25
   slack: fertiglab:LZeNoaS7BnLfqP6XuHiDYdts
... ...
@@ -1,4 +1,9 @@
1
+#include "GapsParameters.h"
2
+#include "GapsResult.h"
1 3
 #include "GapsRunner.h"
4
+#include "data_structures/Matrix.h"
5
+#include "file_parser/FileParser.h"
6
+#include "math/Random.h"
2 7
 #include "utils/GlobalConfig.h"
3 8
 
4 9
 #include <Rcpp.h>
... ...
@@ -1,5 +1,7 @@
1 1
 #include "GapsParameters.h"
2
+#include "utils/Archive.h"
2 3
 #include "utils/GapsPrint.h"
4
+#include "file_parser/FileParser.h"
3 5
 
4 6
 void GapsParameters::print() const
5 7
 {
... ...
@@ -2,11 +2,13 @@
2 2
 #define __COGAPS_GAPS_PARAMETERS_H__
3 3
 
4 4
 #include "data_structures/Matrix.h"
5
-#include "file_parser/FileParser.h"
6 5
 
6
+#include <stdint.h>
7 7
 #include <string>
8 8
 #include <vector>
9 9
 
10
+class Archive;
11
+
10 12
 enum PumpThreshold
11 13
 {
12 14
     PUMP_UNIQUE=1,
... ...
@@ -16,23 +18,17 @@ enum PumpThreshold
16 18
 struct GapsParameters
17 19
 {
18 20
 public:
19
-
20 21
     template <class DataType>
21 22
     explicit GapsParameters(const DataType &data, bool t_transposeData=false,
22 23
         bool t_subsetData=false, bool t_subsetGenes=false,
23 24
         const std::vector<unsigned> &t_dataIndicesSubset=std::vector<unsigned>());
24
-
25 25
     void print() const;
26 26
 
27 27
     Matrix fixedPatterns;
28
-
29 28
     std::vector<unsigned> dataIndicesSubset;
30
-    
31 29
     std::string checkpointFile;
32 30
     std::string checkpointOutFile;
33
-
34 31
     uint32_t seed;
35
-
36 32
     unsigned nGenes;
37 33
     unsigned nSamples;
38 34
     unsigned nPatterns;
... ...
@@ -40,14 +36,11 @@ public:
40 36
     unsigned maxThreads;
41 37
     unsigned outputFrequency;
42 38
     unsigned checkpointInterval;
43
-
44 39
     float alphaA;
45 40
     float alphaP;
46 41
     float maxGibbsMassA;
47 42
     float maxGibbsMassP;
48
-
49 43
     PumpThreshold pumpThreshold;
50
-
51 44
     bool useFixedPatterns;
52 45
     bool subsetData;
53 46
     bool useCheckPoint;
... ...
@@ -59,13 +52,10 @@ public:
59 52
     bool useSparseOptimization;
60 53
     bool takePumpSamples;
61 54
     bool asynchronousUpdates;
62
-
63 55
     char whichMatrixFixed;
64 56
     unsigned workerID;
65 57
     bool runningDistributed;
66
-
67 58
 private:
68
-
69 59
     void calculateDataDimensions(const std::string &file);
70 60
     void calculateDataDimensions(const Matrix &mat);
71 61
 };
... ...
@@ -113,4 +103,4 @@ asynchronousUpdates(true)
113 103
     calculateDataDimensions(data);
114 104
 }
115 105
 
116
-#endif
117 106
\ No newline at end of file
107
+#endif // __COGAPS_GAPS_PARAMETERS_H__
118 108
\ No newline at end of file
... ...
@@ -18,21 +18,7 @@ Psd(stat.Psd()), meanChiSq(0.f), seed(0), chisqHistory(stat.chisqHistory()),
18 18
 atomHistoryA(stat.atomHistory('A')), atomHistoryP(stat.atomHistory('P'))
19 19
 {}
20 20
 
21
-void GapsResult::writeToFile(const std::string &fullPath)
22
-{
23
-    std::size_t pos = fullPath.find_last_of('.');
24
-    std::string base = fullPath.substr(0, pos);
25
-
26
-    switch (FileParser::fileType(fullPath))
27
-    {
28
-        case GAPS_CSV: return writeCsv(base);
29
-        case GAPS_TSV: return writeTsv(base);
30
-        case GAPS_GCT: return writeGct(base);
31
-        default : GAPS_ERROR("Invalid File Type");
32
-    }
33
-}
34
-
35
-void GapsResult::writeCsv(const std::string &path)
21
+void GapsResult::writeToFile(const std::string &path)
36 22
 {
37 23
     unsigned nPatterns = Amean.nCol();
38 24
     std::string label("_" + to_string(nPatterns) + "_");
... ...
@@ -40,24 +26,4 @@ void GapsResult::writeCsv(const std::string &path)
40 26
     FileParser::writeToCsv(path + label + "Pmean.csv", Pmean);
41 27
     FileParser::writeToCsv(path + label + "Asd.csv", Asd);
42 28
     FileParser::writeToCsv(path + label + "Psd.csv", Psd);
43
-}
44
-
45
-void GapsResult::writeTsv(const std::string &path)
46
-{
47
-    unsigned nPatterns = Amean.nCol();
48
-    std::string label("_" + to_string(nPatterns) + "_");
49
-    FileParser::writeToCsv(path + label + "Amean.tsv", Amean);
50
-    FileParser::writeToCsv(path + label + "Pmean.tsv", Pmean);
51
-    FileParser::writeToCsv(path + label + "Asd.tsv", Asd);
52
-    FileParser::writeToCsv(path + label + "Psd.tsv", Psd);
53
-}
54
-
55
-void GapsResult::writeGct(const std::string &path)
56
-{
57
-    unsigned nPatterns = Amean.nCol();
58
-    std::string label("_" + to_string(nPatterns) + "_");
59
-    FileParser::writeToCsv(path + label + "Amean.gct", Amean);
60
-    FileParser::writeToCsv(path + label + "Pmean.gct", Pmean);
61
-    FileParser::writeToCsv(path + label + "Asd.gct", Asd);
62
-    FileParser::writeToCsv(path + label + "Psd.gct", Psd);
63 29
 }
64 30
\ No newline at end of file
... ...
@@ -1,39 +1,34 @@
1 1
 #ifndef __COGAPS_GAPS_RESULT__
2 2
 #define __COGAPS_GAPS_RESULT__
3 3
 
4
-#include "GapsStatistics.h"
5 4
 #include "data_structures/Matrix.h"
6 5
 
7 6
 #include <stdint.h>
8 7
 #include <string>
8
+#include <vector>
9
+
10
+class GapsStatistics;
9 11
 
10 12
 struct GapsResult
11 13
 {
14
+    explicit GapsResult(const GapsStatistics &stat);
15
+    void writeToFile(const std::string &path);
16
+
12 17
     Matrix Amean;
13 18
     Matrix Asd;
14 19
     Matrix Pmean;
15 20
     Matrix Psd;
16 21
     Matrix pumpMatrix;
17 22
     Matrix meanPatternAssignment;
18
-
19 23
     std::vector<float> chisqHistory;
20 24
     std::vector<unsigned> atomHistoryA;
21 25
     std::vector<unsigned> atomHistoryP;
22
-
23 26
     uint64_t totalUpdates;
24 27
     uint32_t seed;
25 28
     unsigned totalRunningTime;
26
-
27 29
     float meanChiSq;
28 30
     float averageQueueLengthA;
29 31
     float averageQueueLengthP;
30
-
31
-    explicit GapsResult(const GapsStatistics &stat);
32
-
33
-    void writeToFile(const std::string &fullPath);
34
-    void writeCsv(const std::string &path);
35
-    void writeTsv(const std::string &path);
36
-    void writeGct(const std::string &path);
37 32
 };
38 33
 
39 34
 #endif // __COGAPS_GAPS_RESULT__
40 35
\ No newline at end of file
... ...
@@ -1,6 +1,10 @@
1 1
 #include "GapsRunner.h"
2
-
2
+#include "GapsResult.h"
3
+#include "GapsParameters.h"
4
+#include "GapsStatistics.h"
5
+#include "math/Random.h"
3 6
 #include "utils/Archive.h"
7
+#include "utils/GlobalConfig.h"
4 8
 #include "gibbs_sampler/AsynchronousGibbsSampler.h"
5 9
 #include "gibbs_sampler/SingleThreadedGibbsSampler.h"
6 10
 #include "gibbs_sampler/DenseNormalModel.h"
... ...
@@ -39,7 +43,7 @@ struct GapsTime
39 43
     unsigned hours;
40 44
     unsigned minutes;
41 45
     unsigned seconds;
42
-    GapsTime(unsigned totalSeconds)
46
+    explicit GapsTime(unsigned totalSeconds)
43 47
     {
44 48
         hours = totalSeconds / 3600;
45 49
         totalSeconds -= hours * 3600;
... ...
@@ -1,10 +1,12 @@
1 1
 #ifndef __COGAPS_GAPS_RUNNER_H__
2 2
 #define __COGAPS_GAPS_RUNNER_H__
3 3
 
4
-#include "GapsResult.h"
5
-#include "GapsParameters.h"
6
-#include "data_structures/Matrix.h"
7
-#include "math/Random.h"
4
+class GapsResult;
5
+class Matrix;
6
+class GapsParameters;
7
+class GapsRandomState;
8
+
9
+#include <string>
8 10
 
9 11
 // these two functions are the top-level functions exposed to the C++
10 12
 // code that is being wrapped by any given language
... ...
@@ -1,4 +1,5 @@
1 1
 #include "GapsStatistics.h"
2
+#include "utils/Archive.h"
2 3
 #include "math/Math.h"
3 4
 
4 5
 GapsStatistics::GapsStatistics(unsigned nGenes, unsigned nSamples, unsigned nPatterns)
... ...
@@ -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
 
... ...
@@ -1,4 +1,25 @@
1 1
 #include "Atom.h"
2
+#include "../utils/Archive.h"
3
+
4
+#include <cstddef>
5
+
6
+AtomNeighborhood::AtomNeighborhood()
7
+    : center(NULL), left(NULL), right(NULL)
8
+{}
9
+
10
+AtomNeighborhood::AtomNeighborhood(Atom *l, Atom *c, Atom *r)
11
+    : center(c), left(l), right(r)
12
+{}
13
+
14
+bool AtomNeighborhood::hasLeft() const
15
+{
16
+    return left != NULL;
17
+}
18
+
19
+bool AtomNeighborhood::hasRight() const
20
+{
21
+    return right != NULL;
22
+}
2 23
 
3 24
 Atom::Atom(uint64_t p, float m)
4 25
     : mPos(p), mIterator(), mLeftIndex(-1), mRightIndex(-1), mIndex(-1), mMass(m)
... ...
@@ -1,9 +1,8 @@
1 1
 #ifndef __COGAPS_ATOM_H__
2 2
 #define __COGAPS_ATOM_H__
3 3
 
4
-#include "../utils/Archive.h"
5
-
6 4
 struct Atom;
5
+class Archive;
7 6
 class AtomicDomain;
8 7
 
9 8
 // this is the map used internally by the atomic domain
... ...
@@ -12,44 +11,34 @@ typedef MutableMap<uint64_t, unsigned> AtomMapType;
12 11
 
13 12
 struct AtomNeighborhood
14 13
 {
14
+    AtomNeighborhood();
15
+    AtomNeighborhood(Atom *l, Atom *c, Atom *r);
16
+    bool hasLeft() const;
17
+    bool hasRight() const;
18
+
15 19
     Atom* center;
16 20
     Atom* left;
17 21
     Atom* right;
18
-
19
-    AtomNeighborhood() : center(NULL), left(NULL), right(NULL) {}
20
-    AtomNeighborhood(Atom *l, Atom *c, Atom *r)
21
-        : center(c), left(l), right(r)
22
-    {}
23
-
24
-    bool hasLeft() const { return left != NULL; }
25
-    bool hasRight() const { return right != NULL; }
26 22
 };
27 23
 
28 24
 struct Atom
29 25
 {
30 26
 public:
31
-
32 27
     Atom(uint64_t p, float m);
33
-
34 28
     uint64_t pos() const;
35 29
     float mass() const;
36 30
     void updateMass(float newMass);
37
-
38 31
     friend Archive& operator<<(Archive& ar, const Atom &a);
39 32
     friend Archive& operator>>(Archive& ar, Atom &a);
40
-
41 33
 private:
42
-
43 34
     // only the atomic domain can change the position of an atom, since it is
44 35
     // responsible for keeping them ordered
45 36
     friend class AtomicDomain;
46 37
     void updatePos(uint64_t newPos);
47
-
48 38
     void setLeftIndex(int index);
49 39
     void setRightIndex(int index);
50 40
     void setIndex(int index);
51 41
     void setIterator(AtomMapType::iterator it);    
52
-
53 42
     bool hasLeft() const;
54 43
     bool hasRight() const;
55 44
     int leftIndex() const;
... ...
@@ -57,15 +46,12 @@ private:
57 46
     int index() const;
58 47
     AtomMapType::iterator iterator() const;
59 48
 
60
-    uint64_t mPos;
61 49
     AtomMapType::iterator mIterator; // iterator to position in map
62
-
63
-    // storing the index allows vector lookup once found in map
64
-    int mLeftIndex;
65
-    int mRightIndex;
66
-    int mIndex;
67
-
68
-    float mMass;
50
+    uint64_t mPos; // position of the atom
51
+    int mLeftIndex; // index of left neighbor
52
+    int mRightIndex; // index of right neighbor
53
+    int mIndex; // storing the index allows vector lookup once found in map
54
+    float mMass; // mass of the atom
69 55
 };
70 56
 
71 57
 #endif // __COGAPS_ATOM_H__
72 58
\ No newline at end of file
... ...
@@ -1,5 +1,7 @@
1 1
 #include "AtomicDomain.h"
2 2
 #include "../utils/GapsAssert.h"
3
+#include "../utils/Archive.h"
4
+#include "../math/Random.h"
3 5
 
4 6
 #include <algorithm>
5 7
 #include <limits>
... ...
@@ -39,7 +41,7 @@ AtomNeighborhood AtomicDomain::randomAtomWithNeighbors(GapsRng *rng)
39 41
 uint64_t AtomicDomain::randomFreePosition(GapsRng *rng) const
40 42
 {
41 43
     uint64_t pos = rng->uniform64(1, mDomainLength);
42
-    while (mAtomMap.count(pos))
44
+    while (mAtomMap.count(pos) != 0u)
43 45
     {
44 46
         pos = rng->uniform64(1, mDomainLength);
45 47
     }
... ...
@@ -97,9 +99,13 @@ void AtomicDomain::erase(Atom *atom)
97 99
         mAtoms[index].setIndex(index);
98 100
         mAtoms[index].iterator()->second = index;
99 101
         if (leftIndex >= 0)
102
+        {
100 103
             mAtoms[leftIndex].setRightIndex(index);
104
+        }
101 105
         if (rightIndex >= 0)
106
+        {
102 107
             mAtoms[rightIndex].setLeftIndex(index);
108
+        }
103 109
     }
104 110
     mAtoms.pop_back();
105 111
 }
... ...
@@ -1,54 +1,44 @@
1 1
 #ifndef __COGAPS_ATOMIC_DOMAIN_H__
2 2
 #define __COGAPS_ATOMIC_DOMAIN_H__
3 3
 
4
-#include "../data_structures/HashSets.h"
5
-#include "../math/Random.h"
6
-#include "../utils/Archive.h"
7 4
 #include "Atom.h"
8 5
 
9 6
 #include <vector>
7
+#include <cstddef>
10 8
 
11
-// needed for friend declarations
12 9
 template <class StoragePolicy>
13 10
 class SingleThreadedGibbsSampler;
14 11
 
12
+class Archive;
13
+class GapsRng;
14
+
15 15
 class AtomicDomain
16 16
 {
17 17
 public:
18
-
19
-    AtomicDomain(uint64_t nBins);
20
-
21
-    // access atoms
18
+    explicit AtomicDomain(uint64_t nBins);
22 19
     Atom* front();
23 20
     Atom* randomAtom(GapsRng *rng);
24 21
     AtomNeighborhood randomAtomWithNeighbors(GapsRng *rng);
25
-
26 22
     uint64_t randomFreePosition(GapsRng *rng) const;
27 23
     uint64_t size() const;
28
-
29
-    // serialization
30 24
     friend Archive& operator<<(Archive &ar, const AtomicDomain &domain);
31 25
     friend Archive& operator>>(Archive &ar, AtomicDomain &domain);
32
-   
33 26
 private:
34
-
35
-    template <class StoragePolicy>
36
-    friend class SingleThreadedGibbsSampler;
37
-
38 27
     // both insert and erase will invalidate pointers to atoms and neither
39 28
     // of these functions should be considered thread-safe in any way. This
40 29
     // class can provide increased performance when thread-safety and 
41 30
     // concurrency is not needed. Friend classes are declared here as a way
42 31
     // to enforce this contract between classes.
32
+    template <class StoragePolicy>
33
+    friend class SingleThreadedGibbsSampler;
34
+
43 35
     Atom* insert(uint64_t pos, float mass);
44 36
     void erase(Atom *atom);
45 37
     void move(Atom *atom, uint64_t newPos);
46 38
 
47 39
     AtomMapType mAtomMap; // sorted, used when inserting atoms to find neighbors
48 40
     std::vector<Atom> mAtoms; // unsorted, used for reads
49
-
50
-    // size of atomic domain to ensure all bins are equal length
51
-    uint64_t mDomainLength;
41
+    uint64_t mDomainLength; // size of atomic domain to ensure all bins are equal length
52 42
 };
53 43
 
54 44
 #endif // __COGAPS_ATOMIC_DOMAIN_H__
... ...
@@ -1,4 +1,25 @@
1 1
 #include "ConcurrentAtom.h"
2
+#include "../utils/Archive.h"
3
+
4
+#include <cstddef>
5
+
6
+ConcurrentAtomNeighborhood::ConcurrentAtomNeighborhood()
7
+    : center(NULL), left(NULL), right(NULL)
8
+{}
9
+
10
+ConcurrentAtomNeighborhood::ConcurrentAtomNeighborhood(ConcurrentAtom *l, ConcurrentAtom *c, ConcurrentAtom *r)
11
+    : center(c), left(l), right(r)
12
+{}
13
+    
14
+bool ConcurrentAtomNeighborhood::hasLeft() const
15
+{
16
+    return left != NULL;
17
+}
18
+
19
+bool ConcurrentAtomNeighborhood::hasRight() const
20
+{
21
+    return right != NULL;
22
+}
2 23
 
3 24
 ConcurrentAtom::ConcurrentAtom(uint64_t p, float m)
4 25
     : mPos(p), mLeft(NULL), mRight(NULL), mIterator(), mIndex(0), mMass(m)
... ...
@@ -44,14 +65,14 @@ void ConcurrentAtom::updatePos(uint64_t newPos)
44 65
     mPos = newPos;
45 66
 }
46 67
 
47
-void ConcurrentAtom::setLeft(ConcurrentAtom* ConcurrentAtom)
68
+void ConcurrentAtom::setLeft(ConcurrentAtom* atom)
48 69
 {
49
-    mLeft = ConcurrentAtom;
70
+    mLeft = atom;
50 71
 }
51 72
 
52
-void ConcurrentAtom::setRight(ConcurrentAtom* ConcurrentAtom)
73
+void ConcurrentAtom::setRight(ConcurrentAtom* atom)
53 74
 {
54
-    mRight = ConcurrentAtom;
75
+    mRight = atom;
55 76
 }
56 77
 
57 78
 void ConcurrentAtom::setIndex(unsigned index)
... ...
@@ -1,10 +1,9 @@
1 1
 #ifndef __COGAPS_CONCURRENT_ATOM_H__
2 2
 #define __COGAPS_CONCURRENT_ATOM_H__
3 3
 
4
-#include "../utils/Archive.h"
5
-
6 4
 struct ConcurrentAtom;
7 5
 class ConcurrentAtomicDomain;
6
+class Archive;
8 7
 
9 8
 // this is the map used internally by the atomic domain
10 9
 #include "../data_structures/MutableMap.h"
... ...
@@ -12,44 +11,33 @@ typedef MutableMap<uint64_t, ConcurrentAtom*> ConcurrentAtomMapType;
12 11
 
13 12
 struct ConcurrentAtomNeighborhood
14 13
 {
14
+    ConcurrentAtomNeighborhood();
15
+    ConcurrentAtomNeighborhood(ConcurrentAtom *l, ConcurrentAtom *c, ConcurrentAtom *r);
16
+    bool hasLeft() const;
17
+    bool hasRight() const;
15 18
     ConcurrentAtom *center;
16 19
     ConcurrentAtom *left;
17 20
     ConcurrentAtom *right;
18
-
19
-    ConcurrentAtomNeighborhood() : center(NULL), left(NULL), right(NULL) {}
20
-    ConcurrentAtomNeighborhood(ConcurrentAtom *l, ConcurrentAtom *c, ConcurrentAtom *r)
21
-        : center(c), left(l), right(r)
22
-    {}
23
-
24
-    bool hasLeft() const { return left != NULL; }
25
-    bool hasRight() const { return right != NULL; }
26 21
 };
27 22
 
28 23
 struct ConcurrentAtom
29 24
 {
30 25
 public:
31
-
32 26
     ConcurrentAtom(uint64_t p, float m);
33
-
34 27
     uint64_t pos() const;
35 28
     float mass() const;
36 29
     void updateMass(float newMass);
37
-
38 30
     friend Archive& operator<<(Archive& ar, const ConcurrentAtom &a);
39 31
     friend Archive& operator>>(Archive& ar, ConcurrentAtom &a);
40
-
41
-//private:
42
-
32
+//private: // TODO
43 33
     // only the atomic domain can change the position of an atom, since it is
44 34
     // responsible for keeping them ordered
45 35
     friend class ConcurrentAtomicDomain;
46 36
     void updatePos(uint64_t newPos);
47
-
48 37
     void setLeft(ConcurrentAtom *atom);
49 38
     void setRight(ConcurrentAtom *atom);
50 39
     void setIndex(unsigned index);
51 40
     void setIterator(ConcurrentAtomMapType::iterator it);    
52
-
53 41
     bool hasLeft() const;
54 42
     bool hasRight() const;
55 43
     ConcurrentAtom* left() const;
... ...
@@ -58,12 +46,9 @@ public:
58 46
     ConcurrentAtomMapType::iterator iterator() const;
59 47
 
60 48
     uint64_t mPos;
61
-
62 49
     ConcurrentAtom *mLeft;
63 50
     ConcurrentAtom *mRight;
64
-
65 51
     ConcurrentAtomMapType::iterator mIterator; // iterator to position in map
66
-
67 52
     unsigned mIndex; // storing the index allows vector lookup once found in map
68 53
     float mMass;
69 54
 };
... ...
@@ -1,4 +1,6 @@
1 1
 #include "ConcurrentAtomicDomain.h"
2
+#include "../math/Random.h"
3
+#include "../utils/Archive.h"
2 4
 #include "../utils/GapsAssert.h"
3 5
 
4 6
 #include <algorithm>
... ...
@@ -1,47 +1,36 @@
1 1
 #ifndef __COGAPS_CONCURRENT_ATOMIC_DOMAIN_H__
2 2
 #define __COGAPS_CONCURRENT_ATOMIC_DOMAIN_H__
3 3
 
4
-#include "../data_structures/HashSets.h"
5
-#include "../math/Random.h"
6
-#include "../utils/Archive.h"
7 4
 #include "ConcurrentAtom.h"
8 5
 
9 6
 #include <vector>
10 7
 
11
-// needed for friend declarations
12
-template <class StoragePolicy>
13
-class AsynchronousGibbsSampler;
14 8
 template <class StoragePolicy>
15 9
 class SingleThreadedGibbsSampler;
10
+
11
+struct ConcurrentAtom;
12
+class Archive;
13
+class GapsRng;
16 14
 class ProposalQueue;
17 15
 
18 16
 class ConcurrentAtomicDomain
19 17
 {
20 18
 public:
21
-
22
-    ConcurrentAtomicDomain(uint64_t nBins);
23
-
24
-    // access atoms
19
+    explicit ConcurrentAtomicDomain(uint64_t nBins);
25 20
     ConcurrentAtom* front();
26 21
     const ConcurrentAtom* front() const;
27 22
     ConcurrentAtom* randomAtom(GapsRng *rng);
28 23
     ConcurrentAtomNeighborhood randomAtomWithNeighbors(GapsRng *rng);
29
-
30 24
     uint64_t randomFreePosition(GapsRng *rng) const;
31 25
     uint64_t size() const;
32
-
33 26
     void cacheErase(ConcurrentAtom *atom); // OpenMP thread safe
34 27
     void move(ConcurrentAtom *atom, uint64_t newPos); // OpenMP thread safe
35 28
     void flushEraseCache();
36
-
37
-    // serialization
38 29
     friend Archive& operator<<(Archive &ar, const ConcurrentAtomicDomain &domain);
39 30
     friend Archive& operator>>(Archive &ar, ConcurrentAtomicDomain &domain);
40
-
41 31
 #ifdef GAPS_DEBUG
42 32
     bool isSorted() const;
43 33
 #endif
44
-   
45 34
 private:
46 35
 
47 36
     // only these classes can call the non-thread safe insert and erase functions
... ...
@@ -54,11 +43,9 @@ private:
54 43
     void erase(ConcurrentAtom *atom);
55 44
 
56 45
     ConcurrentAtomMapType mAtomMap; // sorted, used when inserting atoms to find neighbors
57
-    std::vector<ConcurrentAtom*> mAtoms; // unsorted, used for reads
46
+    std::vector<ConcurrentAtom*> mAtoms; // unsorted, used for random selection of atoms
58 47
     std::vector<ConcurrentAtom*> mEraseCache;
59
-
60
-    // size of atomic domain to ensure all bins are equal length
61
-    uint64_t mDomainLength;
48
+    uint64_t mDomainLength; // size of atomic domain to ensure all bins are equal length
62 49
 };
63 50
 
64 51
 #endif // __COGAPS_CONCURRENT_ATOMIC_DOMAIN_H__
... ...
@@ -1,6 +1,8 @@
1 1
 #include "ProposalQueue.h"
2
+#include "../atomic/ConcurrentAtomicDomain.h"
2 3
 #include "../math/Math.h"
3 4
 #include "../math/Random.h"
5
+#include "../utils/Archive.h"
4 6
 #include "../utils/GapsAssert.h"
5 7
 
6 8
 #include <limits>
... ...
@@ -1,95 +1,75 @@
1 1
 #ifndef __COGAPS_PROPOSAL_QUEUE_H__
2 2
 #define __COGAPS_PROPOSAL_QUEUE_H__
3 3
 
4
-#include "ConcurrentAtomicDomain.h"
5
-#include "../data_structures/HashSets.h"
6 4
 #include "../math/Random.h"
7
-#include "../utils/Archive.h"
5
+#include "../data_structures/HashSets.h"
8 6
 
9 7
 #include <cstddef>
10 8
 #include <stdint.h>
11 9
 #include <vector>
12 10
 
11
+struct ConcurrentAtom;
12
+class Archive;
13
+class ConcurrentAtomicDomain;
14
+
13 15
 struct AtomicProposal
14 16
 {
17
+    AtomicProposal(char t, GapsRandomState *randState);
18
+
15 19
     mutable GapsRng rng; // used for consistency no matter number of threads 
16
- 
17 20
     uint64_t pos; // used for move
18 21
     ConcurrentAtom *atom1; // used for birth/death/move/exchange
19 22
     ConcurrentAtom *atom2; // used for exchange
20
-
21 23
     uint32_t r1; // row of atom1
22 24
     uint32_t c1; // col of atom1
23 25
     uint32_t r2; // row of pos (move) or atom2 (exchange)
24 26
     uint32_t c2; // col of pos (move) or atom2 (exchange)
25
-
26 27
     char type; // birth (B), death (D), move (M), exchange (E)
27
-
28
-    AtomicProposal(char t, GapsRandomState *randState);
29 28
 };
30 29
 
31 30
 class ProposalQueue
32 31
 {
33 32
 public:
34
-
35
-    // initialize
36 33
     ProposalQueue(uint64_t nElements, uint64_t nPatterns, GapsRandomState *randState);
37 34
     void setAlpha(float alpha);
38 35
     void setLambda(float lambda);
39
-
40
-    // modify/access queue
41 36
     void populate(ConcurrentAtomicDomain &domain, unsigned limit);
42 37
     void clear();
43 38
     unsigned size() const;
44 39
     AtomicProposal& operator[](int n);
45
-
46
-    // update min/max atoms
47 40
     void acceptDeath();
48 41
     void rejectDeath();
49 42
     void acceptBirth();
50 43
     void rejectBirth();
51
-
52 44
     unsigned nProcessed() const;
53
-
54
-    // serialization
55 45
     friend Archive& operator<<(Archive &ar, const ProposalQueue &queue);
56 46
     friend Archive& operator>>(Archive &ar, ProposalQueue &queue);
57
-
58 47
 private:
48
+    float deathProb(double nAtoms) const;
49
+    bool makeProposal(ConcurrentAtomicDomain &domain);
50
+    bool birth(ConcurrentAtomicDomain &domain);
51
+    bool death(ConcurrentAtomicDomain &domain);
52
+    bool move(ConcurrentAtomicDomain &domain);
53
+    bool exchange(ConcurrentAtomicDomain &domain);
59 54
 
60 55
     std::vector<AtomicProposal> mQueue; // not really a queue for now
61
-    
62 56
     FixedHashSetU32 mUsedMatrixIndices;
63 57
     SmallHashSetU64 mUsedAtoms;
64 58
     SmallPairedHashSetU64 mProposedMoves;
65
-
66 59
     GapsRandomState *mRandState;
67 60
     mutable GapsRng mRng;
68
-
69 61
     uint64_t mMinAtoms;
70 62
     uint64_t mMaxAtoms;
71 63
     uint64_t mBinLength; // length of single bin
72 64
     uint64_t mNumCols;
73
-
74 65
     double mAlpha;
75 66
     double mDomainLength; // length of entire atomic domain
76 67
     double mNumBins; // number of matrix elements
77
-
78 68
     float mLambda;
79 69
     float mU1;
80 70
     float mU2;
81
-    
82 71
     unsigned mNumProcessed;
83
-
84 72
     bool mUseCachedRng;
85
-
86
-    float deathProb(double nAtoms) const;
87
-
88
-    bool makeProposal(ConcurrentAtomicDomain &domain);
89
-    bool birth(ConcurrentAtomicDomain &domain);
90
-    bool death(ConcurrentAtomicDomain &domain);
91
-    bool move(ConcurrentAtomicDomain &domain);
92
-    bool exchange(ConcurrentAtomicDomain &domain);
93 73
 };
94 74
 
95 75
 #endif // __COGAPS_PROPOSAL_QUEUE_H__
96 76
\ No newline at end of file
... ...
@@ -1,15 +1,9 @@
1 1
 #ifndef __COGAPS_HASH_SETS_H__
2 2
 #define __COGAPS_HASH_SETS_H__
3 3
 
4
-#include "../utils/GlobalConfig.h"
5
-
6 4
 #include <stdint.h>
7 5
 #include <vector>
8 6
 
9
-#ifdef __GAPS_OPENMP__ // defined in global config
10
-#include <omp.h>
11
-#endif
12
-
13 7
 class FixedHashSetU32
14 8
 {
15 9
 public:
... ...
@@ -37,9 +31,9 @@ private:
37 31
 
38 32
 struct PositionPair
39 33
 {
34
+    PositionPair(uint64_t inA, uint64_t inB) : a(inA), b(inB) {}
40 35
     uint64_t a;
41 36
     uint64_t b;
42
-    PositionPair(uint64_t inA, uint64_t inB) : a(inA), b(inB) {}
43 37
 };
44 38
 
45 39
 class SmallPairedHashSetU64
... ...
@@ -1,4 +1,8 @@
1 1
 #include "HybridMatrix.h"
2
+#include "HybridVector.h"
3
+#include "Matrix.h"
4
+#include "../utils/Archive.h"
5
+#include "Vector.h"
2 6
 
3 7
 HybridMatrix::HybridMatrix(unsigned nrow, unsigned ncol)
4 8
     :
... ...
@@ -1,14 +1,14 @@
1 1
 #ifndef __COGAPS_HYBRID_MATRIX_H__
2 2
 #define __COGAPS_HYBRID_MATRIX_H__
3 3
 
4
-#include "HybridVector.h"
5
-#include "Matrix.h"
6 4
 #include "Vector.h"
7
-
8
-#include "../utils/GapsAssert.h"
5
+#include "HybridVector.h"
9 6
 
10 7
 #include <vector>
11 8
 
9
+class Archive;
10
+class Matrix;
11
+
12 12
 // Stores data in row major order with dense vectors, and column major order
13 13
 // with hybrid vectors - which provide a sparse vector interface but have
14 14
 // a dense vector implementation. This data structure provides row-wise
... ...
@@ -18,26 +18,18 @@
18 18
 class HybridMatrix
19 19
 {
20 20
 public:
21
-
22 21
     HybridMatrix(unsigned nrow, unsigned ncol);
23
-    
24 22
     unsigned nRow() const;
25 23
     unsigned nCol() const;
26
-
27 24
     void add(unsigned i, unsigned j, float v);
28 25
     void set(unsigned i, unsigned j, float v);
29 26
     float operator()(unsigned i, unsigned j) const;
30
-    
31 27
     const Vector& getRow(unsigned n) const;
32 28
     const HybridVector& getCol(unsigned n) const;
33
-
34 29
     void operator=(const Matrix &mat);
35
-
36 30
     friend Archive& operator<<(Archive &ar, const HybridMatrix &vec);
37 31
     friend Archive& operator>>(Archive &ar, HybridMatrix &vec);
38
-
39 32
 private:
40
-
41 33
     std::vector<Vector> mRows;
42 34
     std::vector<HybridVector> mCols;
43 35
     unsigned mNumRows;
... ...
@@ -1,6 +1,8 @@
1 1
 #include "HybridVector.h"
2 2
 #include "../math/Math.h"
3 3
 #include "../math/SIMD.h"
4
+#include "../utils/Archive.h"
5
+#include "../utils/GapsAssert.h"
4 6
 
5 7
 #define SIMD_PAD(x) (gaps::simd::Index::increment() + \
6 8
     gaps::simd::Index::increment() * ((x) / gaps::simd::Index::increment())) 
... ...
@@ -49,6 +51,7 @@ unsigned HybridVector::size() const
49 51
     return mSize;
50 52
 }
51 53
 
54
+// can be called from multiple concurrent OpenMP threads
52 55
 bool HybridVector::add(unsigned i, float v)
53 56
 {
54 57
     GAPS_ASSERT(i < mSize);
... ...
@@ -65,6 +68,7 @@ bool HybridVector::add(unsigned i, float v)
65 68
     return false;
66 69
 }
67 70
 
71
+// can be called from multiple concurrent OpenMP threads
68 72
 bool HybridVector::set(unsigned i, float v)
69 73
 {
70 74
     GAPS_ASSERT(i < mSize);
... ...
@@ -1,49 +1,38 @@
1 1
 #ifndef __COGAPS_HYBRID_VECTOR_H__
2 2
 #define __COGAPS_HYBRID_VECTOR_H__
3 3
 
4
-#include "../utils/Archive.h"
5
-
6
-#include "boost/align/aligned_allocator.hpp"
7 4
 #include <vector>
5
+#include <stdint.h>
6
+#pragma GCC diagnostic push
7
+#pragma GCC diagnostic ignored "-Wdeprecated-declarations"
8
+#include <boost/align/aligned_allocator.hpp>
9
+#pragma GCC diagnostic pop
8 10
 
9 11
 // need to align data for SIMD
10 12
 namespace bal = boost::alignment;
11 13
 typedef std::vector<float, bal::aligned_allocator<float,32> > aligned_vector;
12 14
 
13
-class SparseIteratorTwo;
14
-class SparseIteratorThree;
15
+class Archive;
15 16
 
16 17
 // stored as a dense vector (efficient setting of values) but maintains
17 18
 // index bit flags of non-zeros so it can be used with SparseIterator
18 19
 class HybridVector
19 20
 {
20 21
 public:
21
-
22
-    friend class SparseIteratorTwo;
23
-    friend class SparseIteratorThree;
24
-    
25
-    template <unsigned N>
26
-    friend class SparseIterator;
27
-
28 22
     explicit HybridVector(unsigned sz);
29 23
     explicit HybridVector(const std::vector<float> &v);
30
-
31
-    bool empty() const;
32
-    unsigned size() const;
33
-
34 24
     bool add(unsigned i, float v); // true if zeros out data
35 25
     bool set(unsigned i, float v); // true if zeros out data
36
-    float operator[](unsigned i) const;
37
-
26
+    bool empty() const;
27
+    unsigned size() const;
38 28
     const float* ptr() const;
39
-
29
+    float operator[](unsigned i) const;
40 30
     const std::vector<uint64_t>& getBitFlags() const { return mIndexBitFlags; }
41
-
42 31
     friend Archive& operator<<(Archive &ar, const HybridVector &vec);
43 32
     friend Archive& operator>>(Archive &ar, HybridVector &vec);
44
-
45 33
 private:
46
-
34
+    template <unsigned N>
35
+    friend class SparseIterator;
47 36
     std::vector<uint64_t> mIndexBitFlags;
48 37
     aligned_vector mData;
49 38
     unsigned mSize;
... ...
@@ -1,11 +1,13 @@
1
-#include "SparseVector.h"
2 1
 #include "Matrix.h"
2
+#include "SparseVector.h"
3 3
 #include "../file_parser/FileParser.h"
4 4
 #include "../file_parser/MatrixElement.h"
5
+#include "../utils/Archive.h"
5 6
 #include "../utils/GapsAssert.h"
7
+#include "Vector.h"
6 8
 
7
-#include <iterator>
8 9
 #include <algorithm>
10
+#include <iterator>
9 11
 
10 12
 Matrix::Matrix() : mNumRows(0), mNumCols(0) {}
11 13
 
... ...
@@ -1,11 +1,13 @@
1 1
 #ifndef __COGAPS_MATRIX_H__
2 2
 #define __COGAPS_MATRIX_H__
3 3
 
4
-#include "../utils/Archive.h"
5 4
 #include "Vector.h"
6 5
 
6
+#include <string>
7 7
 #include <vector>
8 8
 
9
+class Archive;
10
+
9 11
 class Matrix
10 12
 {
11 13
 public:
... ...
@@ -16,25 +18,17 @@ public:
16 18
         std::vector<unsigned> indices);
17 19
     Matrix(const std::string &path, bool genesInCols, bool subsetGenes,
18 20
         std::vector<unsigned> indices);
19
-
20 21
     unsigned nRow() const;
21 22
     unsigned nCol() const;
22
-
23 23
     void pad(float val);
24
-
25 24
     float operator()(unsigned i, unsigned j) const;
26 25
     float& operator()(unsigned i, unsigned j);
27
-
28 26
     Vector& getCol(unsigned col);
29 27
     const Vector& getCol(unsigned col) const;
30
-    
31 28
     bool empty() const;
32
-
33 29
     friend Archive& operator<<(Archive &ar, const Matrix &mat);
34 30
     friend Archive& operator>>(Archive &ar, Matrix &mat);
35
-
36 31
 private:
37
-
38 32
     std::vector<Vector> mCols;
39 33
     unsigned mNumRows;
40 34
     unsigned mNumCols;
... ...
@@ -2,24 +2,13 @@
2 2
 #define __COGAPS_MUTABLE_MAP_H__
3 3
 
4 4
 #include <stdint.h>
5
-
6
-#define ACTIVE_VERSION 1
7
-
8
-// VERSION 0: std::map
9
-// VERSION 1: cast away const to modify the key directly in std::map
10
-// VERSION 2: wrap keys in struct with mutable members
11
-// VERSION 3: boost multi-index
12
-// VERSION 4: custom implementation
5
+#include <map>
13 6
 
14 7
 // This map should provide the same functionality as std::map, with one additional
15 8
 // feature: it allows the key of an element to be changed, and if that change does
16 9
 // not re-order the map, then the change should happen in O(1), otherwise the change
17 10
 // should happen in O(logN) where N is the number of elements in the map
18 11
 
19
-#if ACTIVE_VERSION == 0
20
-
21
-#include <map>
22
-
23 12
 template <class K, class V>
24 13
 class MutableMap
25 14
 {
... ...
@@ -30,7 +19,6 @@ public:
30 19
     class iterator
31 20
     {
32 21
     public:
33
-    
34 22
         iterator() : mIt() {}
35 23
         iterator(const iterator& it) : mIt(it.mIt) {}
36 24
         ~iterator() {}
... ...
@@ -41,95 +29,15 @@ public:
41 29
         std::pair<const K, V>* operator->() { return &(operator *()); }
42 30
         bool operator==(const iterator &it) { return mIt == it.mIt; }
43 31
         bool operator!=(const iterator &it) { return !(*this == it); }
44
-
45 32
     private:
46
-
47 33
         friend class MutableMap;
48
-        iterator(typename std::map<K, V>::iterator it) : mIt(it) {}
49
-        typename std::map<K, V>::iterator mIt;
50
-    };
51
-
52
-    unsigned count(const K &key) const
53
-    {
54
-        return mMap.count(key);
55
-    }
56
-    
57
-    std::pair<iterator, bool> insert(const std::pair<K, V> &val)
58
-    {
59
-        std::pair< typename std::map<K, V>::iterator, bool> result = mMap.insert(val);
60
-        iterator it(result.first);
61
-        return std::pair<iterator, bool>(it, result.second);
62
-    }
63
-
64
-    void erase(iterator it)
65
-    {
66
-        mMap.erase(it.mIt);
67
-    }
68
-
69
-    iterator find(const K &key)
70
-    {
71
-        return mMap.find(key);
72
-    }
73
-
74
-    void updateKey(iterator it, const K &newKey)
75
-    {
76
-        V val = (*it).second;
77
-        mMap.erase(it.mIt);
78
-        insert(std::pair<K, V>(newKey, val));
79
-    }
80
-
81
-    iterator begin()
82
-    {
83
-        return iterator(mMap.begin());
84
-    }
85
-
86
-    iterator end()
87
-    {
88
-        return iterator(mMap.end());
89
-    }
90
-    
91
-private:
92
-
93
-    std::map<K, V> mMap;
94
-};
95
-
96
-#elif ACTIVE_VERSION == 1
97
-
98
-#include <map>
99
-
100
-template <class K, class V>
101
-class MutableMap
102
-{
103
-public:
104
-
105
-    MutableMap() : mMap() {}
106
-
107
-    class iterator
108
-    {
109
-    public:
110
-    
111
-        iterator() : mIt() {}
112
-        iterator(const iterator& it) : mIt(it.mIt) {}
113
-        ~iterator() {}
114
-        iterator& operator=(const iterator &it) { mIt = it.mIt; return *this; }
115
-        iterator& operator++() { ++mIt; return *this; }
116
-        iterator& operator--() { --mIt; return *this; }
117
-        std::pair<const K, V>& operator*() { return *mIt; }
118
-        std::pair<const K, V>* operator->() { return &(operator *()); }
119
-        bool operator==(const iterator &it) { return mIt == it.mIt; }
120
-        bool operator!=(const iterator &it) { return !(*this == it); }
121
-
122
-    private:
123
-
124
-        friend class MutableMap;
125
-        iterator(typename std::map<K, V>::iterator it) : mIt(it) {}
34
+        explicit iterator(typename std::map<K, V>::iterator it) : mIt(it) {}
126 35
         typename std::map<K, V>::iterator mIt;
127 36
     };
128 37
 
129 38
     class const_iterator
130 39
     {
131 40
     public:
132
-    
133 41
         const_iterator() : mIt() {}
134 42
         const_iterator(const const_iterator& it) : mIt(it.mIt) {}
135 43
         ~const_iterator() {}
... ...
@@ -140,11 +48,9 @@ public:
140 48
         const std::pair<const K, V>* operator->() { return &(operator *()); }
141 49
         bool operator==(const const_iterator &it) { return mIt == it.mIt; }
142 50
         bool operator!=(const const_iterator &it) { return !(*this == it); }
143
-
144 51
     private:
145
-
146 52
         friend class MutableMap;
147
-        const_iterator(typename std::map<K, V>::const_iterator it) : mIt(it) {}
53
+        explicit const_iterator(typename std::map<K, V>::const_iterator it) : mIt(it) {}
148 54
         typename std::map<K, V>::const_iterator mIt;
149 55
     };
150 56
 
... ...
@@ -172,7 +78,7 @@ public:
172 78
 
173 79
     void updateKey(iterator it, const K &newKey)
174 80
     {
175
-        const_cast<uint64_t&>((*it).first) = newKey;
81
+        const_cast<uint64_t&>((*it).first) = newKey; // TODO cleaner solution that this
176 82
     }
177 83
 
178 84
     iterator begin()
... ...
@@ -194,16 +100,12 @@ public:
194 100
     {
195 101
         return const_iterator(mMap.end());
196 102
     }
197
-    
198
-private:
199 103
 
104
+private:
200 105
     std::map<K, V> mMap;
201 106
 };
202 107
 
203
-#elif ACTIVE_VERSION == 2
204
-
205
-#include <map>
206
-
108
+#if 0
207 109
 template <class K>
208 110
 class KeyWrapper
209 111
 {
... ...
@@ -225,80 +127,6 @@ private:
225 127
 
226 128
     mutable K mKey;
227 129
 };
228
-
229
-template <class K, class V>
230
-class MutableMap
231
-{
232
-public:
233
-
234
-    MutableMap() : mMap() {}
235
-
236
-    class iterator
237
-    {
238
-    public:
239
-    
240
-        iterator() : mIt() {}
241
-        iterator(const iterator& it) : mIt(it.mIt) {}
242
-        ~iterator() {}
243
-        iterator& operator=(const iterator &it) { mIt = it.mIt; return *this; }
244
-        iterator& operator++() { ++mIt; return *this; }
245
-        iterator& operator--() { --mIt; return *this; }
246
-        std::pair<K, V> operator*() { return *mIt; }
247
-        //std::pair<K, V>* operator->() { return &(operator *()); }
248
-        bool operator==(const iterator &it) { return mIt == it.mIt; }
249
-        bool operator!=(const iterator &it) { return !(*this == it); }
250
-
251
-    private:
252
-
253
-        friend class MutableMap;
254
-        iterator(typename std::map< KeyWrapper<K>, V>::iterator it) : mIt(it) {}
255
-        typename std::map< KeyWrapper<K>, V>::iterator mIt;
256
-    };
257
-
258
-    unsigned count(const K &key) const
259
-    {
260
-        return mMap.count(KeyWrapper(key));
261
-    }
262
-    
263
-    std::pair<iterator, bool> insert(const std::pair<K, V> &val)
264
-    {
265
-        std::pair< typename std::map< KeyWrapper<K>, V>::iterator, bool> result = mMap.insert(val);
266
-        iterator it(result.first);
267
-        return std::pair<iterator, bool>(it, result.second);
268
-    }
269
-
270
-    void erase(iterator it)
271
-    {
272
-        mMap.erase(it.mIt);
273
-    }
274
-
275
-    iterator find(const K &key)
276
-    {
277
-        return mMap.find(KeyWrapper(key));
278
-    }
279
-
280
-    void updateKey(iterator it, const K &newKey)
281
-    {
282
-        V val = (*it).second;
283
-        mMap.erase(it.mIt);
284
-        insert(std::pair< KeyWrapper<K>, V>(KeyWrapper(newKey), val));
285
-    }
286
-
287
-    iterator begin()
288
-    {
289
-        return iterator(mMap.begin());
290
-    }
291
-
292
-    iterator end()
293
-    {
294
-        return iterator(mMap.end());
295
-    }
296
-    
297
-private:
298
-
299
-    std::map< KeyWrapper<K>, V> mMap;
300
-};
301
-
302 130
 #endif
303 131
 
304 132
 #endif // __COGAPS_MUTABLE_MAP_H__
... ...
@@ -1,4 +1,5 @@
1 1
 #include "SparseIterator.h"
2
+#include "../utils/GapsAssert.h"
2 3
 
3 4
 // counts number of bits set below position
4 5
 // internal expression clears all bits as high as pos or higher
... ...
@@ -111,7 +112,7 @@ bool SparseIterator<1>::atEnd() const
111 112
 void SparseIterator<1>::next()
112 113
 {
113 114
     ++mSparseIndex;
114
-    while (!mSparseFlags)
115
+    while (mSparseFlags == 0u)
115 116
     {
116 117
         // advance to next chunk
117 118
         if (++mBigIndex == mTotalIndices)
... ...
@@ -1,8 +1,11 @@
1 1
 #include "SparseMatrix.h"
2
+#include "Matrix.h"
2 3
 #include "../file_parser/FileParser.h"
4
+#include "../utils/Archive.h"
5
+#include "../utils/GapsAssert.h"
3 6
 
4
-#include <iterator>
5 7
 #include <algorithm>
8
+#include <iterator>
6 9
 
7 10
 // constructor from data set read in as a matrix
8 11
 SparseMatrix::SparseMatrix(const Matrix &mat, bool genesInCols,
... ...
@@ -1,36 +1,30 @@
1 1
 #ifndef __COGAPS_SPARSE_MATRIX_H__
2 2
 #define __COGAPS_SPARSE_MATRIX_H__
3 3
 
4
-#include "Matrix.h"
5 4
 #include "SparseVector.h"
6
-#include "../utils/Archive.h"
7 5
 
6
+#include <string>
8 7
 #include <vector>
9 8
 
9
+class Archive;
10
+class Matrix;
11
+
10 12
 // no random access, all data is const, can only access with iterator
11 13
 // over a given column
12 14
 class SparseMatrix
13 15
 {
14 16
 public:
15
-
16 17
     SparseMatrix(const Matrix &mat, bool genesInCols, bool subsetGenes,
17 18
         std::vector<unsigned> indices);
18
-
19 19
     SparseMatrix(const std::string &path, bool genesInCols, bool subsetGenes,
20 20
         std::vector<unsigned> indices);
21
-
22 21
     unsigned nRow() const;
23 22
     unsigned nCol() const;
24
-
25 23
     const SparseVector& getCol(unsigned n) const;
26
-
27 24
     void operator=(const Matrix &mat);
28
-
29 25
     friend Archive& operator<<(Archive &ar, const SparseMatrix &vec);
30 26
     friend Archive& operator>>(Archive &ar, SparseMatrix &vec);
31
-
32 27
 private:
33
-
34 28
     std::vector<SparseVector> mCols;
35 29
     unsigned mNumRows;
36 30
     unsigned mNumCols;
... ...
@@ -1,5 +1,7 @@
1 1
 #include "SparseVector.h"
2 2
 #include "Vector.h"
3
+#include "../utils/Archive.h"
4
+#include "../utils/GapsAssert.h"
3 5
 
4 6
 // counts number of bits set below position
5 7
 // internal expression clears all bits as high as pos or higher
... ...
@@ -53,17 +55,13 @@ unsigned SparseVector::size() const
53 55
 void SparseVector::insert(unsigned i, float v)
54 56
 {
55 57
     GAPS_ASSERT(v > 0.f);
56
-
57
-    // this data should not exist
58
-    GAPS_ASSERT(!(mIndexBitFlags[i / 64] & (1ull << (i % 64))));
59
-
58
+    GAPS_ASSERT(!(mIndexBitFlags[i / 64] & (1ull << (i % 64)))); // this data should not exist
60 59
     unsigned dataIndex = 0;
61 60
     for (unsigned j = 0; j < i / 64; ++j)
62 61
     {
63 62
         dataIndex += __builtin_popcountll(mIndexBitFlags[j]);
64 63
     }
65 64
     dataIndex += countLowerBits(mIndexBitFlags[i / 64], i % 64);
66
-
67 65
     mData.insert(mData.begin() + dataIndex, v);
68 66
     mIndexBitFlags[i / 64] |= (1ull << (i % 64));
69 67
 }
... ...
@@ -88,11 +86,10 @@ Vector SparseVector::getDense() const
88 86
 
89 87
 float SparseVector::at(unsigned n) const
90 88
 {
91
-    if (!(mIndexBitFlags[n / 64] & (1ull << (n % 64))))
89
+    if ((mIndexBitFlags[n / 64] & (1ull << (n % 64))) == 0u)
92 90
     {
93 91
         return 0.f;
94 92
     }
95
-
96 93
     unsigned sparseNdx = 0;
97 94
     for (unsigned i = 0; i < n / 64; ++i)
98 95
     {
... ...
@@ -1,6 +1,7 @@
1 1
 #ifndef __COGAPS_SPARSE_VECTOR_H__
2 2
 #define __COGAPS_SPARSE_VECTOR_H__
3 3
 
4
+#include <stdint.h>
4 5
 #include <vector>
5 6
 
6 7
 class Archive;
... ...
@@ -22,7 +23,6 @@ public:
22 23
     unsigned nElements() const;
23 24
     friend Archive& operator<<(Archive &ar, const SparseVector &vec);
24 25
     friend Archive& operator>>(Archive &ar, SparseVector &vec);
25
-
26 26
 private:
27 27
     template <unsigned N>
28 28
     friend class SparseIterator;
... ...
@@ -1,5 +1,7 @@
1 1
 #include "Vector.h"
2 2
 #include "../math/SIMD.h"
3
+#include "../utils/Archive.h"
4
+#include "../utils/GapsAssert.h"
3 5
 
4 6
 #define SIMD_PAD(x) (gaps::simd::Index::increment() + \
5 7
     gaps::simd::Index::increment() * ((x) / gaps::simd::Index::increment())) 
... ...
@@ -9,7 +11,7 @@ Vector::Vector(unsigned sz)
9 11
 mData(SIMD_PAD(sz), 0.f),
10 12
 mSize(sz)
11 13
 {
12
-    GAPS_ASSERT(mData.size() % gaps::simd::Index::increment() == 0);
14
+    GAPS_ASSERT((mData.size() % gaps::simd::Index::increment()) == 0);
13 15
 }
14 16
 
15 17
 Vector::Vector(const std::vector<float> &v)
... ...
@@ -17,7 +19,7 @@ Vector::Vector(const std::vector<float> &v)
17 19
 mData(SIMD_PAD(v.size()), 0.f),
18 20
 mSize(v.size())
19 21
 {
20
-    GAPS_ASSERT(mData.size() % gaps::simd::Index::increment() == 0);
22
+    GAPS_ASSERT((mData.size() % gaps::simd::Index::increment()) == 0);
21 23
     for (unsigned i = 0; i < mSize; ++i)
22 24
     {
23 25
         mData[i] = v[i];
... ...
@@ -1,10 +1,13 @@
1 1
 #ifndef __COGAPS_VECTOR_H__
2 2
 #define __COGAPS_VECTOR_H__
3 3
 
4
-#include "../utils/Archive.h"
5
-
6
-#include <boost/align/aligned_allocator.hpp>
7 4
 #include <vector>
5
+#pragma GCC diagnostic push
6
+#pragma GCC diagnostic ignored "-Wdeprecated-declarations"
7
+#include <boost/align/aligned_allocator.hpp>
8
+#pragma GCC diagnostic pop
9
+
10
+class Archive;
8 11
 
9 12
 // need to align data for SIMD
10 13
 namespace bal = boost::alignment;
... ...
@@ -14,29 +17,20 @@ typedef std::vector<float, bal::aligned_allocator<float,32> > aligned_vector;
14 17
 class Vector
15 18
 {
16 19
 public:
17
-
18 20
     explicit Vector(unsigned sz);
19 21
     explicit Vector(const std::vector<float> &v);
20
-
21 22
     float operator[](unsigned i) const;
22 23
     float& operator[](unsigned i);
23
-
24 24
     const float* ptr() const;
25 25
     float* ptr();
26
-
27 26
     unsigned size() const;
28 27
     void pad(float val);
29
-
30 28
     void operator+=(const Vector &v);
31
-
32 29
     void operator*=(float f);
33 30
     void operator/=(float f);
34
-    
35 31
     friend Archive& operator<<(Archive &ar, const Vector &vec);
36 32
     friend Archive& operator>>(Archive &ar, Vector &vec);
37
-
38 33
 private:
39
-
40 34
     aligned_vector mData;
41 35
     unsigned mSize;
42 36
 };
... ...
@@ -28,7 +28,9 @@ static std::string trimNewline(const std::string &s)
28 28
 {
29 29
     std::size_t pos;
30 30
     if ((pos = s.find('\n')) == std::string::npos)
31
+    {
31 32
         return s;
33
+    }
32 34
     return s.substr(0, pos);
33 35
 }
34 36
 
... ...
@@ -37,8 +39,10 @@ static std::vector<std::string> split(const std::string &s, char delimiter)
37 39
     std::vector<std::string> tokens;
38 40
     std::string temp;
39 41
     std::stringstream ss(s);
40
-    while (std::getline(ss, temp, delimiter))
42
+    while (std::getline(ss, temp, delimiter) != 0)
43
+    {
41 44
         tokens.push_back(trim(temp));
45
+    }
42 46
     return tokens;    
43 47
 }
44 48
 
... ...
@@ -94,8 +98,10 @@ mDelimiter(delimiter), mGctFormat(gctFormat)
94 98
         while ((pos = line.find('\n')) == std::string::npos);
95 99
 
96 100
         // find number of rows
97
-        while (std::getline(mFile, line))
101
+        while (std::getline(mFile, line) != 0)
102
+        {
98 103
             ++mNumRows;
104
+        }
99 105
     }
100 106
 
101 107
     // reset file stream to beginning
... ...
@@ -131,9 +137,13 @@ void CharacterDelimitedParser::parseNextLine()
131 137
     std::getline(mFile, fullLine);
132 138
     mCurrentLine = split(fullLine, mDelimiter);
133 139
     if (mRowNamesPresent)
140
+    {
134 141
         mCurrentLine.erase(mCurrentLine.begin());
142
+    }
135 143
     if (mGctFormat)
144
+    {
136 145
         mCurrentLine.erase(mCurrentLine.begin(), mCurrentLine.begin() + 2);
146
+    }
137 147
 }
138 148
 
139 149
 MatrixElement CharacterDelimitedParser::getNext()
... ...
@@ -2,50 +2,41 @@
2 2
 #define __COGAPS_CHARACTER_DELIMITED_PARSER_H__
3 3
 
4 4
 #include "FileParser.h"
5
-#include "MatrixElement.h"
6 5
 
7 6
 #include <fstream>
8
-#include <vector>
9 7
 #include <string>
8
+#include <vector>
9
+
10
+struct MatrixElement;
10 11
 
11 12
 class CharacterDelimitedParser : public AbstractFileParser
12 13
 {
13 14
 public:
14
-
15 15
     explicit CharacterDelimitedParser(const std::string &path, char delimiter, bool gctFormat=false);
16 16
     ~CharacterDelimitedParser();
17
-
18 17
     std::vector<std::string> rowNames() const;
19 18
     std::vector<std::string> colNames() const;
20 19
     unsigned nRow() const;
21 20
     unsigned nCol() const;
22 21
     bool hasNext();
23 22
     MatrixElement getNext();
24
-
25 23
 private:
24
+    CharacterDelimitedParser(const CharacterDelimitedParser &p); // don't allow copies
25
+    CharacterDelimitedParser& operator=(const CharacterDelimitedParser &p); // don't allow copies
26
+    void parseNextLine();
27
+    void checkFileState() const;
26 28
 
27 29
     std::ifstream mFile;
28
-
29 30
     std::vector<std::string> mRowNames;
30 31
     std::vector<std::string> mColNames;
31
-
32 32
     std::vector<std::string> mCurrentLine;
33
-
34 33
     unsigned mNumRows;
35 34
     unsigned mNumCols;
36
-
37 35
     unsigned mCurrentRow;
38 36
     unsigned mCurrentCol;
39
-
40 37
     bool mRowNamesPresent;
41 38
     char mDelimiter;
42 39
     bool mGctFormat;
43
-
44
-    void parseNextLine();
45
-    void checkFileState() const;
46
-
47
-    CharacterDelimitedParser(const CharacterDelimitedParser &p); // don't allow copies
48
-    CharacterDelimitedParser& operator=(const CharacterDelimitedParser &p); // don't allow copies
49 40
 };
50 41
 
51 42
 #endif // __COGAPS_CHARACTER_DELIMITED_PARSER_H__
52 43
\ No newline at end of file
... ...
@@ -1,8 +1,8 @@
1 1
 #ifndef __COGAPS_FILE_PARSER_H__
2 2
 #define __COGAPS_FILE_PARSER_H__
3 3
 
4
-#include "../utils/GapsAssert.h"
5 4
 #include "MatrixElement.h"
5
+#include "../utils/GapsAssert.h"
6 6
 
7 7
 #include <fstream>
8 8
 #include <vector>
... ...
@@ -20,21 +20,16 @@ enum GapsFileType
20 20
 class AbstractFileParser
21 21
 {
22 22
 public:
23
-
24 23
     static AbstractFileParser* create(const std::string &path);
25
-
26 24
     AbstractFileParser();
27 25
     virtual ~AbstractFileParser() = 0;
28
-
29 26
     virtual std::vector<std::string> rowNames() const;
30 27
     virtual std::vector<std::string> colNames() const;
31 28
     virtual unsigned nRow() const = 0;
32 29
     virtual unsigned nCol() const = 0;
33 30
     virtual bool hasNext() = 0;
34 31
     virtual MatrixElement getNext() = 0;
35
-
36 32
 private:
37
-
38 33
     AbstractFileParser(const AbstractFileParser &p); // don't allow copies
39 34
     AbstractFileParser& operator=(const AbstractFileParser &p); // don't allow copies
40 35
 };
... ...
@@ -43,38 +38,31 @@ private:
43 38
 class FileParser
44 39
 {
45 40
 public:
46
-
47 41
     explicit FileParser(const std::string &path);
48 42
     ~FileParser();
49
-
50 43
     std::vector<std::string> rowNames() const;
51 44
     std::vector<std::string> colNames() const;
52 45
     unsigned nRow() const;
53 46
     unsigned nCol() const;
54 47
     bool hasNext();
55 48
     MatrixElement getNext();
56
-
57 49
     static GapsFileType fileType(const std::string &path);
58
-
59 50
     template <class MatrixType>
60 51
     static void writeToCsv(const std::string &path, const MatrixType &mat);
61
-
62 52
 private:
63
-
64
-    AbstractFileParser *mParser;
65
-
66 53
     FileParser(const FileParser &p); // don't allow copies
67 54
     FileParser& operator=(const FileParser &p); // don't allow copies
68
-};
69 55
 
70
-// temporary solution - should be moved into specific file parsers, ok for now
71
-// since writing to file not exposed to user, only used for internal testing
56
+    AbstractFileParser *mParser;
57
+};
72 58
 
73 59
 template <class MatrixType>
74 60
 void FileParser::writeToCsv(const std::string &path, const MatrixType &mat)
75 61
 {
76
-    GAPS_ASSERT(FileParser::fileType(path) == GAPS_CSV);
77
-
62
+    if (FileParser::fileType(path) != GAPS_CSV)
63
+    {
64
+        GAPS_ERROR("output file must be a csv");
65
+    }
78 66
     std::ofstream outputFile;
79 67
     outputFile.open(path.c_str());
80 68
     outputFile << "\"\"";
... ...
@@ -85,7 +73,6 @@ void FileParser::writeToCsv(const std::string &path, const MatrixType &mat)
85 73
         outputFile << ",\"Col" << i << "\"";
86 74
     }
87 75
     outputFile << "\n";
88
-
89 76
     for (unsigned i = 0; i < mat.nRow(); ++i)
90 77
     {
91 78
         // write row names
... ...
@@ -5,12 +5,12 @@
5 5
 
6 6
 struct MatrixElement
7 7
 {
8
+    MatrixElement(unsigned r, unsigned c, float v);
9
+    MatrixElement(unsigned r, unsigned c, const std::string &s);
10
+
8 11
     unsigned row;
9 12
     unsigned col;
10 13
     float value;
11
-
12
-    MatrixElement(unsigned r, unsigned c, float v);
13
-    MatrixElement(unsigned r, unsigned c, const std::string &s);
14 14
 };
15 15
 
16 16
 #endif // __COGAPS_MATRIX_ELEMENT_H__
... ...
@@ -2,35 +2,29 @@
2 2
 #define __COGAPS_MTX_PARSER_H__
3 3
 
4 4
 #include "FileParser.h"
5
-#include "MatrixElement.h"
6 5
 
7 6
 #include <fstream>
8 7
 #include <string>
9 8
 
9
+struct MatrixElement;
10
+
10 11
 class MtxParser : public AbstractFileParser
11 12
 {
12 13
 public:
13
-
14 14
     explicit MtxParser(const std::string &path);
15 15
     ~MtxParser();
16
-
17 16
     unsigned nRow() const;
18 17
     unsigned nCol() const;
19
-
20 18
     bool hasNext();
21 19
     MatrixElement getNext();
22
-
23 20
 private:
21
+    MtxParser(const MtxParser &p); // don't allow copies
22
+    MtxParser& operator=(const MtxParser &p); // don't allow copies
23
+    void checkFileState() const;
24 24
 
25 25
     std::ifstream mFile;
26
-
27 26
     unsigned mNumRows;
28 27
     unsigned mNumCols;
29
-
30
-    MtxParser(const MtxParser &p); // don't allow copies
31
-    MtxParser& operator=(const MtxParser &p); // don't allow copies
32
-
33
-    void checkFileState() const;
34 28
 };
35 29
 
36
-#endif
37 30
\ No newline at end of file
31
+#endif // __COGAPS_MTX_PARSER_H__
38 32
\ No newline at end of file
... ...
@@ -1,5 +1,6 @@
1 1
 #include "AlphaParameters.h"
2 2
 #include "../math/Math.h"
3
+#include "../math/Random.h"
3 4
 
4 5
 #include <cmath>
5 6
 
... ...
@@ -9,7 +10,7 @@ AlphaParameters::AlphaParameters(float t_s, float t_smu)
9 10
 
10 11
 AlphaParameters AlphaParameters::operator+(const AlphaParameters &other) const
11 12
 {
12
-    return AlphaParameters(s + other.s, s_mu - other.s_mu); // not a typo
13
+    return AlphaParameters(s + other.s, s_mu - other.s_mu); // minus sign not a typo
13 14
 }
14 15
 
15 16
 AlphaParameters AlphaParameters::operator*(float v) const
... ...
@@ -19,7 +20,8 @@ AlphaParameters AlphaParameters::operator*(float v) const
19 20
 
20 21
 void AlphaParameters::operator*=(float v)
21 22
 {
22
-    s *= v; s_mu *= v;
23
+    s *= v;
24
+    s_mu *= v;
23 25
 }
24 26
 
25 27
 OptionalFloat gibbsMass(AlphaParameters alpha, float a, float b, GapsRng *rng)
... ...
@@ -1,18 +1,18 @@
1 1
 #ifndef __COGAPS_ALPHA_PARAMETERS_H__
2 2
 #define __COGAPS_ALPHA_PARAMETERS_H__
3 3
 
4
-#include "../math/Random.h"
4
+class OptionalFloat;
5
+class GapsRng;
5 6
 
6 7
 struct AlphaParameters
7 8
 {
8
-    float s;
9
-    float s_mu;
10
-    
11 9
     AlphaParameters(float t_s, float t_smu);
12
-
13 10
     AlphaParameters operator+(const AlphaParameters &other) const;
14 11
     AlphaParameters operator*(float v) const;
15 12
     void operator*=(float v);
13
+
14
+    float s;
15
+    float s_mu;
16 16
 };
17 17
 
18 18
 OptionalFloat gibbsMass(AlphaParameters alpha, float a, float b, GapsRng *rng);
... ...
@@ -7,13 +7,13 @@
7 7
 #include "../math/Math.h"
8 8
 #include "../math/VectorMath.h"
9 9
 #include "../math/MatrixMath.h"
10
-#include "../math/Random.h"
11 10
 #include "../GapsParameters.h"
11
+#include "../math/Random.h"
12 12
 
13
-#include <vector>
14
-#include <limits>
15 13
 #include <cmath>
14
+#include <limits>
16 15
 #include <stdint.h>
16
+#include <vector>
17 17
 
18 18
 //////////////////////////// AsynchronousGibbsSampler Interface ////////////////////////////
19 19
 
... ...
@@ -23,10 +23,10 @@ template <class DataModel>
23 23
 class AsynchronousGibbsSampler;
24 24
 
25 25
 template <class DataModel>
26
-Archive& operator<<(Archive &ar, const AsynchronousGibbsSampler<DataModel> &sampler);
26
+Archive& operator<<(Archive &ar, const AsynchronousGibbsSampler<DataModel> &s);
27 27
 
28 28
 template <class DataModel>
29
-Archive& operator>>(Archive &ar, AsynchronousGibbsSampler<DataModel> &sampler);
29
+Archive& operator>>(Archive &ar, AsynchronousGibbsSampler<DataModel> &s);
30 30
 
31 31
 template <class DataModel>
32 32
 class AsynchronousGibbsSampler : public DataModel
... ...
@@ -37,31 +37,24 @@ public:
37 37
     AsynchronousGibbsSampler(const DataType &data, bool transpose, bool subsetRows,
38 38
         float alpha, float maxGibbsMass, const GapsParameters &params,
39 39
         GapsRandomState *randState);
40
-
41 40
     unsigned nAtoms() const;
42 41
     float getAverageQueueLength() const;
43
-
44 42
     void update(unsigned nSteps, unsigned nThreads);
45
-
46 43
     friend Archive& operator<< <DataModel> (Archive &ar, const AsynchronousGibbsSampler &s);
47 44
     friend Archive& operator>> <DataModel> (Archive &ar, AsynchronousGibbsSampler &s);
48
-
49 45
 private:
50
-
51
-    ConcurrentAtomicDomain mDomain; // data structure providing access to atoms
52
-    ProposalQueue mQueue; // creates queue of proposals that get evaluated by sampler
53
-
54
-    float mAvgQueueLength;
55
-    float mNumQueueSamples;
56
-
57 46
     void birth(const AtomicProposal &prop);
58 47
     void death(const AtomicProposal &prop);
59 48
     void move(const AtomicProposal &prop);
60 49
     void exchange(const AtomicProposal &prop);
61
-
62 50
 #ifdef GAPS_DEBUG
63 51
     float maximumDrift() const;
64 52
 #endif
53
+
54
+    ConcurrentAtomicDomain mDomain; // data structure providing access to atoms
55
+    ProposalQueue mQueue; // creates queue of proposals that get evaluated by sampler
56
+    float mAvgQueueLength;
57
+    float mNumQueueSamples;
65 58
 };
66 59
 
67 60
 //////////////////// AsynchronousGibbsSampler - templated functions ////////////////////////
... ...
@@ -72,7 +65,7 @@ AsynchronousGibbsSampler<DataModel>::AsynchronousGibbsSampler(const DataType &da
72 65
 bool transpose, bool subsetRows, float alpha, float maxGibbsMass,
73 66
 const GapsParameters &params, GapsRandomState *randState)
74 67
     :
75
-DataModel(data, transpose, subsetRows, params, alpha),
68
+DataModel(data, transpose, subsetRows, params, alpha, maxGibbsMass),
76 69
 mDomain(DataModel::nElements()),
77 70
 mQueue(DataModel::nElements(), DataModel::nPatterns(), randState),
78 71
 mAvgQueueLength(0),
... ...
@@ -201,7 +194,6 @@ template <class DataModel>
201 194
 void AsynchronousGibbsSampler<DataModel>::exchange(const AtomicProposal &prop)
202 195
 {
203 196
     GAPS_ASSERT(prop.r1 != prop.r2 || prop.c1 != prop.c2);
204
-
205 197
     if (DataModel::canUseGibbs(prop.c1, prop.c2))
206 198
     {
207 199
         OptionalFloat mass = DataModel::sampleExchange(prop.r1, prop.c1, prop.atom1->mass(),
... ...
@@ -1,5 +1,8 @@
1 1
 #include "DenseNormalModel.h"
2 2
 #include "../math/Math.h"
3
+#include "../math/Random.h"
4
+#include "../utils/Archive.h"
5
+#include "../utils/GapsAssert.h"
3 6
 
4 7
 #define GAPS_SQ(x) ((x) * (x))
5 8
 
... ...
@@ -18,10 +21,8 @@ void DenseNormalModel::sync(const DenseNormalModel &model, unsigned nThreads)
18 21
 {
19 22
     GAPS_ASSERT(model.mAPMatrix.nRow() == mAPMatrix.nCol());
20 23
     GAPS_ASSERT(model.mAPMatrix.nCol() == mAPMatrix.nRow());
21
-
22 24
     unsigned nc = model.mAPMatrix.nCol();
23 25
     unsigned nr = model.mAPMatrix.nRow();
24
-
25 26
     #pragma omp parallel for num_threads(nThreads)
26 27
     for (unsigned j = 0; j < nc; ++j)