@@ -15,7 +15,7 @@ matrix_data_t gaps::algo::sum(const Vector &vec)

     matrix_data_t sum = 0.0;

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

     {

-        sum += vec(i);

+        sum += vec[i];

     }

     return sum;

 }

@@ -72,13 +72,13 @@ ColMatrix gaps::algo::computeStdDev(const ColMatrix &stdMat,

 const ColMatrix &meanMat, unsigned nUpdates)

 {

     ColMatrix retMat(stdMat.nRow(), stdMat.nCol());

-    double meanTerm = 0.0;

+    float meanTerm = 0.0;

     for (unsigned r = 0; r < retMat.nRow(); ++r)

     {

         for (unsigned c = 0; c < retMat.nCol(); ++c)

         {

-            meanTerm = meanMat(r,c) * meanMat(r,c) / (double)nUpdates;

-            retMat(r,c) = std::sqrt((stdMat(r,c) - meanTerm) / ((double)nUpdates - 1.0));

+            meanTerm = meanMat(r,c) * meanMat(r,c) / (float)nUpdates;

+            retMat(r,c) = std::sqrt((stdMat(r,c) - meanTerm) / ((float)nUpdates - 1.0));

         }

     }

     return retMat;

@@ -88,13 +88,13 @@ RowMatrix gaps::algo::computeStdDev(const RowMatrix &stdMat,

 const RowMatrix &meanMat, unsigned nUpdates)

 {

     RowMatrix retMat(stdMat.nRow(), stdMat.nCol());

-    double meanTerm = 0.0;

+    float meanTerm = 0.0;

     for (unsigned r = 0; r < retMat.nRow(); ++r)

     {

         for (unsigned c = 0; c < retMat.nCol(); ++c)

         {

-            meanTerm = meanMat(r,c) * meanMat(r,c) / (double)nUpdates;

-            retMat(r,c) = std::sqrt((stdMat(r,c) - meanTerm) / ((double)nUpdates - 1.0));

+            meanTerm = meanMat(r,c) * meanMat(r,c) / (float)nUpdates;

+            retMat(r,c) = std::sqrt((stdMat(r,c) - meanTerm) / ((float)nUpdates - 1.0));

         }

     }

     return retMat;

@@ -108,7 +108,7 @@ const RowMatrix &B)

     {

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

         {

-            double sum = 0.0;

+            float sum = 0.0;

             for (unsigned k = 0; k < A.nCol(); ++k)

             {

                 sum += A(i,k) * B(k,j);

@@ -118,37 +118,37 @@ const RowMatrix &B)

     }

 }

-Vector gaps::algo::scalarMultiple(const Vector &vec, double n)

+Vector gaps::algo::scalarMultiple(const Vector &vec, float n)

 {

     Vector retVec(vec.size());

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

     {

-        retVec(i) = vec(i) * n;

+        retVec[i] = vec[i] * n;

     }

     return retVec;

 }

-Vector gaps::algo::squaredScalarMultiple(const Vector &vec, double n)

+Vector gaps::algo::squaredScalarMultiple(const Vector &vec, float n)

 {

     Vector retVec(vec.size());

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

     {

-        retVec(i) = GAPS_SQ(vec(i) * n);

+        retVec[i] = GAPS_SQ(vec[i] * n);

     }

     return retVec;

 }

-Vector gaps::algo::squaredScalarDivision(const Vector &vec, double n)

+Vector gaps::algo::squaredScalarDivision(const Vector &vec, float n)

 {

     Vector retVec(vec.size());

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

     {

-        retVec(i) = GAPS_SQ(vec(i) / n);

+        retVec[i] = GAPS_SQ(vec[i] / n);

     }

     return retVec;

 }

-ColMatrix gaps::algo::scalarDivision(const ColMatrix &mat, double n)

+ColMatrix gaps::algo::scalarDivision(const ColMatrix &mat, float n)

 {

     ColMatrix retMat(mat.nRow(), mat.nCol());

     for (unsigned r = 0; r < retMat.nRow(); ++r)

@@ -166,12 +166,12 @@ Vector gaps::algo::scalarDivision(const Vector &vec, matrix_data_t n)

     Vector temp(vec.size());

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

     {

-        temp(i) = vec(i) / n;

+        temp[i] = vec[i] / n;

     }

     return temp;

 }

-RowMatrix gaps::algo::scalarDivision(const RowMatrix &mat, double n)

+RowMatrix gaps::algo::scalarDivision(const RowMatrix &mat, float n)

 {

     RowMatrix retMat(mat.nRow(), mat.nCol());

     for (unsigned r = 0; r < retMat.nRow(); ++r)

@@ -187,7 +187,7 @@ RowMatrix gaps::algo::scalarDivision(const RowMatrix &mat, double n)

 matrix_data_t gaps::algo::loglikelihood(const TwoWayMatrix &D,

 const TwoWayMatrix &S, const TwoWayMatrix &AP)

 {

-    double chi2 = 0.0;

+    float chi2 = 0.0;

     for (unsigned r = 0; r < D.nRow(); ++r)

     {

         for (unsigned c = 0; c < D.nCol(); ++c)

@@ -208,12 +208,12 @@ bool gaps::algo::isColZero(const ColMatrix &mat, unsigned col)

     return gaps::algo::sum(mat.getCol(col)) == 0;

 }

-static double deltaLL_A(const TwoWayMatrix &D, const TwoWayMatrix &S,

+static float deltaLL_A(const TwoWayMatrix &D, const TwoWayMatrix &S,

 const RowMatrix &P, const TwoWayMatrix &AP, unsigned row,

-unsigned col1, double delta1, unsigned col2=0, double delta2=0.0,

+unsigned col1, float delta1, unsigned col2=0, float delta2=0.0,

 bool twoChanges=false)

 {

-    double numer = 0.0, denom = 0.0, delLL = 0.0, d1 = 0.0, d2 = 0.0;

+    float numer = 0.0, denom = 0.0, delLL = 0.0, d1 = 0.0, d2 = 0.0;

     const Vector &Drow = D.getRow(row);

     const Vector &AProw = AP.getRow(row);

     const Vector &Srow = S.getRow(row);

@@ -221,21 +221,21 @@ bool twoChanges=false)

     const Vector &Prow2 = P.getRow(col2);

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

     {

-        d1 = delta1 * Prow1(j);

-        d2 = twoChanges ? delta2 * Prow2(j) : 0.0;

-        numer = 2.0 * (Drow(j) - AProw(j)) * (d1 + d2) - GAPS_SQ(d1 + d2);

-        denom = 2.0 * GAPS_SQ(Srow(j));

+        d1 = delta1 * Prow1[j];

+        d2 = twoChanges ? delta2 * Prow2[j] : 0.0;

+        numer = 2.0 * (Drow[j] - AProw[j]) * (d1 + d2) - GAPS_SQ(d1 + d2);

+        denom = 2.0 * GAPS_SQ(Srow[j]);

         delLL += numer / denom;

     }

     return delLL;

 }

-static double deltaLL_P(const TwoWayMatrix &D, const TwoWayMatrix &S,

+static float deltaLL_P(const TwoWayMatrix &D, const TwoWayMatrix &S,

 const ColMatrix &A, const TwoWayMatrix &AP, unsigned col,

-unsigned row1, double delta1, unsigned row2=0, double delta2=0.0,

+unsigned row1, float delta1, unsigned row2=0, float delta2=0.0,

 bool twoChanges=false)

 {

-    double numer = 0.0, denom = 0.0, delLL = 0.0, d1 = 0.0, d2 = 0.0;

+    float numer = 0.0, denom = 0.0, delLL = 0.0, d1 = 0.0, d2 = 0.0;

     const Vector &Dcol = D.getCol(col);

     const Vector &APcol = AP.getCol(col);

     const Vector &Scol = S.getCol(col);

@@ -243,16 +243,16 @@ bool twoChanges=false)

     const Vector &Acol2 = A.getCol(row2);

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

     {

-        d1 = delta1 * Acol1(i);

-        d2 = twoChanges ? delta2 * Acol2(i) : 0.0;

-        numer = 2.0 * (Dcol(i) - APcol(i)) * (d1 + d2) - GAPS_SQ(d1 + d2);

-        denom = 2.0 * GAPS_SQ(Scol(i));

+        d1 = delta1 * Acol1[i];

+        d2 = twoChanges ? delta2 * Acol2[i] : 0.0;

+        numer = 2.0 * (Dcol[i] - APcol[i]) * (d1 + d2) - GAPS_SQ(d1 + d2);

+        denom = 2.0 * GAPS_SQ(Scol[i]);

         delLL += numer / denom;

     }

     return delLL;

 }

-double gaps::algo::deltaLL(const MatrixChange &ch, const TwoWayMatrix &D,

+float gaps::algo::deltaLL(const MatrixChange &ch, const TwoWayMatrix &D,

 const TwoWayMatrix &S, const ColMatrix &A, const RowMatrix &P,

 const TwoWayMatrix &AP)

 {

@@ -285,7 +285,7 @@ static AlphaParameters alphaParameters_A(const TwoWayMatrix &D,

 const TwoWayMatrix &S, const RowMatrix &P, const TwoWayMatrix &AP,

 unsigned row, unsigned col1, unsigned col2=0, bool twoChanges=false)

 {

-    double s = 0.0, su = 0.0, ratio = 0.0, p2 = 0.0;

+    float s = 0.0, su = 0.0, ratio = 0.0, p2 = 0.0;

     const Vector &Drow = D.getRow(row);

     const Vector &AProw = AP.getRow(row);

     const Vector &Srow = S.getRow(row);

@@ -293,10 +293,10 @@ unsigned row, unsigned col1, unsigned col2=0, bool twoChanges=false)

     const Vector &Prow2 = P.getRow(col2);

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

     {

-        p2 = twoChanges ? Prow2(j) : 0.0;

-        ratio = (Prow1(j) - p2) / Srow(j);

+        p2 = twoChanges ? Prow2[j] : 0.0;

+        ratio = (Prow1[j] - p2) / Srow[j];

         s += GAPS_SQ(ratio);

-        su += ratio * (Drow(j) - AProw(j)) / Srow(j);

+        su += ratio * (Drow[j] - AProw[j]) / Srow[j];

     }

     return AlphaParameters(s, su);

 }

@@ -305,7 +305,7 @@ static AlphaParameters alphaParameters_P(const TwoWayMatrix &D,

 const TwoWayMatrix &S, const ColMatrix &A, const TwoWayMatrix &AP,

 unsigned col, unsigned row1, unsigned row2=0, bool twoChanges=false)

 {

-    double s = 0.0, su = 0.0, ratio = 0.0, a2 = 0.0;

+    float s = 0.0, su = 0.0, ratio = 0.0, a2 = 0.0;

     const Vector &Dcol = D.getCol(col);

     const Vector &APcol = AP.getCol(col);

     const Vector &Scol = S.getCol(col);

@@ -313,10 +313,10 @@ unsigned col, unsigned row1, unsigned row2=0, bool twoChanges=false)

     const Vector &Acol2 = A.getCol(row2);

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

     {

-        a2 = twoChanges ? Acol2(i) : 0.0;

-        ratio = (Acol1(i) - a2) / Scol(i);

+        a2 = twoChanges ? Acol2[i] : 0.0;

+        ratio = (Acol1[i] - a2) / Scol[i];

         s += GAPS_SQ(ratio);

-        su += ratio * (Dcol(i) - APcol(i)) / Scol(i);

+        su += ratio * (Dcol[i] - APcol[i]) / Scol[i];

     }

     return AlphaParameters(s, su);

 }

@@ -2,13 +2,14 @@

 #define __COGAPS_ALGORITHMS_H__

 #include "Matrix.h"

+#include <xmmintrin.h>

 struct AlphaParameters

 {

-    double s;

-    double su;

+    float s;

+    float su;

-    AlphaParameters(double inS, double inSU)

+    AlphaParameters(float inS, float inSU)

         : s(inS), su(inSU)

     {}

@@ -20,6 +21,41 @@ struct AlphaParameters

     }

 };

+class vec4f

+{

+private:

+    

+    __m128 mValue;

+

+public:

+

+    inline vec4f() {}

+    inline vec4f(float f) : mValue(_mm_set1_ps(f)) {}

+    inline vec4f(float f0, float f1, float f2, float f3)

+        : mValue(_mm_setr_ps(f0,f1,f2,f3))

+    {}

+    inline vec4f(const __m128 &rhs) : mValue(rhs) {}

+    

+    inline vec4f& operator=(const __m128 &rhs)

+    {

+        mValue = rhs;

+        return *this;

+    }

+

+    inline operator __m128() const {return mValue;}

+

+    inline vec4f& operator+=(const vec4f &rhs)

+    {

+        *this = *this + rhs;

+        return *this;

+    }

+};

+

+inline vec4f operator+(const vec4f &lhs, const vec4f &rhs)

+{

+    return _mm_add_ps(lhs,rhs);

+}

+

 namespace gaps

 {

@@ -56,7 +92,7 @@ namespace algo

     matrix_data_t loglikelihood(const TwoWayMatrix &D, const TwoWayMatrix &S,

         const TwoWayMatrix &AP);

-    double deltaLL(const MatrixChange &ch, const TwoWayMatrix &D,

+    float deltaLL(const MatrixChange &ch, const TwoWayMatrix &D,

         const TwoWayMatrix &S, const ColMatrix &A,

         const RowMatrix &P, const TwoWayMatrix &AP);

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

 #include "AtomicSupport.h"

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

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

 AtomicSupport::AtomicSupport(char label, uint64_t nrow, uint64_t ncol,

-double alpha, double lambda)

+float alpha, float lambda)

     :

 mNumRows(nrow), mNumCols(ncol), mNumBins(nrow * ncol),

 mNumAtoms(0), mTotalMass(0.0), mLabel(label), mAlpha(alpha), 

@@ -36,7 +36,7 @@ uint64_t AtomicSupport::randomFreePosition() const

 uint64_t AtomicSupport::randomAtomPosition() const

 {

     uint64_t num = gaps::random::uniform64(0, mNumAtoms - 1);

-    std::map<uint64_t, double>::const_iterator it = mAtomicDomain.begin();

+    std::map<uint64_t, float>::const_iterator it = mAtomicDomain.begin();

     std::advance(it, num);

     return it->first;

 }

@@ -44,14 +44,14 @@ uint64_t AtomicSupport::randomAtomPosition() const

 AtomicProposal AtomicSupport::proposeDeath() const

 {

     uint64_t location = randomAtomPosition();

-    double mass = mAtomicDomain.at(location);

+    float mass = mAtomicDomain.at(location);

     return AtomicProposal(mLabel, 'D', location, -mass);

 }

 AtomicProposal AtomicSupport::proposeBirth() const

 {

     uint64_t location = randomFreePosition();

-    double mass = gaps::random::exponential(mLambda);

+    float mass = gaps::random::exponential(mLambda);

     return AtomicProposal(mLabel, 'B', location, mass);

 }

@@ -60,10 +60,10 @@ AtomicProposal AtomicSupport::proposeMove() const

 {

     // get random atom

     uint64_t location = randomAtomPosition();

-    double mass = mAtomicDomain.at(location);

+    float mass = mAtomicDomain.at(location);

     // get an iterator to this atom

-    std::map<uint64_t, double>::const_iterator it, left, right;

+    std::map<uint64_t, float>::const_iterator it, left, right;

     it = mAtomicDomain.find(location);

     left = it; right = it;

     uint64_t rbound = mMaxNumAtoms - 1, lbound = 0;

@@ -86,10 +86,10 @@ AtomicProposal AtomicSupport::proposeExchange() const

 {

     // get random atom

     uint64_t pos1 = randomAtomPosition();

-    double mass1 = mAtomicDomain.at(pos1);

+    float mass1 = mAtomicDomain.at(pos1);

     // find atom to the right, wrap around the end

-    std::map<uint64_t, double>::const_iterator it;

+    std::map<uint64_t, float>::const_iterator it;

     it = ++(mAtomicDomain.find(pos1));

     if (it == mAtomicDomain.end())

     {

@@ -98,22 +98,22 @@ AtomicProposal AtomicSupport::proposeExchange() const

     // get adjacent atom info

     uint64_t pos2 = it->first;

-    double mass2 = it->second;

+    float mass2 = it->second;

     // calculate new mass

-    double pupper = gaps::random::p_gamma(mass1 + mass2, 2.0, 1.0 / mLambda);

-    double newMass = gaps::random::q_gamma(gaps::random::uniform(0.0, pupper),

+    float pupper = gaps::random::p_gamma(mass1 + mass2, 2.0, 1.0 / mLambda);

+    float newMass = gaps::random::q_gamma(gaps::random::uniform(0.0, pupper),

         2.0, 1.0 / mLambda);

     // calculate mass changes

-    double delta1 = mass1 > mass2 ? newMass - mass1 : mass2 - newMass;

-    double delta2 = mass2 > mass1 ? newMass - mass2 : mass1 - newMass;

+    float delta1 = mass1 > mass2 ? newMass - mass1 : mass2 - newMass;

+    float delta2 = mass2 > mass1 ? newMass - mass2 : mass1 - newMass;

     // preserve total mass

     return AtomicProposal(mLabel, 'E', pos1, delta1, pos2, delta2);

 }

-double AtomicSupport::updateAtomMass(char type, uint64_t pos, double delta)

+float AtomicSupport::updateAtomMass(char type, uint64_t pos, float delta)

 {

     if (mAtomicDomain.count(pos)) // update atom if it exists

     {

@@ -121,7 +121,7 @@ double AtomicSupport::updateAtomMass(char type, uint64_t pos, double delta)

     }

     else // create a new atom if it doesn't exist

     {

-        mAtomicDomain.insert(std::pair<uint64_t, double>(pos, delta));

+        mAtomicDomain.insert(std::pair<uint64_t, float>(pos, delta));

         mNumAtoms++;

     }

@@ -142,7 +142,7 @@ AtomicProposal AtomicSupport::makeProposal() const

         return proposeBirth();

     }

-    double unif = gaps::random::uniform();

+    float unif = gaps::random::uniform();

     if ((mNumAtoms < 2 && unif <= 0.6667) || unif <= 0.5) // birth/death

     {

         if (mNumAtoms >= mMaxNumAtoms)

@@ -150,19 +150,19 @@ AtomicProposal AtomicSupport::makeProposal() const

             return proposeDeath();

         }

-        double pDelete = 0.5; // default uniform prior for mAlpha == 0

+        float pDelete = 0.5; // default uniform prior for mAlpha == 0

         if (mAlpha > 0) // poisson prior

         {

-            double dMax = (double)mMaxNumAtoms;

-            double dNum = (double)mNumAtoms;

-            double maxTerm = (dMax - dNum) / dMax;

+            float dMax = (float)mMaxNumAtoms;

+            float dNum = (float)mNumAtoms;

+            float maxTerm = (dMax - dNum) / dMax;

-            pDelete = dNum / (dNum + mAlpha * (double)mNumBins * maxTerm);

+            pDelete = dNum / (dNum + mAlpha * (float)mNumBins * maxTerm);

         }

         else if (mAlpha < 0) // geometric prior

         {

-            double c = -mAlpha * mNumBins / (-mAlpha * mNumBins + 1.0);

-            pDelete = (double)mNumAtoms / ((mNumAtoms + 1) * c + (double)mNumAtoms);

+            float c = -mAlpha * mNumBins / (-mAlpha * mNumBins + 1.0);

+            pDelete = (float)mNumAtoms / ((mNumAtoms + 1) * c + (float)mNumAtoms);

         }

         return gaps::random::uniform() < pDelete ? proposeDeath() : proposeBirth();

     }

@@ -16,16 +16,16 @@ struct AtomicProposal

     unsigned nChanges;

     uint64_t pos1;

-    double delta1;

+    float delta1;

     uint64_t pos2;

-    double delta2;

+    float delta2;

-    AtomicProposal(char l, char t, uint64_t p1, double m1)

+    AtomicProposal(char l, char t, uint64_t p1, float m1)

         : label(l), type(t), nChanges(1), pos1(p1), delta1(m1), pos2(0), delta2(0.0)

     {}

-    AtomicProposal(char l, char t, uint64_t p1, double m1, uint64_t p2, double m2)

+    AtomicProposal(char l, char t, uint64_t p1, float m1, uint64_t p2, float m2)

         : label(l), type(t), nChanges(2), pos1(p1), delta1(m1), pos2(p2), delta2(m2)

     {}

 };

@@ -42,12 +42,12 @@ public:

     char mLabel;

     // storage of the atomic domain

-    std::map<uint64_t, double> mAtomicDomain;

+    std::map<uint64_t, float> mAtomicDomain;

     uint64_t mNumAtoms;

     uint64_t mMaxNumAtoms;

     // total mass of all atoms

-    double mTotalMass;

+    float mTotalMass;

     // matrix information

     uint64_t mNumRows;

@@ -56,10 +56,10 @@ public:

     uint64_t mBinSize;

     // average number of atoms per bin (must be > 0)

-    double mAlpha;

+    float mAlpha;

     // expected magnitude of each atom (must be > 0)

-    double mLambda;

+    float mLambda;

     // convert atomic position to row/col of the matrix

     uint64_t getRow(uint64_t pos) const;

@@ -76,13 +76,13 @@ public:

     AtomicProposal proposeExchange() const;

     // update the mass of an atom, return the total amount changed

-    double updateAtomMass(char type, uint64_t pos, double delta);

+    float updateAtomMass(char type, uint64_t pos, float delta);

 public:

     // constructor

-    AtomicSupport(char label, uint64_t nrow, uint64_t ncol, double alpha=1.0,

-        double lambda=1.0);

+    AtomicSupport(char label, uint64_t nrow, uint64_t ncol, float alpha=1.0,

+        float lambda=1.0);

     // create and accept a proposal

     AtomicProposal makeProposal() const;

@@ -92,15 +92,15 @@ public:

     MatrixChange getMatrixChange(const AtomicProposal &prop) const;

     // getters

-    double alpha() const {return mAlpha;}

-    double lambda() const {return mLambda;}

-    double totalMass() const {return mTotalMass;}

+    float alpha() const {return mAlpha;}

+    float lambda() const {return mLambda;}

+    float totalMass() const {return mTotalMass;}

     uint64_t numAtoms() const {return mNumAtoms;}

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

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

     // setters

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

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

+    void setAlpha(float alpha) {mAlpha = alpha;}

+    void setLambda(float lambda) {mLambda = lambda;}

     //void setMaxNumAtoms(uint64_t max) {mMaxNumAtoms = max;}

 };

@@ -21,9 +21,9 @@ unsigned numSnapshots);

 // [[Rcpp::export]]

 Rcpp::List cogaps(Rcpp::NumericMatrix DMatrix, Rcpp::NumericMatrix SMatrix,

-unsigned nFactor, double alphaA, double alphaP, unsigned nEquil,

-unsigned nEquilCool, unsigned nSample, double maxGibbsMassA,

-double maxGibbsMassP, Rcpp::NumericMatrix fixedPatterns,

+unsigned nFactor, float alphaA, float alphaP, unsigned nEquil,

+unsigned nEquilCool, unsigned nSample, float maxGibbsMassA,

+float maxGibbsMassP, Rcpp::NumericMatrix fixedPatterns,

 char whichMatrixFixed, int seed, bool messages, bool singleCellRNASeq,

 unsigned numOutputs, unsigned numSnapshots)

 {

@@ -92,13 +92,13 @@ unsigned& nIterA, unsigned& nIterP, Vector& chi2Vec, Vector& aAtomVec,

 Vector& pAtomVec, GapsPhase phase, unsigned numOutputs, bool messages,

 unsigned numSnapshots)

 {

-    double tempDenom = (double)nIterTotal / 2.0;

+    float tempDenom = (float)nIterTotal / 2.f;

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

     {

         if (phase == GAPS_CALIBRATION)

         {

-            sampler.setAnnealingTemp(std::min(((double)i + 2.0) / tempDenom, 1.0));

+            sampler.setAnnealingTemp(std::min(((float)i + 2.f) / tempDenom, 1.f));

         }

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

@@ -121,11 +121,11 @@ unsigned numSnapshots)

             }

         }

-        double numAtomsA = sampler.totalNumAtoms('A');

-        double numAtomsP = sampler.totalNumAtoms('P');

-        aAtomVec(i) = numAtomsA;

-        pAtomVec(i) = numAtomsP;

-        chi2Vec(i) = sampler.chi2();

+        float numAtomsA = sampler.totalNumAtoms('A');

+        float numAtomsP = sampler.totalNumAtoms('P');

+        aAtomVec[i] = numAtomsA;

+        pAtomVec[i] = numAtomsP;

+        chi2Vec[i] = sampler.chi2();

         if (phase != GAPS_COOLING)

         {

@@ -137,8 +137,8 @@ unsigned numSnapshots)

                     << sampler.chi2() << '\n';

             }

