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

 #define __COGAPS_GAPS_DISPATCHER_H__

 #include "GapsRunner.h"

+#include "math/SIMD.h"

 #include <string>

@@ -13,8 +13,8 @@

 GapsRunner::GapsRunner(const RowMatrix &data, unsigned nPatterns, float alphaA,

 float alphaP, float maxGibbsMassA, float maxGibbsMassP, bool singleCell)

     :

-mASampler(data, data * 0.1f, nPatterns, alphaA, maxGibbsMassA),

-mPSampler(data, data * 0.1f, nPatterns, alphaP, maxGibbsMassP),

+mASampler(data, data.pmax(0.1f), nPatterns, alphaA, maxGibbsMassA),

+mPSampler(data, data.pmax(0.1f), nPatterns, alphaP, maxGibbsMassP),

 mStatistics(data.nRow(), data.nCol(), nPatterns),

 mSamplePhase(false), mNumUpdatesA(0), mNumUpdatesP(0),

 mNumRows(data.nRow()), mNumCols(data.nCol())

@@ -76,6 +76,7 @@ unsigned nCores)

         {

             displayStatus(outputFreq, i, nIter);

         }

+

         if (mSamplePhase)

         {

             mStatistics.update(mASampler, mPSampler);

@@ -103,6 +103,19 @@ RowMatrix RowMatrix::operator/(float val) const

     return mat;

 }

+RowMatrix RowMatrix::pmax(float scale) const

+{

+    RowMatrix mat(mNumRows, mNumCols);

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

+    {

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

+        {

+            mat(i,j) = std::max(this->operator()(i,j) * scale, scale);

+        }

+    }

+    return mat;

+}

+

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

 {

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

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

     RowMatrix operator*(float val) const;

     RowMatrix operator/(float val) const;

     RowMatrix& operator=(const ColMatrix &mat);

+    RowMatrix pmax(float scale) const;

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

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