@@ -33,7 +33,7 @@ public:

             *this >> magicNum;

             if (magicNum != ARCHIVE_MAGIC_NUM)

             {

-                Rcpp::Rcout << "warning: invalid/old checkpoint file" << std::endl;

+                Rcpp::Rcout << "warning: invalid checkpoint file" << std::endl;

             }

         }*/

     }

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

 #include "GapsAssert.h"

 #include "AtomicDomain.h"

-#include "Random.h"

+#include "math/Random.h"

 #include <stdint.h>

 #include <utility>

@@ -1,4 +1,4 @@

-#include "SIMD.h"

+#include "math/SIMD.h"

 #include "GapsRunner.h"

 #include <Rcpp.h>

@@ -2,7 +2,7 @@

 #define __COGAPS_GAPS_RUNNER_H__

 #include "Archive.h"

-#include "Matrix.h"

+#include "math/Matrix.h"

 #include "GibbsSampler.h"

 #include "GapsStatistics.h"

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

 #include "GapsStatistics.h"

-#include "Algorithms.h"

+#include "math/Algorithms.h"

 GapsStatistics::GapsStatistics(unsigned nRow, unsigned nCol, unsigned nFactor)

     : mAMeanMatrix(nRow, nFactor), mAStdMatrix(nRow, nFactor),

@@ -2,7 +2,7 @@

 #define __COGAPS_GAPS_STATISTICS_H__

 #include "GibbsSampler.h"

-#include "Matrix.h"

+#include "math/Matrix.h"

 class GapsStatistics

 {

@@ -42,15 +42,9 @@ void AmplitudeGibbsSampler::sync(PatternGibbsSampler &sampler)

 void AmplitudeGibbsSampler::updateAPMatrix(unsigned row, unsigned col, float delta)

 {

-    //Rprintf("ap_update - %d %d %.12f\n", row, col, delta);

-    //Rprintf("ap_update - %d %d %x\n", row, col, *((uint32_t*)(&delta)));

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

     {

-        float temp = (*mOtherMatrix)(col,j);

-        //Rprintf("P[j] - %x\n", *((uint32_t*)(&temp)));

         mAPMatrix(row,j) += delta * (*mOtherMatrix)(col,j);

-        //float temp = mAPMatrix(row,j);

-        //Rprintf("ap_update_fin - %x\n", *((uint32_t*)(&temp)));

     }

 }

@@ -141,14 +135,9 @@ void PatternGibbsSampler::sync(AmplitudeGibbsSampler &sampler)

 void PatternGibbsSampler::updateAPMatrix(unsigned row, unsigned col, float delta)

 {

-    //Rprintf("ap_update - %d %d %x\n", row, col, *((uint32_t*)(&delta)));

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

     {

-        float temp = (*mOtherMatrix)(i,row);

-        //Rprintf("A[i] - %x\n", *((uint32_t*)(&temp)));

         mAPMatrix(i,col) += delta * (*mOtherMatrix)(i,row);

-        //float temp = mAPMatrix(i,col);

-        //Rprintf("ap_update_fin - %x\n", *((uint32_t*)(&temp)));

     }

 }

@@ -3,9 +3,9 @@

 #include "GapsAssert.h"

 #include "Archive.h"

-#include "Matrix.h"

-#include "Random.h"

-#include "Algorithms.h"

+#include "math/Matrix.h"

+#include "math/Random.h"

+#include "math/Algorithms.h"

 #include "ProposalQueue.h"

 #include "AtomicDomain.h"

@@ -230,7 +230,6 @@ unsigned col)

         : temp;

     if (mass >= gaps::algo::epsilon)

     {

-        //Rprintf("B - %x\n", *((uint32_t*)&mass));

         addMass(pos, mass, row, col);

     }

 }

@@ -247,7 +246,6 @@ unsigned col)

     float deltaLL = impl()->computeDeltaLL(row, col, mass);

     if (deltaLL * mAnnealingTemp >= std::log(gaps::random::uniform()))

     {

-        //Rprintf("D - %x\n", *((uint32_t*)&mass));

         addMass(pos, mass, row, col);

     }

 }

@@ -262,8 +260,6 @@ unsigned r1, unsigned c1, unsigned r2, unsigned c2)

         float deltaLL = impl()->computeDeltaLL(r1, c1, -mass, r2, c2, mass);

         if (deltaLL * mAnnealingTemp > std::log(gaps::random::uniform()))

         {

-            float temp = -mass;

-            //Rprintf("%x %x\n", *((uint32_t*)(&temp)), *((uint32_t*)(&mass)));

             removeMass(src, mass, r1, c1);

             addMass(dest, mass, r2, c2);

         }

