@@ -63,7 +63,7 @@ std::string getBuildReport_cpp()

     std::string simd = "SIMD not enabled\n";

 #endif

-#if defined(_GAPS_OPENMP)

+#if defined(__GAPS_OPENMP__)

     std::string openmp = "Compiled with OpenMP\n";

 #else

     std::string openmp = "Compiler did not support OpenMP\n";

@@ -207,11 +207,6 @@ void GapsRunner::displayStatus(const std::string &type, unsigned nIterTotal)

 {

     if (mNumOutputs > 0 && (mCurrentIter + 1) % mNumOutputs == 0 && mPrintMessages)

     {

-        // print algo status

-        Rprintf("%s %d of %d, Atoms:%lu(%lu) Chi2 = %.2f\n", type.c_str(),

-            mCurrentIter + 1, nIterTotal, mASampler.nAtoms(),

-            mPSampler.nAtoms(), mASampler.chi2());

-

         // compute time

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

         double elapsed = diff.total_seconds();

@@ -221,6 +216,11 @@ void GapsRunner::displayStatus(const std::string &type, unsigned nIterTotal)

         printTime("Elapsed Time", elapsed);

         printTime("Estimated Total Time", estTotal);

         printTime("Estimated Remaining Time", estTotal - elapsed);

+

+        // print algo status

+        Rprintf("%s %d of %d, Atoms:%lu(%lu) Chi2 = %.2f\n\n", type.c_str(),

+            mCurrentIter + 1, nIterTotal, mASampler.nAtoms(),

+            mPSampler.nAtoms(), mASampler.chi2());

     }

 }

@@ -3,6 +3,7 @@

 #include "Archive.h"

 #include "data_structures/Matrix.h"

+#include "data_structures/Vector.h"

 #include "GibbsSampler.h"

 #include "GapsStatistics.h"

@@ -11,18 +11,16 @@ void GapsStatistics::update(const AmplitudeGibbsSampler &ASampler,

 const PatternGibbsSampler &PSampler)

 {

     mStatUpdates++;

-

-    Vector normVec(mNumPatterns);

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

     {

-        normVec[j] = gaps::algo::sum(PSampler.mMatrix.getRow(j));

-        normVec[j] = normVec[j] == 0 ? 1.f : normVec[j];

+        float norm = gaps::algo::sum(PSampler.mMatrix.getRow(j));

+        norm = norm == 0.f ? 1.f : norm;

-        Vector quot(PSampler.mMatrix.getRow(j) / normVec[j]);

+        Vector quot(PSampler.mMatrix.getRow(j) / norm);

         mPMeanMatrix.getRow(j) += quot;

         mPStdMatrix.getRow(j) += gaps::algo::elementSq(quot);

-        Vector prod(ASampler.mMatrix.getCol(j) * normVec[j]);

+        Vector prod(ASampler.mMatrix.getCol(j) * norm);

         mAMeanMatrix.getCol(j) += prod;

         mAStdMatrix.getCol(j) += gaps::algo::elementSq(prod); 

     }

@@ -9,6 +9,7 @@ OBJECTS =   AtomicDomain.o \

             test-runner.o \

             math/Algorithms.o \

             data_structures/Matrix.o \

+            data_structures/Vector.o \

             math/Random.o \

             cpp_tests/testAlgorithms.o \

             cpp_tests/testAtomicDomain.o \

@@ -48,27 +48,5 @@ TEST_CASE("Test Matrix.h")

         REQUIRE(cmS.nRow() == 1363);

         REQUIRE(cmS.nCol() == 9);

     }

-

-    SECTION("Matrix Update")

-    {

-        RowMatrix rm(10, 25);

-        ColMatrix cm(10, 25);

-

-        rm.update(MatrixChange('A', 3, 4, 3.0));

-        cm.update(MatrixChange('P', 1, 2, 13.0));

-

-        REQUIRE(rm(3,4) == 3.0);

-        REQUIRE(cm(1,2) == 13.0);

-        REQUIRE(rm.getRow(3)[4] == 3.0);

-        REQUIRE(cm.getCol(2)[1] == 13.0);

-

-        rm.update(MatrixChange('A', 3, 4, 3.0));

-        cm.update(MatrixChange('P', 1, 2, 5.0));

-        

-        REQUIRE(rm(3,4) == 6.0);

-        REQUIRE(cm(1,2) == 18.0);

-        REQUIRE(rm.getRow(3)[4] == 6.0);

-        REQUIRE(cm.getCol(2)[1] == 18.0);

-    }

 }

