Browse code

tidied up notation around pattern markers

Tom Sherman authored on 05/09/2019 18:14:53
Showing1 changed files
... ...
@@ -133,13 +133,13 @@ static double estimatedNumUpdates(double current, double total, float nAtoms)
133 133
 
134 134
 template <class Sampler>
135 135
 static double estimatedPercentComplete(const GapsParameters &params,
136
-const Sampler &ASampler, const Sampler &PSampler, char phase, unsigned iter)
136
+const Sampler &ASampler, const Sampler &PSampler, GapsAlgorithmPhase phase, unsigned iter)
137 137
 {
138 138
     double nIter = static_cast<double>(iter);
139 139
     double nAtomsA = static_cast<double>(ASampler.nAtoms());
140 140
     double nAtomsP = static_cast<double>(PSampler.nAtoms());
141 141
     
142
-    if (phase == 'S')
142
+    if (phase == GAPS_SAMPLING_PHASE)
143 143
     {
144 144
         nIter += params.nIterations;
145 145
     }
... ...
@@ -155,7 +155,7 @@ const Sampler &ASampler, const Sampler &PSampler, char phase, unsigned iter)
155 155
 template <class Sampler>
156 156
 static void displayStatus(const GapsParameters &params,
157 157
 const Sampler &ASampler, const Sampler &PSampler, bpt::ptime startTime,
158
-char phase, unsigned iter, GapsStatistics &stats)
158
+GapsAlgorithmPhase phase, unsigned iter, GapsStatistics &stats)
159 159
 {
160 160
     if (params.outputFrequency > 0 && ((iter + 1) % params.outputFrequency) == 0)
161 161
     {
... ...
@@ -208,7 +208,7 @@ Sampler &PSampler, unsigned nA, unsigned nP)
208 208
 template <class Sampler>
209 209
 static void createCheckpoint(const GapsParameters &params,
210 210
 Sampler &ASampler, Sampler &PSampler, const GapsRandomState *randState,
211
-const GapsStatistics &stats, const GapsRng &rng, char phase, unsigned iter)
211
+const GapsStatistics &stats, const GapsRng &rng, GapsAlgorithmPhase phase, unsigned iter)
212 212
 {
213 213
     if (params.checkpointInterval > 0 && ((iter + 1) % params.checkpointInterval) == 0
214 214
     && !params.subsetData)
Browse code

Snapshots working for both phases

Tom Sherman authored on 05/09/2019 00:20:30
Showing1 changed files
... ...
@@ -221,7 +221,7 @@ const GapsStatistics &stats, const GapsRng &rng, char phase, unsigned iter)
221 221
         Archive ar(params.checkpointOutFile, ARCHIVE_WRITE);
222 222
         ar << params;
223 223
         ar << *randState;
224
-        ar << ASampler << PSampler << stats << phase << iter << rng;
224
+        ar << ASampler << PSampler << stats << static_cast<int>(phase) << iter << rng;
225 225
         
226 226
         // delete backup file
227 227
         std::remove((params.checkpointOutFile + ".backup").c_str());
... ...
@@ -239,22 +239,24 @@ const GapsStatistics &stats, const GapsRng &rng, char phase, unsigned iter)
239 239
 template <class Sampler>
240 240
 static void processCheckpoint(GapsParameters &params, Sampler &ASampler,
241 241
 Sampler &PSampler, GapsRandomState *randState, GapsStatistics &stats,
242
-GapsRng &rng, char &phase, unsigned &currentIter)
242
+GapsRng &rng, GapsAlgorithmPhase &phase, unsigned &currentIter)
243 243
 {    
244 244
     // check if running from checkpoint, get all saved data
245 245
     if (params.useCheckPoint)
246 246
     {
247
+        int iPhase;
247 248
         Archive ar(params.checkpointFile, ARCHIVE_READ);
248 249
         ar >> params;
249 250
         ar >> *randState;
250
-        ar >> ASampler >> PSampler >> stats >> phase >> currentIter >> rng;
251
+        ar >> ASampler >> PSampler >> stats >> iPhase >> currentIter >> rng;
252
+        phase = static_cast<GapsAlgorithmPhase>(iPhase);
251 253
     }
252 254
 }
253 255
 
254 256
 template <class Sampler>
255 257
 static uint64_t runOnePhase(const GapsParameters &params, Sampler &ASampler,
256 258
 Sampler &PSampler, GapsStatistics &stats, const GapsRandomState *randState,
257
-GapsRng &rng, bpt::ptime startTime, char phase, unsigned &currentIter)
259
+GapsRng &rng, bpt::ptime startTime, GapsAlgorithmPhase phase, unsigned &currentIter)
258 260
 {
259 261
     uint64_t totalUpdates = 0;
260 262
     for (; currentIter < params.nIterations; ++currentIter)
... ...
@@ -264,7 +266,7 @@ GapsRng &rng, bpt::ptime startTime, char phase, unsigned &currentIter)
264 266
             rng, phase, currentIter);
265 267
         
266 268
         // set annealing temperature in Equilibration phase
267
-        if (phase == 'C')
269
+        if (phase == GAPS_EQUILIBRATION_PHASE)
268 270
         {        
269 271
             float temp = static_cast<float>(2 * currentIter)
270 272
                 / static_cast<float>(params.nIterations);
... ...
@@ -278,16 +280,19 @@ GapsRng &rng, bpt::ptime startTime, char phase, unsigned &currentIter)
278 280
         updateSampler(params, ASampler, PSampler, nA, nP);
279 281
         totalUpdates += nA + nP;
280 282
 
281
-        if (phase == 'S')
283
+        if (phase == GAPS_SAMPLING_PHASE)
282 284
         {
283 285
             stats.update(ASampler, PSampler);
284 286
             if (params.takePumpSamples)
285 287
             {
286 288
                 stats.updatePump(ASampler);
287 289
             }
290
+        }
291
+        if (params.snapshotPhase == phase || params.snapshotPhase == GAPS_ALL_PHASES)
292
+        {
288 293
             if (params.snapshotFrequency > 0 && ((currentIter + 1) % params.snapshotFrequency) == 0)
289 294
             {
290
-                stats.takeSnapshot();
295
+                stats.takeSnapshot(phase, ASampler, PSampler);
291 296
             }
292 297
         }
293 298
         displayStatus(params, ASampler, PSampler, startTime, phase,
... ...
@@ -405,7 +410,7 @@ const DataType &uncertainty, GapsRandomState *randState)
405 410
     // these variables will get overwritten by checkpoint if provided
406 411
     GapsStatistics stats(params.nGenes, params.nSamples, params.nPatterns);
407 412
     GapsRng rng(randState);
408
-    char phase = 'C';
413
+    GapsAlgorithmPhase phase(GAPS_EQUILIBRATION_PHASE);
409 414
     unsigned currentIter = 0;
410 415
     processCheckpoint(params, ASampler, PSampler, randState, stats, rng, phase, currentIter);
411 416
     calculateNumberOfThreads(params);
... ...
@@ -420,18 +425,18 @@ const DataType &uncertainty, GapsRandomState *randState)
420 425
     bpt::ptime startTime = bpt_now();
421 426
 
422 427
     // fallthrough through phases, allows algorithm to be resumed in either phase
423
-    GAPS_ASSERT(phase == 'C' || phase == 'S');
428
+    GAPS_ASSERT(phase == GAPS_EQUILIBRATION_PHASE || phase == GAPS_SAMPLING_PHASE);
424 429
     uint64_t totalUpdates = 0;
425 430
     switch (phase)
