@@ -133,13 +133,13 @@ static double estimatedNumUpdates(double current, double total, float nAtoms)

 template <class Sampler>

 static double estimatedPercentComplete(const GapsParameters &params,

-const Sampler &ASampler, const Sampler &PSampler, char phase, unsigned iter)

+const Sampler &ASampler, const Sampler &PSampler, GapsAlgorithmPhase phase, unsigned iter)

 {

     double nIter = static_cast<double>(iter);

     double nAtomsA = static_cast<double>(ASampler.nAtoms());

     double nAtomsP = static_cast<double>(PSampler.nAtoms());

-    if (phase == 'S')

+    if (phase == GAPS_SAMPLING_PHASE)

     {

         nIter += params.nIterations;

     }

@@ -155,7 +155,7 @@ const Sampler &ASampler, const Sampler &PSampler, char phase, unsigned iter)

 template <class Sampler>

 static void displayStatus(const GapsParameters &params,

 const Sampler &ASampler, const Sampler &PSampler, bpt::ptime startTime,

-char phase, unsigned iter, GapsStatistics &stats)

+GapsAlgorithmPhase phase, unsigned iter, GapsStatistics &stats)

 {

     if (params.outputFrequency > 0 && ((iter + 1) % params.outputFrequency) == 0)

     {

@@ -208,7 +208,7 @@ Sampler &PSampler, unsigned nA, unsigned nP)

 template <class Sampler>

 static void createCheckpoint(const GapsParameters &params,

 Sampler &ASampler, Sampler &PSampler, const GapsRandomState *randState,

-const GapsStatistics &stats, const GapsRng &rng, char phase, unsigned iter)

+const GapsStatistics &stats, const GapsRng &rng, GapsAlgorithmPhase phase, unsigned iter)

 {

     if (params.checkpointInterval > 0 && ((iter + 1) % params.checkpointInterval) == 0

     && !params.subsetData)

@@ -221,7 +221,7 @@ const GapsStatistics &stats, const GapsRng &rng, char phase, unsigned iter)

         Archive ar(params.checkpointOutFile, ARCHIVE_WRITE);

         ar << params;

         ar << *randState;

-        ar << ASampler << PSampler << stats << phase << iter << rng;

+        ar << ASampler << PSampler << stats << static_cast<int>(phase) << iter << rng;

         // delete backup file

         std::remove((params.checkpointOutFile + ".backup").c_str());

@@ -239,22 +239,24 @@ const GapsStatistics &stats, const GapsRng &rng, char phase, unsigned iter)

 template <class Sampler>

 static void processCheckpoint(GapsParameters &params, Sampler &ASampler,

 Sampler &PSampler, GapsRandomState *randState, GapsStatistics &stats,

-GapsRng &rng, char &phase, unsigned &currentIter)

+GapsRng &rng, GapsAlgorithmPhase &phase, unsigned &currentIter)

 {    

     // check if running from checkpoint, get all saved data

     if (params.useCheckPoint)

     {

+        int iPhase;

         Archive ar(params.checkpointFile, ARCHIVE_READ);

         ar >> params;

         ar >> *randState;

-        ar >> ASampler >> PSampler >> stats >> phase >> currentIter >> rng;

+        ar >> ASampler >> PSampler >> stats >> iPhase >> currentIter >> rng;

+        phase = static_cast<GapsAlgorithmPhase>(iPhase);

     }

 }

 template <class Sampler>

 static uint64_t runOnePhase(const GapsParameters &params, Sampler &ASampler,

 Sampler &PSampler, GapsStatistics &stats, const GapsRandomState *randState,

-GapsRng &rng, bpt::ptime startTime, char phase, unsigned &currentIter)

+GapsRng &rng, bpt::ptime startTime, GapsAlgorithmPhase phase, unsigned &currentIter)

 {

     uint64_t totalUpdates = 0;

     for (; currentIter < params.nIterations; ++currentIter)

@@ -264,7 +266,7 @@ GapsRng &rng, bpt::ptime startTime, char phase, unsigned &currentIter)

             rng, phase, currentIter);

         // set annealing temperature in Equilibration phase

-        if (phase == 'C')

+        if (phase == GAPS_EQUILIBRATION_PHASE)

         {        

             float temp = static_cast<float>(2 * currentIter)

                 / static_cast<float>(params.nIterations);

@@ -278,16 +280,19 @@ GapsRng &rng, bpt::ptime startTime, char phase, unsigned &currentIter)

         updateSampler(params, ASampler, PSampler, nA, nP);

         totalUpdates += nA + nP;

-        if (phase == 'S')

+        if (phase == GAPS_SAMPLING_PHASE)

         {

             stats.update(ASampler, PSampler);

             if (params.takePumpSamples)

             {

                 stats.updatePump(ASampler);

             }

+        }

+        if (params.snapshotPhase == phase || params.snapshotPhase == GAPS_ALL_PHASES)

