@@ -8,9 +8,7 @@

 // this file contains the blueprint for creating a wrapper around the C++

 // interface used for running CoGAPS. It exposes some functions to R, has a

 // method for converting the R parameters to the standard GapsParameters

-// struct, and creates a GapsRunner object. The GapsRunner class manages

-// all information from a CoGAPS run and is used to set off the run

-// and get return data.

+// and calls gaps::run

 ////////////////// functions for converting matrix types ///////////////////////

@@ -53,26 +51,20 @@ const Rcpp::Nullable<Rcpp::IntegerVector> &indices)

 {

     // check if subsetting data

     const Rcpp::S4 &gapsParams(allParams["gaps"]);

-    bool subsetData = false;

-    bool printThreadUsage = true;

     bool subsetGenes = false;

-    char whichFixedMatrix = 'N';

     std::vector<unsigned> subset;

     if (indices.isNotNull())

     {

-        subsetData = true;

-        printThreadUsage = false;

         std::string d(Rcpp::as<std::string>(gapsParams.slot("distributed")));

         subsetGenes = (d == "genome-wide");

-        whichFixedMatrix = (d == "genome-wide") ? 'P' : 'A';

         subset = Rcpp::as< std::vector<unsigned> >(Rcpp::IntegerVector(indices));

     }

     // create standard CoGAPS parameters struct

-    GapsParameters params(data, allParams["transposeData"], subsetData,

+    GapsParameters params(data, allParams["transposeData"], indices.isNotNull(),

         subsetGenes, subset);

-    params.printThreadUsage = printThreadUsage;

-    params.whichFixedMatrix = whichFixedMatrix;

+    params.printThreadUsage = !indices.isNotNull();

+    params.whichFixedMatrix = indices.isNotNull() ? (subsetGenes ? 'P' : 'A') : 'N';

     // get configuration parameters

     params.maxThreads = allParams["nThreads"];

@@ -104,13 +96,12 @@ const Rcpp::Nullable<Rcpp::IntegerVector> &indices)

     {

         params.checkpointFile = Rcpp::as<std::string>(allParams["checkpointInFile"]);

         params.useCheckPoint = true;

-        params.peekCheckpoint(params.checkpointFile);

     }

     return params;

 }

-////////// main function that creates a GapsRunner and runs CoGAPS /////////////

+////////////////////// main function that runs CoGAPS //////////////////////////

 // note uncertainty matrix gets special treatment since it's the same size as

 // the data (potentially large), so we want to avoid copying it into the 

@@ -122,23 +113,12 @@ const DataType &uncertainty, const Rcpp::Nullable<Rcpp::IntegerVector> &indices,

 const Rcpp::Nullable<Rcpp::NumericMatrix> &fixedMatrix, bool isMaster)

 {

     // convert R parameters to GapsParameters struct

-    GapsParameters gapsParams(getGapsParameters(data, allParams, isMaster,

+    GapsParameters params(getGapsParameters(data, allParams, isMaster,

         fixedMatrix, indices));

     // create GapsRunner, note we must first initialize the random generator

-    GapsRng::setSeed(gapsParams.seed);

-    gaps_printf("Loading Data...");

-    GapsRunner runner(data, gapsParams);

-

-    // set uncertainty

-    if (!uncertainty.empty())

-    {

-        runner.setUncertainty(uncertainty, gapsParams);

-    }

-    gaps_printf("Done!\n");

-    

-    // run cogaps

-    GapsResult result(runner.run());

+    GapsRandomState randState(params.seed);

+    GapsResult result(gaps::run(data, params, uncertainty, &randState));

     // write result to file if requested

     if (allParams["outputToFile"] != R_NilValue)

@@ -152,7 +132,7 @@ const Rcpp::Nullable<Rcpp::NumericMatrix> &fixedMatrix, bool isMaster)

         Rcpp::Named("Pmean") = createRMatrix(result.Pmean),

         Rcpp::Named("Asd") = createRMatrix(result.Asd),

         Rcpp::Named("Psd") = createRMatrix(result.Psd),

-        Rcpp::Named("seed") = gapsParams.seed,

+        Rcpp::Named("seed") = params.seed,

         Rcpp::Named("meanChiSq") = result.meanChiSq,

         Rcpp::Named("geneNames") = allParams["geneNames"],

         Rcpp::Named("sampleNames") = allParams["sampleNames"],

@@ -1,9 +1,19 @@

 #include "GapsParameters.h"

-void GapsParameters::peekCheckpoint(const std::string &file)

+Archive& operator<<(Archive &ar, const GapsParameters &p)

 {

-    Archive ar(file, ARCHIVE_READ);

-    ar >> nPatterns >> seed >> nIterations >> whichFixedMatrix;

+    ar << p.seed << p.nGenes << p.nSamples << p.nPatterns << p.nIterations

+        << p.alphaA << p.alphaP << p.maxGibbsMassA << p.maxGibbsMassP

+        << p.singleCell << p.useSparseOptimization;

+    return ar;

+}

+

+Archive& operator>>(Archive &ar, GapsParameters &p)

+{

+    ar >> p.seed >> p.nGenes >> p.nSamples >> p.nPatterns >> p.nIterations

+        >> p.alphaA >> p.alphaP >> p.maxGibbsMassA >> p.maxGibbsMassP

+        >> p.singleCell >> p.useSparseOptimization;

+    return ar;

 }

 void GapsParameters::calculateDataDimensions(const std::string &file)

@@ -16,8 +16,6 @@ public:

         bool t_subsetData=false, bool t_subsetGenes=false,

         const std::vector<unsigned> t_dataIndicesSubset=std::vector<unsigned>());

-    void peekCheckpoint(const std::string &file);

-

     Matrix fixedMatrix;

     std::vector<unsigned> dataIndicesSubset;

@@ -58,6 +56,9 @@ private:

     void calculateDataDimensions(const Matrix &mat);

 };