426 431
     {
427
-        case 'C':
432
+        case GAPS_EQUILIBRATION_PHASE:
428 433
             GAPS_MESSAGE(params.printMessages, "-- Equilibration Phase --\n");
429 434
             totalUpdates += runOnePhase(params, ASampler, PSampler, stats, randState,
430 435
                 rng, startTime, phase, currentIter);
431
-            phase = 'S';
436
+            phase = GAPS_SAMPLING_PHASE;
432 437
             currentIter = 0;
433 438
         // fall through
434
-        case 'S':
439
+        case GAPS_SAMPLING_PHASE:
435 440
             GAPS_MESSAGE(params.printMessages, "-- Sampling Phase --\n");
436 441
             totalUpdates += runOnePhase(params, ASampler, PSampler, stats, randState,
437 442
                 rng, startTime, phase, currentIter);
Browse code

Added option to take snapshots during sampling

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

Tom Sherman authored on 03/09/2019 15:23:59
Showing1 changed files
... ...
@@ -285,6 +285,10 @@ GapsRng &rng, bpt::ptime startTime, char phase, unsigned &currentIter)
285 285
             {
286 286
                 stats.updatePump(ASampler);
287 287
             }
288
+            if (params.snapshotFrequency > 0 && ((currentIter + 1) % params.snapshotFrequency) == 0)
289
+            {
290
+                stats.takeSnapshot();
291
+            }
288 292
         }
289 293
         displayStatus(params, ASampler, PSampler, startTime, phase,
290 294
             currentIter, stats);
Browse code

fixed bug where not all workers would print start/stop status during distributed runs

Tom Sherman authored on 28/06/2019 17:12:00
Showing1 changed files
... ...
@@ -392,7 +392,7 @@ const DataType &uncertainty, GapsRandomState *randState)
392 392
     }
393 393
 
394 394
     // if we are running distributed, each worker needs to print when it's started
395
-    if (params.runningDistributed && params.printMessages)
395
+    if (params.runningDistributed)
396 396
     {
397 397
         gaps_printf("    worker %d is starting!\n", params.workerID);
398 398
         gaps_flush();
... ...
@@ -449,7 +449,7 @@ const DataType &uncertainty, GapsRandomState *randState)
449 449
     }
450 450
 
451 451
     // if we are running distributed, each worker needs to print when it's done
452
-    if (params.runningDistributed && params.printMessages)
452
+    if (params.runningDistributed)
453 453
     {
454 454
         GapsTime elapsed(static_cast<unsigned>((bpt_now() - startTime).total_seconds()));
455 455
         gaps_printf("    worker %d is finished! Time: %02d:%02d:%02d\n",
sherman5 authored on 26/06/2019 18:58:07
Showing0 changed files
Browse code

fixed incorrect renaming of Equilibration phase

sherman5 authored on 26/06/2019 18:37:19
Showing1 changed files
... ...
@@ -229,7 +229,7 @@ GapsRng &rng, bpt::ptime startTime, char phase, unsigned &currentIter)
229 229
         createCheckpoint(params, ASampler, PSampler, randState, stats,
230 230
             rng, phase, currentIter);
231 231
         
232
-        // set annealing temperature in calibration phase
232
+        // set annealing temperature in Equilibration phase
233 233
         if (phase == 'C')
234 234
         {        
235 235
             float temp = static_cast<float>(2 * currentIter)
... ...
@@ -386,7 +386,7 @@ const DataType &uncertainty, GapsRandomState *randState)
386 386
     switch (phase)
387 387
     {
388 388
         case 'C':
389
-            GAPS_MESSAGE(params.printMessages, "-- Calibration Phase --\n");
389
+            GAPS_MESSAGE(params.printMessages, "-- Equilibration Phase --\n");
390 390
             runOnePhase(params, ASampler, PSampler, stats, randState, rng,
391 391
                 startTime, phase, currentIter);
392 392
             phase = 'S';
Browse code

fix compiler warnings

sherman5 authored on 25/06/2019 22:21:06
Showing1 changed files
... ...
@@ -426,7 +426,7 @@ const DataType &uncertainty, GapsRandomState *randState)
426 426
                 rng, startTime, phase, currentIter);
427 427
             phase = 'S';
428 428
             currentIter = 0;
429
-
429
+        // fall through
430 430
         case 'S':
431 431
             GAPS_MESSAGE(params.printMessages, "-- Sampling Phase --\n");
432 432
             totalUpdates += runOnePhase(params, ASampler, PSampler, stats, randState,
Browse code

clean up linter warnings

Tom Sherman authored on 24/06/2019 19:44:02
Showing1 changed files
... ...
@@ -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;
Browse code

passing debug tests locally

Tom Sherman authored on 24/06/2019 16:15:51
Showing1 changed files
... ...
@@ -39,7 +39,6 @@ struct GapsTime
39 39
     unsigned hours;
40 40
     unsigned minutes;
41 41
     unsigned seconds;
42
-
43 42
     GapsTime(unsigned totalSeconds)
44 43
     {
45 44
         hours = totalSeconds / 3600;
Browse code

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

Tom Sherman authored on 21/06/2019 19:44:45
Showing1 changed files
... ...
@@ -154,28 +154,28 @@ static void displayStatus(const GapsParameters &params,
154 154
 const Sampler &ASampler, const Sampler &PSampler, bpt::ptime startTime,
155 155
 char phase, unsigned iter, GapsStatistics &stats)
156 156
 {
157
-    if (params.printMessages && params.outputFrequency > 0
158
-    && ((iter + 1) % params.outputFrequency) == 0)
157
+    if (params.outputFrequency > 0 && ((iter + 1) % params.outputFrequency) == 0)
159 158
     {
160
-        bpt::time_duration diff = bpt_now() - startTime;
161
-        double perComplete = estimatedPercentComplete(params, ASampler,
162
-            PSampler, phase, iter);
163
-
164
-        GapsTime elapsedTime(static_cast<unsigned>(diff.total_seconds()));
165
-        GapsTime totalTime(static_cast<unsigned>(diff.total_seconds() / perComplete));
166
-
167 159
         float cs = PSampler.chiSq();
168 160
         unsigned nA = ASampler.nAtoms();
169 161
         unsigned nP = PSampler.nAtoms();
170
-
171 162
         stats.addChiSq(cs);
172 163
         stats.addAtomCount(nA, nP);
173
-
174
-        gaps_printf("%d of %d, Atoms: %d(A), %d(P), ChiSq: %.0f, Time: %02d:%02d:%02d / %02d:%02d:%02d\n",
175
-            iter + 1, params.nIterations, nA, nP, cs, elapsedTime.hours,
176
-            elapsedTime.minutes, elapsedTime.seconds, totalTime.hours,
177
-            totalTime.minutes, totalTime.seconds);
178
-        gaps_flush();
164
+        if (params.printMessages)
165
+        {
166
+            bpt::time_duration diff = bpt_now() - startTime;
167
+            double perComplete = estimatedPercentComplete(params, ASampler,
168
+                PSampler, phase, iter);
169
+
170
+            GapsTime elapsedTime(static_cast<unsigned>(diff.total_seconds()));
171
+            GapsTime totalTime(static_cast<unsigned>(diff.total_seconds() / perComplete));
172
+
173
+            gaps_printf("%d of %d, Atoms: %d(A), %d(P), ChiSq: %.0f, Time: %02d:%02d:%02d / %02d:%02d:%02d\n",
174
+                iter + 1, params.nIterations, nA, nP, cs, elapsedTime.hours,
175
+                elapsedTime.minutes, elapsedTime.seconds, totalTime.hours,
176
+                totalTime.minutes, totalTime.seconds);
177
+            gaps_flush();
178
+        }
179 179
     }
180 180
 }
181 181
 
Browse code

check if printing is requested before displaying messages about sampler type

sherman5 authored on 20/06/2019 20:47:54
Showing1 changed files
... ...
@@ -63,11 +63,11 @@ const DataType &uncertainty, GapsRandomState *randState)
63 63
 {
64 64
     if (params.asynchronousUpdates)
65 65
     {
66
-        gaps_printf("Sampler Type: Asynchronous\n");
66
+        GAPS_MESSAGE(params.printMessages, "Sampler Type: Asynchronous\n");
67 67
         return runCoGAPSAlgorithm< AsynchronousGibbsSampler<DataModel> >(data,
68 68
             params, uncertainty, randState);
69 69
     }
70
-    gaps_printf("Sampler Type: Sequential\n");
70
+    GAPS_MESSAGE(params.printMessages, "Sampler Type: Sequential\n");
71 71
     return runCoGAPSAlgorithm< SingleThreadedGibbsSampler<DataModel> >(data,
72 72
         params, uncertainty, randState);
73 73
 }
... ...
@@ -78,10 +78,10 @@ const DataType &uncertainty, GapsRandomState *randState)
78 78
 {
79 79
     if (params.useSparseOptimization)
80 80
     {
81
-        gaps_printf("Data Model: Sparse, Normal\n");
81
+        GAPS_MESSAGE(params.printMessages, "Data Model: Sparse, Normal\n");
82 82
         return chooseSampler<SparseNormalModel>(data, params, uncertainty, randState);
83 83
     }
84
-    gaps_printf("Data Model: Dense, Normal\n");        
84
+    GAPS_MESSAGE(params.printMessages, "Data Model: Dense, Normal\n");        
85 85
     return chooseSampler<DenseNormalModel>(data, params, uncertainty, randState);
86 86
 }