@@ -20,7 +20,7 @@ public:

     void insert(uint64_t n) {mSet[n] = mCurrentKey;}

 };

-// have sorted vector with at least some % of holes

+// TODO have sorted vector with at least some % of holes

 // even distribute entries along it

 // when shift happens, should be minimal

@@ -37,7 +37,7 @@ public:

     void insert(uint64_t p) {mVec.push_back(p);}

     void clear() {mVec.clear();}

-    // inclusive of a and b

+    // inclusive of a and b, TODO improve performance

     bool isEmptyInterval(uint64_t a, uint64_t b)

     {

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

@@ -1,26 +1,7 @@

 #include "Matrix.h"

-#include <stdexcept>

-

-static const float EPSILON = 1.e-10;

-

 /******************************** HELPER *******************************/

-static void updateHelper2(float& val, float delta)

-{

-    val = std::abs(val + delta) < EPSILON ? 0.0 : val + delta;

-}

-

-template<class GenericMatrix>

-static void updateHelper(GenericMatrix &mat, const MatrixChange &change)

-{

-    updateHelper2(mat(change.row1, change.col1), change.delta1);

-    if (change.nChanges > 1)

-    {

-        updateHelper2(mat(change.row2, change.col2), change.delta2);

-    }

-}

-

 template<class GenericMatrix>

 static Rcpp::NumericMatrix convertToRMatrix(const GenericMatrix &mat)

 {

@@ -35,118 +16,16 @@ static Rcpp::NumericMatrix convertToRMatrix(const GenericMatrix &mat)

     return rmat;

 }

-/******************************** VECTOR *******************************/

-

-void Vector::concat(const Vector& vec)

-{

-    mValues.insert(mValues.end(), vec.mValues.begin(), vec.mValues.end());

-}

-

-void Vector::operator+=(const Vector &vec)

-{

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

-    {

-        mValues[i] += vec[i];

-    }

-}

-

-void Vector::operator=(const Vector &vec)

-{

-    mValues = vec.mValues;

-}

-

-Vector Vector::operator+(Vector vec) const

-{

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

-    {

-        vec[i] += mValues[i];

-    }

-    return vec;

-}

-

-Vector Vector::operator-(Vector vec) const

-{

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

-    {

-        vec[i] -= mValues[i];

-    }

-    return vec;

-}

-

-Vector Vector::operator*(Vector vec) const

-{

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

-    {

-        vec[i] *= mValues[i];

-    }

-    return vec;

-}

-

-Vector Vector::operator/(Vector vec) const

-{

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

-    {

-        vec[i] /= mValues[i];

-    }

-    return vec;

-}

-

-Vector Vector::operator+(float val) const

-{

-    Vector vec(*this);

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

-    {

-        vec[i] += val;

-    }

-    return vec;

-}

-

-Vector Vector::operator-(float val) const

-{

-    Vector vec(*this);

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

-    {

-        vec[i] -= val;

-    }

-    return vec;

-}

-

-Vector Vector::operator*(float val) const

-{

-    Vector vec(*this);

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

-    {

-        vec[i] *= val;

-    }

-    return vec;

-}

-

-Vector Vector::operator/(float val) const

-{

-    Vector vec(*this);

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

-    {

-        vec[i] /= val;

-    }

-    return vec;

-}

-

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

+template<class MatA, class MatB>

+static void copyMatrix(MatA &dest, const MatB &source)

 {

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

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

     {

-        ar << vec[i];

-    }

-    return ar;

-}

-

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

-{

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

-    {

-        ar >> vec.mValues[i];

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

+        {

+            dest(i,j) = source(i,j);

+        }

     }

-    return ar;

 }

 /****************************** ROW MATRIX *****************************/

@@ -173,36 +52,19 @@ RowMatrix::RowMatrix(const Rcpp::NumericMatrix &rmat)

     }

 }

-RowMatrix::RowMatrix(const RowMatrix &mat)

-: mNumRows(mat.nRow()), mNumCols(mat.nCol())

+void RowMatrix::operator=(const RowMatrix &mat)

 {

-    for (unsigned i = 0; i < mNumRows; ++i)

-    {

-        mRows.push_back(mat.getRow(i));

-    }

-}

-

-RowMatrix::RowMatrix(const ColMatrix &mat)

-: mNumRows(mat.nRow()), mNumCols(mat.nCol())

