From 9aa9f76cb16a0e71fea19dd125761034aaca4fd4 Mon Sep 17 00:00:00 2001 From: xiangzhou Date: Mon, 22 Sep 2014 11:12:46 -0400 Subject: version 0.95alpha --- bslmm.cpp | 1928 ------------------------------------------------------------- 1 file changed, 1928 deletions(-) delete mode 100644 bslmm.cpp (limited to 'bslmm.cpp') diff --git a/bslmm.cpp b/bslmm.cpp deleted file mode 100644 index 55a05ca..0000000 --- a/bslmm.cpp +++ /dev/null @@ -1,1928 +0,0 @@ -/* - Genome-wide Efficient Mixed Model Association (GEMMA) - Copyright (C) 2011 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 - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - This program is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - 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 . - */ - -#include -#include -#include - -#include -#include -#include -#include -#include -#include -#include -#include - -#include "gsl/gsl_vector.h" -#include "gsl/gsl_matrix.h" -#include "gsl/gsl_linalg.h" -#include "gsl/gsl_blas.h" -#include "gsl/gsl_eigen.h" -#include "gsl/gsl_randist.h" -#include "gsl/gsl_cdf.h" -#include "gsl/gsl_roots.h" - - - - -#include "lapack.h" - -#ifdef FORCE_FLOAT -#include "param_float.h" -#include "bslmm_float.h" -#include "lmm_float.h" //for class FUNC_PARAM and MatrixCalcLR -#include "lm_float.h" -#include "mathfunc_float.h" //for function CenterVector -#else -#include "param.h" -#include "bslmm.h" -#include "lmm.h" -#include "lm.h" -#include "mathfunc.h" -#endif - -using namespace std; - - - - -void BSLMM::CopyFromParam (PARAM &cPar) -{ - a_mode=cPar.a_mode; - d_pace=cPar.d_pace; - - file_bfile=cPar.file_bfile; - file_geno=cPar.file_geno; - file_out=cPar.file_out; - path_out=cPar.path_out; - - l_min=cPar.h_min; - l_max=cPar.h_max; - n_region=cPar.n_region; - pve_null=cPar.pve_null; - pheno_mean=cPar.pheno_mean; - - time_UtZ=0.0; - time_Omega=0.0; - n_accept=0; - - h_min=cPar.h_min; - h_max=cPar.h_max; - h_scale=cPar.h_scale; - rho_min=cPar.rho_min; - rho_max=cPar.rho_max; - rho_scale=cPar.rho_scale; - logp_min=cPar.logp_min; - logp_max=cPar.logp_max; - logp_scale=cPar.logp_scale; - - s_min=cPar.s_min; - s_max=cPar.s_max; - w_step=cPar.w_step; - s_step=cPar.s_step; - r_pace=cPar.r_pace; - w_pace=cPar.w_pace; - n_mh=cPar.n_mh; - geo_mean=cPar.geo_mean; - randseed=cPar.randseed; - trace_G=cPar.trace_G; - - ni_total=cPar.ni_total; - ns_total=cPar.ns_total; - ni_test=cPar.ni_test; - ns_test=cPar.ns_test; - n_cvt=cPar.n_cvt; - - indicator_idv=cPar.indicator_idv; - indicator_snp=cPar.indicator_snp; - snpInfo=cPar.snpInfo; - - return; -} - - -void BSLMM::CopyToParam (PARAM &cPar) -{ - cPar.time_UtZ=time_UtZ; - cPar.time_Omega=time_Omega; - cPar.time_Proposal=time_Proposal; - cPar.cHyp_initial=cHyp_initial; - cPar.n_accept=n_accept; - cPar.pheno_mean=pheno_mean; - cPar.randseed=randseed; - - return; -} - - - -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) -{ - 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) -{ - size_t pos; - for (size_t i=0; isize2, UtXgamma->size2); - gsl_vector *Xty=gsl_vector_alloc (UtXgamma->size2); - gsl_vector *OiXty=gsl_vector_alloc (UtXgamma->size2); - - gsl_matrix_set_identity (Omega); - gsl_matrix_scale (Omega, 1.0/sigma_a2); - -#ifdef WITH_LAPACK - lapack_dgemm ((char *)"T", (char *)"N", 1.0, UtXgamma, UtXgamma, 1.0, Omega); -#else - gsl_blas_dgemm (CblasTrans, CblasNoTrans, 1.0, UtXgamma, UtXgamma, 1.0, Omega); -#endif - gsl_blas_dgemv (CblasTrans, 1.0, UtXgamma, Uty, 0.0, Xty); - - CholeskySolve(Omega, Xty, OiXty); - - gsl_blas_ddot (Xty, OiXty, &pve); - gsl_blas_ddot (Uty, Uty, &var_y); - - pve/=var_y; - - gsl_matrix_free (Omega); - gsl_vector_free (Xty); - gsl_vector_free (OiXty); - - return pve; -} - - -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; - for (size_t i=0; iq_genome) {cHyp.n_gamma++;} - } - if (cHyp.n_gamma<10) {cHyp.n_gamma=10;} - - if (cHyp.n_gamma>s_max) {cHyp.n_gamma=s_max;} - if (cHyp.n_gamma1.0) {cHyp.rho=1.0;} - - if (cHyp.hh_max) {cHyp.h=h_max;} - if (cHyp.rhorho_max) {cHyp.rho=rho_max;} - if (cHyp.logplogp_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 = "<size; i++) - { - d=gsl_ran_gaussian(gsl_r, 1); - gsl_vector_set(beta, i, d); - } - gsl_blas_dtrsv(CblasUpper, CblasNoTrans, CblasNonUnit, Omega, beta); - - - //it compuates inv(L^T(Omega)) %*% beta; - gsl_vector_scale(beta, sqrt(sigma_a2/tau)); - gsl_vector_add(beta, beta_hat); - gsl_blas_dgemv (CblasNoTrans, 1.0, UtXgamma, beta, 0.0, UtXb); - - //sample alpha - gsl_vector_memcpy (alpha_prime, Uty); - gsl_vector_sub (alpha_prime, UtXb); - gsl_vector_mul (alpha_prime, weight_Hi); - gsl_vector_scale (alpha_prime, sigma_b2); - - //sample u - gsl_vector_memcpy (Utu, alpha_prime); - gsl_vector_mul (Utu, K_eval); - - if (a_mode==11) {gsl_vector_scale (Utu_rand, sqrt(1.0/tau));} - gsl_vector_add (Utu, Utu_rand); - - - //for quantitative traits, calculate pve and pge - if (a_mode==11) { - gsl_blas_ddot (UtXb, UtXb, &d); - cHyp.pge=d/(double)ni_test; - - gsl_blas_ddot (Utu, Utu, &d); - cHyp.pve=cHyp.pge+d/(double)ni_test; - - if (cHyp.pve==0) {cHyp.pge=0.0;} - else {cHyp.pge/=cHyp.pve;} - cHyp.pve/=cHyp.pve+1.0/tau; - } - - - gsl_matrix_free (UtXgamma_eval); - gsl_matrix_free (Omega); - gsl_vector_free (XtHiy); - gsl_vector_free (beta_hat); - gsl_vector_free (Utu_rand); - gsl_vector_free (weight_Hi); - - logpost=-0.5*logdet_H-0.5*logdet_O; - if (a_mode==11) {logpost-=0.5*(double)ni_test*log(P_yy);} - else {logpost-=0.5*P_yy;} -// else {logpost+=-0.5*P_yy*tau+0.5*(double)ni_test*log(tau);} - logpost+=((double)cHyp.n_gamma-1.0)*cHyp.logp+((double)ns_test-(double)cHyp.n_gamma)*log(1.0-exp(cHyp.logp)); - - return logpost; -} - - - -//calculate pve and pge, and calculate z_hat for case-control data -void BSLMM::CalcCC_PVEnZ (const gsl_matrix *U, const gsl_vector *Utu, gsl_vector *z_hat, class HYPBSLMM &cHyp) -{ - double d; - - gsl_blas_ddot (Utu, Utu, &d); - cHyp.pve=d/(double)ni_test; - - gsl_blas_dgemv (CblasNoTrans, 1.0, U, Utu, 0.0, z_hat); - - cHyp.pve/=cHyp.pve+1.0; - cHyp.pge=0.0; - - return; -} - - -//calculate pve and pge, and calculate z_hat for case-control data -void BSLMM::CalcCC_PVEnZ (const gsl_matrix *U, const gsl_vector *UtXb, const gsl_vector *Utu, gsl_vector *z_hat, class HYPBSLMM &cHyp) -{ - double d; - gsl_vector *UtXbU=gsl_vector_alloc (Utu->size); - - gsl_blas_ddot (UtXb, UtXb, &d); - cHyp.pge=d/(double)ni_test; - - gsl_blas_ddot (Utu, Utu, &d); - cHyp.pve=cHyp.pge+d/(double)ni_test; - - gsl_vector_memcpy (UtXbU, Utu); - gsl_vector_add (UtXbU, UtXb); - gsl_blas_dgemv (CblasNoTrans, 1.0, U, UtXbU, 0.0, z_hat); - - if (cHyp.pve==0) {cHyp.pge=0.0;} - else {cHyp.pge/=cHyp.pve;} - - cHyp.pve/=cHyp.pve+1.0; - - gsl_vector_free(UtXbU); - return; -} - - - - -void BSLMM::SampleZ (const gsl_vector *y, const gsl_vector *z_hat, gsl_vector *z) -{ - double d1, d2, z_rand=0.0; - for (size_t i=0; isize; ++i) { - d1=gsl_vector_get (y, i); - d2=gsl_vector_get (z_hat, i); - //y is centerred for case control studies - if (d1<=0.0) { - //control, right truncated - do { - z_rand=d2+gsl_ran_gaussian(gsl_r, 1.0); - } while (z_rand>0.0); - } - else { - do { - z_rand=d2+gsl_ran_gaussian(gsl_r, 1.0); - } while (z_rand<0.0); - } - - gsl_vector_set (z, i, z_rand); - } - - return; -} - - - - - -double BSLMM::ProposeHnRho (const class HYPBSLMM &cHyp_old, class HYPBSLMM &cHyp_new, const size_t &repeat) -{ - - double h=cHyp_old.h, rho=cHyp_old.rho; - - double d_h=(h_max-h_min)*h_scale, d_rho=(rho_max-rho_min)*rho_scale; - - for (size_t i=0; ih_max) {h=2*h_max-h;} - - rho=rho+(gsl_rng_uniform(gsl_r)-0.5)*d_rho; - if (rhorho_max) {rho=2*rho_max-rho;} - } - /* - //Grid Sampling - for (size_t i=0; ih_max) {h=h_min;} - } - - for (size_t i=0; irho_max) {rho=rho_min;} - } - */ - cHyp_new.h=h; - cHyp_new.rho=rho; - return 0.0; -} - - -double BSLMM::ProposePi (const class HYPBSLMM &cHyp_old, class HYPBSLMM &cHyp_new, const size_t &repeat) -{ - double logp_old=cHyp_old.logp, logp_new=cHyp_old.logp; - double log_ratio=0.0; - - double d_logp=min(0.1, (logp_max-logp_min)*logp_scale); - - for (size_t i=0; ilogp_max) {logp_new=2*logp_max-logp_new;} - - log_ratio+=logp_new-logp_old; - logp_old=logp_new; - } - /* - //Grid Sampling - for (size_t i=0; ilogp_max) {logp_new=logp_min;} - - log_ratio+=logp_new-logp_old; - logp_old=logp_new; - } - */ - cHyp_new.logp=logp_new; - - return log_ratio; -} - -bool comp_vec (size_t a, size_t b) -{ - return (a < b); -} - - -double BSLMM::ProposeGamma (const vector &rank_old, vector &rank_new, const double *p_gamma, const class HYPBSLMM &cHyp_old, class HYPBSLMM &cHyp_new, const size_t &repeat) -{ - map mapRank2in; - size_t r; - double unif, logp=0.0; - int flag_gamma; - size_t r_add, r_remove, col_id; - - rank_new.clear(); - if (cHyp_old.n_gamma!=rank_old.size()) {cout<<"size wrong"<=0.40 && unif < 0.80 && cHyp_new.n_gamma>s_min) {flag_gamma=2;} - else if (unif>=0.80 && cHyp_new.n_gamma>0 && cHyp_new.n_gamma a, pair b) -{ - return (a.second > b.second); -} - - - - - - - -//if a_mode==13 then Uty==y -void BSLMM::MCMC (const gsl_matrix *U, const gsl_matrix *UtX, const gsl_vector *Uty, const gsl_vector *K_eval, const gsl_vector *y) { - clock_t time_start; - - class HYPBSLMM cHyp_old, cHyp_new; - - gsl_matrix *Result_hyp=gsl_matrix_alloc (w_pace, 6); - gsl_matrix *Result_gamma=gsl_matrix_alloc (w_pace, s_max); - - gsl_vector *alpha_prime=gsl_vector_alloc (ni_test); - gsl_vector *alpha_new=gsl_vector_alloc (ni_test); - gsl_vector *alpha_old=gsl_vector_alloc (ni_test); - gsl_vector *Utu=gsl_vector_alloc (ni_test); - gsl_vector *Utu_new=gsl_vector_alloc (ni_test); - gsl_vector *Utu_old=gsl_vector_alloc (ni_test); - - gsl_vector *UtXb_new=gsl_vector_alloc (ni_test); - gsl_vector *UtXb_old=gsl_vector_alloc (ni_test); - - gsl_vector *z_hat=gsl_vector_alloc (ni_test); - gsl_vector *z=gsl_vector_alloc (ni_test); - gsl_vector *Utz=gsl_vector_alloc (ni_test); - - gsl_vector_memcpy (Utz, Uty); - - double logPost_new, logPost_old; - double logMHratio; - double mean_z=0.0; - - gsl_matrix_set_zero (Result_gamma); - gsl_vector_set_zero (Utu); - gsl_vector_set_zero (alpha_prime); - if (a_mode==13) { - pheno_mean=0.0; - } - - vector > beta_g; - for (size_t i=0; i rank_new, rank_old; - vector beta_new, beta_old; - - vector > pos_loglr; - - time_start=clock(); - MatrixCalcLR (U, UtX, Utz, K_eval, l_min, l_max, n_region, pos_loglr); - time_Proposal=(clock()-time_start)/(double(CLOCKS_PER_SEC)*60.0); - - stable_sort (pos_loglr.begin(), pos_loglr.end(), comp_lr); - for (size_t i=0; itm_hour%24*3600+ptm->tm_min*60+ptm->tm_sec); - } - gsl_r = gsl_rng_alloc(gslType); - gsl_rng_set(gsl_r, randseed); - - double *p_gamma = new double[ns_test]; - CalcPgamma (p_gamma); - - gsl_t=gsl_ran_discrete_preproc (ns_test, p_gamma); - - //initial parameters - InitialMCMC (UtX, Utz, rank_old, cHyp_old, pos_loglr); -// if (fix_sigma>=0) { -// rho_max=1-fix_sigma; -// cHyp_old.h=fix_sigma/(1-cHyp_old.rho); -// } - - cHyp_initial=cHyp_old; - - 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)); - } - gsl_matrix_free (UtXgamma); - gsl_vector_free (beta); - } - - //calculate centered z_hat, and pve - if (a_mode==13) { - time_start=clock(); - if (cHyp_old.n_gamma==0 || cHyp_old.rho==0) { - CalcCC_PVEnZ (U, Utu_old, z_hat, cHyp_old); - } - else { - CalcCC_PVEnZ (U, UtXb_old, Utu_old, z_hat, cHyp_old); - } - time_UtZ+=(clock()-time_start)/(double(CLOCKS_PER_SEC)*60.0); - } - - //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 (a_mode==13) { - SampleZ (y, z_hat, z); - mean_z=CenterVector (z); - - 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); - - //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)); - } - gsl_matrix_free (UtXgamma); - gsl_vector_free (beta); - } - } - - //MH 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); - beta_new.clear(); - for (size_t i=0; isize; ++i) { - beta_new.push_back(gsl_vector_get(beta, i)); - } - gsl_matrix_free (UtXgamma); - gsl_vector_free (beta); - } - - logMHratio+=logPost_new-logPost_old; - - if (logMHratio>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); - - gsl_vector_memcpy (H_eval, eval); - gsl_vector_scale (H_eval, lambda); - gsl_vector_add_constant (H_eval, 1.0); - - 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_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); - - WriteParam (beta); - WriteBV(bv); - - gsl_vector_free (H_eval); - gsl_vector_free (beta); - gsl_vector_free (bv); - - 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) -{ - time_t time_start=clock(); - 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); - -#ifdef WITH_LAPACK - lapack_dgemm ((char *)"T", (char *)"N", 1.0, &X_sub.matrix, &X_sub.matrix, 0.0, &XtX_sub.matrix); -#else - gsl_blas_dgemm (CblasTrans, CblasNoTrans, 1.0, &X_sub.matrix, &X_sub.matrix, 0.0, &XtX_sub.matrix); -#endif - gsl_blas_dgemv(CblasTrans, 1.0, &X_sub.matrix, y, 0.0, &Xty_sub.vector); - - time_Omega+=(clock()-time_start)/(double(CLOCKS_PER_SEC)*60.0); - - 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) -{ - 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 - size_t v_size=20; - 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()); - 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()); - 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()); - rank_union.resize(it-rank_union.begin()); - - //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 - /* - 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); - gsl_vector_memcpy(&Xtynew_sub.vector, &Xtyold_sub.vector); - } else { - gsl_matrix *X_temp=gsl_matrix_alloc(X_old->size1, rank_old.size()-rank_remove.size() ); - gsl_matrix *XtX_temp=gsl_matrix_alloc(X_temp->size2, X_temp->size2); - gsl_vector *Xty_temp=gsl_vector_alloc(X_temp->size2); - - if (rank_remove.size()==0) { - gsl_matrix_memcpy (X_temp, &Xold_sub.matrix); - gsl_matrix_memcpy (XtX_temp, &XtXold_sub.matrix); - gsl_vector_memcpy (Xty_temp, &Xtyold_sub.vector); - } else { - size_t i_temp=0, j_temp; - for (size_t i=0; isize1, rank_add.size() ); - gsl_matrix *XtX_aa=gsl_matrix_alloc(X_add->size2, X_add->size2); - gsl_matrix *XtX_at=gsl_matrix_alloc(X_add->size2, X_temp->size2); - gsl_vector *Xty_add=gsl_vector_alloc(X_add->size2); - - //get X_add - SetXgamma (X_add, X, rank_add); - - //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 - //#ifdef WITH_LAPACK - //lapack_dgemm ((char *)"T", (char *)"N", 1.0, X_add, X_add, 0.0, XtX_aa); - //lapack_dgemm ((char *)"T", (char *)"N", 1.0, X_add, X_temp, 0.0, XtX_at); - - //#else - 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_temp, 0.0, XtX_at); - //#endif - gsl_blas_dgemv(CblasTrans, 1.0, X_add, y, 0.0, Xty_add); - - time_Omega+=(clock()-time_start)/(double(CLOCKS_PER_SEC)*60.0); - - //save to X_new, XtX_new and Xty_new - size_t i_temp=0, j_temp, i_flag=0, j_flag=0; - for (size_t i=0; isize1, rank_add.size() ); - gsl_matrix *XtX_aa=gsl_matrix_alloc(X_add->size2, X_add->size2); - 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 - SetXgamma (X_add, X, rank_add); - - //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 - //#ifdef WITH_LAPACK - //lapack_dgemm ((char *)"T", (char *)"N", 1.0, X_add, X_add, 0.0, XtX_aa); - //lapack_dgemm ((char *)"T", (char *)"N", 1.0, X_add, X_temp, 0.0, XtX_at); - - //#else - 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, 0.0, XtX_ao); - //#endif - gsl_blas_dgemv(CblasTrans, 1.0, X_add, y, 0.0, Xty_add); - - time_Omega+=(clock()-time_start)/(double(CLOCKS_PER_SEC)*60.0); - - //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 *Omega=gsl_matrix_alloc (s_size, s_size); - gsl_matrix *M_temp=gsl_matrix_alloc (s_size, s_size); - gsl_vector *beta_hat=gsl_vector_alloc (s_size); - gsl_vector *Xty_temp=gsl_vector_alloc (s_size); - - gsl_vector_memcpy (Xty_temp, &Xty_sub.vector); - - //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 - 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 - double tau=1.0; - if (a_mode==11) {tau =gsl_ran_gamma (gsl_r, (double)ni_test/2.0, 2.0/P_yy); } - - //sample beta - for (size_t i=0; i > beta_g; - for (size_t i=0; i rank_new, rank_old; - vector > pos_loglr; - - time_start=clock(); - MatrixCalcLmLR (X, z, pos_loglr); - time_Proposal=(clock()-time_start)/(double(CLOCKS_PER_SEC)*60.0); - - stable_sort (pos_loglr.begin(), pos_loglr.end(), comp_lr); - for (size_t i=0; itm_hour%24*3600+ptm->tm_min*60+ptm->tm_sec); - } - gsl_r = gsl_rng_alloc(gslType); - gsl_rng_set(gsl_r, randseed); - - double *p_gamma = new double[ns_test]; - CalcPgamma (p_gamma); - - gsl_t=gsl_ran_discrete_preproc (ns_test, p_gamma); - - //initial parameters - InitialMCMC (X, z, rank_old, cHyp_old, pos_loglr); - - cHyp_initial=cHyp_old; - - if (cHyp_old.n_gamma==0) { - logPost_old=CalcPosterior (ztz, cHyp_old); - } - 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); - } - - //calculate centered z_hat, and pve - if (a_mode==13) { - if (cHyp_old.n_gamma==0) { - CalcCC_PVEnZ (z_hat, cHyp_old); - } - else { - CalcCC_PVEnZ (Xb_old, z_hat, cHyp_old); - } - } - - //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 (a_mode==13) { - SampleZ (y, z_hat, z); - mean_z=CenterVector (z); - gsl_blas_ddot(z,z,&ztz); - - //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); - } - } - - //MH steps - for (size_t i=0; i0 || log(gsl_rng_uniform(gsl_r))