87 87
 
Browse code

fixed cleanup script and removed incorrect assert statements

sherman5 authored on 19/06/2019 18:28:25
Showing1 changed files
... ...
@@ -296,9 +296,9 @@ Sampler &PSampler)
296 296
     // check if we're fixing a matrix
297 297
     if (params.useFixedPatterns)
298 298
     {
299
+        GAPS_ASSERT(params.fixedPatterns.nCol() == params.nPatterns);
299 300
         switch (params.whichMatrixFixed)
300 301
         {
301
-            GAPS_ASSERT(params.fixedPatterns.nCol() == params.nPatterns);
302 302
             case 'A' :
303 303
                 GAPS_ASSERT(params.fixedPatterns.nRow() == params.nGenes);
304 304
                 ASampler.setMatrix(params.fixedPatterns);
Browse code

return total running time in result object

sherman5 authored on 19/06/2019 03:29:24
Showing1 changed files
... ...
@@ -400,6 +400,7 @@ const DataType &uncertainty, GapsRandomState *randState)
400 400
     
401 401
     // get result
402 402
     GapsResult result(stats);
403
+    result.totalRunningTime = static_cast<unsigned>((bpt_now() - startTime).total_seconds());
403 404
     result.meanChiSq = stats.meanChiSq(PSampler);
404 405
     result.averageQueueLengthA = ASampler.getAverageQueueLength();
405 406
     result.averageQueueLengthP = PSampler.getAverageQueueLength();
Browse code

working implementation of concurrent and single-threaded atomic domain

sherman5 authored on 19/06/2019 03:09:22
Showing1 changed files
... ...
@@ -7,7 +7,10 @@
7 7
 #include "gibbs_sampler/SparseNormalModel.h"
8 8
 
9 9
 #ifdef __GAPS_R_BUILD__
10
+#pragma GCC diagnostic push
11
+#pragma GCC diagnostic ignored "-Wdeprecated-declarations"
10 12
 #include <Rcpp.h>
13
+#pragma GCC diagnostic pop
11 14
 #endif
12 15
 
13 16
 #ifdef __GAPS_OPENMP__
... ...
@@ -15,7 +18,10 @@
15 18
 #endif
16 19
 
17 20
 // boost time helpers
21
+#pragma GCC diagnostic push
22
+#pragma GCC diagnostic ignored "-Wdeprecated-declarations"
18 23
 #include <boost/date_time/posix_time/posix_time.hpp>
24
+#pragma GCC diagnostic pop
19 25
 namespace bpt = boost::posix_time;
20 26
 #define bpt_now() bpt::microsec_clock::local_time()
21 27
 
Browse code

record number of updates in the result object

Tom Sherman authored on 18/06/2019 18:15:06
Showing1 changed files
... ...
@@ -243,10 +243,11 @@ GapsRng &rng, char &phase, unsigned &currentIter)
243 243
 }
244 244
 
245 245
 template <class Sampler>
246
-static void runOnePhase(const GapsParameters &params, Sampler &ASampler,
246
+static uint64_t runOnePhase(const GapsParameters &params, Sampler &ASampler,
247 247
 Sampler &PSampler, GapsStatistics &stats, const GapsRandomState *randState,
248 248
 GapsRng &rng, bpt::ptime startTime, char phase, unsigned &currentIter)
249 249
 {
250
+    uint64_t totalUpdates = 0;
250 251
     for (; currentIter < params.nIterations; ++currentIter)
251 252
     {
252 253
         gaps_check_interrupt();
... ...
@@ -266,6 +267,7 @@ GapsRng &rng, bpt::ptime startTime, char phase, unsigned &currentIter)
266 267
         unsigned nA = rng.poisson(gaps::max(ASampler.nAtoms(), 10));
267 268
         unsigned nP = rng.poisson(gaps::max(PSampler.nAtoms(), 10));
268 269
         updateSampler(params, ASampler, PSampler, nA, nP);
270
+        totalUpdates += nA + nP;
269 271
 
270 272
         if (phase == 'S')
271 273
         {
... ...
@@ -278,6 +280,7 @@ GapsRng &rng, bpt::ptime startTime, char phase, unsigned &currentIter)
278 280
         displayStatus(params, ASampler, PSampler, startTime, phase,
279 281
             currentIter, stats);
280 282
     }
283
+    return totalUpdates;
281 284
 }
282 285
 
283 286
 template <class Sampler>
... ...
@@ -405,19 +408,20 @@ const DataType &uncertainty, GapsRandomState *randState)
405 408
 
406 409
     // fallthrough through phases, allows algorithm to be resumed in either phase
407 410
     GAPS_ASSERT(phase == 'C' || phase == 'S');
411
+    uint64_t totalUpdates = 0;
408 412
     switch (phase)
409 413
     {
410 414
         case 'C':
411 415
             GAPS_MESSAGE(params.printMessages, "-- Calibration Phase --\n");
412
-            runOnePhase(params, ASampler, PSampler, stats, randState, rng,
413
-                startTime, phase, currentIter);
416
+            totalUpdates += runOnePhase(params, ASampler, PSampler, stats, randState,
417
+                rng, startTime, phase, currentIter);
414 418
             phase = 'S';
415 419
             currentIter = 0;
416 420
 
417 421
         case 'S':
418 422
             GAPS_MESSAGE(params.printMessages, "-- Sampling Phase --\n");
419
-            runOnePhase(params, ASampler, PSampler, stats, randState, rng,
420
-                startTime, phase, currentIter);
423
+            totalUpdates += runOnePhase(params, ASampler, PSampler, stats, randState,
424
+                rng, startTime, phase, currentIter);
421 425
     }
422 426
     
423 427
     // get result
... ...
@@ -426,6 +430,7 @@ const DataType &uncertainty, GapsRandomState *randState)
426 430
     result.meanChiSq = stats.meanChiSq(PSampler);
427 431
     result.averageQueueLengthA = ASampler.getAverageQueueLength();
428 432
     result.averageQueueLengthP = PSampler.getAverageQueueLength();
433
+    result.totalUpdates = totalUpdates;
429 434
 
430 435
     // handle pump statistics
431 436
     if (params.takePumpSamples)
... ...
@@ -442,7 +447,6 @@ const DataType &uncertainty, GapsRandomState *randState)
442 447
             params.workerID, elapsed.hours, elapsed.minutes, elapsed.seconds);
443 448
         gaps_flush();
444 449
     }
445
-
446 450
     return result;
447 451
 }
448 452
 
Tom Sherman authored on 18/06/2019 18:02:06
Showing0 changed files
Browse code

return running time in result object's metadata

Tom Sherman authored on 18/06/2019 18:00:21
Showing1 changed files
... ...
@@ -422,6 +422,7 @@ const DataType &uncertainty, GapsRandomState *randState)
422 422
     
423 423
     // get result
424 424
     GapsResult result(stats);
425
+    result.totalRunningTime = static_cast<unsigned>((bpt_now() - startTime).total_seconds());
425 426
     result.meanChiSq = stats.meanChiSq(PSampler);
426 427
     result.averageQueueLengthA = ASampler.getAverageQueueLength();
427 428
     result.averageQueueLengthP = PSampler.getAverageQueueLength();