-{

-    for (unsigned i = 0; i < mNumRows; ++i)

-    {

-        mRows.push_back(Vector(mNumCols));

-        for (unsigned j = 0; j < mNumCols; ++j)

-        {

-            this->operator()(i,j) = mat(i,j);

-        }

-    }

+    copyMatrix(*this, mat);

 }

-void RowMatrix::update(const MatrixChange &change)

+void RowMatrix::operator=(const ColMatrix &mat)

 {

-    updateHelper<RowMatrix>(*this, change);

+    copyMatrix(*this, mat);

 }

 Rcpp::NumericMatrix RowMatrix::rMatrix() const

 {

-    return convertToRMatrix<RowMatrix>(*this);

+    return convertToRMatrix(*this);

 }

 RowMatrix RowMatrix::operator/(float val) const

@@ -215,17 +77,6 @@ RowMatrix RowMatrix::operator/(float val) const

     return mat;

 }

-void RowMatrix::operator=(const ColMatrix &mat)

-{

-    for (unsigned i = 0; i < mNumRows; ++i)

-    {

-        for (unsigned j = 0; j < mNumCols; ++j)

-        {

-            this->operator()(i,j) = mat(i,j);

-        }

-    }

-}

-

 Archive& operator<<(Archive &ar, RowMatrix &mat)

 {

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

@@ -268,20 +119,6 @@ ColMatrix::ColMatrix(const Rcpp::NumericMatrix &rmat)

     }

 }

-ColMatrix::ColMatrix(const ColMatrix &mat)

-: mNumRows(mat.nRow()), mNumCols(mat.nCol())

-{

-    for (unsigned j = 0; j < mNumCols; ++j)

-    {

-        mCols.push_back(mat.getCol(j));

-    }

-}

-

-void ColMatrix::update(const MatrixChange &change)

-{

-    updateHelper<ColMatrix>(*this, change);

-}

-

 ColMatrix ColMatrix::operator/(float val) const

 {

     ColMatrix mat(mNumRows, mNumCols);

@@ -292,20 +129,19 @@ ColMatrix ColMatrix::operator/(float val) const

     return mat;

 }

+void ColMatrix::operator=(const ColMatrix &mat)

+{

+    copyMatrix(*this, mat);

+}

+

 void ColMatrix::operator=(const RowMatrix &mat)

 {

-    for (unsigned i = 0; i < mNumRows; ++i)

-    {

-        for (unsigned j = 0; j < mNumCols; ++j)

-        {

-            this->operator()(i,j) = mat(i,j);

-        }

-    }

+    copyMatrix(*this, mat);

 }

 Rcpp::NumericMatrix ColMatrix::rMatrix() const

 {

-    return convertToRMatrix<ColMatrix>(*this);

+    return convertToRMatrix(*this);

 }

 Archive& operator<<(Archive &ar, ColMatrix &mat)

@@ -2,9 +2,9 @@

 #define __COGAPS_MATRIX_H__

 #include "../Archive.h"

+#include "Vector.h"

 #include <Rcpp.h>

-#include <boost/align/aligned_allocator.hpp>

 #include <vector>

 struct MatrixChange

@@ -32,46 +32,6 @@ struct MatrixChange

     {}

 };

-namespace bal = boost::alignment;

-typedef std::vector<float, bal::aligned_allocator<float,32> > aligned_vector;

-

-class Vector

-{

-private:

-

-    aligned_vector mValues;

-

-public:

-

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

-    Vector(unsigned size, float val) : mValues(aligned_vector(size, val)) {}

-    Vector(const aligned_vector &v) : mValues(v) {}

-    Vector(const Vector &vec) : mValues(vec.mValues) {}

-

-    const float* ptr() const {return &mValues[0];}

-    float* ptr() {return &mValues[0];}

-

-    float& operator[](unsigned i) {return mValues[i];}

-    float operator[](unsigned i) const {return mValues[i];}

-    unsigned size() const {return mValues.size();}

-

-    Rcpp::NumericVector rVec() const {return Rcpp::wrap(mValues);}

-    void concat(const Vector& vec);

-    void operator+=(const Vector &vec);

-    void operator=(const Vector &vec);

-

-    Vector operator+(Vector vec) const;

-    Vector operator-(Vector vec) const;

-    Vector operator*(Vector vec) const;

-    Vector operator/(Vector vec) const;

-    Vector operator+(float val) const;

-    Vector operator-(float val) const;

-    Vector operator*(float val) const;

-    Vector operator/(float val) const;

-

-    friend Archive& operator<<(Archive &ar, Vector &vec);

-    friend Archive& operator>>(Archive &ar, Vector &vec);

-};

 class ColMatrix;

