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 +++++++++++++++++++++++++++++--------------------------
src/bslmm.h | 158 +++++++----
src/bslmmdap.cpp | 260 +++++++++--------
src/bslmmdap.h | 62 ++--
src/gemma.h | 17 +-
src/lmm.h | 113 +++++---
src/logistic.cpp | 46 +--
src/mvlmm.h | 81 +++---
8 files changed, 874 insertions(+), 702 deletions(-)
(limited to 'src')
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)) indicator_idv; //indicator for individuals (phenotypes), 0 missing, 1 available for analysis
- vector indicator_snp; //sequence indicator for SNPs: 0 ignored because of (a) maf, (b) miss, (c) non-poly; 1 available for analysis
-
- vector snpInfo; //record SNP information
+ double time_Omega; // Time spent on optimization iterations.
+
+ // Time spent on constructing the proposal distribution for
+ // gamma (i.e. lmm or lm analysis).
+ double time_Proposal;
+
+ // Indicator for individuals (phenotypes): 0 missing, 1
+ // available for analysis.
+ vector indicator_idv;
+
+ // Sequence indicator for SNPs: 0 ignored because of (a) maf,
+ // (b) miss, (c) non-poly; 1 available for analysis.
+ vector indicator_snp;
+
+ // Record SNP information.
+ vector snpInfo;
- // Not included in PARAM
+ // Not included in PARAM.
gsl_rng *gsl_r;
gsl_ran_discrete_t *gsl_t;
map mapRank2pos;
- // Main Functions
+ // Main functions.
void CopyFromParam (PARAM &cPar);
void CopyToParam (PARAM &cPar);
- void RidgeR(const gsl_matrix *U, const gsl_matrix *UtX, const gsl_vector *Uty, const gsl_vector *eval, const double lambda);
+ void RidgeR(const gsl_matrix *U, const gsl_matrix *UtX,
+ const gsl_vector *Uty, const gsl_vector *eval,
+ const double lambda);
- void MCMC (const gsl_matrix *U, const gsl_matrix *UtX, const gsl_vector *Uty, const gsl_vector *K_eval, const gsl_vector *y);
+ void MCMC (const gsl_matrix *U, const gsl_matrix *UtX,
+ const gsl_vector *Uty, const gsl_vector *K_eval,
+ const gsl_vector *y);
void WriteLog ();
void WriteLR ();
void WriteBV (const gsl_vector *bv);
- void WriteParam (vector > &beta_g, const gsl_vector *alpha, const size_t w);
+ void WriteParam (vector > &beta_g,
+ const gsl_vector *alpha, const size_t w);
void WriteParam (const gsl_vector *alpha);
- void WriteResult (const int flag, const gsl_matrix *Result_hyp, const gsl_matrix *Result_gamma, const size_t w_col);
+ void WriteResult (const int flag, const gsl_matrix *Result_hyp,
+ const gsl_matrix *Result_gamma, const size_t w_col);
- //Subfunctions inside MCMC
+ // Subfunctions inside MCMC.
void CalcPgamma (double *p_gammar);
- double CalcPveLM (const gsl_matrix *UtXgamma, const gsl_vector *Uty, const double sigma_a2);
- void InitialMCMC (const gsl_matrix *UtX, const gsl_vector *Uty, vector &rank_old, class HYPBSLMM &cHyp, vector > &pos_loglr);
- double CalcPosterior (const gsl_vector *Uty, const gsl_vector *K_eval, gsl_vector *Utu, gsl_vector *alpha_prime, class HYPBSLMM &cHyp);
- double 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);
- void CalcCC_PVEnZ (const gsl_matrix *U, const gsl_vector *Utu, gsl_vector *z_hat, class HYPBSLMM &cHyp);
- void CalcCC_PVEnZ (const gsl_matrix *U, const gsl_vector *UtXb, const gsl_vector *Utu, gsl_vector *z_hat, class HYPBSLMM &cHyp);
- double CalcREMLE (const gsl_matrix *Utw, const gsl_vector *Uty, const gsl_vector *K_eval);
- double CalcLR (const gsl_matrix *U, const gsl_matrix *UtX, const gsl_vector *Uty, const gsl_vector *K_eval, vector > &loglr_sort); //calculate the maximum marginal likelihood ratio for each analyzed SNPs with gemma, use it to rank SNPs
- void SampleZ (const gsl_vector *y, const gsl_vector *z_hat, gsl_vector *z);
- double ProposeHnRho (const class HYPBSLMM &cHyp_old, class HYPBSLMM &cHyp_new, const size_t &repeat);
- double ProposePi (const class HYPBSLMM &cHyp_old, class HYPBSLMM &cHyp_new, const size_t &repeat);
- double 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);
- void SetXgamma (gsl_matrix *Xgamma, const gsl_matrix *X, vector &rank);
-
- void CalcXtX (const gsl_matrix *X_new, const gsl_vector *y, const size_t s_size, gsl_matrix *XtX_new, gsl_vector *Xty_new);
- void 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 CalcPveLM (const gsl_matrix *UtXgamma, const gsl_vector *Uty,
+ const double sigma_a2);
+ void InitialMCMC (const gsl_matrix *UtX, const gsl_vector *Uty,
+ vector &rank_old, class HYPBSLMM &cHyp,
+ vector > &pos_loglr);
+ double CalcPosterior (const gsl_vector *Uty, const gsl_vector *K_eval,
+ gsl_vector *Utu, gsl_vector *alpha_prime,
+ class HYPBSLMM &cHyp);
+ double 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);
+ void CalcCC_PVEnZ (const gsl_matrix *U, const gsl_vector *Utu,
+ gsl_vector *z_hat, class HYPBSLMM &cHyp);
+ void CalcCC_PVEnZ (const gsl_matrix *U, const gsl_vector *UtXb,
+ const gsl_vector *Utu, gsl_vector *z_hat,
+ class HYPBSLMM &cHyp);
+ double CalcREMLE (const gsl_matrix *Utw, const gsl_vector *Uty,
+ const gsl_vector *K_eval);
+
+ // Calculate the maximum marginal likelihood ratio for each
+ // analyzed SNPs with gemma, use it to rank SNPs.
+ double CalcLR (const gsl_matrix *U, const gsl_matrix *UtX,
+ const gsl_vector *Uty, const gsl_vector *K_eval,
+ vector > &loglr_sort);
+ void SampleZ (const gsl_vector *y, const gsl_vector *z_hat,
+ gsl_vector *z);
+ double ProposeHnRho (const class HYPBSLMM &cHyp_old,
+ class HYPBSLMM &cHyp_new, const size_t &repeat);
+ double ProposePi (const class HYPBSLMM &cHyp_old,
+ class HYPBSLMM &cHyp_new,
+ const size_t &repeat);
+ double 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);
+ void SetXgamma (gsl_matrix *Xgamma, const gsl_matrix *X,
+ vector &rank);
+
+ void CalcXtX (const gsl_matrix *X_new, const gsl_vector *y,
+ const size_t s_size, gsl_matrix *XtX_new,
+ gsl_vector *Xty_new);
+ void 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 CalcPosterior (const double yty, class HYPBSLMM &cHyp);
- double 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 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);
void CalcCC_PVEnZ (gsl_vector *z_hat, class HYPBSLMM &cHyp);
- void CalcCC_PVEnZ (const gsl_vector *Xb, gsl_vector *z_hat, class HYPBSLMM &cHyp);
+ void CalcCC_PVEnZ (const gsl_vector *Xb, gsl_vector *z_hat,
+ class HYPBSLMM &cHyp);
void MCMC (const gsl_matrix *X, const gsl_vector *y);
-
- //utility functions
-// double vec_sum (gsl_vector *v);
-// void vec_center (gsl_vector *v);
-// double calc_var (gsl_vector *v);
-// void calc_sigma (MCMC &cMcmc);
-// bool comp_lr (pair a, pair b);
};
-
-
#endif
diff --git a/src/bslmmdap.cpp b/src/bslmmdap.cpp
index e1b20d8..ebbff70 100644
--- a/src/bslmmdap.cpp
+++ b/src/bslmmdap.cpp
@@ -49,8 +49,7 @@
using namespace std;
-void BSLMMDAP::CopyFromParam (PARAM &cPar)
-{
+void BSLMMDAP::CopyFromParam (PARAM &cPar) {
file_out=cPar.file_out;
path_out=cPar.path_out;
@@ -83,9 +82,7 @@ void BSLMMDAP::CopyFromParam (PARAM &cPar)
return;
}
-
-void BSLMMDAP::CopyToParam (PARAM &cPar)
-{
+void BSLMMDAP::CopyToParam (PARAM &cPar) {
cPar.time_UtZ=time_UtZ;
cPar.time_Omega=time_Omega;
@@ -94,13 +91,16 @@ void BSLMMDAP::CopyToParam (PARAM &cPar)
-//read hyp file
-void ReadFile_hyb (const string &file_hyp, vector &vec_sa2, vector &vec_sb2, vector &vec_wab)
-{
+// Read hyp file.
+void ReadFile_hyb (const string &file_hyp, vector &vec_sa2,
+ vector &vec_sb2, vector &vec_wab) {
vec_sa2.clear(); vec_sb2.clear(); vec_wab.clear();
igzstream infile (file_hyp.c_str(), igzstream::in);
- if (!infile) {cout<<"error! fail to open hyp file: "< &vec_sa2, vector &vec_rs, vector > > &BF)
-{
+// Read bf file.
+void ReadFile_bf (const string &file_bf, vector &vec_rs,
+ vector > > &BF) {
BF.clear(); vec_rs.clear();
igzstream infile (file_bf.c_str(), igzstream::in);
@@ -172,7 +171,9 @@ void ReadFile_bf (const string &file_bf, vector &vec_rs, vector &vec_rs, vector &vec_rs, gsl_matrix *Ac, gsl_matrix_int *Ad, gsl_vector_int *dlevel, size_t &kc, size_t &kd)
-{
+// Read category files.
+// Read both continuous and discrete category file, record mapRS2catc.
+void ReadFile_cat (const string &file_cat, const vector &vec_rs,
+ gsl_matrix *Ac, gsl_matrix_int *Ad, gsl_vector_int *dlevel,
+ size_t &kc, size_t &kd) {
igzstream infile (file_cat.c_str(), igzstream::in);
- if (!infile) {cout<<"error! fail to open category file: "< &vec_rs, gsl_mat
vector catc;
vector catd;
- //read the following lines to record mapRS2cat
+ // Read the following lines to record mapRS2cat.
while (!safeGetline(infile, line).eof()) {
ch_ptr=strtok ((char *)line.c_str(), " , \t");
@@ -255,7 +260,7 @@ void ReadFile_cat (const string &file_cat, const vector &vec_rs, gsl_mat
if (mapRS2catd.count(rs)==0 && kd>0) {mapRS2catd[rs]=catd;}
}
- //load into Ad and Ac
+ // Load into Ad and Ac.
if (kc>0) {
Ac=gsl_matrix_alloc(vec_rs.size(), kc);
for (size_t i=0; i &vec_rs, gsl_mat
return;
}
-
-
-
-
-
-
-
-void BSLMMDAP::WriteResult (const gsl_matrix *Hyper, const gsl_matrix *BF)
-{
+void BSLMMDAP::WriteResult (const gsl_matrix *Hyper, const gsl_matrix *BF) {
string file_bf, file_hyp;
file_bf=path_out+"/"+file_out;
file_bf+=".bf.txt";
@@ -325,10 +322,17 @@ void BSLMMDAP::WriteResult (const gsl_matrix *Hyper, const gsl_matrix *BF)
outfile_bf.open (file_bf.c_str(), ofstream::out);
outfile_hyp.open (file_hyp.c_str(), ofstream::out);
- if (!outfile_bf) {cout<<"error writing file: "<size1; i++) {
for (size_t j=0; jsize2; j++) {
@@ -366,10 +370,9 @@ void BSLMMDAP::WriteResult (const gsl_matrix *Hyper, const gsl_matrix *BF)
return;
}
-
-
-void BSLMMDAP::WriteResult (const vector &vec_rs, const gsl_matrix *Hyper, const gsl_vector *pip, const gsl_vector *coef)
-{
+void BSLMMDAP::WriteResult (const vector &vec_rs,
+ const gsl_matrix *Hyper, const gsl_vector *pip,
+ const gsl_vector *coef) {
string file_gamma, file_hyp, file_coef;
file_gamma=path_out+"/"+file_out;
file_gamma+=".gamma.txt";
@@ -384,11 +387,21 @@ void BSLMMDAP::WriteResult (const vector &vec_rs, const gsl_matrix *Hype
outfile_hyp.open (file_hyp.c_str(), ofstream::out);
outfile_coef.open (file_coef.c_str(), ofstream::out);
- if (!outfile_gamma) {cout<<"error writing file: "<size1; i++) {
for (size_t j=0; jsize2; j++) {
@@ -397,10 +410,10 @@ void BSLMMDAP::WriteResult (const vector &vec_rs, const gsl_matrix *Hype
outfile_hyp< &vec_rs, const gsl_matrix *Hype
}
-
-
-/*
-void BSLMMDAP::SetXgamma (gsl_matrix *Xgamma, const gsl_matrix *X, vector &rank)
-{
- size_t pos;
- for (size_t i=0; isize);
double logm=0.0;
@@ -452,7 +449,7 @@ double BSLMMDAP::CalcMarginal (const gsl_vector *Uty, const gsl_vector *K_eval,
Hi_yy+=d*uy*uy;
}
- //calculate likelihood
+ // Calculate likelihood.
logm=-0.5*logdet_H-0.5*tau*Hi_yy+0.5*log(tau)*(double)ni_test;
gsl_vector_free (weight_Hi);
@@ -460,14 +457,17 @@ double BSLMMDAP::CalcMarginal (const gsl_vector *Uty, const gsl_vector *K_eval,
return logm;
}
-
-double BSLMMDAP::CalcMarginal (const gsl_matrix *UtXgamma, const gsl_vector *Uty, const gsl_vector *K_eval, const double sigma_a2, const double sigma_b2, const double tau)
-{
+double BSLMMDAP::CalcMarginal (const gsl_matrix *UtXgamma,
+ const gsl_vector *Uty,
+ const gsl_vector *K_eval,
+ const double sigma_a2,
+ const double sigma_b2, const double tau) {
clock_t time_start;
double logm=0.0;
double d, 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 *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);
@@ -477,7 +477,7 @@ double BSLMMDAP::CalcMarginal (const gsl_matrix *UtXgamma, const gsl_vector *Uty
logdet_H=0.0; P_yy=0.0;
for (size_t i=0; i sum_pip;
map sum;
@@ -623,27 +626,26 @@ void single_ct_regression(const gsl_matrix_int *Xd, const gsl_vector_int *dlevel
gsl_vector_set(prior_vec,i,sum_pip[cat]/sum[cat]);
}
- //double baseline=0;
for(int i=0;i &vec_rs, const vector &vec_sa2, const vector &vec_sb2, const vector &wab, const vector > > &BF, gsl_matrix *Ac, gsl_matrix_int *Ad, gsl_vector_int *dlevel) {
+// Where A is the ni_test by n_cat matrix of annotations.
+void BSLMMDAP::DAP_EstimateHyper (const size_t kc, const size_t kd,
+ const vector &vec_rs,
+ const vector &vec_sa2,
+ const vector &vec_sb2,
+ const vector &wab,
+ const vector > > &BF,
+ gsl_matrix *Ac, gsl_matrix_int *Ad,
+ gsl_vector_int *dlevel) {
clock_t time_start;
- //set up BF
+ // Set up BF.
double h, rho, sigma_a2, sigma_b2, d, s, logm, logm_save;
size_t t1, t2;
size_t n_grid=wab.size(), ns_test=vec_rs.size();
@@ -653,10 +655,10 @@ void BSLMMDAP::DAP_EstimateHyper (const size_t kc, const size_t kd, const vector
gsl_vector *pip=gsl_vector_alloc(ns_test);
gsl_vector *coef=gsl_vector_alloc(kc+kd+1);
- //perform the EM algorithm
+ // Perform the EM algorithm.
vector vec_wab, vec_wab_new;
- //initial values
+ // Initial values.
for (size_t t=0; t1e-3) {
- //update E_gamma
+
+ // Update E_gamma.
t1=0, t2=0;
for (size_t b=0; bsize; t++) {
s+=gsl_vector_get(pip, t);
@@ -733,22 +738,28 @@ void BSLMMDAP::DAP_EstimateHyper (const size_t kc, const size_t kd, const vector
}
gsl_vector_set (coef, 0, log(s/(1-s)));
- } else if(kc==0 && kd!=0){//only discrete annotations
+ } else if(kc==0 && kd!=0){
+
+ // Only discrete annotations.
if(kd == 1){
single_ct_regression(Ad, dlevel, pip, coef, prior_vec);
}else{
logistic_cat_fit(coef, Ad, dlevel, pip, 0, 0);
logistic_cat_pred(coef, Ad, dlevel, prior_vec);
}
- } else if (kc!=0 && kd==0) {//only continuous annotations
+ } else if (kc!=0 && kd==0) {
+
+ // Only continuous annotations.
logistic_cont_fit(coef, Ac, pip, 0, 0);
logistic_cont_pred(coef, Ac, prior_vec);
- } else if (kc!=0 && kd!=0) {//both continuous and categorical annotations
+ } else if (kc!=0 && kd!=0) {
+
+ // Both continuous and categorical annotations.
logistic_mixed_fit(coef, Ad, dlevel, Ac, pip, 0, 0);
logistic_mixed_pred(coef, Ad, dlevel, Ac, prior_vec);
}
- //compute marginal likelihood
+ // Compute marginal likelihood.
logm=0;
t1=0;
@@ -757,7 +768,8 @@ void BSLMMDAP::DAP_EstimateHyper (const size_t kc, const size_t kd, const vector
for (size_t m=0; m.
+ along with this program. If not, see .
*/
#ifndef __LMM_H__
@@ -21,22 +21,12 @@
#include "gsl/gsl_vector.h"
#include "gsl/gsl_matrix.h"
-
-
-#ifdef FORCE_FLOAT
-#include "param_float.h"
-#include "io_float.h"
-#else
#include "param.h"
#include "io.h"
-#endif
using namespace std;
-
-
-class FUNC_PARAM
-{
+class FUNC_PARAM {
public:
bool calc_null;
@@ -48,15 +38,12 @@ public:
size_t e_mode;
};
-
-
-
class LMM {
public:
- // IO related parameters
- int a_mode; //analysis mode, 1/2/3/4 for Frequentist tests
- size_t d_pace; //display pace
+ // IO-related parameters
+ int a_mode; // Analysis mode: 1/2/3/4 for Frequentist tests.
+ size_t d_pace; // Display pace.
string file_bfile;
string file_geno;
@@ -75,42 +62,76 @@ public:
double logl_mle_H0;
// Summary statistics
- size_t ni_total, ni_test; //number of individuals
- size_t ns_total, ns_test; //number of snps
- size_t ng_total, ng_test; //number of genes
+ size_t ni_total, ni_test; // Number of individuals.
+ size_t ns_total, ns_test; // Number of SNPs.
+ size_t ng_total, ng_test; // Number of genes.
size_t n_cvt;
- double time_UtX; //time spent on optimization iterations
- double time_opt; //time spent on optimization iterations
+ double time_UtX; // Time spent on optimization iterations.
+ double time_opt; // Time spent on optimization iterations.
- vector indicator_idv; //indicator for individuals (phenotypes), 0 missing, 1 available for analysis
- vector indicator_snp; //sequence indicator for SNPs: 0 ignored because of (a) maf, (b) miss, (c) non-poly; 1 available for analysis
+ // Indicator for individuals (phenotypes): 0 missing, 1
+ // available for analysis.
+ vector indicator_idv;
- vector snpInfo; //record SNP information
+ // Sequence indicator for SNPs: 0 ignored because of (a) maf,
+ // (b) miss, (c) non-poly; 1 available for analysis.
+ vector indicator_snp;
- // Not included in PARAM
- vector sumStat; //Output SNPSummary Data
+ vector snpInfo; // Record SNP information.
- // Main functions
+ // Not included in PARAM.
+ vector sumStat; // Output SNPSummary Data.
+
+ // Main functions.
void CopyFromParam (PARAM &cPar);
void CopyToParam (PARAM &cPar);
- void AnalyzeGene (const gsl_matrix *U, const gsl_vector *eval, const gsl_matrix *UtW, const gsl_vector *Utx, const gsl_matrix *W, const gsl_vector *x);
- void AnalyzePlink (const gsl_matrix *U, const gsl_vector *eval, const gsl_matrix *UtW, const gsl_vector *Uty, const gsl_matrix *W, const gsl_vector *y);
- // WJA added
- void Analyzebgen (const gsl_matrix *U, const gsl_vector *eval, const gsl_matrix *UtW, const gsl_vector *Uty, const gsl_matrix *W, const gsl_vector *y);
- void AnalyzeBimbam (const gsl_matrix *U, const gsl_vector *eval, const gsl_matrix *UtW, const gsl_vector *Uty, const gsl_matrix *W, const gsl_vector *y);
- void AnalyzePlinkGXE (const gsl_matrix *U, const gsl_vector *eval, const gsl_matrix *UtW, const gsl_vector *Uty, const gsl_matrix *W, const gsl_vector *y, const gsl_vector *env);
- void AnalyzeBimbamGXE (const gsl_matrix *U, const gsl_vector *eval, const gsl_matrix *UtW, const gsl_vector *Uty, const gsl_matrix *W, const gsl_vector *y, const gsl_vector *env);
+ void AnalyzeGene (const gsl_matrix *U, const gsl_vector *eval,
+ const gsl_matrix *UtW, const gsl_vector *Utx,
+ const gsl_matrix *W, const gsl_vector *x);
+ void AnalyzePlink (const gsl_matrix *U, const gsl_vector *eval,
+ const gsl_matrix *UtW, const gsl_vector *Uty,
+ const gsl_matrix *W, const gsl_vector *y);
+ // WJA added.
+ void Analyzebgen (const gsl_matrix *U, const gsl_vector *eval,
+ const gsl_matrix *UtW, const gsl_vector *Uty,
+ const gsl_matrix *W, const gsl_vector *y);
+ void AnalyzeBimbam (const gsl_matrix *U, const gsl_vector *eval,
+ const gsl_matrix *UtW, const gsl_vector *Uty,
+ const gsl_matrix *W, const gsl_vector *y);
+ void AnalyzePlinkGXE (const gsl_matrix *U, const gsl_vector *eval,
+ const gsl_matrix *UtW, const gsl_vector *Uty,
+ const gsl_matrix *W, const gsl_vector *y,
+ const gsl_vector *env);
+ void AnalyzeBimbamGXE (const gsl_matrix *U, const gsl_vector *eval,
+ const gsl_matrix *UtW, const gsl_vector *Uty,
+ const gsl_matrix *W, const gsl_vector *y,
+ const gsl_vector *env);
void WriteFiles ();
- void CalcRLWald (const double &lambda, const FUNC_PARAM ¶ms, double &beta, double &se, double &p_wald);
- void CalcRLScore (const double &l, const FUNC_PARAM ¶ms, double &beta, double &se, double &p_score);
+ void CalcRLWald (const double &lambda, const FUNC_PARAM ¶ms,
+ double &beta, double &se, double &p_wald);
+ void CalcRLScore (const double &l, const FUNC_PARAM ¶ms,
+ double &beta, double &se, double &p_score);
};
-void MatrixCalcLR (const gsl_matrix *U, const gsl_matrix *UtX, const gsl_vector *Uty, const gsl_vector *K_eval, const double l_min, const double l_max, const size_t n_region, vector > &pos_loglr);
-void CalcLambda (const char func_name, FUNC_PARAM ¶ms, const double l_min, const double l_max, const size_t n_region, double &lambda, double &logf);
-void CalcLambda (const char func_name, const gsl_vector *eval, const gsl_matrix *UtW, const gsl_vector *Uty, const double l_min, const double l_max, const size_t n_region, double &lambda, double &logl_H0);
-void CalcPve (const gsl_vector *eval, const gsl_matrix *UtW, const gsl_vector *Uty, const double lambda, const double trace_G, double &pve, double &pve_se);
-void CalcLmmVgVeBeta (const gsl_vector *eval, const gsl_matrix *UtW, const gsl_vector *Uty, const double lambda, double &vg, double &ve, gsl_vector *beta, gsl_vector *se_beta);
+void MatrixCalcLR (const gsl_matrix *U, const gsl_matrix *UtX,
+ const gsl_vector *Uty, const gsl_vector *K_eval,
+ const double l_min, const double l_max,
+ const size_t n_region,
+ vector > &pos_loglr);
+void CalcLambda (const char func_name, FUNC_PARAM ¶ms,
+ const double l_min, const double l_max,
+ const size_t n_region, double &lambda, double &logf);
+void CalcLambda (const char func_name, const gsl_vector *eval,
+ const gsl_matrix *UtW, const gsl_vector *Uty,
+ const double l_min, const double l_max,
+ const size_t n_region, double &lambda, double &logl_H0);
+void CalcPve (const gsl_vector *eval, const gsl_matrix *UtW,
+ const gsl_vector *Uty, const double lambda,
+ const double trace_G, double &pve, double &pve_se);
+void CalcLmmVgVeBeta (const gsl_vector *eval, const gsl_matrix *UtW,
+ const gsl_vector *Uty, const double lambda, double &vg,
+ double &ve, gsl_vector *beta, gsl_vector *se_beta);
#endif
diff --git a/src/logistic.cpp b/src/logistic.cpp
index 002ce98..3f9d6ff 100644
--- a/src/logistic.cpp
+++ b/src/logistic.cpp
@@ -13,7 +13,7 @@ typedef struct{
gsl_matrix_int *X;
gsl_vector_int *nlev;
gsl_vector *y;
- gsl_matrix *Xc; // continuous covariates matrix Nobs x Kc (NULL if not used)
+ gsl_matrix *Xc; // Continuous covariates matrix Nobs x Kc (NULL if not used).
double lambdaL1;
double lambdaL2;
} fix_parm_mixed_T;
@@ -29,8 +29,9 @@ double fLogit_mixed(gsl_vector *beta,
int npar = beta->size;
double total = 0;
double aux = 0;
+
// Changed loop start at 1 instead of 0 to avoid regularization of
- // beta_0*\/ */
+ // beta_0*\/
// #pragma omp parallel for reduction (+:total)
for(int i = 1; i < npar; ++i)
total += beta->data[i]*beta->data[i];
@@ -57,8 +58,9 @@ double fLogit_mixed(gsl_vector *beta,
}
-void logistic_mixed_pred(gsl_vector *beta // Vector of parameters length = 1 + Sum_k(C_k - 1)
- ,gsl_matrix_int *X //Matrix Nobs x K
+void logistic_mixed_pred(gsl_vector *beta // Vector of parameters
+ // length = 1 + Sum_k(C_k -1)
+ ,gsl_matrix_int *X //Matrix Nobs x K
,gsl_vector_int *nlev // Vector with number categories
,gsl_matrix *Xc // continuous covariates Matrix Nobs x Kc (NULL if not used)
,gsl_vector *yhat //Vector of prob. predicted by the logistic
@@ -287,39 +289,40 @@ int logistic_mixed_fit(gsl_vector *beta
/* Categorical */
/***************/
-// I need to bundle all the data that goes to the function to optimze together.
+// I need to bundle all the data that goes to the function to optimze
+// together.
typedef struct{
gsl_matrix_int *X;
gsl_vector_int *nlev;
gsl_vector *y;
double lambdaL1;
double lambdaL2;
-}fix_parm_cat_T;
-
-
-double fLogit_cat(gsl_vector *beta
- ,gsl_matrix_int *X
- ,gsl_vector_int *nlev
- ,gsl_vector *y
- ,double lambdaL1
- ,double lambdaL2)
-{
+} fix_parm_cat_T;
+
+double fLogit_cat(gsl_vector *beta,
+ gsl_matrix_int *X,
+ gsl_vector_int *nlev,
+ gsl_vector *y,
+ double lambdaL1,
+ double lambdaL2) {
int n = y->size;
- // int k = X->size2;
int npar = beta->size;
double total = 0;
double aux = 0;
- /* omp_set_num_threads(ompthr); */
- /* /\* Changed loop start at 1 instead of 0 to avoid regularization of beta 0*\/ */
- /* /\*#pragma omp parallel for reduction (+:total)*\/ */
+ // omp_set_num_threads(ompthr); /\* Changed loop start at 1 instead
+ // of 0 to avoid regularization of beta 0*\/ /\*#pragma omp parallel
+ // for reduction (+:total)*\/
for(int i = 1; i < npar; ++i)
total += beta->data[i]*beta->data[i];
total = (-total*lambdaL2/2);
- /* /\*#pragma omp parallel for reduction (+:aux)*\/ */
+
+ // /\*#pragma omp parallel for reduction (+:aux)*\/
for(int i = 1; i < npar; ++i)
aux += (beta->data[i]>0 ? beta->data[i] : -beta->data[i]);
total = total-aux*lambdaL1;
- /* #pragma omp parallel for schedule(static) shared(n,beta,X,nlev,y) reduction (+:total) */
+
+ // #pragma omp parallel for schedule(static) shared(n,beta,X,nlev,y)
+ // #reduction (+:total)
for(int i = 0; i < n; ++i) {
double Xbetai=beta->data[0];
int iParm=1;
@@ -333,7 +336,6 @@ double fLogit_cat(gsl_vector *beta
return -total;
}
-
void logistic_cat_pred(gsl_vector *beta // Vector of parameters length = 1 + Sum_k(C_k - 1)
,gsl_matrix_int *X //Matrix Nobs x K
,gsl_vector_int *nlev // Vector with number categories
diff --git a/src/mvlmm.h b/src/mvlmm.h
index 9ff567c..d495c26 100644
--- a/src/mvlmm.h
+++ b/src/mvlmm.h
@@ -1,6 +1,6 @@
/*
Genome-wide Efficient Mixed Model Association (GEMMA)
- Copyright (C) 2011 Xiang Zhou
+ Copyright (C) 2011-2017, Xiang Zhou
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
@@ -13,36 +13,25 @@
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
- along with this program. If not, see .
- */
+ along with this program. If not, see .
+*/
#ifndef __MVLMM_H__
#define __MVLMM_H__
#include "gsl/gsl_vector.h"
#include "gsl/gsl_matrix.h"
-
-
-#ifdef FORCE_FLOAT
-#include "param_float.h"
-#include "io_float.h"
-#else
#include "param.h"
#include "io.h"
-#endif
using namespace std;
-
-
-
-
class MVLMM {
public:
- // IO related parameters
- int a_mode; //analysis mode, 1/2/3/4 for Frequentist tests
- size_t d_pace; //display pace
+ // IO-related parameters.
+ int a_mode; // Analysis mode: 1/2/3/4 for Frequentist tests.
+ size_t d_pace; // Display pace.
string file_bfile;
string file_geno;
@@ -50,48 +39,68 @@ public:
string file_out;
string path_out;
- // MVLMM related parameters
+ // MVLMM-related parameters.
double l_min;
double l_max;
size_t n_region;
double logl_remle_H0, logl_mle_H0;
vector Vg_remle_null, Ve_remle_null, Vg_mle_null, Ve_mle_null;
- vector VVg_remle_null, VVe_remle_null, VVg_mle_null, VVe_mle_null;
- vector beta_remle_null, se_beta_remle_null, beta_mle_null, se_beta_mle_null;
+ vector VVg_remle_null, VVe_remle_null, VVg_mle_null;
+ vector VVe_mle_null;
+ vector beta_remle_null, se_beta_remle_null, beta_mle_null;
+ vector se_beta_mle_null;
double p_nr;
size_t em_iter, nr_iter;
double em_prec, nr_prec;
size_t crt;
- // Summary statistics
- size_t ni_total, ni_test; //number of individuals
- size_t ns_total, ns_test; //number of snps
+ // Summary statistics.
+ size_t ni_total, ni_test; // Number of individuals.
+ size_t ns_total, ns_test; // Number of SNPs.
size_t n_cvt;
size_t n_ph;
- double time_UtX; //time spent on optimization iterations
- double time_opt; //time spent on optimization iterations
+ double time_UtX; // Time spent on optimization iterations.
+ double time_opt; // Time spent on optimization iterations.
+
+ // Indicator for individuals (phenotypes): 0 missing, 1
+ // available for analysis.
+ vector indicator_idv;
- vector indicator_idv; //indicator for individuals (phenotypes), 0 missing, 1 available for analysis
- vector indicator_snp; //sequence indicator for SNPs: 0 ignored because of (a) maf, (b) miss, (c) non-poly; 1 available for analysis
+ // Sequence indicator for SNPs: 0 ignored because of (a) maf,
+ // (b) miss, (c) non-poly; 1 available for analysis.
+ vector indicator_snp;
- vector snpInfo; //record SNP information
+ vector snpInfo; // Record SNP information.
- // Not included in PARAM
- vector sumStat; //Output SNPSummary Data
+ // Not included in PARAM.
+ vector sumStat; // Output SNPSummary Data.
// Main functions
void CopyFromParam (PARAM &cPar);
void CopyToParam (PARAM &cPar);
- void AnalyzeBimbam (const gsl_matrix *U, const gsl_vector *eval, const gsl_matrix *UtW, const gsl_matrix *UtY);
- void AnalyzePlink (const gsl_matrix *U, const gsl_vector *eval, const gsl_matrix *UtW, const gsl_matrix *UtY);
- void Analyzebgen (const gsl_matrix *U, const gsl_vector *eval, const gsl_matrix *UtW, const gsl_matrix *UtY);
- void AnalyzeBimbamGXE (const gsl_matrix *U, const gsl_vector *eval, const gsl_matrix *UtW, const gsl_matrix *UtY, const gsl_vector *env);
- void AnalyzePlinkGXE (const gsl_matrix *U, const gsl_vector *eval, const gsl_matrix *UtW, const gsl_matrix *UtY, const gsl_vector *env);
+ void AnalyzeBimbam (const gsl_matrix *U, const gsl_vector *eval,
+ const gsl_matrix *UtW, const gsl_matrix *UtY);
+ void AnalyzePlink (const gsl_matrix *U, const gsl_vector *eval,
+ const gsl_matrix *UtW, const gsl_matrix *UtY);
+ void Analyzebgen (const gsl_matrix *U, const gsl_vector *eval,
+ const gsl_matrix *UtW, const gsl_matrix *UtY);
+ void AnalyzeBimbamGXE (const gsl_matrix *U, const gsl_vector *eval,
+ const gsl_matrix *UtW, const gsl_matrix *UtY,
+ const gsl_vector *env);
+ void AnalyzePlinkGXE (const gsl_matrix *U, const gsl_vector *eval,
+ const gsl_matrix *UtW, const gsl_matrix *UtY,
+ const gsl_vector *env);
void WriteFiles ();
};
-void CalcMvLmmVgVeBeta (const gsl_vector *eval, const gsl_matrix *UtW, const gsl_matrix *UtY, const size_t em_iter, const size_t nr_iter, const double em_prec, const double nr_prec, const double l_min, const double l_max, const size_t n_region, gsl_matrix *V_g, gsl_matrix *V_e, gsl_matrix *B, gsl_matrix *se_B);
+void CalcMvLmmVgVeBeta (const gsl_vector *eval, const gsl_matrix *UtW,
+ const gsl_matrix *UtY, const size_t em_iter,
+ const size_t nr_iter, const double em_prec,
+ const double nr_prec, const double l_min,
+ const double l_max, const size_t n_region,
+ gsl_matrix *V_g, gsl_matrix *V_e, gsl_matrix *B,
+ gsl_matrix *se_B);
#endif
--
cgit v1.2.3