-            nIterA = gaps::random::poisson(std::max(numAtomsA, 10.0));

-            nIterP = gaps::random::poisson(std::max(numAtomsP, 10.0));

+            nIterA = gaps::random::poisson(std::max(numAtomsA, 10.f));

+            nIterP = gaps::random::poisson(std::max(numAtomsP, 10.f));

         }

     }

 }

\ No newline at end of file

@@ -3,11 +3,11 @@

 #include <Rcpp.h>

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

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

 GibbsSampler::GibbsSampler(Rcpp::NumericMatrix D, Rcpp::NumericMatrix S,

-unsigned int nFactor, double alphaA, double alphaP, double maxGibbsMassA,

-double maxGibbsMassP, bool singleCellRNASeq, Rcpp::NumericMatrix fixedPat,

+unsigned int nFactor, float alphaA, float alphaP, float maxGibbsMassA,

+float maxGibbsMassP, bool singleCellRNASeq, Rcpp::NumericMatrix fixedPat,

 char whichMat)

     :

 mDMatrix(D), mSMatrix(S), mAMatrix(D.nrow(), nFactor),

@@ -19,7 +19,7 @@ mAMeanMatrix(D.nrow(), nFactor), mAStdMatrix(D.nrow(), nFactor),

 mPMeanMatrix(nFactor, D.ncol()), mPStdMatrix(nFactor, D.ncol()),

 mStatUpdates(0), mFixedMat(whichMat)

 {

-    double meanD = mSingleCellRNASeq ? gaps::algo::nonZeroMean(mDMatrix)

+    float meanD = mSingleCellRNASeq ? gaps::algo::nonZeroMean(mDMatrix)

         : gaps::algo::mean(mDMatrix);

     mADomain.setAlpha(alphaA);

@@ -59,10 +59,10 @@ Rcpp::NumericMatrix GibbsSampler::getNormedMatrix(char mat)

     Vector normVec(mPMatrix.nRow());

     for (unsigned r = 0; r < mPMatrix.nRow(); ++r)

     {

-        normVec(r) = gaps::algo::sum(mPMatrix.getRow(r));

-        if (normVec(r) == 0)

+        normVec[r] = gaps::algo::sum(mPMatrix.getRow(r));

+        if (normVec[r] == 0)

         {

-            normVec(r) = 1.0;

+            normVec[r] = 1.0;

         }

     }

@@ -72,7 +72,7 @@ Rcpp::NumericMatrix GibbsSampler::getNormedMatrix(char mat)

         for (unsigned c = 0; c < normedA.nCol(); ++c)

         {

             normedA.getCol(c) += gaps::algo::scalarMultiple(normedA.getCol(c),

-                normVec(c));

+                normVec[c]);

         }

         return normedA.rMatrix();

     }

