@@ -31,4 +31,5 @@ update_remotes.sh

 inst/profiling/callgrind.out.*

 inst/profiling/vgcore.*

 inst/profiling/valgrind_out.txt

->>>>>>> develop

+

+*.temp

@@ -66,7 +66,6 @@ gapsRun <- function(D, S, ABins = data.frame(), PBins = data.frame(),

                     alphaP = 0.01, nMaxP = 100000, max_gibbmass_paraP = 100.0,

                     seed=-1, messages=TRUE, singleCellRNASeq=FALSE,

                     fixedPatterns = matrix(0), whichMatrixFixed = 'N',

-                    checkpoint_file_name = "gaps_checkpoint.out",

                     checkpoint_interval = 0)

 {

     # Floor the parameters that are integers to prevent allowing doubles.

@@ -13,8 +13,7 @@ gapsRun(D, S, ABins = data.frame(), PBins = data.frame(), nFactor = 7,

   max_gibbmass_paraA = 100, alphaP = 0.01, nMaxP = 1e+05,

   max_gibbmass_paraP = 100, seed = -1, messages = TRUE,

   singleCellRNASeq = FALSE, fixedPatterns = matrix(0),

-  whichMatrixFixed = "N", checkpoint_file_name = "gaps_checkpoint.out",

-  checkpoint_interval = 0)

+  whichMatrixFixed = "N", checkpoint_interval = 0)

 }

 \arguments{

 \item{D}{data matrix}

@@ -71,9 +70,9 @@ domain for relative probabilities}

 \item{whichMatrixFixed}{character to indicate whether A or P matrix

 contains the fixed patterns}

-\item{checkpoint_file_name}{name of file to store checkpoint}

-

 \item{checkpoint_interval}{time (in seconds) between cogaps checkpoints}

+

+\item{checkpoint_file_name}{name of file to store checkpoint}

 }

 \description{

 \code{gapsRun} calls the C++ MCMC code and performs Bayesian

@@ -98,13 +98,12 @@ public:

     double totalMass() const {return mTotalMass;}

     uint64_t numAtoms() const {return mNumAtoms;}

     double at(uint64_t loc) const {return mAtomicDomain.at(loc);}

+    double test_at(uint64_t loc) const {return mAtomicDomain.count(loc) ? mAtomicDomain.at(loc) : 0.0;}

     // setters

     void setAlpha(double alpha) {mAlpha = alpha;}

     void setLambda(double lambda) {mLambda = lambda;}

-    void serializeAndWrite(const std::ofstream &file);

-

     friend void operator<<(Archive &ar, AtomicSupport &sampler);

     friend void operator>>(Archive &ar, AtomicSupport &sampler);

 };

@@ -10,6 +10,11 @@

 #include <boost/archive/text_oarchive.hpp>

 #include <boost/archive/text_iarchive.hpp>

+// no C++11 std::to_string

+#include <sstream>

+#define SSTR( x ) static_cast< std::ostringstream & >( \

+        ( std::ostringstream() << std::dec << x ) ).str()

+

 #define ARCHIVE_MAGIC_NUM 0xCE45D32A

 namespace bpt = boost::posix_time;

