Browse code

tests passing

Tom Sherman authored on 25/10/2018 22:31:44
Showing38 changed files

... ...
@@ -4,7 +4,7 @@ Archive& operator<<(Archive &ar, const GapsParameters &p)
4 4
 {
5 5
     ar << p.seed << p.nGenes << p.nSamples << p.nPatterns << p.nIterations
6 6
         << p.alphaA << p.alphaP << p.maxGibbsMassA << p.maxGibbsMassP
7
-        << p.singleCell << p.useSparseOptimization;
7
+        << p.singleCell << p.useSparseOptimization << p.checkpointInterval;
8 8
     return ar;
9 9
 }
10 10
 
... ...
@@ -12,7 +12,7 @@ Archive& operator>>(Archive &ar, GapsParameters &p)
12 12
 {
13 13
     ar >> p.seed >> p.nGenes >> p.nSamples >> p.nPatterns >> p.nIterations
14 14
         >> p.alphaA >> p.alphaP >> p.maxGibbsMassA >> p.maxGibbsMassP
15
-        >> p.singleCell >> p.useSparseOptimization;
15
+        >> p.singleCell >> p.useSparseOptimization >> p.checkpointInterval;
16 16
     return ar;
17 17
 }
18 18
     
... ...
@@ -15,6 +15,14 @@
15 15
 namespace bpt = boost::posix_time;
16 16
 #define bpt_now() bpt::microsec_clock::local_time()
17 17
 
18
+// for conditionally printing status messages
19
+#define GAPS_MESSAGE(b, m) \
20
+    do { \
21
+        if (b) { \
22
+            gaps_printf(m); \
23
+        } \
24
+    } while(0)
25
+
18 26
 // forward declaration
19 27
 template <class Sampler, class DataType>
20 28
 static GapsResult runCoGAPSAlgorithm(const DataType &data, GapsParameters &params,
... ...
@@ -160,11 +168,12 @@ Sampler &PSampler, unsigned nA, unsigned nP)
160 168
 
161 169
 template <class Sampler>
162 170
 static void createCheckpoint(const GapsParameters &params,
163
-const Sampler &ASampler, const Sampler &PSampler, const GapsRandomState *randState,
164
-const GapsRng &rng, char phase, unsigned iter)
171
+Sampler &ASampler, Sampler &PSampler, const GapsRandomState *randState,
172
+const GapsStatistics &stats, const GapsRng &rng, char phase, unsigned iter)
165 173
 {
166 174
     if (params.checkpointInterval > 0
167
-    && ((iter + 1) % params.checkpointInterval) == 0)
175
+    && ((iter + 1) % params.checkpointInterval) == 0
176
+    && !params.subsetData)
168 177
     {
169 178
         // create backup file
170 179
         std::rename(params.checkpointOutFile.c_str(),
... ...
@@ -174,10 +183,28 @@ const GapsRng &rng, char phase, unsigned iter)
174 183
         Archive ar(params.checkpointOutFile, ARCHIVE_WRITE);
175 184
         ar << params;
176 185
         ar << *randState;
177
-        ar << ASampler << PSampler << phase << iter << rng;
186
+        ar << ASampler << PSampler << stats << phase << iter << rng;
178 187
         
179 188
         // delete backup file
180 189
         std::remove((params.checkpointOutFile + ".backup").c_str());
190
+
191
+        ASampler.extraInitialization();
192
+        PSampler.extraInitialization();
193
+    }
194
+}
195
+
196
+template <class Sampler>
197
+static void processCheckpoint(GapsParameters &params, Sampler &ASampler,
198
+Sampler &PSampler, GapsRandomState *randState, GapsStatistics &stats,
199
+GapsRng &rng, char &phase, unsigned &currentIter)
200
+{    
201
+    // check if running from checkpoint, get all saved data
202
+    if (params.useCheckPoint)
203
+    {
204
+        Archive ar(params.checkpointFile, ARCHIVE_READ);
205
+        ar >> params;
206
+        ar >> *randState;
207
+        ar >> ASampler >> PSampler >> stats >> phase >> currentIter >> rng;
181 208
     }
182 209
 }
183 210
 
... ...
@@ -192,8 +219,8 @@ GapsRng &rng, bpt::ptime startTime, char phase, unsigned &currentIter)
192 219
         Rcpp::checkUserInterrupt();
193 220
         #endif
194 221
 
195
-        createCheckpoint(params, ASampler, PSampler, randState, rng, phase,
196
-            currentIter);
222
+        createCheckpoint(params, ASampler, PSampler, randState, stats,
223
+            rng, phase, currentIter);
197 224
         
198 225
         // set annealing temperature in calibration phase
199 226
         if (phase == 'C')
... ...
@@ -217,51 +244,10 @@ GapsRng &rng, bpt::ptime startTime, char phase, unsigned &currentIter)
217 244
     }
218 245
 }
219 246
 
220
-// here is the CoGAPS algorithm
221
-template <class Sampler, class DataType>
222
-static GapsResult runCoGAPSAlgorithm(const DataType &data, GapsParameters &params,
223
-const DataType &uncertainty, GapsRandomState *randState)
247
+template <class Sampler>
248
+static void processFixedMatrix(const GapsParameters &params, Sampler &ASampler,
249
+Sampler &PSampler)
224 250
 {
225
-    // check if running in debug mode
226
-    #ifdef GAPS_DEBUG
227
-    if (params.printMessages)
228
-    {
229
-        gaps_printf("Running in debug mode\n");
230
-    }
231
-    #endif
232
-
233
-    // load data into gibbs samplers
234
-    // we transpose the data in the A sampler so that the update step
235
-    // is symmetrical for each sampler, this simplifies the code 
236
-    // within the sampler, note the subsetting genes/samples flag must be
237
-    // flipped if we are flipping the transpose flag
238
-    if (params.printMessages)
239
-    {
240
-        gaps_printf("Loading Data...");
241
-    }
242
-    Sampler ASampler(data, !params.transposeData, !params.subsetGenes,
243
-        params.alphaA, params.maxGibbsMassA, params, randState);
244
-    Sampler PSampler(data, params.transposeData, params.subsetGenes,
245
-        params.alphaP, params.maxGibbsMassP, params, randState);
246
-
247
-    // read in the uncertainty matrix if one is provided
248
-    if (!uncertainty.empty())
249
-    {
250
-        ASampler.setUncertainty(uncertainty, !params.transposeData,
251
-            !params.subsetGenes, params);
252
-        PSampler.setUncertainty(uncertainty, params.transposeData,
253
-            params.subsetGenes, params);
254
-    }
255
-    if (params.printMessages)
256
-    {
257
-        gaps_printf("Done!\n");
258
-    }
259
-
260
-    // these variables will get overwritten by checkpoint if provided
261
-    GapsRng rng(randState);
262
-    char phase = 'C';
263
-    unsigned currentIter = 0;
264
-
265 251
     // check if we're fixing a matrix
266 252
     if (params.useFixedMatrix)
267 253
     {
... ...
@@ -279,24 +265,10 @@ const DataType &uncertainty, GapsRandomState *randState)
279 265
             default: break; // 'N' for none
280 266
         }
281 267
     }
282
-     
283
-    // check if running from checkpoint, get all saved data
284
-    if (params.useCheckPoint)
285
-    {
286
-        Archive ar(params.checkpointFile, ARCHIVE_READ);
287
-        ar >> params;
288
-        ar >> *randState;
289
-        ar >> ASampler >> PSampler >> phase >> currentIter >> rng;
290
-    }
291
-
292
-    // sync samplers, second parameter indicates this should be a full sync
293
-    ASampler.sync(PSampler);
294
-    PSampler.sync(ASampler);
295
-
296
-    // sampler may need to run additional initialization after parameters set
297
-    ASampler.extraInitialization();
298
-    PSampler.extraInitialization();
268
+}
299 269
 
270
+static void calculateNumberOfThreads(GapsParameters params)
271
+{
300 272
     // calculate appropiate number of threads if compiled with openmp
301 273
     #ifdef __GAPS_OPENMP__
302 274
     if (params.printMessages && params.printThreadUsage)
... ...
@@ -307,36 +279,79 @@ const DataType &uncertainty, GapsRandomState *randState)
307 279
             params.maxThreads, availableThreads);
308 280
     }
309 281
     #endif