... ...
@@ -436,8 +437,7 @@ const DataType &uncertainty, GapsRandomState *randState)
436 437
     // if we are running distributed, each worker needs to print when it's done
437 438
     if (params.runningDistributed)
438 439
     {
439
-        bpt::time_duration diff = bpt_now() - startTime;
440
-        GapsTime elapsed(static_cast<unsigned>(diff.total_seconds()));
440
+        GapsTime elapsed(static_cast<unsigned>((bpt_now() - startTime).total_seconds()));
441 441
         gaps_printf("    worker %d is finished! Time: %02d:%02d:%02d\n",
442 442
             params.workerID, elapsed.hours, elapsed.minutes, elapsed.seconds);
443 443
         gaps_flush();
Browse code

added a clearer display message for the number of atoms in each domain

Tom Sherman authored on 04/06/2019 14:12:13
Showing1 changed files
... ...
@@ -145,7 +145,7 @@ char phase, unsigned iter, GapsStatistics &stats)
145 145
         stats.addChiSq(cs);
146 146
         stats.addAtomCount(nA, nP);
147 147
 
148
-        gaps_printf("%d of %d, Atoms: %d(%d), ChiSq: %.0f, Time: %02d:%02d:%02d / %02d:%02d:%02d\n",
148
+        gaps_printf("%d of %d, Atoms: %d(A), %d(P), ChiSq: %.0f, Time: %02d:%02d:%02d / %02d:%02d:%02d\n",
149 149
             iter + 1, params.nIterations, nA, nP, cs, elapsedTime.hours,
150 150
             elapsedTime.minutes, elapsedTime.seconds, totalTime.hours,
151 151
             totalTime.minutes, totalTime.seconds);
Browse code

all sampler/model combinations seem to be working

Tom Sherman authored on 03/06/2019 17:59:25
Showing1 changed files
... ...
@@ -4,7 +4,6 @@
4 4
 #include "gibbs_sampler/AsynchronousGibbsSampler.h"
5 5
 #include "gibbs_sampler/SingleThreadedGibbsSampler.h"
6 6
 #include "gibbs_sampler/DenseNormalModel.h"
7
-#include "gibbs_sampler/DenseNormalWithUncertaintyModel.h"
8 7
 #include "gibbs_sampler/SparseNormalModel.h"
9 8
 
10 9
 #ifdef __GAPS_R_BUILD__
... ...
@@ -58,9 +57,11 @@ const DataType &uncertainty, GapsRandomState *randState)
58 57
 {
59 58
     if (params.asynchronousUpdates)
60 59
     {
60
+        gaps_printf("Sampler Type: Asynchronous\n");
61 61
         return runCoGAPSAlgorithm< AsynchronousGibbsSampler<DataModel> >(data,
62 62
             params, uncertainty, randState);
63 63
     }
64
+    gaps_printf("Sampler Type: Sequential\n");
64 65
     return runCoGAPSAlgorithm< SingleThreadedGibbsSampler<DataModel> >(data,
65 66
         params, uncertainty, randState);
66 67
 }
... ...
@@ -71,13 +72,11 @@ const DataType &uncertainty, GapsRandomState *randState)
71 72
 {
72 73
     if (params.useSparseOptimization)
73 74
     {
75
+        gaps_printf("Data Model: Sparse, Normal\n");
74 76
         return chooseSampler<SparseNormalModel>(data, params, uncertainty, randState);
75 77
     }
76
-    if (uncertainty.empty())
77
-    {
78
-        return chooseSampler<DenseNormalModel>(data, params, uncertainty, randState);
79
-    }
80
-    return chooseSampler<DenseNormalWithUncertaintyModel>(data, params, uncertainty, randState);
78
+    gaps_printf("Data Model: Dense, Normal\n");        
79
+    return chooseSampler<DenseNormalModel>(data, params, uncertainty, randState);
81 80
 }
82 81
 
83 82
 // helper function, this dispatches the correct run function depending
Browse code

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

Tom Sherman authored on 23/05/2019 18:55:10
Showing1 changed files
... ...
@@ -4,6 +4,7 @@
4 4
 #include "gibbs_sampler/AsynchronousGibbsSampler.h"
5 5
 #include "gibbs_sampler/SingleThreadedGibbsSampler.h"
6 6
 #include "gibbs_sampler/DenseNormalModel.h"
7
+#include "gibbs_sampler/DenseNormalWithUncertaintyModel.h"
7 8
 #include "gibbs_sampler/SparseNormalModel.h"
8 9
 
9 10
 #ifdef __GAPS_R_BUILD__
... ...
@@ -72,7 +73,11 @@ const DataType &uncertainty, GapsRandomState *randState)
72 73
     {
73 74
         return chooseSampler<SparseNormalModel>(data, params, uncertainty, randState);
74 75
     }
75
-    return chooseSampler<DenseNormalModel>(data, params, uncertainty, randState);
76
+    if (uncertainty.empty())
77
+    {
78
+        return chooseSampler<DenseNormalModel>(data, params, uncertainty, randState);
79
+    }
80
+    return chooseSampler<DenseNormalWithUncertaintyModel>(data, params, uncertainty, randState);
76 81
 }
77 82
 
78 83
 // helper function, this dispatches the correct run function depending
Browse code

format and comments

Tom Sherman authored on 23/05/2019 18:00:31
Showing1 changed files
... ...
@@ -27,6 +27,7 @@ namespace bpt = boost::posix_time;
27 27
         } \
28 28
     } while(0)
29 29
 
30
+// for converting seconds to h:m:s
30 31
 struct GapsTime
31 32
 {
32 33
     unsigned hours;
... ...
@@ -75,7 +76,7 @@ const DataType &uncertainty, GapsRandomState *randState)
75 76
 }
76 77
 
77 78
 // helper function, this dispatches the correct run function depending
78
-// on the type of GibbsSampler needed for the given parameters
79
+// on the type of GibbsSampler and DataModel needed for the given parameters
79 80
 template <class DataType>
80 81
 static GapsResult run_helper(const DataType &data, GapsParameters &params,
81 82
 const DataType &uncertainty, GapsRandomState *randState)
... ...
@@ -131,13 +132,10 @@ const Sampler &ASampler, const Sampler &PSampler, char phase, unsigned iter)
131 132
     }
132 133
 
133 134
     double totalIter = 2.0 * static_cast<double>(params.nIterations);
134
-
135 135
     double estimatedCompleted = estimatedNumUpdates(nIter, nIter, nAtomsA) + 
136 136
         estimatedNumUpdates(nIter, nIter, nAtomsP);
137
-
138 137
     double estimatedTotal = estimatedNumUpdates(nIter, totalIter, nAtomsA) + 
139 138
         estimatedNumUpdates(nIter, totalIter, nAtomsP);
140
-
141 139
     return estimatedCompleted / estimatedTotal;
142 140
 }
143 141
 
... ...
@@ -199,8 +197,7 @@ static void createCheckpoint(const GapsParameters &params,
199 197
 Sampler &ASampler, Sampler &PSampler, const GapsRandomState *randState,
200 198
 const GapsStatistics &stats, const GapsRng &rng, char phase, unsigned iter)
201 199
 {
202
-    if (params.checkpointInterval > 0
203
-    && ((iter + 1) % params.checkpointInterval) == 0
200
+    if (params.checkpointInterval > 0 && ((iter + 1) % params.checkpointInterval) == 0
204 201
     && !params.subsetData)
205 202
     {
206 203
         // create backup file
... ...
@@ -216,6 +213,11 @@ const GapsStatistics &stats, const GapsRng &rng, char phase, unsigned iter)
216 213
         // delete backup file
217 214
         std::remove((params.checkpointOutFile + ".backup").c_str());
218 215
 
216
+        // running the extra initialization here allows for consistency with runs
217
+        // started from a checkpoint. This initialization phase will be run first 
218
+        // thing once a checkpoint is loaded since large matrices which aren't stored
219
+        // need to be initialized. By running it here we make sure that the algorithm
220
+        // is in the same state it will be when started from a checkpoint
219 221
         ASampler.extraInitialization();
220 222
         PSampler.extraInitialization();
221 223
     }
... ...
@@ -332,7 +334,7 @@ const DataType &uncertainty, GapsRandomState *randState)
332 334
 {
333 335
     // check if running in debug mode
334 336
     #ifdef GAPS_DEBUG
335
-    GAPS_MESSAGE(params.printMessages, "Running in debug mode\n");
337
+    gaps_printf("Running in debug mode\n");
336 338
     #endif
337 339
 
338 340
     // load data into gibbs samplers
... ...
@@ -346,7 +348,6 @@ const DataType &uncertainty, GapsRandomState *randState)
346 348
     // the process or waiting for it to finish
