version 0.95alpha

15 files changed, 43 insertions, 22 deletions

diff --git a/bslmm.cpp b/bslmm.cpp
index ff9618d..55a05ca 100644
--- a/bslmm.cpp
+++ b/bslmm.cpp
@@ -70,6 +70,7 @@ void BSLMM::CopyFromParam (PARAM &cPar)
 	file_bfile=cPar.file_bfile;
 	file_geno=cPar.file_geno;
 	file_out=cPar.file_out;
+	path_out=cPar.path_out;
 	
 	l_min=cPar.h_min;	
 	l_max=cPar.h_max;  
@@ -134,7 +135,7 @@ void BSLMM::CopyToParam (PARAM &cPar)
 void BSLMM::WriteBV (const gsl_vector *bv) 
 {
 	string file_str;
-	file_str="./output/"+file_out;
+	file_str=path_out+"/"+file_out;
 	file_str+=".bv.txt";
 
 	ofstream outfile (file_str.c_str(), ofstream::out);
@@ -162,7 +163,7 @@ void BSLMM::WriteBV (const gsl_vector *bv)
 void BSLMM::WriteParam (vector<pair<double, double> > &beta_g, const gsl_vector *alpha, const size_t w) 
 {
 	string file_str;
-	file_str="./output/"+file_out;
+	file_str=path_out+"/"+file_out;
 	file_str+=".param.txt";
 
 	ofstream outfile (file_str.c_str(), ofstream::out);
@@ -198,7 +199,7 @@ void BSLMM::WriteParam (vector<pair<double, double> > &beta_g, const gsl_vector
 void BSLMM::WriteParam (const gsl_vector *alpha) 
 {
 	string file_str;
-	file_str="./output/"+file_out;
+	file_str=path_out+"/"+file_out;
 	file_str+=".param.txt";
 
 	ofstream outfile (file_str.c_str(), ofstream::out);
@@ -228,9 +229,9 @@ void BSLMM::WriteParam (const gsl_vector *alpha)
 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="./output/"+file_out;
+	file_gamma=path_out+"/"+file_out;
 	file_gamma+=".gamma.txt";
-	file_hyp="./output/"+file_out;
+	file_hyp=path_out+"/"+file_out;
 	file_hyp+=".hyp.txt";
 
 	ofstream outfile_gamma, outfile_hyp;
@@ -1448,7 +1449,7 @@ void BSLMM::SetXgamma (const gsl_matrix *X, const gsl_matrix *X_old, const gsl_m
       //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_old, 0.0, XtX_ao);
+      //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);
diff --git a/bslmm.h b/bslmm.h
index f407794..8b5edc7 100644
--- a/bslmm.h
+++ b/bslmm.h
@@ -49,6 +49,7 @@ public:
 	string file_bfile;
 	string file_geno;
 	string file_out;
+	string path_out;
 	
 	// LMM related parameters
 	double l_min;
diff --git a/gemma.cpp b/gemma.cpp
index 093cd05..b8693a8 100644
--- a/gemma.cpp
+++ b/gemma.cpp
@@ -61,7 +61,7 @@ using namespace std;
 
 
 GEMMA::GEMMA(void):	
-version("0.95"), date("08/08/2014"), year("2011")
+version("0.95alpha"), date("08/08/2014"), year("2011")
 {}
 
 void GEMMA::PrintHeader (void)
@@ -238,6 +238,7 @@ void GEMMA::PrintHelp(size_t option)
 		cout<<"                   2: \"id  id  value\" format"<<endl;
 		cout<<" -n        [num]          "<<" specify phenotype column in the phenotype/*.fam file (optional; default 1)"<<endl;	
 		cout<<" -pace     [num]          "<<" specify terminal display update pace (default 100000 SNPs or 100000 iterations)."<<endl;
+		cout<<" -outdir   [path]         "<<" specify output directory path (default \"./output/\")"<<endl; 
 		cout<<" -o        [prefix]       "<<" specify output file prefix (default \"result\")"<<endl;  
 		cout<<"          output: prefix.cXX.txt or prefix.sXX.txt from kinship/relatedness matrix estimation"<<endl;	
 		cout<<"          output: prefix.assoc.txt and prefix.log.txt form association tests"<<endl;	
@@ -519,6 +520,13 @@ void GEMMA::Assign(int argc, char ** argv, PARAM &cPar)
 			str.assign(argv[i]);
 			cPar.d_pace=atoi(str.c_str());
 		}
+		else if (strcmp(argv[i], "-outdir")==0) {
+			if(argv[i+1] == NULL || argv[i+1][0] == '-') {continue;}
+			++i;
+			str.clear();
+			str.assign(argv[i]);
+			cPar.path_out=str;
+		}
 		else if (strcmp(argv[i], "-o")==0) {
 			if(argv[i+1] == NULL || argv[i+1][0] == '-') {continue;}
 			++i;
@@ -1567,7 +1575,7 @@ void GEMMA::BatchRun (PARAM &cPar)
 void GEMMA::WriteLog (int argc, char ** argv, PARAM &cPar) 
 {
 	string file_str;
-	file_str="./output/"+cPar.file_out;
+	file_str=cPar.path_out+"/"+cPar.file_out;
 	file_str+=".log.txt";
 	
 	ofstream outfile (file_str.c_str(), ofstream::out);
diff --git a/lm.cpp b/lm.cpp
index c983253..7577d0a 100644
--- a/lm.cpp
+++ b/lm.cpp
@@ -65,6 +65,7 @@ void LM::CopyFromParam (PARAM &cPar)
 	file_bfile=cPar.file_bfile;
 	file_geno=cPar.file_geno;
 	file_out=cPar.file_out;
+	path_out=cPar.path_out;
 	file_gene=cPar.file_gene;
 	
 	time_opt=0.0;
@@ -100,7 +101,7 @@ void LM::CopyToParam (PARAM &cPar)
 void LM::WriteFiles () 
 {
 	string file_str;
-	file_str="./output/"+file_out;
+	file_str=path_out+"/"+file_out;
 	file_str+=".assoc.txt";
 
 	ofstream outfile (file_str.c_str(), ofstream::out);
diff --git a/lm.h b/lm.h
index 84a0322..ceec060 100644
--- a/lm.h
+++ b/lm.h
@@ -44,6 +44,7 @@ public:
 	string file_bfile;
 	string file_geno;
 	string file_out;
+	string path_out;
 	
 	string file_gene;
 	
diff --git a/lmm.cpp b/lmm.cpp
index fed94ee..e0b4160 100644
--- a/lmm.cpp
+++ b/lmm.cpp
@@ -1,5 +1,5 @@
 /*
-    Genome-wide Efficient Mixed Model Association (GEMMA)
+	Genome-wide Efficient Mixed Model Association (GEMMA)
     Copyright (C) 2011  Xiang Zhou
 
     This program is free software: you can redistribute it and/or modify
@@ -66,6 +66,7 @@ void LMM::CopyFromParam (PARAM &cPar)
 	file_bfile=cPar.file_bfile;
 	file_geno=cPar.file_geno;
 	file_out=cPar.file_out;
+	path_out=cPar.path_out;
 	file_gene=cPar.file_gene;
 	
 	l_min=cPar.l_min;
@@ -109,7 +110,7 @@ void LMM::CopyToParam (PARAM &cPar)
 void LMM::WriteFiles () 
 {
 	string file_str;
-	file_str="./output/"+file_out;
+	file_str=path_out+"/"+file_out;
 	file_str+=".assoc.txt";
 
 	ofstream outfile (file_str.c_str(), ofstream::out);
diff --git a/lmm.h b/lmm.h
index d65b785..45f9b72 100644
--- a/lmm.h
+++ b/lmm.h
@@ -61,6 +61,7 @@ public:
 	string file_bfile;
 	string file_geno;
 	string file_out;
+	string path_out;
 	
 	string file_gene;
 	
diff --git a/main.cpp b/main.cpp
index 9ab98ea..e1fb336 100644
--- a/main.cpp
+++ b/main.cpp
@@ -56,14 +56,14 @@ int main(int argc, char * argv[])
 	if (argc==2 && argv[1][0] == '-' && argv[1][1] == 'l') {
 		cGemma.PrintLicense();
 		return EXIT_SUCCESS;
-	}
-	
-	ifstream check_dir("output/");
-	if (!check_dir) {
-		mkdir("output", S_IRWXU|S_IRGRP|S_IROTH);
 	}	
 	
 	cGemma.Assign(argc, argv, cPar); 
+
+	ifstream check_dir((cPar.path_out).c_str());
+	if (!check_dir) {
+	  mkdir((cPar.path_out).c_str(), S_IRWXU|S_IRGRP|S_IROTH);
+	}	
 		
 	if (cPar.error==true) {return EXIT_FAILURE;}
 	     
diff --git a/mvlmm.cpp b/mvlmm.cpp
index 56540d8..4b910ee 100644
--- a/mvlmm.cpp
+++ b/mvlmm.cpp
@@ -68,6 +68,7 @@ void MVLMM::CopyFromParam (PARAM &cPar)
 	file_bfile=cPar.file_bfile;
 	file_geno=cPar.file_geno;
 	file_out=cPar.file_out;
+	path_out=cPar.path_out;
 	
 	l_min=cPar.l_min;
 	l_max=cPar.l_max;
@@ -132,7 +133,7 @@ void MVLMM::CopyToParam (PARAM &cPar)
 void MVLMM::WriteFiles () 
 {
 	string file_str;
-	file_str="./output/"+file_out;
+	file_str=path_out+"/"+file_out;
 	file_str+=".assoc.txt";
 	
 	ofstream outfile (file_str.c_str(), ofstream::out);
diff --git a/mvlmm.h b/mvlmm.h
index 5aa9bbf..129879c 100644
--- a/mvlmm.h
+++ b/mvlmm.h
@@ -47,6 +47,7 @@ public:
 	string file_bfile;
 	string file_geno;
 	string file_out;
+	string path_out;
 	
 	// MVLMM related parameters
 	double l_min;
diff --git a/param.cpp b/param.cpp
index edacc42..7a89ff8 100644
--- a/param.cpp
+++ b/param.cpp
@@ -41,7 +41,7 @@ using namespace std;
 
 PARAM::PARAM(void):	
 mode_silence (false), a_mode (0), k_mode(1), d_pace (100000),
-file_out("result"), 
+file_out("result"), path_out("./output/"),
 miss_level(0.05), maf_level(0.01), hwe_level(0), r2_level(0.9999),
 l_min(1e-5), l_max(1e5), n_region(10),p_nr(0.001),em_prec(0.0001),nr_prec(0.0001),em_iter(10000),nr_iter(100),crt(0),
 pheno_mean(0),
@@ -614,7 +614,7 @@ void PARAM::CalcKin (gsl_matrix *matrix_kin)  {
 void PARAM::WriteMatrix (const gsl_matrix *matrix_U, const string suffix) 
 {
 	string file_str;
-	file_str="./output/"+file_out;
+	file_str=path_out+"/"+file_out;
 	file_str+=".";
 	file_str+=suffix;
 	file_str+=".txt";	
@@ -640,7 +640,7 @@ void PARAM::WriteMatrix (const gsl_matrix *matrix_U, const string suffix)
 void PARAM::WriteVector (const gsl_vector *vector_D, const string suffix) 
 {
 	string file_str;
-	file_str="./output/"+file_out;
+	file_str=path_out+"/"+file_out;
 	file_str+=".";
 	file_str+=suffix;
 	file_str+=".txt";
diff --git a/param.h b/param.h
index 5d58c4c..fa18181 100644
--- a/param.h
+++ b/param.h
@@ -100,6 +100,7 @@ public:
 	string file_ku, file_kd;
 	string file_mk;
 	string file_out;
+	string path_out;
 	
 	string file_epm;		//estimated parameter file
 	string file_ebv;		//estimated breeding value file
diff --git a/prdt.cpp b/prdt.cpp
index 7570d36..2875119 100644
--- a/prdt.cpp
+++ b/prdt.cpp
@@ -61,6 +61,7 @@ void PRDT::CopyFromParam (PARAM &cPar)
 	file_bfile=cPar.file_bfile;
 	file_geno=cPar.file_geno;
 	file_out=cPar.file_out;
+	path_out=cPar.path_out;
 	
 	indicator_pheno=cPar.indicator_pheno;	
 	indicator_cvt=cPar.indicator_cvt;
@@ -94,7 +95,7 @@ void PRDT::CopyToParam (PARAM &cPar)
 void PRDT::WriteFiles (gsl_vector *y_prdt) 
 {
 	string file_str;
-	file_str="./output/"+file_out;
+	file_str=path_out+"/"+file_out;
 	file_str+=".";
 	file_str+="prdt";
 	file_str+=".txt";
@@ -121,7 +122,7 @@ void PRDT::WriteFiles (gsl_vector *y_prdt)
 void PRDT::WriteFiles (gsl_matrix *Y_full) 
 {
 	string file_str;
-	file_str="./output/"+file_out;
+	file_str=path_out+"/"+file_out;
 	file_str+=".prdt.txt";
 	
 	ofstream outfile (file_str.c_str(), ofstream::out);
diff --git a/prdt.h b/prdt.h
index 69043df..8af2cee 100644
--- a/prdt.h
+++ b/prdt.h
@@ -44,6 +44,7 @@ public:
 	string file_bfile;
 	string file_geno;
 	string file_out;
+	string path_out;
 	
 	vector<vector<int> > indicator_pheno;
 	vector<int> indicator_cvt;
diff --git a/vc.h b/vc.h
index aed5247..f34d72e 100644
--- a/vc.h
+++ b/vc.h
@@ -57,6 +57,7 @@ class VC {
 public:
 	// IO related parameters
 	string file_out;
+	string path_out;
 
 	vector<double> v_sigma2;
 	vector<double> v_se_sigma2;
@@ -72,6 +73,7 @@ public:
 	// Main functions
 	void CopyFromParam (PARAM &cPar);
 	void CopyToParam (PARAM &cPar);
+	void CalcVChe (const gsl_matrix *K, const gsl_matrix *W, const gsl_vector *y);
 	void CalcVCreml (const gsl_matrix *K, const gsl_matrix *W, const gsl_vector *y);
 };