@@ -82,13 +82,13 @@ Rcpp::NumericMatrix GibbsSampler::getNormedMatrix(char mat)

         for (unsigned r = 0; r < normedP.nRow(); ++r)

         {

             normedP.getRow(r) += gaps::algo::scalarDivision(normedP.getRow(r),

-                normVec(r));

+                normVec[r]);

         }

         return normedP.rMatrix();

     }

 }

-double GibbsSampler::getGibbsMass(const MatrixChange &change)

+float GibbsSampler::getGibbsMass(const MatrixChange &change)

 {

     // check if this change is death (only called in birth/death)

     bool death = change.delta1 < 0;

@@ -100,25 +100,25 @@ double GibbsSampler::getGibbsMass(const MatrixChange &change)

     // calculate mean and standard deviation

     alphaParam.s *= mAnnealingTemp / 2.0;

     alphaParam.su *= mAnnealingTemp / 2.0;

-    double lambda = change.label == 'A' ? mADomain.lambda() : mPDomain.lambda();

-    double mean  = (2.0 * alphaParam.su - lambda) / (2.0 * alphaParam.s);

-    double sd = 1.0 / std::sqrt(2.0 * alphaParam.s);

+    float lambda = change.label == 'A' ? mADomain.lambda() : mPDomain.lambda();

+    float mean  = (2.0 * alphaParam.su - lambda) / (2.0 * alphaParam.s);

+    float sd = 1.0 / std::sqrt(2.0 * alphaParam.s);

     // note: is bounded below by zero so have to use inverse sampling!

     // based upon algorithm in DistScalarRmath.cc (scalarRandomSample)

-    double plower = gaps::random::p_norm(0.0, mean, sd);

-    double u = gaps::random::uniform(plower, 1.0);

+    float plower = gaps::random::p_norm(0.f, mean, sd);

+    float u = gaps::random::uniform(plower, 1.f);

     // if the likelihood is flat and nonzero, sample strictly from the prior

-    double newMass = 0.0;

-    if (plower == 1 || alphaParam.s < 1.e-5)

+    float newMass = 0.f;

+    if (plower == 1.f || alphaParam.s < 0.00001f)

     {

-        newMass = death ? std::abs(change.delta1) : 0.0;

+        newMass = death ? std::abs(change.delta1) : 0.f;

     }

-    else if (plower >= 0.99)

+    else if (plower >= 0.99f)

     {

-        double tmp1 = gaps::random::d_norm(0, mean, sd);

-        double tmp2 = gaps::random::d_norm(10 * lambda, mean, sd);

+        float tmp1 = gaps::random::d_norm(0.f, mean, sd);

+        float tmp2 = gaps::random::d_norm(10.f * lambda, mean, sd);

         if (tmp1 > EPSILON && std::abs(tmp1 - tmp2) < EPSILON)

         {

@@ -133,10 +133,10 @@ double GibbsSampler::getGibbsMass(const MatrixChange &change)

     newMass = (change.label == 'A' ? std::min(newMass, mMaxGibbsMassA)

         : std::min(newMass, mMaxGibbsMassP));

-    return std::max(newMass, 0.0);

+    return std::max(newMass, 0.f);

 }

-double GibbsSampler::computeDeltaLL(const MatrixChange &change)

+float GibbsSampler::computeDeltaLL(const MatrixChange &change)

 {

     return gaps::algo::deltaLL(change, mDMatrix, mSMatrix, mAMatrix,

         mPMatrix, mAPMatrix);

@@ -178,26 +178,26 @@ uint64_t GibbsSampler::totalNumAtoms(char matrixLabel) const

     return matrixLabel == 'A' ? mADomain.numAtoms() : mPDomain.numAtoms();

 }

-void GibbsSampler::setChi2(double chi2)

+void GibbsSampler::setChi2(float chi2)

 {

     mChi2 = chi2;

 }

-double GibbsSampler::chi2() const

+float GibbsSampler::chi2() const

 {

     return mChi2;

 }

-void GibbsSampler::setAnnealingTemp(double temp)

+void GibbsSampler::setAnnealingTemp(float temp)

 {

     mAnnealingTemp = temp;

 }

 bool GibbsSampler::evaluateChange(AtomicSupport &domain,

-const AtomicProposal &proposal, double threshold, bool accept)

+const AtomicProposal &proposal, float threshold, bool accept)

 {

     MatrixChange change = domain.getMatrixChange(proposal);

-    double delLL = computeDeltaLL(change);

+    float delLL = computeDeltaLL(change);

     if (accept || delLL * mAnnealingTemp >= threshold)

     {

         change = domain.acceptProposal(proposal);

@@ -209,23 +209,23 @@ const AtomicProposal &proposal, double threshold, bool accept)

     return false;

 }

-void GibbsSampler::updateAPMatrix_A(unsigned row, unsigned col, double delta)

+void GibbsSampler::updateAPMatrix_A(unsigned row, unsigned col, float delta)

 {

     const Vector &APvec = mAPMatrix.getRow(row);

     const Vector &Pvec = mPMatrix.getRow(col);

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

     {

-        mAPMatrix.set(row, j, APvec(j) + delta * Pvec(j));

+        mAPMatrix.set(row, j, APvec[j] + delta * Pvec[j]);

     }

 }

-void GibbsSampler::updateAPMatrix_P(unsigned row, unsigned col, double delta)

+void GibbsSampler::updateAPMatrix_P(unsigned row, unsigned col, float delta)

 {

     const Vector &APvec = mAPMatrix.getCol(col);

     const Vector &Avec = mAMatrix.getCol(row);

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

     {

-        mAPMatrix.set(i, col, APvec(i) + delta * Avec(i));

+        mAPMatrix.set(i, col, APvec[i] + delta * Avec[i]);

     }

 }

@@ -271,7 +271,7 @@ bool GibbsSampler::death(AtomicSupport &domain, AtomicProposal &proposal)

     // rebirth

     MatrixChange ch = domain.getMatrixChange(proposal);

-    double newMass = canUseGibbs(ch) ? getGibbsMass(ch) : 0.0;

+    float newMass = canUseGibbs(ch) ? getGibbsMass(ch) : 0.0;

     proposal.delta1 = newMass < EPSILON ? -proposal.delta1 : newMass;    

     // attempt to accept rebirth

@@ -309,11 +309,11 @@ bool GibbsSampler::exchange(AtomicSupport &domain, AtomicProposal &proposal)

     uint64_t pos1 = proposal.delta1 > 0 ? proposal.pos1 : proposal.pos2;

     uint64_t pos2 = proposal.delta1 > 0 ? proposal.pos2 : proposal.pos1;

-    double mass1 = domain.at(pos1);

-    double mass2 = domain.at(pos2);

+    float mass1 = domain.at(pos1);

+    float mass2 = domain.at(pos2);

-    double newMass1 = mass1 + std::max(proposal.delta1, proposal.delta2);

-    double newMass2 = mass2 + std::min(proposal.delta1, proposal.delta2);

+    float newMass1 = mass1 + std::max(proposal.delta1, proposal.delta2);

+    float newMass2 = mass2 + std::min(proposal.delta1, proposal.delta2);

     unsigned row1 = change.delta1 > 0 ? change.row1 : change.row2;

     unsigned row2 = change.delta1 > 0 ? change.row2 : change.row1;

@@ -337,14 +337,14 @@ bool GibbsSampler::exchange(AtomicSupport &domain, AtomicProposal &proposal)

         if (alphaParam.s > EPSILON)

         {

-            double mean = alphaParam.su / alphaParam.s;

-            double sd = 1.0 / std::sqrt(alphaParam.s);

+            float mean = alphaParam.su / alphaParam.s;

+            float sd = 1.f / std::sqrt(alphaParam.s);

-            double plower = gaps::random::p_norm(-mass1, mean, sd);

-            double pupper = gaps::random::p_norm(mass2, mean, sd);

-            double u = gaps::random::uniform(plower, pupper);

+            float plower = gaps::random::p_norm(-mass1, mean, sd);

+            float pupper = gaps::random::p_norm(mass2, mean, sd);

+            float u = gaps::random::uniform(plower, pupper);

-            if (!(plower >  0.95 || pupper < 0.05))

+            if (!(plower >  0.95f || pupper < 0.05f))

             {

                 proposal.delta1 = gaps::random::q_norm(u, mean, sd);

                 proposal.delta1 = std::max(proposal.delta1, -mass1);

@@ -358,11 +358,11 @@ bool GibbsSampler::exchange(AtomicSupport &domain, AtomicProposal &proposal)

         }

     }

-    double pnewMass = mass1 > mass2 ? newMass1 : newMass2;

-    double poldMass = newMass1 > newMass2 ? mass1 : mass2;

+    float pnewMass = mass1 > mass2 ? newMass1 : newMass2;

+    float poldMass = newMass1 > newMass2 ? mass1 : mass2;

-    double pnew = gaps::random::d_gamma(pnewMass, 2.0, 1.0 / domain.lambda());

-    double pold = gaps::random::d_gamma(poldMass, 2.0, 1.0 / domain.lambda());

+    float pnew = gaps::random::d_gamma(pnewMass, 2.0, 1.f / domain.lambda());

+    float pold = gaps::random::d_gamma(poldMass, 2.0, 1.f / domain.lambda());

     proposal.pos1 = pos1;

     proposal.delta1 = newMass1 - mass1;

@@ -375,7 +375,7 @@ bool GibbsSampler::exchange(AtomicSupport &domain, AtomicProposal &proposal)

     }

     else

     {

-        double priorLL = pold == 0.0 ? 0.0 : log(pnew / pold);

+        float priorLL = pold == 0.0 ? 0.0 : log(pnew / pold);

         return evaluateChange(domain, proposal, std::log(gaps::random::uniform())

             - priorLL);

     }

