@@ -198,22 +198,21 @@ float delta2)

     return delLL;

 }

-/*

 // horribly slow, don't call often

-void gaps::algo::matrixMultiplication(TwoWayMatrix &C, const ColMatrix &A,

-const RowMatrix &B)

+RowMatrix gaps::algo::matrixMultiplication(const ColMatrix &A, const RowMatrix &B)

 {

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

+    RowMatrix temp(A.nRow(), B.nCol());

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

     {

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

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

         {

             float sum = 0.0;

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

             {

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

             }

-            C.set(i, j, sum);

+            temp(i, j) = sum;

         }

     }

-}

-*/

+    return temp;

+}

\ No newline at end of file

@@ -16,9 +16,7 @@ struct AlphaParameters

     AlphaParameters operator+(const AlphaParameters &other) const

     {

-        float rs = s + other.s;

-        float rsu = su - other.su; // weird

-        return AlphaParameters(rs, rsu);

+        return AlphaParameters(s + other.s, su - other.su);

     }

 };

@@ -54,8 +52,7 @@ namespace algo

     // specific matrix algorithms

     bool isRowZero(const RowMatrix &mat, unsigned row); // TODO take pointer

     bool isColZero(const ColMatrix &mat, unsigned col);

-    //void matrixMultiplication(TwoWayMatrix &C, const ColMatrix &A,

-    //    const RowMatrix &B);

+    RowMatrix matrixMultiplication(const ColMatrix &A, const RowMatrix &B);

     // chiSq / 2

     template <class Matrix>

@@ -72,7 +72,7 @@ Rcpp::List GapsRunner::run()

     chi2Vec.concat(mChiSqSample);

     // print final chi-sq value

-    float meanChiSq = mStatistics.meanChiSq();

+    float meanChiSq = mStatistics.meanChiSq(mASampler);

     if (mPrintMessages)

     {

         Rprintf("Chi-Squared of Mean: %.2f\n", meanChiSq);

@@ -50,11 +50,12 @@ Rcpp::NumericMatrix GapsStatistics::PStd() const

         mStatUpdates).rMatrix();

 }

-float GapsStatistics::meanChiSq() const

+float GapsStatistics::meanChiSq(const AmplitudeGibbsSampler &ASampler) const

 {

-    //ColMatrix A = mAMeanMatrix / mStatUpdates;

-    //RowMatrix P = mPMeanMatrix / mStatUpdates;

-    //RowMatrix M(A * P);

-    return 0.f;

+    ColMatrix A = mAMeanMatrix / mStatUpdates;

+    RowMatrix P = mPMeanMatrix / mStatUpdates;

+    RowMatrix M(gaps::algo::matrixMultiplication(A, P));

+    return 2.f * gaps::algo::loglikelihood(ASampler.mDMatrix, ASampler.mSMatrix,

+        M);

 }

@@ -25,7 +25,7 @@ public:

     Rcpp::NumericMatrix PMean() const;

     Rcpp::NumericMatrix PStd() const;

-    float meanChiSq() const;

+    float meanChiSq(const AmplitudeGibbsSampler &ASampler) const;

     void update(const AmplitudeGibbsSampler &ASampler,

         const PatternGibbsSampler &PSampler);

@@ -154,7 +154,8 @@ mAnnealingTemp(0.f), mNumRows(nrow), mNumCols(ncol)

         / static_cast<uint64_t>(mNumRows * mNumCols);

     uint64_t remain = std::numeric_limits<uint64_t>::max()

         % static_cast<uint64_t>(mNumRows * mNumCols);

-    mQueue.setDomainSize(std::numeric_limits<uint64_t>::max() - remain);

+    //mQueue.setDomainSize(std::numeric_limits<uint64_t>::max() - remain);

+    mQueue.setDomainSize(std::numeric_limits<uint64_t>::max());

 }

 template <class T, class MatA, class MatB>

@@ -35,7 +35,7 @@ AtomicProposal ProposalQueue::makeProposal(const AtomicDomain &domain)

     float bdProb = domain.size() < 2 ? 0.6667f : 0.5f;

     float u = gaps::random::uniform();

-    if (u < bdProb)

+    if (u <= bdProb)

     {

         return gaps::random::uniform() < deathProb(domain.size()) ? 

             death(domain) : birth(domain);

@@ -68,7 +68,7 @@ AtomicProposal ProposalQueue::move(const AtomicDomain &domain)

     Atom a = domain.randomAtom();

     uint64_t lbound = domain.hasLeft(a) ? domain.left(a).pos : 0;

     uint64_t rbound = domain.hasRight(a) ? domain.right(a).pos : mDomainSize;

-    uint64_t newLocation = gaps::random::uniform64(lbound + 1, rbound - 1);

+    uint64_t newLocation = gaps::random::uniform64(lbound, rbound - 1);

     GAPS_ASSERT(domain.test(a.pos));

     return AtomicProposal('M', a.pos, a.mass, newLocation);

 }