347 349
     GAPS_MESSAGE(params.printMessages, "Loading Data...");
348 350
     bpt::ptime readStart = bpt_now();
349
-
350 351
     gaps_check_interrupt();
351 352
     Sampler ASampler(data, !params.transposeData, !params.subsetGenes,
352 353
         params.alphaA, params.maxGibbsMassA, params, randState);
... ...
@@ -359,12 +360,6 @@ const DataType &uncertainty, GapsRandomState *randState)
359 360
     processFixedMatrix(params, ASampler, PSampler);
360 361
     gaps_check_interrupt();
361 362
 
362
-    // check if data is sparse and sparseOptimization is not enabled
363
-    if (params.printMessages && !params.useSparseOptimization && ASampler.dataSparsity() > 0.80f)
364
-    {
365
-        gaps_printf("\nWarning: data is more than 80%% sparse and sparseOptimization is not enabled\n");
366
-    }
367
-
368 363
     // elapsed time for reading data
369 364
     bpt::time_duration readDiff = bpt_now() - readStart;
370 365
     GapsTime elapsed(static_cast<unsigned>(readDiff.total_seconds()));
... ...
@@ -374,6 +369,12 @@ const DataType &uncertainty, GapsRandomState *randState)
374 369
             elapsed.seconds);
375 370
     }
376 371
 
372
+    // check if data is sparse and sparseOptimization is not enabled
373
+    if (params.printMessages && !params.useSparseOptimization && ASampler.dataSparsity() > 0.80f)
374
+    {
375
+        gaps_printf("\nWarning: data is more than 80%% sparse and sparseOptimization is not enabled\n");
376
+    }
377
+
377 378
     // if we are running distributed, each worker needs to print when it's started
378 379
     if (params.runningDistributed)
379 380
     {
... ...
@@ -386,8 +387,7 @@ const DataType &uncertainty, GapsRandomState *randState)
386 387
     GapsRng rng(randState);
387 388
     char phase = 'C';
388 389
     unsigned currentIter = 0;
389
-    processCheckpoint(params, ASampler, PSampler, randState, stats, rng,
390
-        phase, currentIter);
390
+    processCheckpoint(params, ASampler, PSampler, randState, stats, rng, phase, currentIter);
391 391
     calculateNumberOfThreads(params);
392 392
 
393 393
     // sync samplers and run any additional initialization needed
... ...
@@ -422,6 +422,7 @@ const DataType &uncertainty, GapsRandomState *randState)
422 422
     result.averageQueueLengthA = ASampler.getAverageQueueLength();
423 423
     result.averageQueueLengthP = PSampler.getAverageQueueLength();
424 424
 
425
+    // handle pump statistics
425 426
     if (params.takePumpSamples)
426 427
     {
427 428
         result.pumpMatrix = stats.pumpMatrix();
Browse code

infrastructure in place for extending samplers and models

Tom Sherman authored on 22/05/2019 21:27:14
Showing1 changed files
... ...
@@ -1,10 +1,10 @@
1 1
 #include "GapsRunner.h"
2 2
 
3 3
 #include "utils/Archive.h"
4
-#include "gibbs_sampler/GibbsSampler.h"
4
+#include "gibbs_sampler/AsynchronousGibbsSampler.h"
5 5
 #include "gibbs_sampler/SingleThreadedGibbsSampler.h"
6
-#include "gibbs_sampler/DenseStoragePolicy.h"
7
-#include "gibbs_sampler/SparseStoragePolicy.h"
6
+#include "gibbs_sampler/DenseNormalModel.h"
7
+#include "gibbs_sampler/SparseNormalModel.h"
8 8
 
9 9
 #ifdef __GAPS_R_BUILD__
10 10
 #include <Rcpp.h>
... ...
@@ -50,28 +50,28 @@ static GapsResult runCoGAPSAlgorithm(const DataType &data, GapsParameters &param
50 50
 
51 51
 ////////////////////////////////////////////////////////////////////////////////
52 52
 
53
-template <class StorageType, class DataType>
53
+template <class DataModel, class DataType>
54 54
 static GapsResult chooseSampler(const DataType &data, GapsParameters &params,
55 55
 const DataType &uncertainty, GapsRandomState *randState)
56 56
 {
57 57
     if (params.asynchronousUpdates)
58 58
     {
59
-        return runCoGAPSAlgorithm< GibbsSampler<StorageType> >(data,
59
+        return runCoGAPSAlgorithm< AsynchronousGibbsSampler<DataModel> >(data,
60 60
             params, uncertainty, randState);
61 61
     }
62
-    return runCoGAPSAlgorithm< SingleThreadedGibbsSampler<StorageType> >(data,
62
+    return runCoGAPSAlgorithm< SingleThreadedGibbsSampler<DataModel> >(data,
63 63
         params, uncertainty, randState);
64 64
 }
65 65
 
66 66
 template <class DataType>
67
-static GapsResult chooseStorage(const DataType &data, GapsParameters &params,
67
+static GapsResult chooseDataModel(const DataType &data, GapsParameters &params,
68 68
 const DataType &uncertainty, GapsRandomState *randState)
69 69
 {
70 70
     if (params.useSparseOptimization)
71 71
     {
72
-        return chooseSampler<SparseStorage>(data, params, uncertainty, randState);
72
+        return chooseSampler<SparseNormalModel>(data, params, uncertainty, randState);
73 73
     }
74
-    return chooseSampler<DenseStorage>(data, params, uncertainty, randState);
74
+    return chooseSampler<DenseNormalModel>(data, params, uncertainty, randState);
75 75
 }
76 76
 
77 77
 // helper function, this dispatches the correct run function depending
... ...
@@ -87,7 +87,7 @@ const DataType &uncertainty, GapsRandomState *randState)
87 87
         ar >> params;
88 88
         ar >> *randState;
89 89
     }
90
-    return chooseStorage(data, params, uncertainty, randState);
90
+    return chooseDataModel(data, params, uncertainty, randState);
91 91
 }
92 92
 
93 93
 // these two functions are the top-level functions exposed to the C++
Browse code

single-threaded sampler working

Tom Sherman authored on 22/05/2019 17:38:53
Showing1 changed files
... ...
@@ -2,6 +2,7 @@
2 2
 
3 3
 #include "utils/Archive.h"
4 4
 #include "gibbs_sampler/GibbsSampler.h"
5
+#include "gibbs_sampler/SingleThreadedGibbsSampler.h"
5 6
 #include "gibbs_sampler/DenseStoragePolicy.h"
6 7
 #include "gibbs_sampler/SparseStoragePolicy.h"
7 8
 
... ...
@@ -49,6 +50,30 @@ static GapsResult runCoGAPSAlgorithm(const DataType &data, GapsParameters &param
49 50
 
50 51
 ////////////////////////////////////////////////////////////////////////////////
51 52
 
53
+template <class StorageType, class DataType>
54
+static GapsResult chooseSampler(const DataType &data, GapsParameters &params,
55
+const DataType &uncertainty, GapsRandomState *randState)
56
+{
57
+    if (params.asynchronousUpdates)
58
+    {
59
+        return runCoGAPSAlgorithm< GibbsSampler<StorageType> >(data,
60
+            params, uncertainty, randState);
61
+    }
62
+    return runCoGAPSAlgorithm< SingleThreadedGibbsSampler<StorageType> >(data,
63
+        params, uncertainty, randState);
64
+}
65
+
66
+template <class DataType>
67
+static GapsResult chooseStorage(const DataType &data, GapsParameters &params,
68
+const DataType &uncertainty, GapsRandomState *randState)
69
+{
70
+    if (params.useSparseOptimization)
71
+    {
72
+        return chooseSampler<SparseStorage>(data, params, uncertainty, randState);
73
+    }
74
+    return chooseSampler<DenseStorage>(data, params, uncertainty, randState);
75
+}
76
+
52 77
 // helper function, this dispatches the correct run function depending
53 78
 // on the type of GibbsSampler needed for the given parameters
54 79
 template <class DataType>
... ...
@@ -62,14 +87,7 @@ const DataType &uncertainty, GapsRandomState *randState)
62 87
         ar >> params;
63 88
         ar >> *randState;
64 89
     }
65
-
66
-    if (params.useSparseOptimization)
67
-    {
68
-        return runCoGAPSAlgorithm< GibbsSampler<SparseStorage> >(data, params,
69
-            uncertainty, randState);
70
-    }
71
-    return runCoGAPSAlgorithm< GibbsSampler<DenseStorage> >(data, params,
72
-        uncertainty, randState);
90
+    return chooseStorage(data, params, uncertainty, randState);
73 91
 }