+        {

             if (params.snapshotFrequency > 0 && ((currentIter + 1) % params.snapshotFrequency) == 0)

             {

-                stats.takeSnapshot();

+                stats.takeSnapshot(phase, ASampler, PSampler);

             }

         }

         displayStatus(params, ASampler, PSampler, startTime, phase,

@@ -405,7 +410,7 @@ const DataType &uncertainty, GapsRandomState *randState)

     // these variables will get overwritten by checkpoint if provided

     GapsStatistics stats(params.nGenes, params.nSamples, params.nPatterns);

     GapsRng rng(randState);

-    char phase = 'C';

+    GapsAlgorithmPhase phase(GAPS_EQUILIBRATION_PHASE);

     unsigned currentIter = 0;

     processCheckpoint(params, ASampler, PSampler, randState, stats, rng, phase, currentIter);

     calculateNumberOfThreads(params);

@@ -420,18 +425,18 @@ const DataType &uncertainty, GapsRandomState *randState)

     bpt::ptime startTime = bpt_now();

     // fallthrough through phases, allows algorithm to be resumed in either phase

-    GAPS_ASSERT(phase == 'C' || phase == 'S');

+    GAPS_ASSERT(phase == GAPS_EQUILIBRATION_PHASE || phase == GAPS_SAMPLING_PHASE);

     uint64_t totalUpdates = 0;

     switch (phase)

     {

-        case 'C':

+        case GAPS_EQUILIBRATION_PHASE:

             GAPS_MESSAGE(params.printMessages, "-- Equilibration Phase --\n");

             totalUpdates += runOnePhase(params, ASampler, PSampler, stats, randState,

                 rng, startTime, phase, currentIter);

-            phase = 'S';

+            phase = GAPS_SAMPLING_PHASE;

             currentIter = 0;

         // fall through

-        case 'S':

+        case GAPS_SAMPLING_PHASE:

             GAPS_MESSAGE(params.printMessages, "-- Sampling Phase --\n");

             totalUpdates += runOnePhase(params, ASampler, PSampler, stats, randState,

                 rng, startTime, phase, currentIter);

@@ -285,6 +285,10 @@ GapsRng &rng, bpt::ptime startTime, char phase, unsigned &currentIter)

             {

                 stats.updatePump(ASampler);

             }

+            if (params.snapshotFrequency > 0 && ((currentIter + 1) % params.snapshotFrequency) == 0)

+            {

+                stats.takeSnapshot();

+            }

         }

         displayStatus(params, ASampler, PSampler, startTime, phase,

             currentIter, stats);

@@ -392,7 +392,7 @@ const DataType &uncertainty, GapsRandomState *randState)

     }

     // if we are running distributed, each worker needs to print when it's started

-    if (params.runningDistributed && params.printMessages)

+    if (params.runningDistributed)

     {

         gaps_printf("    worker %d is starting!\n", params.workerID);

         gaps_flush();

@@ -449,7 +449,7 @@ const DataType &uncertainty, GapsRandomState *randState)

     }

     // if we are running distributed, each worker needs to print when it's done

-    if (params.runningDistributed && params.printMessages)

+    if (params.runningDistributed)

     {

         GapsTime elapsed(static_cast<unsigned>((bpt_now() - startTime).total_seconds()));

         gaps_printf("    worker %d is finished! Time: %02d:%02d:%02d\n",

@@ -229,7 +229,7 @@ GapsRng &rng, bpt::ptime startTime, char phase, unsigned &currentIter)

         createCheckpoint(params, ASampler, PSampler, randState, stats,

             rng, phase, currentIter);

-        // set annealing temperature in calibration phase

+        // set annealing temperature in Equilibration phase

         if (phase == 'C')

         {        

             float temp = static_cast<float>(2 * currentIter)

@@ -386,7 +386,7 @@ const DataType &uncertainty, GapsRandomState *randState)

     switch (phase)

     {

         case 'C':

-            GAPS_MESSAGE(params.printMessages, "-- Calibration Phase --\n");

+            GAPS_MESSAGE(params.printMessages, "-- Equilibration Phase --\n");

             runOnePhase(params, ASampler, PSampler, stats, randState, rng,

                 startTime, phase, currentIter);

             phase = 'S';

@@ -426,7 +426,7 @@ const DataType &uncertainty, GapsRandomState *randState)

                 rng, startTime, phase, currentIter);

             phase = 'S';

             currentIter = 0;

-

