From a5598b09d5022ec7b9f75a12098068ca40bc4d58 Mon Sep 17 00:00:00 2001 From: Peter Carbonetto Date: Thu, 29 Jun 2017 08:30:55 -0500 Subject: Remove FORCE_FLOAT from vc.cpp. --- src/logistic.cpp | 389 ++++++++--------- src/vc.cpp | 1237 +++++++++++++++++++++++++++--------------------------- 2 files changed, 778 insertions(+), 848 deletions(-) (limited to 'src') diff --git a/src/logistic.cpp b/src/logistic.cpp index 3f9d6ff..c1ddac1 100644 --- a/src/logistic.cpp +++ b/src/logistic.cpp @@ -57,15 +57,14 @@ double fLogit_mixed(gsl_vector *beta, return -total; } - -void logistic_mixed_pred(gsl_vector *beta // Vector of parameters +void logistic_mixed_pred(gsl_vector *beta, // Vector of parameters // length = 1 + Sum_k(C_k -1) - ,gsl_matrix_int *X //Matrix Nobs x K - ,gsl_vector_int *nlev // Vector with number categories - ,gsl_matrix *Xc // continuous covariates Matrix Nobs x Kc (NULL if not used) - ,gsl_vector *yhat //Vector of prob. predicted by the logistic - ) -{ + gsl_matrix_int *X, // Matrix Nobs x K + gsl_vector_int *nlev, // Vector with number categories + gsl_matrix *Xc, // Continuous covariates matrix: + // obs x Kc (NULL if not used). + gsl_vector *yhat){ // Vector of prob. predicted by + // the logistic for(int i = 0; i < X->size1; ++i) { double Xbetai=beta->data[0]; int iParm=1; @@ -74,7 +73,7 @@ void logistic_mixed_pred(gsl_vector *beta // Vector of parameters Xbetai+=beta->data[gsl_matrix_int_get(X,i,k)-1+iParm]; iParm+=nlev->data[k]-1; } - // Adding the continuous + // Adding the continuous. for(int k = 0; k < (Xc->size2); ++k) Xbetai+= gsl_matrix_get(Xc,i,k)*beta->data[iParm++]; yhat->data[i]=1/(1 + gsl_sf_exp(-Xbetai)); @@ -82,11 +81,9 @@ void logistic_mixed_pred(gsl_vector *beta // Vector of parameters } -/* The gradient of f, df = (df/dx, df/dy). */ -void -wgsl_mixed_optim_df (const gsl_vector *beta, void *params, - gsl_vector *out) -{ +// The gradient of f, df = (df/dx, df/dy). +void wgsl_mixed_optim_df (const gsl_vector *beta, void *params, + gsl_vector *out) { fix_parm_mixed_T *p = (fix_parm_mixed_T *)params; int n = p->y->size; int K = p->X->size2; @@ -95,7 +92,8 @@ wgsl_mixed_optim_df (const gsl_vector *beta, void *params, // Intitialize gradient out necessary? for(int i = 0; i < npar; ++i) - out->data[i]= 0; + out->data[i]= 0; + // Changed loop start at 1 instead of 0 to avoid regularization of beta 0. for(int i = 1; i < npar; ++i) out->data[i]= p->lambdaL2*beta->data[i]; @@ -142,14 +140,17 @@ void wgsl_mixed_optim_hessian (const gsl_vector *beta, void *params, int K = p->X->size2; int Kc = p->Xc->size2; int npar = beta->size; - gsl_vector *gn = gsl_vector_alloc(npar); /* gn */ + gsl_vector *gn = gsl_vector_alloc(npar); // gn + // Intitialize Hessian out necessary ??? gsl_matrix_set_zero(out); + /* Changed loop start at 1 instead of 0 to avoid regularization of beta 0*/ for(int i = 1; i < npar; ++i) - gsl_matrix_set(out,i,i,(p->lambdaL2)); // Double check this - // L1 penalty not working yet, as not differentiable, I may need to do coordinate descent (as in glm_net) + gsl_matrix_set(out,i,i,(p->lambdaL2)); // Double-check this. + // L1 penalty not working yet, as not differentiable, I may need to + // do coordinate descent (as in glm_net) for(int i = 0; i < n; ++i) { double pn=0; double aux=0; @@ -160,15 +161,17 @@ void wgsl_mixed_optim_hessian (const gsl_vector *beta, void *params, Xbetai+=beta->data[gsl_matrix_int_get(p->X,i,k)-1+iParm1]; iParm1+=p->nlev->data[k]-1; //-1? } - // Adding the continuous + + // Adding the continuous. for(int k = 0; k < Kc; ++k) Xbetai+= gsl_matrix_get(p->Xc,i,k)*beta->data[iParm1++]; pn= 1/(1 + gsl_sf_exp(-Xbetai)); + // Add a protection for pn very close to 0 or 1? aux=pn*(1-pn); - // Calculate sub-gradient vector gn + // Calculate sub-gradient vector gn. gsl_vector_set_zero(gn); gn->data[0]= 1; iParm1=1; @@ -177,7 +180,8 @@ void wgsl_mixed_optim_hessian (const gsl_vector *beta, void *params, gn->data[gsl_matrix_int_get(p->X,i,k)-1+iParm1]=1; iParm1+=p->nlev->data[k]-1; } - // Adding the continuous + + // Adding the continuous. for(int k = 0; k < Kc; ++k) { gn->data[iParm1++] = gsl_matrix_get(p->Xc,i,k); } @@ -191,42 +195,31 @@ void wgsl_mixed_optim_hessian (const gsl_vector *beta, void *params, gsl_vector_free(gn); } - -double wgsl_mixed_optim_f(gsl_vector *v, void *params) -{ +double wgsl_mixed_optim_f(gsl_vector *v, void *params) { double mLogLik=0; fix_parm_mixed_T *p = (fix_parm_mixed_T *)params; mLogLik = fLogit_mixed(v,p->X,p->nlev,p->Xc,p->y,p->lambdaL1,p->lambdaL2); return mLogLik; } - -/* Compute both f and df together. */ +// Compute both f and df together. void -wgsl_mixed_optim_fdf (gsl_vector *x, void *params, - double *f, gsl_vector *df) -{ +wgsl_mixed_optim_fdf (gsl_vector *x, void *params, double *f, gsl_vector *df) { *f = wgsl_mixed_optim_f(x, params); wgsl_mixed_optim_df(x, params, df); } - -int logistic_mixed_fit(gsl_vector *beta - ,gsl_matrix_int *X - ,gsl_vector_int *nlev - ,gsl_matrix *Xc // continuous covariates Matrix Nobs x Kc (NULL if not used) - ,gsl_vector *y - ,double lambdaL1 - ,double lambdaL2) -{ - +// Xc is the matrix of continuous covariates, Nobs x Kc (NULL if not used). +int logistic_mixed_fit(gsl_vector *beta, gsl_matrix_int *X, + gsl_vector_int *nlev, gsl_matrix *Xc, + gsl_vector *y, double lambdaL1, double lambdaL2) { double mLogLik=0; fix_parm_mixed_T p; int npar = beta->size; int iter=0; double maxchange=0; - //Intializing fix parameters + // Intializing fix parameters. p.X=X; p.Xc=Xc; p.nlev=nlev; @@ -234,49 +227,35 @@ int logistic_mixed_fit(gsl_vector *beta p.lambdaL1=lambdaL1; p.lambdaL2=lambdaL2; - //Initial fit - //#ifdef _RPR_DEBUG_ + // Initial fit. mLogLik = wgsl_mixed_optim_f(beta,&p); - //fprintf(stderr,"#Initial -log(Lik(0))=%lf\n",mLogLik); - //#endif //_RPR_DEBUG - gsl_matrix *myH = gsl_matrix_alloc(npar,npar); /* Hessian matrix*/ - gsl_vector *stBeta = gsl_vector_alloc(npar); /* Direction to move */ + gsl_matrix *myH = gsl_matrix_alloc(npar,npar); // Hessian matrix. + gsl_vector *stBeta = gsl_vector_alloc(npar); // Direction to move. - gsl_vector *myG = gsl_vector_alloc(npar); /* Gradient*/ - gsl_vector *tau = gsl_vector_alloc(npar); /* tau for QR*/ + gsl_vector *myG = gsl_vector_alloc(npar); // Gradient. + gsl_vector *tau = gsl_vector_alloc(npar); // tau for QR. for(iter=0;iter<100;iter++){ - wgsl_mixed_optim_hessian(beta,&p,myH); //Calculate Hessian - wgsl_mixed_optim_df(beta,&p,myG); //Calculate Gradient - gsl_linalg_QR_decomp(myH,tau); //Calculate next beta + wgsl_mixed_optim_hessian(beta,&p,myH); // Calculate Hessian. + wgsl_mixed_optim_df(beta,&p,myG); // Calculate Gradient. + gsl_linalg_QR_decomp(myH,tau); // Calculate next beta. gsl_linalg_QR_solve(myH,tau,myG,stBeta); gsl_vector_sub(beta,stBeta); - //Monitor convergence + // Monitor convergence. maxchange=0; for(int i=0;idata[i])) maxchange=fabs(stBeta->data[i]); -#ifdef _RPR_DEBUG_ - //mLogLik = wgsl_mixed_optim_f(beta,&p); - //fprintf(stderr,"#iter %d, -log(Lik(0))=%lf,%lf\n",(int)iter,mLogLik,maxchange); -#endif //_RPR_DEBUG - if(maxchange<1E-4) break; } -#ifdef _RPR_DEBUG_ - //for (int i = 0; i < npar; i++) - // fprintf(stderr,"#par_%d= %lf\n",i,beta->data[i]); -#endif //_RPR_DEBUG - - //Final fit + // Final fit. mLogLik = wgsl_mixed_optim_f(beta,&p); - //fprintf(stderr,"#Final %d) -log(Lik(0))=%lf, maxchange %lf\n",iter,mLogLik,maxchange); - + gsl_vector_free (tau); gsl_vector_free (stBeta); gsl_vector_free (myG); @@ -291,7 +270,7 @@ int logistic_mixed_fit(gsl_vector *beta // I need to bundle all the data that goes to the function to optimze // together. -typedef struct{ +typedef struct { gsl_matrix_int *X; gsl_vector_int *nlev; gsl_vector *y; @@ -299,16 +278,13 @@ typedef struct{ double lambdaL2; } fix_parm_cat_T; -double fLogit_cat(gsl_vector *beta, - gsl_matrix_int *X, - gsl_vector_int *nlev, - gsl_vector *y, - double lambdaL1, - double lambdaL2) { +double fLogit_cat (gsl_vector *beta, gsl_matrix_int *X, gsl_vector_int *nlev, + gsl_vector *y, double lambdaL1, double lambdaL2) { int n = y->size; int npar = beta->size; double total = 0; double aux = 0; + // omp_set_num_threads(ompthr); /\* Changed loop start at 1 instead // of 0 to avoid regularization of beta 0*\/ /\*#pragma omp parallel // for reduction (+:total)*\/ @@ -336,12 +312,12 @@ double fLogit_cat(gsl_vector *beta, return -total; } -void logistic_cat_pred(gsl_vector *beta // Vector of parameters length = 1 + Sum_k(C_k - 1) - ,gsl_matrix_int *X //Matrix Nobs x K - ,gsl_vector_int *nlev // Vector with number categories - ,gsl_vector *yhat //Vector of prob. predicted by the logistic - ) -{ +void logistic_cat_pred (gsl_vector *beta, // Vector of parameters + // length = 1 + Sum_k(C_k-1). + gsl_matrix_int *X, // Matrix Nobs x K + gsl_vector_int *nlev, // Vector with #categories + gsl_vector *yhat){ // Vector of prob. predicted by + // the logistic. for(int i = 0; i < X->size1; ++i) { double Xbetai=beta->data[0]; int iParm=1; @@ -354,20 +330,19 @@ void logistic_cat_pred(gsl_vector *beta // Vector of parameters length = 1 + Su } } - -/* The gradient of f, df = (df/dx, df/dy). */ -void -wgsl_cat_optim_df (const gsl_vector *beta, void *params, - gsl_vector *out) -{ +// The gradient of f, df = (df/dx, df/dy). +void wgsl_cat_optim_df (const gsl_vector *beta, void *params, + gsl_vector *out) { fix_parm_cat_T *p = (fix_parm_cat_T *)params; int n = p->y->size; int K = p->X->size2; - int npar = beta->size; + int npar = beta->size; + // Intitialize gradient out necessary? for(int i = 0; i < npar; ++i) - out->data[i]= 0; - /* Changed loop start at 1 instead of 0 to avoid regularization of beta 0 */ + out->data[i]= 0; + + // Changed loop start at 1 instead of 0 to avoid regularization of beta 0. for(int i = 1; i < npar; ++i) out->data[i]= p->lambdaL2*beta->data[i]; for(int i = 1; i < npar; ++i) @@ -382,7 +357,7 @@ wgsl_cat_optim_df (const gsl_vector *beta, void *params, Xbetai+=beta->data[gsl_matrix_int_get(p->X,i,k)-1+iParm]; iParm+=p->nlev->data[k]-1; } - // total += y->data[i]*Xbetai-log(1+gsl_sf_exp(Xbetai)); + pn= -( p->y->data[i] - 1/(1 + gsl_sf_exp(-Xbetai)) ); out->data[0]+= pn; @@ -395,24 +370,23 @@ wgsl_cat_optim_df (const gsl_vector *beta, void *params, } } - -/* The Hessian of f */ -void -wgsl_cat_optim_hessian (const gsl_vector *beta, void *params, - gsl_matrix *out) -{ +// The Hessian of f. +void wgsl_cat_optim_hessian (const gsl_vector *beta, void *params, + gsl_matrix *out) { fix_parm_cat_T *p = (fix_parm_cat_T *)params; int n = p->y->size; int K = p->X->size2; - int npar = beta->size; - // Intitialize Hessian out necessary ??? + int npar = beta->size; + + // Intitialize Hessian out necessary. gsl_matrix_set_zero(out); - /* Changed loop start at 1 instead of 0 to avoid regularization of beta 0*/ + + // Changed loop start at 1 instead of 0 to avoid regularization of beta. for(int i = 1; i < npar; ++i) - gsl_matrix_set(out,i,i,(p->lambdaL2)); // Double check this - // L1 penalty not working yet, as not differentiable, I may need to do coordinate descent (as in glm_net) - + gsl_matrix_set(out,i,i,(p->lambdaL2)); // Double-check this. + // L1 penalty not working yet, as not differentiable, I may need to + // do coordinate descent (as in glm_net). for(int i = 0; i < n; ++i) { double pn=0; double aux=0; @@ -424,8 +398,9 @@ wgsl_cat_optim_hessian (const gsl_vector *beta, void *params, Xbetai+=beta->data[gsl_matrix_int_get(p->X,i,k)-1+iParm1]; iParm1+=p->nlev->data[k]-1; //-1? } - // total += y->data[i]*Xbetai-log(1+gsl_sf_exp(Xbetai)); + pn= 1/(1 + gsl_sf_exp(-Xbetai)); + // Add a protection for pn very close to 0 or 1? aux=pn*(1-pn); *gsl_matrix_ptr(out,0,0)+=aux; @@ -441,7 +416,8 @@ wgsl_cat_optim_hessian (const gsl_vector *beta, void *params, *gsl_matrix_ptr(out,gsl_matrix_int_get(p->X,i,k1)-1+iParm1,0)+=aux; iParm2=1; for(int k2 = 0; k2 < K; ++k2) { - if((gsl_matrix_int_get(p->X,i,k1)>0) && (gsl_matrix_int_get(p->X,i,k2)>0)) + if((gsl_matrix_int_get(p->X,i,k1)>0) && + (gsl_matrix_int_get(p->X,i,k2)>0)) *gsl_matrix_ptr(out ,gsl_matrix_int_get(p->X,i,k1)-1+iParm1 ,gsl_matrix_int_get(p->X,i,k2)-1+iParm2 @@ -453,66 +429,56 @@ wgsl_cat_optim_hessian (const gsl_vector *beta, void *params, } } - -double wgsl_cat_optim_f(gsl_vector *v, void *params) -{ +double wgsl_cat_optim_f(gsl_vector *v, void *params) { double mLogLik=0; fix_parm_cat_T *p = (fix_parm_cat_T *)params; mLogLik = fLogit_cat(v,p->X,p->nlev,p->y,p->lambdaL1,p->lambdaL2); return mLogLik; } - -/* Compute both f and df together. */ -void -wgsl_cat_optim_fdf (gsl_vector *x, void *params, - double *f, gsl_vector *df) -{ +// Compute both f and df together. +void wgsl_cat_optim_fdf (gsl_vector *x, void *params, double *f, + gsl_vector *df) { *f = wgsl_cat_optim_f(x, params); wgsl_cat_optim_df(x, params, df); } - -int logistic_cat_fit(gsl_vector *beta - ,gsl_matrix_int *X - ,gsl_vector_int *nlev - ,gsl_vector *y - ,double lambdaL1 - ,double lambdaL2) -{ +int logistic_cat_fit(gsl_vector *beta, + gsl_matrix_int *X, + gsl_vector_int *nlev, + gsl_vector *y, + double lambdaL1, + double lambdaL2) { double mLogLik=0; fix_parm_cat_T p; int npar = beta->size; int iter=0; double maxchange=0; - //Intializing fix parameters + // Intializing fix parameters. p.X=X; p.nlev=nlev; p.y=y; p.lambdaL1=lambdaL1; p.lambdaL2=lambdaL2; - //Initial fit - //#ifdef _RPR_DEBUG_ + // Initial fit. mLogLik = wgsl_cat_optim_f(beta,&p); - //fprintf(stderr,"#Initial -log(Lik(0))=%lf\n",mLogLik); - //#endif //_RPR_DEBUG - gsl_matrix *myH = gsl_matrix_alloc(npar,npar); /* Hessian matrix*/ - gsl_vector *stBeta = gsl_vector_alloc(npar); /* Direction to move */ + gsl_matrix *myH = gsl_matrix_alloc(npar,npar); // Hessian matrix. + gsl_vector *stBeta = gsl_vector_alloc(npar); // Direction to move. - gsl_vector *myG = gsl_vector_alloc(npar); /* Gradient*/ - gsl_vector *tau = gsl_vector_alloc(npar); /* tau for QR*/ + gsl_vector *myG = gsl_vector_alloc(npar); // Gradient. + gsl_vector *tau = gsl_vector_alloc(npar); // tau for QR. for(iter=0;iter<100;iter++){ - wgsl_cat_optim_hessian(beta,&p,myH); //Calculate Hessian - wgsl_cat_optim_df(beta,&p,myG); //Calculate Gradient - gsl_linalg_QR_decomp(myH,tau); //Calculate next beta + wgsl_cat_optim_hessian(beta,&p,myH); // Calculate Hessian. + wgsl_cat_optim_df(beta,&p,myG); // Calculate Gradient. + gsl_linalg_QR_decomp(myH,tau); // Calculate next beta. gsl_linalg_QR_solve(myH,tau,myG,stBeta); gsl_vector_sub(beta,stBeta); - //Monitor convergence + // Monitor convergence. maxchange=0; for(int i=0;idata[i])) @@ -520,21 +486,14 @@ int logistic_cat_fit(gsl_vector *beta #ifdef _RPR_DEBUG_ mLogLik = wgsl_cat_optim_f(beta,&p); - //fprintf(stderr,"#iter %d, -log(Lik(0))=%lf,%lf\n",(int)iter,mLogLik,maxchange); -#endif //_RPR_DEBUG +#endif if(maxchange<1E-4) break; } -#ifdef _RPR_DEBUG_ - //for (int i = 0; i < npar; i++) - // fprintf(stderr,"#par_%d= %lf\n",i,beta->data[i]); -#endif //_RPR_DEBUG - - //Final fit + // Final fit. mLogLik = wgsl_cat_optim_f(beta,&p); - //fprintf(stderr,"#Final %d) -log(Lik(0))=%lf, maxchange %lf\n",iter,mLogLik,maxchange); gsl_vector_free (tau); gsl_vector_free (stBeta); @@ -544,41 +503,40 @@ int logistic_cat_fit(gsl_vector *beta return 0; } - /***************/ /* Continuous */ /***************/ -// I need to bundle all the data that goes to the function to optimze together. +// I need to bundle all the data that goes to the function to optimze +// together. typedef struct{ - gsl_matrix *Xc; // continuous covariates Matrix Nobs x Kc + gsl_matrix *Xc; // continuous covariates; Matrix Nobs x Kc gsl_vector *y; double lambdaL1; double lambdaL2; }fix_parm_cont_T; - -double fLogit_cont(gsl_vector *beta - ,gsl_matrix *Xc - ,gsl_vector *y - ,double lambdaL1 - ,double lambdaL2) -{ +double fLogit_cont(gsl_vector *beta, gsl_matrix *Xc, gsl_vector *y, + double lambdaL1, double lambdaL2) { int n = y->size; int npar = beta->size; double total = 0; double aux = 0; - /* omp_set_num_threads(ompthr); */ - /* /\* Changed loop start at 1 instead of 0 to avoid regularization of beta_0*\/ */ - /* /\*#pragma omp parallel for reduction (+:total)*\/ */ + + // omp_set_num_threads(ompthr); /\* Changed loop start at 1 instead + // of 0 to avoid regularization of beta_0*\/ /\*#pragma omp parallel + // for reduction (+:total)*\/ for(int i = 1; i < npar; ++i) total += beta->data[i]*beta->data[i]; total = (-total*lambdaL2/2); - /* /\*#pragma omp parallel for reduction (+:aux)*\/ */ + + // /\*#pragma omp parallel for reduction (+:aux)*\/ for(int i = 1; i < npar; ++i) aux += (beta->data[i]>0 ? beta->data[i] : -beta->data[i]); total = total-aux*lambdaL1; - /* #pragma omp parallel for schedule(static) shared(n,beta,X,nlev,y) reduction (+:total) */ + + // #pragma omp parallel for schedule(static) shared(n,beta,X,nlev,y) + // #reduction (+:total) for(int i = 0; i < n; ++i) { double Xbetai=beta->data[0]; int iParm=1; @@ -589,12 +547,12 @@ double fLogit_cont(gsl_vector *beta return -total; } - -void logistic_cont_pred(gsl_vector *beta // Vector of parameters length = 1 + Sum_k(C_k - 1) - ,gsl_matrix *Xc // continuous covariates Matrix Nobs x Kc (NULL if not used) - ,gsl_vector *yhat //Vector of prob. predicted by the logistic - ) -{ +void logistic_cont_pred(gsl_vector *beta, // Vector of parameters + // length = 1 + Sum_k(C_k-1). + gsl_matrix *Xc, // Continuous covariates matrix, + // Nobs x Kc (NULL if not used). + ,gsl_vector *yhat) { // Vector of prob. predicted by + // the logistic. for(int i = 0; i < Xc->size1; ++i) { double Xbetai=beta->data[0]; int iParm=1; @@ -604,20 +562,19 @@ void logistic_cont_pred(gsl_vector *beta // Vector of parameters length = 1 + S } } - -/* The gradient of f, df = (df/dx, df/dy). */ -void -wgsl_cont_optim_df (const gsl_vector *beta, void *params, - gsl_vector *out) -{ +// The gradient of f, df = (df/dx, df/dy). +void wgsl_cont_optim_df (const gsl_vector *beta, void *params, + gsl_vector *out) { fix_parm_cont_T *p = (fix_parm_cont_T *)params; int n = p->y->size; int Kc = p->Xc->size2; - int npar = beta->size; + int npar = beta->size; + // Intitialize gradient out necessary? for(int i = 0; i < npar; ++i) - out->data[i]= 0; - /* Changed loop start at 1 instead of 0 to avoid regularization of beta 0 */ + out->data[i]= 0; + + // Changed loop start at 1 instead of 0 to avoid regularization of beta 0. for(int i = 1; i < npar; ++i) out->data[i]= p->lambdaL2*beta->data[i]; for(int i = 1; i < npar; ++i) @@ -634,31 +591,34 @@ wgsl_cont_optim_df (const gsl_vector *beta, void *params, out->data[0]+= pn; iParm=1; - // Adding the continuous + + // Adding the continuous. for(int k = 0; k < Kc; ++k) { out->data[iParm++] += gsl_matrix_get(p->Xc,i,k)*pn; } } } - -/* The Hessian of f */ -void -wgsl_cont_optim_hessian (const gsl_vector *beta, void *params, - gsl_matrix *out) -{ +// The Hessian of f. +void wgsl_cont_optim_hessian (const gsl_vector *beta, void *params, + gsl_matrix *out) { fix_parm_cont_T *p = (fix_parm_cont_T *)params; int n = p->y->size; int Kc = p->Xc->size2; int npar = beta->size; - gsl_vector *gn = gsl_vector_alloc(npar); /* gn */ + gsl_vector *gn = gsl_vector_alloc(npar); // gn. + // Intitialize Hessian out necessary ??? + gsl_matrix_set_zero(out); - /* Changed loop start at 1 instead of 0 to avoid regularization of beta 0*/ + + // Changed loop start at 1 instead of 0 to avoid regularization of + // beta 0. for(int i = 1; i < npar; ++i) - gsl_matrix_set(out,i,i,(p->lambdaL2)); // Double check this - // L1 penalty not working yet, as not differentiable, I may need to do coordinate descent (as in glm_net) + gsl_matrix_set(out,i,i,(p->lambdaL2)); // Double-check this. + // L1 penalty not working yet, as not differentiable, I may need to + // do coordinate descent (as in glm_net). for(int i = 0; i < n; ++i) { double pn=0; double aux=0; @@ -668,10 +628,11 @@ wgsl_cont_optim_hessian (const gsl_vector *beta, void *params, Xbetai+= gsl_matrix_get(p->Xc,i,k)*beta->data[iParm1++]; pn= 1/(1 + gsl_sf_exp(-Xbetai)); + // Add a protection for pn very close to 0 or 1? aux=pn*(1-pn); - // Calculate sub-gradient vector gn + // Calculate sub-gradient vector gn. gsl_vector_set_zero(gn); gn->data[0]= 1; iParm1=1; @@ -688,32 +649,26 @@ wgsl_cont_optim_hessian (const gsl_vector *beta, void *params, gsl_vector_free(gn); } - -double wgsl_cont_optim_f(gsl_vector *v, void *params) -{ +double wgsl_cont_optim_f(gsl_vector *v, void *params) { double mLogLik=0; fix_parm_cont_T *p = (fix_parm_cont_T *)params; mLogLik = fLogit_cont(v,p->Xc,p->y,p->lambdaL1,p->lambdaL2); return mLogLik; } - -/* Compute both f and df together. */ -void -wgsl_cont_optim_fdf (gsl_vector *x, void *params, - double *f, gsl_vector *df) -{ +// Compute both f and df together. +void wgsl_cont_optim_fdf (gsl_vector *x, void *params, + double *f, gsl_vector *df) { *f = wgsl_cont_optim_f(x, params); wgsl_cont_optim_df(x, params, df); } - -int logistic_cont_fit(gsl_vector *beta - ,gsl_matrix *Xc // continuous covariates Matrix Nobs x Kc (NULL if not used) - ,gsl_vector *y - ,double lambdaL1 - ,double lambdaL2) -{ +int logistic_cont_fit (gsl_vector *beta, + gsl_matrix *Xc, // Continuous covariates matrix, + // Nobs x Kc (NULL if not used). + gsl_vector *y, + double lambdaL1, + double lambdaL2) { double mLogLik=0; fix_parm_cont_T p; @@ -721,32 +676,29 @@ int logistic_cont_fit(gsl_vector *beta int iter=0; double maxchange=0; - //Intializing fix parameters + // Initializing fix parameters. p.Xc=Xc; p.y=y; p.lambdaL1=lambdaL1; p.lambdaL2=lambdaL2; - //Initial fit - //#ifdef _RPR_DEBUG_ + // Initial fit. mLogLik = wgsl_cont_optim_f(beta,&p); - //fprintf(stderr,"#Initial -log(Lik(0))=%lf\n",mLogLik); - //#endif //_RPR_DEBUG - gsl_matrix *myH = gsl_matrix_alloc(npar,npar); /* Hessian matrix*/ - gsl_vector *stBeta = gsl_vector_alloc(npar); /* Direction to move */ + gsl_matrix *myH = gsl_matrix_alloc(npar,npar); // Hessian matrix. + gsl_vector *stBeta = gsl_vector_alloc(npar); // Direction to move. - gsl_vector *myG = gsl_vector_alloc(npar); /* Gradient*/ - gsl_vector *tau = gsl_vector_alloc(npar); /* tau for QR*/ + gsl_vector *myG = gsl_vector_alloc(npar); // Gradient. + gsl_vector *tau = gsl_vector_alloc(npar); // tau for QR. for(iter=0;iter<100;iter++){ - wgsl_cont_optim_hessian(beta,&p,myH); //Calculate Hessian - wgsl_cont_optim_df(beta,&p,myG); //Calculate Gradient - gsl_linalg_QR_decomp(myH,tau); //Calculate next beta + wgsl_cont_optim_hessian(beta,&p,myH); // Calculate Hessian. + wgsl_cont_optim_df(beta,&p,myG); // Calculate Gradient. + gsl_linalg_QR_decomp(myH,tau); // Calculate next beta. gsl_linalg_QR_solve(myH,tau,myG,stBeta); gsl_vector_sub(beta,stBeta); - //Monitor convergence + // Monitor convergence. maxchange=0; for(int i=0;idata[i])) @@ -754,22 +706,15 @@ int logistic_cont_fit(gsl_vector *beta #ifdef _RPR_DEBUG_ mLogLik = wgsl_cont_optim_f(beta,&p); - //fprintf(stderr,"#iter %d, -log(Lik(0))=%lf,%lf\n",(int)iter,mLogLik,maxchange); -#endif //_RPR_DEBUG +#endif if(maxchange<1E-4) break; } -#ifdef _RPR_DEBUG_ - //for (int i = 0; i < npar; i++) - // fprintf(stderr,"#par_%d= %lf\n",i,beta->data[i]); -#endif //_RPR_DEBUG - - //Final fit + // Final fit. mLogLik = wgsl_cont_optim_f(beta,&p); - //fprintf(stderr,"#Final %d) -log(Lik(0))=%lf, maxchange %lf\n",iter,mLogLik,maxchange); - + gsl_vector_free (tau); gsl_vector_free (stBeta); gsl_vector_free (myG); diff --git a/src/vc.cpp b/src/vc.cpp index 1e0c562..9fe1894 100644 --- a/src/vc.cpp +++ b/src/vc.cpp @@ -1,6 +1,6 @@ /* Genome-wide Efficient Mixed Model Association (GEMMA) - Copyright (C) 2011 Xiang Zhou + Copyright (C) 2011-2017, Xiang Zhou This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by @@ -13,10 +13,8 @@ GNU General Public License for more details. You should have received a copy of the GNU General Public License - along with this program. If not, see . - */ - - + along with this program. If not, see . +*/ #include #include @@ -51,25 +49,14 @@ #include "eigenlib.h" #include "gzstream.h" #include "mathfunc.h" - -#ifdef FORCE_FLOAT -#include "lmm_float.h" -#include "vc_float.h" -#else #include "lmm.h" #include "vc.h" -#endif - - using namespace std; using namespace Eigen; -//in this file, X, Y are already transformed (i.e. UtX and UtY) - - -void VC::CopyFromParam (PARAM &cPar) -{ +// In this file, X, Y are already transformed (i.e. UtX and UtY). +void VC::CopyFromParam (PARAM &cPar) { a_mode=cPar.a_mode; file_cat=cPar.file_cat; @@ -81,8 +68,6 @@ void VC::CopyFromParam (PARAM &cPar) file_out=cPar.file_out; path_out=cPar.path_out; - //v_sigma2=cPar.v_sigma2; - time_UtX=0.0; time_opt=0.0; @@ -102,9 +87,7 @@ void VC::CopyFromParam (PARAM &cPar) return; } - -void VC::CopyToParam (PARAM &cPar) -{ +void VC::CopyToParam (PARAM &cPar) { cPar.time_UtX=time_UtX; cPar.time_opt=time_opt; @@ -128,16 +111,18 @@ void VC::CopyToParam (PARAM &cPar) return; } - - -void VC::WriteFile_qs (const gsl_vector *s_vec, const gsl_vector *q_vec, const gsl_vector *qvar_vec, const gsl_matrix *S_mat, const gsl_matrix *Svar_mat) -{ +void VC::WriteFile_qs (const gsl_vector *s_vec, const gsl_vector *q_vec, + const gsl_vector *qvar_vec, const gsl_matrix *S_mat, + const gsl_matrix *Svar_mat) { string file_str; file_str=path_out+"/"+file_out; file_str+=".qvec.txt"; ofstream outfile_q (file_str.c_str(), ofstream::out); - if (!outfile_q) {cout<<"error writing file: "<size; i++) { outfile_q<size1; i++) { for (size_t j=0; jsize2; j++) { @@ -177,15 +165,7 @@ void VC::WriteFile_qs (const gsl_vector *s_vec, const gsl_vector *q_vec, const g return; } - - - - - - - -void UpdateParam (const gsl_vector *log_sigma2, VC_PARAM *p) -{ +void UpdateParam (const gsl_vector *log_sigma2, VC_PARAM *p) { size_t n1=(p->K)->size1, n_vc=log_sigma2->size-1, n_cvt=(p->W)->size2; gsl_matrix *K_temp=gsl_matrix_alloc(n1, n1); @@ -195,17 +175,19 @@ void UpdateParam (const gsl_vector *log_sigma2, VC_PARAM *p) gsl_matrix *WtHiWiWtHi=gsl_matrix_alloc(n_cvt, n1); double sigma2; - //calculate H=\sum_i^{k+1} \sigma_i^2 K_i + + // Calculate H = \sum_i^{k+1} \sigma_i^2 K_i. gsl_matrix_set_zero (p->P); for (size_t i=0; iK, 0, n1*i, n1, n1); + gsl_matrix_const_view K_sub= + gsl_matrix_const_submatrix (p->K, 0, n1*i, n1, n1); gsl_matrix_memcpy (K_temp, &K_sub.matrix); } - //when unconstrained, update on sigma2 instead of log_sigma2 + // When unconstrained, update on sigma2 instead of log_sigma2. if (p->noconstrain) { sigma2=gsl_vector_get (log_sigma2, i); } else { @@ -214,42 +196,22 @@ void UpdateParam (const gsl_vector *log_sigma2, VC_PARAM *p) gsl_matrix_scale(K_temp, sigma2); gsl_matrix_add (p->P, K_temp); } - - //calculate H^{-1} - /* - int sig; - gsl_permutation * pmt1=gsl_permutation_alloc (n1); - LUDecomp (p->P, pmt1, &sig); - LUInvert (p->P, pmt1, K_temp); - gsl_permutation_free(pmt1); - - gsl_matrix_memcpy (p->P, K_temp); - */ + + // Calculate H^{-1}. eigenlib_invert(p->P); - //calculate P=H^{-1}-H^{-1}W(W^TH^{-1}W)^{-1}W^TH^{-1} - //gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, p->P, p->W, 0.0, HiW); - //gsl_blas_dgemm (CblasTrans, CblasNoTrans, 1.0, p->W, HiW, 0.0, WtHiW); - eigenlib_dgemm ("N", "N", 1.0, p->P, p->W, 0.0, HiW); eigenlib_dgemm ("T", "N", 1.0, p->W, HiW, 0.0, WtHiW); - //gsl_permutation * pmt2=gsl_permutation_alloc (n_cvt); - //LUDecomp (WtHiW, pmt2, &sig); - //LUInvert (WtHiW, pmt2, WtHiWi); - //gsl_permutation_free(pmt2); eigenlib_invert(WtHiW); gsl_matrix_memcpy(WtHiWi, WtHiW); - //gsl_blas_dgemm (CblasNoTrans, CblasTrans, 1.0, WtHiWi, HiW, 0.0, WtHiWiWtHi); - //gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, -1.0, HiW, WtHiWiWtHi, 1.0, p->P); eigenlib_dgemm ("N", "T", 1.0, WtHiWi, HiW, 0.0, WtHiWiWtHi); eigenlib_dgemm ("N", "N", -1.0, HiW, WtHiWiWtHi, 1.0, p->P); - //calculate Py, KPy, PKPy + // Calculate Py, KPy, PKPy. gsl_blas_dgemv(CblasNoTrans, 1.0, p->P, p->y, 0.0, p->Py); - //eigenlib_dgemv("N", 1.0, p->P, p->y, 0.0, p->Py); - + double d; for (size_t i=0; iKPy_mat, i); @@ -259,20 +221,29 @@ void UpdateParam (const gsl_vector *log_sigma2, VC_PARAM *p) gsl_vector_memcpy (&KPy.vector, p->Py); } else { gsl_matrix_const_view K_sub=gsl_matrix_const_submatrix (p->K, 0, n1*i, n1, n1); - //seems to be important to use gsl dgemv here instead of eigenlib_dgemv; otherwise - gsl_blas_dgemv(CblasNoTrans, 1.0, &K_sub.matrix, p->Py, 0.0, &KPy.vector); - //eigenlib_dgemv("N", 1.0, &K_sub.matrix, p->Py, 0.0, &KPy.vector); + + // Seems to be important to use gsl dgemv here instead of + // eigenlib_dgemv; otherwise. + gsl_blas_dgemv(CblasNoTrans, 1.0, &K_sub.matrix, p->Py, 0.0, + &KPy.vector); } gsl_blas_dgemv(CblasNoTrans, 1.0, p->P, &KPy.vector, 0.0, &PKPy.vector); - //eigenlib_dgemv("N", 1.0, p->P, &KPy.vector, 0.0, &PKPy.vector); - //when phenotypes are not normalized well, then some values in the following matrix maybe nan; change that to 0; this seems to only happen when eigenlib_dgemv was used above + // When phenotypes are not normalized well, then some values in + // the following matrix maybe NaN; change that to 0; this seems to + // only happen when eigenlib_dgemv was used above. for (size_t j=0; jKPy_mat->size1; j++) { d=gsl_matrix_get (p->KPy_mat, j, i); - if (std::isnan(d)) {gsl_matrix_set (p->KPy_mat, j, i, 0); cout<<"nan appears in "<KPy_mat, j, i, 0); + cout<<"nan appears in "<PKPy_mat, j, i); - if (std::isnan(d)) {gsl_matrix_set (p->PKPy_mat, j, i, 0); cout<<"nan appears in "<PKPy_mat, j, i, 0); + cout<<"nan appears in "<K)->size1, n_vc=log_sigma2->size-1; double tr, d; - //update parameters + // Update parameters. UpdateParam (log_sigma2, p); - //calculate dev1=-0.5*trace(PK_i)+0.5*yPKPy + // Calculate dev1=-0.5*trace(PK_i)+0.5*yPKPy. for (size_t i=0; isize-1; double d, sigma2_i, sigma2_j; - //update parameters + // Update parameters. UpdateParam (log_sigma2, p); - //calculate dev2=0.5(yPKPKPy) + // Calculate dev2 = 0.5(yPKPKPy). for (size_t i=0; iKPy_mat, i); if (p->noconstrain) { @@ -370,32 +336,20 @@ int LogRL_dev2 (const gsl_vector *log_sigma2, void *params, gsl_matrix *dev2) } gsl_matrix_memcpy (p->Hessian, dev2); - /* - for (size_t i=0; isize1; i++) { - for (size_t j=0; jsize2; j++) { - cout<K)->size1, n_vc=log_sigma2->size-1; double tr, d, sigma2_i, sigma2_j; - //update parameters + // Update parameters. UpdateParam (log_sigma2, p); - //calculate dev1=(-0.5*trace(PK_i)+0.5*yPK_iPy)*sigma2_i - //calculate dev2=0.5(yPK_iPK_jPy)*sigma2_i*sigma2_j for (size_t i=0; i rs_set(rs_ptr, rs_ptr+10); string chr_ptr[]={"chr","CHR"}; set chr_set(chr_ptr, chr_ptr+2); - string pos_ptr[]={"ps","PS","pos","POS","base_position","BASE_POSITION", "bp", "BP"}; + string pos_ptr[]={"ps","PS","pos","POS","base_position","BASE_POSITION", + "bp", "BP"}; set pos_set(pos_ptr, pos_ptr+8); string cm_ptr[]={"cm","CM"}; set cm_set(cm_ptr, cm_ptr+2); @@ -486,7 +437,8 @@ bool ReadHeader_vc (const string &line, HEADER &header) string nobs_ptr[]={"nobs","NOBS","n_obs","N_OBS"}; set nobs_set(nobs_ptr, nobs_ptr+4); - string af_ptr[]={"af","AF","maf","MAF","f","F","allele_freq","ALLELE_FREQ","allele_frequency","ALLELE_FREQUENCY"}; + string af_ptr[]={"af","AF","maf","MAF","f","F","allele_freq", + "ALLELE_FREQ","allele_frequency","ALLELE_FREQUENCY"}; set af_set(af_ptr, af_ptr+10); string var_ptr[]={"var","VAR"}; set var_set(var_ptr, var_ptr+2); @@ -496,7 +448,11 @@ bool ReadHeader_vc (const string &line, HEADER &header) string cor_ptr[]={"cor","COR","r","R"}; set cor_set(cor_ptr, cor_ptr+4); - header.rs_col=0; header.chr_col=0; header.pos_col=0; header.a1_col=0; header.a0_col=0; header.z_col=0; header.beta_col=0; header.sebeta_col=0; header.chisq_col=0; header.p_col=0; header.n_col=0; header.nmis_col=0; header.nobs_col=0; header.af_col=0; header.var_col=0; header.ws_col=0; header.cor_col=0; header.coln=0; + header.rs_col=0; header.chr_col=0; header.pos_col=0; header.a1_col=0; + header.a0_col=0; header.z_col=0; header.beta_col=0; header.sebeta_col=0; + header.chisq_col=0; header.p_col=0; header.n_col=0; header.nmis_col=0; + header.nobs_col=0; header.af_col=0; header.var_col=0; header.ws_col=0; + header.cor_col=0; header.coln=0; char *ch_ptr; string type; @@ -506,46 +462,134 @@ bool ReadHeader_vc (const string &line, HEADER &header) while (ch_ptr!=NULL) { type=ch_ptr; if (rs_set.count(type)!=0) { - if (header.rs_col==0) {header.rs_col=header.coln+1;} else {cout<<"error! more than two rs columns in the file."< &setSnps, vector &vec_rs, vector &vec_n, vector &vec_cm, vector &vec_bp, map &mapRS2in, map &mapRS2var) -{ +// Read cov file the first time, record mapRS2in, mapRS2var (in case +// var is not provided in the z file), store vec_n and vec_rs. +void ReadFile_cor (const string &file_cor, const set &setSnps, + vector &vec_rs, vector &vec_n, + vector &vec_cm, vector &vec_bp, + map &mapRS2in, map &mapRS2var) { vec_rs.clear(); vec_n.clear(); mapRS2in.clear(); mapRS2var.clear(); igzstream infile (file_cor.c_str(), igzstream::in); - if (!infile) {cout<<"error! fail to open cov file: "< &setSnps, vector &setSnps, vector &setSnps, vector0) {vec_cm.push_back(d_cm);} else {vec_cm.push_back(d_cm);} if (d_pos>0) {vec_bp.push_back(d_pos);} else {vec_bp.push_back(d_pos);} - //record mapRS2in and mapRS2var + // Record mapRS2in and mapRS2var. if (setSnps.size()==0 || setSnps.count(rs)!=0) { if (mapRS2in.count(rs)==0) { mapRS2in[rs]=1; @@ -648,37 +699,39 @@ void ReadFile_cor (const string &file_cor, const set &setSnps, vector1) { - cout<<"error! more than one snp has the same id "< &vec_rs, const vector &vec_n, const vector &vec_cm, const vector &vec_bp, const map &mapRS2cat, const map &mapRS2in, const map &mapRS2var, const map &mapRS2nsamp, const size_t crt, const double &window_cm, const double &window_bp, const double &window_ns, gsl_matrix *S_mat, gsl_matrix *Svar_mat, gsl_vector *qvar_vec, size_t &ni_total, size_t &ns_total, size_t &ns_test, size_t &ns_pair) -{ +// Read covariance file the second time. +// Look for rs, n_mis+n_obs, var, window_size, cov. +// If window_cm/bp/ns is provided, then use these max values to +// calibrate estimates. +void ReadFile_cor (const string &file_cor, const vector &vec_rs, + const vector &vec_n, const vector &vec_cm, + const vector &vec_bp, + const map &mapRS2cat, + const map &mapRS2in, + const map &mapRS2var, + const map &mapRS2nsamp, + const size_t crt, const double &window_cm, + const double &window_bp, const double &window_ns, + gsl_matrix *S_mat, gsl_matrix *Svar_mat, + gsl_vector *qvar_vec, size_t &ni_total, + size_t &ns_total, size_t &ns_test, size_t &ns_pair) { igzstream infile (file_cor.c_str(), igzstream::in); - if (!infile) {cout<<"error! fail to open cov file: "< &vec_rs, const v double d_pos1, d_pos2, d_pos, d_cm1, d_cm2, d_cm; ns_test=0; ns_total=0; ns_pair=0; ni_total=0; - //header + // Header. HEADER header; !safeGetline(infile, line).eof(); ReadHeader_vc (line, header); while (!safeGetline(infile, line).eof()) { - //do not read cor values this time; upto col_n-1 + + // Do not read cor values this time; upto col_n-1. d_pos1=0; d_cm1=0; ch_ptr=strtok ((char *)line.c_str(), " , \t"); for (size_t i=0; i &vec_rs, const v rs1=rs; - if ( (mapRS2cat.size()==0 || mapRS2cat.count(rs1)!=0) && mapRS2in.count(rs1)!=0 && mapRS2in.at(rs1)==2) { + if ( (mapRS2cat.size()==0 || mapRS2cat.count(rs1)!=0) && + mapRS2in.count(rs1)!=0 && mapRS2in.at(rs1)==2) { var1=mapRS2var.at(rs1); nsamp1=mapRS2nsamp.at(rs1); d2=var1*var1; if (mapRS2cat.size()!=0) { - mat_S[mapRS2cat.at(rs1) ][mapRS2cat.at(rs1) ]+=(1-1.0/(double)vec_n[ns_total])*d2; - mat_Svar[mapRS2cat.at(rs1) ][mapRS2cat.at(rs1) ]+=d2*d2/((double)vec_n[ns_total]*(double)vec_n[ns_total]); + mat_S[mapRS2cat.at(rs1) ][mapRS2cat.at(rs1) ]+= + (1-1.0/(double)vec_n[ns_total])*d2; + mat_Svar[mapRS2cat.at(rs1) ][mapRS2cat.at(rs1) ]+= + d2*d2/((double)vec_n[ns_total]*(double)vec_n[ns_total]); if (crt==1) { - mat3d_Sbin[mapRS2cat.at(rs1) ][mapRS2cat.at(rs1) ][0]+=(1-1.0/(double)vec_n[ns_total])*d2; + mat3d_Sbin[mapRS2cat.at(rs1) ][mapRS2cat.at(rs1) ][0]+= + (1-1.0/(double)vec_n[ns_total])*d2; } } else { - //mat_S[0][0]+=(1-1.0/(double)vec_n[ns_total])*d2; mat_S[0][0]+=(1-1.0/(double)vec_n[ns_total])*d2; - mat_Svar[0][0]+=d2*d2/((double)vec_n[ns_total]*(double)vec_n[ns_total]); + mat_Svar[0][0]+= + d2*d2/((double)vec_n[ns_total]*(double)vec_n[ns_total]); if (crt==1) { mat3d_Sbin[0][0][0]+=(1-1.0/(double)vec_n[ns_total])*d2; } @@ -923,7 +1004,8 @@ void ReadFile_cor (const string &file_cor, const vector &vec_rs, const v n_nb=0; while(ch_ptr!=NULL) { type=ch_ptr; - if (type.compare("NA")!=0 && type.compare("na")!=0 && type.compare("nan")!=0 && type.compare("-nan")!=0) { + if (type.compare("NA")!=0 && type.compare("na")!=0 && + type.compare("nan")!=0 && type.compare("-nan")!=0) { cor=atof(ch_ptr); rs2=vec_rs[ns_total+n_nb+1]; d_pos2=vec_bp[ns_total+n_nb+1]; @@ -931,40 +1013,45 @@ void ReadFile_cor (const string &file_cor, const vector &vec_rs, const v d_pos=abs(d_pos2-d_pos1); d_cm=abs(d_cm2-d_cm1); - if ( (mapRS2cat.size()==0 || mapRS2cat.count(rs2)!=0) && mapRS2in.count(rs2)!=0 && mapRS2in.at(rs2)==2) { + if ( (mapRS2cat.size()==0 || mapRS2cat.count(rs2)!=0) && + mapRS2in.count(rs2)!=0 && mapRS2in.at(rs2)==2) { var2=mapRS2var.at(rs2); nsamp2=mapRS2nsamp.at(rs2); - d1=cor*cor-1.0/(double)min(vec_n[ns_total], vec_n[ns_total+n_nb+1]); + d1=cor*cor-1.0/(double)min(vec_n[ns_total], + vec_n[ns_total+n_nb+1]); d2=var1*var2; d3=cor*cor/((double)nsamp1*(double)nsamp2); n12=min(vec_n[ns_total], vec_n[ns_total+n_nb+1]); - //compute bin + // Compute bin. if (crt==1) { if (window_cm!=0 && d_cm1!=0 && d_cm2!=0) { bin=min( (int)floor(d_cm/window_cm*bin_size), (int)bin_size); } else if (window_bp!=0 && d_pos1!=0 && d_pos2!=0) { bin=min( (int)floor(d_pos/window_bp*bin_size), (int)bin_size); } else if (window_ns!=0) { - bin=min( (int)floor(((double)n_nb+1)/window_ns*bin_size), (int)bin_size); + bin=min( (int)floor(((double)n_nb+1)/window_ns*bin_size), + (int)bin_size); } } - //if (mat_S[0][0]!=mat_S[0][0] && flag_nan==0) { - //if (rs1.compare("rs10915560")==0 || rs1.compare("rs241273")==0) {cout< &vec_rs, const v ns_total++; } - //use S_bin to fit a rational function y=1/(a+bx)^2, where x=seq(0.5,bin_size-0.5,by=1) - //and then compute a correlation factor as a percentage + // Use S_bin to fit a rational function y=1/(a+bx)^2, where + // x=seq(0.5,bin_size-0.5,by=1) and then compute a correlation + // factor as a percentage. double a, b, x, y, n, var_y, var_x, mean_y, mean_x, cov_xy, crt_factor; if (crt==1) { for (size_t i=0; isize1; i++) { for (size_t j=i; jsize2; j++) { - //correct mat_S + // Correct mat_S. n=0; var_y=0; var_x=0; mean_y=0; mean_x=0; cov_xy=0; for (size_t k=0; k &vec_rs, const v mat_S[i][j]*=crt_factor; mat_S[j][i]*=crt_factor; } cout<size1; i++) { d1=gsl_vector_get(qvar_vec, i)+2*vec_qvar[i]; gsl_vector_set(qvar_vec, i, d1); for (size_t j=0; jsize2; j++) { if (i==j) { gsl_matrix_set(S_mat, i, j, mat_S[i][i]); - gsl_matrix_set(Svar_mat, i, j, 2.0*mat_Svar[i][i]*ns_test*ns_test/(2.0*ns_pair) ); + gsl_matrix_set(Svar_mat, i, j, + 2.0*mat_Svar[i][i]*ns_test*ns_test/(2.0*ns_pair) ); } else { gsl_matrix_set(S_mat, i, j, mat_S[i][j]+mat_S[j][i]); - gsl_matrix_set(Svar_mat, i, j, 2.0*(mat_Svar[i][j]+mat_Svar[j][i])*ns_test*ns_test/(2.0*ns_pair) ); + gsl_matrix_set(Svar_mat, i, j, + 2.0*(mat_Svar[i][j]+mat_Svar[j][i])* + ns_test*ns_test/(2.0*ns_pair) ); } } } @@ -1059,9 +1151,18 @@ void ReadFile_cor (const string &file_cor, const vector &vec_rs, const v return; } -//use the new method to calculate variance components with summary statistics -//first, use a function CalcS to compute S matrix (where the diagonal elements are part of V(q) ), and then use bootstrap to compute the variance for S, use a set of genotypes, phenotypes, and individual ids, and snp category label -void CalcVCss(const gsl_matrix *Vq, const gsl_matrix *S_mat, const gsl_matrix *Svar_mat, const gsl_vector *q_vec, const gsl_vector *s_vec, const double df, vector &v_pve, vector &v_se_pve, double &pve_total, double &se_pve_total, vector &v_sigma2, vector &v_se_sigma2, vector &v_enrich, vector &v_se_enrich) { +// Use the new method to calculate variance components with summary +// statistics first, use a function CalcS to compute S matrix (where +// the diagonal elements are part of V(q) ), and then use bootstrap to +// compute the variance for S, use a set of genotypes, phenotypes, and +// individual ids, and snp category label. +void CalcVCss(const gsl_matrix *Vq, const gsl_matrix *S_mat, + const gsl_matrix *Svar_mat, const gsl_vector *q_vec, + const gsl_vector *s_vec, const double df, + vector &v_pve, vector &v_se_pve, + double &pve_total, double &se_pve_total, + vector &v_sigma2, vector &v_se_sigma2, + vector &v_enrich, vector &v_se_enrich) { size_t n_vc=S_mat->size1; gsl_matrix *Si_mat=gsl_matrix_alloc (n_vc, n_vc); @@ -1082,69 +1183,38 @@ void CalcVCss(const gsl_matrix *Vq, const gsl_matrix *S_mat, const gsl_matrix *S double d; - //calculate S^{-1}q + // Calculate S^{-1}q. gsl_matrix_memcpy (tmp_mat, S_mat); int sig; gsl_permutation * pmt=gsl_permutation_alloc (n_vc); LUDecomp (tmp_mat, pmt, &sig); LUInvert (tmp_mat, pmt, Si_mat); - //calculate sigma2snp and pve + // Calculate sigma2snp and pve. gsl_blas_dgemv (CblasNoTrans, 1.0, Si_mat, q_vec, 0.0, pve); gsl_vector_memcpy(sigma2persnp, pve); gsl_vector_div(sigma2persnp, s_vec); - //get qvar_mat - /* - if (n_block==0 || n_block==1) { - double s=1.0; - for (size_t i=0; isize1, n2=K->size2; size_t n_vc=n2/n1; @@ -1287,14 +1355,14 @@ void VC::CalcVChe (const gsl_matrix *K, const gsl_matrix *W, const gsl_vector *y double d, tr, s, v; vector traceG_new; - //new matrices/vectors + // New matrices/vectors. gsl_matrix *K_scale=gsl_matrix_alloc (n1, n2); gsl_vector *y_scale=gsl_vector_alloc (n1); gsl_matrix *Kry=gsl_matrix_alloc (n1, n_vc); gsl_matrix *yKrKKry=gsl_matrix_alloc (n_vc, n_vc*(n_vc+1) ); gsl_vector *KKry=gsl_vector_alloc (n1); - //old matrices/vectors + // Old matrices/vectors. gsl_vector *pve=gsl_vector_alloc (n_vc); gsl_vector *se_pve=gsl_vector_alloc (n_vc); gsl_vector *q_vec=gsl_vector_alloc (n_vc); @@ -1304,10 +1372,12 @@ void VC::CalcVChe (const gsl_matrix *K, const gsl_matrix *W, const gsl_vector *y gsl_matrix *Si_mat=gsl_matrix_alloc (n_vc, n_vc); gsl_matrix *Var_mat=gsl_matrix_alloc (n_vc, n_vc); - //center and scale K by W + // Center and scale K by W. for (size_t i=0; isize1, n2=K->size2; size_t n_vc=n2/n1; gsl_vector *log_sigma2=gsl_vector_alloc (n_vc+1); double d, s; - /* - //compare eigenlib vs lapack - //dgemm - gsl_matrix *K2=gsl_matrix_alloc(K->size1, K->size1); - - clock_t time_start=clock(); - gsl_blas_dgemm(CblasNoTrans, CblasTrans, 1.0, K, K, 0.0, K2); - cout<<"standard time: "<<(clock()-time_start)/(double(CLOCKS_PER_SEC)*60.0)<size1); - time_start=clock(); - for (size_t i=0; i<1000; i++) { - gsl_blas_dgemv(CblasNoTrans, 1.0, K2, &W_col.vector, 0.0, v); - } - cout<<"standard time: "<<(clock()-time_start)/(double(CLOCKS_PER_SEC)*60.0)<size1, K->size1); - gsl_matrix *K3=gsl_matrix_alloc(K->size1, K->size1); - - gsl_matrix_memcpy(K2copy, K2); - time_start=clock(); - EigenDecomp(K2copy, K3, v, 0); - cout<<"standard time 0: "<<(clock()-time_start)/(double(CLOCKS_PER_SEC)*60.0)<size1); - LUDecomp (K2copy, pmt1, &sigcopy); - LUInvert (K2copy, pmt1, K3); - gsl_permutation_free(pmt1); - cout<<"standard time: "<<(clock()-time_start)/(double(CLOCKS_PER_SEC)*60.0)<f, 1e-3); } while (status==GSL_CONTINUE && iterx, ¶ms, dev1, dev2); - /* - for (size_t i=0; isize1; i++) { - for (size_t j=0; jsize2; j++) { - cout<size1, n2=K->size2; size_t n_vc=n2/n1; double d, y2_sum, tau_inv, se_tau_inv; - //new matrices/vectors + // New matrices/vectors. gsl_matrix *K_scale=gsl_matrix_alloc (n1, n2); gsl_vector *y_scale=gsl_vector_alloc (n1); gsl_vector *y2=gsl_vector_alloc (n1); @@ -1860,7 +1775,7 @@ void VC::CalcVCacl (const gsl_matrix *K, const gsl_matrix *W, const gsl_vector * gsl_matrix *K_tmp=gsl_matrix_alloc (n1, n1); gsl_matrix *V_mat=gsl_matrix_alloc (n1, n1); - //old matrices/vectors + // Old matrices/vectors. gsl_vector *pve=gsl_vector_alloc (n_vc); gsl_vector *se_pve=gsl_vector_alloc (n_vc); gsl_vector *q_vec=gsl_vector_alloc (n_vc); @@ -1874,23 +1789,24 @@ void VC::CalcVCacl (const gsl_matrix *K, const gsl_matrix *W, const gsl_vector * int sig; gsl_permutation * pmt=gsl_permutation_alloc (n_vc); - //center and scale K by W - //and standardize K further so that all diagonal elements are 1 + // Center and scale K by W, and standardize K further so that all + // diagonal elements are 1 for (size_t i=0; isize; t++) { d+=gsl_vector_get (n1_vec, t); @@ -1939,7 +1860,7 @@ void VC::CalcVCacl (const gsl_matrix *K, const gsl_matrix *W, const gsl_vector * gsl_matrix_set (S2, i, j, d); if (i!=j) {gsl_matrix_set (S2, j, i, d);} - //save information to compute J + // Save information to compute J. gsl_vector_view K2col1=gsl_matrix_column (K2, n_vc*i+j); gsl_vector_view K2col2=gsl_matrix_column (K2, n_vc*j+i); @@ -1948,57 +1869,61 @@ void VC::CalcVCacl (const gsl_matrix *K, const gsl_matrix *W, const gsl_vector * } } - //iterate to solve tau and h's + // Iterate to solve tau and h's. size_t it=0; double s=1; while (abs(s)>1e-3 && it<100) { - //update tau_inv + + // Update tau_inv. gsl_blas_ddot (q_vec, pve, &d); if (it>0) {s=y2_sum/(double)n1-d/((double)n1*((double)n1-1))-tau_inv;} tau_inv=y2_sum/(double)n1-d/((double)n1*((double)n1-1)); if (it>0) {s/=tau_inv;} - //update S + // Update S. gsl_matrix_memcpy (S_mat, S2); gsl_matrix_scale (S_mat, -1*tau_inv); gsl_matrix_add (S_mat, S1); - //update h=S^{-1}q + // Update h=S^{-1}q. int sig; gsl_permutation * pmt=gsl_permutation_alloc (n_vc); LUDecomp (S_mat, pmt, &sig); LUInvert (S_mat, pmt, Si_mat); gsl_blas_dgemv (CblasNoTrans, 1.0, Si_mat, q_vec, 0.0, pve); - //cout< &indicator_idv, vector &indicator_snp, const vector &vec_cat, const gsl_vector *w, const gsl_vector *z, size_t ns_test, gsl_matrix *XWz) -{ +// Read bimbam mean genotype file and compute XWz. +bool BimbamXwz (const string &file_geno, const int display_pace, + vector &indicator_idv, vector &indicator_snp, + const vector &vec_cat, const gsl_vector *w, + const gsl_vector *z, size_t ns_test, gsl_matrix *XWz) { igzstream infile (file_geno.c_str(), igzstream::in); - //ifstream infile (file_geno.c_str(), ifstream::in); - if (!infile) {cout<<"error reading genotype file:"< &in for (size_t t=0; t &in for (size_t i=0; i &in geno_var+=geno_mean*geno_mean*(double)n_miss; geno_var/=(double)ni_test; geno_var-=geno_mean*geno_mean; -// geno_var=geno_mean*(1-geno_mean*0.5); for (size_t i=0; i &in return true; } - - - - - -//read plink bed file and compute XWz -bool PlinkXwz (const string &file_bed, const int display_pace, vector &indicator_idv, vector &indicator_snp, const vector &vec_cat, const gsl_vector *w, const gsl_vector *z, size_t ns_test, gsl_matrix *XWz) -{ +// Read PLINK bed file and compute XWz. +bool PlinkXwz (const string &file_bed, const int display_pace, + vector &indicator_idv, vector &indicator_snp, + const vector &vec_cat, const gsl_vector *w, + const gsl_vector *z, size_t ns_test, gsl_matrix *XWz) { ifstream infile (file_bed.c_str(), ios::binary); - if (!infile) {cout<<"error reading bed file:"< b; @@ -2216,38 +2143,60 @@ bool PlinkXwz (const string &file_bed, const int display_pace, vector &indi gsl_vector_mul(wz, w); int n_bit; - //calculate n_bit and c, the number of bit for each snp + + // Calculate n_bit and c, the number of bit for each snp. if (ni_total%4==0) {n_bit=ni_total/4;} else {n_bit=ni_total/4+1; } - //print the first three majic numbers + // Print the first three magic numbers. for (int i=0; i<3; ++i) { infile.read(ch,1); b=ch[0]; } for (size_t t=0; t &indi } } - geno_mean/=(double)(ni_test-n_miss); geno_var+=geno_mean*geno_mean*(double)n_miss; geno_var/=(double)ni_test; geno_var-=geno_mean*geno_mean; -// geno_var=geno_mean*(1-geno_mean*0.5); for (size_t i=0; i &indi gsl_vector_add_constant (geno, -1.0*geno_mean); - gsl_vector_view XWz_col=gsl_matrix_column(XWz, vec_cat[ns_test]); + gsl_vector_view XWz_col= + gsl_matrix_column(XWz, vec_cat[ns_test]); d=gsl_vector_get (wz, ns_test); gsl_blas_daxpy (d/sqrt(geno_var), geno, &XWz_col.vector); @@ -2286,15 +2234,19 @@ bool PlinkXwz (const string &file_bed, const int display_pace, vector &indi return true; } - - -//read multiple genotype files and compute XWz -bool MFILEXwz (const size_t mfile_mode, const string &file_mfile, const int display_pace, vector &indicator_idv, vector > &mindicator_snp, const vector &vec_cat, const gsl_vector *w, const gsl_vector *z, gsl_matrix *XWz) -{ +// Read multiple genotype files and compute XWz. +bool MFILEXwz (const size_t mfile_mode, const string &file_mfile, + const int display_pace, vector &indicator_idv, + vector > &mindicator_snp, + const vector &vec_cat, const gsl_vector *w, + const gsl_vector *z, gsl_matrix *XWz) { gsl_matrix_set_zero(XWz); igzstream infile (file_mfile.c_str(), igzstream::in); - if (!infile) {cout<<"error! fail to open mfile file: "< &indicator_idv, vector &indicator_snp, const gsl_matrix *XWz, size_t ns_test, gsl_matrix *XtXWz) -{ +// Read bimbam mean genotype file and compute X_i^TX_jWz. +bool BimbamXtXwz (const string &file_geno, const int display_pace, + vector &indicator_idv, vector &indicator_snp, + const gsl_matrix *XWz, size_t ns_test, gsl_matrix *XtXWz) { igzstream infile (file_geno.c_str(), igzstream::in); - //ifstream infile (file_geno.c_str(), ifstream::in); - if (!infile) {cout<<"error reading genotype file:"< & for (size_t t=0; t & for (size_t i=0; i & geno_var+=geno_mean*geno_mean*(double)n_miss; geno_var/=(double)ni_test; geno_var-=geno_mean*geno_mean; -// geno_var=geno_mean*(1-geno_mean*0.5); for (size_t i=0; isize2; i++) { - gsl_vector_const_view XWz_col=gsl_matrix_const_column(XWz, i); + gsl_vector_const_view XWz_col= + gsl_matrix_const_column(XWz, i); gsl_blas_ddot (geno, &XWz_col.vector, &d); gsl_matrix_set (XtXWz, ns_test, i, d/sqrt(geno_var)); } @@ -2398,16 +2356,15 @@ bool BimbamXtXwz (const string &file_geno, const int display_pace, vector & return true; } - - - - - -//read plink bed file and compute XWz -bool PlinkXtXwz (const string &file_bed, const int display_pace, vector &indicator_idv, vector &indicator_snp, const gsl_matrix *XWz, size_t ns_test, gsl_matrix *XtXWz) -{ +// Read PLINK bed file and compute XWz. +bool PlinkXtXwz (const string &file_bed, const int display_pace, + vector &indicator_idv, vector &indicator_snp, + const gsl_matrix *XWz, size_t ns_test, gsl_matrix *XtXWz) { ifstream infile (file_bed.c_str(), ios::binary); - if (!infile) {cout<<"error reading bed file:"< b; @@ -2421,11 +2378,11 @@ bool PlinkXtXwz (const string &file_bed, const int display_pace, vector &in int n_bit; - //calculate n_bit and c, the number of bit for each snp + // Calculate n_bit and c, the number of bit for each snp. if (ni_total%4==0) {n_bit=ni_total/4;} else {n_bit=ni_total/4+1; } - //print the first three majic numbers + // Print the first three magic numbers. for (int i=0; i<3; ++i) { infile.read(ch,1); b=ch[0]; @@ -2435,24 +2392,45 @@ bool PlinkXtXwz (const string &file_bed, const int display_pace, vector &in if (t%display_pace==0 || t==(indicator_snp.size()-1)) {ProgressBar ("Reading SNPs ", t, indicator_snp.size()-1);} if (indicator_snp[t]==0) {continue;} - infile.seekg(t*n_bit+3); //n_bit, and 3 is the number of magic numbers + // n_bit, and 3 is the number of magic numbers. + infile.seekg(t*n_bit+3); - //read genotypes + // Read genotypes. geno_mean=0.0; n_miss=0; ci_total=0; geno_var=0.0; ci_test=0; for (int i=0; i &in } } - geno_mean/=(double)(ni_test-n_miss); geno_var+=geno_mean*geno_mean*(double)n_miss; geno_var/=(double)ni_test; geno_var-=geno_mean*geno_mean; -// geno_var=geno_mean*(1-geno_mean*0.5); for (size_t i=0; i &in gsl_vector_add_constant (geno, -1.0*geno_mean); for (size_t i=0; isize2; i++) { - gsl_vector_const_view XWz_col=gsl_matrix_const_column(XWz, i); + gsl_vector_const_view XWz_col= + gsl_matrix_const_column(XWz, i); gsl_blas_ddot (geno, &XWz_col.vector, &d); gsl_matrix_set (XtXWz, ns_test, i, d/sqrt(geno_var)); } ns_test++; - } + } cout< &in return true; } - - -//read multiple genotype files and compute XWz -bool MFILEXtXwz (const size_t mfile_mode, const string &file_mfile, const int display_pace, vector &indicator_idv, vector > &mindicator_snp, const gsl_matrix *XWz, gsl_matrix *XtXWz) -{ +// Read multiple genotype files and compute XWz. +bool MFILEXtXwz (const size_t mfile_mode, const string &file_mfile, + const int display_pace, vector &indicator_idv, + vector > &mindicator_snp, const gsl_matrix *XWz, + gsl_matrix *XtXWz) { gsl_matrix_set_zero(XtXWz); igzstream infile (file_mfile.c_str(), igzstream::in); - if (!infile) {cout<<"error! fail to open mfile file: "< &vec_cat, const vector &v_pve, vector &v_se_pve, double &pve_total, double &se_pve_total, vector &v_sigma2, vector &v_se_sigma2, vector &v_enrich, vector &v_se_enrich) { +// Compute confidence intervals from summary statistics. +void CalcCIss(const gsl_matrix *Xz, const gsl_matrix *XWz, + const gsl_matrix *XtXWz, const gsl_matrix *S_mat, + const gsl_matrix *Svar_mat, const gsl_vector *w, + const gsl_vector *z, const gsl_vector *s_vec, + const vector &vec_cat, const vector &v_pve, + vector &v_se_pve, double &pve_total, + double &se_pve_total, vector &v_sigma2, + vector &v_se_sigma2, vector &v_enrich, + vector &v_se_enrich) { size_t n_vc=XWz->size2, ns_test=w->size, ni_test=XWz->size1; - //set up matrices + // Set up matrices. gsl_vector *w_pve=gsl_vector_alloc (ns_test); gsl_vector *wz=gsl_vector_alloc (ns_test); gsl_vector *zwz=gsl_vector_alloc (n_vc); @@ -2544,7 +2533,7 @@ void CalcCIss(const gsl_matrix *Xz, const gsl_matrix *XWz, const gsl_matrix *XtX double d, s0, s1, s, s_pve, s_snp; - //compute wz and zwz + // Compute wz and zwz. gsl_vector_memcpy (wz, z); gsl_vector_mul (wz, w); @@ -2560,7 +2549,7 @@ void CalcCIss(const gsl_matrix *Xz, const gsl_matrix *XWz, const gsl_matrix *XtX gsl_vector_set (zz, vec_cat[i], d); } - //compute wz, ve and Xz_pve + // Compute wz, ve and Xz_pve. gsl_vector_set_zero (Xz_pve); s_pve=0; s_snp=0; for (size_t i=0; isize; i++) { d=v_pve[vec_cat[i]]/gsl_vector_get(s_vec, vec_cat[i]); gsl_vector_set (w_pve, i, d); } - //compute Vq (in qvar_mat) + // Compute Vq (in qvar_mat). s0=1-s_pve; for (size_t i=0; i