74 92
 
75 93
 // these two functions are the top-level functions exposed to the C++
Browse code

improve unit tests

Tom Sherman authored on 04/04/2019 21:09:16
Showing1 changed files
... ...
@@ -357,7 +357,7 @@ const DataType &uncertainty, GapsRandomState *randState)
357 357
     }
358 358
 
359 359
     // if we are running distributed, each worker needs to print when it's started
360
-    if (params.runningDistributed)
360
+    if (params.runningDistributed && params.printMessages)
361 361
     {
362 362
         gaps_printf("    worker %d is starting!\n", params.workerID);
363 363
         gaps_flush();
... ...
@@ -411,7 +411,7 @@ const DataType &uncertainty, GapsRandomState *randState)
411 411
     }
412 412
 
413 413
     // if we are running distributed, each worker needs to print when it's done
414
-    if (params.runningDistributed)
414
+    if (params.runningDistributed && params.printMessages)
415 415
     {
416 416
         bpt::time_duration diff = bpt_now() - startTime;
417 417
         GapsTime elapsed(static_cast<unsigned>(diff.total_seconds()));
Browse code

added warning for when data is sparse and sparseOptimization is not enabled

Tom Sherman authored on 28/03/2019 17:38:34
Showing1 changed files
... ...
@@ -341,6 +341,12 @@ const DataType &uncertainty, GapsRandomState *randState)
341 341
     processFixedMatrix(params, ASampler, PSampler);
342 342
     gaps_check_interrupt();
343 343
 
344
+    // check if data is sparse and sparseOptimization is not enabled
345
+    if (params.printMessages && !params.useSparseOptimization && ASampler.dataSparsity() > 0.80f)
346
+    {
347
+        gaps_printf("\nWarning: data is more than 80%% sparse and sparseOptimization is not enabled\n");
348
+    }
349
+
344 350
     // elapsed time for reading data
345 351
     bpt::time_duration readDiff = bpt_now() - readStart;
346 352
     GapsTime elapsed(static_cast<unsigned>(readDiff.total_seconds()));
Browse code

print start message for workers

Tom Sherman authored on 19/02/2019 03:18:22
Showing1 changed files
... ...
@@ -350,6 +350,13 @@ const DataType &uncertainty, GapsRandomState *randState)
350 350
             elapsed.seconds);
351 351
     }
352 352
 
353
+    // if we are running distributed, each worker needs to print when it's started
354
+    if (params.runningDistributed)
355
+    {
356
+        gaps_printf("    worker %d is starting!\n", params.workerID);
357
+        gaps_flush();
358
+    }
359
+
353 360
     // these variables will get overwritten by checkpoint if provided
354 361
     GapsStatistics stats(params.nGenes, params.nSamples, params.nPatterns);
355 362
     GapsRng rng(randState);
Tom Sherman authored on 07/02/2019 21:34:56
Showing0 changed files
Browse code

check for interrupt more often when reading data

Tom Sherman authored on 22/01/2019 23:29:45
Showing1 changed files
... ...
@@ -222,10 +222,7 @@ GapsRng &rng, bpt::ptime startTime, char phase, unsigned &currentIter)
222 222
 {
223 223
     for (; currentIter < params.nIterations; ++currentIter)
224 224
     {
225
-        #ifdef __GAPS_R_BUILD__
226
-        Rcpp::checkUserInterrupt();
227
-        #endif
228
-
225
+        gaps_check_interrupt();
229 226
         createCheckpoint(params, ASampler, PSampler, randState, stats,
230 227
             rng, phase, currentIter);
231 228
         
... ...
@@ -323,12 +320,17 @@ const DataType &uncertainty, GapsRandomState *randState)
323 320
     // within the sampler, note the subsetting genes/samples flag must be
324 321
     // flipped if we are flipping the transpose flag
325 322
     GAPS_MESSAGE(params.printMessages, "Loading Data...");
323
+    gaps_check_interrupt();
326 324
     Sampler ASampler(data, !params.transposeData, !params.subsetGenes,
327 325
         params.alphaA, params.maxGibbsMassA, params, randState);
326
+    gaps_check_interrupt();
328 327
     Sampler PSampler(data, params.transposeData, params.subsetGenes,
329 328
         params.alphaP, params.maxGibbsMassP, params, randState);
329
+    gaps_check_interrupt();
330 330
     processUncertainty(params, ASampler, PSampler, uncertainty);
331
+    gaps_check_interrupt();
331 332
     processFixedMatrix(params, ASampler, PSampler);
333
+    gaps_check_interrupt();
332 334
     GAPS_MESSAGE(params.printMessages, "Done!\n");
333 335
 
334 336
     // these variables will get overwritten by checkpoint if provided
Browse code

return average queue length in metadata

Tom Sherman authored on 09/01/2019 22:01:27
Showing1 changed files
... ...
@@ -369,6 +369,8 @@ const DataType &uncertainty, GapsRandomState *randState)
369 369
     // get result
370 370
     GapsResult result(stats);
371 371
     result.meanChiSq = stats.meanChiSq(PSampler);
372
+    result.averageQueueLengthA = ASampler.getAverageQueueLength();
373
+    result.averageQueueLengthP = PSampler.getAverageQueueLength();
372 374
 
373 375
     if (params.takePumpSamples)
374 376
     {
Browse code

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

Tom Sherman authored on 08/01/2019 22:16:52
Showing1 changed files
... ...
@@ -1,6 +1,9 @@
1 1
 #include "GapsRunner.h"
2 2
 
3 3
 #include "utils/Archive.h"
4
+#include "gibbs_sampler/GibbsSampler.h"
5
+#include "gibbs_sampler/DenseStoragePolicy.h"
6
+#include "gibbs_sampler/SparseStoragePolicy.h"
4 7
 
5 8
 #ifdef __GAPS_R_BUILD__
6 9
 #include <Rcpp.h>
... ...
@@ -62,10 +65,10 @@ const DataType &uncertainty, GapsRandomState *randState)
62 65
 
63 66
     if (params.useSparseOptimization)
64 67
     {
65
-        return runCoGAPSAlgorithm<SparseGibbsSampler>(data, params,
68
+        return runCoGAPSAlgorithm< GibbsSampler<SparseStorage> >(data, params,
66 69
             uncertainty, randState);
67 70
     }
68
-    return runCoGAPSAlgorithm<DenseGibbsSampler>(data, params,
71
+    return runCoGAPSAlgorithm< GibbsSampler<DenseStorage> >(data, params,
69 72
         uncertainty, randState);
70 73
 }
71 74
 
Browse code

link against Rhdf5 library

Tom Sherman authored on 08/01/2019 15:10:46
Showing1 changed files
... ...
@@ -323,13 +323,20 @@ const DataType &uncertainty, GapsRandomState *randState)
323 323
     // within the sampler, note the subsetting genes/samples flag must be
324 324
     // flipped if we are flipping the transpose flag
325 325
     GAPS_MESSAGE(params.printMessages, "Loading Data...");