+        // fall through

         case 'S':

             GAPS_MESSAGE(params.printMessages, "-- Sampling Phase --\n");

             totalUpdates += runOnePhase(params, ASampler, PSampler, stats, randState,

@@ -1,6 +1,10 @@

 #include "GapsRunner.h"

-

+#include "GapsResult.h"

+#include "GapsParameters.h"

+#include "GapsStatistics.h"

+#include "math/Random.h"

 #include "utils/Archive.h"

+#include "utils/GlobalConfig.h"

 #include "gibbs_sampler/AsynchronousGibbsSampler.h"

 #include "gibbs_sampler/SingleThreadedGibbsSampler.h"

 #include "gibbs_sampler/DenseNormalModel.h"

@@ -39,7 +43,7 @@ struct GapsTime

     unsigned hours;

     unsigned minutes;

     unsigned seconds;

-    GapsTime(unsigned totalSeconds)

+    explicit GapsTime(unsigned totalSeconds)

     {

         hours = totalSeconds / 3600;

         totalSeconds -= hours * 3600;

@@ -39,7 +39,6 @@ struct GapsTime

     unsigned hours;

     unsigned minutes;

     unsigned seconds;

-

     GapsTime(unsigned totalSeconds)

     {

         hours = totalSeconds / 3600;

@@ -154,28 +154,28 @@ static void displayStatus(const GapsParameters &params,

 const Sampler &ASampler, const Sampler &PSampler, bpt::ptime startTime,

 char phase, unsigned iter, GapsStatistics &stats)

 {

-    if (params.printMessages && params.outputFrequency > 0

-    && ((iter + 1) % params.outputFrequency) == 0)

+    if (params.outputFrequency > 0 && ((iter + 1) % params.outputFrequency) == 0)

     {

-        bpt::time_duration diff = bpt_now() - startTime;

-        double perComplete = estimatedPercentComplete(params, ASampler,

-            PSampler, phase, iter);

-

-        GapsTime elapsedTime(static_cast<unsigned>(diff.total_seconds()));

-        GapsTime totalTime(static_cast<unsigned>(diff.total_seconds() / perComplete));

-

         float cs = PSampler.chiSq();

         unsigned nA = ASampler.nAtoms();

         unsigned nP = PSampler.nAtoms();

-

         stats.addChiSq(cs);

         stats.addAtomCount(nA, nP);

-

-        gaps_printf("%d of %d, Atoms: %d(A), %d(P), ChiSq: %.0f, Time: %02d:%02d:%02d / %02d:%02d:%02d\n",

-            iter + 1, params.nIterations, nA, nP, cs, elapsedTime.hours,

-            elapsedTime.minutes, elapsedTime.seconds, totalTime.hours,

-            totalTime.minutes, totalTime.seconds);

-        gaps_flush();

+        if (params.printMessages)

+        {

+            bpt::time_duration diff = bpt_now() - startTime;

+            double perComplete = estimatedPercentComplete(params, ASampler,

+                PSampler, phase, iter);

+

+            GapsTime elapsedTime(static_cast<unsigned>(diff.total_seconds()));

+            GapsTime totalTime(static_cast<unsigned>(diff.total_seconds() / perComplete));

+

+            gaps_printf("%d of %d, Atoms: %d(A), %d(P), ChiSq: %.0f, Time: %02d:%02d:%02d / %02d:%02d:%02d\n",

+                iter + 1, params.nIterations, nA, nP, cs, elapsedTime.hours,

+                elapsedTime.minutes, elapsedTime.seconds, totalTime.hours,

+                totalTime.minutes, totalTime.seconds);

+            gaps_flush();

+        }

     }

 }

@@ -63,11 +63,11 @@ const DataType &uncertainty, GapsRandomState *randState)

 {

     if (params.asynchronousUpdates)

     {

-        gaps_printf("Sampler Type: Asynchronous\n");

+        GAPS_MESSAGE(params.printMessages, "Sampler Type: Asynchronous\n");

         return runCoGAPSAlgorithm< AsynchronousGibbsSampler<DataModel> >(data,

             params, uncertainty, randState);

     }

-    gaps_printf("Sampler Type: Sequential\n");

+    GAPS_MESSAGE(params.printMessages, "Sampler Type: Sequential\n");

     return runCoGAPSAlgorithm< SingleThreadedGibbsSampler<DataModel> >(data,

         params, uncertainty, randState);

 }

@@ -78,10 +78,10 @@ const DataType &uncertainty, GapsRandomState *randState)

 {

     if (params.useSparseOptimization)

     {

-        gaps_printf("Data Model: Sparse, Normal\n");

+        GAPS_MESSAGE(params.printMessages, "Data Model: Sparse, Normal\n");

         return chooseSampler<SparseNormalModel>(data, params, uncertainty, randState);

     }

-    gaps_printf("Data Model: Dense, Normal\n");        

+    GAPS_MESSAGE(params.printMessages, "Data Model: Dense, Normal\n");        

     return chooseSampler<DenseNormalModel>(data, params, uncertainty, randState);

 }

@@ -296,9 +296,9 @@ Sampler &PSampler)

     // check if we're fixing a matrix

     if (params.useFixedPatterns)

     {

+        GAPS_ASSERT(params.fixedPatterns.nCol() == params.nPatterns);

         switch (params.whichMatrixFixed)

         {

-            GAPS_ASSERT(params.fixedPatterns.nCol() == params.nPatterns);

             case 'A' :

                 GAPS_ASSERT(params.fixedPatterns.nRow() == params.nGenes);

                 ASampler.setMatrix(params.fixedPatterns);

@@ -400,6 +400,7 @@ const DataType &uncertainty, GapsRandomState *randState)

     // get result

     GapsResult result(stats);

+    result.totalRunningTime = static_cast<unsigned>((bpt_now() - startTime).total_seconds());

     result.meanChiSq = stats.meanChiSq(PSampler);

     result.averageQueueLengthA = ASampler.getAverageQueueLength();

     result.averageQueueLengthP = PSampler.getAverageQueueLength();

@@ -7,7 +7,10 @@

 #include "gibbs_sampler/SparseNormalModel.h"

 #ifdef __GAPS_R_BUILD__

+#pragma GCC diagnostic push

+#pragma GCC diagnostic ignored "-Wdeprecated-declarations"

 #include <Rcpp.h>

+#pragma GCC diagnostic pop

 #endif

 #ifdef __GAPS_OPENMP__

@@ -15,7 +18,10 @@

 #endif

 // boost time helpers

+#pragma GCC diagnostic push

+#pragma GCC diagnostic ignored "-Wdeprecated-declarations"

 #include <boost/date_time/posix_time/posix_time.hpp>

+#pragma GCC diagnostic pop

 namespace bpt = boost::posix_time;

 #define bpt_now() bpt::microsec_clock::local_time()

@@ -243,10 +243,11 @@ GapsRng &rng, char &phase, unsigned &currentIter)

 }

 template <class Sampler>