+Archive& operator<<(Archive &ar, const GapsParameters &p);

+Archive& operator>>(Archive &ar, GapsParameters &p);

+

 template <class DataType>

 GapsParameters::GapsParameters(const DataType &data, bool t_transposeData,

 bool t_subsetData, bool t_subsetGenes,

@@ -1,5 +1,7 @@

 #include "GapsRunner.h"

+#include "utils/Archive.h"

+

 #ifdef __GAPS_R_BUILD__

 #include <Rcpp.h>

 #endif

@@ -8,117 +10,61 @@

 #include <omp.h>

 #endif

-///////////////////////////// RAII wrapper /////////////////////////////////////

+// boost time helpers

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

+namespace bpt = boost::posix_time;

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

-GapsRunner::~GapsRunner()

-{

-    delete mRunner;

-}

+// forward declaration

+template <class Sampler, class DataType>

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

+    const DataType &uncertainty, GapsRandomState *randState);

-GapsResult GapsRunner::run()

-{

-    return mRunner->run();

-}

+////////////////////////////////////////////////////////////////////////////////

-///////////////////////// Abstract Interface ///////////////////////////////////

-

-AbstractGapsRunner::AbstractGapsRunner(const GapsParameters &params)

-    :

-mStatistics(params.nGenes, params.nSamples, params.nPatterns),

-mCheckpointOutFile(params.checkpointOutFile),

-mCurrentIteration(0),

-mMaxIterations(params.nIterations),

-mMaxThreads(params.maxThreads),

-mOutputFrequency(params.outputFrequency),

-mCheckpointInterval(params.checkpointInterval),

-mNumPatterns(params.nPatterns),

-mNumUpdatesA(0),

-mNumUpdatesP(0),

-mSeed(params.seed),

-mPrintMessages(params.printMessages),

-mPrintThreadUsage(params.printThreadUsage),

-mPhase('C'),

-mFixedMatrix(params.whichFixedMatrix)

-{}

-

-GapsResult AbstractGapsRunner::run()

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

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

+template <class DataType>

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