@@ -82,14 +42,11 @@ private:

     std::vector<Vector> mRows;

     unsigned mNumRows, mNumCols;

-    RowMatrix() {}

-

 public:

     RowMatrix(unsigned nrow, unsigned ncol);

     RowMatrix(const Rcpp::NumericMatrix &rmat);

-    RowMatrix(const RowMatrix &mat);

-    RowMatrix(const ColMatrix &mat);

+    RowMatrix(const std::string &path);

     unsigned nRow() const {return mNumRows;}

     unsigned nCol() const {return mNumCols;}

@@ -105,9 +62,9 @@ public:

     RowMatrix operator/(float val) const;

+    void operator=(const RowMatrix &mat);

     void operator=(const ColMatrix &mat);

-    void update(const MatrixChange &change);

     Rcpp::NumericMatrix rMatrix() const;

     friend Archive& operator<<(Archive &ar, RowMatrix &mat);

@@ -121,14 +78,11 @@ private:

     std::vector<Vector> mCols;

     unsigned mNumRows, mNumCols;

-    ColMatrix() {}

-

 public:

     ColMatrix(unsigned nrow, unsigned ncol);

     ColMatrix(const Rcpp::NumericMatrix &rmat);

-    ColMatrix(const ColMatrix &mat);

-    ColMatrix(const RowMatrix &mat);

+    ColMatrix(const std::string &path);

     unsigned nRow() const {return mNumRows;}

     unsigned nCol() const {return mNumCols;}

@@ -144,9 +98,9 @@ public:

     ColMatrix operator/(float val) const;

+    void operator=(const ColMatrix &mat);

     void operator=(const RowMatrix &mat);

-    void update(const MatrixChange &change);

     Rcpp::NumericMatrix rMatrix() const;

     friend Archive& operator<<(Archive &ar, ColMatrix &mat);

new file mode 100644

@@ -0,0 +1,52 @@

+#include "Vector.h"

+

+void Vector::concat(const Vector& vec)

+{

+    mValues.insert(mValues.end(), vec.mValues.begin(), vec.mValues.end());

+}

+

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

+{

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

+    {

+        mValues[i] += vec[i];

+    }

+}

+

+Vector Vector::operator*(float val) const

+{

+    Vector vec(*this);

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

+    {

+        vec[i] *= val;

+    }

+    return vec;

+}

+

+Vector Vector::operator/(float val) const

+{

+    Vector vec(*this);

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

+    {

+        vec[i] /= val;

+    }

+    return vec;

+}

+

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

+{

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

+    {

+        ar << vec[i];

+    }

+    return ar;

+}

+

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

+{

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

+    {

+        ar >> vec.mValues[i];

+    }

+    return ar;

+}

new file mode 100644

@@ -0,0 +1,41 @@

+#ifndef __COGAPS_VECTOR_H__

+#define __COGAPS_VECTOR_H__

+

+#include "../Archive.h"

+

+#include <boost/align/aligned_allocator.hpp>

+#include <vector>

+

+namespace bal = boost::alignment;

+typedef std::vector<float, bal::aligned_allocator<float,32> > aligned_vector;

+

+class Vector

+{

+private:

+

+    aligned_vector mValues;

+

+public:

+

+    Vector(unsigned size) : mValues(aligned_vector(size, 0.f)) {}

+

+    const float* ptr() const {return &mValues[0];}

+    float* ptr() {return &mValues[0];}

+

+    float& operator[](unsigned i) {return mValues[i];}

+    float operator[](unsigned i) const {return mValues[i];}

+

+    unsigned size() const {return mValues.size();}

+

+    Rcpp::NumericVector rVec() const {return Rcpp::wrap(mValues);}

+    void concat(const Vector& vec);

+    void operator+=(const Vector &vec);

+

+    Vector operator*(float val) const;

+    Vector operator/(float val) const;

+

+    friend Archive& operator<<(Archive &ar, Vector &vec);

+    friend Archive& operator>>(Archive &ar, Vector &vec);

+};

+

+#endif

\ No newline at end of file

@@ -21,7 +21,7 @@

 #include <stdint.h>

-#if defined(__GAPS_OPENMP__)

+#ifdef __GAPS_OPENMP__

     #include <omp.h>

 #else

     typedef int omp_int_t;