[Genabel-commits] r1537 - branches/ProbABEL-0.50/src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Jan 10 23:56:32 CET 2014
Author: maartenk
Date: 2014-01-10 23:56:32 +0100 (Fri, 10 Jan 2014)
New Revision: 1537
Modified:
branches/ProbABEL-0.50/src/reg1.cpp
Log:
removed double copied code (1line) and remove verbose ,dead code, and outcommented lines
Modified: branches/ProbABEL-0.50/src/reg1.cpp
===================================================================
--- branches/ProbABEL-0.50/src/reg1.cpp 2014-01-10 22:41:44 UTC (rev 1536)
+++ branches/ProbABEL-0.50/src/reg1.cpp 2014-01-10 22:56:32 UTC (rev 1537)
@@ -313,27 +313,13 @@
invvarmatrixin.update_mask(rdatain.masked_data);
// invvarmatrixin.masked_data->print();
}
- if (verbose)
- {
- cout << rdata.is_interaction_excluded
- << " <-rdata.is_interaction_excluded\n";
- // std::cout << "invvarmatrix:\n";
- // invvarmatrixin.masked_data->print();
- std::cout << "rdata.X:\n";
- rdata.X.print();
- }
+
mematrix<double> X = apply_model(rdata.X, model, interaction, ngpreds,
rdata.is_interaction_excluded, false,
nullmodel);
- if (verbose)
- {
- std::cout << "X:\n";
- X.print();
- std::cout << "Y:\n";
- rdata.Y.print();
- }
+
int length_beta = X.ncol;
beta.reinit(length_beta, 1);
sebeta.reinit(length_beta, 1);
@@ -359,10 +345,9 @@
//X.data=X.data.transpose()*invvarmatrixin.masked_data->data;
X.data =MatrixXd(X.data.cols(),X.data.cols()).setZero().selfadjointView<Lower>().rankUpdate(X.data.adjoint());
- X.data =MatrixXd(X.data.cols(),X.data.cols()).setZero().selfadjointView<Lower>().rankUpdate(X.data.adjoint());
- X.nrow=X.data.rows();
- X.ncol=X.data.cols();
- X.nelements=X.ncol*X.nrow;
+ X.nrow=X.data.rows();
+ X.ncol=X.data.cols();
+ X.nelements=X.ncol*X.nrow;
#else
tX = transpose(X) * invvarmatrixin.masked_data;
#endif
@@ -394,35 +379,7 @@
#endif
double N = X.nrow;
-#if EIGEN_COMMENTEDOUT
- MatrixXd Xeigen = X.data;
- MatrixXd tXeigen = Xeigen.transpose();
- MatrixXd tXXeigen = tXeigen * Xeigen;
- VectorXd Yeigen = rdata.Y.data;
- VectorXd tXYeigen = tXeigen * Yeigen;
- // Solve X^T * X * beta = X^T * Y for beta:
- VectorXd betaeigen = tXXeigen.fullPivLu().solve(tXYeigen);
- beta.data = betaeigen;
- if (verbose)
- {
- std::cout << setprecision(9) << "Xeigen:\n" << Xeigen << endl;
- std::cout << setprecision(9) << "tX:\n" << tXeigen << endl;
- std::cout << setprecision(9) << "tXX:\n" << tXXeigen << endl;
- std::cout << setprecision(9) << "tXY:\n" << tXYeigen << endl;
- std::cout << setprecision(9) << "beta:\n"<< betaeigen << endl;
- printf("----\n");
- printf("beta[0] = %e\n", betaeigen.data()[0]);
- printf("----\n");
-// (beta).print();
- double relative_error = (tXXeigen * betaeigen - tXYeigen).norm() /
- tXYeigen.norm(); // norm() is L2 norm
- cout << "The relative error is:\n" << relative_error << endl;
- }
-
- // This one is needed later on in this function
- mematrix<double> tXX_i = invert(tXX);
-#else
//
// use cholesky to invert
//
@@ -442,22 +399,6 @@
beta = tXX * tXY;
#endif
- if (verbose)
- {
- std::cout << "tX:\n";
- //tX.print();
- std::cout << "tXX:\n";
- //tXX.print();
- std::cout << "chole tXX:\n";
- tXX_i.print();
- std::cout << "tXX-1:\n";
- tXX_i.print();
- std::cout << "tXY:\n";
- //tXY.print();
- std::cout << "beta:\n";
- (beta).print();
- }
-#endif
// now compute residual variance
sigma2 = 0.;
#if EIGEN
@@ -469,64 +410,22 @@
mematrix<double> ttX = transpose(tX);
mematrix<double> sigma2_matrix = rdata.Y;
mematrix<double> sigma2_matrix1 = ttX * beta;
- // std::cout << "sigma2_matrix\n";
- // sigma2_matrix.print();
- //
- // std::cout << "sigma2_matrix1\n";
- // sigma2_matrix1.print();
sigma2_matrix = sigma2_matrix - sigma2_matrix1;
- // std::cout << "sigma2_matrix\n";
- // sigma2_matrix.print();
- // 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++)
{
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
- // sigma2 = 0.;
- // for (int i =0;i<(rdata.Y).nrow;i++)
- // sigma2 += ((rdata.Y).get(i,0))*((rdata.Y).get(i,0));
- // for (int i=0;i<length_beta;i++)
- // sigma2 -= 2. * (beta.get(i,0)) * tXY.get(i,0);
- // for (int i=0;i<(length_beta);i++)
- // for (int j=0;j<(length_beta);j++)
- // sigma2 += (beta.get(i,0)) * (beta.get(j,0)) * tXX.get(i,j);
- // std::cout<<"sigma2="<<sigma2<<"\n";
- // std::cout<<"sigma2_internal="<<sigma2_internal<<"\n";
- // replaced for ML
- // sigma2_internal = sigma2/(N - double(length_beta) - 1);
- // std::cout << "sigma2/=N = "<< sigma2 << "\n";
+
sigma2 /= N;
- // std::cout<<"N="<<N<<", length_beta="<<length_beta<<"\n";
- if (verbose)
- {
- std::cout << "sigma2 = " << sigma2 << "\n";
- }
- /*
- loglik = 0.;
- double ss=0;
- for (int i=0;i<rdata.nids;i++) {
- double resid = rdata.Y[i] - beta.get(0,0); // intercept
- for (int j=1;j<beta.nrow;j++) resid -= beta.get(j,0)*X.get(i,j);
- // residuals[i] = resid;
- ss += resid*resid;
- }
- sigma2 = ss/N;
- */
- //cout << "estimate " << rdata.nids << "\n";
- //(rdata.X).print();
- //for (int i=0;i<rdata.nids;i++) cout << rdata.masked_data[i] << " ";
- //cout << endl;
+
loglik = 0.;
double halfrecsig2 = .5 / sigma2;
#if EIGEN
@@ -550,13 +449,11 @@
}
#endif
loglik -= static_cast<double>(rdata.nids) * log(sqrt(sigma2));
- // cout << "estimate " << rdata.nids << "\n";
//
// Ugly fix to the fact that if we do mmscore, sigma2 is already
// in the matrix...
// YSA, 2009.07.20
//
- //cout << "estimate 0\n";
if (invvarmatrixin.length_of_mask != 0)
sigma2_internal = 1.0;
More information about the Genabel-commits
mailing list