-static void runOnePhase(const GapsParameters &params, Sampler &ASampler,

+static uint64_t runOnePhase(const GapsParameters &params, Sampler &ASampler,

 Sampler &PSampler, GapsStatistics &stats, const GapsRandomState *randState,

 GapsRng &rng, bpt::ptime startTime, char phase, unsigned &currentIter)

 {

+    uint64_t totalUpdates = 0;

     for (; currentIter < params.nIterations; ++currentIter)

     {

         gaps_check_interrupt();

@@ -266,6 +267,7 @@ GapsRng &rng, bpt::ptime startTime, char phase, unsigned &currentIter)

         unsigned nA = rng.poisson(gaps::max(ASampler.nAtoms(), 10));

         unsigned nP = rng.poisson(gaps::max(PSampler.nAtoms(), 10));

         updateSampler(params, ASampler, PSampler, nA, nP);

+        totalUpdates += nA + nP;

         if (phase == 'S')

         {

@@ -278,6 +280,7 @@ GapsRng &rng, bpt::ptime startTime, char phase, unsigned &currentIter)

         displayStatus(params, ASampler, PSampler, startTime, phase,

             currentIter, stats);

     }

+    return totalUpdates;

 }

 template <class Sampler>

@@ -405,19 +408,20 @@ const DataType &uncertainty, GapsRandomState *randState)

     // fallthrough through phases, allows algorithm to be resumed in either phase

     GAPS_ASSERT(phase == 'C' || phase == 'S');

+    uint64_t totalUpdates = 0;

     switch (phase)

     {

         case 'C':

             GAPS_MESSAGE(params.printMessages, "-- Calibration Phase --\n");

-            runOnePhase(params, ASampler, PSampler, stats, randState, rng,

-                startTime, phase, currentIter);

+            totalUpdates += runOnePhase(params, ASampler, PSampler, stats, randState,

+                rng, startTime, phase, currentIter);

             phase = 'S';

             currentIter = 0;

         case 'S':

             GAPS_MESSAGE(params.printMessages, "-- Sampling Phase --\n");

-            runOnePhase(params, ASampler, PSampler, stats, randState, rng,

-                startTime, phase, currentIter);

+            totalUpdates += runOnePhase(params, ASampler, PSampler, stats, randState,

+                rng, startTime, phase, currentIter);

     }

     // get result

@@ -426,6 +430,7 @@ const DataType &uncertainty, GapsRandomState *randState)

     result.meanChiSq = stats.meanChiSq(PSampler);

     result.averageQueueLengthA = ASampler.getAverageQueueLength();

     result.averageQueueLengthP = PSampler.getAverageQueueLength();

+    result.totalUpdates = totalUpdates;

     // handle pump statistics

     if (params.takePumpSamples)

@@ -442,7 +447,6 @@ const DataType &uncertainty, GapsRandomState *randState)

             params.workerID, elapsed.hours, elapsed.minutes, elapsed.seconds);

         gaps_flush();

     }

-

     return result;

 }

@@ -422,6 +422,7 @@ const DataType &uncertainty, GapsRandomState *randState)

     // get result

     GapsResult result(stats);

+    result.totalRunningTime = static_cast<unsigned>((bpt_now() - startTime).total_seconds());

     result.meanChiSq = stats.meanChiSq(PSampler);

     result.averageQueueLengthA = ASampler.getAverageQueueLength();

     result.averageQueueLengthP = PSampler.getAverageQueueLength();

