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

 Package: CoGAPS

-Version: 3.3.32

+Version: 3.3.33

 Date: 2018-04-24

 Title: Coordinated Gene Activity in Pattern Sets

 Author: Thomas Sherman, Wai-shing Lee, Conor Kelton, Ondrej Maxian, Jacob Carey,

@@ -76,31 +76,34 @@ void IndexFlagSet::release(uint64_t n)

 AtomPool::AtomPool()

 {

-    mPool = new Atom[POOL_SIZE];

+    //mPool = new Atom[POOL_SIZE];

 }

 AtomPool::~AtomPool()

 {

-    delete[] mPool;

+    //delete[] mPool;

 }

 Atom* AtomPool::alloc()

 {

-    unsigned n = mIndexFlags.getFirstFree();

-    mIndexFlags.set(n);

-    Atom *a = &(mPool[n]);

-    a->poolIndex = n;

-    return a;

+    //unsigned n = mIndexFlags.getFirstFree();

+    //mIndexFlags.set(n);

+    //Atom *a = &(mPool[n]);

+    //a->poolIndex = n;

+    //return a;

+    return new Atom();

 }

 void AtomPool::free(Atom* a)

 {

-    mIndexFlags.release(a->poolIndex);

+    //mIndexFlags.release(a->poolIndex);

+    delete a;

 }

 bool AtomPool::depleted() const

 {

-    return !mIndexFlags.isAnyFree();

+    //return !mIndexFlags.isAnyFree();

+    return false;

 }

 ////////////////////////////// AtomAllocator ///////////////////////////////////

@@ -30,6 +30,7 @@ public:

     AtomicDomain(uint64_t nBins);

+    // TODO can we have internal rng since these are always called sequentially

     // access atoms

     Atom* front();

     Atom* randomAtom(GapsRng *rng);

@@ -242,6 +242,7 @@ void GapsRunner::createCheckpoint()

         // create checkpoint file

         Archive ar(mCheckpointOutFile, ARCHIVE_WRITE);

+        GapsRng::save(ar);

         ar << mNumPatterns << mSeed << mASampler << mPSampler << mStatistics

             << mFixedMatrix << mMaxIterations << mPhase << mCurrentIteration

             << mNumUpdatesA << mNumUpdatesP << mRng;

@@ -84,7 +84,7 @@ template <class DataType>

 GapsRunner::GapsRunner(const DataType &data, bool transposeData,

 unsigned nPatterns, bool partitionRows, const std::vector<unsigned> &indices)

     :

-mASampler(data, !transposeData, nPatterns,!partitionRows, indices),

+mASampler(data, !transposeData, nPatterns, !partitionRows, indices),

 mPSampler(data, transposeData, nPatterns, partitionRows, indices),

 mStatistics(mPSampler.dataRows(), mPSampler.dataCols(), nPatterns),

 mFixedMatrix('N'), mMaxIterations(1000), mMaxThreads(1), mPrintMessages(true),

