[Genabel-commits] r1671 - in branches/ProbABEL-0.50: examples src

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

Author: maartenk
Date: 2014-04-02 22:25:43 +0200 (Wed, 02 Apr 2014)
New Revision: 1671

Modified:
   branches/ProbABEL-0.50/examples/mmscore.R
   branches/ProbABEL-0.50/src/probabel
   branches/ProbABEL-0.50/src/reg1.cpp
Log:
reg1.cpp : simplified code of mmscore(palinear) and fixed failing test_mms.sh check
probabel: introduce check that verifies phenotype file does exists before running pa* (prevents errors while running the script)  
mmscore.R: fixed small typo

Modified: branches/ProbABEL-0.50/examples/mmscore.R
===================================================================
--- branches/ProbABEL-0.50/examples/mmscore.R	2014-04-02 15:57:31 UTC (rev 1670)
+++ branches/ProbABEL-0.50/examples/mmscore.R	2014-04-02 20:25:43 UTC (rev 1671)
@@ -88,5 +88,5 @@
 ## 2) residuals of the phenotype, which will be the new phenotype that
 ## ProbABEL will analyse.
 
-## Mow, go to ProbABEL and start analysis
+## Now, go to ProbABEL and start analysis
 

Modified: branches/ProbABEL-0.50/src/probabel
===================================================================
--- branches/ProbABEL-0.50/src/probabel	2014-04-02 15:57:31 UTC (rev 1670)
+++ branches/ProbABEL-0.50/src/probabel	2014-04-02 20:25:43 UTC (rev 1671)
@@ -169,6 +169,11 @@
 
 
 my $phename = $ARGV[5];
+if (! -e $phename.".PHE"){
+die "Phenotype file $phename.PHE does not exists. The phenotype file should be specified without the .PHE extension.\n";
+}
+
+
 # By default the output file prefix is the same as the name of the
 # phenotype file (minus the .PHE extension and any paths)
 use File::Basename;

Modified: branches/ProbABEL-0.50/src/reg1.cpp
===================================================================
--- branches/ProbABEL-0.50/src/reg1.cpp	2014-04-02 15:57:31 UTC (rev 1670)
+++ branches/ProbABEL-0.50/src/reg1.cpp	2014-04-02 20:25:43 UTC (rev 1671)
@@ -313,35 +313,21 @@
 void linear_reg::mmscore_regression(const mematrix<double>& X,
         const masked_matrix& W_masked, LDLT<MatrixXd>& Ch) {
 
-
     VectorXd Y = reg_data.Y.data.col(0);
-    if (X.data.cols() == 3)
-    {
-        Matrix<double, Dynamic, 3> tXW = W_masked.masked_data->data * X.data;
-        Matrix2d xWx = tXW.transpose() * X.data;
-        Ch = LDLT<MatrixXd>(xWx);
-        Vector3d beta_3f = Ch.solve(tXW.transpose() * Y);
-        sigma2 = (Y - tXW * beta_3f).squaredNorm();
-        beta.data = beta_3f;
-    }
-    else if (X.data.cols() == 2)
-    {
-        Matrix<double,  Dynamic,2> tXW =  W_masked.masked_data->data*X.data;
-        Matrix2d xWx = tXW.transpose() * X.data;
-        Ch = LDLT<MatrixXd> (xWx);
-        Vector2d beta_2f = Ch.solve(tXW.transpose() * Y);
-        sigma2 = (Y - tXW * beta_2f).squaredNorm();
-        beta.data = beta_2f;
-    }
-    else
-    {
-        // next line is  5997000 flops
-        MatrixXd tXW = X.data.transpose() * W_masked.masked_data->data;
-        Ch = LDLT<MatrixXd>(tXW * X.data); // 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();
-    }
+    /*
+     in ProbABEL <0.50 this calculation was performed like t(X)*W
+     This changed to W*X since this is better vectorized since the left hand
+     side has more rows: this introduces an additional transpose, but can be
+     neglected compared to the speedup this brings(about a factor 2 for the
+     palinear with 1 predictor)
+     */
+    MatrixXd tXW = W_masked.masked_data->data * X.data;
+    MatrixXd xWx = tXW.transpose() * X.data;
+    Ch = LDLT<MatrixXd>(xWx);
+    VectorXd beta_vec = Ch.solve(tXW.transpose() * Y);
+    sigma2 = (Y - tXW * beta_vec).squaredNorm();
+    beta.data = beta_vec;
+
 }
 
 void linear_reg::logLikelihood(const mematrix<double>& X) {