[Genabel-commits] r907 - branches/ProbABEL-refactoring/ProbABEL/src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed May 16 17:29:00 CEST 2012
Author: maartenk
Date: 2012-05-16 17:29:00 +0200 (Wed, 16 May 2012)
New Revision: 907
Modified:
branches/ProbABEL-refactoring/ProbABEL/src/main.cpp
branches/ProbABEL-refactoring/ProbABEL/src/reg1.h
Log:
-changes casting style (eg. double -> static_cast<double>)
-removed whitespace
-shorten some lines
Modified: branches/ProbABEL-refactoring/ProbABEL/src/main.cpp
===================================================================
--- branches/ProbABEL-refactoring/ProbABEL/src/main.cpp 2012-05-14 19:16:26 UTC (rev 906)
+++ branches/ProbABEL-refactoring/ProbABEL/src/main.cpp 2012-05-16 15:29:00 UTC (rev 907)
@@ -1,4 +1,4 @@
-//=====================================================================================
+//=============================================================================
//
// Filename: src/main.cpp
//
@@ -11,28 +11,28 @@
// Author: Yurii S. Aulchenko (cox, log, lin regressions)
// Modified by: L.C. Karssen,
// Maksim V. Struchalin
-//
-// modified 14-May-2009 by MVS: interaction with SNP, interaction with SNP with exclusion of interacted covariates,
-// mmscore implemented (poor me)
+//
+// modified 14-May-2009 by MVS: interaction with SNP, interaction with SNP
+// with exclusion of interacted covariates,
+// mmscore implemented (poor me)
// modified 20-Jul-2009 by YSA: small changes, bug fix in mmscore option
// modified 22-Sep-2009 by YSA: "robust" option added
//
// Modified by Han Chen (hanchen at bu.edu) on Nov 9, 2009
-// to extract the covariance between the estimate of beta(SNP) and the estimate of beta(interaction)
-// based on src/main.cpp version 0.1-0 as of Oct 19, 2009
+// to extract the covariance between the estimate of beta(SNP) and the estimate
+// of beta(interaction) based on src/main.cpp version 0.1-0 as of Oct 19, 2009
//
-// Company: Department of Epidemiology, ErasmusMC Rotterdam, The Netherlands.
-// Email: i.aoultchenko at erasmusmc.nl, m.struchalin at erasmusmc.nl
+// Company: Department of Epidemiology, ErasmusMC Rotterdam, The Netherlands.
+// Email: i.aoultchenko at erasmusmc.nl, m.struchalin at erasmusmc.nl
//
-//=====================================================================================
-
+//=============================================================================
+#include <stdio.h>
#include <iostream>
#include <cstdlib>
#include <fstream>
#include <sstream>
#include <string>
#include <iomanip>
-#include <stdio.h>
#include <vector>
#if EIGEN
@@ -57,12 +57,12 @@
if (csnp == 0)
{
fprintf(stdout, "Analysis: %6.2f ...",
- 100. * double(csnp) / double(nsnps));
+ 100. * static_cast<double>(csnp) / static_cast<double>(nsnps));
}
else
{
fprintf(stdout, "\b\b\b\b\b\b\b\b\b\b%6.2f ...",
- 100. * double(csnp) / double(nsnps));
+ 100. * static_cast<double>(csnp) / static_cast<double>(nsnps));
}
std::cout.flush();
}
@@ -90,12 +90,10 @@
exit(1);
}
}
-
}
int create_phenoytype(phedata& phd, cmdvars& input_var)
{
-
phd.set_is_interaction_excluded(input_var.isIsInteractionExcluded());
phd.setphedata(input_var.getPhefilename(), input_var.getNoutcomes(),
input_var.getNpeople(), input_var.getInteraction(),
@@ -173,7 +171,6 @@
//TODO(unknown): compare in create_header_1 and create_header_2 the next lines.
if (input_var.getInteraction() != 0)
{
-
//Han Chen
*outfile[0] << input_var.getSep() << "beta_SNP_A1A2_"
<< phd.model_terms[interaction_cox] << input_var.getSep()
@@ -213,18 +210,16 @@
*outfile[2] << input_var.getSep() << "loglik\n"; //"chi2_SNP_domA1\n";
*outfile[3] << input_var.getSep() << "loglik\n"; //"chi2_SNP_recA1\n";
*outfile[4] << input_var.getSep() << "loglik\n"; //"chi2_SNP_odom\n";
-
}
void create_header2(std::vector<std::ofstream*>& outfile, cmdvars& input_var,
phedata phd, int interaction_cox)
{
-
create_start_of_header(outfile, input_var, phd);
*outfile[0] << input_var.getSep() << "beta_SNP_add" << input_var.getSep()
<< "sebeta_SNP_add";
- //TODO(unknown): compare in create_header_1 and create_header_2 the next lines.
+//TODO(unknown): compare in create_header_1 and create_header_2 the next lines.
if (input_var.getInteraction() != 0)
{
@@ -311,18 +306,18 @@
gendata gtd (input_var.getGenfilename(),nsnps,input_var.getNgpreds(),phd.nids_all,phd.nids,phd.allmeasured,skipd,phd.idnames);
**/
// estimate null model
- //TODO: remove this unused variable if there is not a reason to keep it
+ //TODO(maarten): remove this unused variable if there is not a reason to keep it
//double null_loglik = 0.;
#if COXPH
- coxph_data nrgd=coxph_data(phd,gtd,-1);
+ coxph_data nrgd = coxph_data(phd, gtd, -1);
#else
regdata nrgd = regdata(phd, gtd, -1, input_var.isIsInteractionExcluded());
#endif
std::cout << " loaded null data ...";
#if LOGISTIC
- logistic_reg nrd=logistic_reg(nrgd);
- nrd.estimate(nrgd,0,MAXITER,EPS,CHOLTOL,0, input_var.getInteraction(), input_var.getNgpreds(), invvarmatrix, input_var.getRobust(), 1);
+ logistic_reg nrd = logistic_reg(nrgd);
+ nrd.estimate(nrgd, 0, MAXITER, EPS, CHOLTOL, 0, input_var.getInteraction(), input_var.getNgpreds(), invvarmatrix, input_var.getRobust(), 1);
#elif LINEAR
linear_reg nrd = linear_reg(nrgd);
@@ -334,13 +329,13 @@
#elif COXPH
coxph_reg nrd(nrgd);
- nrd.estimate(nrgd,0,MAXITER,EPS,CHOLTOL,0, input_var.getInteraction(), input_var.getNgpreds(), 1);
+ nrd.estimate(nrgd, 0, MAXITER, EPS, CHOLTOL, 0, input_var.getInteraction(), input_var.getNgpreds(), 1);
#endif
std::cout << " estimated null model ...";
// end null
#if COXPH
- coxph_data rgd(phd,gtd,0);
+ coxph_data rgd(phd, gtd, 0);
#else
regdata rgd(phd, gtd, 0, input_var.isIsInteractionExcluded());
#endif
@@ -355,8 +350,8 @@
std::string outfilename_str(input_var.getOutfilename());
std::vector<std::ofstream*> outfile;
-
- if (input_var.getNgpreds() == 2) //All models output. One file per each model
+ //All models output.One file per each model
+ if (input_var.getNgpreds() == 2)
{
open_files_for_output(outfile, outfilename_str);
if (input_var.getNohead() != 1)
@@ -366,7 +361,6 @@
}
else //Only additive model. Only one output file
{
-
outfile.push_back(
new std::ofstream((outfilename_str + "_add.out.txt").c_str()));
@@ -437,7 +431,6 @@
for (int csnp = 0; csnp < nsnps; csnp++)
{
-
rgd.update_snp(gtd, csnp);
double freq = 0.;
int gcount = 0;
@@ -445,7 +438,7 @@
float snpdata2[gtd.nids];
if (input_var.getNgpreds() == 2)
{
- // freq = ((gtd.G).column_mean(csnp*2)*2.+(gtd.G).column_mean(csnp*2+1))/2.;
+ //freq = ((gtd.G).column_mean(csnp*2)*2.+(gtd.G).column_mean(csnp*2+1))/2.;
gtd.get_var(csnp * 2, snpdata1);
gtd.get_var(csnp * 2 + 1, snpdata2);
for (int ii = 0; ii < gtd.nids; ii++)
@@ -466,15 +459,15 @@
freq += snpdata1[ii] * 0.5;
}
}
- freq /= (double) (gcount);
+ freq /= static_cast<double> (gcount);
int poly = 1;
if (fabs(freq) < 1.e-16 || fabs(1. - freq) < 1.e-16)
poly = 0;
if (fabs(mli.Rsq[csnp]) < 1.e-16)
poly = 0;
-
- if (input_var.getNgpreds() == 2) //All models output. One file per each model
+ //All models output. One file per each model
+ if (input_var.getNgpreds() == 2)
{
//Write mlinfo to output:
for (int file = 0; file < outfile.size(); file++)
@@ -501,9 +494,9 @@
#if LOGISTIC
logistic_reg rd(rgd);
if (input_var.getScore())
- rd.score(nrd.residuals,rgd,0,CHOLTOL,model,input_var.getInteraction(), input_var.getNgpreds(), invvarmatrix);
+ rd.score(nrd.residuals, rgd, 0, CHOLTOL, model, input_var.getInteraction(), input_var.getNgpreds(), invvarmatrix);
else
- rd.estimate(rgd,0,MAXITER,EPS,CHOLTOL,model,input_var.getInteraction(), input_var.getNgpreds(), invvarmatrix, input_var.getRobust());
+ rd.estimate(rgd, 0, MAXITER, EPS, CHOLTOL, model,input_var.getInteraction(), input_var.getNgpreds(), invvarmatrix, input_var.getRobust());
#elif LINEAR
linear_reg rd(rgd);
if (input_var.getScore())
@@ -520,7 +513,7 @@
}
#elif COXPH
coxph_reg rd(rgd);
- rd.estimate(rgd,0,MAXITER,EPS,CHOLTOL,model,input_var.getInteraction(), true, input_var.getNgpreds());
+ rd.estimate(rgd, 0, MAXITER, EPS, CHOLTOL, model, input_var.getInteraction(), true, input_var.getNgpreds());
#endif
if (!input_var.getAllcov() && model == 0
@@ -584,7 +577,6 @@
*chi2[model] << "nan";
}
//________________________________
-
}
else //beta, sebeta = nan
{
@@ -687,7 +679,6 @@
#endif
*outfile[4] << chi2[4]->str() << "\n";
//Oct 26, 2009
-
}
else //Only additive model. Only one output file
{
@@ -709,9 +700,9 @@
#if LOGISTIC
logistic_reg rd(rgd);
if (input_var.getScore())
- rd.score(nrd.residuals,rgd,0,CHOLTOL,model, input_var.getInteraction(), input_var.getNgpreds(), invvarmatrix);
+ rd.score(nrd.residuals, rgd, 0, CHOLTOL, model, input_var.getInteraction(), input_var.getNgpreds(), invvarmatrix);
else
- rd.estimate(rgd,0,MAXITER,EPS,CHOLTOL,model, input_var.getInteraction(), input_var.getNgpreds(), invvarmatrix, input_var.getRobust());
+ rd.estimate( rgd, 0, MAXITER, EPS, CHOLTOL, model, input_var.getInteraction(), input_var.getNgpreds(), invvarmatrix, input_var.getRobust());
#elif LINEAR
//cout << (rgd.get_unmasked_data()).nids << " 1\n";
#if DEBUG
@@ -797,7 +788,7 @@
}
#elif COXPH
coxph_reg rd(rgd);
- rd.estimate(rgd,0,MAXITER,EPS,CHOLTOL,model, input_var.getInteraction(), true, input_var.getNgpreds());
+ rd.estimate(rgd, 0, MAXITER, EPS, CHOLTOL, model, input_var.getInteraction(), true, input_var.getNgpreds());
#endif
if (!input_var.getAllcov() && input_var.getInteraction() == 0)
@@ -812,7 +803,6 @@
else
{
start_pos = 0;
-
}
#if DEBUG
cout << " start_pos;" << start_pos << "\n";
@@ -912,8 +902,8 @@
*outfile[0] << beta_sebeta[0]->str() << "\n";
#if DEBUG
cout << "Se beta" << beta_sebeta[0] << "\n";
-#endif
- }
+#endif
+ }
}
//clean chi2
for (int i = 0; i < 5; i++)
@@ -925,7 +915,6 @@
chi2[i]->str("");
}
update_progress_to_cmd_line(csnp, nsnps);
-
}
fprintf(stdout, "\b\b\b\b\b\b\b\b\b\b%6.2f", 100.);
Modified: branches/ProbABEL-refactoring/ProbABEL/src/reg1.h
===================================================================
--- branches/ProbABEL-refactoring/ProbABEL/src/reg1.h 2012-05-14 19:16:26 UTC (rev 906)
+++ branches/ProbABEL-refactoring/ProbABEL/src/reg1.h 2012-05-16 15:29:00 UTC (rev 907)
@@ -50,7 +50,7 @@
std::cout << "apply_model():interaction" << interaction << std::endl;
std::cout << "apply_model():ngpreds" << ngpreds << std::endl;
std::cout << "apply_model():is_interaction_excluded"
- << is_interaction_excluded << std::endl;
+ << is_interaction_excluded << std::endl;
std::cout << "apply_model():iscox" << iscox << std::endl;
std::cout << "apply_model():nullmodel" << nullmodel << std::endl;
#endif
@@ -74,17 +74,19 @@
{
if (iscox)
{
+ //Maksim: interaction with SNP;;
nX[i * nX.ncol + c1] = X[i * X.ncol + csnp_p1]
- * X[i * X.ncol + interaction - 1]; //Maksim: interaction with SNP;;
+ * X[i * X.ncol + interaction - 1];
nX[i * nX.ncol + c2] = X[i * X.ncol + csnp_p2]
- * X[i * X.ncol + interaction - 1]; //Maksim: interaction with SNP;;
+ * X[i * X.ncol + interaction - 1];
}
else
{
+ //Maksim: interaction with SNP;;
nX[i * nX.ncol + c1] = X[i * X.ncol + csnp_p1]
- * X[i * X.ncol + interaction]; //Maksim: interaction with SNP;;
+ * X[i * X.ncol + interaction];
nX[i * nX.ncol + c2] = X[i * X.ncol + csnp_p2]
- * X[i * X.ncol + interaction]; //Maksim: interaction with SNP;;
+ * X[i * X.ncol + interaction];
}
}
//________________________
@@ -187,7 +189,6 @@
{
return (X);
}
-
}
mematrix<double> nX;
if (interaction != 0)
@@ -242,7 +243,6 @@
nX_without_interact_phe[row * nX_without_interact_phe.ncol
+ col_new] = nX[row * nX.ncol + col];
}
-
}
}
return nX_without_interact_phe;
@@ -301,92 +301,7 @@
// delete sebeta;
// delete residuals;
}
- // void mmscore(regdata &rdata,int verbose, double tol_chol, int model, int interaction, int ngpreds, mematrix<double> & invvarmatrix, int nullmodel=0)
- // {
- // int length_beta = (rdata.X).ncol;
- // beta.reinit(length_beta,1);
- // sebeta.reinit(length_beta,1);
- // mematrix<double> G;
- // mematrix<double> Y = rdata.Y;
- //
- // G.reinit((rdata.X).nrow,1);
- // for (int i=0;i<(rdata.X).nrow;i++) G[i] = (rdata.X).get(i,1);
- //// double sumgt = column_sum(G).get(0,0);
- // double meangt = G.column_mean(0);
- // double meany = Y.column_mean(0);
- //
- // G = G - meangt;
- // Y = Y - meany;
- // double U, V;
- // mematrix<double> U_mema, V_mema;
- //
- //
- // U_mema = transpose(G)*invvarmatrix;
- // U = (U_mema * Y).get(0, 0);
- //
- // V_mema = transpose(G)*invvarmatrix;
- // V = (V_mema*G).get(0, 0);
- //
- //
- // chi2_score = U*U/V;
- //
- // beta.put(-100, 0,0);
- // double bb = U/V;
- // beta.put(bb,1,0);
- // sebeta.put(-100.,0,0);
- // double sebe=sqrt(1./V);
- // sebeta.put(sebe,1,0);
- //
- //
- //
- //
- // /*
- // for(snp in snps)
- // {
- // print(snp)
- // ns <- which(mm$snpnames==snp)
- // X <- gt[,snp]
- // V <- h2ht$InvSigma #*var(h2ht$residualY,na.rm=T)
- // UU <- ((X-mean(X))%*%V%*%(Y-mean(Y)))
- // VV <- ((X-mean(X))%*%V%*%(X-mean(X)))
- // beta <- UU/VV
- // se <- sqrt(1./VV)
- // chi2 <- UU*UU/VV
- // print(c(re$chi2.1df[ns],qt$chi2.1df[ns],gr$chi2.1df[ns],mm$chi2.1df[ns],(beta/se)^2))
- // print(c(re$effB[ns],qt$effB[ns],gr$effB[ns],mm$effB[ns],beta))
- // }
- //
- // */
- //
- //
- //
- ///*
- // int length_beta = (rdata.X).ncol;
- // beta.reinit(length_beta,1);
- // sebeta.reinit(length_beta,1);
- //
- // mematrix<double> G;
- // G.reinit((rdata.X).nrow,1);
- // for (int i=0;i<(rdata.X).nrow;i++) G[i] = (rdata.X).get(i,1);
- // double sumgt = column_sum(G).get(0,0);
- // double meangt = G.column_mean(0);
- // G = G - meangt;
- // mematrix<double> YY = rdata.Y;
- // YY = YY - YY.column_mean(0);
- // mematrix<double> tG = transpose(G);
- // tG = tG*invvarmatrix;
- // double U = (tG*YY).get(0,0);
- // double V = (tG*G).get(0,0);
- //
- // beta.put(-100, 0,0);
- // double bb = (U/V)*var(YY);
- // beta.put(bb,1,0);
- // sebeta.put(-100.,0,0);
- // sebeta.put(abs(bb)/sqrt(U*U/V),1,0);
- //
- //*/
- //
- // }
+
void estimate(regdata &rdatain, int verbose, double tol_chol, int model,
int interaction, int ngpreds, mematrix<double> invvarmatrixin,
int robust, int nullmodel = 0)
@@ -398,7 +313,6 @@
{
cout << rdata.is_interaction_excluded
<< " <-irdata.is_interaction_excluded\n";
-
}
mematrix<double> invvarmatrix;
if (invvarmatrixin.nrow != 0 && invvarmatrixin.ncol != 0)
@@ -417,7 +331,6 @@
}
i1++;
}
-
}
if (verbose)
{
@@ -457,13 +370,12 @@
}
//Oct 26, 2009
mematrix<double> tX = transpose(X);
-
if (invvarmatrix.nrow != 0 && invvarmatrix.ncol != 0)
{
tX = tX * invvarmatrix;
//!check if quicker
//tX = productXbySymM(tX,invvarmatrix);
- // X = invvarmatrix*X; std::cout<<"new tX.nrow="<<X.nrow<<" tX.ncol="<<X.ncol<<"\n";
+ // = invvarmatrix*X; std::cout<<"new tX.nrow="<<X.nrow<<" tX.ncol="<<X.ncol<<"\n";
}
if (verbose)
{
@@ -538,7 +450,7 @@
// printf("sigma2+= = %f\n", sigma2);
}
- double sigma2_internal = sigma2 / (N - double(length_beta));
+ double sigma2_internal = sigma2 / (N - static_cast<double>(length_beta));
// now compute residual variance
// sigma2 = 0.;
@@ -591,7 +503,7 @@
residuals[i] = resid;
loglik -= halfrecsig2 * resid * resid;
}
- loglik -= double(rdata.nids) * log(sqrt(sigma2));
+ loglik -= static_cast<double>(rdata.nids) * log(sqrt(sigma2));
//cout << "estimate " << rdata.nids << "\n";
//
@@ -697,7 +609,7 @@
rdata.is_interaction_excluded, false, nullmodel);
beta.reinit(X.ncol, 1);
sebeta.reinit(X.ncol, 1);
- double N = double(resid.nrow);
+ double N = static_cast<double>(resid.nrow);
mematrix<double> tX = transpose(X);
@@ -734,7 +646,6 @@
chi2 = chi2 * (1. / sigma2_internal);
chi2_score = chi2[0];
}
-
};
class logistic_reg
@@ -804,7 +715,6 @@
}
i1++;
}
-
}
mematrix<double> X = apply_model(rdata.X, model, interaction, ngpreds,
@@ -857,7 +767,7 @@
fprintf(stdout,"tXInv %f %f %f\n",tX.get(0,0),tX.get(1,0),tX.get(2,0));
if (X.ncol==4) fprintf(stdout,"X[4] %f\n",tX.get(3,0));
*/
- //TODO: remove this unused variable if there is not a reason to keep it
+ //TODO(maarten): remove this unused variable if there is not a reason to keep it
//double N;
niter = 0;
double delta = 1.;
@@ -1025,7 +935,7 @@
rdata.is_interaction_excluded, false, nullmodel);
beta.reinit(X.ncol, 1);
sebeta.reinit(X.ncol, 1);
- double N = double(resid.nrow);
+ double N = static_cast<double>(resid.nrow);
mematrix<double> tX = transpose(X);
if (invvarmatrix.nrow != 0 && invvarmatrix.ncol != 0)
More information about the Genabel-commits
mailing list