326
+    bpt::ptime readStart = bpt_now();
326 327
     Sampler ASampler(data, !params.transposeData, !params.subsetGenes,
327 328
         params.alphaA, params.maxGibbsMassA, params, randState);
328 329
     Sampler PSampler(data, params.transposeData, params.subsetGenes,
329 330
         params.alphaP, params.maxGibbsMassP, params, randState);
330 331
     processUncertainty(params, ASampler, PSampler, uncertainty);
331 332
     processFixedMatrix(params, ASampler, PSampler);
332
-    GAPS_MESSAGE(params.printMessages, "Done!\n");
333
+    bpt::time_duration readDiff = bpt_now() - readStart;
334
+    GapsTime elapsed(static_cast<unsigned>(readDiff.total_seconds()));
335
+    if (params.printMessages)
336
+    {
337
+        gaps_printf("Done! (%02d:%02d:%02d)\n", elapsed.hours, elapsed.minutes,
338
+            elapsed.seconds);
339
+    }
333 340
 
334 341
     // these variables will get overwritten by checkpoint if provided
335 342
     GapsStatistics stats(params.nGenes, params.nSamples, params.nPatterns);
Browse code

added PUMP back to COGAPS

Tom Sherman authored on 21/12/2018 00:52:37
Showing1 changed files
... ...
@@ -246,6 +246,10 @@ GapsRng &rng, bpt::ptime startTime, char phase, unsigned &currentIter)
246 246
         if (phase == 'S')
247 247
         {
248 248
             stats.update(ASampler, PSampler);
249
+            if (params.takePumpSamples)
250
+            {
251
+                stats.updatePump(ASampler);
252
+            }
249 253
         }
250 254
         displayStatus(params, ASampler, PSampler, startTime, phase,
251 255
             currentIter, stats);
... ...
@@ -366,6 +370,12 @@ const DataType &uncertainty, GapsRandomState *randState)
366 370
     GapsResult result(stats);
367 371
     result.meanChiSq = stats.meanChiSq(PSampler);
368 372
 
373
+    if (params.takePumpSamples)
374
+    {
375
+        result.pumpMatrix = stats.pumpMatrix();
376
+        result.meanPatternAssignment = stats.meanPattern();
377
+    }
378
+
369 379
     // if we are running distributed, each worker needs to print when it's done
370 380
     if (params.runningDistributed)
371 381
     {
Browse code

added fixedPatterns option for normal CoGAPS

Tom Sherman authored on 19/12/2018 14:59:41
Showing1 changed files
... ...
@@ -257,18 +257,18 @@ static void processFixedMatrix(const GapsParameters &params, Sampler &ASampler,
257 257
 Sampler &PSampler)
258 258
 {
259 259
     // check if we're fixing a matrix
260
-    if (params.useFixedMatrix)
260
+    if (params.useFixedPatterns)
261 261
     {
262
-        switch (params.whichFixedMatrix)
262
+        switch (params.whichMatrixFixed)
263 263
         {
264
-            GAPS_ASSERT(params.fixedMatrix.nCol() == params.nPatterns);
264
+            GAPS_ASSERT(params.fixedPatterns.nCol() == params.nPatterns);
265 265
             case 'A' :
266
-                GAPS_ASSERT(params.fixedMatrix.nRow() == params.nGenes);
267
-                ASampler.setMatrix(params.fixedMatrix);
266
+                GAPS_ASSERT(params.fixedPatterns.nRow() == params.nGenes);
267
+                ASampler.setMatrix(params.fixedPatterns);
268 268
                 break;
269 269
             case 'P' :
270
-                GAPS_ASSERT(params.fixedMatrix.nRow() == params.nSamples);
271
-                PSampler.setMatrix(params.fixedMatrix);
270
+                GAPS_ASSERT(params.fixedPatterns.nRow() == params.nSamples);
271
+                PSampler.setMatrix(params.fixedPatterns);
272 272
                 break;
273 273
             default: break; // 'N' for none
274 274
         }
Browse code

get started on fixed matrix option

Tom Sherman authored on 18/12/2018 18:00:13
Showing1 changed files
... ...
@@ -154,19 +154,19 @@ template <class Sampler>
154 154
 static void updateSampler(const GapsParameters &params, Sampler &ASampler,
155 155
 Sampler &PSampler, unsigned nA, unsigned nP)
156 156
 {
157
-    if (!params.useFixedMatrix || params.whichFixedMatrix != 'A')
157
+    if (params.whichMatrixFixed != 'A')
158 158
     {
159 159
         ASampler.update(nA, params.maxThreads);
160
-        if (!params.useFixedMatrix || params.whichFixedMatrix != 'P')
160
+        if (params.whichMatrixFixed != 'P')
161 161
         {
162 162
             PSampler.sync(ASampler, params.maxThreads);
163 163
         }
164 164
     }
165 165
 
166
-    if (!params.useFixedMatrix || params.whichFixedMatrix != 'P')
166
+    if (params.whichMatrixFixed != 'P')
167 167
     {
168 168
         PSampler.update(nP, params.maxThreads);
169
-        if (!params.useFixedMatrix || params.whichFixedMatrix != 'A')
169
+        if (params.whichMatrixFixed != 'A')
170 170
         {
171 171
             ASampler.sync(PSampler, params.maxThreads);
172 172
         }
Browse code

flush print

Tom Sherman authored on 23/11/2018 02:01:15
Showing1 changed files
... ...
@@ -373,6 +373,7 @@ const DataType &uncertainty, GapsRandomState *randState)
373 373
         GapsTime elapsed(static_cast<unsigned>(diff.total_seconds()));
374 374
         gaps_printf("    worker %d is finished! Time: %02d:%02d:%02d\n",
375 375
             params.workerID, elapsed.hours, elapsed.minutes, elapsed.seconds);
376
+        gaps_flush();
376 377
     }
377 378
 
378 379
     return result;
Browse code

better parameter use

Tom Sherman authored on 23/11/2018 00:35:31
Showing1 changed files
... ...
@@ -367,7 +367,7 @@ const DataType &uncertainty, GapsRandomState *randState)
367 367
     result.meanChiSq = stats.meanChiSq(PSampler);
368 368
 
369 369
     // if we are running distributed, each worker needs to print when it's done
370
-    if (!params.printThreadUsage)
370
+    if (params.runningDistributed)
371 371
     {
372 372
         bpt::time_duration diff = bpt_now() - startTime;
373 373
         GapsTime elapsed(static_cast<unsigned>(diff.total_seconds()));
Browse code

added time to worker message

Tom Sherman authored on 16/11/2018 22:39:31
Showing1 changed files
... ...
@@ -23,6 +23,22 @@ namespace bpt = boost::posix_time;
23 23
         } \
24 24
     } while(0)
25 25
 
26
+struct GapsTime
27
+{
28
+    unsigned hours;
29
+    unsigned minutes;
30
+    unsigned seconds;
31
+
32
+    GapsTime(unsigned totalSeconds)
33
+    {
34
+        hours = totalSeconds / 3600;
35
+        totalSeconds -= hours * 3600;
36
+        minutes = totalSeconds / 60;
37
+        totalSeconds -= minutes * 60;
38
+        seconds = totalSeconds;
39
+    }
40
+};
41
+
26 42
 // forward declaration
27 43
 template <class Sampler, class DataType>
28 44
 static GapsResult runCoGAPSAlgorithm(const DataType &data, GapsParameters &params,
... ...
@@ -115,21 +131,9 @@ char phase, unsigned iter, GapsStatistics &stats)
115 131
         bpt::time_duration diff = bpt_now() - startTime;
116 132
         double perComplete = estimatedPercentComplete(params, ASampler,
117 133
             PSampler, phase, iter);
118
-        double nSecondsCurrent = diff.total_seconds();
119
-        double nSecondsTotal = nSecondsCurrent / perComplete;
120
-
121
-        unsigned elapsedSeconds = static_cast<unsigned>(nSecondsCurrent);
122
-        unsigned totalSeconds = static_cast<unsigned>(nSecondsTotal);
123 134
 
124
-        unsigned elapsedHours = elapsedSeconds / 3600;
125
-        elapsedSeconds -= elapsedHours * 3600;
126
-        unsigned elapsedMinutes = elapsedSeconds / 60;
127
-        elapsedSeconds -= elapsedMinutes * 60;
128
-
129
-        unsigned totalHours = totalSeconds / 3600;
130
-        totalSeconds -= totalHours * 3600;
131
-        unsigned totalMinutes = totalSeconds / 60;
132
-        totalSeconds -= totalMinutes * 60;
135
+        GapsTime elapsedTime(static_cast<unsigned>(diff.total_seconds()));
136
+        GapsTime totalTime(static_cast<unsigned>(diff.total_seconds() / perComplete));
133 137
 
134 138
         float cs = PSampler.chiSq();
135 139
         unsigned nA = ASampler.nAtoms();
... ...
@@ -139,9 +143,9 @@ char phase, unsigned iter, GapsStatistics &stats)
139 143
         stats.addAtomCount(nA, nP);
