about summary refs log tree commit diff
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;
 }
generated by cgit-pink 1.4.1 (git 2.36.1) at 2025-11-20 18:41:44 +0000