diff options
Diffstat (limited to 'src/ldr.cpp')
-rw-r--r-- | src/ldr.cpp | 207 |
1 files changed, 16 insertions, 191 deletions
diff --git a/src/ldr.cpp b/src/ldr.cpp index e28490a..a1e5791 100644 --- a/src/ldr.cpp +++ b/src/ldr.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,8 +13,8 @@ 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/>. - */ + along with this program. If not, see <http://www.gnu.org/licenses/>. +*/ #include <iostream> #include <fstream> @@ -39,29 +39,16 @@ #include "gsl/gsl_roots.h" #include "Eigen/Dense" - - #include "lapack.h" - -#ifdef FORCE_FLOAT -#include "param_float.h" -#include "ldr_float.h" -#include "lm_float.h" -#include "mathfunc_float.h" //for function CenterVector -#else #include "param.h" #include "ldr.h" #include "lm.h" #include "mathfunc.h" -#endif using namespace std; using namespace Eigen; - - -void LDR::CopyFromParam (PARAM &cPar) -{ +void LDR::CopyFromParam (PARAM &cPar) { a_mode=cPar.a_mode; d_pace=cPar.d_pace; @@ -83,176 +70,15 @@ void LDR::CopyFromParam (PARAM &cPar) return; } - -void LDR::CopyToParam (PARAM &cPar) -{ - //cPar.pheno_mean=pheno_mean; - //cPar.randseed=randseed; - - return; -} - - -/* -void BSLMM::WriteBV (const gsl_vector *bv) -{ - string file_str; - file_str=path_out+"/"+file_out; - file_str+=".bv.txt"; - - ofstream outfile (file_str.c_str(), ofstream::out); - if (!outfile) {cout<<"error writing file: "<<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(); +void LDR::CopyToParam (PARAM &cPar) { return; } - - - -void BSLMM::WriteParam (vector<pair<double, double> > &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: "<<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=path_out+"/"+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=path_out+"/"+file_out; - file_gamma+=".gamma.txt"; - file_hyp=path_out+"/"+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; -} - -*/ - -//X is a p by n matrix -void LDR::VB (const vector<vector<unsigned char> > &Xt, const gsl_matrix *W_gsl, const gsl_vector *y_gsl) -{ - //save gsl_vector and gsl_matrix into eigen library formats +//X is a p by n matrix. +void LDR::VB (const vector<vector<unsigned char> > &Xt, + const gsl_matrix *W_gsl, const gsl_vector *y_gsl) { + + // Save gsl_vector and gsl_matrix into Eigen library formats. MatrixXd W(W_gsl->size1, W_gsl->size2); VectorXd y(y_gsl->size); VectorXd x_col(y_gsl->size); @@ -266,7 +92,7 @@ void LDR::VB (const vector<vector<unsigned char> > &Xt, const gsl_matrix *W_gsl, } } - //initial VB values by lm + // Initial VB values by lm. cout<<indicator_snp[0]<<" "<<indicator_snp[1]<<" "<<indicator_snp[2]<<endl; uchar_matrix_get_row (Xt, 0, x_col); @@ -274,12 +100,11 @@ void LDR::VB (const vector<vector<unsigned char> > &Xt, const gsl_matrix *W_gsl, cout<<x_col(j)<<endl; } + // Run VB iterations. + // TO DO. - //run VB iterations - - - - //save results - + // Save results. + // TO DO. + return; } |