+const DataType &uncertainty, GapsRandomState *randState)

 {

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

-

-    mStartTime = bpt_now();

-

-    // check if running in debug mode

-    #ifdef GAPS_DEBUG

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

-    #endif

-

-    // calculate appropiate number of threads if compiled with openmp

-    #ifdef __GAPS_OPENMP__

-    if (mPrintMessages && mPrintThreadUsage)

+    // fetch parameters from checkpoint - some are used in initialization

+    if (params.useCheckPoint)

     {

-        unsigned availableThreads = omp_get_max_threads();

-        mMaxThreads = gaps::min(availableThreads, mMaxThreads);

-        gaps_printf("Running on %d out of %d available threads\n",

-            mMaxThreads, availableThreads);

+        Archive ar(params.checkpointFile, ARCHIVE_READ);

+        ar >> params;

+        ar >> *randState;

     }

-    #endif

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

-    switch (mPhase)

+    if (params.useSparseOptimization)

     {

-        case 'C':

-            if (mPrintMessages)

-            {

-                gaps_printf("-- Calibration Phase --\n");

-            }

-            runOnePhase();

-            mPhase = 'S';

-            mCurrentIteration = 0;

-

-        case 'S':

-            if (mPrintMessages)

-            {

-                gaps_printf("-- Sampling Phase --\n");

-            }

-            runOnePhase();

-            break;

+        return runCoGAPSAlgorithm<SparseGibbsSampler>(data, params,

+            uncertainty, randState);

+    }

+    else

+    {

+        return runCoGAPSAlgorithm<DenseGibbsSampler>(data, params,

+            uncertainty, randState);

     }

-    GapsResult result(mStatistics);

-    result.meanChiSq = meanChiSq();

-    return result;    

 }

-void AbstractGapsRunner::runOnePhase()

-{

-    for (; mCurrentIteration < mMaxIterations; ++mCurrentIteration)

-    {

-        createCheckpoint();

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

+// code that is being wrapped by any given language

-        #ifdef __GAPS_R_BUILD__

-        Rcpp::checkUserInterrupt();

-        #endif

-        

-        // set annealing temperature in calibration phase

-        if (mPhase == 'C')

-        {        

-            float temp = static_cast<float>(2 * mCurrentIteration)

-                / static_cast<float>(mMaxIterations);

-            setAnnealingTemp(gaps::min(1.f, temp));

-        }

-    

-        // number of updates per iteration is poisson 

-        unsigned nA = mRng.poisson(gaps::max(nAtoms('A'), 10));

-        unsigned nP = mRng.poisson(gaps::max(nAtoms('P'), 10));

-        updateSampler(nA, nP);

+GapsResult gaps::run(const Matrix &data, GapsParameters &params,

+const Matrix &uncertainty, GapsRandomState *randState)

+{

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

+}

-        if (mPhase == 'S')

-        {

-            updateStatistics();

-        }

-        displayStatus();

-    }

+GapsResult gaps::run(const std::string &data, GapsParameters &params,

+const std::string &uncertainty, GapsRandomState *randState)

+{

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

 }

+////////////////////////////////////////////////////////////////////////////////

+

 // sum coef * log(i) for i = 1 to total, fit coef from number of atoms

 // approximates sum of number of atoms (stirling approx to factorial)

 // this should be proportional to total number of updates

@@ -129,19 +75,21 @@ static double estimatedNumUpdates(double current, double total, float nAtoms)

         total * coef * std::log(total) - total * coef;

 }

-

-double AbstractGapsRunner::estimatedPercentComplete() const

+template <class Sampler>

+static double estimatedPercentComplete(const GapsParameters &params,

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

+char phase, unsigned iter)

 {

-    double nIter = static_cast<double>(mCurrentIteration);

-    double nAtomsA = static_cast<double>(nAtoms('A'));

-    double nAtomsP = static_cast<double>(nAtoms('P'));

+    double nIter = static_cast<double>(iter);

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

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

-    if (mPhase == 'S')

+    if (phase == 'S')

     {

-        nIter += mMaxIterations;

+        nIter += params.nIterations;

     }

-    double totalIter = 2.0 * static_cast<double>(mMaxIterations);

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

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

         estimatedNumUpdates(nIter, nIter, nAtomsP);

@@ -152,13 +100,19 @@ double AbstractGapsRunner::estimatedPercentComplete() const

     return estimatedCompleted / estimatedTotal;

 }

-void AbstractGapsRunner::displayStatus()

+template <class Sampler>