@@ -436,8 +437,7 @@ const DataType &uncertainty, GapsRandomState *randState)

     // if we are running distributed, each worker needs to print when it's done

     if (params.runningDistributed)

     {

-        bpt::time_duration diff = bpt_now() - startTime;

-        GapsTime elapsed(static_cast<unsigned>(diff.total_seconds()));

+        GapsTime elapsed(static_cast<unsigned>((bpt_now() - startTime).total_seconds()));

         gaps_printf("    worker %d is finished! Time: %02d:%02d:%02d\n",

             params.workerID, elapsed.hours, elapsed.minutes, elapsed.seconds);

         gaps_flush();

@@ -145,7 +145,7 @@ char phase, unsigned iter, GapsStatistics &stats)

         stats.addChiSq(cs);

         stats.addAtomCount(nA, nP);

-        gaps_printf("%d of %d, Atoms: %d(%d), ChiSq: %.0f, Time: %02d:%02d:%02d / %02d:%02d:%02d\n",

+        gaps_printf("%d of %d, Atoms: %d(A), %d(P), ChiSq: %.0f, Time: %02d:%02d:%02d / %02d:%02d:%02d\n",

             iter + 1, params.nIterations, nA, nP, cs, elapsedTime.hours,

             elapsedTime.minutes, elapsedTime.seconds, totalTime.hours,

             totalTime.minutes, totalTime.seconds);

@@ -4,7 +4,6 @@

 #include "gibbs_sampler/AsynchronousGibbsSampler.h"

 #include "gibbs_sampler/SingleThreadedGibbsSampler.h"

 #include "gibbs_sampler/DenseNormalModel.h"

-#include "gibbs_sampler/DenseNormalWithUncertaintyModel.h"

 #include "gibbs_sampler/SparseNormalModel.h"

 #ifdef __GAPS_R_BUILD__

@@ -58,9 +57,11 @@ const DataType &uncertainty, GapsRandomState *randState)

 {

     if (params.asynchronousUpdates)

     {

+        gaps_printf("Sampler Type: Asynchronous\n");

         return runCoGAPSAlgorithm< AsynchronousGibbsSampler<DataModel> >(data,

             params, uncertainty, randState);

     }

+    gaps_printf("Sampler Type: Sequential\n");

     return runCoGAPSAlgorithm< SingleThreadedGibbsSampler<DataModel> >(data,

         params, uncertainty, randState);

 }

@@ -71,13 +72,11 @@ const DataType &uncertainty, GapsRandomState *randState)

 {

     if (params.useSparseOptimization)

     {

+        gaps_printf("Data Model: Sparse, Normal\n");

         return chooseSampler<SparseNormalModel>(data, params, uncertainty, randState);

     }

-    if (uncertainty.empty())

-    {

-        return chooseSampler<DenseNormalModel>(data, params, uncertainty, randState);

-    }

-    return chooseSampler<DenseNormalWithUncertaintyModel>(data, params, uncertainty, randState);

+    gaps_printf("Data Model: Dense, Normal\n");        

+    return chooseSampler<DenseNormalModel>(data, params, uncertainty, randState);

 }

 // helper function, this dispatches the correct run function depending

@@ -4,6 +4,7 @@

 #include "gibbs_sampler/AsynchronousGibbsSampler.h"

 #include "gibbs_sampler/SingleThreadedGibbsSampler.h"

 #include "gibbs_sampler/DenseNormalModel.h"

+#include "gibbs_sampler/DenseNormalWithUncertaintyModel.h"

 #include "gibbs_sampler/SparseNormalModel.h"

 #ifdef __GAPS_R_BUILD__

@@ -72,7 +73,11 @@ const DataType &uncertainty, GapsRandomState *randState)

     {

         return chooseSampler<SparseNormalModel>(data, params, uncertainty, randState);

     }

-    return chooseSampler<DenseNormalModel>(data, params, uncertainty, randState);

+    if (uncertainty.empty())

+    {

+        return chooseSampler<DenseNormalModel>(data, params, uncertainty, randState);

+    }

+    return chooseSampler<DenseNormalWithUncertaintyModel>(data, params, uncertainty, randState);

 }

 // helper function, this dispatches the correct run function depending

@@ -27,6 +27,7 @@ namespace bpt = boost::posix_time;

         } \

     } while(0)

+// for converting seconds to h:m:s

 struct GapsTime

 {

     unsigned hours;

@@ -75,7 +76,7 @@ const DataType &uncertainty, GapsRandomState *randState)

 }

 // helper function, this dispatches the correct run function depending

-// on the type of GibbsSampler needed for the given parameters

+// on the type of GibbsSampler and DataModel needed for the given parameters

 template <class DataType>

 static GapsResult run_helper(const DataType &data, GapsParameters &params,

 const DataType &uncertainty, GapsRandomState *randState)

@@ -131,13 +132,10 @@ const Sampler &ASampler, const Sampler &PSampler, char phase, unsigned iter)

     }

     double totalIter = 2.0 * static_cast<double>(params.nIterations);