282
+}
283
+
284
+template <class Sampler, class DataType>
285
+static void processUncertainty(const GapsParameters params, Sampler &ASampler,
286
+Sampler &PSampler, const DataType &uncertainty)
287
+{
288
+    // read in the uncertainty matrix if one is provided
289
+    if (!uncertainty.empty())
290
+    {
291
+        ASampler.setUncertainty(uncertainty, !params.transposeData,
292
+            !params.subsetGenes, params);
293
+        PSampler.setUncertainty(uncertainty, params.transposeData,
294
+            params.subsetGenes, params);
295
+    }
296
+}
297
+
298
+// here is the CoGAPS algorithm
299
+template <class Sampler, class DataType>
300
+static GapsResult runCoGAPSAlgorithm(const DataType &data, GapsParameters &params,
301
+const DataType &uncertainty, GapsRandomState *randState)
302
+{
303
+    // check if running in debug mode
304
+    #ifdef GAPS_DEBUG
305
+    GAPS_MESSAGE(params.printMessages, "Running in debug mode\n");
306
+    #endif
307
+
308
+    // load data into gibbs samplers
309
+    // we transpose the data in the A sampler so that the update step
310
+    // is symmetrical for each sampler, this simplifies the code 
311
+    // within the sampler, note the subsetting genes/samples flag must be
312
+    // flipped if we are flipping the transpose flag
313
+    GAPS_MESSAGE(params.printMessages, "Loading Data...");
314
+    Sampler ASampler(data, !params.transposeData, !params.subsetGenes,
315
+        params.alphaA, params.maxGibbsMassA, params, randState);
316
+    Sampler PSampler(data, params.transposeData, params.subsetGenes,
317
+        params.alphaP, params.maxGibbsMassP, params, randState);
318
+    processUncertainty(params, ASampler, PSampler, uncertainty);
319
+    processFixedMatrix(params, ASampler, PSampler);
320
+    GAPS_MESSAGE(params.printMessages, "Done!\n");
321
+
322
+    // these variables will get overwritten by checkpoint if provided
323
+    GapsStatistics stats(params.nGenes, params.nSamples, params.nPatterns);
324
+    GapsRng rng(randState);
325
+    char phase = 'C';
326
+    unsigned currentIter = 0;
327
+    processCheckpoint(params, ASampler, PSampler, randState, stats, rng,
328
+        phase, currentIter);
329
+    calculateNumberOfThreads(params);
330
+
331
+    // sync samplers and run any additional initialization needed
332
+    ASampler.sync(PSampler);
333
+    PSampler.sync(ASampler);
334
+    ASampler.extraInitialization();
335
+    PSampler.extraInitialization();
310 336
 
311 337
     // record start time
312 338
     bpt::ptime startTime = bpt_now();
313 339
 
314
-    // cascade through phases, allows algorithm to be resumed in either phase
315
-    GapsStatistics stats(params.nGenes, params.nSamples, params.nPatterns);
340
+    // fallthrough through phases, allows algorithm to be resumed in any phase
316 341
     GAPS_ASSERT(phase == 'C' || phase == 'S');
317 342
     switch (phase)
318 343
     {
319 344
         case 'C':
320
-            if (params.printMessages)
321
-            {
322
-                gaps_printf("-- Calibration Phase --\n");
323
-            }
345
+            GAPS_MESSAGE(params.printMessages, "-- Calibration Phase --\n");
324 346
             runOnePhase(params, ASampler, PSampler, stats, randState, rng,
325 347
                 startTime, phase, currentIter);
326 348
             phase = 'S';
327 349
             currentIter = 0;
328 350
 
329
-
330
-
331 351
         case 'S':
332
-            if (params.printMessages)
333
-            {
334
-                gaps_printf("-- Sampling Phase --\n");
335
-            }
352
+            GAPS_MESSAGE(params.printMessages, "-- Sampling Phase --\n");
336 353
             runOnePhase(params, ASampler, PSampler, stats, randState, rng,
337 354
                 startTime, phase, currentIter);
338
-
339
-        default: break;
340 355
     }
341 356
     
342 357
     // get result
... ...
@@ -75,7 +75,23 @@ float GapsStatistics::meanChiSq(const SparseGibbsSampler &PSampler) const
75 75
 {
76 76
     Matrix A(mAMeanMatrix / static_cast<float>(mStatUpdates));
77 77
     Matrix P(mPMeanMatrix / static_cast<float>(mStatUpdates));
78
-    return 0.f; // TODO
78
+    
79
+    float chisq = 0.f;
80
+    for (unsigned i = 0; i < PSampler.mDMatrix.nRow(); ++i)
81
+    {
82
+        for (unsigned j = 0; j < PSampler.mDMatrix.nCol(); ++j)
83
+        {
84
+            float m = 0.f;
85
+            for (unsigned k = 0; k < A.nCol(); ++k)
86
+            {
87
+                m += A(i,k) * P(j,k);
88
+            }
89
+            float d = PSampler.mDMatrix.getCol(j).at(i);
90
+            float s = gaps::max(d * 0.1f, 0.1f);
91
+            chisq += GAPS_SQ((d - m) / s);
92
+        }
93
+    }
94
+    return chisq;
79 95
 }
80 96
 
81 97
 Archive& operator<<(Archive &ar, const GapsStatistics &stat)
... ...
@@ -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>
... ...
@@ -32,14 +32,17 @@ OBJECTS =   Cogaps.o \
32 32
             math/Random.o \
33 33
             math/VectorMath.o \
34 34
             cpp_tests/testAtomicDomain.o \
35
+            cpp_tests/testDenseGibbsSampler.o \
36
+            cpp_tests/testFileParsers.o \
35 37
             cpp_tests/testHashSets.o \
36 38
             cpp_tests/testHybridMatrix.o \
37 39
             cpp_tests/testHybridVector.o \
38 40
             cpp_tests/testMatrix.o \
41
+            cpp_tests/testSerialization.o \
42
+            cpp_tests/testSparseGibbsSampler.o \
39 43
             cpp_tests/testSparseIterator.o \
40 44
             cpp_tests/testSparseMatrix.o \
41 45
             cpp_tests/testSparseVector.o \
42
-            cpp_tests/testVector.o \
43
-            cpp_tests/testFileParsers.o \
44
-            cpp_tests/testDenseGibbsSampler.o \
45
-            cpp_tests/testSparseGibbsSampler.o
46 46
\ No newline at end of file
47
+            cpp_tests/testVector.o
48
+
49
+
... ...
@@ -30,7 +30,9 @@ public:
30 30
     friend Archive& operator<<(Archive& ar, const Atom &a);
31 31
     friend Archive& operator>>(Archive& ar, Atom &a);
32 32
 
33
+#ifndef GAPS_INTERNAL_TESTS
33 34
 private:
35
+#endif
34 36
 
35 37
     // used by the allocator
36 38
     friend class AtomAllocator;
... ...
@@ -51,7 +53,9 @@ public:
51 53
     void set(uint64_t n);
52 54
     void release(uint64_t n);
53 55
 
56
+#ifndef GAPS_INTERNAL_TESTS
54 57
 private:
58
+#endif
55 59
 
56 60
     uint64_t parts[NUM_INDEX_CHUNKS];
57 61
 };
... ...
@@ -67,7 +71,9 @@ public:
67 71
     void free(Atom* a);
68 72
     bool depleted() const;
69 73
 
74
+#ifndef GAPS_INTERNAL_TESTS
70 75
 private:
76
+#endif
71 77
     
72 78
     Atom* mPool;
73 79
     IndexFlagSet mIndexFlags;
... ...
@@ -86,7 +92,9 @@ public:
86 92
     Atom* createAtom(uint64_t pos, float mass);
87 93
     void destroyAtom(Atom *a);
88 94
 
95
+#ifndef GAPS_INTERNAL_TESTS
89 96
 private:
97
+#endif
90 98
 
91 99
     std::vector<AtomPool*> mPools;
92 100
     unsigned mIndex;
... ...
@@ -56,7 +56,7 @@ public:
56 56
 private:
57 57
 #endif
58 58
 
59
-    // only the proposal queue can insert
59
+    // only the proposal queue can insert, insert not thread-safe
60 60
     friend class ProposalQueue;
61 61
     Atom* insert(uint64_t pos, float mass);
62 62
 
... ...
@@ -6,7 +6,8 @@
6 6
 //////////////////////////////// AtomicProposal ////////////////////////////////
7 7
 
8 8
 AtomicProposal::AtomicProposal(char t, GapsRandomState *randState)
9
-    : rng(randState), pos(0), atom1(NULL), atom2(NULL), type(t)
9
+    : rng(randState), pos(0), atom1(NULL), atom2(NULL), r1(0), c1(0), r2(0),
10
+    c2(0), type(t)
10 11
 {}
11 12
     
12 13
 //////////////////////////////// ProposalQueue /////////////////////////////////
... ...
@@ -46,7 +47,7 @@ void ProposalQueue::populate(AtomicDomain &domain, unsigned limit)
46 47
     GAPS_ASSERT(mUsedMatrixIndices.isEmpty());
47 48
     GAPS_ASSERT(mProposedMoves.isEmpty());
48 49
     GAPS_ASSERT(mMinAtoms == mMaxAtoms);
49
-    GAPS_ASSERT(mMaxAtoms == domain.size());
50
+    GAPS_ASSERT_MSG(mMaxAtoms == domain.size(), mMaxAtoms << " != " << domain.size());
50 51
 
51 52
     unsigned nIter = 0;
52 53
     bool success = true;
... ...
@@ -115,6 +116,10 @@ float ProposalQueue::deathProb(double nAtoms) const
115 116
 
116 117
 bool ProposalQueue::makeProposal(AtomicDomain &domain)
