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

 Package: CoGAPS

-Version: 3.3.33

+Version: 3.3.50

 Date: 2018-04-24

 Title: Coordinated Gene Activity in Pattern Sets

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

@@ -103,6 +103,8 @@ outputToFile=NULL, ...)

         stop("uncertainty must be a matrix unless data is a file path")

     if (!is(data, "character"))

         checkDataMatrix(data, uncertainty, allParams$gaps)

+    if (!is.null(uncertainty) & allParams$gaps@useSparseOptimization)

+        stop("must use default uncertainty when enabling useSparseOptimization")

     # check single cell parameter

     if (!is.null(allParams$gaps@distributed))

@@ -11,6 +11,8 @@

 #' @slot maxGibbsMassP atomic mass restriction for sample matrix

 #' @slot seed random number generator seed

 #' @slot singleCell is the data single cell?

+#' @slot useSparseOptimization speeds up performance with sparse data, note

+#' this can only be used with the default uncertainty

 #' @slot distributed either "genome-wide" or "single-cell" indicating which

 #' distributed algorithm should be used

 #' @slot nSets [distributed parameter] number of sets to break data into

@@ -35,6 +37,7 @@ setClass("CogapsParams", slots = c(

     maxGibbsMassP = "numeric",

     seed = "numeric",

     singleCell = "logical",

+    useSparseOptimization = "logical",

     distributed = "character_OR_NULL",

     nSets = "numeric",

     cut = "numeric",

@@ -63,6 +66,7 @@ setMethod("initialize", "CogapsParams",

         .Object@maxGibbsMassP <- 100

         .Object@seed <- getMilliseconds(as.POSIXlt(Sys.time()))

         .Object@singleCell <- FALSE

+        .Object@useSparseOptimization <- FALSE

         .Object@distributed <- NULL

         .Object@cut <- .Object@nPatterns

         .Object@nSets <- 4

@@ -72,6 +72,7 @@ const Rcpp::Nullable<Rcpp::IntegerVector> &indices)

     params.maxGibbsMassA = gapsParams.slot("maxGibbsMassA");

     params.maxGibbsMassP = gapsParams.slot("maxGibbsMassP");

     params.singleCell = gapsParams.slot("singleCell");

+    params.useSparseOptimization = gapsParams.slot("useSparseOptimization");

     // check if using fixed matrix

     if (fixedMatrix.isNotNull())

@@ -47,6 +47,11 @@ GapsResult AbstractGapsRunner::run()

     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)

@@ -47,4 +47,5 @@ OBJECTS =   Cogaps.o \

             cpp_tests/testSparseVector.o \

             cpp_tests/testVector.o \

             cpp_tests/testFileParsers.o \

-            cpp_tests/testDenseGibbsSampler.o

\ No newline at end of file

+            cpp_tests/testDenseGibbsSampler.o \

+            cpp_tests/testSparseGibbsSampler.o

\ No newline at end of file

new file mode 100644

@@ -0,0 +1,247 @@

+#include "catch.h"

+#include "../gibbs_sampler/DenseGibbsSampler.h"

+#include "../gibbs_sampler/SparseGibbsSampler.h"

+

+#define TEST_APPROX(x) Approx(x).epsilon(0.001f)

+

+TEST_CASE("Test SparseGibbsSampler")

+{

+    SECTION("Construct from data matrix")

+    {

+        Matrix data(25, 50);

+        for (unsigned i = 0; i < data.nRow(); ++i)

+        {

+            for (unsigned j = 0; j < data.nCol(); ++j)

+            {

+                data(i,j) = i + j + 1.f;

+            }

+        }

+

+        GapsParameters params(data);

+        SparseGibbsSampler ASampler(data, true, false, params.alphaA,

+            params.maxGibbsMassA, params);

+        SparseGibbsSampler PSampler(data, false, false, params.alphaP,

+            params.maxGibbsMassP, params);

+    

+        ASampler.sync(PSampler);

+        PSampler.sync(ASampler);

+

+        REQUIRE(ASampler.chiSq() == 100.f * data.nRow() * data.nCol());

+        REQUIRE(PSampler.chiSq() == 100.f * data.nRow() * data.nCol());

+

+    #ifdef GAPS_DEBUG

+        REQUIRE(ASampler.internallyConsistent());

+        REQUIRE(PSampler.internallyConsistent());

+    #endif

+    }

+

+#ifdef GAPS_INTERNAL_TESTS

+    SECTION("Test consistency between alpha parameters calculations")

+    {

+        // create the "data"

+        Matrix data(100, 75);

+        GapsRng::setSeed(123);

+        GapsRng rng;

+        for (unsigned i = 0; i < data.nRow(); ++i)

+        {

+            for (unsigned j = 0; j < data.nCol(); ++j)

+            {

+                data(i,j) = rng.uniform32(1,14) * (rng.uniform() < 0.5f ? 0.f : 1.f);

+            }

+        }

+

+        // create pair of sparse gibbs samplers

+        GapsParameters params(data);

+        SparseGibbsSampler sparse_ASampler(data, true, false, params.alphaA,

+            params.maxGibbsMassA, params);

+        SparseGibbsSampler sparse_PSampler(data, false, false, params.alphaP,

+            params.maxGibbsMassP, params);

+        sparse_ASampler.sync(sparse_PSampler);

+        sparse_PSampler.sync(sparse_ASampler);

+

+        // create pair of dense gibbs samplers

+        DenseGibbsSampler dense_ASampler(data, true, false, params.alphaA,

+            params.maxGibbsMassA, params);

+        DenseGibbsSampler dense_PSampler(data, false, false, params.alphaP,

+            params.maxGibbsMassP, params);

+        dense_ASampler.sync(dense_PSampler);

+        dense_PSampler.sync(dense_ASampler);

+

+        // set the A and P matrix to the same thing

+        for (unsigned i = 0; i < data.nRow(); ++i)

+        {

+            for (unsigned k = 0; k < params.nPatterns; ++k)

+            {

+                float val = rng.uniform(0.f, 10.f) * (rng.uniform() < 0.2f ? 0.f : 1.f);

+                dense_ASampler.mMatrix(i,k) = val;

+                sparse_ASampler.mMatrix.add(i, k, val);

+            }

+        }        

+        REQUIRE(gaps::sum(dense_ASampler.mMatrix) == gaps::sum(sparse_ASampler.mMatrix));

+

+        for (unsigned j = 0; j < data.nCol(); ++j)

+        {

+            for (unsigned k = 0; k < params.nPatterns; ++k)

+            {

+                float val = rng.uniform(0.f, 10.f) * (rng.uniform() < 0.2f ? 0.f : 1.f);

+                dense_PSampler.mMatrix(j,k) = val;

+                sparse_PSampler.mMatrix.add(j, k, val);

+            }

+        }        

+        REQUIRE(gaps::sum(dense_PSampler.mMatrix) == gaps::sum(sparse_PSampler.mMatrix));

+

+        // sync them back up

+        sparse_ASampler.sync(sparse_PSampler);

+        sparse_PSampler.sync(sparse_ASampler);

+        dense_ASampler.sync(dense_PSampler);

+        dense_PSampler.sync(dense_ASampler);

+        dense_ASampler.recalculateAPMatrix();

+        dense_PSampler.recalculateAPMatrix();

+

+///////////////// test that alphaParameters are the same ///////////////////////

+        for (unsigned i = 0; i < data.nRow(); ++i)

+        {

+            for (unsigned k = 0; k < params.nPatterns; ++k)

+            {

+                AlphaParameters sa = sparse_ASampler.alphaParameters(i,k);

+                AlphaParameters da = dense_ASampler.alphaParameters(i,k);

+                REQUIRE(sa.s >= 0.f);

+                REQUIRE(da.s >= 0.f);

+                if (sa.s <= gaps::epsilon || da.s <= gaps::epsilon)

+                {

+                    REQUIRE(sa.s <= gaps::epsilon);

+                    REQUIRE(da.s <= gaps::epsilon);

+                }

+                REQUIRE(sa.s == TEST_APPROX(da.s));

+                REQUIRE(sa.s_mu == TEST_APPROX(da.s_mu));

+            }

+        }

+

+        for (unsigned j = 0; j < data.nCol(); ++j)

+        {

+            for (unsigned k = 0; k < params.nPatterns; ++k)

+            {

+                AlphaParameters sa = sparse_PSampler.alphaParameters(j,k);

+                AlphaParameters da = dense_PSampler.alphaParameters(j,k);

+                REQUIRE(sa.s >= 0.f);

+                REQUIRE(da.s >= 0.f);

+                if (sa.s <= gaps::epsilon || da.s <= gaps::epsilon)

+                {

+                    REQUIRE(sa.s <= gaps::epsilon);

+                    REQUIRE(da.s <= gaps::epsilon);

+                }

+                REQUIRE(sa.s == TEST_APPROX(da.s));

+                REQUIRE(sa.s_mu == TEST_APPROX(da.s_mu));

+            }

+        }

+

+///////////// test two dimensional alphaParameters are the same ////////////////

+        for (unsigned i = 0; i < data.nRow(); ++i)

+        {

+            for (unsigned k1 = 0; k1 < params.nPatterns; ++k1)

+            {

+                for (unsigned k2 = k1+1; k2 < params.nPatterns; ++k2)

+                {

+                    AlphaParameters sa = sparse_ASampler.alphaParameters(i,k1,i,k2);

+                    AlphaParameters da = dense_ASampler.alphaParameters(i,k1,i,k2);

+                    REQUIRE(sa.s >= 0.f);

+                    REQUIRE(da.s >= 0.f);

+                    if (sa.s <= gaps::epsilon || da.s <= gaps::epsilon)

+                    {

+                        REQUIRE(sa.s <= gaps::epsilon);

+                        REQUIRE(da.s <= gaps::epsilon);

+                    }

+                    REQUIRE(sa.s == TEST_APPROX(da.s));

+                    REQUIRE(sa.s_mu == TEST_APPROX(da.s_mu));

+

+                    // symmetry

+                    sa = sparse_ASampler.alphaParameters(i,k2,i,k1);

+                    da = dense_ASampler.alphaParameters(i,k2,i,k1);

+                    REQUIRE(sa.s >= 0.f);

+                    REQUIRE(da.s >= 0.f);

+                    if (sa.s <= gaps::epsilon || da.s <= gaps::epsilon)

+                    {

+                        REQUIRE(sa.s <= gaps::epsilon);

+                        REQUIRE(da.s <= gaps::epsilon);

+                    }

+                    REQUIRE(sa.s == TEST_APPROX(da.s));

+                    REQUIRE(sa.s_mu == TEST_APPROX(da.s_mu));

+                }

+            }

+        }

+

+        for (unsigned j = 0; j < data.nCol(); ++j)

+        {

+            for (unsigned k1 = 0; k1 < params.nPatterns; ++k1)

+            {

+                for (unsigned k2 = k1+1; k2 < params.nPatterns; ++k2)

+                {

+                    AlphaParameters sa = sparse_PSampler.alphaParameters(j,k1,j,k2);

+                    AlphaParameters da = dense_PSampler.alphaParameters(j,k1,j,k2);

+                    REQUIRE(sa.s >= 0.f);

+                    REQUIRE(da.s >= 0.f);

+                    if (sa.s <= gaps::epsilon || da.s <= gaps::epsilon)

+                    {

+                        REQUIRE(sa.s <= gaps::epsilon);

+                        REQUIRE(da.s <= gaps::epsilon);

+                    }

+                    REQUIRE(sa.s == TEST_APPROX(da.s));

+                    REQUIRE(sa.s_mu == TEST_APPROX(da.s_mu));

+

+                    // symmetry

+                    sa = sparse_PSampler.alphaParameters(j,k2,j,k1);

+                    da = dense_PSampler.alphaParameters(j,k2,j,k1);

+                    REQUIRE(sa.s >= 0.f);

+                    REQUIRE(da.s >= 0.f);

+                    if (sa.s <= gaps::epsilon || da.s <= gaps::epsilon)

+                    {

+                        REQUIRE(sa.s <= gaps::epsilon);

+                        REQUIRE(da.s <= gaps::epsilon);

+                    }

+                    REQUIRE(sa.s == TEST_APPROX(da.s));

+                    REQUIRE(sa.s_mu == TEST_APPROX(da.s_mu));

+                }

+            }

+        }

+

+///////////// test alphaParameters with change are the same ////////////////////

+        for (unsigned i = 0; i < data.nRow(); ++i)

+        {

+            for (unsigned k = 0; k < params.nPatterns; ++k)

+            {

+                float ch = rng.uniform(0.f, 25.f);

+                AlphaParameters sa = sparse_ASampler.alphaParametersWithChange(i,k,ch);

+                AlphaParameters da = dense_ASampler.alphaParametersWithChange(i,k,ch);

+                REQUIRE(sa.s >= 0.f);

+                REQUIRE(da.s >= 0.f);

+                if (sa.s <= gaps::epsilon || da.s <= gaps::epsilon)

+                {

+                    REQUIRE(sa.s <= gaps::epsilon);

+                    REQUIRE(da.s <= gaps::epsilon);

+                }

+                REQUIRE(sa.s == TEST_APPROX(da.s));

+                REQUIRE(sa.s_mu == TEST_APPROX(da.s_mu));

+            }

+        }

+

+        for (unsigned j = 0; j < data.nCol(); ++j)

+        {

+            for (unsigned k = 0; k < params.nPatterns; ++k)

+            {

+                float ch = rng.uniform(0.f, 25.f);

+                AlphaParameters sa = sparse_PSampler.alphaParametersWithChange(j,k,ch);

+                AlphaParameters da = dense_PSampler.alphaParametersWithChange(j,k,ch);

+                REQUIRE(sa.s >= 0.f);

+                REQUIRE(da.s >= 0.f);

+                if (sa.s <= gaps::epsilon || da.s <= gaps::epsilon)

+                {

+                    REQUIRE(sa.s <= gaps::epsilon);

+                    REQUIRE(da.s <= gaps::epsilon);

+                }

+                REQUIRE(sa.s == TEST_APPROX(da.s));

+                REQUIRE(sa.s_mu == TEST_APPROX(da.s_mu));

+            }

+        }

+    }

+#endif

+}

\ No newline at end of file

@@ -7,6 +7,8 @@

 #include "../data_structures/HybridMatrix.h"

 #include "../data_structures/SparseIterator.h"

+#include <bitset>

+

 TEST_CASE("Test SparseIterator.h - One Dimensional")

 {

 #ifdef GAPS_INTERNAL_TESTS

@@ -86,6 +88,112 @@ TEST_CASE("Test SparseIterator.h - Two Dimensional")

         it.next();

         REQUIRE(it.atEnd());

     }

+

+    SECTION("First overlap happens after 64 entries")

+    {

+        SparseVector sv(100);

+        sv.insert(1, 1.f);

+        sv.insert(2, 2.f);

+        sv.insert(3, 3.f);

+        sv.insert(4, 4.f);

+        sv.insert(5, 5.f);

+        sv.insert(74, 74.f);

+        sv.insert(75, 75.f);

+        sv.insert(76, 76.f);

+

+        HybridVector hv(100);

+        hv.add(6, 7.f);

+        hv.add(7, 8.f);

+        hv.add(8, 9.f);

+        hv.add(75, 76.f);

+

+        SparseIteratorTwo it(sv, hv);

+        REQUIRE(it.getValue_1() == 75.f);

+        REQUIRE(it.getValue_2() == 76.f);

+        it.next();

+        REQUIRE(it.atEnd());

+    }

+

+    SECTION("Test Dot Product with gap")

+    {

+        SparseVector sv(300);

+        HybridVector hv(300);

+        Vector dv1(300), dv2(300);

+    

+        // fill vectors

+        GapsRng::setSeed(123);

+        GapsRng rng;

+

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

+        {

+            float val = rng.uniform(50.f,500.f);

+            sv.insert(i, val);

+            dv1[i] = val;

+        }

+

+        for (unsigned i = 32; i < 60; ++i)

+        {

+            float val = rng.uniform(50.f,500.f);

+            hv.add(i, val);

+            dv2[i] = val;

+        }

+

+        for (unsigned i = 70; i < 120; i+=3)

+        {

+            float v1 = rng.uniform(50.f,500.f);

+            sv.insert(i, v1);

+            dv1[i] = v1;   

+

+            float v2 = rng.uniform(50.f,500.f);

+            hv.add(i, v2);

+            dv2[i] = v2;

+        }

+

+        // this part needs to be accounted for

+        for (unsigned i = 128; i < 196; ++i)

+        {

+            float val = rng.uniform(5.f,10.f);

+            sv.insert(i, val);

+            dv1[i] = val;

+        }

+

+        for (unsigned i = 200; i < 300; ++i)

+        {   

+            float v1 = rng.uniform(50.f,500.f);

+            sv.insert(i, v1);

+            dv1[i] = v1;   

+

+            float v2 = rng.uniform(50.f,500.f);

+            hv.add(i, v2);

+            dv2[i] = v2;

+        }

+

+        // calculate dot product

+        float sdot = 0.f, ddot = 0.f;

+        SparseIteratorTwo it(sv, hv);

+        unsigned i = 0;

+        while (!it.atEnd())

+        {

+            while (dv1[i] == 0.f || dv2[i] == 0.f)

+            {

+                ++i;

+            }

+

+            if (i < dv1.size())

+            {

+                ddot += dv1[i] * dv2[i];

+                REQUIRE(dv1[i] == it.getValue_1());

+                REQUIRE(dv2[i] == it.getValue_2());

+                ++i;

+            }

+

+            sdot += it.getValue_1() * it.getValue_2();

+

+            it.next();

+        }

+        REQUIRE(ddot == gaps::dot(dv1, dv2));

+        REQUIRE(sdot == ddot);        

+    }

 #endif

     // could this fail because of SIMD?

@@ -1,14 +1,21 @@

 #include "HybridVector.h"

 #include "../math/Math.h"

+#include "../math/SIMD.h"

+

+#define PAD_SIZE_FOR_SIMD(x) (gaps::simd::Index::increment() * (1 + ((x) - 1) / gaps::simd::Index::increment()))

 HybridVector::HybridVector(unsigned size)

-    : mIndexBitFlags(size / 64 + 1, 0), mData(size, 0.f)

+    :

+mIndexBitFlags(size / 64 + 1, 0),

+mData(PAD_SIZE_FOR_SIMD(size), 0.f),

+mSize(size)

 {}

 HybridVector::HybridVector(const std::vector<float> &v)

     :

 mIndexBitFlags(v.size() / 64 + 1, 0),

-mData(v.size(), 0.f)

+mData(PAD_SIZE_FOR_SIMD(v.size()), 0.f),

+mSize(v.size())

 {

     for (unsigned i = 0; i < v.size(); ++i)

     {

@@ -34,7 +41,7 @@ bool HybridVector::empty() const

 unsigned HybridVector::size() const

 {

-    return mData.size();

+    return mSize;

 }

 bool HybridVector::add(unsigned i, float v)

@@ -65,13 +72,13 @@ const float* HybridVector::densePtr() const

 Archive& operator<<(Archive &ar, HybridVector &vec)

 {

-    ar << vec.mData.size();

+    ar << vec.mSize;

     for (unsigned i = 0; i < vec.mIndexBitFlags.size(); ++i)

     {

         ar << vec.mIndexBitFlags[i];

     }

-    for (unsigned i = 0; i < vec.mData.size(); ++i)

+    for (unsigned i = 0; i < vec.mSize; ++i)

     {

         ar << vec.mData[i];

     }

@@ -82,14 +89,14 @@ Archive& operator>>(Archive &ar, HybridVector &vec)

 {

     unsigned sz = 0;

     ar >> sz;

-    GAPS_ASSERT(sz == vec.mData.size());

+    GAPS_ASSERT(sz == vec.size());

     for (unsigned i = 0; i < vec.mIndexBitFlags.size(); ++i)

     {

         ar >> vec.mIndexBitFlags[i];

     }

-    for (unsigned i = 0; i < vec.mData.size(); ++i)

+    for (unsigned i = 0; i < vec.mSize; ++i)

     {

         ar >> vec.mData[i];

     }

@@ -40,6 +40,7 @@ private:

     std::vector<uint64_t> mIndexBitFlags;

     aligned_vector mData;

+    unsigned mSize;

 };

 #endif // __COGAPS_HYBRID_VECTOR_H__

\ No newline at end of file

@@ -11,7 +11,11 @@ static unsigned countLowerBits(uint64_t u, unsigned pos)

 // clears all bits as low as pos or lower

 static uint64_t clearLowerBits(uint64_t u, unsigned pos)

 {

-    return u & ~((1ull << (pos + 1ull)) - 1ull);    

+    if (pos == 63)

+    {

+        return 0; // TODO understand this bug - should happen automatically

+    }

+    return u & ~((1ull << (pos + 1ull)) - 1ull);

 }

 SparseIterator::SparseIterator(const SparseVector &v)

@@ -49,8 +53,9 @@ mSparseIndex(0),

 mAtEnd(false)

 {   

     GAPS_ASSERT(v1.size() == v2.size());

+

     next();

-    mSparseIndex -= 1; // this gets advanced one too far

+    mSparseIndex -= 1; // next puts us at position 1, this resets to 0

 }

 bool SparseIteratorTwo::atEnd() const

@@ -60,17 +65,23 @@ bool SparseIteratorTwo::atEnd() const

 void SparseIteratorTwo::next()

 {

+    // get the common indices in this chunk

+    mCommon = mFlags_1 & mFlags_2;

+

+    // if nothing common in this chunk, find a chunk that has common indices

     while (!mCommon)

     {

-        // no common values in this chunk, go to next one

-        ++mBigIndex;

-        if (mBigIndex == mTotalIndices)

+        // first count how many sparse indices we are skipping

+        mSparseIndex += __builtin_popcountll(mFlags_1);

+        

+        // advance to next chunk

+        if (++mBigIndex == mTotalIndices)

         {   

             mAtEnd = true;

             return;

         }

-        // check if we have any common values here

+        // update the flags

         mFlags_1 = mSparse.mIndexBitFlags[mBigIndex];

         mFlags_2 = mHybrid.mIndexBitFlags[mBigIndex];

         mCommon = mFlags_1 & mFlags_2;

@@ -83,11 +94,7 @@ void SparseIteratorTwo::next()

     mSparseIndex += 1 + countLowerBits(mFlags_1, mSmallIndex);

     // clear out all skipped indices and the current index from the bitflags

-    // this is needed so that countLowerBits is accurate in the next iteration

     mFlags_1 = clearLowerBits(mFlags_1, mSmallIndex);

-

-    // clear out this bit from common

-    mCommon ^= (1ull << mSmallIndex);

 }

 float SparseIteratorTwo::getValue_1() const

@@ -123,8 +130,9 @@ mAtEnd(false)

 {   

     GAPS_ASSERT(v1.size() == v2.size());

     GAPS_ASSERT(v2.size() == v3.size());

+

     next();

-    mSparseIndex -= 1; // this gets advanced one too far

+    mSparseIndex -= 1;

 }

 bool SparseIteratorThree::atEnd() const

@@ -134,17 +142,23 @@ bool SparseIteratorThree::atEnd() const

 void SparseIteratorThree::next()

 {

+    // get the common indices in this chunk

+    mCommon = mFlags_1 & mFlags_2 & mFlags_3;

+

+    // if nothing common in this chunk, find a chunk that has common indices

     while (!mCommon)

     {

-        // no common values in this chunk, go to next one

-        ++mBigIndex;

-        if (mBigIndex == mTotalIndices)

+        // first count how many sparse indices we are skipping

+        mSparseIndex += __builtin_popcountll(mFlags_1);

+

+        // advance to next chunk

+        if (++mBigIndex == mTotalIndices)

         {   

             mAtEnd = true;

             return;

         }

-        // check if we have any common values here

+        // update the flags

         mFlags_1 = mSparse.mIndexBitFlags[mBigIndex];

         mFlags_2 = mHybrid_1.mIndexBitFlags[mBigIndex];

         mFlags_3 = mHybrid_2.mIndexBitFlags[mBigIndex];

@@ -158,11 +172,7 @@ void SparseIteratorThree::next()

     mSparseIndex += 1 + countLowerBits(mFlags_1, mSmallIndex);

     // clear out all skipped indices and the current index from the bitflags

-    // this is needed so that countLowerBits is accurate in the next iteration

     mFlags_1 = clearLowerBits(mFlags_1, mSmallIndex);

-

-    // clear out this bit from common

-    mCommon ^= (1ull << mSmallIndex);

 }

 float SparseIteratorThree::getValue_1() const

@@ -35,7 +35,9 @@ public:

     float getValue_2() const;

     unsigned getIndex() const;

+#ifndef GAPS_INTERNAL_TESTS

 private:

+#endif

     const SparseVector &mSparse;

     const HybridVector &mHybrid;

@@ -4,13 +4,15 @@

 SparseVector::SparseVector(unsigned size)

     :

 mSize(size),

-mIndexBitFlags(size / 64 + 1, 0)

+mIndexBitFlags(size / 64 + 1, 0),

+mIndexStart(mIndexBitFlags.size(), 0)

 {}

 SparseVector::SparseVector(const std::vector<float> &v)

     :

 mSize(v.size()),

-mIndexBitFlags(v.size() / 64 + 1, 0)

+mIndexBitFlags(v.size() / 64 + 1, 0),

+mIndexStart(mIndexBitFlags.size(), 0)

 {

     for (unsigned i = 0; i < v.size(); ++i)

     {

@@ -18,6 +20,7 @@ mIndexBitFlags(v.size() / 64 + 1, 0)

         {

             mData.push_back(v[i]);

             mIndexBitFlags[i / 64] ^= (1ull << (i % 64));

+            propogate(i / 64);

         }

     }

 }

@@ -25,7 +28,8 @@ mIndexBitFlags(v.size() / 64 + 1, 0)

 SparseVector::SparseVector(const Vector &v)

     :

 mSize(v.size()),

-mIndexBitFlags(v.size() / 64 + 1, 0)

+mIndexBitFlags(v.size() / 64 + 1, 0),

+mIndexStart(mIndexBitFlags.size(), 0)

 {

     for (unsigned i = 0; i < v.size(); ++i)

     {

@@ -33,6 +37,7 @@ mIndexBitFlags(v.size() / 64 + 1, 0)

         {

             mData.push_back(v[i]);

             mIndexBitFlags[i / 64] ^= (1ull << (i % 64));

+            propogate(i / 64);

         }

     }

 }

@@ -64,6 +69,15 @@ void SparseVector::insert(unsigned i, float v)

     mData.insert(mData.begin() + dataIndex, v);

     mIndexBitFlags[i / 64] |= (1ull << (i % 64));

+    propogate(i / 64);

+}

+

+void SparseVector::propogate(unsigned ndx)

+{

+    for (unsigned i = ndx + 1; i < mIndexStart.size(); ++i)

+    {

+        ++mIndexStart[i];

+    }

 }

 Vector SparseVector::getDense() const

@@ -32,11 +32,13 @@ public:

 private:

 #endif

+    unsigned mSize;

     std::vector<uint64_t> mIndexBitFlags;

+    std::vector<unsigned> mIndexStart;

     std::vector<float> mData;

-    unsigned mSize;

     void insert(unsigned i, float v);

+    void propogate(unsigned ndx);

 };

 #endif // __COGAPS_SPARSE_VECTOR_H__

\ No newline at end of file

@@ -1,8 +1,18 @@

 #include "Vector.h"

+#include "../math/SIMD.h"

-Vector::Vector(unsigned size) : mData(size, 0.f) {}

+#define PAD_SIZE_FOR_SIMD(x) (gaps::simd::Index::increment() * (1 + ((x) - 1) / gaps::simd::Index::increment()))

-Vector::Vector(const std::vector<float> &v) : mData(v.size(), 0.f)

+Vector::Vector(unsigned size)

+    :

+mData(PAD_SIZE_FOR_SIMD(size), 0.f),

+mSize(size)

+{}

+

+Vector::Vector(const std::vector<float> &v)

+    :

+mData(PAD_SIZE_FOR_SIMD(v.size()), 0.f),

+mSize(v.size())

 {

     for (unsigned i = 0; i < v.size(); ++i)

     {

@@ -32,12 +42,12 @@ const float* Vector::ptr() const

 unsigned Vector::size() const

 {

-    return mData.size();

+    return mSize;

 }

 void Vector::operator+=(const Vector &v)

 {

-    for (unsigned i = 0; i < mData.size(); ++i)

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

     {

         mData[i] += v[i];

     }

@@ -45,7 +55,7 @@ void Vector::operator+=(const Vector &v)

 Archive& operator<<(Archive &ar, Vector &vec)

 {

-    for (unsigned i = 0; i < vec.mData.size(); ++i)

+    for (unsigned i = 0; i < vec.mSize; ++i)

     {

         ar << vec.mData[i];

     }

@@ -54,7 +64,7 @@ Archive& operator<<(Archive &ar, Vector &vec)

 Archive& operator>>(Archive &ar, Vector &vec)

 {

-    for (unsigned i = 0; i < vec.mData.size(); ++i)

+    for (unsigned i = 0; i < vec.mSize; ++i)

     {

         ar >> vec.mData[i];

     }

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

 private:

     aligned_vector mData;

+    unsigned mSize;

 };

 #endif

\ No newline at end of file

@@ -80,27 +80,19 @@ AlphaParameters DenseGibbsSampler::alphaParameters(unsigned row, unsigned col)

     const float *mat = mOtherMatrix->getCol(col).ptr();

     gaps::simd::PackedFloat pMat, pD, pAP, pS;

-    gaps::simd::PackedFloat partialS(0.f), partialS_mu(0.f);

+    gaps::simd::PackedFloat s(0.f), s_mu(0.f);

     gaps::simd::Index i(0);

-    for (; i <= size - i.increment(); ++i)

+    for (; i < size; ++i)

     {   

         pMat.load(mat + i);

         pD.load(D + i);

         pAP.load(AP + i);

         pS.load(S + i);

         gaps::simd::PackedFloat ratio(pMat / (pS * pS));

-        partialS += pMat * ratio;

-        partialS_mu += ratio * (pD - pAP);

+        s += pMat * ratio;

+        s_mu += ratio * (pD - pAP);

     }

-

-    float s = partialS.scalar(), s_mu = partialS_mu.scalar();

-    for (unsigned j = i.value(); j < size; ++j)

-    {

-        float ratio = mat[j] / (S[j] * S[j]);

-        s += mat[j] * ratio;

-        s_mu += ratio * (D[j] - AP[j]);

-    }

-    return AlphaParameters(s, s_mu);

+    return AlphaParameters(s.scalar(), s_mu.scalar());

 }

 // PERFORMANCE_CRITICAL

@@ -117,28 +109,20 @@ unsigned r2, unsigned c2)

         const float *mat2 = mOtherMatrix->getCol(c2).ptr();

         gaps::simd::PackedFloat pMat1, pMat2, pD, pAP, pS;

-        gaps::simd::PackedFloat partialS(0.f), partialS_mu(0.f);

+        gaps::simd::PackedFloat s(0.f), s_mu(0.f);

         gaps::simd::Index i(0);

-        for (; i <= size - i.increment(); ++i)

-        {   

+        for (; i < size; ++i)

+        {

             pMat1.load(mat1 + i);

             pMat2.load(mat2 + i);

             pD.load(D + i);

             pAP.load(AP + i);

             pS.load(S + i);

             gaps::simd::PackedFloat ratio((pMat1 - pMat2) / (pS * pS));

-            partialS += (pMat1 - pMat2) * ratio;

-            partialS_mu += ratio * (pD - pAP);

+            s += (pMat1 - pMat2) * ratio;

+            s_mu += ratio * (pD - pAP);

         }

-

-        float s = partialS.scalar(), s_mu = partialS_mu.scalar();

-        for (unsigned j = i.value(); j < size; ++j)

-        {

-            float ratio = (mat1[j] - mat2[j]) / (S[j] * S[j]);

-            s += (mat1[j] - mat2[j]) * ratio;

-            s_mu += ratio * (D[j] - AP[j]);

-        }

-        return AlphaParameters(s, s_mu);

+        return AlphaParameters(s.scalar(), s_mu.scalar());

     }

     return alphaParameters(r1, c1) + alphaParameters(r2, c2);

 }

@@ -155,28 +139,19 @@ unsigned col, float ch)

     gaps::simd::PackedFloat pCh(ch);

     gaps::simd::PackedFloat pMat, pD, pAP, pS;

-    gaps::simd::PackedFloat partialS(0.f), partialS_mu(0.f);

+    gaps::simd::PackedFloat s(0.f), s_mu(0.f);

     gaps::simd::Index i(0);

-    for (; i <= size - i.increment(); ++i)

+    for (; i < size; ++i)

     {   

         pMat.load(mat + i);

         pD.load(D + i);

         pAP.load(AP + i);

         pS.load(S + i);

         gaps::simd::PackedFloat ratio(pMat / (pS * pS));

-        partialS += pMat * ratio;

-        partialS_mu += ratio * (pD - (pAP + pCh * pMat));

-    }