-

     double estimatedCompleted = estimatedNumUpdates(nIter, nIter, nAtomsA) + 

         estimatedNumUpdates(nIter, nIter, nAtomsP);

-

     double estimatedTotal = estimatedNumUpdates(nIter, totalIter, nAtomsA) + 

         estimatedNumUpdates(nIter, totalIter, nAtomsP);

-

     return estimatedCompleted / estimatedTotal;

 }

@@ -199,8 +197,7 @@ static void createCheckpoint(const GapsParameters &params,

 Sampler &ASampler, Sampler &PSampler, const GapsRandomState *randState,

 const GapsStatistics &stats, const GapsRng &rng, char phase, unsigned iter)

 {

-    if (params.checkpointInterval > 0

-    && ((iter + 1) % params.checkpointInterval) == 0

+    if (params.checkpointInterval > 0 && ((iter + 1) % params.checkpointInterval) == 0

     && !params.subsetData)

     {

         // create backup file

@@ -216,6 +213,11 @@ const GapsStatistics &stats, const GapsRng &rng, char phase, unsigned iter)

         // delete backup file

         std::remove((params.checkpointOutFile + ".backup").c_str());

+        // running the extra initialization here allows for consistency with runs

+        // started from a checkpoint. This initialization phase will be run first 

+        // thing once a checkpoint is loaded since large matrices which aren't stored

+        // need to be initialized. By running it here we make sure that the algorithm

+        // is in the same state it will be when started from a checkpoint

         ASampler.extraInitialization();

         PSampler.extraInitialization();

     }

@@ -332,7 +334,7 @@ const DataType &uncertainty, GapsRandomState *randState)

 {

     // check if running in debug mode

     #ifdef GAPS_DEBUG

-    GAPS_MESSAGE(params.printMessages, "Running in debug mode\n");

+    gaps_printf("Running in debug mode\n");

     #endif

     // load data into gibbs samplers

@@ -346,7 +348,6 @@ const DataType &uncertainty, GapsRandomState *randState)

     // the process or waiting for it to finish

     GAPS_MESSAGE(params.printMessages, "Loading Data...");

     bpt::ptime readStart = bpt_now();

-

     gaps_check_interrupt();

     Sampler ASampler(data, !params.transposeData, !params.subsetGenes,

         params.alphaA, params.maxGibbsMassA, params, randState);

@@ -359,12 +360,6 @@ const DataType &uncertainty, GapsRandomState *randState)

     processFixedMatrix(params, ASampler, PSampler);

     gaps_check_interrupt();

-    // check if data is sparse and sparseOptimization is not enabled

-    if (params.printMessages && !params.useSparseOptimization && ASampler.dataSparsity() > 0.80f)

-    {

-        gaps_printf("\nWarning: data is more than 80%% sparse and sparseOptimization is not enabled\n");

-    }

-

     // elapsed time for reading data

     bpt::time_duration readDiff = bpt_now() - readStart;

     GapsTime elapsed(static_cast<unsigned>(readDiff.total_seconds()));

@@ -374,6 +369,12 @@ const DataType &uncertainty, GapsRandomState *randState)

             elapsed.seconds);

     }

+    // check if data is sparse and sparseOptimization is not enabled

+    if (params.printMessages && !params.useSparseOptimization && ASampler.dataSparsity() > 0.80f)

+    {

+        gaps_printf("\nWarning: data is more than 80%% sparse and sparseOptimization is not enabled\n");

+    }

+

     // if we are running distributed, each worker needs to print when it's started

     if (params.runningDistributed)

     {

@@ -386,8 +387,7 @@ const DataType &uncertainty, GapsRandomState *randState)

     GapsRng rng(randState);

     char phase = 'C';

     unsigned currentIter = 0;

-    processCheckpoint(params, ASampler, PSampler, randState, stats, rng,

-        phase, currentIter);

+    processCheckpoint(params, ASampler, PSampler, randState, stats, rng, phase, currentIter);

     calculateNumberOfThreads(params);

     // sync samplers and run any additional initialization needed

@@ -422,6 +422,7 @@ const DataType &uncertainty, GapsRandomState *randState)

     result.averageQueueLengthA = ASampler.getAverageQueueLength();

     result.averageQueueLengthP = PSampler.getAverageQueueLength();