+static void displayStatus(const GapsParameters &params,

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

+char phase, unsigned iter)

 {

-    if (mPrintMessages && mOutputFrequency > 0 && ((mCurrentIteration + 1) % mOutputFrequency) == 0)

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

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

     {

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

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

+        double perComplete = estimatedPercentComplete(params, ASampler,

+            PSampler, startTime, phase, iter);

         double nSecondsCurrent = diff.total_seconds();

-        double nSecondsTotal = nSecondsCurrent / estimatedPercentComplete();

+        double nSecondsTotal = nSecondsCurrent / perComplete;

         unsigned elapsedSeconds = static_cast<unsigned>(nSecondsCurrent);

         unsigned totalSeconds = static_cast<unsigned>(nSecondsTotal);

@@ -174,180 +128,201 @@ void AbstractGapsRunner::displayStatus()

         totalSeconds -= totalMinutes * 60;

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

-            mCurrentIteration + 1, mMaxIterations, nAtoms('A'),

-            nAtoms('P'), chiSq(), elapsedHours, elapsedMinutes,

+            iter + 1, params.nIterations, ASampler.nAtoms(),

+            PSampler.nAtoms(), PSampler.chiSq(), elapsedHours, elapsedMinutes,

             elapsedSeconds, totalHours, totalMinutes, totalSeconds);

         gaps_flush();

     }

 }

-void AbstractGapsRunner::createCheckpoint()

+template <class Sampler>

+static void updateSampler(const GapsParameters &params, Sampler &ASampler,

+Sampler &PSampler, unsigned nA, unsigned nP)

 {

-    if (mCheckpointInterval > 0 && ((mCurrentIteration + 1) % mCheckpointInterval) == 0)

+    if (params.whichFixedMatrix != 'A')

     {

-        // create backup file

-        std::rename(mCheckpointOutFile.c_str(), (mCheckpointOutFile + ".backup").c_str());

-    

-        // create checkpoint file

-        Archive ar(mCheckpointOutFile, ARCHIVE_WRITE);

-        ar << mNumPatterns << mSeed << mMaxIterations << mFixedMatrix << mPhase

-            << mCurrentIteration << mNumUpdatesA << mNumUpdatesP << mRng;

-        writeSamplers(ar);

-        GapsRng::save(ar);

-

-        // delete backup file

-        std::remove((mCheckpointOutFile + ".backup").c_str());

-    }

-}

-

-///////////////////// DenseGapsRunner Implementation ///////////////////////////

-

-float DenseGapsRunner::chiSq() const

-{

-    // doesn't matter which sampler is called

-    return mPSampler.chiSq();

-}

-

-float DenseGapsRunner::meanChiSq() const

-{

-    // need to pass P sampler (due to configuration of internal data)

-    return mStatistics.meanChiSq(mPSampler);

-}

-

-unsigned DenseGapsRunner::nAtoms(char which) const

-{

-    return which == 'A' ? mASampler.nAtoms() : mPSampler.nAtoms();

-}

-

-void DenseGapsRunner::setAnnealingTemp(float temp)

-{

-    mASampler.setAnnealingTemp(temp);

-    mPSampler.setAnnealingTemp(temp);

-}

-

-void DenseGapsRunner::updateStatistics()

-{

-    mStatistics.update(mASampler, mPSampler);

-}

-

-Archive& DenseGapsRunner::readSamplers(Archive &ar)

-{

-    ar >> mASampler >> mPSampler;

-    return ar;

-}

-

-Archive& DenseGapsRunner::writeSamplers(Archive &ar)

-{

-    ar << mASampler << mPSampler;

-    return ar;

-}

-

-void DenseGapsRunner::updateSampler(unsigned nA, unsigned nP)

-{

-    if (mFixedMatrix != 'A')

-    {

-        mNumUpdatesA += nA;

-        mASampler.update(nA, mMaxThreads);

-        if (mFixedMatrix != 'P')

+        ASampler.update(nA, params.maxThreads);

+        if (params.whichFixedMatrix != 'P')

         {

-            mPSampler.sync(mASampler, mMaxThreads);

+            PSampler.sync(ASampler, params.maxThreads);

         }

     }

-    if (mFixedMatrix != 'P')

+    if (params.whichFixedMatrix != 'P')

     {

-        mNumUpdatesP += nP;

-        mPSampler.update(nP, mMaxThreads);

-        if (mFixedMatrix != 'A')

+        PSampler.update(nP, params.maxThreads);

+        if (params.whichFixedMatrix != 'A')

         {

-            mASampler.sync(mPSampler, mMaxThreads);

+            ASampler.sync(PSampler, params.maxThreads);

         }

     }

 }