@@ -368,7 +364,6 @@ template <class T, class MatA, class MatB>

 float GibbsSampler<T, MatA, MatB>::gibbsMass(unsigned row, unsigned col, float mass)

 {        

     AlphaParameters alpha = impl()->alphaParameters(row, col);

-    //Rprintf("alpha param - %f %f\n", alpha.s, alpha.su);

     alpha.s *= mAnnealingTemp / 2.0;

     alpha.su *= mAnnealingTemp / 2.0;

     float mean = (2.0 * alpha.su - mLambda) / (2.0 * alpha.s);

@@ -402,7 +397,6 @@ float GibbsSampler<T, MatA, MatB>::gibbsMass(unsigned r1, unsigned c1, float m1,

 unsigned r2, unsigned c2, float m2)

 {

     AlphaParameters alpha = impl()->alphaParameters(r1, c1, r2, c2);

-    //Rprintf("alpha param - %f %f\n", alpha.s, alpha.su);

     alpha.s *= mAnnealingTemp;

     alpha.su *= mAnnealingTemp;

@@ -1,16 +1,15 @@

-PKG_CPPFLAGS = -DGAPS_DEBUG -DBOOST_MATH_PROMOTE_DOUBLE_POLICY=0

-PKG_CPPFLAGS = -DBOOST_MATH_PROMOTE_DOUBLE_POLICY=0

-OBJECTS =   Algorithms.o \

-            AtomicDomain.o \

+PKG_CPPFLAGS =  -DBOOST_MATH_PROMOTE_DOUBLE_POLICY=0

+OBJECTS =   AtomicDomain.o \

             Cogaps.o \

             GapsRunner.o \

             GapsStatistics.o \

             GibbsSampler.o \

-            Matrix.o \

             ProposalQueue.o \

-            Random.o \

             RcppExports.o \

             test-runner.o \

+            math/Algorithms.o \

+            math/Matrix.o \

+            math/Random.o \

             cpp_tests/testAlgorithms.o \

             cpp_tests/testAtomicDomain.o \

             cpp_tests/testGapsRunner.o \

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

 #include "GapsAssert.h"

 #include "ProposalQueue.h"

-#include "Random.h"

+#include "math/Random.h"

 void ProposalQueue::setNumBins(unsigned nBins)

 {

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

 #include "catch.h"

-#include "../Matrix.h"

-#include "../Algorithms.h"

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

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

 #define MAT_SUM(nR, nC) ((nR + nC - 2) * nR * nC / 2.f)

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

 #include "catch.h"

 #include "../GibbsSampler.h"

-#include "../Algorithms.h"

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

 #include <Rcpp.h>

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

 #include "catch.h"

-#include "../Matrix.h"

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

 TEST_CASE("Test Matrix.h")

 {

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

 #include "catch.h"

-#include "../Random.h"

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

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

@@ -1,10 +1,10 @@

 #include "catch.h"

 #include "../Archive.h"

-#include "../Matrix.h"

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

 #include "../AtomicSupport.h"

 #include "../GibbsSampler.h"

 #include "../InternalState.h"

-#include "../Random.h"

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

 #if 0

similarity index 87%

rename from src/Algorithms.cpp

rename to src/math/Algorithms.cpp

@@ -98,7 +98,6 @@ bool gaps::algo::isColZero(const ColMatrix &mat, unsigned col)

 AlphaParameters gaps::algo::alphaParameters(unsigned size, const float *D,

 const float *S, const float *AP, const float *mat)

 {

-/*

     gaps::simd::packedFloat ratio, pMat, pD, pAP, pS;

     gaps::simd::packedFloat partialS = 0.f, partialSU = 0.f;

     gaps::simd::Index i = 0;

@@ -120,25 +119,12 @@ const float *S, const float *AP, const float *mat)

         su += (fratio * (D[j] - AP[j])) / S[j];

     }

     return AlphaParameters(s,su);

-    float fratio, s = partialS.scalar(), su = partialSU.scalar();

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

-*/

-    float fratio = 0.f, s = 0.f, su = 0.f;

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

-    {

-        fratio = mat[j] / S[j];

-        s += fratio * fratio;

-        su += fratio * (D[j] - AP[j]) / S[j];

-        //Rprintf("AP[j] - %x\n", *((uint32_t*)(&AP[j])));

-    }

-    return AlphaParameters(s,su);

 }

 //

 AlphaParameters gaps::algo::alphaParameters(unsigned size, const float *D,

 const float *S, const float *AP, const float *mat1, const float *mat2)

 {

-/*

     gaps::simd::packedFloat ratio, pMat1, pMat2, pD, pAP, pS;

     gaps::simd::packedFloat partialS = 0.f, partialSU = 0.f;

     gaps::simd::Index i = 0;

@@ -153,16 +139,13 @@ const float *S, const float *AP, const float *mat1, const float *mat2)

         partialS += ratio * ratio;

         partialSU += ratio * (pD - pAP) / pS;

     }

+

     float fratio, s = partialS.scalar(), su = partialSU.scalar();

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

-*/

-    float fratio = 0.f, s = 0.f, su = 0.f;

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

     {

         fratio = (mat1[j] - mat2[j]) / S[j];

         s += fratio * fratio;

         su += fratio * (D[j] - AP[j]) / S[j];

-        //Rprintf("AP[j] - %x\n", *((uint32_t*)(&AP[j])));

     }

     return AlphaParameters(s,su);

 }

similarity index 100%

rename from src/Algorithms.h

rename to src/math/Algorithms.h

similarity index 100%

rename from src/Matrix.cpp

rename to src/math/Matrix.cpp

similarity index 99%

rename from src/Matrix.h

rename to src/math/Matrix.h

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

 #ifndef __COGAPS_MATRIX_H__

 #define __COGAPS_MATRIX_H__

-#include "Archive.h"

+#include "../Archive.h"

 #include <Rcpp.h>

 #include <boost/align/aligned_allocator.hpp>

similarity index 99%

rename from src/Random.cpp

rename to src/math/Random.cpp

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

 // [[Rcpp::depends(BH)]]

 #include "Random.h"

-#include "GapsAssert.h"

+#include "../GapsAssert.h"

 #include <boost/random/uniform_01.hpp>

 #include <boost/random/uniform_real_distribution.hpp>

similarity index 96%

rename from src/Random.h

rename to src/math/Random.h

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

 #ifndef __COGAPS_RANDOM_H__

 #define __COGAPS_RANDOM_H__

-#include "Archive.h"

+#include "../Archive.h"

 #include <stdint.h>

 #include <vector>

similarity index 100%

rename from src/SIMD.h

rename to src/math/SIMD.h