140 144
 
141 145
         gaps_printf("%d of %d, Atoms: %d(%d), ChiSq: %.0f, Time: %02d:%02d:%02d / %02d:%02d:%02d\n",
142
-            iter + 1, params.nIterations, nA, nP, cs, elapsedHours,
143
-            elapsedMinutes, elapsedSeconds, totalHours, totalMinutes,
144
-            totalSeconds);
146
+            iter + 1, params.nIterations, nA, nP, cs, elapsedTime.hours,
147
+            elapsedTime.minutes, elapsedTime.seconds, totalTime.hours,
148
+            totalTime.minutes, totalTime.seconds);
145 149
         gaps_flush();
146 150
     }
147 151
 }
... ...
@@ -361,6 +365,16 @@ const DataType &uncertainty, GapsRandomState *randState)
361 365
     // get result
362 366
     GapsResult result(stats);
363 367
     result.meanChiSq = stats.meanChiSq(PSampler);
368
+
369
+    // if we are running distributed, each worker needs to print when it's done
370
+    if (!params.printThreadUsage)
371
+    {
372
+        bpt::time_duration diff = bpt_now() - startTime;
373
+        GapsTime elapsed(static_cast<unsigned>(diff.total_seconds()));
374
+        gaps_printf("    worker %d is finished! Time: %02d:%02d:%02d\n",
375
+            params.workerID, elapsed.hours, elapsed.minutes, elapsed.seconds);
376
+    }
377
+
364 378
     return result;
365 379
 }
366 380
 
Browse code

store atom and chisq history

Tom Sherman authored on 15/11/2018 21:58:54
Showing1 changed files
... ...
@@ -107,7 +107,7 @@ const Sampler &ASampler, const Sampler &PSampler, char phase, unsigned iter)
107 107
 template <class Sampler>
108 108
 static void displayStatus(const GapsParameters &params,
109 109
 const Sampler &ASampler, const Sampler &PSampler, bpt::ptime startTime,
110
-char phase, unsigned iter)
110
+char phase, unsigned iter, GapsStatistics &stats)
111 111
 {
112 112
     if (params.printMessages && params.outputFrequency > 0
113 113
     && ((iter + 1) % params.outputFrequency) == 0)
... ...
@@ -131,10 +131,17 @@ char phase, unsigned iter)
131 131
         unsigned totalMinutes = totalSeconds / 60;
132 132
         totalSeconds -= totalMinutes * 60;
133 133
 
134
+        float cs = PSampler.chiSq();
135
+        unsigned nA = ASampler.nAtoms();
136
+        unsigned nP = PSampler.nAtoms();
137
+
138
+        stats.addChiSq(cs);
139
+        stats.addAtomCount(nA, nP);
140
+
134 141
         gaps_printf("%d of %d, Atoms: %d(%d), ChiSq: %.0f, Time: %02d:%02d:%02d / %02d:%02d:%02d\n",
135
-            iter + 1, params.nIterations, ASampler.nAtoms(),
136
-            PSampler.nAtoms(), PSampler.chiSq(), elapsedHours, elapsedMinutes,
137
-            elapsedSeconds, totalHours, totalMinutes, totalSeconds);
142
+            iter + 1, params.nIterations, nA, nP, cs, elapsedHours,
143
+            elapsedMinutes, elapsedSeconds, totalHours, totalMinutes,
144
+            totalSeconds);
138 145
         gaps_flush();
139 146
     }
140 147
 }
... ...
@@ -236,7 +243,8 @@ GapsRng &rng, bpt::ptime startTime, char phase, unsigned &currentIter)
236 243
         {
237 244
             stats.update(ASampler, PSampler);
238 245
         }
239
-        displayStatus(params, ASampler, PSampler, startTime, phase, currentIter);
246
+        displayStatus(params, ASampler, PSampler, startTime, phase,
247
+            currentIter, stats);
240 248
     }
241 249
 }
242 250
 
Browse code

cleaned up version and regenerated vignette

Tom Sherman authored on 02/11/2018 20:05:05
Showing1 changed files
... ...
@@ -49,11 +49,8 @@ const DataType &uncertainty, GapsRandomState *randState)
49 49
         return runCoGAPSAlgorithm<SparseGibbsSampler>(data, params,
50 50
             uncertainty, randState);
51 51
     }
52
-    else
53
-    {
54
-        return runCoGAPSAlgorithm<DenseGibbsSampler>(data, params,
55
-            uncertainty, randState);
56
-    }
52
+    return runCoGAPSAlgorithm<DenseGibbsSampler>(data, params,
53
+        uncertainty, randState);
57 54
 }
58 55
 
59 56
 // these two functions are the top-level functions exposed to the C++
... ...
@@ -85,8 +82,7 @@ static double estimatedNumUpdates(double current, double total, float nAtoms)
85 82
 
86 83
 template <class Sampler>
87 84
 static double estimatedPercentComplete(const GapsParameters &params,
88
-const Sampler &ASampler, const Sampler &PSampler, bpt::ptime startTime,
89
-char phase, unsigned iter)
85
+const Sampler &ASampler, const Sampler &PSampler, char phase, unsigned iter)
90 86
 {
91 87
     double nIter = static_cast<double>(iter);
92 88
     double nAtomsA = static_cast<double>(ASampler.nAtoms());
... ...
@@ -118,7 +114,7 @@ char phase, unsigned iter)
118 114
     {
119 115
         bpt::time_duration diff = bpt_now() - startTime;
120 116
         double perComplete = estimatedPercentComplete(params, ASampler,
121
-            PSampler, startTime, phase, iter);
117
+            PSampler, phase, iter);
122 118
         double nSecondsCurrent = diff.total_seconds();
123 119
         double nSecondsTotal = nSecondsCurrent / perComplete;
124 120
 
... ...
@@ -135,7 +131,7 @@ char phase, unsigned iter)
135 131
         unsigned totalMinutes = totalSeconds / 60;
136 132
         totalSeconds -= totalMinutes * 60;
137 133
 
138
-        gaps_printf("%d of %d, Atoms: %lu(%lu), ChiSq: %.0f, Time: %02d:%02d:%02d / %02d:%02d:%02d\n",
134
+        gaps_printf("%d of %d, Atoms: %d(%d), ChiSq: %.0f, Time: %02d:%02d:%02d / %02d:%02d:%02d\n",
139 135
             iter + 1, params.nIterations, ASampler.nAtoms(),
140 136
             PSampler.nAtoms(), PSampler.chiSq(), elapsedHours, elapsedMinutes,
141 137
             elapsedSeconds, totalHours, totalMinutes, totalSeconds);
Browse code

updated config to commit file permissions

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

tests passing

Tom Sherman authored on 25/10/2018 22:31:44
Showing1 changed files
... ...
@@ -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
Browse code

all tests passing besides checkpoints

Tom Sherman authored on 23/10/2018 23:49:14
Showing1 changed files
... ...
@@ -139,19 +139,19 @@ template <class Sampler>
139 139
 static void updateSampler(const GapsParameters &params, Sampler &ASampler,
140 140
 Sampler &PSampler, unsigned nA, unsigned nP)
141 141
 {
142
-    if (params.whichFixedMatrix != 'A')