From 93a7a2adb03f61e80badf6a5004fa4850dbb7d48 Mon Sep 17 00:00:00 2001 From: Peter Carbonetto Date: Wed, 7 Jun 2017 23:23:35 -0500 Subject: Removed FORCE_FLOAT from a few more files. --- src/bslmm.cpp | 839 +++++++++++++++++++++++++++++++--------------------------- 1 file changed, 450 insertions(+), 389 deletions(-) (limited to 'src/bslmm.cpp') diff --git a/src/bslmm.cpp b/src/bslmm.cpp index d295fd8..92762e2 100644 --- a/src/bslmm.cpp +++ b/src/bslmm.cpp @@ -14,7 +14,7 @@ You should have received a copy of the GNU General Public License along with this program. If not, see . - */ +*/ #include #include @@ -47,8 +47,7 @@ using namespace std; -void BSLMM::CopyFromParam (PARAM &cPar) -{ +void BSLMM::CopyFromParam (PARAM &cPar) { a_mode=cPar.a_mode; d_pace=cPar.d_pace; @@ -101,9 +100,7 @@ void BSLMM::CopyFromParam (PARAM &cPar) return; } - -void BSLMM::CopyToParam (PARAM &cPar) -{ +void BSLMM::CopyToParam (PARAM &cPar) { cPar.time_UtZ=time_UtZ; cPar.time_Omega=time_Omega; cPar.time_Proposal=time_Proposal; @@ -115,16 +112,16 @@ void BSLMM::CopyToParam (PARAM &cPar) return; } - - -void BSLMM::WriteBV (const gsl_vector *bv) -{ +void BSLMM::WriteBV (const gsl_vector *bv) { string file_str; file_str=path_out+"/"+file_out; file_str+=".bv.txt"; ofstream outfile (file_str.c_str(), ofstream::out); - if (!outfile) {cout<<"error writing file: "< > &beta_g, const gsl_vector *alpha, const size_t w) -{ +void BSLMM::WriteParam (vector > &beta_g, + const gsl_vector *alpha, const size_t w) { string file_str; file_str=path_out+"/"+file_out; file_str+=".param.txt"; ofstream outfile (file_str.c_str(), ofstream::out); - if (!outfile) {cout<<"error writing file: "< > &beta_g, const gsl_vector if (indicator_snp[i]==0) {continue;} outfile< > &beta_g, const gsl_vector return; } - -void BSLMM::WriteParam (const gsl_vector *alpha) -{ +void BSLMM::WriteParam (const gsl_vector *alpha) { string file_str; file_str=path_out+"/"+file_out; file_str+=".param.txt"; ofstream outfile (file_str.c_str(), ofstream::out); - if (!outfile) {cout<<"error writing file: "< &rank) -{ +void BSLMM::SetXgamma (gsl_matrix *Xgamma, const gsl_matrix *X, + vector &rank) { size_t pos; for (size_t i=0; i & return; } - - -double BSLMM::CalcPveLM (const gsl_matrix *UtXgamma, const gsl_vector *Uty, const double sigma_a2) -{ +double BSLMM::CalcPveLM (const gsl_matrix *UtXgamma, const gsl_vector *Uty, + const double sigma_a2) { double pve, var_y; gsl_matrix *Omega=gsl_matrix_alloc (UtXgamma->size2, UtXgamma->size2); @@ -333,9 +342,9 @@ double BSLMM::CalcPveLM (const gsl_matrix *UtXgamma, const gsl_vector *Uty, cons return pve; } - -void BSLMM::InitialMCMC (const gsl_matrix *UtX, const gsl_vector *Uty, vector &rank, class HYPBSLMM &cHyp, vector > &pos_loglr) -{ +void BSLMM::InitialMCMC (const gsl_matrix *UtX, const gsl_vector *Uty, + vector &rank, class HYPBSLMM &cHyp, + vector > &pos_loglr) { double q_genome=gsl_cdf_chisq_Qinv(0.05/(double)ns_test, 1); cHyp.n_gamma=0; @@ -362,7 +371,8 @@ void BSLMM::InitialMCMC (const gsl_matrix *UtX, const gsl_vector *Uty, vectorlogp_max) {cHyp.logp=logp_max;} - -// if (fix_sigma>=0) { -// fix_sigma=cHyp.h; -// rho_max=1-cHyp.h; -// cHyp.rho=rho_max/2.0; -// } - - //Initial for grid sampling: -// cHyp.h=0.225; -// cHyp.rho=1.0; -// cHyp.logp=-4.835429; - cout<<"initial value of h = "<10) {break;} + if (t%d_pace==0 || t==total_step-1) { + ProgressBar ("Running MCMC ", t, total_step-1, + (double)n_accept/(double)(t*n_mh+1)); + } if (a_mode==13) { SampleZ (y, z_hat, z); @@ -969,60 +948,75 @@ void BSLMM::MCMC (const gsl_matrix *U, const gsl_matrix *UtX, const gsl_vector * time_start=clock(); gsl_blas_dgemv (CblasTrans, 1.0, U, z, 0.0, Utz); - time_UtZ+=(clock()-time_start)/(double(CLOCKS_PER_SEC)*60.0); + time_UtZ+=(clock()-time_start)/ + (double(CLOCKS_PER_SEC)*60.0); - //First proposal - if (cHyp_old.n_gamma==0 || cHyp_old.rho==0) { - logPost_old=CalcPosterior(Utz, K_eval, Utu_old, alpha_old, cHyp_old); + // First proposal. + if (cHyp_old.n_gamma==0 || cHyp_old.rho==0) { + logPost_old= + CalcPosterior(Utz, K_eval, Utu_old, + alpha_old, cHyp_old); beta_old.clear(); for (size_t i=0; isize; ++i) { - beta_old.push_back(gsl_vector_get(beta, i)); + beta_old.push_back(gsl_vector_get(beta, i)); } gsl_matrix_free (UtXgamma); gsl_vector_free (beta); } } - //MH steps + // M-H steps. for (size_t i=0; i=0) { -// cHyp_new.h=fix_sigma/(1-cHyp_new.rho); -// } - if (cHyp_new.n_gamma==0 || cHyp_new.rho==0) { - logPost_new=CalcPosterior(Utz, K_eval, Utu_new, alpha_new, cHyp_new); + logPost_new=CalcPosterior(Utz, K_eval, Utu_new, + alpha_new, cHyp_new); beta_new.clear(); for (size_t i=0; isize; ++i) { - beta_new.push_back(gsl_vector_get(beta, i)); + beta_new.push_back(gsl_vector_get(beta, i)); } gsl_matrix_free (UtXgamma); gsl_vector_free (beta); @@ -1030,17 +1024,20 @@ void BSLMM::MCMC (const gsl_matrix *U, const gsl_matrix *UtX, const gsl_vector * logMHratio+=logPost_new-logPost_old; - if (logMHratio>0 || log(gsl_rng_uniform(gsl_r))0 || + log(gsl_rng_uniform(gsl_r))size2); gsl_vector *H_eval=gsl_vector_alloc (Uty->size); gsl_vector *bv=gsl_vector_alloc (Uty->size); @@ -1169,7 +1173,8 @@ void BSLMM::RidgeR(const gsl_matrix *U, const gsl_matrix *UtX, const gsl_vector gsl_vector_memcpy (bv, Uty); gsl_vector_div (bv, H_eval); - gsl_blas_dgemv (CblasTrans, lambda/(double)UtX->size2, UtX, bv, 0.0, beta); + gsl_blas_dgemv (CblasTrans, lambda/(double)UtX->size2, + UtX, bv, 0.0, beta); gsl_vector_add_constant (H_eval, -1.0); gsl_vector_mul (H_eval, bv); gsl_blas_dgemv (CblasNoTrans, 1.0, U, H_eval, 0.0, bv); @@ -1183,28 +1188,13 @@ void BSLMM::RidgeR(const gsl_matrix *U, const gsl_matrix *UtX, const gsl_vector return; } - - - - - - - - - - - - - - - - -//below fits MCMC for rho=1 -void BSLMM::CalcXtX (const gsl_matrix *X, const gsl_vector *y, const size_t s_size, gsl_matrix *XtX, gsl_vector *Xty) -{ +// Below fits MCMC for rho=1. +void BSLMM::CalcXtX (const gsl_matrix *X, const gsl_vector *y, + const size_t s_size, gsl_matrix *XtX, gsl_vector *Xty) { time_t time_start=clock(); - gsl_matrix_const_view X_sub=gsl_matrix_const_submatrix(X, 0, 0, X->size1, s_size); + gsl_matrix_const_view X_sub=gsl_matrix_const_submatrix(X, 0, 0, X->size1, + s_size); gsl_matrix_view XtX_sub=gsl_matrix_submatrix(XtX, 0, 0, s_size, s_size); gsl_vector_view Xty_sub=gsl_vector_subvector(Xty, 0, s_size); @@ -1217,29 +1207,34 @@ void BSLMM::CalcXtX (const gsl_matrix *X, const gsl_vector *y, const size_t s_si return; } - -void BSLMM::SetXgamma (const gsl_matrix *X, const gsl_matrix *X_old, const gsl_matrix *XtX_old, const gsl_vector *Xty_old, const gsl_vector *y, const vector &rank_old, const vector &rank_new, gsl_matrix *X_new, gsl_matrix *XtX_new, gsl_vector *Xty_new) -{ +void BSLMM::SetXgamma (const gsl_matrix *X, const gsl_matrix *X_old, + const gsl_matrix *XtX_old, const gsl_vector *Xty_old, + const gsl_vector *y, const vector &rank_old, + const vector &rank_new, gsl_matrix *X_new, + gsl_matrix *XtX_new, gsl_vector *Xty_new) { double d; - //rank_old and rank_new are sorted already inside PorposeGamma - //calculate vectors rank_remove and rank_add - // size_t v_size=max(rank_old.size(), rank_new.size()); - //make sure that v_size is larger than repeat + // rank_old and rank_new are sorted already inside PorposeGamma + // calculate vectors rank_remove and rank_add. + // make sure that v_size is larger than repeat. size_t v_size=20; - vector rank_remove(v_size), rank_add(v_size), rank_union(s_max+v_size); + vector rank_remove(v_size), rank_add(v_size), + rank_union(s_max+v_size); vector::iterator it; - it=set_difference (rank_old.begin(), rank_old.end(), rank_new.begin(), rank_new.end(), rank_remove.begin()); + it=set_difference(rank_old.begin(), rank_old.end(), rank_new.begin(), + rank_new.end(), rank_remove.begin()); rank_remove.resize(it-rank_remove.begin()); - it=set_difference (rank_new.begin(), rank_new.end(), rank_old.begin(), rank_old.end(), rank_add.begin()); + it=set_difference (rank_new.begin(), rank_new.end(), rank_old.begin(), + rank_old.end(), rank_add.begin()); rank_add.resize(it-rank_add.begin()); - it=set_union (rank_new.begin(), rank_new.end(), rank_old.begin(), rank_old.end(), rank_union.begin()); + it=set_union (rank_new.begin(), rank_new.end(), rank_old.begin(), + rank_old.end(), rank_union.begin()); rank_union.resize(it-rank_union.begin()); - //map rank_remove and rank_add + // Map rank_remove and rank_add. map mapRank2in_remove, mapRank2in_add; for (size_t i=0; isize1, rank_old.size()); - gsl_matrix_const_view XtXold_sub=gsl_matrix_const_submatrix(XtX_old, 0, 0, rank_old.size(), rank_old.size()); - gsl_vector_const_view Xtyold_sub=gsl_vector_const_subvector(Xty_old, 0, rank_old.size()); - - gsl_matrix_view Xnew_sub=gsl_matrix_submatrix(X_new, 0, 0, X_new->size1, rank_new.size()); - gsl_matrix_view XtXnew_sub=gsl_matrix_submatrix(XtX_new, 0, 0, rank_new.size(), rank_new.size()); - gsl_vector_view Xtynew_sub=gsl_vector_subvector(Xty_new, 0, rank_new.size()); - - //get X_new and calculate XtX_new + // Obtain the subset of matrix/vector. + gsl_matrix_const_view Xold_sub= + gsl_matrix_const_submatrix(X_old, 0, 0, X_old->size1, rank_old.size()); + gsl_matrix_const_view XtXold_sub= + gsl_matrix_const_submatrix(XtX_old, 0, 0, rank_old.size(), + rank_old.size()); + gsl_vector_const_view Xtyold_sub= + gsl_vector_const_subvector(Xty_old, 0, rank_old.size()); + + gsl_matrix_view Xnew_sub= + gsl_matrix_submatrix(X_new, 0, 0, X_new->size1, rank_new.size()); + gsl_matrix_view XtXnew_sub= + gsl_matrix_submatrix(XtX_new, 0, 0, rank_new.size(), rank_new.size()); + gsl_vector_view Xtynew_sub= + gsl_vector_subvector(Xty_new, 0, rank_new.size()); + + // Get X_new and calculate XtX_new. if (rank_remove.size()==0 && rank_add.size()==0) { gsl_matrix_memcpy(&Xnew_sub.matrix, &Xold_sub.matrix); gsl_matrix_memcpy(&XtXnew_sub.matrix, &XtXold_sub.matrix); @@ -1295,13 +1297,13 @@ void BSLMM::SetXgamma (const gsl_matrix *X, const gsl_matrix *X_old, const gsl_m gsl_matrix *XtX_ao=gsl_matrix_alloc(X_add->size2, X_old->size2); gsl_vector *Xty_add=gsl_vector_alloc(X_add->size2); - //get X_add + // Get X_add. SetXgamma (X_add, X, rank_add); - //get t(X_add)X_add and t(X_add)X_temp + // Get t(X_add)X_add and t(X_add)X_temp. clock_t time_start=clock(); - //somehow the lapack_dgemm does not work here + // Somehow the lapack_dgemm does not work here. gsl_blas_dgemm (CblasTrans, CblasNoTrans, 1.0, X_add, X_add, 0.0, XtX_aa); gsl_blas_dgemm (CblasTrans, CblasNoTrans, 1.0, X_add, X_old, @@ -1310,18 +1312,26 @@ void BSLMM::SetXgamma (const gsl_matrix *X, const gsl_matrix *X_old, const gsl_m time_Omega+=(clock()-time_start)/(double(CLOCKS_PER_SEC)*60.0); - //save to X_new, XtX_new and Xty_new + // Save to X_new, XtX_new and Xty_new. i_old=0; i_new=0; i_add=0; for (size_t i=0; isize1, s_size); - gsl_matrix_const_view XtX_sub=gsl_matrix_const_submatrix (XtX, 0, 0, s_size, s_size); - gsl_vector_const_view Xty_sub=gsl_vector_const_subvector (Xty, 0, s_size); + gsl_matrix_const_view Xgamma_sub= + gsl_matrix_const_submatrix (Xgamma, 0, 0, Xgamma->size1, s_size); + gsl_matrix_const_view XtX_sub= + gsl_matrix_const_submatrix (XtX, 0, 0, s_size, s_size); + gsl_vector_const_view Xty_sub= + gsl_vector_const_subvector (Xty, 0, s_size); gsl_matrix *Omega=gsl_matrix_alloc (s_size, s_size); gsl_matrix *M_temp=gsl_matrix_alloc (s_size, s_size); @@ -1411,38 +1431,42 @@ double BSLMM::CalcPosterior (const gsl_matrix *Xgamma, const gsl_matrix *XtX, co gsl_vector_memcpy (Xty_temp, &Xty_sub.vector); - //calculate Omega + // Calculate Omega. gsl_matrix_memcpy (Omega, &XtX_sub.matrix); gsl_matrix_scale (Omega, sigma_a2); gsl_matrix_set_identity (M_temp); gsl_matrix_add (Omega, M_temp); - //calculate beta_hat + // Calculate beta_hat. logdet_O=CholeskySolve(Omega, Xty_temp, beta_hat); gsl_vector_scale (beta_hat, sigma_a2); gsl_blas_ddot (Xty_temp, beta_hat, &d); P_yy-=d; - //sample tau + // Sample tau. double tau=1.0; - if (a_mode==11) {tau =gsl_ran_gamma (gsl_r, (double)ni_test/2.0, 2.0/P_yy); } + if (a_mode==11) { + tau = gsl_ran_gamma (gsl_r, (double)ni_test/2.0, 2.0/P_yy); + } - //sample beta + // Sample beta. for (size_t i=0; itm_hour%24*3600+ptm->tm_min*60+ptm->tm_sec); + randseed = (unsigned) (ptm->tm_hour%24*3600+ + ptm->tm_min*60+ptm->tm_sec); } gsl_r = gsl_rng_alloc(gslType); gsl_rng_set(gsl_r, randseed); @@ -1569,7 +1593,7 @@ void BSLMM::MCMC (const gsl_matrix *X, const gsl_vector *y) { gsl_t=gsl_ran_discrete_preproc (ns_test, p_gamma); - //initial parameters + // Initial parameters. InitialMCMC (X, z, rank_old, cHyp_old, pos_loglr); cHyp_initial=cHyp_old; @@ -1580,10 +1604,12 @@ void BSLMM::MCMC (const gsl_matrix *X, const gsl_vector *y) { else { SetXgamma (Xgamma_old, X, rank_old); CalcXtX (Xgamma_old, z, rank_old.size(), XtX_old, Xtz_old); - logPost_old=CalcPosterior (Xgamma_old, XtX_old, Xtz_old, ztz, rank_old.size(), Xb_old, beta_old, cHyp_old); + logPost_old=CalcPosterior (Xgamma_old, XtX_old, Xtz_old, ztz, + rank_old.size(), Xb_old, beta_old, + cHyp_old); } - //calculate centered z_hat, and pve + // Calculate centered z_hat, and pve. if (a_mode==13) { if (cHyp_old.n_gamma==0) { CalcCC_PVEnZ (z_hat, cHyp_old); @@ -1593,65 +1619,94 @@ void BSLMM::MCMC (const gsl_matrix *X, const gsl_vector *y) { } } - //start MCMC + // Start MCMC. int accept; size_t total_step=w_step+s_step; size_t w=0, w_col, pos; size_t repeat=0; for (size_t t=0; t10) {break;} + if (t%d_pace==0 || t==total_step-1) { + ProgressBar ("Running MCMC ", t, total_step-1, + (double)n_accept/(double)(t*n_mh+1)); + } + if (a_mode==13) { SampleZ (y, z_hat, z); mean_z=CenterVector (z); gsl_blas_ddot(z,z,&ztz); - //First proposal + // First proposal. if (cHyp_old.n_gamma==0) { logPost_old=CalcPosterior (ztz, cHyp_old); } else { - gsl_matrix_view Xold_sub=gsl_matrix_submatrix(Xgamma_old, 0, 0, ni_test, rank_old.size()); - gsl_vector_view Xtz_sub=gsl_vector_subvector(Xtz_old, 0, rank_old.size()); - gsl_blas_dgemv (CblasTrans, 1.0, &Xold_sub.matrix, z, 0.0, &Xtz_sub.vector); - logPost_old=CalcPosterior (Xgamma_old, XtX_old, Xtz_old, ztz, rank_old.size(), Xb_old, beta_old, cHyp_old); + gsl_matrix_view Xold_sub= + gsl_matrix_submatrix(Xgamma_old, 0, 0, ni_test, + rank_old.size()); + gsl_vector_view Xtz_sub= + gsl_vector_subvector(Xtz_old, 0, rank_old.size()); + gsl_blas_dgemv (CblasTrans, 1.0, &Xold_sub.matrix, + z, 0.0, &Xtz_sub.vector); + logPost_old= + CalcPosterior (Xgamma_old, XtX_old, Xtz_old, ztz, + rank_old.size(), Xb_old, beta_old, + cHyp_old); } } - //MH steps + // M-H steps. for (size_t i=0; i0 || log(gsl_rng_uniform(gsl_r))0 || + log(gsl_rng_uniform(gsl_r))