-void DenseGapsRunner::setUncertainty(const Matrix &unc, const GapsParameters &params)

+template <class Sampler>

+static void createCheckpoint(const GapsParameters &params,

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

+const GapsRng &rng, char phase, unsigned iter)

 {

-    mASampler.setUncertainty(unc, !params.transposeData, !params.subsetGenes, params);

-    mPSampler.setUncertainty(unc, params.transposeData, params.subsetGenes, params);

+    if (params.checkpointInterval > 0

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

+    {

+        // create backup file

+        std::rename(params.checkpointOutFile.c_str(),

+            (params.checkpointOutFile + ".backup").c_str());

+    

+        // create checkpoint file

+        Archive ar(params.checkpointOutFile, ARCHIVE_WRITE);

+        ar << params;

+        ar << *randState;

+        ar << ASampler << PSampler << phase << iter << rng;

+        

+        // delete backup file

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

+    }

 }

-void DenseGapsRunner::setUncertainty(const std::string &unc, const GapsParameters &params)

+template <class Sampler>

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

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

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

 {

-    mASampler.setUncertainty(unc, !params.transposeData, !params.subsetGenes, params);

-    mPSampler.setUncertainty(unc, params.transposeData, params.subsetGenes, params);

-}

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

+    {

+        #ifdef __GAPS_R_BUILD__

+        Rcpp::checkUserInterrupt();

+        #endif

-///////////////////// SparseGapsRunner Implementation //////////////////////////

+        createCheckpoint(params, ASampler, PSampler, randState, rng, phase,

+            currentIter);

+        

+        // set annealing temperature in calibration phase

+        if (phase == 'C')

+        {        

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

+                / static_cast<float>(params.nIterations);

+            ASampler.setAnnealingTemp(gaps::min(1.f, temp));

+            PSampler.setAnnealingTemp(gaps::min(1.f, temp));

+        }

+    

+        // number of updates per iteration is poisson 

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

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

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

-float SparseGapsRunner::chiSq() const

-{

-    // doesn't matter which sampler is called

-    return mPSampler.chiSq();

+        if (phase == 'S')

+        {

+            stats.update(ASampler, PSampler);

+        }

+        displayStatus(params, ASampler, PSampler, startTime, phase, currentIter);

+    }

 }

-float SparseGapsRunner::meanChiSq() const

+// here is the CoGAPS algorithm

+template <class Sampler, class DataType>

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

+const DataType &uncertainty, GapsRandomState *randState)

 {

-    // need to pass P sampler (due to configuration of internal data)

-    return mStatistics.meanChiSq(mPSampler);

-}

+    // check if running in debug mode

+    #ifdef GAPS_DEBUG

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

+    #endif

-unsigned SparseGapsRunner::nAtoms(char which) const

-{

-    return which == 'A' ? mASampler.nAtoms() : mPSampler.nAtoms();

-}

+    // load data into gibbs samplers

+    // we transpose the data in the A sampler so that the update step

+    // is symmetrical for each sampler, this simplifies the code 

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

+    // flipped if we are flipping the transpose flag

+    gaps_printf("Loading Data...");

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

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

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

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

+

+    // read in the uncertainty matrix if one is provided

+    if (!uncertainty.empty())

+    {

+        ASampler.setUncertainty(uncertainty, !params.transposeData,

+            !params.subsetGenes, params);

+        PSampler.setUncertainty(uncertainty, params.transposeData,

+            params.subsetGenes, params);

+    }

+    gaps_printf("Done!\n");

-void SparseGapsRunner::setAnnealingTemp(float temp)

-{

-    mASampler.setAnnealingTemp(temp);

-    mPSampler.setAnnealingTemp(temp);

-}

+    // these variables will get overwritten by checkpoint if provided

+    GapsRng rng(randState);

+    char phase = 'C';

+    unsigned currentIter = 0;

-void SparseGapsRunner::updateStatistics()

-{

-    mStatistics.update(mASampler, mPSampler);

-}

+    // check if we're fixing a matrix

+    switch (params.whichFixedMatrix)

+    {

+        case 'A' : ASampler.setMatrix(params.fixedMatrix); break;

+        case 'P' : PSampler.setMatrix(params.fixedMatrix); break;

+        default: break; // 'N' for none

+    }