+    // handle pump statistics

     if (params.takePumpSamples)

     {

         result.pumpMatrix = stats.pumpMatrix();

@@ -1,10 +1,10 @@

 #include "GapsRunner.h"

 #include "utils/Archive.h"

-#include "gibbs_sampler/GibbsSampler.h"

+#include "gibbs_sampler/AsynchronousGibbsSampler.h"

 #include "gibbs_sampler/SingleThreadedGibbsSampler.h"

-#include "gibbs_sampler/DenseStoragePolicy.h"

-#include "gibbs_sampler/SparseStoragePolicy.h"

+#include "gibbs_sampler/DenseNormalModel.h"

+#include "gibbs_sampler/SparseNormalModel.h"

 #ifdef __GAPS_R_BUILD__

 #include <Rcpp.h>

@@ -50,28 +50,28 @@ static GapsResult runCoGAPSAlgorithm(const DataType &data, GapsParameters &param

 ////////////////////////////////////////////////////////////////////////////////

-template <class StorageType, class DataType>

+template <class DataModel, class DataType>

 static GapsResult chooseSampler(const DataType &data, GapsParameters &params,

 const DataType &uncertainty, GapsRandomState *randState)

 {

     if (params.asynchronousUpdates)

     {

-        return runCoGAPSAlgorithm< GibbsSampler<StorageType> >(data,

+        return runCoGAPSAlgorithm< AsynchronousGibbsSampler<DataModel> >(data,

             params, uncertainty, randState);

     }

-    return runCoGAPSAlgorithm< SingleThreadedGibbsSampler<StorageType> >(data,

+    return runCoGAPSAlgorithm< SingleThreadedGibbsSampler<DataModel> >(data,

         params, uncertainty, randState);

 }

 template <class DataType>

-static GapsResult chooseStorage(const DataType &data, GapsParameters &params,

+static GapsResult chooseDataModel(const DataType &data, GapsParameters &params,

 const DataType &uncertainty, GapsRandomState *randState)

 {

     if (params.useSparseOptimization)

     {

-        return chooseSampler<SparseStorage>(data, params, uncertainty, randState);

+        return chooseSampler<SparseNormalModel>(data, params, uncertainty, randState);

     }

-    return chooseSampler<DenseStorage>(data, params, uncertainty, randState);

+    return chooseSampler<DenseNormalModel>(data, params, uncertainty, randState);

 }

 // helper function, this dispatches the correct run function depending

@@ -87,7 +87,7 @@ const DataType &uncertainty, GapsRandomState *randState)

         ar >> params;

         ar >> *randState;

     }

-    return chooseStorage(data, params, uncertainty, randState);

+    return chooseDataModel(data, params, uncertainty, randState);

 }

 // these two functions are the top-level functions exposed to the C++

@@ -2,6 +2,7 @@

 #include "utils/Archive.h"

 #include "gibbs_sampler/GibbsSampler.h"

+#include "gibbs_sampler/SingleThreadedGibbsSampler.h"

 #include "gibbs_sampler/DenseStoragePolicy.h"

 #include "gibbs_sampler/SparseStoragePolicy.h"

@@ -49,6 +50,30 @@ static GapsResult runCoGAPSAlgorithm(const DataType &data, GapsParameters &param

 ////////////////////////////////////////////////////////////////////////////////

+template <class StorageType, class DataType>

+static GapsResult chooseSampler(const DataType &data, GapsParameters &params,

+const DataType &uncertainty, GapsRandomState *randState)

+{

+    if (params.asynchronousUpdates)

+    {

+        return runCoGAPSAlgorithm< GibbsSampler<StorageType> >(data,

+            params, uncertainty, randState);

+    }

+    return runCoGAPSAlgorithm< SingleThreadedGibbsSampler<StorageType> >(data,

+        params, uncertainty, randState);

+}

+

+template <class DataType>

+static GapsResult chooseStorage(const DataType &data, GapsParameters &params,

+const DataType &uncertainty, GapsRandomState *randState)

+{

+    if (params.useSparseOptimization)

+    {

+        return chooseSampler<SparseStorage>(data, params, uncertainty, randState);

+    }

+    return chooseSampler<DenseStorage>(data, params, uncertainty, randState);

+}

+

 // helper function, this dispatches the correct run function depending

 // on the type of GibbsSampler needed for the given parameters

 template <class DataType>

@@ -62,14 +87,7 @@ const DataType &uncertainty, GapsRandomState *randState)

         ar >> params;

         ar >> *randState;

     }

-

-    if (params.useSparseOptimization)

-    {

-        return runCoGAPSAlgorithm< GibbsSampler<SparseStorage> >(data, params,

-            uncertainty, randState);

-    }

-    return runCoGAPSAlgorithm< GibbsSampler<DenseStorage> >(data, params,

-        uncertainty, randState);

+    return chooseStorage(data, params, uncertainty, randState);

 }

 // these two functions are the top-level functions exposed to the C++

@@ -357,7 +357,7 @@ const DataType &uncertainty, GapsRandomState *randState)

     }

     // if we are running distributed, each worker needs to print when it's started

-    if (params.runningDistributed)

+    if (params.runningDistributed && params.printMessages)

     {

         gaps_printf("    worker %d is starting!\n", params.workerID);

         gaps_flush();

@@ -411,7 +411,7 @@ const DataType &uncertainty, GapsRandomState *randState)

     }

     // if we are running distributed, each worker needs to print when it's done

-    if (params.runningDistributed)

