aboutsummaryrefslogtreecommitdiff
path: root/src/ldr.cpp
diff options
context:
space:
mode:
Diffstat (limited to 'src/ldr.cpp')
-rw-r--r--src/ldr.cpp207
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;
}