@@ -18,9 +23,11 @@ static bpt::ptime lastCheckpoint;

 static void createCheckpoint(GapsInternalState &state)

 {

+    state.numCheckpoints++;

     std::cout << "creating gaps checkpoint...";

     bpt::ptime start = bpt::microsec_clock::local_time();

-    Archive ar("gaps_checkpoint.out", ARCHIVE_WRITE);

+    Archive ar("gaps_checkpoint_" + SSTR(state.numCheckpoints) + ".out",

+        ARCHIVE_WRITE);

     ar << ARCHIVE_MAGIC_NUM;

     gaps::random::save(ar);

     ar << state.nEquil;

@@ -41,6 +41,7 @@ struct GapsInternalState

     uint32_t seed;

     long checkpointInterval;

+    unsigned numCheckpoints;

     GibbsSampler sampler;

@@ -60,7 +61,7 @@ struct GapsInternalState

         nIterA(10), nIterP(10), nEquil(nE), nEquilCool(nEC), nSample(nS),

         nSnapshots(numSnapshots), nOutputs(numOutputs), messages(msg),

         iter(0), phase(GAPS_CALIBRATION), seed(in_seed),

-        checkpointInterval(cptInterval),

+        checkpointInterval(cptInterval), numCheckpoints(0),

         sampler(DMatrix, SMatrix, nFactor, alphaA, alphaP,

             maxGibbsMassA, maxGibbsMassP, singleCellRNASeq, fixedPatterns,

             whichMatrixFixed)

@@ -75,7 +76,7 @@ struct GapsInternalState

     {}

 };

-void operator<<(Archive &ar, GapsInternalState &state)

+inline void operator<<(Archive &ar, GapsInternalState &state)

 {

     ar << state.chi2VecEquil;

     ar << state.nAtomsAEquil;

@@ -95,10 +96,11 @@ void operator<<(Archive &ar, GapsInternalState &state)

     ar << state.phase;

     ar << state.seed;

     ar << state.checkpointInterval;

+    ar << state.numCheckpoints;

     ar << state.sampler;

 }

-void operator>>(Archive &ar, GapsInternalState &state)

+inline void operator>>(Archive &ar, GapsInternalState &state)

 {

     ar >> state.chi2VecEquil;

     ar >> state.nAtomsAEquil;

@@ -117,7 +119,8 @@ void operator>>(Archive &ar, GapsInternalState &state)

     ar >> state.iter;

     ar >> state.phase;

     ar >> state.seed;

-    ar << state.checkpointInterval;

+    ar >> state.checkpointInterval;

+    ar >> state.numCheckpoints;

     ar >> state.sampler;

 }

@@ -11,5 +11,6 @@ OBJECTS =   Algorithms.o \

             cpp_tests/testAtomicSupport.o \

             cpp_tests/testGibbsSampler.o \

             cpp_tests/testMatrix.o \

-            cpp_tests/testRandom.o

+            cpp_tests/testRandom.o \

+            cpp_tests/testSerialization.o

@@ -126,14 +126,10 @@ TEST_CASE("Internal AtomicSupport Tests")

     {

         for (unsigned i = 0; i < 10000; ++i)

         {

-            REQUIRE(Adomain.getRow(gaps::random::uniform64()) >= 0);

             REQUIRE(Adomain.getRow(gaps::random::uniform64()) < nrow);

-            REQUIRE(Adomain.getCol(gaps::random::uniform64()) >= 0);

             REQUIRE(Adomain.getCol(gaps::random::uniform64()) < ncol);

-            REQUIRE(Pdomain.getRow(gaps::random::uniform64()) >= 0);

             REQUIRE(Pdomain.getRow(gaps::random::uniform64()) < nrow);

-            REQUIRE(Pdomain.getCol(gaps::random::uniform64()) >= 0);

             REQUIRE(Pdomain.getCol(gaps::random::uniform64()) < ncol);            

         }

     }

