77 sparseoptimization pattern discrepancy

src/gibbs_sampler/DenseNormalModel.h

@@ -53,9 +53,9 @@ class DenseNormalModel
     AlphaParameters alphaParametersWithChange(unsigned row, unsigned col, float ch);
     void updateAPMatrix(unsigned row, unsigned col, float delta);
 
-    Matrix mDMatrix; // samples by genes for A, genes by samples for P
-    Matrix mMatrix; // genes by patterns for A, samples by patterns for P
-    const Matrix *mOtherMatrix; // pointer to P if this is A, and vice versa
+    Matrix mDMatrix; // Data Matrix, samples x genes or genes x samples
+    Matrix mMatrix; // A (left mult) or P (right mult) based on mDMatrix
+    const Matrix *mOtherMatrix; // pointer to vis-a-vis of mMartix
     Matrix mSMatrix; // uncertainty values for each data point
     Matrix mAPMatrix; // cached product of A and P
     float mMaxGibbsMass;

src/gibbs_sampler/SparseNormalModel.cpp

@@ -36,29 +36,32 @@ void SparseNormalModel::extraInitialization()
     // nop - not needed
 }
 
+
 float SparseNormalModel::chiSq() const
 {
     float chisq = 0.f;
-    for (unsigned j = 0; j < mDMatrix.nCol(); ++j)
+    for (unsigned j = 0; j < mDMatrix.nCol(); ++j)                                          // for each col of Data Matrix
     {
-        for (unsigned i = 0; i < mDMatrix.nRow(); ++i)
+        for (unsigned i = 0; i < mDMatrix.nRow(); ++i)                                      // for each row of Data Matrix
         {
-            float dot = gaps::dot(mMatrix.getRow(j), mOtherMatrix->getRow(i));
-            chisq += dot * dot;
+            float dot = gaps::dot(mMatrix.getRow(j), mOtherMatrix->getRow(i));              // dot = A * P
+            chisq += dot * dot;                                                             // chisq = chisq+dot^2
         }
 
-        SparseIterator<1> it(mDMatrix.getCol(j));
+        SparseIterator<1> it(mDMatrix.getCol(j));                                           // iterate non-zero elements of col j
         while (!it.atEnd())
         {
-            float dot = gaps::dot(mMatrix.getRow(j), mOtherMatrix->getRow(it.getIndex()));
-            float dsq = get<1>(it) * get<1>(it);
-            chisq += 1 + dot * (dot - 2 * get<1>(it) - dsq * dot) / dsq;
+            float val = get<1>(it);                                                         // val = Data
+            float dot = gaps::dot(mMatrix.getRow(j), mOtherMatrix->getRow(it.getIndex()));  // dot = A * P again
+            float dsq = val * val;                                                          // dsq = Data squared  
+            chisq += (1 + dot * (dot - 2 * val - dsq * dot)) / (1 + dsq);                   // chisq
             it.next();
         }
     }
-    return chisq * mBeta;
+    return chisq * mBeta;                                                                   // mBeta is 100
 }
 
+
 float SparseNormalModel::dataSparsity() const
 {
     return gaps::sparsity(mDMatrix);

src/gibbs_sampler/SparseNormalModel.h

@@ -54,9 +54,9 @@ class SparseNormalModel
     AlphaParameters alphaParameters(unsigned r1, unsigned c1, unsigned r2, unsigned c2);
     AlphaParameters alphaParametersWithChange(unsigned row, unsigned col, float ch);
 
-    SparseMatrix mDMatrix; // samples by genes for A, genes by samples for P
-    HybridMatrix mMatrix; // genes by patterns for A, samples by patterns for P
-    const HybridMatrix *mOtherMatrix; // pointer to P if this is A, and vice versa
+    SparseMatrix mDMatrix; // Data Matrix D, samp x genes or genes x samp
+    HybridMatrix mMatrix; // A (left mult) or P (right mult) based on mDMatrix
+    const HybridMatrix *mOtherMatrix; // pointer to vis-a-vis of mMartix
     Matrix mZ2;
     Vector mZ1;
     float mBeta;