+    if (params.runningDistributed && params.printMessages)

     {

         bpt::time_duration diff = bpt_now() - startTime;

         GapsTime elapsed(static_cast<unsigned>(diff.total_seconds()));

@@ -341,6 +341,12 @@ const DataType &uncertainty, GapsRandomState *randState)

     processFixedMatrix(params, ASampler, PSampler);

     gaps_check_interrupt();

+    // check if data is sparse and sparseOptimization is not enabled

+    if (params.printMessages && !params.useSparseOptimization && ASampler.dataSparsity() > 0.80f)

+    {

+        gaps_printf("\nWarning: data is more than 80%% sparse and sparseOptimization is not enabled\n");

+    }

+

     // elapsed time for reading data

     bpt::time_duration readDiff = bpt_now() - readStart;

     GapsTime elapsed(static_cast<unsigned>(readDiff.total_seconds()));

@@ -350,6 +350,13 @@ const DataType &uncertainty, GapsRandomState *randState)

             elapsed.seconds);

     }

+    // if we are running distributed, each worker needs to print when it's started

+    if (params.runningDistributed)

+    {

+        gaps_printf("    worker %d is starting!\n", params.workerID);

+        gaps_flush();

+    }

+

     // these variables will get overwritten by checkpoint if provided

     GapsStatistics stats(params.nGenes, params.nSamples, params.nPatterns);

     GapsRng rng(randState);

@@ -222,10 +222,7 @@ GapsRng &rng, bpt::ptime startTime, char phase, unsigned &currentIter)

 {

     for (; currentIter < params.nIterations; ++currentIter)

     {

-        #ifdef __GAPS_R_BUILD__

-        Rcpp::checkUserInterrupt();

-        #endif

-

+        gaps_check_interrupt();

         createCheckpoint(params, ASampler, PSampler, randState, stats,

             rng, phase, currentIter);

@@ -323,12 +320,17 @@ const DataType &uncertainty, GapsRandomState *randState)

     // within the sampler, note the subsetting genes/samples flag must be

     // flipped if we are flipping the transpose flag

     GAPS_MESSAGE(params.printMessages, "Loading Data...");

+    gaps_check_interrupt();

     Sampler ASampler(data, !params.transposeData, !params.subsetGenes,

         params.alphaA, params.maxGibbsMassA, params, randState);

+    gaps_check_interrupt();

     Sampler PSampler(data, params.transposeData, params.subsetGenes,

         params.alphaP, params.maxGibbsMassP, params, randState);

+    gaps_check_interrupt();

     processUncertainty(params, ASampler, PSampler, uncertainty);

+    gaps_check_interrupt();

     processFixedMatrix(params, ASampler, PSampler);

+    gaps_check_interrupt();

     GAPS_MESSAGE(params.printMessages, "Done!\n");

     // these variables will get overwritten by checkpoint if provided

@@ -369,6 +369,8 @@ const DataType &uncertainty, GapsRandomState *randState)

     // get result

     GapsResult result(stats);

     result.meanChiSq = stats.meanChiSq(PSampler);

+    result.averageQueueLengthA = ASampler.getAverageQueueLength();

+    result.averageQueueLengthP = PSampler.getAverageQueueLength();

     if (params.takePumpSamples)

     {

@@ -1,6 +1,9 @@

 #include "GapsRunner.h"

 #include "utils/Archive.h"

+#include "gibbs_sampler/GibbsSampler.h"

+#include "gibbs_sampler/DenseStoragePolicy.h"

+#include "gibbs_sampler/SparseStoragePolicy.h"

 #ifdef __GAPS_R_BUILD__

 #include <Rcpp.h>

@@ -62,10 +65,10 @@ const DataType &uncertainty, GapsRandomState *randState)

     if (params.useSparseOptimization)

     {

-        return runCoGAPSAlgorithm<SparseGibbsSampler>(data, params,

+        return runCoGAPSAlgorithm< GibbsSampler<SparseStorage> >(data, params,

             uncertainty, randState);

     }

-    return runCoGAPSAlgorithm<DenseGibbsSampler>(data, params,

+    return runCoGAPSAlgorithm< GibbsSampler<DenseStorage> >(data, params,

         uncertainty, randState);

 }

@@ -323,13 +323,20 @@ const DataType &uncertainty, GapsRandomState *randState)

     // within the sampler, note the subsetting genes/samples flag must be

     // flipped if we are flipping the transpose flag

     GAPS_MESSAGE(params.printMessages, "Loading Data...");

+    bpt::ptime readStart = bpt_now();

     Sampler ASampler(data, !params.transposeData, !params.subsetGenes,

         params.alphaA, params.maxGibbsMassA, params, randState);

     Sampler PSampler(data, params.transposeData, params.subsetGenes,

         params.alphaP, params.maxGibbsMassP, params, randState);

     processUncertainty(params, ASampler, PSampler, uncertainty);

     processFixedMatrix(params, ASampler, PSampler);

-    GAPS_MESSAGE(params.printMessages, "Done!\n");

+    bpt::time_duration readDiff = bpt_now() - readStart;