+     

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

+    if (params.useCheckPoint)

+    {

+        Archive ar(params.checkpointFile, ARCHIVE_READ);

+        ar >> params;

+        ar >> *randState;

+        ar >> ASampler >> PSampler >> phase >> currentIter >> rng;

+    }

-Archive& SparseGapsRunner::readSamplers(Archive &ar)

-{

-    ar >> mASampler >> mPSampler;

-    return ar;

-}

+    // sync samplers, second parameter indicates this should be a full sync

+    ASampler.sync(PSampler);

+    PSampler.sync(ASampler);

-Archive& SparseGapsRunner::writeSamplers(Archive &ar)

-{

-    ar << mASampler << mPSampler;

-    return ar;

-}

+    // sampler may need to run additional initialization after parameters set

+    ASampler.extraInitialization();

+    PSampler.extraInitialization();

-void SparseGapsRunner::updateSampler(unsigned nA, unsigned nP)

-{

-    if (mFixedMatrix != 'A')

+    // calculate appropiate number of threads if compiled with openmp

+    #ifdef __GAPS_OPENMP__

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

     {

-        mNumUpdatesA += nA;

-        mASampler.update(nA, mMaxThreads);

-        if (mFixedMatrix != 'P')

-        {

-            mPSampler.sync(mASampler, mMaxThreads);

-        }

+        unsigned availableThreads = omp_get_max_threads();

+        params.maxThreads = gaps::min(availableThreads, params.maxThreads);

+        gaps_printf("Running on %d out of %d available threads\n",

+            params.maxThreads, availableThreads);

     }

+    #endif

+

+    // record start time

+    bpt::ptime startTime = bpt_now();

-    if (mFixedMatrix != 'P')

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

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

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

+    switch (phase)

     {

-        mNumUpdatesP += nP;

-        mPSampler.update(nP, mMaxThreads);

-        if (mFixedMatrix != 'A')

-        {

-            mASampler.sync(mPSampler, mMaxThreads);

-        }

-    }

-}

+        case 'C':

+            if (params.printMessages)

+            {

+                gaps_printf("-- Calibration Phase --\n");

+            }

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

+                startTime, phase, currentIter);

+            phase = 'S';

+            currentIter = 0;

-void SparseGapsRunner::setUncertainty(const Matrix &unc, const GapsParameters &params)

-{

-    // nothing happens - SparseGibbsSampler assumes default uncertainty always

-}

-void SparseGapsRunner::setUncertainty(const std::string &unc, const GapsParameters &params)

-{

-    // nothing happens - SparseGibbsSampler assumes default uncertainty always

+

+        case 'S':

+            if (params.printMessages)

+            {

+                gaps_printf("-- Sampling Phase --\n");

+            }

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

+                startTime, phase, currentIter);

+

+        default: break;

+    }

+    

+    // get result

+    GapsResult result(stats);

+    result.meanChiSq = stats.meanChiSq(PSampler);

+    return result;

 }

+

@@ -1,265 +1,23 @@

 #ifndef __COGAPS_GAPS_RUNNER_H__

 #define __COGAPS_GAPS_RUNNER_H__

-#include "GapsParameters.h"

 #include "GapsResult.h"

-#include "GapsStatistics.h"

-#include "gibbs_sampler/GibbsSampler.h"

-#include "gibbs_sampler/DenseGibbsSampler.h"

-#include "gibbs_sampler/SparseGibbsSampler.h"

-

-#include <string>

-

-// boost time helpers

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

-namespace bpt = boost::posix_time;

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

-

-// forward declarations

-class AbstractGapsRunner;

-

-///////////////////////////// RAII wrapper /////////////////////////////////////

-

-// This is the class that is exposed to the top-level CoGAPS routine - all 

-// aspects of CoGAPS can be managed through this class. The class itself is

-// just a lightweight wrapper around an abstract interface, which allows for

-// multiple types of GapsRunner to be declared. Which implementation is used

-// depends on the parameters passed to the GapsRunner constructor.

-class GapsRunner

