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

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

Author: maartenk
Date: 2014-03-01 01:26:29 +0100 (Sat, 01 Mar 2014)
New Revision: 1627

Modified:
   branches/ProbABEL-0.50/src/reg1.cpp
   branches/ProbABEL-0.50/src/reg1.h
Log:
using predefined matrices in mmscore method are now in a seperate function. This still does not live up to my expectations.

Modified: branches/ProbABEL-0.50/src/reg1.cpp
===================================================================
--- branches/ProbABEL-0.50/src/reg1.cpp	2014-02-27 22:20:11 UTC (rev 1626)
+++ branches/ProbABEL-0.50/src/reg1.cpp	2014-03-01 00:26:29 UTC (rev 1627)
@@ -310,6 +310,38 @@
     chi2_score = chi2[0];
 }
 
+void linear_reg::mmscore_regression(const MatrixXd& X,
+        const MatrixXd& W, LDLT<MatrixXd>& Ch) {
+    MatrixXd::ConstColXpr Y(reg_data.Y.data.col(0));
+    //VectorXd Y = reg_data.Y.data.col(0);
+    if (X.cols() == 3)
+    {
+        Matrix<double, 3, Dynamic> tXW = X.transpose()* W;
+        Matrix3d xWx = tXW * X;
+        Ch = LDLT<MatrixXd>(xWx);
+        Vector3d beta_3f = Ch.solve(tXW * Y);
+        sigma2 = (Y - tXW.transpose() * beta_3f).squaredNorm();
+        beta.data = beta_3f;
+    }
+    else if (X.cols() == 2)
+    {
+        Matrix<double, 2, Dynamic> tXW = X.transpose()* W;
+        Matrix2d xWx = tXW * X;
+        Ch = LDLT<MatrixXd>(xWx);
+        Vector2d beta_2f = Ch.solve(tXW * Y);
+        sigma2 = (Y - tXW.transpose() * beta_2f).squaredNorm();
+        beta.data = beta_2f;
+    }
+    else
+    {
+        // next line is  5997000 flops
+        MatrixXd tXW = X.transpose() * W;
+        Ch = LDLT<MatrixXd>(tXW * X); // 17991 flops
+        beta.data = Ch.solve(tXW * Y); //5997 flops
+        //next line is: 1000+5000+3000= 9000 flops
+        sigma2 = (Y - tXW.transpose() * beta.data).squaredNorm();
+    }
+}
 
 void linear_reg::estimate(int verbose, double tol_chol,
         int model, int interaction, int ngpreds, masked_matrix& invvarmatrixin,
@@ -374,29 +406,7 @@
         //Oct 26, 2009
 
 #if EIGEN
-        if (X.data.cols()== 3){
-            Matrix<double,3,Dynamic> tXW = X.data.transpose()*invvarmatrixin.masked_data->data;
-            Matrix3d xWx = tXW * X.data;
-            Ch = LDLT <MatrixXd>  (xWx );
-            Vector3d beta_3f = Ch.solve(tXW * reg_data.Y.data);
-            sigma2 = (reg_data.Y.data - tXW.transpose() * beta_3f).squaredNorm();
-            beta.data = beta_3f;
-        }
-       else if(X.data.cols()== 2){
-            Matrix<double,2,Dynamic> tXW = X.data.transpose()*invvarmatrixin.masked_data->data;
-            Matrix2d xWx = tXW * X.data;
-            Ch = LDLT <MatrixXd>  (xWx );
-            Vector2d beta_2f = Ch.solve(tXW * reg_data.Y.data);
-            sigma2 = (reg_data.Y.data - tXW.transpose() * beta_2f).squaredNorm();
-            beta.data = beta_2f;
-        }else{
-            // next line is  5997000 flops
-            MatrixXd tXW = X.data.transpose() * invvarmatrixin.masked_data->data;
-            Ch = LDLT <MatrixXd>(tXW * X.data); // 17991 flops
-            beta.data = Ch.solve(tXW * reg_data.Y.data);//5997 flops
-            //next line is: 1000+5000+3000= 9000 flops
-            sigma2 = (reg_data.Y.data - tXW.transpose() * beta.data).squaredNorm();
-        }
+        mmscore_regression(X.data, invvarmatrixin.masked_data->data, Ch);
 #else
         // next line is  5997000 flops
         mematrix<double> tXW = transpose(X) * invvarmatrixin.masked_data;

Modified: branches/ProbABEL-0.50/src/reg1.h
===================================================================
--- branches/ProbABEL-0.50/src/reg1.h	2014-02-27 22:20:11 UTC (rev 1626)
+++ branches/ProbABEL-0.50/src/reg1.h	2014-03-01 00:26:29 UTC (rev 1627)
@@ -99,6 +99,10 @@
     void score(mematrix<double>& resid,
                double tol_chol, int model, int interaction, int ngpreds,
                const masked_matrix& invvarmatrix, int nullmodel = 0);
+
+private:
+    void mmscore_regression(const MatrixXd& X,
+            const MatrixXd& W, LDLT<MatrixXd>& Ch);
 };
 
 class logistic_reg: public base_reg {