@@ -120,24 +120,18 @@ void GibbsSampler::processProposal(const AtomicProposal &prop)

 // exponential distribution or with the gibbs mass distribution

 void GibbsSampler::birth(const AtomicProposal &prop)

 {

-    // calculate proposed mass

-    float mass = 0.f;

-    if (canUseGibbs(prop.c1))

-    {

-        mass = gibbsMass(alphaParameters(prop.r1, prop.c1) * mAnnealingTemp,

-            0.f, mMaxGibbsMass, mLambda, &(prop.rng)).value();

-    }

-    else

-    {

-        mass = prop.rng.exponential(mLambda);

-    }

-

-    // accept mass as long as it's non-zero

-    if (mass >= gaps::epsilon)

+    // try to get mass using gibbs, resort to exponential if needed

+    OptionalFloat mass = canUseGibbs(prop.c1)

+        ? gibbsMass(alphaParameters(prop.r1, prop.c1) * mAnnealingTemp, 0.f,

+            mMaxGibbsMass, mLambda, &(prop.rng))

+        : prop.rng.exponential(mLambda);

+

+    // accept mass as long as gibbs succeded or it's non-zero

+    if (mass.hasValue() && mass.value() >= gaps::epsilon)

     {

         mQueue.acceptBirth();

-        prop.atom1->mass = mass;

-        changeMatrix(prop.r1, prop.c1, mass);

+        prop.atom1->mass = mass.value();

+        changeMatrix(prop.r1, prop.c1, mass.value());

     }

     else

     {

@@ -154,7 +148,6 @@ void GibbsSampler::death(const AtomicProposal &prop)

     AlphaParameters alpha = alphaParametersWithChange(prop.r1, prop.c1,

         -prop.atom1->mass);

-    // try to calculate rebirth mass using gibbs distribution, otherwise use original

     float rebirthMass = prop.atom1->mass;

     if (canUseGibbs(prop.c1))

     {

@@ -223,7 +216,7 @@ AlphaParameters alpha)

 {

     // compute amount of mass to be exchanged

     float totalMass = prop.atom1->mass + prop.atom2->mass;

-    float newMass = prop.rng.truncGammaUpper(totalMass, 2.f, 1.f / mLambda);

+    float newMass = prop.rng.truncGammaUpper(totalMass, 1.f / mLambda);

     // compute amount to change atom1 by - always change larger mass to newMass

     float delta = (prop.atom1->mass > prop.atom2->mass)

@@ -266,8 +259,9 @@ void GibbsSampler::acceptExchange(const AtomicProposal &prop, float delta)

 void GibbsSampler::changeMatrix(unsigned row, unsigned col, float delta)

 {

     mMatrix(row, col) += delta;

-    GAPS_ASSERT(mMatrix(row, col) >= 0.f);

     updateAPMatrix(row, col, delta);

+

+    GAPS_ASSERT(mMatrix(row, col) >= 0.f);

 }

 // delta could be negative, this is needed to prevent negative values in matrix

@@ -340,13 +334,17 @@ unsigned col, float ch)

 Archive& operator<<(Archive &ar, GibbsSampler &s)

 {

-    // TODO

+    ar << s.mMatrix << s.mDomain << s.mQueue << s.mAlpha

+        << s.mLambda << s.mMaxGibbsMass << s.mAnnealingTemp << s.mNumPatterns

+        << s.mNumBins << s.mBinLength << s.mDomainLength;

     return ar;

 }

 Archive& operator>>(Archive &ar, GibbsSampler &s)

 {

-    // TODO

+    ar >> s.mMatrix >> s.mDomain >> s.mQueue >> s.mAlpha

+        >> s.mLambda >> s.mMaxGibbsMass >> s.mAnnealingTemp >> s.mNumPatterns

+        >> s.mNumBins >> s.mBinLength >> s.mDomainLength;

     return ar;

 }

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

 #ifndef __COGAPS_GIBBS_SAMPLER_H__

 #define __COGAPS_GIBBS_SAMPLER_H__

+#define DEFAULT_ALPHA           0.01f

+#define DEFAULT_MAX_GIBBS_MASS  100.f

+

 #include "AtomicDomain.h"

 #include "ProposalQueue.h"

 #include "data_structures/Matrix.h"

@@ -112,7 +115,7 @@ mBinLength(std::numeric_limits<uint64_t>::max() / mNumBins),

 mDomainLength(mBinLength * mNumBins)

 {

     // default sparsity parameters

-    setSparsity(0.01, 100.f, false);

+    setSparsity(DEFAULT_ALPHA, DEFAULT_MAX_GIBBS_MASS, false);

 }

 template <class DataType>

@@ -251,14 +251,14 @@ bool ProposalQueue::exchange(AtomicDomain &domain)

     if (prop.r1 == prop.r2 && prop.c1 == prop.c2)

     {

-        float newMass = prop.rng.truncGammaUpper(prop.atom1->mass + prop.atom2->mass, 2.f, 1.f / mLambda);

+        float newMass = prop.rng.truncGammaUpper(prop.atom1->mass + prop.atom2->mass, 1.f / mLambda);

         float delta = (prop.atom1->mass > prop.atom2->mass) ? newMass - prop.atom1->mass : prop.atom2->mass - newMass;

         if (prop.atom1->mass + delta > gaps::epsilon && prop.atom2->mass - delta > gaps::epsilon)

         {

             prop.atom1->mass += delta;

             prop.atom2->mass -= delta;

         }        

-        return true;

+        return true; // automatically accept exchanges in same bin

     }

     mQueue.push_back(prop);

@@ -4,6 +4,9 @@

 TEST_CASE("AtomicDomain")

 {

+    GapsRng::setSeed(123);

+    GapsRng rng;

+

     SECTION("Construction")

     {

         AtomicDomain domain(10);

@@ -47,7 +50,7 @@ TEST_CASE("AtomicDomain")

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

         {

             domain.insert(i, static_cast<float>(i));

-            REQUIRE_NOTHROW(domain.randomFreePosition());

+            REQUIRE_NOTHROW(domain.randomFreePosition(&rng));

         }

     }

@@ -60,10 +63,10 @@ TEST_CASE("AtomicDomain")

             domain.insert(i, static_cast<float>(i));

             // single random atom

-            REQUIRE(domain.randomAtom()->pos < i + 1);

+            REQUIRE(domain.randomAtom(&rng)->pos < i + 1);

             // random atom for exchange

-            AtomNeighborhood hood = domain.randomAtomWithRightNeighbor();

+            AtomNeighborhood hood = domain.randomAtomWithRightNeighbor(&rng);

             REQUIRE(hood.center->pos < i + 1);

             REQUIRE(!hood.hasLeft());

@@ -78,7 +81,7 @@ TEST_CASE("AtomicDomain")

             }

             // random atom for move

-            hood = domain.randomAtomWithNeighbors();

+            hood = domain.randomAtomWithNeighbors(&rng);

             REQUIRE(hood.center->pos < i + 1);

             if (hood.center->pos == 0)

@@ -4,10 +4,11 @@

 #include "../GibbsSampler.h"

 #include "../math/Random.h"

-#if 0

-

 TEST_CASE("Test utils/Archive.h")

 {

+    GapsRng::setSeed(123);

+    GapsRng rng;

+

     SECTION("Reading/Writing to an Archive")

     {

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

@@ -59,13 +60,13 @@ TEST_CASE("Test utils/Archive.h")

         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::uniform(0.0, 2.0);

+            vec_write[i] = rng.uniform(0.f, 2.f);

         }

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

@@ -86,30 +87,24 @@ TEST_CASE("Test utils/Archive.h")

     SECTION("Matrix Serialization")

     {

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

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

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

         {

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

             {

-                rMat_write(i,j) = gaps::random::uniform(0.0, 2.0);

-                cMat_write(i,j) = gaps::random::uniform(0.0, 2.0);

+                cMat_write(i,j) = rng.uniform(0.f, 2.f);

             }

         }

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

-        arWrite << rMat_write;

         arWrite << cMat_write;

         arWrite.close();

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

-        arRead >> rMat_read;

         arRead >> cMat_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());

@@ -117,7 +112,6 @@ TEST_CASE("Test utils/Archive.h")

         {

             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));

             }

         }