@@ -412,28 +412,28 @@ void GibbsSampler::updateStatistics()

     for (unsigned r = 0; r < mPMatrix.nRow(); ++r)

     {

-        normVec(r) = gaps::algo::sum(mPMatrix.getRow(r));

-        if (normVec(r) == 0)

+        normVec[r] = gaps::algo::sum(mPMatrix.getRow(r));

+        if (normVec[r] == 0)

         {

-            normVec(r) = 1.0;

+            normVec[r] = 1.0;

         }

     }

     for (unsigned c = 0; c < mAMeanMatrix.nCol(); ++c)

     {

         mAMeanMatrix.getCol(c) += gaps::algo::scalarMultiple(mAMatrix.getCol(c),

-            normVec(c));

+            normVec[c]);

         mAStdMatrix.getCol(c) += gaps::algo::squaredScalarMultiple(mAMatrix.getCol(c),

-            normVec(c));

+            normVec[c]);

     }

     for (unsigned r = 0; r < mPMeanMatrix.nRow(); ++r)

     {

         mPMeanMatrix.getRow(r) += gaps::algo::scalarDivision(mPMatrix.getRow(r),

-            normVec(r));

+            normVec[r]);

         mPStdMatrix.getRow(r) += gaps::algo::squaredScalarDivision(mPMatrix.getRow(r),

-            normVec(r));

+            normVec[r]);

     }

 }

\ No newline at end of file

@@ -23,11 +23,11 @@ public:

     ColMatrix mAMeanMatrix, mAStdMatrix;

     unsigned mStatUpdates;

-    double mMaxGibbsMassA;

-    double mMaxGibbsMassP;

+    float mMaxGibbsMassA;

+    float mMaxGibbsMassP;

-    double mAnnealingTemp;

-    double mChi2;

+    float mAnnealingTemp;

+    float mChi2;

     bool mSingleCellRNASeq;

@@ -40,30 +40,30 @@ public:

     bool exchange(AtomicSupport &domain, AtomicProposal &proposal);

     bool evaluateChange(AtomicSupport &domain, const AtomicProposal &proposal,

-        double threshold, bool accept=false);

+        float threshold, bool accept=false);