117 118
 {
119
+    mU1 = mUseCachedRng ? mU1 : mRng.uniform();
120
+    mU2 = mUseCachedRng ? mU2: mRng.uniform();
121
+    mUseCachedRng = false;
122
+
118 123
     if (mMinAtoms < 2 && mMaxAtoms >= 2)
119 124
     {
120 125
         return false; // special indeterminate case
... ...
@@ -125,10 +130,6 @@ bool ProposalQueue::makeProposal(AtomicDomain &domain)
125 130
         return birth(domain); // always birth when no atoms exist
126 131
     }
127 132
 
128
-    mU1 = mUseCachedRng ? mU1 : mRng.uniform();
129
-    mU2 = mUseCachedRng ? mU2: mRng.uniform();
130
-    mUseCachedRng = false;
131
-
132 133
     float lowerBound = deathProb(static_cast<double>(mMinAtoms));
133 134
     float upperBound = deathProb(static_cast<double>(mMaxAtoms));
134 135
 
... ...
@@ -273,13 +274,15 @@ bool ProposalQueue::exchange(AtomicDomain &domain)
273 274
 Archive& operator<<(Archive &ar, const ProposalQueue &q)
274 275
 {
275 276
     ar << q.mRng << q.mMinAtoms << q.mMaxAtoms << q.mBinLength << q.mNumCols
276
-        << q.mAlpha << q.mDomainLength << q.mNumBins;
277
+        << q.mAlpha << q.mDomainLength << q.mNumBins << q.mLambda
278
+        << q.mUseCachedRng << q.mU1 << q.mU2;
277 279
     return ar;
278 280
 }
279 281
 
280 282
 Archive& operator>>(Archive &ar, ProposalQueue &q)
281 283
 {
282 284
     ar >> q.mRng >> q.mMinAtoms >> q.mMaxAtoms >> q.mBinLength >> q.mNumCols
283
-        >> q.mAlpha >> q.mDomainLength >> q.mNumBins;
285
+        >> q.mAlpha >> q.mDomainLength >> q.mNumBins >> q.mLambda
286
+        >> q.mUseCachedRng >> q.mU1 >> q.mU2;
284 287
     return ar;
285 288
 }
286 289
\ No newline at end of file
... ...
@@ -53,7 +53,9 @@ public:
53 53
     friend Archive& operator<<(Archive &ar, const ProposalQueue &queue);
54 54
     friend Archive& operator>>(Archive &ar, ProposalQueue &queue);
55 55
 
56
+#ifndef GAPS_INTERNAL_TESTS
56 57
 private:
58
+#endif
57 59
 
58 60
     std::vector<AtomicProposal> mQueue; // not really a queue for now
59 61
     
... ...
@@ -1,38 +1,42 @@
1 1
 #include "catch.h"
2 2
 #include "../utils/Archive.h"
3 3
 #include "../data_structures/Matrix.h"
4
-#include "../GibbsSampler.h"
5 4
 #include "../math/Random.h"
5
+#include "../atomic/AtomicDomain.h"
6
+#include "../atomic/ProposalQueue.h"
6 7
 
7
-TEST_CASE("Test utils/Archive.h")
8
-{
9
-    GapsRandomState randState(123);
10
-    GapsRng rng(&randState);
8
+// put Archive in it's own scope so it gets destructed (file stream closed)
11 9
 
12
-    SECTION("Reading/Writing to an Archive")
10
+TEST_CASE("Reading/Writing to an Archive")
11
+{
13 12
     {
14 13
         Archive ar1("test_ar.temp", ARCHIVE_WRITE);
15 14
         ar1 << 3;
16
-        ar1.close();
15
+    }
17 16
 
17
+    {
18 18
         Archive ar2("test_ar.temp", ARCHIVE_READ);
19 19
         unsigned i = 0;
20 20
         ar2 >> i;
21 21
         REQUIRE(i == 3);
22
-        ar2.close();
23 22
     }
24 23
 
25
-    SECTION("Serialization of primitive types")
26
-    {
27
-        // test values
28
-        unsigned u_read = 0, u_write = 456;
29
-        uint32_t u32_read = 0, u32_write = 512;
30
-        uint64_t u64_read = 0, u64_write = 0xAABBCCDDEE;
31
-        float f_read = 0.f, f_write = 0.123542f;
32
-        double d_read = 0., d_write = 0.54362;
33
-        bool b_read = false, b_write = true;
24
+    // cleanup directory
25
+    std::remove("test_ar.temp");
26
+}
27
+
28
+TEST_CASE("Serialization of primitive types")
29
+{
30
+    // test values
31
+    unsigned u_read = 0, u_write = 456;
32
+    uint32_t u32_read = 0, u32_write = 512;
33
+    uint64_t u64_read = 0, u64_write = 0xAABBCCDDEE;
34
+    float f_read = 0.f, f_write = 0.123542f;
35
+    double d_read = 0., d_write = 0.54362;
36
+    bool b_read = false, b_write = true;
34 37
 
35
-        // write to archive
38
+    // write to archive
39
+    {
36 40
         Archive arWrite("test_ar.temp", ARCHIVE_WRITE);
37 41
         arWrite << u_write;
38 42
         arWrite << u32_write;
... ...
@@ -40,9 +44,10 @@ TEST_CASE("Test utils/Archive.h")
40 44
         arWrite << f_write;
41 45
         arWrite << d_write;
42 46
         arWrite << b_write;
43
-        arWrite.close();
47
+    }
44 48
 
45
-        // read from archive
49
+    // read from archive
50
+    {
46 51
         Archive arRead("test_ar.temp", ARCHIVE_READ);
47 52
         arRead >> u_read;
48 53
         arRead >> u32_read;
... ...
@@ -50,139 +55,387 @@ TEST_CASE("Test utils/Archive.h")
50 55
         arRead >> f_read;
51 56
         arRead >> d_read;
52 57
         arRead >> b_read;
53
-        arRead.close();
54
-
55
-        // test that values are the same
56
-        REQUIRE(u_read == u_write);
57
-        REQUIRE(u32_read == u32_write);
58
-        REQUIRE(u64_read == u64_write);
59
-        REQUIRE(f_read == f_write);
60
-        REQUIRE(d_read == d_write);
61
-        REQUIRE(b_read == b_write);
62 58
     }
63 59
 
64
-    SECTION("Vector Serialization")
60
+    // test that values are the same
61
+    REQUIRE(u_read == u_write);
62
+    REQUIRE(u32_read == u32_write);
63
+    REQUIRE(u64_read == u64_write);
64
+    REQUIRE(f_read == f_write);
65
+    REQUIRE(d_read == d_write);
66
+    REQUIRE(b_read == b_write);
67
+
68
+    // cleanup directory
69
+    std::remove("test_ar.temp");
70
+}
71
+
72
+TEST_CASE("Vector Serialization")
73
+{
74
+    GapsRandomState randState(123);
75
+    GapsRng rng(&randState);
76
+
77
+    Vector vec_read(100), vec_write(100);
78
+    for (unsigned i = 0; i < 100; ++i)
65 79
     {
66
-        Vector vec_read(100), vec_write(100);
67
-        for (unsigned i = 0; i < 100; ++i)
68
-        {
69
-            vec_write[i] = rng.uniform(0.f, 2.f);
70
-        }
80
+        vec_write[i] = rng.uniform(0.f, 2.f);
81
+    }
71 82
 
83
+    {
72 84
         Archive arWrite("test_ar.temp", ARCHIVE_WRITE);
73 85
         arWrite << vec_write;
74
-        arWrite.close();
86
+    }
87
+
75 88
 
89
+    {
76 90
         Archive arRead("test_ar.temp", ARCHIVE_READ);
77 91
         arRead >> vec_read;
78
-        arRead.close();
92
+    }
79 93
 
80
-        REQUIRE(vec_read.size() == vec_write.size());
94
+    REQUIRE(vec_read.size() == vec_write.size());
81 95
 
82
-        for (unsigned i = 0; i < 100; ++i)
83
-        {
84
-            REQUIRE(vec_read[i] == vec_write[i]);
85
-        }
96
+    for (unsigned i = 0; i < 100; ++i)
97
+    {
98
+        REQUIRE(vec_read[i] == vec_write[i]);
86 99
     }
87 100
 
88
-    SECTION("Matrix Serialization")
89
-    {
90
-        ColMatrix cMat_read(100,100), cMat_write(100,100);
101
+    // cleanup directory
102
+    std::remove("test_ar.temp");
103
+}
104
+
105
+TEST_CASE("HybridVector Serialization")
106
+{
91 107
 
92
-        for (unsigned i = 0; i < 100; ++i)
108
+}
109
+
110
+TEST_CASE("SparseVector Serialization")
111
+{
112
+
113
+}
114
+
115
+TEST_CASE("Matrix Serialization")
116
+{
117
+    GapsRandomState randState(123);
118
+    GapsRng rng(&randState);
119
+
120
+    Matrix matRead(100,100), matWrite(100,100);
121
+
122
+    for (unsigned i = 0; i < 100; ++i)
123
+    {
124
+        for (unsigned j = 0; j < 100; ++j)
93 125
         {
94
-            for (unsigned j = 0; j < 100; ++j)
95
-            {
96
-                cMat_write(i,j) = rng.uniform(0.f, 2.f);
97
-            }
126
+            matWrite(i,j) = rng.uniform(0.f, 2.f);
98 127
         }
128
+    }
99 129
 
130
+    {
100 131
         Archive arWrite("test_ar.temp", ARCHIVE_WRITE);
101
-        arWrite << cMat_write;
102
-        arWrite.close();
132
+        arWrite << matWrite;
133
+    }
103 134
 
135
+    {
104 136
         Archive arRead("test_ar.temp", ARCHIVE_READ);
105
-        arRead >> cMat_read;
106
-        arRead.close();
137
+        arRead >> matRead;
138
+    }
107 139
 
108
-        REQUIRE(cMat_read.nRow() == cMat_write.nRow());
109
-        REQUIRE(cMat_read.nCol() == cMat_write.nCol());
110
-    
111
-        for (unsigned i = 0; i < 100; ++i)
140
+    REQUIRE(matRead.nRow() == matWrite.nRow());
141
+    REQUIRE(matRead.nCol() == matWrite.nCol());
142
+
143
+    for (unsigned i = 0; i < 100; ++i)
144
+    {
145
+        for (unsigned j = 0; j < 100; ++j)
112 146
         {
113
-            for (unsigned j = 0; j < 100; ++j)
114
-            {
115
-                REQUIRE(cMat_read(i,j) == cMat_write(i,j));
116
-            }
147
+            REQUIRE(matRead(i,j) == matWrite(i,j));
117 148
         }
118 149
     }
119 150
 
120
-    SECTION("GibbsSampler Serialization")
121
-    {
122
-        Rcpp::Environment env = Rcpp::Environment::global_env();
123
-        std::string csvPath = Rcpp::as<std::string>(env["gistCsvPath"]);
151
+    // cleanup directory
152
+    std::remove("test_ar.temp");
153
+}
154
+
155
+TEST_CASE("HybridMatrix Serialization")
156
+{
157
+
158
+}
159
+
160
+TEST_CASE("SparseMatrix Serialization")
161
+{
162
+
163
+}
164
+
165
+TEST_CASE("Random Generator Serialization")
166
+{
167
+    std::vector<float> randSequence;
124 168
 
125
-        GibbsSampler Asampler(csvPath, false, 7, false, std::vector<unsigned>());
126
-        GibbsSampler Psampler(csvPath, true, 7, false, std::vector<unsigned>());
127
-        Asampler.sync(Psampler);
128
-        Psampler.sync(Asampler);
129
-        
130
-        Asampler.update(10000, 1);
169
+    GapsRandomState randStateWrite(123);
170
+    GapsRng rngWrite(&randStateWrite);
131 171
 
172
+    volatile float burn_in = 0.0;
173
+    for (unsigned i = 0; i < 1000; ++i)
174
+    {
175
+        burn_in = rngWrite.uniform(0,1);
176
+    }
177
+    REQUIRE(burn_in < 1);
178
+
179
+    {
132 180
         Archive arWrite("test_ar.temp", ARCHIVE_WRITE);
133
-        arWrite << Asampler;
134
-        arWrite.close();
181
+        arWrite << rngWrite;
182
+    }
183
+    
184
+    for (unsigned i = 0; i < 1000; ++i)
185
+    {
186
+        randSequence.push_back(rngWrite.uniform());
187
+        randSequence.push_back(rngWrite.uniform(0.3, 6.4));
188
+        randSequence.push_back(rngWrite.exponential(5.5));
189
+    }
190
+
191
+    GapsRandomState randStateRead(456);
192
+    GapsRng rngRead(&randStateRead);
135 193
 
136
-        GibbsSampler savedAsampler(csvPath, false, 7, false, std::vector<unsigned>());
194
+    {
137 195
         Archive arRead("test_ar.temp", ARCHIVE_READ);
138
-        arRead >> savedAsampler;
139
-        arRead.close();
196
+        arRead >> rngRead;
140 197
     }
141 198
 
142
-#ifdef GAPS_INTERNAL_TESTS    
143
-    SECTION("AtomicDomain Serialization")
199
+    for (unsigned i = 0; i < 1000; ++i)
144 200
     {
145
-
201
+        REQUIRE(rngRead.uniform() == randSequence[i++]);
202
+        REQUIRE(rngRead.uniform(0.3, 6.4) == randSequence[i++]);
203
+        REQUIRE(rngRead.exponential(5.5) == randSequence[i]);
146 204
     }
205
+
206
+    // cleanup directory
207
+    std::remove("test_ar.temp");
208
+}
209
+
210
+TEST_CASE("GibbsSampler Serialization")
211
+{
212
+#if 0
213
+    Rcpp::Environment env = Rcpp::Environment::global_env();
214
+    std::string csvPath = Rcpp::as<std::string>(env["gistCsvPath"]);
215
+
216
+    GibbsSampler Asampler(csvPath, false, 7, false, std::vector<unsigned>());
217
+    GibbsSampler Psampler(csvPath, true, 7, false, std::vector<unsigned>());
218
+    Asampler.sync(Psampler);
219
+    Psampler.sync(Asampler);
220
+    
221
+    Asampler.update(10000, 1);
222
+
223
+    Archive arWrite("test_ar.temp", ARCHIVE_WRITE);
224
+    arWrite << Asampler;
225
+    arWrite.close();
226
+
227
+    GibbsSampler savedAsampler(csvPath, false, 7, false, std::vector<unsigned>());
228
+    Archive arRead("test_ar.temp", ARCHIVE_READ);
229
+    arRead >> savedAsampler;
230
+    arRead.close();
231
+
232
+    // cleanup directory
233
+    std::remove("test_ar.temp");
147 234
 #endif
235
+}
236
+
237
+TEST_CASE("GapsParameters Serialization")
238
+{
148 239
 
149
-/*
150
-    SECTION("Random Generator Serialization")
240
+}
241
+
242
+TEST_CASE("GapsStatistics Serialization")
243
+{
244
+
245
+}
246
+
247
+#ifdef GAPS_INTERNAL_TESTS
248
+TEST_CASE("AtomicDomain Serialization")
249
+{
250
+    GapsRandomState randState(123);
251
+    GapsRng rng(&randState);
252
+    
253
+    AtomicDomain domainWrite(100000);
254
+
255
+    for (unsigned i = 0; i < 1000; ++i)
256
+    {
257
+        domainWrite.insert(rng.uniform64(), rng.uniform(0.f, 100.f));
258
+    }
259
+
260
+    {
261
+        Archive arWrite("test_ar.temp", ARCHIVE_WRITE);
262
+        arWrite << domainWrite;
263
+    }
264
+
265
+    AtomicDomain domainRead(1);
266
+    
267
+    {
268
+        Archive arRead("test_ar.temp", ARCHIVE_READ);
269
+        arRead >> domainRead;
270
+    }
271
+
272
+    REQUIRE(domainWrite.front()->pos == domainRead.front()->pos);
273
+    REQUIRE(domainWrite.front()->mass == domainRead.front()->mass);
274
+    REQUIRE(domainWrite.size() == domainRead.size());
275
+    REQUIRE(domainWrite.mDomainLength == domainRead.mDomainLength);
276
+
277
+    for (unsigned i = 0; i < domainWrite.size(); ++i)
151 278
     {
152
-        std::vector<float> randSequence;
279
+        REQUIRE(domainWrite.mAtoms[i]->pos == domainRead.mAtoms[i]->pos);
280
+        REQUIRE(domainWrite.mAtoms[i]->mass == domainRead.mAtoms[i]->mass);
281
+    }
282
+}
283
+#endif
153 284
 
154
-        gaps::random::setSeed(123);
155
-        volatile float burn_in = 0.0;
156
-        for (unsigned i = 0; i < 1000; ++i)
285
+TEST_CASE("ProposalQueue Serialization")
286
+{
287
+    const unsigned nGenes = 10000;
288
+    const unsigned nPatterns = 100;
289
+    const unsigned nIterations = 1000;
290
+
291
+    GapsRandomState randStateWrite(123);
292
+    AtomicDomain domainWrite(nGenes * nPatterns);
293
+    ProposalQueue queueWrite(nGenes, nPatterns, &randStateWrite);
294
+    queueWrite.setAlpha(0.01f);
295
+    queueWrite.setLambda(0.01f);
296
+
297
+    for (unsigned i = 0; i < nIterations; ++i)
298
+    {
299
+        queueWrite.populate(domainWrite, nIterations);
300
+        for (unsigned j = 0; j < queueWrite.size(); ++j)
157 301
         {
158
-            burn_in = gaps::random::uniform(0,1);
302
+            switch (queueWrite[j].type)
303
+            {
304
+                case 'B':
305
+                    queueWrite.acceptBirth();
306
+                    queueWrite[j].atom1->mass = 3.f;
307
+                    break;
308
+                case 'D':
309
+                    queueWrite.acceptDeath();
310
+                    domainWrite.erase(queueWrite[j].atom1->pos);
311
+                    break;
312
+                case 'M':
313
+                    queueWrite[j].atom1->pos = queueWrite[j].pos;
314
+                    break;
315
+                case 'E':
316
+                    float mass1 = queueWrite[j].atom1->mass;
317
+                    queueWrite[j].atom1->mass = queueWrite[j].atom2->mass;
318
+                    queueWrite[j].atom2->mass = mass1;
319
+                    break;
320
+            }
159 321
         }
160
-        REQUIRE(burn_in < 1);
322
+        queueWrite.clear();
323
+    }
161 324
 
325
+    {
162 326
         Archive arWrite("test_ar.temp", ARCHIVE_WRITE);
163
-        gaps::random::save(arWrite);
164
-        arWrite.close();
327
+        arWrite << randStateWrite << queueWrite << domainWrite;
328
+    }
165 329
 
166
-        for (unsigned i = 0; i < 1000; ++i)
167
-        {
168
-            randSequence.push_back(gaps::random::uniform());
169
-            randSequence.push_back(gaps::random::uniform(0.3, 6.4));
170
-            randSequence.push_back(gaps::random::exponential(5.5));
171
-        }
330
+    GapsRandomState randStateRead(456);
331
+    AtomicDomain domainRead(nGenes * nPatterns);
332
+    ProposalQueue queueRead(nGenes, nPatterns, &randStateRead);
333
+    queueRead.setAlpha(100.f);
334
+    queueRead.setLambda(100.f);
172 335
     
173
-        gaps::random::setSeed(11);
336
+    {
174 337
         Archive arRead("test_ar.temp", ARCHIVE_READ);
175
-        gaps::random::load(arRead);
176
-        arRead.close();
338
+        arRead >> randStateRead >> queueRead >> domainRead;
339
+    }
340
+
341
+    GapsRng rngWriteTest(&randStateWrite);
342
+    GapsRng rngReadTest(&randStateRead);
343
+
344
+    REQUIRE(rngWriteTest.uniform() == rngReadTest.uniform());
345
+    REQUIRE(rngWriteTest.uniform() == rngReadTest.uniform());
346
+    REQUIRE(domainWrite.size() == domainRead.size());
347
+    REQUIRE(queueWrite.size() == 0);
348
+    REQUIRE(queueRead.size() == 0);
349
+
350
+#ifdef GAPS_INTERNAL_TESTS
351
+    REQUIRE(queueWrite.mRng.uniform() == queueRead.mRng.uniform());
352
+    REQUIRE(queueWrite.mRng.uniform() == queueRead.mRng.uniform());
353
+    REQUIRE(queueWrite.mRng.uniform() == queueRead.mRng.uniform());
354
+    REQUIRE(queueWrite.mRng.uniform() == queueRead.mRng.uniform());
355
+    REQUIRE(queueWrite.mMinAtoms == queueRead.mMinAtoms);
356
+    REQUIRE(queueWrite.mMaxAtoms == queueRead.mMaxAtoms);
357
+    REQUIRE(queueWrite.mBinLength == queueRead.mBinLength);
358
+    REQUIRE(queueWrite.mNumCols == queueRead.mNumCols);
359
+    REQUIRE(queueWrite.mAlpha == queueRead.mAlpha);
360
+    REQUIRE(queueWrite.mDomainLength == queueRead.mDomainLength);
361
+    REQUIRE(queueWrite.mNumBins == queueRead.mNumBins);
362
+    REQUIRE(queueWrite.mUseCachedRng == queueRead.mUseCachedRng);
363
+    for (unsigned i = 0; i < domainWrite.size(); ++i)
364
+    {
365
+        REQUIRE(domainWrite.mAtoms[i]->pos == domainRead.mAtoms[i]->pos);
366
+        REQUIRE(domainWrite.mAtoms[i]->mass == domainRead.mAtoms[i]->mass);
367
+    }
368
+#endif
369
+
370
+    for (unsigned i = 0; i < nIterations; ++i)
371
+    {
372
+        queueWrite.populate(domainWrite, nIterations);
373
+        queueRead.populate(domainRead, nIterations);
374
+        REQUIRE(queueWrite.size() == queueRead.size());
177 375
 
178
-        for (unsigned i = 0; i < 1000; ++i)
376
+        for (unsigned j = 0; j < queueWrite.size(); ++j)
179 377
         {
180
-            REQUIRE(gaps::random::uniform() == randSequence[i++]);
181
-            REQUIRE(gaps::random::uniform(0.3, 6.4) == randSequence[i++]);
182
-            REQUIRE(gaps::random::exponential(5.5) == randSequence[i]);
378
+            REQUIRE(queueWrite[j].type == queueRead[j].type);
379
+            REQUIRE(queueWrite[j].pos == queueRead[j].pos);
380
+            REQUIRE(queueWrite[j].r1 == queueRead[j].r1);
381
+            REQUIRE(queueWrite[j].c1 == queueRead[j].c1);
382
+            REQUIRE(queueWrite[j].r2 == queueRead[j].r2);
383
+            REQUIRE(queueWrite[j].c2 == queueRead[j].c2);
384
+
385
+            if (queueWrite[j].atom1 == NULL)
386
+            {
387
+                REQUIRE(queueRead[j].atom1 == NULL);
388
+            }
389
+            else
390
+            {
391
+                REQUIRE(queueWrite[j].atom1->pos == queueRead[j].atom1->pos);
392
+                REQUIRE(queueWrite[j].atom1->mass == queueRead[j].atom1->mass);
393
+            }
394
+
395
+            if (queueWrite[j].atom2 == NULL)
396
+            {
397
+                REQUIRE(queueRead[j].atom2 == NULL);
398
+            }
399
+            else
400
+            {
401
+                REQUIRE(queueWrite[j].atom2->pos == queueRead[j].atom2->pos);
402
+                REQUIRE(queueWrite[j].atom2->mass == queueRead[j].atom2->mass);
403
+            }
404
+
405
+            // process proposal
406
+            switch (queueWrite[j].type)
407
+            {
408
+                case 'B':
409
+                    queueWrite.acceptBirth();
410
+                    queueRead.acceptBirth();
411
+                    queueWrite[j].atom1->mass = 3.f;
412
+                    queueRead[j].atom1->mass = 3.f;
413
+                    break;
414
+                case 'D':
415
+                    queueWrite.acceptDeath();
416
+                    queueRead.acceptDeath();
417
+                    domainWrite.erase(queueWrite[j].atom1->pos);
418
+                    domainRead.erase(queueRead[j].atom1->pos);
419
+                    break;
420
+                case 'M':
421
+                    queueWrite[j].atom1->pos = queueWrite[j].pos;
422
+                    queueRead[j].atom1->pos = queueRead[j].pos;
423
+                    break;
424
+                case 'E':
425
+                    float mass1w = queueWrite[j].atom1->mass;
426
+                    queueWrite[j].atom1->mass = queueWrite[j].atom2->mass;
427
+                    queueWrite[j].atom2->mass = mass1w;
428
+
429
+                    float mass1r = queueRead[j].atom1->mass;
430
+                    queueRead[j].atom1->mass = queueRead[j].atom2->mass;
431
+                    queueRead[j].atom2->mass = mass1r;
432
+                    break;
433
+            }
183 434
         }
184
-    }
185
-*/
186
-    // cleanup directory
187
-    std::remove("test_ar.temp");
188
-}
189 435
\ No newline at end of file
436
+
437
+        queueWrite.clear();
438
+        queueRead.clear();
439
+    }    
440
+}
441
+
442
+
... ...
@@ -21,7 +21,9 @@ public:
21 21
     bool contains(unsigned n);
22 22
     bool isEmpty();
23 23
 
24
+#ifndef GAPS_INTERNAL_TESTS
24 25
 private:
26
+#endif
25 27
 
26 28
     std::vector<uint32_t> mSet;
27 29
     uint64_t mCurrentKey;
... ...
@@ -38,7 +40,9 @@ public:
38 40
     bool contains(uint64_t pos);
39 41
     bool isEmpty();
40 42
 
43
+#ifndef GAPS_INTERNAL_TESTS
41 44
 private:
45
+#endif
42 46
 
43 47
     std::vector<uint64_t> mSet;
44 48
 };
... ...
@@ -63,7 +67,9 @@ public:
63 67
     bool overlap(uint64_t pos); // this position in between pair
64 68
     bool isEmpty();
65 69
 
70
+#ifndef GAPS_INTERNAL_TESTS
66 71
 private:
72
+#endif
67 73
 
68 74
     std::vector<PositionPair> mSet;
69 75
 };
... ...
@@ -34,7 +34,9 @@ public:
34 34
     friend Archive& operator<<(Archive &ar, const HybridMatrix &vec);
35 35
     friend Archive& operator>>(Archive &ar, HybridMatrix &vec);
36 36
 
37
+#ifndef GAPS_INTERNAL_TESTS
37 38
 private:
39
+#endif
38 40
 
39 41
     std::vector<Vector> mRows;
40 42
     std::vector<HybridVector> mCols;
... ...
@@ -2,20 +2,27 @@
2 2
 #include "../math/Math.h"
3 3
 #include "../math/SIMD.h"
4 4
 
5
+#define SIMD_PAD(x) (gaps::simd::Index::increment() + \
6
+    gaps::simd::Index::increment() * ((x) / gaps::simd::Index::increment())) 
7
+
5 8
 HybridVector::HybridVector(unsigned sz)
