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))