initial upload, version 0.95alpha

1 files changed, 1927 insertions, 0 deletions

diff --git a/bslmm.cpp b/bslmm.cpp
new file mode 100644
index 0000000..ff9618d
--- /dev/null
+++ b/bslmm.cpp
@@ -0,0 +1,1927 @@
+/*
+ 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 <http://www.gnu.org/licenses/>.
+ */
+
+#include <iostream>
+#include <fstream>
+#include <sstream>
+
+#include <iomanip>
+#include <cmath>
+#include <iostream>
+#include <stdio.h>
+#include <stdlib.h> 
+#include <ctime>
+#include <cstring>
+#include <algorithm>
+
+#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;
+	
+	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="./output/"+file_out;
+	file_str+=".bv.txt";
+
+	ofstream outfile (file_str.c_str(), ofstream::out);
+	if (!outfile) {cout<<"error writing file: "<<file_str.c_str()<<endl; return;}
+	
+	size_t t=0;
+	for (size_t i=0; i<ni_total; ++i) {
+		if (indicator_idv[i]==0) {
+			outfile<<"NA"<<endl;
+		}		
+		else {
+			outfile<<scientific<<setprecision(6)<<gsl_vector_get(bv, t)<<endl;
+			t++;
+		}
+	}		
+	
+	outfile.clear();	
+	outfile.close();	
+	return;
+}
+
+
+
+
+void BSLMM::WriteParam (vector<pair<double, double> > &beta_g, const gsl_vector *alpha, const size_t w) 
+{
+	string file_str;
+	file_str="./output/"+file_out;
+	file_str+=".param.txt";
+
+	ofstream outfile (file_str.c_str(), ofstream::out);
+	if (!outfile) {cout<<"error writing file: "<<file_str.c_str()<<endl; return;}
+	
+	outfile<<"chr"<<"\t"<<"rs"<<"\t"
+			<<"ps"<<"\t"<<"n_miss"<<"\t"<<"alpha"<<"\t"
+			<<"beta"<<"\t"<<"gamma"<<endl;
+	
+	size_t t=0;
+	for (size_t i=0; i<ns_total; ++i) {
+		if (indicator_snp[i]==0) {continue;}		
+		
+		outfile<<snpInfo[i].chr<<"\t"<<snpInfo[i].rs_number<<"\t"
+		<<snpInfo[i].base_position<<"\t"<<snpInfo[i].n_miss<<"\t";	
+				
+		outfile<<scientific<<setprecision(6)<<gsl_vector_get(alpha, t)<<"\t";
+		if (beta_g[t].second!=0) {
+			outfile<<beta_g[t].first/beta_g[t].second<<"\t"<<beta_g[t].second/(double)w<<endl;
+		}
+		else {
+			outfile<<0.0<<"\t"<<0.0<<endl;
+		}
+		t++;
+	}		
+	
+	outfile.clear();	
+	outfile.close();	
+	return;
+}
+
+
+void BSLMM::WriteParam (const gsl_vector *alpha) 
+{
+	string file_str;
+	file_str="./output/"+file_out;
+	file_str+=".param.txt";
+
+	ofstream outfile (file_str.c_str(), ofstream::out);
+	if (!outfile) {cout<<"error writing file: "<<file_str.c_str()<<endl; return;}
+	
+	outfile<<"chr"<<"\t"<<"rs"<<"\t"
+			<<"ps"<<"\t"<<"n_miss"<<"\t"<<"alpha"<<"\t"
+			<<"beta"<<"\t"<<"gamma"<<endl;
+	
+	size_t t=0;
+	for (size_t i=0; i<ns_total; ++i) {
+		if (indicator_snp[i]==0) {continue;}		
+
+		outfile<<snpInfo[i].chr<<"\t"<<snpInfo[i].rs_number<<"\t"
+				<<snpInfo[i].base_position<<"\t"<<snpInfo[i].n_miss<<"\t";				
+		outfile<<scientific<<setprecision(6)<<gsl_vector_get(alpha, t)<<"\t";
+		outfile<<0.0<<"\t"<<0.0<<endl;
+		t++;
+	}		
+	
+	outfile.clear();	
+	outfile.close();	
+	return;
+}
+
+
+void BSLMM::WriteResult (const int flag, const gsl_matrix *Result_hyp, const gsl_matrix *Result_gamma, const size_t w_col) 
+{
+	string file_gamma, file_hyp;
+	file_gamma="./output/"+file_out;
+	file_gamma+=".gamma.txt";
+	file_hyp="./output/"+file_out;
+	file_hyp+=".hyp.txt";
+
+	ofstream outfile_gamma, outfile_hyp;
+		
+	if (flag==0) {
+		outfile_gamma.open (file_gamma.c_str(), ofstream::out);
+		outfile_hyp.open (file_hyp.c_str(), ofstream::out);
+		if (!outfile_gamma) {cout<<"error writing file: "<<file_gamma<<endl; return;}
+		if (!outfile_hyp) {cout<<"error writing file: "<<file_hyp<<endl; return;}
+		
+		outfile_hyp<<"h \t pve \t rho \t pge \t pi \t n_gamma"<<endl;
+		
+		for (size_t i=0; i<s_max; ++i) {
+			outfile_gamma<<"s"<<i<<"\t";
+		}
+		outfile_gamma<<endl;
+	}
+	else {
+		outfile_gamma.open (file_gamma.c_str(), ofstream::app);
+		outfile_hyp.open (file_hyp.c_str(), ofstream::app);
+		if (!outfile_gamma) {cout<<"error writing file: "<<file_gamma<<endl; return;}
+		if (!outfile_hyp) {cout<<"error writing file: "<<file_hyp<<endl; return;}
+		
+		size_t w;
+		if (w_col==0) {w=w_pace;}
+		else {w=w_col;}
+		
+		for (size_t i=0; i<w; ++i) {
+			outfile_hyp<<scientific;
+			for (size_t j=0; j<4; ++j) {
+				outfile_hyp<<setprecision(6)<<gsl_matrix_get (Result_hyp, i, j)<<"\t";
+			}
+			outfile_hyp<<setprecision(6)<<exp(gsl_matrix_get (Result_hyp, i, 4))<<"\t";
+			outfile_hyp<<(int)gsl_matrix_get (Result_hyp, i, 5)<<"\t";
+			outfile_hyp<<endl;
+		}
+		
+		for (size_t i=0; i<w; ++i) {
+			for (size_t j=0; j<s_max; ++j) {
+				outfile_gamma<<(int)gsl_matrix_get (Result_gamma, i, j)<<"\t";
+			}
+			outfile_gamma<<endl;
+		}
+		
+	}
+	
+	outfile_hyp.close();
+	outfile_hyp.clear();
+	outfile_gamma.close();
+	outfile_gamma.clear();	
+	return;
+}
+
+
+
+void BSLMM::CalcPgamma (double *p_gamma)
+{
+	double p, s=0.0;
+	for (size_t i=0; i<ns_test; ++i) {
+		p=0.7*gsl_ran_geometric_pdf (i+1, 1.0/geo_mean)+0.3/(double)ns_test;
+		p_gamma[i]=p;
+		s+=p;
+	}
+	for (size_t i=0; i<ns_test; ++i) {
+		p=p_gamma[i];
+		p_gamma[i]=p/s;
+	}
+	return;
+}
+
+
+
+void BSLMM::SetXgamma (gsl_matrix *Xgamma, const gsl_matrix *X, vector<size_t> &rank)
+{
+	size_t pos;
+	for (size_t i=0; i<rank.size(); ++i) {
+		pos=mapRank2pos[rank[i]];
+		gsl_vector_view Xgamma_col=gsl_matrix_column (Xgamma, i);
+		gsl_vector_const_view X_col=gsl_matrix_const_column (X, pos);
+		gsl_vector_memcpy (&Xgamma_col.vector, &X_col.vector);
+	}
+	
+	return;
+}
+
+
+
+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);
+	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<size_t> &rank, class HYPBSLMM &cHyp, vector<pair<size_t, double> > &pos_loglr)
+{
+	double q_genome=gsl_cdf_chisq_Qinv(0.05/(double)ns_test, 1);
+	
+	cHyp.n_gamma=0;
+	for (size_t i=0; i<pos_loglr.size(); ++i) {
+		if (2.0*pos_loglr[i].second>q_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_gamma<s_min) {cHyp.n_gamma=s_min;}	
+	
+	rank.clear();
+	for (size_t i=0; i<cHyp.n_gamma; ++i) {
+		rank.push_back(i);
+	}
+	
+	cHyp.logp=log((double)cHyp.n_gamma/(double)ns_test);
+	cHyp.h=pve_null; 
+	
+	if (cHyp.logp==0) {cHyp.logp=-0.000001;}
+	if (cHyp.h==0) {cHyp.h=0.1;}
+
+	gsl_matrix *UtXgamma=gsl_matrix_alloc (ni_test, cHyp.n_gamma);
+	SetXgamma (UtXgamma, UtX, rank);
+	double sigma_a2;
+	if (trace_G!=0) {
+	  sigma_a2=cHyp.h*1.0/(trace_G*(1-cHyp.h)*exp(cHyp.logp)*(double)ns_test);
+	} else {
+	  sigma_a2=cHyp.h*1.0/( (1-cHyp.h)*exp(cHyp.logp)*(double)ns_test);
+	}
+	if (sigma_a2==0) {sigma_a2=0.025;}	
+	cHyp.rho=CalcPveLM (UtXgamma, Uty, sigma_a2)/cHyp.h;
+	gsl_matrix_free (UtXgamma);
+	
+	if (cHyp.rho>1.0) {cHyp.rho=1.0;}
+	
+	if (cHyp.h<h_min) {cHyp.h=h_min;}
+	if (cHyp.h>h_max) {cHyp.h=h_max;}
+	if (cHyp.rho<rho_min) {cHyp.rho=rho_min;}
+	if (cHyp.rho>rho_max) {cHyp.rho=rho_max;}
+	if (cHyp.logp<logp_min) {cHyp.logp=logp_min;}
+	if (cHyp.logp>logp_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 = "<<cHyp.h<<endl;
+	cout<<"initial value of rho = "<<cHyp.rho<<endl;
+	cout<<"initial value of pi = "<<exp(cHyp.logp)<<endl;
+	cout<<"initial value of |gamma| = "<<cHyp.n_gamma<<endl;
+	
+	return;
+}
+
+
+
+double BSLMM::CalcPosterior (const gsl_vector *Uty, const gsl_vector *K_eval, gsl_vector *Utu, gsl_vector *alpha_prime, class HYPBSLMM &cHyp)
+{
+	double sigma_b2=cHyp.h*(1.0-cHyp.rho)/(trace_G*(1-cHyp.h));
+	
+	gsl_vector *Utu_rand=gsl_vector_alloc (Uty->size);	
+	gsl_vector *weight_Hi=gsl_vector_alloc (Uty->size);
+	
+	double logpost=0.0;
+	double d, ds, uy, Hi_yy=0, logdet_H=0.0;
+	for (size_t i=0; i<ni_test; ++i) {
+		d=gsl_vector_get (K_eval, i)*sigma_b2;
+		ds=d/(d+1.0);
+		d=1.0/(d+1.0);		
+		gsl_vector_set (weight_Hi, i, d);
+		
+		logdet_H-=log(d);
+		uy=gsl_vector_get (Uty, i);
+		Hi_yy+=d*uy*uy;
+		
+		gsl_vector_set (Utu_rand, i, gsl_ran_gaussian(gsl_r, 1)*sqrt(ds));
+	}
+	
+	//sample tau
+	double tau=1.0;
+	if (a_mode==11) {tau = gsl_ran_gamma (gsl_r, (double)ni_test/2.0,  2.0/Hi_yy); }
+	
+	//sample alpha
+	gsl_vector_memcpy (alpha_prime, Uty);
+	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 ppe
+	if (a_mode==11) {
+		gsl_blas_ddot (Utu, Utu, &d);
+		cHyp.pve=d/(double)ni_test;	
+		cHyp.pve/=cHyp.pve+1.0/tau;
+		cHyp.pge=0.0;	
+	}
+
+	//calculate likelihood
+	logpost=-0.5*logdet_H;
+	if (a_mode==11) {logpost-=0.5*(double)ni_test*log(Hi_yy);}
+	else {logpost-=0.5*Hi_yy;}
+	
+	logpost+=((double)cHyp.n_gamma-1.0)*cHyp.logp+((double)ns_test-(double)cHyp.n_gamma)*log(1-exp(cHyp.logp));
+	
+	gsl_vector_free (Utu_rand);
+	gsl_vector_free (weight_Hi);
+	
+	return logpost;
+}
+
+
+double BSLMM::CalcPosterior (const gsl_matrix *UtXgamma, const gsl_vector *Uty, const gsl_vector *K_eval, gsl_vector *UtXb, gsl_vector *Utu, gsl_vector *alpha_prime, gsl_vector *beta, class HYPBSLMM &cHyp)
+{
+	clock_t time_start;	
+	
+	double sigma_a2=cHyp.h*cHyp.rho/(trace_G*(1-cHyp.h)*exp(cHyp.logp)*(double)ns_test);
+	double sigma_b2=cHyp.h*(1.0-cHyp.rho)/(trace_G*(1-cHyp.h));
+	
+	double logpost=0.0;
+	double d, ds, uy, P_yy=0, logdet_O=0.0, logdet_H=0.0;
+	
+	gsl_matrix *UtXgamma_eval=gsl_matrix_alloc (UtXgamma->size1, UtXgamma->size2);	
+	gsl_matrix *Omega=gsl_matrix_alloc (UtXgamma->size2, UtXgamma->size2);
+	gsl_vector *XtHiy=gsl_vector_alloc (UtXgamma->size2);
+	gsl_vector *beta_hat=gsl_vector_alloc (UtXgamma->size2);
+	gsl_vector *Utu_rand=gsl_vector_alloc (UtXgamma->size1);	
+	gsl_vector *weight_Hi=gsl_vector_alloc (UtXgamma->size1);
+	
+	gsl_matrix_memcpy (UtXgamma_eval, UtXgamma);
+	
+	logdet_H=0.0; P_yy=0.0;
+	for (size_t i=0; i<ni_test; ++i) {
+		gsl_vector_view UtXgamma_row=gsl_matrix_row (UtXgamma_eval, i);
+		d=gsl_vector_get (K_eval, i)*sigma_b2;
+		ds=d/(d+1.0);
+		d=1.0/(d+1.0);
+		gsl_vector_set (weight_Hi, i, d);
+		
+		logdet_H-=log(d);
+		uy=gsl_vector_get (Uty, i);
+		P_yy+=d*uy*uy;
+		gsl_vector_scale (&UtXgamma_row.vector, d);
+		
+		gsl_vector_set (Utu_rand, i, gsl_ran_gaussian(gsl_r, 1)*sqrt(ds));
+	}
+	
+	//calculate Omega
+	gsl_matrix_set_identity (Omega);
+	
+	time_start=clock();
+#ifdef WITH_LAPACK
+	lapack_dgemm ((char *)"T", (char *)"N", sigma_a2, UtXgamma_eval, UtXgamma, 1.0, Omega);
+#else
+	gsl_blas_dgemm (CblasTrans, CblasNoTrans, sigma_a2, UtXgamma_eval, UtXgamma, 1.0, Omega);
+#endif	
+	time_Omega+=(clock()-time_start)/(double(CLOCKS_PER_SEC)*60.0);
+	
+	
+	//calculate beta_hat
+	gsl_blas_dgemv (CblasTrans, 1.0, UtXgamma_eval, Uty, 0.0, XtHiy);	
+
+	logdet_O=CholeskySolve(Omega, XtHiy, beta_hat);
+	
+	gsl_vector_scale (beta_hat, sigma_a2);
+
+	gsl_blas_ddot (XtHiy, 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->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; i<z->size; ++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; i<repeat; ++i) {
+		h=h+(gsl_rng_uniform(gsl_r)-0.5)*d_h;
+		if (h<h_min) {h=2*h_min-h;}
+		if (h>h_max) {h=2*h_max-h;}
+		
+		rho=rho+(gsl_rng_uniform(gsl_r)-0.5)*d_rho;
+		if (rho<rho_min) {rho=2*rho_min-rho;}
+		if (rho>rho_max) {rho=2*rho_max-rho;}
+	}
+	/*
+	//Grid Sampling
+	for (size_t i=0; i<repeat; ++i) {
+		if (gsl_rng_uniform(gsl_r)<0.66) {continue;}
+		h=h+(gsl_rng_uniform_int(gsl_r, 2)-0.5)*0.1;
+		if (h<h_min) {h=h_max;}
+		if (h>h_max) {h=h_min;}
+	}
+	
+	for (size_t i=0; i<repeat; ++i) {
+		if (gsl_rng_uniform(gsl_r)<0.66) {continue;}
+		rho=rho+(gsl_rng_uniform_int(gsl_r, 2)-0.5)*0.1;
+		if (rho<rho_min) {rho=rho_max;}
+		if (rho>rho_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; i<repeat; ++i) {
+		logp_new=logp_old+(gsl_rng_uniform(gsl_r)-0.5)*d_logp;
+		if (logp_new<logp_min) {logp_new=2*logp_min-logp_new;}
+		if (logp_new>logp_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; i<repeat; ++i) {
+		if (gsl_rng_uniform(gsl_r)<0.66) {continue;}
+		logp_new=logp_old+(gsl_rng_uniform_int(gsl_r, 2)-0.5)*0.5*log(10.0);
+		if (logp_new<logp_min) {logp_new=logp_max;}
+		if (logp_new>logp_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<size_t> &rank_old, vector<size_t> &rank_new, const double *p_gamma, const class HYPBSLMM &cHyp_old, class HYPBSLMM &cHyp_new, const size_t &repeat)
+{
+	map<size_t, int> 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"<<endl;}
+	
+	if (cHyp_old.n_gamma!=0) {
+		for (size_t i=0; i<rank_old.size(); ++i) {
+			r=rank_old[i];
+			rank_new.push_back(r);
+			mapRank2in[r]=1;
+		}
+	}
+	cHyp_new.n_gamma=cHyp_old.n_gamma;	
+	
+	for (size_t i=0; i<repeat; ++i) {
+		unif=gsl_rng_uniform(gsl_r); 
+	
+		if (unif < 0.40 && cHyp_new.n_gamma<s_max) {flag_gamma=1;}
+		else if (unif>=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<ns_test) {flag_gamma=3;}
+		else {flag_gamma=4;}
+	
+		if(flag_gamma==1)  {//add a snp; 
+			do {
+				r_add=gsl_ran_discrete (gsl_r, gsl_t);
+			} while (mapRank2in.count(r_add)!=0); 
+		
+			double prob_total=1.0;
+			for (size_t i=0; i<cHyp_new.n_gamma; ++i) {
+				r=rank_new[i];
+				prob_total-=p_gamma[r];
+			}
+
+			mapRank2in[r_add]=1;
+			rank_new.push_back(r_add);
+			cHyp_new.n_gamma++;
+			logp+=-log(p_gamma[r_add]/prob_total)-log((double)cHyp_new.n_gamma);
+		}
+		else if (flag_gamma==2) {//delete a snp;
+			col_id=gsl_rng_uniform_int(gsl_r, cHyp_new.n_gamma);		
+			r_remove=rank_new[col_id];
+		
+			double prob_total=1.0;
+			for (size_t i=0; i<cHyp_new.n_gamma; ++i) {
+				r=rank_new[i];
+				prob_total-=p_gamma[r];
+			}
+			prob_total+=p_gamma[r_remove];
+		
+			mapRank2in.erase(r_remove);
+			rank_new.erase(rank_new.begin()+col_id);
+			logp+=log(p_gamma[r_remove]/prob_total)+log((double)cHyp_new.n_gamma);
+			cHyp_new.n_gamma--;
+		}
+		else if (flag_gamma==3) {//switch a snp;
+			col_id=gsl_rng_uniform_int(gsl_r, cHyp_new.n_gamma);		
+			r_remove=rank_new[col_id];
+		//careful with the proposal
+			do {
+				r_add=gsl_ran_discrete (gsl_r, gsl_t);
+			} while (mapRank2in.count(r_add)!=0); 
+			
+			double prob_total=1.0;
+			for (size_t i=0; i<cHyp_new.n_gamma; ++i) {
+				r=rank_new[i];
+				prob_total-=p_gamma[r];
+			}
+			
+			logp+=log(p_gamma[r_remove]/(prob_total+p_gamma[r_remove]-p_gamma[r_add]) );
+			logp-=log(p_gamma[r_add]/prob_total);
+			
+			mapRank2in.erase(r_remove);
+			mapRank2in[r_add]=1;
+			rank_new.erase(rank_new.begin()+col_id);
+			rank_new.push_back(r_add);
+		}
+		else {logp+=0;}//do not change
+	}
+	
+	stable_sort (rank_new.begin(), rank_new.end(), comp_vec);
+
+	mapRank2in.clear();
+	return logp;
+}
+
+
+
+
+
+
+bool comp_lr (pair<size_t, double> a, pair<size_t, double> 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<pair<double, double> > beta_g;
+	for (size_t i=0; i<ns_test; i++) {
+		beta_g.push_back(make_pair(0.0, 0.0));
+	}
+	
+	vector<size_t> rank_new, rank_old;
+	vector<double> beta_new, beta_old;	
+
+	vector<pair<size_t, double> > 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; i<ns_test; ++i) {
+		mapRank2pos[i]=pos_loglr[i].first;
+	}
+	
+	//calculate proposal distribution for gamma (unnormalized), and set up gsl_r and gsl_t			
+	gsl_rng_env_setup();                
+	const gsl_rng_type * gslType;                                               
+	gslType = gsl_rng_default; 
+	if (randseed<0)
+	{
+		time_t rawtime;
+		time (&rawtime);
+		tm * ptm = gmtime (&rawtime);
+		
+		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);
+	
+	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; i<cHyp_old.n_gamma; ++i) {
+		  beta_old.push_back(0);
+		}	
+	}
+	else {
+		gsl_matrix *UtXgamma=gsl_matrix_alloc (ni_test, cHyp_old.n_gamma);
+		gsl_vector *beta=gsl_vector_alloc (cHyp_old.n_gamma);
+		SetXgamma (UtXgamma, UtX, rank_old);		
+		logPost_old=CalcPosterior(UtXgamma, Utz, K_eval, UtXb_old, Utu_old, alpha_old, beta, cHyp_old);
+	
+		beta_old.clear();
+		for (size_t i=0; i<beta->size; ++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; t<total_step; ++t) {
+		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 (t>10) {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; i<cHyp_old.n_gamma; ++i) {
+				  beta_old.push_back(0);
+				}	
+			}
+			else {
+				gsl_matrix *UtXgamma=gsl_matrix_alloc (ni_test, cHyp_old.n_gamma);
+				gsl_vector *beta=gsl_vector_alloc (cHyp_old.n_gamma);
+				SetXgamma (UtXgamma, UtX, rank_old);
+				logPost_old=CalcPosterior(UtXgamma, Utz, K_eval, UtXb_old, Utu_old, alpha_old, beta, cHyp_old);
+				
+				beta_old.clear();
+				for (size_t i=0; i<beta->size; ++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<n_mh; ++i) {
+			if (gsl_rng_uniform(gsl_r)<0.33) {repeat = 1+gsl_rng_uniform_int(gsl_r, 20);}
+			else {repeat=1;}
+			
+			logMHratio=0.0;
+			logMHratio+=ProposeHnRho(cHyp_old, cHyp_new, repeat);		
+			logMHratio+=ProposeGamma (rank_old, rank_new, p_gamma, cHyp_old, cHyp_new, repeat);	
+			logMHratio+=ProposePi(cHyp_old, cHyp_new, repeat);
+			
+//			if (fix_sigma>=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; i<cHyp_new.n_gamma; ++i) {
+				  beta_new.push_back(0);
+				}	
+			}
+			else {
+				gsl_matrix *UtXgamma=gsl_matrix_alloc (ni_test, cHyp_new.n_gamma);
+				gsl_vector *beta=gsl_vector_alloc (cHyp_new.n_gamma);
+				SetXgamma (UtXgamma, UtX, rank_new);
+				logPost_new=CalcPosterior(UtXgamma, Utz, K_eval, UtXb_new, Utu_new, alpha_new, beta, cHyp_new);
+				beta_new.clear();
+				for (size_t i=0; i<beta->size; ++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))<logMHratio) {accept=1; n_accept++;}
+			else {accept=0;}
+
+			if (accept==1) {			
+				logPost_old=logPost_new;
+				rank_old.clear(); beta_old.clear();
+				if (rank_new.size()!=0) {
+					for (size_t i=0; i<rank_new.size(); ++i) {
+						rank_old.push_back(rank_new[i]);
+						beta_old.push_back(beta_new[i]);
+					}
+				}
+				cHyp_old=cHyp_new;
+				gsl_vector_memcpy (alpha_old, alpha_new);
+				gsl_vector_memcpy (UtXb_old, UtXb_new);
+				gsl_vector_memcpy (Utu_old, Utu_new);
+			}
+			else {cHyp_new=cHyp_old;}
+		}				
+		
+		//calculate 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);
+			}
+			
+			//sample mu and update z hat
+			gsl_vector_sub (z, z_hat);
+			mean_z+=CenterVector(z);
+			mean_z+=gsl_ran_gaussian(gsl_r, sqrt(1.0/(double) ni_test) );			
+			
+			gsl_vector_add_constant (z_hat, mean_z);
+			
+			time_UtZ+=(clock()-time_start)/(double(CLOCKS_PER_SEC)*60.0);
+		}
+		
+		//Save data
+		if (t<w_step) {continue;}
+		else {		
+			if (t%r_pace==0) {
+				w_col=w%w_pace;
+				if (w_col==0) {
+					if (w==0) {WriteResult (0, Result_hyp, Result_gamma, w_col);}					
+					else {
+						WriteResult (1, Result_hyp, Result_gamma, w_col);
+						gsl_matrix_set_zero (Result_hyp);
+						gsl_matrix_set_zero (Result_gamma);
+					}
+				}
+				
+				gsl_matrix_set (Result_hyp, w_col, 0, cHyp_old.h);
+				gsl_matrix_set (Result_hyp, w_col, 1, cHyp_old.pve);
+				gsl_matrix_set (Result_hyp, w_col, 2, cHyp_old.rho);
+				gsl_matrix_set (Result_hyp, w_col, 3, cHyp_old.pge);
+				gsl_matrix_set (Result_hyp, w_col, 4, cHyp_old.logp);
+				gsl_matrix_set (Result_hyp, w_col, 5, cHyp_old.n_gamma);
+				
+				for (size_t i=0; i<cHyp_old.n_gamma; ++i) {
+					pos=mapRank2pos[rank_old[i]]+1;
+
+					gsl_matrix_set (Result_gamma, w_col, i, pos);
+					
+					beta_g[pos-1].first+=beta_old[i];
+					beta_g[pos-1].second+=1.0;	
+				}
+				
+				gsl_vector_add (alpha_prime, alpha_old);
+				gsl_vector_add (Utu, Utu_old);
+				
+				if (a_mode==13) {
+					pheno_mean+=mean_z;
+				}
+				
+				w++;
+				
+			}
+			
+		}
+	}
+	cout<<endl;
+	
+	w_col=w%w_pace;
+	WriteResult (1, Result_hyp, Result_gamma, w_col);	
+	
+	gsl_matrix_free(Result_hyp);
+	gsl_matrix_free(Result_gamma);	
+	
+	gsl_vector_free(z_hat);
+	gsl_vector_free(z);
+	gsl_vector_free(Utz);	
+	gsl_vector_free(UtXb_new);	
+	gsl_vector_free(UtXb_old);
+	gsl_vector_free(alpha_new);	
+	gsl_vector_free(alpha_old);
+	gsl_vector_free(Utu_new);	
+	gsl_vector_free(Utu_old);	
+	
+	gsl_vector_scale (alpha_prime, 1.0/(double)w);	
+	gsl_vector_scale (Utu, 1.0/(double)w);	
+	if (a_mode==13) {
+		pheno_mean/=(double)w;
+	}
+	
+	gsl_vector *alpha=gsl_vector_alloc (ns_test);
+	gsl_blas_dgemv (CblasTrans, 1.0/(double)ns_test, UtX, alpha_prime, 0.0, alpha);	
+	WriteParam (beta_g, alpha, w);
+	gsl_vector_free(alpha);
+	
+	gsl_blas_dgemv (CblasNoTrans, 1.0, U, Utu, 0.0, alpha_prime);
+	WriteBV(alpha_prime);	
+	
+	gsl_vector_free(alpha_prime);				
+	gsl_vector_free(Utu);	
+		
+	delete [] p_gamma;
+	beta_g.clear();
+	
+	return;
+}
+
+
+
+void BSLMM::RidgeR(const gsl_matrix *U, const gsl_matrix *UtX, const gsl_vector *Uty, const gsl_vector *eval, const double lambda)
+{
+	gsl_vector *beta=gsl_vector_alloc (UtX->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<size_t> &rank_old, const vector<size_t> &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<size_t> rank_remove(v_size), rank_add(v_size), rank_union(s_max+v_size);
+  vector<size_t>::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<size_t, int> mapRank2in_remove, mapRank2in_add;
+  for (size_t i=0; i<rank_remove.size(); i++) {
+    mapRank2in_remove[rank_remove[i]]=1;
+  }
+  for (size_t i=0; i<rank_add.size(); i++) {
+    mapRank2in_add[rank_add[i]]=1;
+  }
+
+  //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);
+    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; i<rank_old.size(); i++) {
+	if (mapRank2in_remove.count(rank_old[i])!=0) {continue;}
+	gsl_vector_const_view Xold_col=gsl_matrix_const_column(X_old, i);	
+	gsl_vector_view Xtemp_col=gsl_matrix_column(X_temp, i_temp);
+	gsl_vector_memcpy (&Xtemp_col.vector, &Xold_col.vector);
+
+	d=gsl_vector_get (Xty_old, i);
+	gsl_vector_set (Xty_temp, i_temp, d);
+	
+	j_temp=i_temp;
+	for (size_t j=i; j<rank_old.size(); j++) {
+	  if (mapRank2in_remove.count(rank_old[j])!=0) {continue;}
+	  d=gsl_matrix_get (XtX_old, i, j);
+	  gsl_matrix_set (XtX_temp, i_temp, j_temp, d);
+	  if (i_temp!=j_temp) {gsl_matrix_set (XtX_temp, j_temp, i_temp, d);}
+	  j_temp++;
+	}
+	i_temp++;
+      }
+    }
+
+    if (rank_add.size()==0) {
+      gsl_matrix_memcpy (&Xnew_sub.matrix, X_temp);
+      gsl_matrix_memcpy (&XtXnew_sub.matrix, XtX_temp);
+      gsl_vector_memcpy (&Xtynew_sub.vector, Xty_temp);
+    } else {
+      gsl_matrix *X_add=gsl_matrix_alloc(X_old->size1, 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; i<rank_new.size(); i++) {
+	if (mapRank2in_add.count(rank_new[i])!=0) {i_flag=1;} else {i_flag=0;}
+	gsl_vector_view Xnew_col=gsl_matrix_column(X_new, i); 
+	if (i_flag==1) {
+	  gsl_vector_view Xcopy_col=gsl_matrix_column(X_add, i-i_temp);
+	  gsl_vector_memcpy (&Xnew_col.vector, &Xcopy_col.vector);
+	} else {
+	  gsl_vector_view Xcopy_col=gsl_matrix_column(X_temp, i_temp);	  
+	  gsl_vector_memcpy (&Xnew_col.vector, &Xcopy_col.vector);
+	}	
+
+	if (i_flag==1) {
+          d=gsl_vector_get (Xty_add, i-i_temp);
+        } else {
+          d=gsl_vector_get (Xty_temp, i_temp);
+        }
+	gsl_vector_set (Xty_new, i, d);
+
+	j_temp=i_temp;
+	for (size_t j=i; j<rank_new.size(); j++) {
+          if (mapRank2in_add.count(rank_new[j])!=0) {j_flag=1;} else {j_flag=0;}
+
+	  if (i_flag==1 && j_flag==1) {
+            d=gsl_matrix_get(XtX_aa, i-i_temp, j-j_temp);
+	  } else if (i_flag==1) {
+	    d=gsl_matrix_get(XtX_at, i-i_temp, j_temp);
+	  } else if (j_flag==1) {
+	    d=gsl_matrix_get(XtX_at, j-j_temp, i_temp);
+	  } else {
+	    d=gsl_matrix_get(XtX_temp, i_temp, j_temp);
+	  }
+
+	  gsl_matrix_set (XtX_new, i, j, d);
+	  if (i!=j) {gsl_matrix_set (XtX_new, j, i, d);}
+
+	  if (j_flag==0) {j_temp++;}
+        }
+	if (i_flag==0) {i_temp++;}
+      }
+
+      gsl_matrix_free(X_add);
+      gsl_matrix_free(XtX_aa);
+      gsl_matrix_free(XtX_at);
+      gsl_vector_free(Xty_add);
+    }
+
+    gsl_matrix_free(X_temp);
+    gsl_matrix_free(XtX_temp);
+    gsl_vector_free(Xty_temp);
+  }
+  */
+
+
+  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 {
+    size_t i_old, j_old, i_new, j_new, i_add, j_add, i_flag, j_flag;
+    if (rank_add.size()==0) {
+      i_old=0; i_new=0;
+      for (size_t i=0; i<rank_union.size(); i++) {
+	if (mapRank2in_remove.count(rank_old[i_old])!=0) {i_old++; continue;}
+
+	gsl_vector_view Xnew_col=gsl_matrix_column(X_new, i_new); 
+	gsl_vector_const_view Xcopy_col=gsl_matrix_const_column(X_old, i_old);
+	gsl_vector_memcpy (&Xnew_col.vector, &Xcopy_col.vector);
+
+	d=gsl_vector_get (Xty_old, i_old);
+	gsl_vector_set (Xty_new, i_new, d);
+
+	j_old=i_old; j_new=i_new;
+	for (size_t j=i; j<rank_union.size(); j++) {
+          if (mapRank2in_remove.count(rank_old[j_old])!=0) {j_old++; continue;}
+
+	  d=gsl_matrix_get(XtX_old, i_old, j_old);
+
+	  gsl_matrix_set (XtX_new, i_new, j_new, d);
+	  if (i_new!=j_new) {gsl_matrix_set (XtX_new, j_new, i_new, d);}
+
+	  j_old++; j_new++;
+        }
+	i_old++; i_new++;
+      }	
+    } else {
+      gsl_matrix *X_add=gsl_matrix_alloc(X_old->size1, 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_old, 0.0, XtX_ao);
+      
+      //#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; i<rank_union.size(); i++) {
+	if (mapRank2in_remove.count(rank_old[i_old])!=0) {i_old++; continue;}
+	if (mapRank2in_add.count(rank_new[i_new])!=0) {i_flag=1;} else {i_flag=0;}
+
+	gsl_vector_view Xnew_col=gsl_matrix_column(X_new, i_new); 
+	if (i_flag==1) {
+	  gsl_vector_view Xcopy_col=gsl_matrix_column(X_add, i_add);
+	  gsl_vector_memcpy (&Xnew_col.vector, &Xcopy_col.vector);
+	} else {
+	  gsl_vector_const_view Xcopy_col=gsl_matrix_const_column(X_old, i_old);	  
+	  gsl_vector_memcpy (&Xnew_col.vector, &Xcopy_col.vector);
+	}	
+
+	if (i_flag==1) {
+          d=gsl_vector_get (Xty_add, i_add);
+        } else {
+          d=gsl_vector_get (Xty_old, i_old);
+        }
+	gsl_vector_set (Xty_new, i_new, d);
+
+	j_old=i_old; j_new=i_new; j_add=i_add;
+	for (size_t j=i; j<rank_union.size(); j++) {
+	  if (mapRank2in_remove.count(rank_old[j_old])!=0) {j_old++; continue;}
+	  if (mapRank2in_add.count(rank_new[j_new])!=0) {j_flag=1;} else {j_flag=0;}
+
+	  if (i_flag==1 && j_flag==1) {
+            d=gsl_matrix_get(XtX_aa, i_add, j_add);
+	  } else if (i_flag==1) {
+	    d=gsl_matrix_get(XtX_ao, i_add, j_old);
+	  } else if (j_flag==1) {
+	    d=gsl_matrix_get(XtX_ao, j_add, i_old);
+	  } else {
+	    d=gsl_matrix_get(XtX_old, i_old, j_old);
+	  }
+
+	  gsl_matrix_set (XtX_new, i_new, j_new, d);
+	  if (i_new!=j_new) {gsl_matrix_set (XtX_new, j_new, i_new, d);}
+
+	  j_new++; if (j_flag==1) {j_add++;} else {j_old++;}
+        }
+	i_new++; if (i_flag==1) {i_add++;} else {i_old++;}
+      }
+
+      gsl_matrix_free(X_add);
+      gsl_matrix_free(XtX_aa);
+      gsl_matrix_free(XtX_ao);
+      gsl_vector_free(Xty_add);
+    }
+
+  }
+
+  rank_remove.clear();
+  rank_add.clear();
+  rank_union.clear();
+  mapRank2in_remove.clear();
+  mapRank2in_add.clear();
+	
+  return;
+}
+
+
+double BSLMM::CalcPosterior (const double yty, class HYPBSLMM &cHyp)
+{	
+	double logpost=0.0;
+	
+	//for quantitative traits, calculate pve and pge
+	//pve and pge for case/control data are calculted in CalcCC_PVEnZ
+	if (a_mode==11) {
+		cHyp.pve=0.0;
+		cHyp.pge=1.0;	
+	}
+
+	//calculate likelihood
+	if (a_mode==11) {logpost-=0.5*(double)ni_test*log(yty);}
+	else {logpost-=0.5*yty;}
+	
+	logpost+=((double)cHyp.n_gamma-1.0)*cHyp.logp+((double)ns_test-(double)cHyp.n_gamma)*log(1-exp(cHyp.logp));
+		
+	return logpost;
+}
+
+
+double BSLMM::CalcPosterior (const gsl_matrix *Xgamma, const gsl_matrix *XtX, const gsl_vector *Xty, const double yty, const size_t s_size, gsl_vector *Xb, gsl_vector *beta, class HYPBSLMM &cHyp)
+{	
+	double sigma_a2=cHyp.h/( (1-cHyp.h)*exp(cHyp.logp)*(double)ns_test);
+	double logpost=0.0;
+	double d, P_yy=yty, logdet_O=0.0;
+
+	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);
+	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<s_size; i++)
+	{
+		d=gsl_ran_gaussian(gsl_r, 1); 
+		gsl_vector_set(beta, i, d); 
+	}
+	gsl_vector_view beta_sub=gsl_vector_subvector(beta, 0, s_size);
+	gsl_blas_dtrsv(CblasUpper, CblasNoTrans, CblasNonUnit, Omega, &beta_sub.vector); 
+		
+	//it compuates inv(L^T(Omega)) %*% beta;  
+	gsl_vector_scale(&beta_sub.vector, sqrt(sigma_a2/tau));
+	gsl_vector_add(&beta_sub.vector, beta_hat); 
+	gsl_blas_dgemv (CblasNoTrans, 1.0, &Xgamma_sub.matrix, &beta_sub.vector, 0.0, Xb);		
+	
+	//for quantitative traits, calculate pve and pge
+	if (a_mode==11) {
+		gsl_blas_ddot (Xb, Xb, &d);
+		cHyp.pve=d/(double)ni_test;
+		cHyp.pve/=cHyp.pve+1.0/tau;
+		cHyp.pge=1.0;	
+	}	
+	
+	logpost=-0.5*logdet_O;
+	if (a_mode==11) {logpost-=0.5*(double)ni_test*log(P_yy);}
+	else {logpost-=0.5*P_yy;}
+
+	logpost+=((double)cHyp.n_gamma-1.0)*cHyp.logp+((double)ns_test-(double)cHyp.n_gamma)*log(1.0-exp(cHyp.logp));
+
+	gsl_matrix_free (Omega);
+	gsl_matrix_free (M_temp);
+	gsl_vector_free (beta_hat);
+	gsl_vector_free (Xty_temp);
+
+	return logpost;
+}
+
+
+
+//calculate pve and pge, and calculate z_hat for case-control data	
+void BSLMM::CalcCC_PVEnZ (gsl_vector *z_hat, class HYPBSLMM &cHyp) 
+{
+  gsl_vector_set_zero(z_hat);
+  cHyp.pve=0.0;
+  cHyp.pge=1.0;		
+  return;
+}
+
+
+//calculate pve and pge, and calculate z_hat for case-control data	
+void BSLMM::CalcCC_PVEnZ (const gsl_vector *Xb, gsl_vector *z_hat, class HYPBSLMM &cHyp) 
+{
+	double d;
+	
+	gsl_blas_ddot (Xb, Xb, &d);
+	cHyp.pve=d/(double)ni_test;
+	cHyp.pve/=cHyp.pve+1.0;
+	cHyp.pge=1.0;
+	
+	gsl_vector_memcpy (z_hat, Xb);
+
+	return;
+}
+
+
+
+//if a_mode==13, then run probit model
+void BSLMM::MCMC (const gsl_matrix *X, const gsl_vector *y) {
+	clock_t time_start;	
+	double time_set=0, time_post=0;
+
+	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 *Xb_new=gsl_vector_alloc (ni_test);
+	gsl_vector *Xb_old=gsl_vector_alloc (ni_test);	
+	gsl_vector *z_hat=gsl_vector_alloc (ni_test);
+	gsl_vector *z=gsl_vector_alloc (ni_test);
+
+	gsl_matrix *Xgamma_old=gsl_matrix_alloc (ni_test, s_max);
+	gsl_matrix *XtX_old=gsl_matrix_alloc (s_max, s_max);
+	gsl_vector *Xtz_old=gsl_vector_alloc (s_max);
+	gsl_vector *beta_old=gsl_vector_alloc (s_max);
+
+	gsl_matrix *Xgamma_new=gsl_matrix_alloc (ni_test, s_max);
+	gsl_matrix *XtX_new=gsl_matrix_alloc (s_max, s_max);
+	gsl_vector *Xtz_new=gsl_vector_alloc (s_max);
+	gsl_vector *beta_new=gsl_vector_alloc (s_max);
+
+	double ztz=0.0;
+	gsl_vector_memcpy (z, y);
+	//for quantitative traits, y is centered already in gemma.cpp, but just in case
+	double mean_z=CenterVector (z);				
+	gsl_blas_ddot(z, z, &ztz);
+
+	double logPost_new, logPost_old;
+	double logMHratio;
+	
+	gsl_matrix_set_zero (Result_gamma);
+	if (a_mode==13) {
+		pheno_mean=0.0;
+	}
+	
+	vector<pair<double, double> > beta_g;
+	for (size_t i=0; i<ns_test; i++) {
+		beta_g.push_back(make_pair(0.0, 0.0));
+	}
+	
+	vector<size_t> rank_new, rank_old;
+	vector<pair<size_t, double> > 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; i<ns_test; ++i) {
+		mapRank2pos[i]=pos_loglr[i].first;
+	}
+	
+	//calculate proposal distribution for gamma (unnormalized), and set up gsl_r and gsl_t		
+	gsl_rng_env_setup();                
+	const gsl_rng_type * gslType;                                               
+	gslType = gsl_rng_default; 
+	if (randseed<0)
+	{
+		time_t rawtime;
+		time (&rawtime);
+		tm * ptm = gmtime (&rawtime);
+		
+		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);
+	
+	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; t<total_step; ++t) {
+		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 (t>10) {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; i<n_mh; ++i) {
+			if (gsl_rng_uniform(gsl_r)<0.33) {repeat = 1+gsl_rng_uniform_int(gsl_r, 20);}
+			else {repeat=1;}
+
+			logMHratio=0.0;
+			logMHratio+=ProposeHnRho(cHyp_old, cHyp_new, repeat);	
+			logMHratio+=ProposeGamma (rank_old, rank_new, p_gamma, cHyp_old, cHyp_new, repeat);	
+			logMHratio+=ProposePi(cHyp_old, cHyp_new, repeat);
+			
+			if (cHyp_new.n_gamma==0) {
+				logPost_new=CalcPosterior (ztz, cHyp_new);
+			} else {
+			  //this if makes sure that rank_old.size()==rank_remove.size() does not happen
+			  if (cHyp_new.n_gamma<=20 || cHyp_old.n_gamma<=20) {
+			    time_start=clock();
+			    SetXgamma (Xgamma_new, X, rank_new);	  
+			    CalcXtX (Xgamma_new, z, rank_new.size(), XtX_new, Xtz_new);	
+			    time_set+=(clock()-time_start)/(double(CLOCKS_PER_SEC)*60.0);
+			  } else {
+			    time_start=clock();
+			    SetXgamma (X, Xgamma_old, XtX_old, Xtz_old, z, rank_old, rank_new, Xgamma_new, XtX_new, Xtz_new);
+			    time_set+=(clock()-time_start)/(double(CLOCKS_PER_SEC)*60.0);
+			  }
+			  time_start=clock();
+			  logPost_new=CalcPosterior (Xgamma_new, XtX_new, Xtz_new, ztz, rank_new.size(), Xb_new, beta_new, cHyp_new);
+			  time_post+=(clock()-time_start)/(double(CLOCKS_PER_SEC)*60.0);
+			}	
+			logMHratio+=logPost_new-logPost_old;	
+		
+			if (logMHratio>0 || log(gsl_rng_uniform(gsl_r))<logMHratio) {accept=1; n_accept++;}
+			else {accept=0;}
+			
+			//cout<<rank_new.size()<<"\t"<<rank_old.size()<<"\t"<<logPost_new<<"\t"<<logPost_old<<endl;
+			if (accept==1) {			
+				logPost_old=logPost_new;
+				cHyp_old=cHyp_new;
+				gsl_vector_memcpy (Xb_old, Xb_new);
+
+				rank_old.clear();
+				if (rank_new.size()!=0) {
+					for (size_t i=0; i<rank_new.size(); ++i) {
+						rank_old.push_back(rank_new[i]);
+					}
+								
+					gsl_matrix_view Xold_sub=gsl_matrix_submatrix(Xgamma_old, 0, 0, ni_test, rank_new.size());
+					gsl_matrix_view XtXold_sub=gsl_matrix_submatrix(XtX_old, 0, 0, rank_new.size(), rank_new.size());
+					gsl_vector_view Xtzold_sub=gsl_vector_subvector(Xtz_old, 0, rank_new.size());
+					gsl_vector_view betaold_sub=gsl_vector_subvector(beta_old, 0, rank_new.size());
+
+					gsl_matrix_view Xnew_sub=gsl_matrix_submatrix(Xgamma_new, 0, 0, ni_test, rank_new.size());
+					gsl_matrix_view XtXnew_sub=gsl_matrix_submatrix(XtX_new, 0, 0, rank_new.size(), rank_new.size());
+					gsl_vector_view Xtznew_sub=gsl_vector_subvector(Xtz_new, 0, rank_new.size());
+					gsl_vector_view betanew_sub=gsl_vector_subvector(beta_new, 0, rank_new.size());
+
+					gsl_matrix_memcpy(&Xold_sub.matrix, &Xnew_sub.matrix);
+					gsl_matrix_memcpy(&XtXold_sub.matrix, &XtXnew_sub.matrix);
+					gsl_vector_memcpy(&Xtzold_sub.vector, &Xtznew_sub.vector);
+					gsl_vector_memcpy(&betaold_sub.vector, &betanew_sub.vector);
+				}
+			} else {
+			  cHyp_new=cHyp_old;
+			}
+			
+		}				
+
+		//calculate 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);
+			}
+			
+			//sample mu and update z hat
+			gsl_vector_sub (z, z_hat);
+			mean_z+=CenterVector(z);
+			mean_z+=gsl_ran_gaussian(gsl_r, sqrt(1.0/(double) ni_test) );			
+			
+			gsl_vector_add_constant (z_hat, mean_z);
+		}
+		
+		//Save data
+		if (t<w_step) {continue;}
+		else {		
+			if (t%r_pace==0) {
+				w_col=w%w_pace;
+				if (w_col==0) {
+					if (w==0) {WriteResult (0, Result_hyp, Result_gamma, w_col);}					
+					else {
+						WriteResult (1, Result_hyp, Result_gamma, w_col);
+						gsl_matrix_set_zero (Result_hyp);
+						gsl_matrix_set_zero (Result_gamma);
+					}
+				}
+
+				gsl_matrix_set (Result_hyp, w_col, 0, cHyp_old.h);
+				gsl_matrix_set (Result_hyp, w_col, 1, cHyp_old.pve);
+				gsl_matrix_set (Result_hyp, w_col, 2, cHyp_old.rho);
+				gsl_matrix_set (Result_hyp, w_col, 3, cHyp_old.pge);
+				gsl_matrix_set (Result_hyp, w_col, 4, cHyp_old.logp);
+				gsl_matrix_set (Result_hyp, w_col, 5, cHyp_old.n_gamma);
+				
+				for (size_t i=0; i<cHyp_old.n_gamma; ++i) {
+					pos=mapRank2pos[rank_old[i]]+1;
+
+					gsl_matrix_set (Result_gamma, w_col, i, pos);
+					
+					beta_g[pos-1].first+=gsl_vector_get(beta_old, i);
+					beta_g[pos-1].second+=1.0;	
+				}
+				
+				if (a_mode==13) {
+					pheno_mean+=mean_z;
+				}
+				
+				w++;
+				
+			}
+			
+		}
+	}
+	cout<<endl;
+
+	cout<<"time on selecting Xgamma: "<<time_set<<endl;
+	cout<<"time on calculating posterior: "<<time_post<<endl;
+
+	w_col=w%w_pace;
+	WriteResult (1, Result_hyp, Result_gamma, w_col);	
+	
+	gsl_vector *alpha=gsl_vector_alloc (ns_test);
+	gsl_vector_set_zero (alpha);
+	WriteParam (beta_g, alpha, w);
+	gsl_vector_free(alpha);
+
+	gsl_matrix_free(Result_hyp);
+	gsl_matrix_free(Result_gamma);	
+	
+	gsl_vector_free(z_hat);
+	gsl_vector_free(z);
+	gsl_vector_free(Xb_new);	
+	gsl_vector_free(Xb_old);
+
+	gsl_matrix_free(Xgamma_old);
+	gsl_matrix_free(XtX_old);
+	gsl_vector_free(Xtz_old);
+	gsl_vector_free(beta_old);
+
+	gsl_matrix_free(Xgamma_new);
+	gsl_matrix_free(XtX_new);
+	gsl_vector_free(Xtz_new);
+	gsl_vector_free(beta_new);
+	
+	delete [] p_gamma;
+	beta_g.clear();
+	
+	return;
+}