From f28b8531d135a948a5858844c6b6ce2bbfcf6813 Mon Sep 17 00:00:00 2001 From: Peter Carbonetto Date: Wed, 31 May 2017 15:12:19 -0500 Subject: Removed WITH_LAPACK compile/build option. --- src/bslmm.cpp | 173 +++---------------------------------------- src/bslmm.h | 16 +--- src/bslmmdap.cpp | 220 +------------------------------------------------------ src/bslmmdap.h | 18 +---- src/io.cpp | 139 +++++++++++++++++++++++------------ src/io.h | 72 +++++++++++++----- src/lapack.cpp | 94 ------------------------ src/param.cpp | 2 +- src/param.h | 2 +- 9 files changed, 172 insertions(+), 564 deletions(-) (limited to 'src') diff --git a/src/bslmm.cpp b/src/bslmm.cpp index 55a05ca..0bd4234 100644 --- a/src/bslmm.cpp +++ b/src/bslmm.cpp @@ -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,7 +13,7 @@ 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 . */ #include @@ -38,30 +38,15 @@ #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; @@ -330,11 +315,8 @@ double BSLMM::CalcPveLM (const gsl_matrix *UtXgamma, const gsl_vector *Uty, cons 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 + lapack_dgemm ((char *)"T", (char *)"N", 1.0, UtXgamma, UtXgamma, + 1.0, Omega); gsl_blas_dgemv (CblasTrans, 1.0, UtXgamma, Uty, 0.0, Xty); CholeskySolve(Omega, Xty, OiXty); @@ -517,11 +499,8 @@ double BSLMM::CalcPosterior (const gsl_matrix *UtXgamma, const gsl_vector *Uty, 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 + lapack_dgemm ((char *)"T", (char *)"N", sigma_a2, UtXgamma_eval, + UtXgamma, 1.0, Omega); time_Omega+=(clock()-time_start)/(double(CLOCKS_PER_SEC)*60.0); @@ -1229,11 +1208,8 @@ void BSLMM::CalcXtX (const gsl_matrix *X, const gsl_vector *y, const size_t s_si 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 + lapack_dgemm ((char *)"T", (char *)"N", 1.0, &X_sub.matrix, + &X_sub.matrix, 0.0, &XtX_sub.matrix); 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); @@ -1282,127 +1258,6 @@ void BSLMM::SetXgamma (const gsl_matrix *X, const gsl_matrix *X_old, const gsl_m gsl_vector_view Xtynew_sub=gsl_vector_subvector(Xty_new, 0, rank_new.size()); //get X_new and calculate XtX_new - /* - if (rank_remove.size()==0 && rank_add.size()==0) { - gsl_matrix_memcpy(&Xnew_sub.matrix, &Xold_sub.matrix); - gsl_matrix_memcpy(&XtXnew_sub.matrix, &XtXold_sub.matrix); - gsl_vector_memcpy(&Xtynew_sub.vector, &Xtyold_sub.vector); - } else { - gsl_matrix *X_temp=gsl_matrix_alloc(X_old->size1, rank_old.size()-rank_remove.size() ); - gsl_matrix *XtX_temp=gsl_matrix_alloc(X_temp->size2, X_temp->size2); - gsl_vector *Xty_temp=gsl_vector_alloc(X_temp->size2); - - if (rank_remove.size()==0) { - gsl_matrix_memcpy (X_temp, &Xold_sub.matrix); - gsl_matrix_memcpy (XtX_temp, &XtXold_sub.matrix); - gsl_vector_memcpy (Xty_temp, &Xtyold_sub.vector); - } else { - size_t i_temp=0, j_temp; - for (size_t i=0; isize1, rank_add.size() ); - gsl_matrix *XtX_aa=gsl_matrix_alloc(X_add->size2, X_add->size2); - gsl_matrix *XtX_at=gsl_matrix_alloc(X_add->size2, X_temp->size2); - gsl_vector *Xty_add=gsl_vector_alloc(X_add->size2); - - //get X_add - SetXgamma (X_add, X, rank_add); - - //get t(X_add)X_add and t(X_add)X_temp - clock_t time_start=clock(); - - //somehow the lapack_dgemm does not work here - //#ifdef WITH_LAPACK - //lapack_dgemm ((char *)"T", (char *)"N", 1.0, X_add, X_add, 0.0, XtX_aa); - //lapack_dgemm ((char *)"T", (char *)"N", 1.0, X_add, X_temp, 0.0, XtX_at); - - //#else - gsl_blas_dgemm (CblasTrans, CblasNoTrans, 1.0, X_add, X_add, 0.0, XtX_aa); - gsl_blas_dgemm (CblasTrans, CblasNoTrans, 1.0, X_add, X_temp, 0.0, XtX_at); - //#endif - gsl_blas_dgemv(CblasTrans, 1.0, X_add, y, 0.0, Xty_add); - - time_Omega+=(clock()-time_start)/(double(CLOCKS_PER_SEC)*60.0); - - //save to X_new, XtX_new and Xty_new - size_t i_temp=0, j_temp, i_flag=0, j_flag=0; - for (size_t i=0; i. + along with this program. If not, see . */ @@ -25,24 +25,14 @@ #include #include -#ifdef FORCE_FLOAT -#include "param_float.h" -#else #include "param.h" -#endif - using namespace std; - - - - - class BSLMM { public: - // IO related parameters + // IO-related parameters. int a_mode; size_t d_pace; diff --git a/src/bslmmdap.cpp b/src/bslmmdap.cpp index f6e9e9c..c66189b 100644 --- a/src/bslmmdap.cpp +++ b/src/bslmmdap.cpp @@ -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,7 +13,7 @@ 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 . */ #include @@ -38,31 +38,17 @@ #include "gsl/gsl_cdf.h" #include "gsl/gsl_roots.h" - - #include "logistic.h" #include "lapack.h" #include "io.h" - -#ifdef FORCE_FLOAT -#include "param_float.h" -#include "bslmmdap_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 "bslmmdap.h" #include "lmm.h" #include "lm.h" #include "mathfunc.h" -#endif using namespace std; - - - void BSLMMDAP::CopyFromParam (PARAM &cPar) { file_out=cPar.file_out; @@ -506,11 +492,8 @@ double BSLMMDAP::CalcMarginal (const gsl_matrix *UtXgamma, const gsl_vector *Uty 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 + lapack_dgemm ((char *)"T", (char *)"N", sigma_a2, UtXgamma_eval, + UtXgamma, 1.0, Omega); time_Omega+=(clock()-time_start)/(double(CLOCKS_PER_SEC)*60.0); //calculate beta_hat @@ -819,198 +802,3 @@ void BSLMMDAP::DAP_EstimateHyper (const size_t kc, const size_t kd, const vector gsl_vector_free(coef); return; } - -/* -//readin the estimated hyper-parameters and perform fine mapping for each region -void BSLMM::DAP_FineMapping (const gsl_matrix *U, const gsl_matrix *UtX, const gsl_matrix *A, const gsl_vector *Uty, const gsl_vector *K_eval, const gsl_vector *y, gsl_matrix *Hyper, gsl_vector *alpha, gsl_vector *pip) { - clock_t time_start; - - //two priority sets: S_1 contains all candidate causal SNPs; S_2 contains the prioritized combintion of them - //two marginal probability sets: P_1 contains marginals for S_1; P_2 contains marginals for S_2; - set S1set, S2set; - vector S1vec; - vector > S2vec; - vector P1, P2; - - //calculate P0 (null) and P1 (for every SNP) - - - - //loop through the number of combinations - for (size_t s=0; s10) {break;} - - if (a_mode==13) { - SampleZ (y, z_hat, z); - mean_z=CenterVector (z); - - time_start=clock(); - gsl_blas_dgemv (CblasTrans, 1.0, U, z, 0.0, Utz); - time_UtZ+=(clock()-time_start)/(double(CLOCKS_PER_SEC)*60.0); - - //First proposal - if (cHyp_old.n_gamma==0 || cHyp_old.rho==0) { - logPost_old=CalcPosterior(Utz, K_eval, Utu_old, alpha_old, cHyp_old); - beta_old.clear(); - for (size_t i=0; isize; ++i) { - beta_old.push_back(gsl_vector_get(beta, i)); - } - gsl_matrix_free (UtXgamma); - gsl_vector_free (beta); - } - } - - - delete [] p_gamma; - beta_g.clear(); - - 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; -} - - - -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. + along with this program. If not, see . */ - #ifndef __BSLMMDAP_H__ #define __BSLMMDAP_H__ @@ -24,25 +23,14 @@ #include #include #include - -#ifdef FORCE_FLOAT -#include "param_float.h" -#else #include "param.h" -#endif - using namespace std; - - - - - class BSLMMDAP { public: - // IO related parameters + // IO-related parameters. int a_mode; size_t d_pace; diff --git a/src/io.cpp b/src/io.cpp index 1a4d466..e8f00ed 100644 --- a/src/io.cpp +++ b/src/io.cpp @@ -1,6 +1,6 @@ /* Genome-wide Efficient Mixed Model Association (GEMMA) - Copyright (C) 2011-2017 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 @@ -194,11 +194,13 @@ bool ReadFile_snps_header (const string &file_snps, set &setSnps) { return true; } -//Read log file -bool ReadFile_log (const string &file_log, double &pheno_mean) -{ +// Read log file. +bool ReadFile_log (const string &file_log, double &pheno_mean) { ifstream infile (file_log.c_str(), ifstream::in); - if (!infile) {cout<<"error! fail to open log file: "< &mapRS2chr, map &mapRS2bp, map &mapRS2cM) -{ +// Read bimbam annotation file. +bool ReadFile_anno (const string &file_anno, map &mapRS2chr, + map &mapRS2bp, + map &mapRS2cM) { mapRS2chr.clear(); mapRS2bp.clear(); ifstream infile (file_anno.c_str(), ifstream::in); - if (!infile) {cout<<"error opening annotation file: "< &mapRS2chr, map ch_ptr=strtok ((char *)line.c_str(), " , \t"); rs=ch_ptr; ch_ptr=strtok (NULL, " , \t"); - if (strcmp(ch_ptr, "NA")==0) {b_pos=-9;} else {b_pos=atol(ch_ptr);} + if (strcmp(ch_ptr, "NA")==0) { + b_pos=-9; + } else { + b_pos=atol(ch_ptr); + } ch_ptr=strtok (NULL, " , \t"); - if (ch_ptr==NULL || strcmp(ch_ptr, "NA")==0) {chr="-9";} else {chr=ch_ptr;} + if (ch_ptr==NULL || strcmp(ch_ptr, "NA")==0) { + chr="-9"; + } else { + chr=ch_ptr; + } ch_ptr=strtok (NULL, " , \t"); - if (ch_ptr==NULL || strcmp(ch_ptr, "NA")==0) {cM=-9;} else {cM=atof(ch_ptr);} + if (ch_ptr==NULL || strcmp(ch_ptr, "NA")==0) { + cM=-9; + } else { + cM=atof(ch_ptr); + } mapRS2chr[rs]=chr; mapRS2bp[rs]=b_pos; @@ -268,15 +285,17 @@ bool ReadFile_anno (const string &file_anno, map &mapRS2chr, map return true; } -//read one column of phenotype -bool ReadFile_column (const string &file_pheno, vector &indicator_idv, vector &pheno, const int &p_column) -{ +// Read 1 column of phenotype. +bool ReadFile_column (const string &file_pheno, vector &indicator_idv, + vector &pheno, const int &p_column) { indicator_idv.clear(); pheno.clear(); igzstream infile (file_pheno.c_str(), igzstream::in); -// ifstream infile (file_pheno.c_str(), ifstream::in); - if (!infile) {cout<<"error! fail to open phenotype file: "< &indicator_idv, vect for (int i=0; i<(p_column-1); ++i) { ch_ptr=strtok (NULL, " , \t"); } - if (strcmp(ch_ptr, "NA")==0) {indicator_idv.push_back(0); pheno.push_back(-9);} //pheno is different from pimass2 - else {p=atof(ch_ptr); indicator_idv.push_back(1); pheno.push_back(p);} + if (strcmp(ch_ptr, "NA")==0) { + indicator_idv.push_back(0); + pheno.push_back(-9); + } // Pheno is different from pimass2. + else { + p=atof(ch_ptr); + indicator_idv.push_back(1); + pheno.push_back(p); + } } infile.close(); @@ -298,17 +324,19 @@ bool ReadFile_column (const string &file_pheno, vector &indicator_idv, vect return true; } - - -//Read bimbam phenotype file, p_column=1, 2 ... -bool ReadFile_pheno (const string &file_pheno, vector > &indicator_pheno, vector > &pheno, const vector &p_column) -{ +// Read bimbam phenotype file, p_column=1, 2,... +bool ReadFile_pheno (const string &file_pheno, + vector > &indicator_pheno, + vector > &pheno, + const vector &p_column) { indicator_pheno.clear(); pheno.clear(); igzstream infile (file_pheno.c_str(), igzstream::in); -// ifstream infile (file_pheno.c_str(), ifstream::in); - if (!infile) {cout<<"error! fail to open phenotype file: "< > &indicator_p size_t i=0; while (i > &indicator_p return true; } - -bool ReadFile_cvt (const string &file_cvt, vector &indicator_cvt, vector > &cvt, size_t &n_cvt) -{ +bool ReadFile_cvt (const string &file_cvt, vector &indicator_cvt, + vector > &cvt, size_t &n_cvt) { indicator_cvt.clear(); ifstream infile (file_cvt.c_str(), ifstream::in); - if (!infile) {cout<<"error! fail to open covariates file: "< &indicator_cvt, vector &indicator_cvt, vector::size_type i=0; i &indicator_cvt, vector &snpInfo) -{ - snpInfo.clear(); +// Read .bim file. +bool ReadFile_bim (const string &file_bim, vector &snpInfo) { + snpInfo.clear(); ifstream infile (file_bim.c_str(), ifstream::in); - if (!infile) {cout<<"error opening .bim file: "< &snpInfo) ch_ptr=strtok (NULL, " \t"); major=ch_ptr; - SNPINFO sInfo={chr, rs, cM, b_pos, minor, major, 0, -9, -9, 0, 0, 0}; + SNPINFO sInfo={chr, rs, cM, b_pos, minor, major, + 0, -9, -9, 0, 0, 0}; snpInfo.push_back(sInfo); } @@ -438,8 +486,7 @@ bool ReadFile_bim (const string &file_bim, vector &snpInfo) return true; } - -//Read .fam file +// Read .fam file. bool ReadFile_fam (const string &file_fam, vector > &indicator_pheno, vector > &pheno, map &mapID2num, const vector &p_column) { indicator_pheno.clear(); diff --git a/src/io.h b/src/io.h index 365127a..9d6f8cc 100644 --- a/src/io.h +++ b/src/io.h @@ -1,6 +1,6 @@ /* Genome-wide Efficient Mixed Model Association (GEMMA) - Copyright (C) 2011-2017 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 @@ -121,22 +121,60 @@ bool ReadFile_gene (const string &file_gene, vector &vec_read, vector &snpInfo, size_t &ng_total); bool ReadHeader_io (const string &line, HEADER &header); -bool ReadFile_cat (const string &file_cat, map &mapRS2cat, size_t &n_vc); -bool ReadFile_mcat (const string &file_mcat, map &mapRS2cat, size_t &n_vc); - -bool ReadFile_catc (const string &file_cat, map > &mapRS2catc, size_t &n_cat); -bool ReadFile_mcatc (const string &file_mcat, map > &mapRS2catc, size_t &n_cat); - -bool BimbamKin (const string &file_geno, const int display_pace, const vector &indicator_idv, const vector &indicator_snp, const map &mapRS2weight, const map &mapRS2cat, const vector &snpInfo, const gsl_matrix *W, gsl_matrix *matrix_kin, gsl_vector *vector_ns); -bool PlinkKin (const string &file_bed, const int display_pace, const vector &indicator_idv, const vector &indicator_snp, const map &mapRS2weight, const map &mapRS2cat, const vector &snpInfo, const gsl_matrix *W, gsl_matrix *matrix_kin, gsl_vector *vector_ns); -bool MFILEKin (const size_t mfile_mode, const string &file_mfile, const int display_pace, const vector &indicator_idv, const vector > &mindicator_snp, const map &mapRS2weight, const map &mapRS2cat, const vector > &msnpInfo, const gsl_matrix *W, gsl_matrix *matrix_kin, gsl_vector *vector_ns); - -bool ReadFile_wsnp (const string &file_wsnp, map &mapRS2double); -bool ReadFile_wsnp (const string &file_wcat, const size_t n_vc, map > &mapRS2vector); - -void ReadFile_beta (const string &file_beta, const map &mapRS2cat, const map &mapRS2wA, vector &vec_cat, vector &vec_ni, vector &vec_weight, vector &vec_z2, size_t &ni_total, size_t &ns_total, size_t &ns_test); -void ReadFile_beta (const string &file_beta, const map &mapRS2wA, map &mapRS2A1, map &mapRS2z); -void Calcq (const size_t n_block, const vector &vec_cat, const vector &vec_ni, const vector &vec_weight, const vector &vec_z2, gsl_matrix *Vq, gsl_vector *q, gsl_vector *s); +bool ReadFile_cat (const string &file_cat, map &mapRS2cat, + size_t &n_vc); +bool ReadFile_mcat (const string &file_mcat, map &mapRS2cat, + size_t &n_vc); + +bool ReadFile_catc (const string &file_cat, + map > &mapRS2catc, + size_t &n_cat); +bool ReadFile_mcatc (const string &file_mcat, map > &mapRS2catc, size_t &n_cat); + +bool BimbamKin (const string &file_geno, const int display_pace, + const vector &indicator_idv, + const vector &indicator_snp, + const map &mapRS2weight, + const map &mapRS2cat, + const vector &snpInfo, const gsl_matrix *W, + gsl_matrix *matrix_kin, gsl_vector *vector_ns); +bool PlinkKin (const string &file_bed, const int display_pace, + const vector &indicator_idv, + const vector &indicator_snp, + const map &mapRS2weight, + const map &mapRS2cat, + const vector &snpInfo, + const gsl_matrix *W, gsl_matrix *matrix_kin, + gsl_vector *vector_ns); +bool MFILEKin (const size_t mfile_mode, const string &file_mfile, + const int display_pace, const vector &indicator_idv, + const vector > &mindicator_snp, + const map &mapRS2weight, + const map &mapRS2cat, + const vector > &msnpInfo, + const gsl_matrix *W, gsl_matrix *matrix_kin, + gsl_vector *vector_ns); + +bool ReadFile_wsnp (const string &file_wsnp, + map &mapRS2double); +bool ReadFile_wsnp (const string &file_wcat, const size_t n_vc, + map > &mapRS2vector); + +void ReadFile_beta (const string &file_beta, + const map &mapRS2cat, + const map &mapRS2wA, + vector &vec_cat, vector &vec_ni, + vector &vec_weight, vector &vec_z2, + size_t &ni_total, size_t &ns_total, size_t &ns_test); +void ReadFile_beta (const string &file_beta, + const map &mapRS2wA, + map &mapRS2A1, + map &mapRS2z); +void Calcq (const size_t n_block, const vector &vec_cat, + const vector &vec_ni, + const vector &vec_weight, const vector &vec_z2, + gsl_matrix *Vq, gsl_vector *q, gsl_vector *s); void ReadFile_study (const string &file_study, gsl_matrix *Vq, gsl_vector *q_vec, gsl_vector *s_vec, size_t &ni); diff --git a/src/lapack.cpp b/src/lapack.cpp index 18399a2..2bbdf62 100644 --- a/src/lapack.cpp +++ b/src/lapack.cpp @@ -398,13 +398,7 @@ void lapack_eigen_symmv (gsl_matrix *A, gsl_vector *eval, gsl_matrix *evec, // DO NOT set eigenvalues to be positive. double EigenDecomp (gsl_matrix *G, gsl_matrix *U, gsl_vector *eval, const size_t flag_largematrix) { -#ifdef WITH_LAPACK lapack_eigen_symmv (G, eval, U, flag_largematrix); -#else - gsl_eigen_symmv_workspace *w=gsl_eigen_symmv_alloc (G->size1); - gsl_eigen_symmv (G, eval, U, w); - gsl_eigen_symmv_free (w); -#endif // Calculate track_G=mean(diag(G)). double d=0.0; @@ -420,48 +414,7 @@ double EigenDecomp (gsl_matrix *G, gsl_matrix *U, gsl_vector *eval, // DO NOT set eigen values to be positive. double EigenDecomp (gsl_matrix_float *G, gsl_matrix_float *U, gsl_vector_float *eval, const size_t flag_largematrix) { -#ifdef WITH_LAPACK lapack_float_eigen_symmv (G, eval, U, flag_largematrix); -#else - - // GSL doesn't provide single precision eigen decomposition; - // plus, float precision eigen decomposition in LAPACK may not - // work on OS 10.4 first change to double precision. - gsl_matrix *G_double=gsl_matrix_alloc (G->size1, G->size2); - gsl_matrix *U_double=gsl_matrix_alloc (U->size1, U->size2); - gsl_vector *eval_double=gsl_vector_alloc (eval->size); - for (size_t i=0; isize1; i++) { - for (size_t j=0; jsize2; j++) { - gsl_matrix_set(G_double, i, j, - gsl_matrix_float_get(G, i, j)); - } - } - gsl_eigen_symmv_workspace *w_space=gsl_eigen_symmv_alloc (G->size1); - gsl_eigen_symmv (G_double, eval_double, U_double, w_space); - gsl_eigen_symmv_free (w_space); - - // Change back to float precision. - for (size_t i=0; isize1; i++) { - for (size_t j=0; jsize2; j++) { - gsl_matrix_float_set(K, i, j, - gsl_matrix_get(G_double, i, j)); - } - } - for (size_t i=0; isize1; i++) { - for (size_t j=0; jsize2; j++) { - gsl_matrix_float_set(U, i, j, - gsl_matrix_get(U_double, i, j)); - } - } - for (size_t i=0; isize; i++) { - gsl_vector_float_set(eval, i, gsl_vector_get(eval_double, i)); - } - - // Delete double-precision matrices. - gsl_matrix_free (G_double); - gsl_matrix_free (U_double); - gsl_vector_free (eval_double); -#endif // Calculate track_G=mean(diag(G)). double d = 0.0; @@ -477,28 +430,12 @@ double EigenDecomp (gsl_matrix_float *G, gsl_matrix_float *U, double CholeskySolve(gsl_matrix *Omega, gsl_vector *Xty, gsl_vector *OiXty) { double logdet_O=0.0; -#ifdef WITH_LAPACK lapack_cholesky_decomp(Omega); for (size_t i=0; isize1; ++i) { logdet_O+=log(gsl_matrix_get (Omega, i, i)); } logdet_O*=2.0; lapack_cholesky_solve(Omega, Xty, OiXty); -#else - int status = gsl_linalg_cholesky_decomp(Omega); - if(status == GSL_EDOM) { - cout << "## non-positive definite matrix" << endl; - } - - for (size_t i=0; isize1; ++i) { - logdet_O+=log(gsl_matrix_get (Omega, i, i)); - } - logdet_O*=2.0; - - gsl_vector_memcpy (OiXty, Xty); - gsl_blas_dtrsv(CblasLower, CblasNoTrans, CblasNonUnit, Omega, OiXty); - gsl_blas_dtrsv(CblasUpper, CblasNoTrans, CblasNonUnit, Omega, OiXty); -#endif return logdet_O; } @@ -508,43 +445,12 @@ double CholeskySolve(gsl_matrix_float *Omega, gsl_vector_float *Xty, gsl_vector_float *OiXty) { double logdet_O=0.0; -#ifdef WITH_LAPACK lapack_float_cholesky_decomp(Omega); for (size_t i=0; isize1; ++i) { logdet_O+=log(gsl_matrix_float_get (Omega, i, i)); } logdet_O*=2.0; lapack_float_cholesky_solve(Omega, Xty, OiXty); -#else - gsl_matrix *Omega_double=gsl_matrix_alloc (Omega->size1, Omega->size2); - double d; - for (size_t i=0; isize1; ++i) { - for (size_t j=0; jsize2; ++j) { - d=(double)gsl_matrix_float_get (Omega, i, j); - gsl_matrix_set (Omega_double, i, j, d); - } - } - - int status = gsl_linalg_cholesky_decomp(Omega_double); - if(status == GSL_EDOM) { - cout << "## non-positive definite matrix" << endl; - } - - for (size_t i=0; isize1; ++i) { - for (size_t j=0; jsize2; ++j) { - d=gsl_matrix_get (Omega_double, i, j); - if (j==i) {logdet_O+=log(d);} - gsl_matrix_float_set (Omega, i, j, (float)d); - } - } - logdet_O*=2.0; - - gsl_vector_float_memcpy (OiXty, Xty); - gsl_blas_strsv(CblasLower, CblasNoTrans, CblasNonUnit, Omega, OiXty); - gsl_blas_strsv(CblasUpper, CblasNoTrans, CblasNonUnit, Omega, OiXty); - - gsl_matrix_free (Omega_double); -#endif return logdet_O; } diff --git a/src/param.cpp b/src/param.cpp index 63bf349..0a1d721 100644 --- a/src/param.cpp +++ b/src/param.cpp @@ -1,6 +1,6 @@ /* Genome-wide Efficient Mixed Model Association (GEMMA) - Copyright (C) 2011-2017 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 diff --git a/src/param.h b/src/param.h index 18d5e36..9707790 100644 --- a/src/param.h +++ b/src/param.h @@ -1,6 +1,6 @@ /* Genome-wide Efficient Mixed Model Association (GEMMA) - Copyright (C) 2011-2017 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 -- cgit v1.2.3