@@ -125,9 +119,34 @@ TEST_CASE("Test utils/Archive.h")

     SECTION("GibbsSampler Serialization")

     {

-        //TODO

+        Rcpp::Environment env = Rcpp::Environment::global_env();

+        std::string csvPath = Rcpp::as<std::string>(env["gistCsvPath"]);

+

+        GibbsSampler Asampler(csvPath, false, 7, false, std::vector<unsigned>());

+        GibbsSampler Psampler(csvPath, true, 7, false, std::vector<unsigned>());

+        Asampler.sync(Psampler);

+        Psampler.sync(Asampler);

+        

+        Asampler.update(10000, 1);

+

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

+        arWrite << Asampler;

+        arWrite.close();

+

+        GibbsSampler savedAsampler(csvPath, false, 7, false, std::vector<unsigned>());

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

+        arRead >> savedAsampler;

+        arRead.close();

+    }

+

+#ifdef GAPS_INTERNAL_TESTS    

+    SECTION("AtomicDomain Serialization")

+    {

+

     }

+#endif

+/*

     SECTION("Random Generator Serialization")

     {

         std::vector<float> randSequence;

@@ -163,9 +182,7 @@ TEST_CASE("Test utils/Archive.h")

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

         }

     }

-

+*/

     // cleanup directory

     std::remove("test_ar.temp");

-}

-

-#endif

\ No newline at end of file

+}

\ No newline at end of file

@@ -181,8 +181,8 @@ Archive& operator>>(Archive &ar, ColMatrix &mat)

     // should already by allocated

     unsigned nr = 0, nc = 0;

     ar >> nr >> nc;