@@ -30,9 +30,6 @@ TEST_CASE("Test Random.h - Random Number Generation")

     SECTION("Test uniform distribution over general interval")

     {

-        // invalid bounds

-        REQUIRE_THROWS(gaps::random::uniform(2.0, 1.4));

-

         // bounds equal

         REQUIRE(gaps::random::uniform(4.3,4.3) == 4.3);

@@ -0,0 +1,219 @@

+#include "catch.h"

+#include "../Archive.h"

+#include "../Matrix.h"

+#include "../AtomicSupport.h"

+#include "../GibbsSampler.h"

+#include "../InternalState.h"

+#include "../Random.h"

+

+TEST_CASE("Test Archive.h")

+{

+    SECTION("Reading/Writing to an Archive")

+    {

+        Archive ar1("test_ar.temp", ARCHIVE_WRITE);

+        ar1 << 3;

+        ar1.close();

+

+        Archive ar2("test_ar.temp", ARCHIVE_READ);

+        unsigned i = 0;

+        ar2 >> i;

+        REQUIRE(i == 3);

+        ar2.close();

+    }

+

+    SECTION("Serialization of primitive types")

+    {

+        // test values

+        unsigned u_read = 0, u_write = 456;

+        uint32_t u32_read = 0, u32_write = 512;

+        uint64_t u64_read = 0, u64_write = 0xAABBCCDDEE;

+        float f_read = 0.f, f_write = 0.123542f;

+        double d_read = 0., d_write = 0.54362;

+        bool b_read = false, b_write = true;

+

+        // write to archive

+        Archive arWrite("test_ar.temp", ARCHIVE_WRITE);

+        arWrite << u_write;

+        arWrite << u32_write;

+        arWrite << u64_write;

+        arWrite << f_write;

+        arWrite << d_write;

+        arWrite << b_write;

+        arWrite.close();

+

+        // read from archive

+        Archive arRead("test_ar.temp", ARCHIVE_READ);

+        arRead >> u_read;

+        arRead >> u32_read;

+        arRead >> u64_read;

+        arRead >> f_read;

+        arRead >> d_read;

+        arRead >> b_read;

+        arRead.close();

+

+        // test that values are the same

+        REQUIRE(u_read == u_write);

+        REQUIRE(u32_read == u32_write);

+        REQUIRE(u64_read == u64_write);

+        REQUIRE(f_read == f_write);

+        REQUIRE(d_read == d_write);

+        REQUIRE(b_read == b_write);

+    }

+    

+    SECTION("Vector Serialization")

+    {

+        Vector vec_read(100), vec_write(100);

+        for (unsigned i = 0; i < 100; ++i)

+        {

+            vec_write[i] = gaps::random::normal(0.0, 2.0);

+        }

+

+        Archive arWrite("test_ar.temp", ARCHIVE_WRITE);

+        arWrite << vec_write;

+        arWrite.close();

+

+        Archive arRead("test_ar.temp", ARCHIVE_READ);

+        arRead >> vec_read;

+        arRead.close();

+

+        REQUIRE(vec_read.size() == vec_write.size());

+

+        for (unsigned i = 0; i < 100; ++i)

+        {

+            REQUIRE(vec_read[i] == vec_write[i]);

+        }

+    }

+

+    SECTION("Matrix Serialization")

+    {

+        RowMatrix rMat_read(100,100), rMat_write(100,100);

+        ColMatrix cMat_read(100,100), cMat_write(100,100);

+        TwoWayMatrix twMat_read(100,100), twMat_write(100,100);

+

+        for (unsigned i = 0; i < 100; ++i)

+        {

+            for (unsigned j = 0; j < 100; ++j)

+            {

+                rMat_write(i,j) = gaps::random::normal(0.0, 2.0);

+                cMat_write(i,j) = gaps::random::normal(0.0, 2.0);

+                twMat_write.set(i,j,gaps::random::normal(0.0, 2.0));

+            }

+        }

+

+        Archive arWrite("test_ar.temp", ARCHIVE_WRITE);

+        arWrite << rMat_write;

+        arWrite << cMat_write;

+        arWrite << twMat_write;

+        arWrite.close();

+

+        Archive arRead("test_ar.temp", ARCHIVE_READ);

+        arRead >> rMat_read;

+        arRead >> cMat_read;

+        arRead >> twMat_read;

+        arRead.close();

+

+        REQUIRE(rMat_read.nRow() == rMat_write.nRow());

+        REQUIRE(rMat_read.nCol() == rMat_write.nCol());

+        REQUIRE(cMat_read.nRow() == cMat_write.nRow());

+        REQUIRE(cMat_read.nCol() == cMat_write.nCol());

+        REQUIRE(twMat_read.nRow() == twMat_write.nRow());

+        REQUIRE(twMat_read.nCol() == twMat_write.nCol());

+    

+        for (unsigned i = 0; i < 100; ++i)

+        {

+            for (unsigned j = 0; j < 100; ++j)

+            {

+                REQUIRE(rMat_read(i,j) == rMat_write(i,j));

+                REQUIRE(cMat_read(i,j) == cMat_write(i,j));

+                REQUIRE(twMat_read.getRow(i)[j] == twMat_write.getRow(i)[j]);

+            }

+        }

+    }

+    

+    SECTION("Atomic Serialization")

+    {

+        AtomicSupport domain_read('A',100,100), domain_write('A',100,100);

+        std::vector<uint64_t> locations;

+        for (unsigned i = 0; i < 10000; ++i)

+        {

+            AtomicProposal prop = domain_write.makeProposal();

+            locations.push_back(prop.pos1);

+            locations.push_back(prop.pos2);

+            domain_write.acceptProposal(prop);

+        }

+

+        Archive arWrite("test_ar.temp", ARCHIVE_WRITE);

+        arWrite << domain_write;

+        arWrite.close();

+    

+        Archive arRead("test_ar.temp", ARCHIVE_READ);

+        arRead >> domain_read;

+        arRead.close();

+

+        REQUIRE(domain_read.alpha() == domain_write.alpha());

+        REQUIRE(domain_read.lambda() == domain_write.lambda());

+        REQUIRE(domain_read.totalMass() == domain_write.totalMass());

+        REQUIRE(domain_read.numAtoms() == domain_write.numAtoms());

+

+#ifdef GAPS_INTERNAL_TESTS

+        std::map<uint64_t, double>::iterator readIter, writeIter;

+        readIter = domain_read.mAtomicDomain.begin();

+        writeIter = domain_write.mAtomicDomain.begin();        

+        while (readIter != domain_read.mAtomicDomain.end())

+        {

+            REQUIRE(readIter->first == writeIter->first);   

+            REQUIRE(readIter->second == writeIter->second);

+            ++readIter;

+            ++writeIter;

+        }

+#endif

+    }

+

+    SECTION("GibbsSampler Serialization")

+    {

+

+    }

+

+    SECTION("GapsInternalState Serialization")

+    {

+

+    }

+

+    SECTION("Random Generator Serialization")

+    {

+        std::vector<double> randSequence;

+

+        gaps::random::setSeed(123);

+        volatile double burn_in = 0.0;

+        for (unsigned i = 0; i < 1000; ++i)

+        {

+            burn_in = gaps::random::uniform(0,1);

+        }

+        REQUIRE(burn_in < 1);

+

+        Archive arWrite("test_ar.temp", ARCHIVE_WRITE);

+        gaps::random::save(arWrite);

+        arWrite.close();

+

+        for (unsigned i = 0; i < 1000; ++i)

+        {

+            randSequence.push_back(gaps::random::uniform());

+            randSequence.push_back(gaps::random::uniform(0.3, 6.4));

+            randSequence.push_back(gaps::random::normal(0.0, 2.0));

+            randSequence.push_back(gaps::random::exponential(5.5));

+        }

+    

+        gaps::random::setSeed(11);

+        Archive arRead("test_ar.temp", ARCHIVE_READ);

+        gaps::random::load(arRead);

+        arRead.close();

+

+        for (unsigned i = 0; i < 1000; ++i)

+        {

+            REQUIRE(gaps::random::uniform() == randSequence[i++]);

+            REQUIRE(gaps::random::uniform(0.3, 6.4) == randSequence[i++]);

+            REQUIRE(gaps::random::normal(0.0, 2.0) == randSequence[i++]);

+            REQUIRE(gaps::random::exponential(5.5) == randSequence[i]);

+        }

+    }

+}

\ No newline at end of file