6 9
     :
7 10
 mIndexBitFlags(sz / 64 + 1, 0),
8
-mData(sz, 0.f),
11
+mData(SIMD_PAD(sz), 0.f),
9 12
 mSize(sz)
10
-{}
13
+{
14
+    GAPS_ASSERT(mData.size() % gaps::simd::Index::increment() == 0);
15
+}
11 16
 
12 17
 HybridVector::HybridVector(const std::vector<float> &v)
13 18
     :
14 19
 mIndexBitFlags(v.size() / 64 + 1, 0),
15
-mData(v.size(), 0.f),
20
+mData(SIMD_PAD(v.size()), 0.f),
16 21
 mSize(v.size())
17 22
 {
18
-    for (unsigned i = 0; i < v.size(); ++i)
23
+    GAPS_ASSERT(mData.size() % gaps::simd::Index::increment() == 0);
24
+
25
+    for (unsigned i = 0; i < mSize; ++i)
19 26
     {
20 27
         mData[i] = v[i];
21 28
         if (v[i] > 0.f)
... ...
@@ -60,6 +67,7 @@ bool HybridVector::add(unsigned i, float v)
60 67
 
61 68
 float HybridVector::operator[](unsigned i) const
62 69
 {
70
+    GAPS_ASSERT(i < mSize);
63 71
     return mData[i];
64 72
 }
65 73
 
... ...
@@ -39,7 +39,9 @@ public:
39 39
     friend Archive& operator<<(Archive &ar, const HybridVector &vec);
40 40
     friend Archive& operator>>(Archive &ar, HybridVector &vec);
41 41
 
42
+#ifndef GAPS_INTERNAL_TESTS
42 43
 private:
44
+#endif
43 45
 
44 46
     std::vector<uint64_t> mIndexBitFlags;
45 47
     aligned_vector mData;
... ...
@@ -15,6 +15,14 @@ mNumRows(nrow),
15 15
 mNumCols(ncol)
16 16
 {}
17 17
 
18
+void Matrix::pad(float val)
19
+{
20
+    for (unsigned j = 0; j < mNumCols; ++j)
21
+    {
22
+        mCols[j].pad(val);
23
+    }
24
+}
25
+
18 26
 // constructor from data set read in as a matrix
19 27
 Matrix::Matrix(const Matrix &mat, bool genesInCols, bool subsetGenes,
20 28
 std::vector<unsigned> indices)
... ...
@@ -162,26 +170,26 @@ bool Matrix::empty() const
162 170
     return mNumRows == 0;
163 171
 }
164 172
 
165
-Archive& operator<<(Archive &ar, const Matrix &vec)
173
+Archive& operator<<(Archive &ar, const Matrix &mat)
166 174
 {
167
-    ar << vec.mNumRows << vec.mNumCols;
168
-    for (unsigned j = 0; j < vec.mNumCols; ++j)
175
+    ar << mat.mNumRows << mat.mNumCols;
176
+    for (unsigned j = 0; j < mat.mNumCols; ++j)
169 177
     {
170
-        ar << vec.mCols[j];
178
+        ar << mat.mCols[j];
171 179
     }
172 180
     return ar;
173 181
 }
174 182
 
175
-Archive& operator>>(Archive &ar, Matrix &vec)
183
+Archive& operator>>(Archive &ar, Matrix &mat)
176 184
 {
177 185
     unsigned nr = 0, nc = 0;
178 186
     ar >> nr >> nc;
179
-    GAPS_ASSERT(nr == vec.mNumRows);
180
-    GAPS_ASSERT(nc == vec.mNumCols);
187
+    GAPS_ASSERT(nr == mat.mNumRows);
188
+    GAPS_ASSERT(nc == mat.mNumCols);
181 189
 
182
-    for (unsigned j = 0; j < vec.mNumCols; ++j)
190
+    for (unsigned j = 0; j < mat.mNumCols; ++j)
183 191
     {
184
-        ar >> vec.mCols[j];
192
+        ar >> mat.mCols[j];
185 193
     }
186 194
     return ar;
187 195
 }
... ...
@@ -20,6 +20,8 @@ public:
20 20
     unsigned nRow() const;
21 21
     unsigned nCol() const;
22 22
 
23
+    void pad(float val);
24
+
23 25
     float operator()(unsigned i, unsigned j) const;
24 26
     float& operator()(unsigned i, unsigned j);
25 27
 
... ...
@@ -31,7 +33,9 @@ public:
31 33
     friend Archive& operator<<(Archive &ar, const Matrix &vec);
32 34
     friend Archive& operator>>(Archive &ar, Matrix &vec);
33 35
 
36
+#ifndef GAPS_INTERNAL_TESTS
34 37
 private:
38
+#endif
35 39
 
36 40
     std::vector<Vector> mCols;
37 41
     unsigned mNumRows;
... ...
@@ -25,7 +25,9 @@ public:
25 25
     bool atEnd() const;
26 26
     void next();
27 27
 
28
+#ifndef GAPS_INTERNAL_TESTS
28 29
 private:
30
+#endif
29 31
 
30 32
     friend float get<1>(const SparseIterator<1> &it);
31 33
 
... ...
@@ -46,7 +48,9 @@ public:
46 48
     void getFlags();
47 49
     unsigned getIndex() const;
48 50
 
51
+#ifndef GAPS_INTERNAL_TESTS
49 52
 private:
53
+#endif
50 54
 
51 55
     template <class Iter>
52 56
     friend void gotoNextCommon(Iter &it);
... ...
@@ -82,7 +86,9 @@ public:
82 86
     void getFlags();
83 87
     unsigned getIndex() const;
84 88
 
89
+#ifndef GAPS_INTERNAL_TESTS
85 90
 private:
91
+#endif
86 92
 
87 93
     template <class Iter>
88 94
     friend void gotoNextCommon(Iter &it);
... ...
@@ -29,7 +29,9 @@ public:
29 29
     friend Archive& operator<<(Archive &ar, const SparseMatrix &vec);
30 30
     friend Archive& operator>>(Archive &ar, SparseMatrix &vec);
31 31
 
32
+#ifndef GAPS_INTERNAL_TESTS
32 33
 private:
34
+#endif
33 35
 
34 36
     std::vector<SparseVector> mCols;
35 37
     unsigned mNumRows;
... ...
@@ -111,6 +111,7 @@ Archive& operator<<(Archive &ar, const SparseVector &vec)
111 111
 Archive& operator>>(Archive &ar, SparseVector &vec)
112 112
 {
113 113
     unsigned sz = 0;
114
+    ar >> sz;
114 115
     GAPS_ASSERT(sz == vec.mSize);
115 116
 
116 117
     for (unsigned i = 0; i < vec.mIndexBitFlags.size(); ++i)
... ...
@@ -1,32 +1,46 @@
1 1
 #include "Vector.h"
2 2
 #include "../math/SIMD.h"
3 3
 
4
+#define SIMD_PAD(x) (gaps::simd::Index::increment() + \
5
+    gaps::simd::Index::increment() * ((x) / gaps::simd::Index::increment())) 
6
+
4 7
 Vector::Vector(unsigned sz)
5 8
     :
6
-mData(sz, 0.f),
9
+mData(SIMD_PAD(sz), 0.f),
7 10
 mSize(sz)
8
-{}
11
+{
12
+    GAPS_ASSERT(mData.size() % gaps::simd::Index::increment() == 0);
13
+}
9 14
 
10 15
 Vector::Vector(const std::vector<float> &v)
11 16
     :
12
-mData(v.size(), 0.f),
17
+mData(SIMD_PAD(v.size()), 0.f),
13 18
 mSize(v.size())
14 19
 {
15
-    for (unsigned i = 0; i < v.size(); ++i)
20
+    GAPS_ASSERT(mData.size() % gaps::simd::Index::increment() == 0);
21
+    for (unsigned i = 0; i < mSize; ++i)
16 22
     {
17 23
         mData[i] = v[i];
18 24
     }
19 25
 }
20 26
 
27
+void Vector::pad(float val)
28
+{
29
+    for (unsigned i = mSize; i < mData.size(); ++i)
30
+    {
31
+        mData[i] = val;
32
+    }
33
+}
34
+
21 35
 float Vector::operator[](unsigned i) const
22 36
 {
23
-    GAPS_ASSERT(i < mData.size());
37
+    GAPS_ASSERT(i < mSize);
24 38
     return mData[i];
25 39
 }
26 40
 
27 41
 float& Vector::operator[](unsigned i)
28 42
 {
29
-    GAPS_ASSERT(i < mData.size());
43
+    GAPS_ASSERT(i < mSize);
30 44
     return mData[i];
31 45
 }
32 46
 
... ...
@@ -55,7 +69,7 @@ void Vector::operator+=(const Vector &v)
55 69
 
56 70
 void Vector::operator*=(float f)
57 71
 {
58
-    for (unsigned i = 0; i < mData.size(); ++i)
72
+    for (unsigned i = 0; i < mSize; ++i)
59 73
     {
60 74
         GAPS_ASSERT_MSG(mData[i] >= 0.f, mData[i]);
61 75
         mData[i] *= f;
... ...
@@ -64,15 +78,16 @@ void Vector::operator*=(float f)
64 78
 
65 79
 void Vector::operator/=(float f)
66 80
 {
67
-    for (unsigned i = 0; i < mData.size(); ++i)
81
+    for (unsigned i = 0; i < mSize; ++i)
68 82
     {
69
-        GAPS_ASSERT_MSG(mData[i] >= 0.f, i << " , " << mData.size() << " : " << mData[i]);
83
+        GAPS_ASSERT_MSG(mData[i] >= 0.f, i << " , " << mSize << " : " << mData[i]);
70 84
         mData[i] /= f;
71 85
     }
72 86
 }
73 87
 
74 88
 Archive& operator<<(Archive &ar, const Vector &vec)
75 89
 {
90
+    ar << vec.mSize;
76 91
     for (unsigned i = 0; i < vec.mSize; ++i)
77 92
     {
78 93
         ar << vec.mData[i];
... ...
@@ -82,6 +97,10 @@ Archive& operator<<(Archive &ar, const Vector &vec)
82 97
 
83 98
 Archive& operator>>(Archive &ar, Vector &vec)
84 99
 {
100
+    unsigned sz = 0;
101
+    ar >> sz;
102
+    GAPS_ASSERT(sz == vec.mSize);
103
+
85 104
     for (unsigned i = 0; i < vec.mSize; ++i)
86 105
     {
87 106
         ar >> vec.mData[i];
... ...
@@ -25,6 +25,7 @@ public:
25 25
     float* ptr();
26 26
 
27 27
     unsigned size() const;
28
+    void pad(float val);
28 29
 
29 30
     void operator+=(const Vector &v);
30 31
 
... ...
@@ -34,7 +35,9 @@ public:
34 35
     friend Archive& operator<<(Archive &ar, const Vector &vec);
35 36
     friend Archive& operator>>(Archive &ar, Vector &vec);
36 37
 
38
+#ifndef GAPS_INTERNAL_TESTS
37 39
 private:
40
+#endif
38 41
 
39 42
     aligned_vector mData;
40 43
     unsigned mSize;
... ...
@@ -9,7 +9,20 @@
9 9
 
10 10
 class CsvParser : public AbstractFileParser
11 11
 {
12
+public:
13
+
14
+    explicit CsvParser(const std::string &path);
15
+    ~CsvParser() { mFile.close(); }
16
+
17
+    unsigned nRow() const { return mNumRows; }
18
+    unsigned nCol() const { return mNumCols; }
19
+
20
+    bool hasNext();
21
+    MatrixElement getNext();
22
+
23
+#ifndef GAPS_INTERNAL_TESTS
12 24
 private:
25
+#endif
13 26
 
14 27
     std::ifstream mFile;
15 28
 
... ...
@@ -21,17 +34,6 @@ private:
21 34
 
22 35
     CsvParser(const CsvParser &p); // don't allow copies
23 36
     CsvParser& operator=(const CsvParser &p); // don't allow copies
24
-
25
-public:
26
-
27
-    explicit CsvParser(const std::string &path);
28
-    ~CsvParser() { mFile.close(); }
29
-
30
-    unsigned nRow() const { return mNumRows; }
31
-    unsigned nCol() const { return mNumCols; }
32
-
33
-    bool hasNext();
34
-    MatrixElement getNext();
35 37
 };
36 38
 
37 39
 #endif
38 40
\ No newline at end of file
... ...
@@ -18,11 +18,6 @@ enum GapsFileType
18 18
 // file parser interface
19 19
 class AbstractFileParser
20 20
 {
21
-private:
22
-
23
-    AbstractFileParser(const AbstractFileParser &p); // don't allow copies
24
-    AbstractFileParser& operator=(const AbstractFileParser &p); // don't allow copies
25
-
26 21
 public:
27 22
 
28 23
     AbstractFileParser() {}
... ...
@@ -36,18 +31,16 @@ public:
36 31
 
37 32
     virtual bool hasNext() = 0;
38 33
     virtual MatrixElement getNext() = 0;
34
+
35
+private:
36
+
37
+    AbstractFileParser(const AbstractFileParser &p); // don't allow copies
38
+    AbstractFileParser& operator=(const AbstractFileParser &p); // don't allow copies
39 39
 };
40 40
 
41 41
 // wrap the pointer to the parser implementation
42 42
 class FileParser
43 43
 {
44
-private:
45
-
46
-    AbstractFileParser *mParser;
47
-
48
-    FileParser(const FileParser &p); // don't allow copies
49
-    FileParser& operator=(const FileParser &p); // don't allow copies
50
-
51 44
 public:
52 45
 
53 46
     explicit FileParser(const std::string &path);
... ...
@@ -72,6 +65,15 @@ public:
72 65
 
73 66
     template <class MatrixType>
74 67
     static void writeToGct(const std::string &path, const MatrixType &mat);
68
+
69
+#ifndef GAPS_INTERNAL_TESTS
70
+private:
71
+#endif
72
+
73
+    AbstractFileParser *mParser;
74
+
75
+    FileParser(const FileParser &p); // don't allow copies
76
+    FileParser& operator=(const FileParser &p); // don't allow copies
75 77
 };
76 78
 
77 79
 // temporary solution - should be moved into specific file parsers, ok for now
... ...
@@ -9,7 +9,20 @@
9 9
 
10 10
 class GctParser : public AbstractFileParser
11 11
 {
12
+public:
13
+
14
+    explicit GctParser(const std::string &path);
15
+    ~GctParser() { mFile.close(); }
16
+
17
+    unsigned nRow() const { return mNumRows; }
18
+    unsigned nCol() const { return mNumCols; }
19
+
20
+    bool hasNext();
21
+    MatrixElement getNext();
22
+
23
+#ifndef GAPS_INTERNAL_TESTS
12 24
 private:
25
+#endif
13 26
 
14 27
     std::ifstream mFile;
15 28
 
... ...
@@ -21,17 +34,6 @@ private:
21 34
 
22 35
     GctParser(const GctParser &p); // don't allow copies
23 36
     GctParser& operator=(const GctParser &p); // don't allow copies
24
-
25
-public:
26
-
27
-    explicit GctParser(const std::string &path);
28
-    ~GctParser() { mFile.close(); }
29
-
30
-    unsigned nRow() const { return mNumRows; }
31
-    unsigned nCol() const { return mNumCols; }
32
-
33
-    bool hasNext();
34
-    MatrixElement getNext();
35 37
 };
36 38
 
37 39
 #endif
38 40
\ No newline at end of file
... ...
@@ -9,16 +9,6 @@
9 9
 
10 10
 class MtxParser : public AbstractFileParser
11 11
 {
12
-private:
13
-
14
-    std::ifstream mFile;
15
-
16
-    unsigned mNumRows;
17
-    unsigned mNumCols;
18
-
19
-    MtxParser(const MtxParser &p); // don't allow copies
20
-    MtxParser& operator=(const MtxParser &p); // don't allow copies
21
-
22 12
 public:
23 13
 
24 14
     explicit MtxParser(const std::string &path);
... ...
@@ -29,6 +19,18 @@ public:
29 19
 
30 20
     bool hasNext();
31 21
     MatrixElement getNext();
22
+
23
+#ifndef GAPS_INTERNAL_TESTS
24
+private:
25
+#endif
26
+
27
+    std::ifstream mFile;
28
+
29
+    unsigned mNumRows;
30
+    unsigned mNumCols;
31
+
32
+    MtxParser(const MtxParser &p); // don't allow copies
33
+    MtxParser& operator=(const MtxParser &p); // don't allow copies
32 34
 };
33 35
 
34 36
 #endif
35 37
\ No newline at end of file
... ...
@@ -9,7 +9,20 @@
9 9
 
10 10
 class TsvParser : public AbstractFileParser
11 11
 {
12
+public:
13
+
14
+    explicit TsvParser(const std::string &path);
15
+    ~TsvParser() { mFile.close(); }
16
+
17
+    unsigned nRow() const { return mNumRows; }
18
+    unsigned nCol() const { return mNumCols; }
19
+
20
+    bool hasNext();
21
+    MatrixElement getNext();
22
+
23
+#ifndef GAPS_INTERNAL_TESTS
12 24
 private:
25
+#endif
13 26
 
14 27
     std::ifstream mFile;
15 28
 
... ...
@@ -21,17 +34,6 @@ private:
21 34
 
22 35
     TsvParser(const TsvParser &p); // don't allow copies
23 36
     TsvParser& operator=(const TsvParser &p); // don't allow copies
24
-
25
-public:
26
-
27
-    explicit TsvParser(const std::string &path);
28
-    ~TsvParser() { mFile.close(); }
29
-
30
-    unsigned nRow() const { return mNumRows; }
31
-    unsigned nCol() const { return mNumCols; }
32
-
33
-    bool hasNext();
34
-    MatrixElement getNext();
35 37
 };
36 38
 
37 39
 #endif
38 40
\ No newline at end of file
... ...
@@ -85,8 +85,7 @@ AlphaParameters DenseGibbsSampler::alphaParameters(unsigned row, unsigned col)
85 85
 
86 86
     gaps::simd::PackedFloat pMat, pD, pAP, pS;
87 87
     gaps::simd::PackedFloat partialS(0.f), partialS_mu(0.f);
88
-    gaps::simd::Index i(0);
89
-    for (; i <= size - gaps::simd::Index::increment(); ++i)
88
+    for (gaps::simd::Index i(0); i < size; ++i)
90 89
     {   
91 90
         pMat.load(mat + i);
92 91
         pD.load(D + i);
... ...
@@ -96,14 +95,9 @@ AlphaParameters DenseGibbsSampler::alphaParameters(unsigned row, unsigned col)
96 95
         partialS += pMat * ratio;
97 96
         partialS_mu += ratio * (pD - pAP);
98 97
     }
99
-
100
-    float s = partialS.scalar(), s_mu = partialS_mu.scalar();
101
-    for (unsigned j = i.value(); j < size; ++j)
102
-    {
103
-        float ratio = mat[j] / (S[j] * S[j]);
104
-        s += mat[j] * ratio;
105
-        s_mu += ratio * (D[j] - AP[j]);
106
-    }
98
+    
99
+    float s = partialS.scalar();
100
+    float s_mu = partialS_mu.scalar();
107 101
 
108 102
 #if 0 // optional debug section, set to 1 to enable, 0 to disable
109 103
 #ifdef GAPS_DEBUG
... ...
@@ -117,19 +111,12 @@ AlphaParameters DenseGibbsSampler::alphaParameters(unsigned row, unsigned col)
117 111
    
118 112
     float s_denom = debug_s < 10.f ? 10.f : debug_s;
119 113
     float smu_denom = std::abs(debug_smu) < 10.f ? 10.f : std::abs(debug_smu);
120
-    const float tolerance = 0.01f;
121 114
 
122
-    if (std::abs(s - debug_s) / s_denom >= tolerance)
123
-    {
124
-        gaps_printf("%f != %f\n", s, debug_s);
125
-    }
126
-    if (std::abs(s_mu - debug_smu) / smu_denom >= tolerance)
127
-    {
128
-        gaps_printf("%f != %f\n", s_mu, debug_smu);
129
-    }
130
-
131
-    GAPS_ASSERT(std::abs(s - debug_s) / s_denom < tolerance);
132
-    GAPS_ASSERT(std::abs(s_mu - debug_smu) / smu_denom < tolerance);
115
+    const float tolerance = 0.01f;
116
+    GAPS_ASSERT_MSG(std::abs(s - debug_s) / s_denom < tolerance,
117
+        s << " != " << debug_s);
118
+    GAPS_ASSERT_MSG(std::abs(s_mu - debug_smu) / smu_denom < tolerance,
119
+        s_mu << " != " << debug_smu);
133 120
 #endif
134 121
 #endif
135 122
 
... ...
@@ -140,7 +127,7 @@ AlphaParameters DenseGibbsSampler::alphaParameters(unsigned row, unsigned col)
140 127
 AlphaParameters DenseGibbsSampler::alphaParameters(unsigned r1, unsigned c1,
141 128
 unsigned r2, unsigned c2)
142 129
 {
143
-    if (r1 == r2) // expect false ?
130
+    if (r1 == r2)
144 131
     {
145 132
         unsigned size = mDMatrix.nRow();
146 133
         const float *D = mDMatrix.getCol(r1).ptr();
... ...
@@ -150,9 +137,8 @@ unsigned r2, unsigned c2)
150 137
         const float *mat2 = mOtherMatrix->getCol(c2).ptr();
151 138
 
152 139
         gaps::simd::PackedFloat pMat1, pMat2, pD, pAP, pS;
153
-        gaps::simd::PackedFloat partialS(0.f), partialS_mu(0.f);
154
-        gaps::simd::Index i(0);
155
-        for (; i <= size - gaps::simd::Index::increment(); ++i)
140
+        gaps::simd::PackedFloat packedS(0.f), packedS_mu(0.f);
141
+        for (gaps::simd::Index i(0); i < size; ++i)
156 142
         {   
157 143
             pMat1.load(mat1 + i);
158 144
             pMat2.load(mat2 + i);
... ...
@@ -160,18 +146,10 @@ unsigned r2, unsigned c2)
160 146
             pAP.load(AP + i);
161 147
             pS.load(S + i);
162 148
             gaps::simd::PackedFloat ratio((pMat1 - pMat2) / (pS * pS));
163
-            partialS += (pMat1 - pMat2) * ratio;
164
-            partialS_mu += ratio * (pD - pAP);
149
+            packedS += (pMat1 - pMat2) * ratio;
150
+            packedS_mu += ratio * (pD - pAP);
165 151
         }
166
-
167
-        float s = partialS.scalar(), s_mu = partialS_mu.scalar();
168
-        for (unsigned j = i.value(); j < size; ++j)
169
-        {
170
-            float ratio = (mat1[j] - mat2[j]) / (S[j] * S[j]);
171
-            s += (mat1[j] - mat2[j]) * ratio;
172
-            s_mu += ratio * (D[j] - AP[j]);
173
-        }
174
-        return AlphaParameters(s, s_mu);
152
+        return AlphaParameters(packedS.scalar(), packedS_mu.scalar());
175 153
     }
176 154
     return alphaParameters(r1, c1) + alphaParameters(r2, c2);
177 155
 }
... ...
@@ -188,28 +166,18 @@ unsigned col, float ch)
188 166
 
189 167
     gaps::simd::PackedFloat pCh(ch);
190 168
     gaps::simd::PackedFloat pMat, pD, pAP, pS;
191
-    gaps::simd::PackedFloat partialS(0.f), partialS_mu(0.f);
192
-    gaps::simd::Index i(0);
193
-    for (; i <= size - gaps::simd::Index::increment(); ++i)
169
+    gaps::simd::PackedFloat packedS(0.f), packedS_mu(0.f);
170
+    for (gaps::simd::Index i(0); i < size; ++i)
194 171
     {   
195 172
         pMat.load(mat + i);
196 173
         pD.load(D + i);
197 174
         pAP.load(AP + i);
198 175
         pS.load(S + i);
199 176
         gaps::simd::PackedFloat ratio(pMat / (pS * pS));
200
-        partialS += pMat * ratio;
201
-        partialS_mu += ratio * (pD - (pAP + pCh * pMat));
202
-    }
203
-
204
-    float s = partialS.scalar(), s_mu = partialS_mu.scalar();
205
-    for (unsigned j = i.value(); j < size; ++j)
206
-    {
207
-        float ratio = mat[j] / (S[j] * S[j]);
208
-        s += mat[j] * ratio;
209
-        s_mu += ratio * (D[j] - (AP[j] + ch * mat[j]));
177
+        packedS += pMat * ratio;
178
+        packedS_mu += ratio * (pD - (pAP + pCh * pMat));
210 179
     }
211
-    return AlphaParameters(s, s_mu);
212
-
180
+    return AlphaParameters(packedS.scalar(), packedS_mu.scalar());
213 181
 }
214 182
 
215 183
 // PERFORMANCE_CRITICAL
... ...
@@ -220,33 +188,29 @@ void DenseGibbsSampler::updateAPMatrix(unsigned row, unsigned col, float delta)
220 188
     unsigned size = mAPMatrix.nRow();
221 189
 
222 190
     gaps::simd::PackedFloat pOther, pAP;
223
-    gaps::simd::Index i(0);
224 191
     gaps::simd::PackedFloat pDelta(delta);
225
-    for (; i <= size - gaps::simd::Index::increment(); ++i)
192
+    for (gaps::simd::Index i(0); i < size; ++i)
226 193
     {
227 194
         pOther.load(other + i);
228 195
         pAP.load(ap + i);
229 196
         pAP += pDelta * pOther;
230 197
         pAP.store(ap + i);
231 198
     }
232
-
233
-    for (unsigned j = i.value(); j < size; ++j)
234
-    {
235
-        ap[j] += delta * other[j];
236
-    }
237 199
 }
238 200
 
239 201
 Archive& operator<<(Archive &ar, const DenseGibbsSampler &s)
240 202
 {
241
-    ar << s.mMatrix << s.mDomain << s.mAlpha << s.mLambda << s.mMaxGibbsMass
242
-        << s.mAnnealingTemp << s.mNumPatterns << s.mNumBins << s.mBinLength;
203
+    ar << s.mMatrix << s.mDomain << s.mQueue << s.mAlpha << s.mLambda
204
+        << s.mMaxGibbsMass << s.mAnnealingTemp << s.mNumPatterns << s.mNumBins
205
+        << s.mBinLength;
243 206
     return ar;
244 207
 }
245 208
 
246 209
 Archive& operator>>(Archive &ar, DenseGibbsSampler &s)
247 210
 {
248
-    ar >> s.mMatrix >> s.mDomain >> s.mAlpha >> s.mLambda >> s.mMaxGibbsMass
249
-        >> s.mAnnealingTemp >> s.mNumPatterns >> s.mNumBins >> s.mBinLength;
211
+    ar >> s.mMatrix >> s.mDomain >> s.mQueue >> s.mAlpha >> s.mLambda
212
+        >> s.mMaxGibbsMass >> s.mAnnealingTemp >> s.mNumPatterns >> s.mNumBins
213
+        >> s.mBinLength;
250 214
     return ar;
251 215
 }
252 216
 
... ...
@@ -29,6 +29,8 @@ public:
29 29
     friend Archive& operator<<(Archive &ar, const DenseGibbsSampler &s);
30 30
     friend Archive& operator>>(Archive &ar, DenseGibbsSampler &s);
31 31
 
32
+    float apSum() const { return gaps::sum(mAPMatrix); }
33
+
32 34
 #ifdef GAPS_DEBUG
33 35
     bool internallyConsistent() const;
34 36
 #endif
... ...
@@ -60,13 +62,16 @@ GapsRandomState *randState)
60 62
 GibbsSampler(data, transpose, subsetRows, alpha, maxGibbsMass, params, randState),
61 63
 mSMatrix(gaps::pmax(mDMatrix, 0.1f)),
62 64
 mAPMatrix(mDMatrix.nRow(), mDMatrix.nCol())
63
-{}
65
+{
66
+    mSMatrix.pad(1.f); // can't divide by zero
67
+}
64 68
 
65 69
 template <class DataType>
66 70
 void DenseGibbsSampler::setUncertainty(const DataType &unc, bool transpose,
67 71
 bool subsetRows, const GapsParameters &params)
68 72
 {
69 73
     mSMatrix = Matrix(unc, transpose, subsetRows, params.dataIndicesSubset);
74
+    mSMatrix.pad(1.f); // can't divide by zero
70 75
 }
71 76
 
72 77
 #endif // __COGAPS_DENSE_GIIBS_SAMPLER_H__
73 78
\ No newline at end of file
... ...
@@ -47,8 +47,7 @@ public:
47 47
 
48 48
     void update(unsigned nSteps, unsigned nThreads);
49 49
 
50
-//#ifndef GAPS_INTERNAL_TESTS
51
-#if 0
50
+#ifndef GAPS_INTERNAL_TESTS
52 51
 protected:
53 52
 #endif
54 53
 
... ...
@@ -242,8 +241,8 @@ void GibbsSampler<Derived, DataMatrix, FactorMatrix>::death(const AtomicProposal
242 241
         if (rebirthMass != prop.atom1->mass)
243 242
         {
244 243
             impl()->safelyChangeMatrix(prop.r1, prop.c1, rebirthMass - prop.atom1->mass);
244
+            prop.atom1->mass = rebirthMass;
245 245
         }
246
-        prop.atom1->mass = rebirthMass;
247 246
     }
248 247
     else // reject rebirth
249 248
     {
... ...
@@ -15,7 +15,8 @@ float SparseGibbsSampler::chiSq() const
15 15
         Vector S(gaps::pmax(D, 0.1f));
16 16
         for (unsigned i = 0; i < D.size(); ++i)
17 17
         {
18
-            float ap = gaps::dot(mMatrix.getRow(j), mOtherMatrix->getRow(i));
18
+            //float ap = gaps::dot(mMatrix.getRow(j), mOtherMatrix->getRow(i));
19
+            float ap = 0.f;
19 20
             chisq += GAPS_SQ(D[i] - ap) / GAPS_SQ(S[i]);
20 21
         }
21 22
     }
... ...
@@ -103,8 +104,7 @@ unsigned col, float ch)
103 104
     SparseIterator<2> it(mDMatrix.getCol(row), mOtherMatrix->getCol(col));