-{

-public:

-

-    template <class DataType>

-    GapsRunner(const DataType &data, const GapsParameters &params);

-

-    ~GapsRunner();

-

-    template <class DataType>

-    void setUncertainty(const DataType &unc, const GapsParameters &params);

-

-    GapsResult run();

-

-private:

-

-    AbstractGapsRunner *mRunner;

-

-    GapsRunner(const GapsRunner &p); // don't allow copies

-    GapsRunner& operator=(const GapsRunner &p); // don't allow copies    

-};

-

-///////////////////////// Abstract Interface ///////////////////////////////////

-

-// This class is the abstract interface that any implementation of GapsRunner

-// must satisfy. It provides a factory method that will create the appropiate

-// derived class depending on the parameters passed in.

-class AbstractGapsRunner

-{

-public:

-

-    AbstractGapsRunner(const GapsParameters &params);

-    virtual ~AbstractGapsRunner() {}

-

-    template <class DataType>

-    static AbstractGapsRunner* create(const DataType &data, const GapsParameters &params);

-

-    // can't use template with virtual function

-    virtual void setUncertainty(const Matrix &unc, const GapsParameters &params) = 0;

-    virtual void setUncertainty(const std::string &unc, const GapsParameters &params) = 0;

-

-    GapsResult run();

-

-protected:

-

-    GapsStatistics mStatistics;

-

-    mutable GapsRng mRng;

-

-    std::string mCheckpointOutFile;

-

-    bpt::ptime mStartTime;

-

-    unsigned mCurrentIteration;

-    unsigned mMaxIterations;

-    unsigned mMaxThreads;

-    unsigned mOutputFrequency;

-    unsigned mCheckpointInterval;

-    unsigned mNumPatterns;

-    unsigned mNumUpdatesA;

-    unsigned mNumUpdatesP;

-    uint32_t mSeed;

-

-    bool mPrintMessages;

-    bool mPrintThreadUsage;

-

-    char mPhase;

-    char mFixedMatrix;

-        

-    void runOnePhase();

-    double estimatedPercentComplete() const;

-    void displayStatus();

-    void createCheckpoint();

-

-    virtual float chiSq() const = 0;

-    virtual float meanChiSq() const = 0;

-    virtual unsigned nAtoms(char which) const = 0;

-    virtual void setAnnealingTemp(float temp) = 0;

-    virtual void updateStatistics() = 0;

-    virtual Archive& readSamplers(Archive &ar) = 0;

-    virtual Archive& writeSamplers(Archive &ar) = 0;

-    virtual void updateSampler(unsigned nA, unsigned nP) = 0;

-};

-

-///////////////////// GapsRunner Implementations ///////////////////////////////

-

-// This implementation uses a DenseGibbsSampler internally

-class DenseGapsRunner : public AbstractGapsRunner

-{

-public:

-

-    ~DenseGapsRunner() {}

-

-    template <class DataType>

-    DenseGapsRunner(const DataType &data, const GapsParameters &params);

-

-    void setUncertainty(const Matrix &unc, const GapsParameters &params);

-    void setUncertainty(const std::string &unc, const GapsParameters &params);

-

-private:

-

-    DenseGibbsSampler mASampler;

-    DenseGibbsSampler mPSampler;

+#include "GapsParameters.h"

+#include "data_structures/Matrix.h"

+#include "math/Random.h"

-    float chiSq() const;

-    float meanChiSq() const;

-    unsigned nAtoms(char which) const;

-    void setAnnealingTemp(float temp);

-    void updateStatistics();

-    Archive& readSamplers(Archive &ar);

-    Archive& writeSamplers(Archive &ar);

-    void updateSampler(unsigned nA, unsigned nP);

-};

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

+// code that is being wrapped by any given language

-// This implementation uses a SparseGibbsSampler internally

-class SparseGapsRunner : public AbstractGapsRunner

+namespace gaps

 {

-public:

-

-    ~SparseGapsRunner() {}

-

-    template <class DataType>

-    SparseGapsRunner(const DataType &data, const GapsParameters &params);

-

-    void setUncertainty(const Matrix &unc, const GapsParameters &params);

-    void setUncertainty(const std::string &unc, const GapsParameters &params);

+    // data stored in matrix

+    GapsResult run(const Matrix &data, GapsParameters &params,

+        const Matrix &uncertainty, GapsRandomState *randState);

-private:

-

-    SparseGibbsSampler mASampler;

-    SparseGibbsSampler mPSampler;