-    GAPS_ASSERT(nr == mat.mNumRows);

-    GAPS_ASSERT(nc == mat.mNumCols);

+    GAPS_ASSERT_MSG(nr == mat.mNumRows, nr << " != " << mat.mNumRows);

+    GAPS_ASSERT_MSG(nc == mat.mNumCols, nc << " != " << mat.mNumCols);

     // read in data

     for (unsigned i = 0; i < mat.mNumCols; ++i)

@@ -26,6 +26,7 @@ namespace gaps

     const float epsilon = 1.0e-5f;

     const float pi = 3.1415926535897932384626433832795f;

     const float pi_double = 3.1415926535897932384626433832795;

+    const float sqrt2 = 1.4142135623730950488016887242097f;

     float min(float a, float b);

     unsigned min(unsigned a, unsigned b);

@@ -20,6 +20,41 @@ const double maxU32AsDouble = static_cast<double>(std::numeric_limits<uint32_t>:

 static Xoroshiro128plus seeder;

+////////////////////////////// Lookup Tables ///////////////////////////////////

+

+static float erf_lookup_table[3001];

+static float erfinv_lookup_table[5001];

+static float qgamma_lookup_table[5001];

+

+static void initLookupTables()

+{

+    // erf

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

+    {

+        float x = static_cast<float>(i) / 1000.f;

+        erf_lookup_table[i] = 2.f * gaps::p_norm(x * gaps::sqrt2, 0.f, 1.f) - 1.f;

+    }

+

+    // erfinv

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

+    {

+        float x = static_cast<float>(i) / 5000.f;

+        erfinv_lookup_table[i] = gaps::q_norm((1.f + x) / 2.f, 0.f, 1.f) / gaps::sqrt2;

+    }

+    erfinv_lookup_table[5000] = gaps::q_norm(1.9998f / 2.f, 0.f, 1.f) / gaps::sqrt2;

+

+    // qgamma

+    qgamma_lookup_table[0] = 0.f;

+    for (unsigned i = 1; i < 5000; ++i)

+    {

+        float x = static_cast<float>(i) / 5000.f;

+        qgamma_lookup_table[i] = gaps::q_gamma(x, 2.f, 1.f);

+    }

+    qgamma_lookup_table[5000] = gaps::q_gamma(0.9998f, 2.f, 1.f);

+

+    GAPS_ASSERT(erf_lookup_table[3000] < 1.f);

+}

+

 /////////////////////////////// OptionalFloat //////////////////////////////////

 OptionalFloat::OptionalFloat() : mValue(0.f), mHasValue(false) {}

@@ -94,6 +129,7 @@ Archive& operator>>(Archive &ar, Xoroshiro128plus &gen)

 void GapsRng::setSeed(uint32_t sd)

 {

+    initLookupTables();

     seeder.seed(sd);

 }

@@ -105,6 +141,7 @@ Archive& GapsRng::save(Archive &ar)

 Archive& GapsRng::load(Archive &ar)

 {

+    initLookupTables();

     ar >> seeder;

     return ar;

 }

@@ -204,14 +241,7 @@ uint64_t GapsRng::uniform64(uint64_t a, uint64_t b)

 int GapsRng::poisson(double lambda)

 {

-    if (lambda <= 5.0)

-    {

-        return poissonSmall(lambda);

-    }

-    else

-    {

-        return poissonLarge(lambda);

-    }

+    return lambda <= 5.0 ? poissonSmall(lambda) : poissonLarge(lambda);

 }

 // lambda <= 5

@@ -263,33 +293,29 @@ float GapsRng::exponential(float lambda)

     return -1.f * std::log(uniform()) / lambda;

 }

+// fails if too far in tail

 OptionalFloat GapsRng::truncNormal(float a, float b, float mean, float sd)

 {

-    float pLower = gaps::p_norm(a, mean, sd);

-    float pUpper = gaps::p_norm(b, mean, sd);

+    float pLower = gaps::p_norm_fast(a, mean, sd);

+    float pUpper = gaps::p_norm_fast(b, mean, sd);

     if (!(pLower > 0.95f || pUpper < 0.05f)) // too far in tail

     {

-        float u = uniform(pLower, pUpper);

-        while (u == 0.f || u == 1.f)

-        {

-            u = uniform(pLower, pUpper);

-        }

-        float m = gaps::q_norm(u, mean, sd); 

-        m = gaps::max(a, gaps::min(m, b));

-        return m;

+        GAPS_ASSERT(pLower > 0.f);

+        GAPS_ASSERT(pUpper < 1.f);

+

+        float z = gaps::q_norm_fast(uniform(pLower, pUpper), mean, sd); 

+        z = gaps::max(a, gaps::min(z, b));

+        return z;

     }

     return OptionalFloat();

 }

