[Genabel-commits] r1532 - branches/ProbABEL-0.50/src

[noreply at r-forge.r-project.org]

Author: maartenk
Date: 2014-01-10 21:29:35 +0100 (Fri, 10 Jan 2014)
New Revision: 1532

Modified:
   branches/ProbABEL-0.50/src/reg1.cpp
Log:
added residual variance calculation in EIGEN specific code. WARNING:this probaly brakes MMSCORE

Modified: branches/ProbABEL-0.50/src/reg1.cpp
===================================================================
--- branches/ProbABEL-0.50/src/reg1.cpp	2014-01-10 20:21:52 UTC (rev 1531)
+++ branches/ProbABEL-0.50/src/reg1.cpp	2014-01-10 20:29:35 UTC (rev 1532)
@@ -423,6 +423,12 @@
 #endif
     // now compute residual variance
     sigma2 = 0.;
+#if EIGEN
+//TODO: fix this part(residual variance) for mmscore this part: eigen part is not equel to non eigen code
+    mematrix<double>  sigma2_matrix;
+    sigma2_matrix.data=rdata.Y.data-(X.data*beta.data);
+    sigma2= sigma2_matrix.data.col(0).array().square().sum();
+#else
     mematrix<double> ttX = transpose(tX);
     mematrix<double> sigma2_matrix = rdata.Y;
     mematrix<double> sigma2_matrix1 = ttX * beta;
@@ -438,14 +444,15 @@
     //  std::cout << "sigma2_matrix.nrow=" << sigma2_matrix.nrow
     //            << "sigma2_matrix.ncol" << sigma2_matrix.ncol
     //            <<"\n";
-
+    static double val;
     for (int i = 0; i < sigma2_matrix.nrow; i++)
-    {
-        double val = sigma2_matrix.get(i, 0);
-        //            std::cout << "val = " << val << "\n";
-        sigma2 += val * val;
-        //            std::cout << "sigma2+= " << sigma2 << "\n";
-    }
+     {
+         val = sigma2_matrix.get(i, 0);
+         //            std::cout << "val = " << val << "\n";
+         sigma2 += val * val;
+         //            std::cout << "sigma2+= " << sigma2 << "\n";
+     }
+#endif
 
     double sigma2_internal = sigma2 / (N - static_cast<double>(length_beta));
     // now compute residual variance