104 105
     while (!it.atEnd())
105 106
     {
106
-        s += (get<2>(it) * get<2>(it))
107
-            / (get<1>(it) * get<1>(it))
107
+        s += (get<2>(it) * get<2>(it)) / (get<1>(it) * get<1>(it))
108 108
             - (get<2>(it) * get<2>(it));
109 109
 
110 110
         s_mu += get<2>(it) / get<1>(it)
... ...
@@ -34,6 +34,8 @@ public:
34 34
     friend Archive& operator<<(Archive &ar, const SparseGibbsSampler &s);
35 35
     friend Archive& operator>>(Archive &ar, SparseGibbsSampler &s);
36 36
 
37
+    float apSum() const { return 0.f; }
38
+
37 39
 #ifdef GAPS_DEBUG
38 40
     bool internallyConsistent() const;
39 41
 #endif
... ...
@@ -42,7 +42,6 @@ uint32_t GapsRng::next()
42 42
 
43 43
 void GapsRng::advance()
44 44
 {
45
-    mPreviousState = mState;
46 45
     mState = mState * 6364136223846793005ull + (54u|1);
47 46
 }
48 47
 
... ...
@@ -53,11 +53,12 @@ public:
53 53
     friend Archive& operator<<(Archive &ar, const GapsRng &gen);
54 54
     friend Archive& operator>>(Archive &ar, GapsRng &gen);
55 55
 
56
+#ifndef GAPS_INTERNAL_TESTS
56 57
 private:
58
+#endif
57 59
 
58 60
     const GapsRandomState *mRandState;
59 61
     uint64_t mState;
60
-    uint32_t mPreviousState;
61 62
 
62 63
     uint32_t next();
63 64
     void advance();
... ...
@@ -71,7 +72,11 @@ private:
71 72
 // private class used by GapsRandomState for seeding individual rngs
72 73
 class Xoroshiro128plus
73 74
 {
75
+#ifndef GAPS_INTERNAL_TESTS
74 76
 private:
77
+#else
78
+public:
79
+#endif
75 80
 
76 81
     friend class GapsRandomState;
77 82
 
... ...
@@ -106,8 +111,10 @@ public:
106 111
 
107 112
     friend Archive& operator<<(Archive &ar, const GapsRandomState &s);
108 113
     friend Archive& operator>>(Archive &ar, GapsRandomState &s);
109
-    
114
+   
115
+#ifndef GAPS_INTERNAL_TESTS
110 116
 private:
117
+#endif
111 118
 
112 119
     friend GapsRng;
113 120
 
... ...
@@ -4,21 +4,14 @@
4 4
 
5 5
 static float dot_helper(const float *v1, const float *v2, unsigned size)
6 6
 {
7
-    gaps::simd::PackedFloat partialDot(0.f), p1, p2;
8
-    gaps::simd::Index i(0);
9
-    for (; i <= size - gaps::simd::Index::increment(); ++i)
7
+    gaps::simd::PackedFloat packedDot(0.f), p1, p2;
8
+    for (gaps::simd::Index i(0); i < size; ++i)
10 9
     {
11 10
         p1.load(v1 + i);
12 11
         p2.load(v2 + i);
13
-        partialDot += p1 * p2;
12
+        packedDot += p1 * p2;
14 13
     }
15
-
16
-    float dot = partialDot.scalar();
17
-    for (unsigned j = i.value(); j < size; ++j)
18
-    {
19
-        dot += v1[j] * v2[j];
20
-    }
21
-    return dot;
14
+    return packedDot.scalar();
22 15
 }
23 16
 
24 17
 float gaps::min(const Vector &v)
... ...
@@ -17,10 +17,6 @@
17 17
 
18 18
 class Archive
19 19
 {
20
-private:
21
-
22
-    std::fstream mStream;
23
-
24 20
 public:
25 21
 
26 22
     Archive(const std::string &path, std::ios_base::openmode flags)
... ...
@@ -83,6 +79,10 @@ public:
83 79
     friend Archive& operator>>(Archive &ar, int64_t &val)  { return readFromArchive(ar, val); }
84 80
     friend Archive& operator>>(Archive &ar, float &val)    { return readFromArchive(ar, val); }
85 81
     friend Archive& operator>>(Archive &ar, double &val)   { return readFromArchive(ar, val); }
82
+
83
+private:
84
+
85
+    std::fstream mStream;
86 86
 };
87 87
 
88 88
 #endif // __COGAPS_ARCHIVE_H__
89 89
new file mode 100644
90 90
Binary files /dev/null and b/tests/testthat/test_ar.temp differ
... ...
@@ -2,11 +2,13 @@ context("CoGAPS")
2 2
 
3 3
 test_that("Checkpoint System",
4 4
 {
5
-    #data(SimpSim)
6
-    #run1 <- CoGAPS(SimpSim.data, checkpointInterval=100,
7
-    #    checkpointOutFile="test.out", messages=FALSE)
8
-    #run2 <- CoGAPS(SimpSim.data, checkpointInFile="test.out", messages=FALSE)
5
+    data(SimpSim)
6
+    run1 <- CoGAPS(SimpSim.data, checkpointInterval=501,
7
+        checkpointOutFile="test.out", messages=FALSE)
8
+    run2 <- CoGAPS(SimpSim.data, checkpointInFile="test.out", messages=FALSE)
9 9
 
10
-    #expect_true(all.equal(run1@featureLoadings, run2@featureLoadings))
11
-    #expect_true(all.equal(run1@sampleFactors, run2@sampleFactors))
10
+    print(max(run1@featureLoadings - run2@featureLoadings))
11
+
12
+    expect_true(all.equal(run1@featureLoadings, run2@featureLoadings))
13
+    expect_true(all.equal(run1@sampleFactors, run2@sampleFactors))
12 14
 })
13 15
\ No newline at end of file
... ...
@@ -108,3 +108,4 @@ test_that("Valid Top-Level CoGAPS Calls",
108 108
 #{
109 109
 
110 110
 #})
111
+