-float GapsRng::truncGammaUpper(float b, float shape, float scale)

+// shape is hardcoded to 2 since it never changes

+float GapsRng::truncGammaUpper(float b, float scale)

 {

-    float upper = gaps::p_gamma(b, shape, scale);

-    float u = uniform(0.f, upper);

-    while (u == 0.f || u == 1.f)

-    {

-        u = uniform(0.f, upper);

-    }

-    return gaps::q_gamma(u, shape, scale);

+    float upper = 1.f - std::exp(-b / scale) * (1.f + b / scale);

+    unsigned ndx = static_cast<unsigned>(uniform(0.f, upper * 5000.f));

+    return qgamma_lookup_table[ndx] * scale;

 }

 Archive& operator<<(Archive &ar, GapsRng &gen)

@@ -346,13 +372,39 @@ float gaps::p_norm(float p, float mean, float sd)

     return cdf(norm, p);

 }

-double gaps::lgamma(double x)

+float gaps::p_norm_fast(float p, float mean, float sd)

 {

-    return boost::math::lgamma(x);

+    float term = (p - mean) / (sd * gaps::sqrt2);

+    float erf = 0.f;

+    if (term < 0.f)

+    {

+        term = gaps::max(term, -3.f);

+        erf = -erf_lookup_table[static_cast<unsigned>(-term * 1000.f)];

+    }

+    else

+    {

+        term = gaps::min(term, 3.f);

+        erf = erf_lookup_table[static_cast<unsigned>(term * 1000.f)];

+    }

+    return 0.5f * (1.f + erf);

 }

-float gaps::d_norm_fast(float d, float mean, float sd)

+float gaps::q_norm_fast(float q, float mean, float sd)

 {

-    return std::exp((d - mean) * (d - mean) / (-2.f * sd * sd))

-        / std::sqrt(2.f * gaps::pi * sd * sd);

+    float term = 2.f * q - 1.f;

+    float erfinv = 0.f;

+    if (term < 0.f)

+    {

+        erfinv = -erfinv_lookup_table[static_cast<unsigned>(-term * 5000.f)];

+    }

+    else

+    {

+        erfinv = erfinv_lookup_table[static_cast<unsigned>(term * 5000.f)];

+    }

+    return mean + sd * gaps::sqrt2 * erfinv;

+}

+

+double gaps::lgamma(double x)

+{

+    return boost::math::lgamma(x);

 }

\ No newline at end of file

@@ -27,19 +27,20 @@ namespace gaps

 {

     double lgamma(double x);

+    // fast enough with default implementation

     float d_gamma(float d, float shape, float scale);

     float p_gamma(float p, float shape, float scale);

+

+    // standard functions

     float q_gamma(float q, float shape, float scale);

     float d_norm(float d, float mean, float sd);

-    float q_norm(float q, float mean, float sd);

     float p_norm(float p, float mean, float sd);

+    float q_norm(float q, float mean, float sd);

-    float d_gamma_fast(float d, float shape, float scale);

-    float p_gamma_fast(float p, float shape, float scale);

-    float q_gamma_fast(float q, float shape, float scale);

-    float d_norm_fast(float d, float mean, float sd);

+    // fast versions, mostly using lookup tables

     float p_norm_fast(float p, float mean, float sd);

     float q_norm_fast(float q, float mean, float sd);

+

 }

 // used for seeding individual rngs

@@ -83,7 +84,7 @@ public:

     float exponential(float lambda);

     OptionalFloat truncNormal(float a, float b, float mean, float sd);

-    float truncGammaUpper(float b, float shape, float scale);

+    float truncGammaUpper(float b, float scale); // shape hardcoded to 2

     static void setSeed(uint32_t sd);

     static Archive& save(Archive &ar);