Removing FORCE_FLOAT from some C++ source files.

5 files changed, 677 insertions, 462 deletions

diff --git a/src/eigenlib.h b/src/eigenlib.h
index e9768c2..493907c 100644
--- a/src/eigenlib.h
+++ b/src/eigenlib.h
@@ -1,6 +1,6 @@
 /*
-	Genome-wide Efficient Mixed Model Association (GEMMA)
-    Copyright (C) 2011  Xiang Zhou
+    Genome-wide Efficient Mixed Model Association (GEMMA)
+    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 <http://www.gnu.org/licenses/>.
+    along with this program. If not, see <http://www.gnu.org/licenses/>.
 */
 
 #ifndef __EIGENLIB_H__                
@@ -23,12 +23,10 @@
 
 using namespace std;
 
-
 void eigenlib_dgemm (const char *TransA, const char *TransB, const double alpha, const gsl_matrix *A, const gsl_matrix *B, const double beta, gsl_matrix *C);
 void eigenlib_dgemv (const char *TransA, const double alpha, const gsl_matrix *A, const gsl_vector *x, const double beta, gsl_vector *y);
 void eigenlib_invert(gsl_matrix *A);
 void eigenlib_dsyr (const double alpha, const gsl_vector *b, gsl_matrix *A);
 void eigenlib_eigensymm (const gsl_matrix *G, gsl_matrix *U, gsl_vector *eval);
 
-
 #endif
diff --git a/src/lapack.cpp b/src/lapack.cpp
index d9e4149..18399a2 100644
--- a/src/lapack.cpp
+++ b/src/lapack.cpp
@@ -1,6 +1,6 @@
 /*
-	Genome-wide Efficient Mixed Model Association (GEMMA)
-    Copyright (C) 2011  Xiang Zhou
+    Genome-wide Efficient Mixed Model Association (GEMMA)
+    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 <http://www.gnu.org/licenses/>.
+    along with this program. If not, see <http://www.gnu.org/licenses/>.
 */
 
 #include <iostream>
@@ -25,82 +25,129 @@
 
 using namespace std;
 
-extern "C" void sgemm_(char *TRANSA, char *TRANSB, int *M, int *N, int *K, float *ALPHA, float *A, int *LDA, float *B, int *LDB, float *BETA, float *C, int *LDC);
+extern "C" void sgemm_(char *TRANSA, char *TRANSB, int *M, int *N, int *K,
+		       float *ALPHA, float *A, int *LDA, float *B, int *LDB,
+		       float *BETA, float *C, int *LDC);
 extern "C" void spotrf_(char *UPLO, int *N, float *A, int *LDA, int *INFO);
-extern "C" void spotrs_(char *UPLO, int *N, int *NRHS, float *A, int *LDA, float *B, int *LDB, int *INFO);
-extern "C" void ssyev_(char* JOBZ, char* UPLO, int *N, float *A, int *LDA, float *W, float *WORK, int *LWORK, int *INFO);
-extern "C" void ssyevr_(char* JOBZ, char *RANGE, char* UPLO, int *N, float *A, int *LDA, float *VL, float *VU, int *IL, int *IU, float *ABSTOL, int *M, float *W, float *Z, int *LDZ, int *ISUPPZ, float *WORK, int *LWORK, int *IWORK, int *LIWORK, int *INFO);
+extern "C" void spotrs_(char *UPLO, int *N, int *NRHS, float *A, int *LDA,
+			float *B, int *LDB, int *INFO);
+extern "C" void ssyev_(char* JOBZ, char* UPLO, int *N, float *A, int *LDA,
+		       float *W, float *WORK, int *LWORK, int *INFO);
+extern "C" void ssyevr_(char* JOBZ, char *RANGE, char* UPLO, int *N,
+			float *A, int *LDA, float *VL, float *VU, int *IL,
+			int *IU, float *ABSTOL, int *M, float *W, float *Z,
+			int *LDZ, int *ISUPPZ, float *WORK, int *LWORK,
+			int *IWORK, int *LIWORK, int *INFO);
 extern "C" double sdot_(int *N, float *DX, int *INCX, float *DY, int *INCY);
 
-extern "C" void dgemm_(char *TRANSA, char *TRANSB, int *M, int *N, int *K, double *ALPHA, double *A, int *LDA, double *B, int *LDB, double *BETA, double *C, int *LDC);
+extern "C" void dgemm_(char *TRANSA, char *TRANSB, int *M, int *N, int *K,
+		       double *ALPHA, double *A, int *LDA, double *B,
+		       int *LDB, double *BETA, double *C, int *LDC);
 extern "C" void dpotrf_(char *UPLO, int *N, double *A, int *LDA, int *INFO);
-extern "C" void dpotrs_(char *UPLO, int *N, int *NRHS, double *A, int *LDA, double *B, int *LDB, int *INFO);
-extern "C" void dsyev_(char* JOBZ, char* UPLO, int *N, double *A, int *LDA, double *W, double *WORK, int *LWORK, int *INFO);
-extern "C" void dsyevr_(char* JOBZ, char *RANGE, char* UPLO, int *N, double *A, int *LDA, double *VL, double *VU, int *IL, int *IU, double *ABSTOL, int *M, double *W, double *Z, int *LDZ, int *ISUPPZ, double *WORK, int *LWORK, int *IWORK, int *LIWORK, int *INFO);
+extern "C" void dpotrs_(char *UPLO, int *N, int *NRHS, double *A, int *LDA,
+			double *B, int *LDB, int *INFO);
+extern "C" void dsyev_(char* JOBZ, char* UPLO, int *N, double *A, int *LDA,
+		       double *W, double *WORK, int *LWORK, int *INFO);
+extern "C" void dsyevr_(char* JOBZ, char *RANGE, char* UPLO, int *N,
+			double *A, int *LDA, double *VL, double *VU,
+			int *IL, int *IU, double *ABSTOL, int *M,
+			double *W, double *Z, int *LDZ, int *ISUPPZ,
+			double *WORK, int *LWORK, int *IWORK,
+			int *LIWORK, int *INFO);
 extern "C" double ddot_(int *N, double *DX, int *INCX, double *DY, int *INCY);
 
-
-//cholesky decomposition, A is distroyed
-void lapack_float_cholesky_decomp (gsl_matrix_float *A)
-{
+// Cholesky decomposition, A is destroyed.
+void lapack_float_cholesky_decomp (gsl_matrix_float *A) {
 	int N=A->size1, LDA=A->size1, INFO;
 	char UPLO='L';
 	
-	if (N!=(int)A->size2) {cout<<"Matrix needs to be symmetric and same dimension in lapack_cholesky_decomp."<<endl; return;}
+	if (N!=(int)A->size2) {
+	  cout << "Matrix needs to be symmetric and same dimension in" <<
+	    "lapack_cholesky_decomp." << endl;
+	  return;
+	}
 	
 	spotrf_(&UPLO, &N, A->data, &LDA, &INFO);
-	if (INFO!=0) {cout<<"Cholesky decomposition unsuccessful in lapack_cholesky_decomp."<<endl; return;}	
+	if (INFO!=0) {
+	  cout << "Cholesky decomposition unsuccessful in" <<
+	    "lapack_cholesky_decomp." << endl;
+	  return;
+	}	
 	
 	return;
 }
 
-//cholesky decomposition, A is distroyed
-void lapack_cholesky_decomp (gsl_matrix *A)
-{
+// Cholesky decomposition, A is destroyed.
+void lapack_cholesky_decomp (gsl_matrix *A) {
 	int N=A->size1, LDA=A->size1, INFO;
 	char UPLO='L';
 	
-	if (N!=(int)A->size2) {cout<<"Matrix needs to be symmetric and same dimension in lapack_cholesky_decomp."<<endl; return;}
+	if (N!=(int)A->size2) {
+	  cout << "Matrix needs to be symmetric and same dimension in" <<
+	    "lapack_cholesky_decomp." << endl;
+	  return;
+	}
 	
 	dpotrf_(&UPLO, &N, A->data, &LDA, &INFO);
-	if (INFO!=0) {cout<<"Cholesky decomposition unsuccessful in lapack_cholesky_decomp."<<endl; return;}	
-	
+	if (INFO!=0) {
+	  cout << "Cholesky decomposition unsuccessful in" <<
+	    "lapack_cholesky_decomp."<<endl;
+	  return;
+	} 
+
 	return;
 }
 
-//cholesky solve, A is decomposed, 
-void lapack_float_cholesky_solve (gsl_matrix_float *A, const gsl_vector_float *b, gsl_vector_float *x)
-{
+// Cholesky solve, A is decomposed.
+void lapack_float_cholesky_solve (gsl_matrix_float *A,
+				  const gsl_vector_float *b,
+				  gsl_vector_float *x) {
 	int N=A->size1, NRHS=1, LDA=A->size1, LDB=b->size, INFO;
 	char UPLO='L';
 	
-	if (N!=(int)A->size2 || N!=LDB) {cout<<"Matrix needs to be symmetric and same dimension in lapack_cholesky_solve."<<endl; return;}
+	if (N!=(int)A->size2 || N!=LDB) {
+	  cout << "Matrix needs to be symmetric and same dimension in" <<
+	    "lapack_cholesky_solve." << endl;
+	  return;
+	}
 	
 	gsl_vector_float_memcpy (x, b);
 	spotrs_(&UPLO, &N, &NRHS, A->data, &LDA, x->data, &LDB, &INFO);
-	if (INFO!=0) {cout<<"Cholesky solve unsuccessful in lapack_cholesky_solve."<<endl; return;}	
+	if (INFO!=0) {
+	  cout << "Cholesky solve unsuccessful in lapack_cholesky_solve." <<
+	    endl;
+	  return;
+	}
 	
 	return;
 }
 
-//cholesky solve, A is decomposed, 
-void lapack_cholesky_solve (gsl_matrix *A, const gsl_vector *b, gsl_vector *x)
-{
+// Cholesky solve, A is decomposed.
+void lapack_cholesky_solve (gsl_matrix *A, const gsl_vector *b,
+			    gsl_vector *x) {
 	int N=A->size1, NRHS=1, LDA=A->size1, LDB=b->size, INFO;
 	char UPLO='L';
 	
-	if (N!=(int)A->size2 || N!=LDB) {cout<<"Matrix needs to be symmetric and same dimension in lapack_cholesky_solve."<<endl; return;}
+	if (N!=(int)A->size2 || N!=LDB) {
+	  cout << "Matrix needs to be symmetric and same dimension in" <<
+	    "lapack_cholesky_solve." << endl;
+	  return;
+	}
 	
 	gsl_vector_memcpy (x, b);
 	dpotrs_(&UPLO, &N, &NRHS, A->data, &LDA, x->data, &LDB, &INFO);
-	if (INFO!=0) {cout<<"Cholesky solve unsuccessful in lapack_cholesky_solve."<<endl; return;}	
+	if (INFO!=0) {
+	  cout << "Cholesky solve unsuccessful in lapack_cholesky_solve." <<
+	    endl;
+	  return;
+	}
 	
 	return;
 }
 
-
-void lapack_sgemm (char *TransA, char *TransB, float alpha, const gsl_matrix_float *A, const gsl_matrix_float *B, float beta, gsl_matrix_float *C)
-{
+void lapack_sgemm (char *TransA, char *TransB, float alpha,
+		   const gsl_matrix_float *A, const gsl_matrix_float *B,
+		   float beta, gsl_matrix_float *C) {
 	int M, N, K1, K2, LDA=A->size1, LDB=B->size1, LDC=C->size2;
 	
 	if (*TransA=='N' || *TransA=='n') {M=A->size1; K1=A->size2;}
@@ -111,8 +158,14 @@ void lapack_sgemm (char *TransA, char *TransB, float alpha, const gsl_matrix_flo
 	else if (*TransB=='T' || *TransB=='t')  {N=B->size1; K2=B->size2;}
 	else {cout<<"need 'N' or 'T' in lapack_sgemm"<<endl;  return;}
 	
-	if (K1!=K2) {cout<<"A and B not compatible in lapack_sgemm"<<endl; return;}
-	if (C->size1!=(size_t)M || C->size2!=(size_t)N) {cout<<"C not compatible in lapack_sgemm"<<endl; return;}
+	if (K1!=K2) {
+	  cout<<"A and B not compatible in lapack_sgemm"<<endl;
+	  return;
+	}
+	if (C->size1!=(size_t)M || C->size2!=(size_t)N) {
+	  cout<<"C not compatible in lapack_sgemm"<<endl;
+	  return;
+	}
 	
 	gsl_matrix_float *A_t=gsl_matrix_float_alloc (A->size2, A->size1);
 	gsl_matrix_float_transpose_memcpy (A_t, A);
@@ -121,7 +174,8 @@ void lapack_sgemm (char *TransA, char *TransB, float alpha, const gsl_matrix_flo
 	gsl_matrix_float *C_t=gsl_matrix_float_alloc (C->size2, C->size1);
 	gsl_matrix_float_transpose_memcpy (C_t, C);
 	
-	sgemm_(TransA, TransB, &M, &N, &K1, &alpha, A_t->data, &LDA, B_t->data, &LDB, &beta, C_t->data, &LDC);
+	sgemm_(TransA, TransB, &M, &N, &K1, &alpha, A_t->data, &LDA,
+	       B_t->data, &LDB, &beta, C_t->data, &LDC);
 	gsl_matrix_float_transpose_memcpy (C, C_t);
 	
 	gsl_matrix_float_free (A_t);
@@ -132,8 +186,9 @@ void lapack_sgemm (char *TransA, char *TransB, float alpha, const gsl_matrix_flo
 
 
 
-void lapack_dgemm (char *TransA, char *TransB, double alpha, const gsl_matrix *A, const gsl_matrix *B, double beta, gsl_matrix *C)
-{
+void lapack_dgemm (char *TransA, char *TransB, double alpha,
+		   const gsl_matrix *A, const gsl_matrix *B,
+		   double beta, gsl_matrix *C) {
 	int M, N, K1, K2, LDA=A->size1, LDB=B->size1, LDC=C->size2;
 	
 	if (*TransA=='N' || *TransA=='n') {M=A->size1; K1=A->size2;}
@@ -144,8 +199,14 @@ void lapack_dgemm (char *TransA, char *TransB, double alpha, const gsl_matrix *A
 	else if (*TransB=='T' || *TransB=='t')  {N=B->size1; K2=B->size2;}
 	else {cout<<"need 'N' or 'T' in lapack_dgemm"<<endl;  return;}
 	
-	if (K1!=K2) {cout<<"A and B not compatible in lapack_dgemm"<<endl; return;}
-	if (C->size1!=(size_t)M || C->size2!=(size_t)N) {cout<<"C not compatible in lapack_dgemm"<<endl; return;}
+	if (K1!=K2) {
+	  cout << "A and B not compatible in lapack_dgemm"<<endl;
+	  return;
+	}
+	if (C->size1!=(size_t)M || C->size2!=(size_t)N) {
+	  cout<<"C not compatible in lapack_dgemm"<<endl;
+	  return;
+	}
 	
 	gsl_matrix *A_t=gsl_matrix_alloc (A->size2, A->size1);
 	gsl_matrix_transpose_memcpy (A_t, A);
@@ -154,7 +215,8 @@ void lapack_dgemm (char *TransA, char *TransB, double alpha, const gsl_matrix *A
 	gsl_matrix *C_t=gsl_matrix_alloc (C->size2, C->size1);
 	gsl_matrix_transpose_memcpy (C_t, C);
 
-	dgemm_(TransA, TransB, &M, &N, &K1, &alpha, A_t->data, &LDA, B_t->data, &LDB, &beta, C_t->data, &LDC);
+	dgemm_(TransA, TransB, &M, &N, &K1, &alpha, A_t->data, &LDA,
+	       B_t->data, &LDB, &beta, C_t->data, &LDC);
 
 	gsl_matrix_transpose_memcpy (C, C_t);
 	
@@ -164,58 +226,79 @@ void lapack_dgemm (char *TransA, char *TransB, double alpha, const gsl_matrix *A
 	return;
 }
 
-
-
-//eigen value decomposition, matrix A is destroyed, float seems to have problem with large matrices (in mac)
-void lapack_float_eigen_symmv (gsl_matrix_float *A, gsl_vector_float *eval, gsl_matrix_float *evec, const size_t flag_largematrix)
-{
+// Eigen value decomposition, matrix A is destroyed, float seems to
+// have problem with large matrices (in mac).
+void lapack_float_eigen_symmv (gsl_matrix_float *A, gsl_vector_float *eval,
+			       gsl_matrix_float *evec,
+			       const size_t flag_largematrix) {
 	if (flag_largematrix==1) {
 		int N=A->size1, LDA=A->size1, INFO, LWORK=-1;
 		char JOBZ='V', UPLO='L';
 				
-		if (N!=(int)A->size2 || N!=(int)eval->size) {cout<<"Matrix needs to be symmetric and same dimension in lapack_eigen_symmv."<<endl; return;}
-		
-		//	float temp[1];
-		//	ssyev_(&JOBZ, &UPLO, &N, A->data, &LDA, eval->data, temp, &LWORK, &INFO);
-		//	if (INFO!=0) {cout<<"Work space estimate unsuccessful in lapack_eigen_symmv."<<endl; return;}
-		//	LWORK=(int)temp[0];
+		if (N!=(int)A->size2 || N!=(int)eval->size) {
+		  cout << "Matrix needs to be symmetric and same" <<
+		    "dimension in lapack_eigen_symmv."<<endl;
+		  return;
+		}
 		
 		LWORK=3*N;
 		float *WORK=new float [LWORK];	
-		ssyev_(&JOBZ, &UPLO, &N, A->data, &LDA, eval->data, WORK, &LWORK, &INFO);
-		if (INFO!=0) {cout<<"Eigen decomposition unsuccessful in lapack_eigen_symmv."<<endl; return;}
+		ssyev_(&JOBZ, &UPLO, &N, A->data, &LDA, eval->data, WORK,
+		       &LWORK, &INFO);
+		if (INFO!=0) {
+		  cout << "Eigen decomposition unsuccessful in" <<
+		    "lapack_eigen_symmv."<<endl;
+		  return;
+		}
 		
-		gsl_matrix_float_view A_sub=gsl_matrix_float_submatrix(A, 0, 0, N, N);
+		gsl_matrix_float_view A_sub =
+		  gsl_matrix_float_submatrix(A, 0, 0, N, N);
 		gsl_matrix_float_memcpy (evec, &A_sub.matrix);
 		gsl_matrix_float_transpose (evec);
 		
 		delete [] WORK;
 	} else {	
-		int N=A->size1, LDA=A->size1, LDZ=A->size1, INFO, LWORK=-1, LIWORK=-1;
+		int N=A->size1, LDA=A->size1, LDZ=A->size1, INFO,
+		  LWORK=-1, LIWORK=-1;
 		char JOBZ='V', UPLO='L', RANGE='A';
 		float ABSTOL=1.0E-7;
 		
-		//VL, VU, IL, IU are not referenced; M equals N if RANGE='A'
+		// VL, VU, IL, IU are not referenced; M equals N if RANGE='A'.
 		float VL=0.0, VU=0.0;
 		int IL=0, IU=0, M;
 		
-		if (N!=(int)A->size2 || N!=(int)eval->size) {cout<<"Matrix needs to be symmetric and same dimension in lapack_float_eigen_symmv."<<endl; return;}
+		if (N!=(int)A->size2 || N!=(int)eval->size) {
+		  cout << "Matrix needs to be symmetric and same" <<
+		    "dimension in lapack_float_eigen_symmv." << endl;
+		  return;
+		}
 		
 		int *ISUPPZ=new int [2*N];
 				
 		float WORK_temp[1];
 		int IWORK_temp[1];
-		ssyevr_(&JOBZ, &RANGE, &UPLO, &N, A->data, &LDA, &VL, &VU, &IL, &IU, &ABSTOL, &M, eval->data, evec->data, &LDZ, ISUPPZ, WORK_temp, &LWORK, IWORK_temp, &LIWORK, &INFO);
-		if (INFO!=0) {cout<<"Work space estimate unsuccessful in lapack_float_eigen_symmv."<<endl; return;}	
+		ssyevr_(&JOBZ, &RANGE, &UPLO, &N, A->data, &LDA, &VL,
+			&VU, &IL, &IU, &ABSTOL, &M, eval->data,
+			evec->data, &LDZ, ISUPPZ, WORK_temp, &LWORK,
+			IWORK_temp, &LIWORK, &INFO);
+		if (INFO!=0) {
+		  cout << "Work space estimate unsuccessful in" <<
+		    "lapack_float_eigen_symmv." << endl;
+		  return;
+		}
 		LWORK=(int)WORK_temp[0]; LIWORK=(int)IWORK_temp[0];	
 		 
-		//LWORK=26*N;
-		//LIWORK=10*N;
 		float *WORK=new float [LWORK];
 		int *IWORK=new int [LIWORK];
 		
-		ssyevr_(&JOBZ, &RANGE, &UPLO, &N, A->data, &LDA, &VL, &VU, &IL, &IU, &ABSTOL, &M, eval->data, evec->data, &LDZ, ISUPPZ, WORK, &LWORK, IWORK, &LIWORK, &INFO);
-		if (INFO!=0) {cout<<"Eigen decomposition unsuccessful in lapack_float_eigen_symmv."<<endl; return;}
+		ssyevr_(&JOBZ, &RANGE, &UPLO, &N, A->data, &LDA, &VL,
+			&VU, &IL, &IU, &ABSTOL, &M, eval->data, evec->data,
+			&LDZ, ISUPPZ, WORK, &LWORK, IWORK, &LIWORK, &INFO);
+		if (INFO!=0) {
+		  cout << "Eigen decomposition unsuccessful in" <<
+		    "lapack_float_eigen_symmv." << endl;
+		  return;
+		}
 		
 		gsl_matrix_float_transpose (evec);
 		
@@ -230,24 +313,28 @@ void lapack_float_eigen_symmv (gsl_matrix_float *A, gsl_vector_float *eval, gsl_
 
 
 
-//eigen value decomposition, matrix A is destroyed
-void lapack_eigen_symmv (gsl_matrix *A, gsl_vector *eval, gsl_matrix *evec, const size_t flag_largematrix)
-{
+// Eigenvalue decomposition, matrix A is destroyed.
+void lapack_eigen_symmv (gsl_matrix *A, gsl_vector *eval, gsl_matrix *evec,
+			 const size_t flag_largematrix) {
 	if (flag_largematrix==1) {
 		int N=A->size1, LDA=A->size1, INFO, LWORK=-1;
 		char JOBZ='V', UPLO='L';		
 		
-		if (N!=(int)A->size2 || N!=(int)eval->size) {cout<<"Matrix needs to be symmetric and same dimension in lapack_eigen_symmv."<<endl; return;}
-		
-		//	double temp[1];
-		//	dsyev_(&JOBZ, &UPLO, &N, A->data, &LDA, eval->data, temp, &LWORK, &INFO);
-		//	if (INFO!=0) {cout<<"Work space estimate unsuccessful in lapack_eigen_symmv."<<endl; return;}		
-		//	LWORK=(int)temp[0];
+		if (N!=(int)A->size2 || N!=(int)eval->size) {
+		  cout << "Matrix needs to be symmetric and same" <<
+		    "dimension in lapack_eigen_symmv." << endl;
+		  return;
+		}
 		
 		LWORK=3*N;
 		double *WORK=new double [LWORK];	
-		dsyev_(&JOBZ, &UPLO, &N, A->data, &LDA, eval->data, WORK, &LWORK, &INFO);
-		if (INFO!=0) {cout<<"Eigen decomposition unsuccessful in lapack_eigen_symmv."<<endl; return;}
+		dsyev_(&JOBZ, &UPLO, &N, A->data, &LDA, eval->data, WORK,
+		       &LWORK, &INFO);
+		if (INFO!=0) {
+		  cout<<"Eigen decomposition unsuccessful in" <<
+		    "lapack_eigen_symmv." << endl;
+		  return;
+		}
 		
 		gsl_matrix_view A_sub=gsl_matrix_submatrix(A, 0, 0, N, N);
 		gsl_matrix_memcpy (evec, &A_sub.matrix);
@@ -255,32 +342,48 @@ void lapack_eigen_symmv (gsl_matrix *A, gsl_vector *eval, gsl_matrix *evec, cons
 		
 		delete [] WORK;
 	} else {	
-		int N=A->size1, LDA=A->size1, LDZ=A->size1, INFO, LWORK=-1, LIWORK=-1;
+   	        int N=A->size1, LDA=A->size1, LDZ=A->size1, INFO;
+		int LWORK=-1, LIWORK=-1;
 		char JOBZ='V', UPLO='L', RANGE='A';
 		double ABSTOL=1.0E-7;
 		
-		//VL, VU, IL, IU are not referenced; M equals N if RANGE='A'
+		// VL, VU, IL, IU are not referenced; M equals N if RANGE='A'.
 		double VL=0.0, VU=0.0;
 		int IL=0, IU=0, M;
 		
-		if (N!=(int)A->size2 || N!=(int)eval->size) {cout<<"Matrix needs to be symmetric and same dimension in lapack_eigen_symmv."<<endl; return;}
+		if (N!=(int)A->size2 || N!=(int)eval->size) {
+		  cout << "Matrix needs to be symmetric and same" <<
+		    "dimension in lapack_eigen_symmv." << endl;
+		  return;
+		}
 		
 		int *ISUPPZ=new int [2*N];
 		
 		double WORK_temp[1];
 		int IWORK_temp[1];
 
-		dsyevr_(&JOBZ, &RANGE, &UPLO, &N, A->data, &LDA, &VL, &VU, &IL, &IU, &ABSTOL, &M, eval->data, evec->data, &LDZ, ISUPPZ, WORK_temp, &LWORK, IWORK_temp, &LIWORK, &INFO);
-		if (INFO!=0) {cout<<"Work space estimate unsuccessful in lapack_eigen_symmv."<<endl; return;}	
+		dsyevr_(&JOBZ, &RANGE, &UPLO, &N, A->data, &LDA, &VL, &VU,
+			&IL, &IU, &ABSTOL, &M, eval->data, evec->data,
+			&LDZ, ISUPPZ, WORK_temp, &LWORK, IWORK_temp,
+			&LIWORK, &INFO);
+		if (INFO!=0) {
+		  cout << "Work space estimate unsuccessful in" <<
+		    "lapack_eigen_symmv." << endl;
+		  return;
+		}	
 		LWORK=(int)WORK_temp[0]; LIWORK=(int)IWORK_temp[0];	
 
-		//LWORK=26*N;
-		//LIWORK=10*N;
 		double *WORK=new double [LWORK];
 		int *IWORK=new int [LIWORK];
 		
-		dsyevr_(&JOBZ, &RANGE, &UPLO, &N, A->data, &LDA, &VL, &VU, &IL, &IU, &ABSTOL, &M, eval->data, evec->data, &LDZ, ISUPPZ, WORK, &LWORK, IWORK, &LIWORK, &INFO);
-		if (INFO!=0) {cout<<"Eigen decomposition unsuccessful in lapack_eigen_symmv."<<endl; return;}
+		dsyevr_(&JOBZ, &RANGE, &UPLO, &N, A->data, &LDA, &VL, &VU,
+			&IL, &IU, &ABSTOL, &M, eval->data, evec->data,
+			&LDZ, ISUPPZ, WORK, &LWORK, IWORK, &LIWORK, &INFO);
+		if (INFO!=0) {
+		  cout << "Eigen decomposition unsuccessful in" <<
+		    "lapack_eigen_symmv." << endl;
+		  return;
+		}
 
 		gsl_matrix_transpose (evec);
 		
@@ -292,9 +395,9 @@ void lapack_eigen_symmv (gsl_matrix *A, gsl_vector *eval, gsl_matrix *evec, cons
 	return;
 }
 
-//DO NOT set eigen values to be positive
-double EigenDecomp (gsl_matrix *G, gsl_matrix *U, gsl_vector *eval, const size_t flag_largematrix)
-{
+// 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
@@ -302,15 +405,8 @@ double EigenDecomp (gsl_matrix *G, gsl_matrix *U, gsl_vector *eval, const size_t
 	gsl_eigen_symmv (G, eval, U, w);
 	gsl_eigen_symmv_free (w);	
 #endif	
-	/*
-	for (size_t i=0; i<eval->size; ++i) {
-		if (gsl_vector_get (eval, i)<1e-10) {
-//			cout<<gsl_vector_get (eval, i)<<endl;
-			gsl_vector_set (eval, i, 0);			
-		}
-	}
-	*/
-	//calculate track_G=mean(diag(G))	
+
+	// Calculate track_G=mean(diag(G)).
 	double d=0.0;
 	for (size_t i=0; i<eval->size; ++i) {
 		d+=gsl_vector_get(eval, i);
@@ -321,55 +417,54 @@ double EigenDecomp (gsl_matrix *G, gsl_matrix *U, gsl_vector *eval, const size_t
 }
 
 
-//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)
-{
+// 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 float precision eigen decomposition; plus, float precision eigen decomposition in lapack may not work on OS 10.4
-	//first change to double precision
+	
+	// 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; i<G->size1; i++) {
 		for (size_t j=0; j<G->size2; j++) {
-			gsl_matrix_set(G_double, i, j, gsl_matrix_float_get(G, i, 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
+	// Change back to float precision.
 	for (size_t i=0; i<G->size1; i++) {
 		for (size_t j=0; j<G->size2; j++) {
-			gsl_matrix_float_set(K, i, j, gsl_matrix_get(G_double, i, j));
+			gsl_matrix_float_set(K, i, j,
+					     gsl_matrix_get(G_double, i, j));
 		}
 	}
 	for (size_t i=0; i<U->size1; i++) {
 		for (size_t j=0; j<U->size2; j++) {
-			gsl_matrix_float_set(U, i, j, gsl_matrix_get(U_double, i, j));
+			gsl_matrix_float_set(U, i, j,
+					     gsl_matrix_get(U_double, i, j));
 		}
 	}
 	for (size_t i=0; i<eval->size; i++) {
 		gsl_vector_float_set(eval, i, gsl_vector_get(eval_double, i));
 	}	
 	
-	//delete double precision matrices
+	// Delete double-precision matrices.
 	gsl_matrix_free (G_double);
 	gsl_matrix_free (U_double);
 	gsl_vector_free (eval_double);
 #endif
-	/*
-	for (size_t i=0; i<eval->size; ++i) {
-		if (gsl_vector_float_get (eval, i)<1e-10) {
-			gsl_vector_float_set (eval, i, 0);
-		}
-	}
-	*/
-	//calculate track_G=mean(diag(G))	
-	double d=0.0;
+
+	// Calculate track_G=mean(diag(G)).
+	double d = 0.0;
 	for (size_t i=0; i<eval->size; ++i) {
 		d+=gsl_vector_float_get(eval, i);
 	}
@@ -379,8 +474,7 @@ double EigenDecomp (gsl_matrix_float *G, gsl_matrix_float *U, gsl_vector_float *
 }
 
 
-double CholeskySolve(gsl_matrix *Omega, gsl_vector *Xty, gsl_vector *OiXty)
-{
+double CholeskySolve(gsl_matrix *Omega, gsl_vector *Xty, gsl_vector *OiXty) {
 	double logdet_O=0.0;
 	
 #ifdef WITH_LAPACK
@@ -394,7 +488,6 @@ double CholeskySolve(gsl_matrix *Omega, gsl_vector *Xty, gsl_vector *OiXty)
 	int status = gsl_linalg_cholesky_decomp(Omega);
 	if(status == GSL_EDOM) {
 		cout << "## non-positive definite matrix" << endl; 
-		//		exit(0); 
 	}
 	
 	for (size_t i=0; i<Omega->size1; ++i) {
@@ -404,16 +497,15 @@ double CholeskySolve(gsl_matrix *Omega, gsl_vector *Xty, gsl_vector *OiXty)
 	
 	gsl_vector_memcpy (OiXty, Xty);
 	gsl_blas_dtrsv(CblasLower, CblasNoTrans, CblasNonUnit, Omega, OiXty); 
-	gsl_blas_dtrsv(CblasUpper, CblasNoTrans, CblasNonUnit, Omega, OiXty); 	
-	//	gsl_linalg_cholesky_solve(XtX, Xty, iXty);
+	gsl_blas_dtrsv(CblasUpper, CblasNoTrans, CblasNonUnit, Omega, OiXty);
 #endif
 	
 	return logdet_O;
 }
 
 
-double CholeskySolve(gsl_matrix_float *Omega, gsl_vector_float *Xty, gsl_vector_float *OiXty)
-{
+double CholeskySolve(gsl_matrix_float *Omega, gsl_vector_float *Xty,
+		     gsl_vector_float *OiXty) {
 	double logdet_O=0.0;
 	
 #ifdef WITH_LAPACK
@@ -436,7 +528,6 @@ double CholeskySolve(gsl_matrix_float *Omega, gsl_vector_float *Xty, gsl_vector_
 	int status = gsl_linalg_cholesky_decomp(Omega_double);
 	if(status == GSL_EDOM) {
 		cout << "## non-positive definite matrix" << endl; 
-		//		exit(0); 
 	}	
 	
 	for (size_t i=0; i<Omega->size1; ++i) {
@@ -450,8 +541,7 @@ double CholeskySolve(gsl_matrix_float *Omega, gsl_vector_float *Xty, gsl_vector_
 	
 	gsl_vector_float_memcpy (OiXty, Xty);
 	gsl_blas_strsv(CblasLower, CblasNoTrans, CblasNonUnit, Omega, OiXty); 
-	gsl_blas_strsv(CblasUpper, CblasNoTrans, CblasNonUnit, Omega, OiXty); 	
-	//	gsl_linalg_cholesky_solve(XtX, Xty, iXty);
+	gsl_blas_strsv(CblasUpper, CblasNoTrans, CblasNonUnit, Omega, OiXty);
 	
 	gsl_matrix_free (Omega_double);
 #endif
@@ -460,129 +550,124 @@ double CholeskySolve(gsl_matrix_float *Omega, gsl_vector_float *Xty, gsl_vector_
 }	
 
 
-//LU decomposition
-void LUDecomp (gsl_matrix *LU, gsl_permutation *p, int *signum)
-{
+// LU decomposition.
+void LUDecomp (gsl_matrix *LU, gsl_permutation *p, int *signum) {
 	gsl_linalg_LU_decomp (LU, p, signum);
 	return;
 }
 
-void LUDecomp (gsl_matrix_float *LU, gsl_permutation *p, int *signum)
-{
+void LUDecomp (gsl_matrix_float *LU, gsl_permutation *p, int *signum) {
 	gsl_matrix *LU_double=gsl_matrix_alloc (LU->size1, LU->size2);
 	
-	//copy float matrix to double	
+	// Copy float matrix to double.	
 	for (size_t i=0; i<LU->size1; i++) {
 		for (size_t j=0; j<LU->size2; j++) {
-			gsl_matrix_set (LU_double, i, j, gsl_matrix_float_get(LU, i, j));
+			gsl_matrix_set (LU_double, i, j,
+					gsl_matrix_float_get(LU, i, j));
 		}
 	}
 	
-	//LU decomposition
+	// LU decomposition.
 	gsl_linalg_LU_decomp (LU_double, p, signum);
 	
-	//copy float matrix to double
+	// Copy float matrix to double.
 	for (size_t i=0; i<LU->size1; i++) {
 		for (size_t j=0; j<LU->size2; j++) {
-			gsl_matrix_float_set (LU, i, j, gsl_matrix_get(LU_double, i, j));
+			gsl_matrix_float_set (LU, i, j,
+					      gsl_matrix_get(LU_double, i, j));
 		}
 	}
 	
-	//free matrix
+	// Free matrix.
 	gsl_matrix_free (LU_double);
 	return;
 }
 
 
-//LU invert
-void LUInvert (const gsl_matrix *LU, const gsl_permutation *p, gsl_matrix *inverse)
-{
+// LU invert.
+void LUInvert (const gsl_matrix *LU, const gsl_permutation *p,
+	       gsl_matrix *inverse) {
 	gsl_linalg_LU_invert (LU, p, inverse);
 	return;
 }
 
-void LUInvert (const gsl_matrix_float *LU, const gsl_permutation *p, gsl_matrix_float *inverse)
-{
+void LUInvert (const gsl_matrix_float *LU, const gsl_permutation *p,
+	       gsl_matrix_float *inverse) {
 	gsl_matrix *LU_double=gsl_matrix_alloc (LU->size1, LU->size2);
-	gsl_matrix *inverse_double=gsl_matrix_alloc (inverse->size1, inverse->size2);
+	gsl_matrix *inverse_double=gsl_matrix_alloc (inverse->size1,
+						     inverse->size2);
 	
-	//copy float matrix to double	
+	// Copy float matrix to double.	
 	for (size_t i=0; i<LU->size1; i++) {
 		for (size_t j=0; j<LU->size2; j++) {
-			gsl_matrix_set (LU_double, i, j, gsl_matrix_float_get(LU, i, j));
+			gsl_matrix_set (LU_double, i, j,
+					gsl_matrix_float_get(LU, i, j));
 		}
 	}
 	
-	//LU decomposition
+	// LU decomposition.
 	gsl_linalg_LU_invert (LU_double, p, inverse_double);
 	
-	//copy float matrix to double
+	// Copy float matrix to double.
 	for (size_t i=0; i<inverse->size1; i++) {
 		for (size_t j=0; j<inverse->size2; j++) {
-			gsl_matrix_float_set (inverse, i, j, gsl_matrix_get(inverse_double, i, j));
+			gsl_matrix_float_set (inverse, i, j,
+					      gsl_matrix_get(inverse_double,
+							     i, j));
 		}
 	}
 	
-	//free matrix
+	// Free matrix.
 	gsl_matrix_free (LU_double);
 	gsl_matrix_free (inverse_double);
 	return;
 }
 
-//LU lndet
-double LULndet (gsl_matrix *LU)
-{
+// LU lndet.
+double LULndet (gsl_matrix *LU) {
 	double d;
 	d=gsl_linalg_LU_lndet (LU);
 	return d;
 }
 
-double LULndet (gsl_matrix_float *LU)
-{
+double LULndet (gsl_matrix_float *LU) {
 	gsl_matrix *LU_double=gsl_matrix_alloc (LU->size1, LU->size2);
 	double d;
 	
-	//copy float matrix to double	
+	// Copy float matrix to double.	
 	for (size_t i=0; i<LU->size1; i++) {
 		for (size_t j=0; j<LU->size2; j++) {
 			gsl_matrix_set (LU_double, i, j, gsl_matrix_float_get(LU, i, j));
 		}
 	}
 	
-	//LU decomposition
+	// LU decomposition.
 	d=gsl_linalg_LU_lndet (LU_double);
 	
-	//copy float matrix to double
-	/*
-	for (size_t i=0; i<LU->size1; i++) {
-		for (size_t j=0; j<LU->size2; j++) {
-			gsl_matrix_float_set (LU, i, j, gsl_matrix_get(LU_double, i, j));
-		}
-	}
-	*/
-	//free matrix
+	// Free matrix
 	gsl_matrix_free (LU_double);
 	return d;
 }
 
 
-//LU solve
-void LUSolve (const gsl_matrix *LU, const gsl_permutation *p, const gsl_vector *b, gsl_vector *x)
-{
+// LU solve.
+void LUSolve (const gsl_matrix *LU, const gsl_permutation *p,
+	      const gsl_vector *b, gsl_vector *x) {
 	gsl_linalg_LU_solve (LU, p, b, x);
 	return;
 }
 
-void LUSolve (const gsl_matrix_float *LU, const gsl_permutation *p, const gsl_vector_float *b, gsl_vector_float *x)
-{
+void LUSolve (const gsl_matrix_float *LU, const gsl_permutation *p,
+	      const gsl_vector_float *b, gsl_vector_float *x) {
 	gsl_matrix *LU_double=gsl_matrix_alloc (LU->size1, LU->size2);
 	gsl_vector *b_double=gsl_vector_alloc (b->size);
 	gsl_vector *x_double=gsl_vector_alloc (x->size);	
 	
-	//copy float matrix to double	
+	// Copy float matrix to double.
 	for (size_t i=0; i<LU->size1; i++) {
 		for (size_t j=0; j<LU->size2; j++) {
-			gsl_matrix_set (LU_double, i, j, gsl_matrix_float_get(LU, i, j));
+			gsl_matrix_set (LU_double, i, j,
+					gsl_matrix_float_get(LU, i, j));
 		}
 	}
 	
@@ -594,15 +679,15 @@ void LUSolve (const gsl_matrix_float *LU, const gsl_permutation *p, const gsl_ve
 		gsl_vector_set (x_double, i, gsl_vector_float_get(x, i));
 	}
 	
-	//LU decomposition
+	// LU decomposition.
 	gsl_linalg_LU_solve (LU_double, p, b_double, x_double);
 	
-	//copy float matrix to double
+	// Copy float matrix to double.
 	for (size_t i=0; i<x->size; i++) {
 		gsl_vector_float_set (x, i, gsl_vector_get(x_double, i));
 	}
 	
-	//free matrix
+	// Free matrix.
 	gsl_matrix_free (LU_double);
 	gsl_vector_free (b_double);
 	gsl_vector_free (x_double);
@@ -610,8 +695,7 @@ void LUSolve (const gsl_matrix_float *LU, const gsl_permutation *p, const gsl_ve
 }
 
 
-bool lapack_ddot(vector<double> &x, vector<double> &y, double &v)
-{
+bool lapack_ddot(vector<double> &x, vector<double> &y, double &v) {
   bool flag=false;
   int incx=1;
   int incy=1;
@@ -625,8 +709,7 @@ bool lapack_ddot(vector<double> &x, vector<double> &y, double &v)
 }
 
 
-bool lapack_sdot(vector<float> &x, vector<float> &y, double &v)
-{
+bool lapack_sdot(vector<float> &x, vector<float> &y, double &v) {
   bool flag=false;
   int incx=1;
   int incy=1;
diff --git a/src/lapack.h b/src/lapack.h
index e5a1d37..5277b2f 100644
--- a/src/lapack.h
+++ b/src/lapack.h
@@ -1,6 +1,6 @@
 /*
-	Genome-wide Efficient Mixed Model Association (GEMMA)
-    Copyright (C) 2011  Xiang Zhou
+    Genome-wide Efficient Mixed Model Association (GEMMA)
+    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 <http://www.gnu.org/licenses/>.
+    along with this program. If not, see <http://www.gnu.org/licenses/>.
 */
 
 #ifndef __LAPACK_H__                
@@ -23,31 +23,47 @@
 
 using namespace std;
 
-
 void lapack_float_cholesky_decomp (gsl_matrix_float *A);
 void lapack_cholesky_decomp (gsl_matrix *A);
-void lapack_float_cholesky_solve (gsl_matrix_float *A, const gsl_vector_float *b, gsl_vector_float *x);
+void lapack_float_cholesky_solve (gsl_matrix_float *A,
+				  const gsl_vector_float *b,
+				  gsl_vector_float *x);
 void lapack_cholesky_solve (gsl_matrix *A, const gsl_vector *b, gsl_vector *x);
-void lapack_sgemm (char *TransA, char *TransB, float alpha, const gsl_matrix_float *A, const gsl_matrix_float *B, float beta, gsl_matrix_float *C);
-void lapack_dgemm (char *TransA, char *TransB, double alpha, const gsl_matrix *A, const gsl_matrix *B, double beta, gsl_matrix *C);
-void lapack_float_eigen_symmv (gsl_matrix_float *A, gsl_vector_float *eval, gsl_matrix_float *evec, const size_t flag_largematrix);
-void lapack_eigen_symmv (gsl_matrix *A, gsl_vector *eval, gsl_matrix *evec, const size_t flag_largematrix);
+void lapack_sgemm (char *TransA, char *TransB, float alpha,
+		   const gsl_matrix_float *A, const gsl_matrix_float *B,
+		   float beta, gsl_matrix_float *C);
+void lapack_dgemm (char *TransA, char *TransB, double alpha,
+		   const gsl_matrix *A, const gsl_matrix *B,
+		   double beta, gsl_matrix *C);
+void lapack_float_eigen_symmv (gsl_matrix_float *A, gsl_vector_float *eval,
+			       gsl_matrix_float *evec,
+			       const size_t flag_largematrix);
+void lapack_eigen_symmv (gsl_matrix *A, gsl_vector *eval, gsl_matrix *evec,
+			 const size_t flag_largematrix);
 
-double EigenDecomp (gsl_matrix *G, gsl_matrix *U, gsl_vector *eval, const size_t flag_largematrix);
-double EigenDecomp (gsl_matrix_float *G, gsl_matrix_float *U, gsl_vector_float *eval, const size_t flag_largematrix);
+double EigenDecomp (gsl_matrix *G, gsl_matrix *U, gsl_vector *eval,
+		    const size_t flag_largematrix);
+double EigenDecomp (gsl_matrix_float *G, gsl_matrix_float *U,
+		    gsl_vector_float *eval, const size_t flag_largematrix);
 
 double CholeskySolve(gsl_matrix *Omega, gsl_vector *Xty, gsl_vector *OiXty);
-double CholeskySolve(gsl_matrix_float *Omega, gsl_vector_float *Xty, gsl_vector_float *OiXty);
+double CholeskySolve(gsl_matrix_float *Omega, gsl_vector_float *Xty,
+		     gsl_vector_float *OiXty);
 
 void LUDecomp (gsl_matrix *LU, gsl_permutation *p, int *signum);
 void LUDecomp (gsl_matrix_float *LU, gsl_permutation *p, int *signum);
-void LUInvert (const gsl_matrix *LU, const gsl_permutation *p, gsl_matrix *inverse);
-void LUInvert (const gsl_matrix_float *LU, const gsl_permutation *p, gsl_matrix_float *inverse);
+void LUInvert (const gsl_matrix *LU, const gsl_permutation *p,
+	       gsl_matrix *inverse);
+void LUInvert (const gsl_matrix_float *LU, const gsl_permutation *p,
+	       gsl_matrix_float *inverse);
 double LULndet (gsl_matrix *LU);
 double LULndet (gsl_matrix_float *LU);
-void LUSolve (const gsl_matrix *LU, const gsl_permutation *p, const gsl_vector *b, gsl_vector *x);
-void LUSolve (const gsl_matrix_float *LU, const gsl_permutation *p, const gsl_vector_float *b, gsl_vector_float *x);
+void LUSolve (const gsl_matrix *LU, const gsl_permutation *p,
+	      const gsl_vector *b, gsl_vector *x);
+void LUSolve (const gsl_matrix_float *LU, const gsl_permutation *p,
+	      const gsl_vector_float *b, gsl_vector_float *x);
 
 bool lapack_ddot(vector<double> &x, vector<double> &y, double &v);
 bool lapack_sdot(vector<float> &x, vector<float> &y, double &v);
+
 #endif
diff --git a/src/param.cpp b/src/param.cpp
index 4b8c3a4..63bf349 100644
--- a/src/param.cpp
+++ b/src/param.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 <http://www.gnu.org/licenses/>.
+    along with this program. If not, see <http://www.gnu.org/licenses/>.
 */
 
 #include <iostream>
@@ -33,29 +33,20 @@
 
 #include "eigenlib.h"
 #include "mathfunc.h"
-
-#ifdef FORCE_FLOAT
-#include "param_float.h"
-#include "io_float.h"
-#else
 #include "param.h"
 #include "io.h"
-#endif
 
 using namespace std;
 
-
-
-
-
 PARAM::PARAM(void):
 mode_silence (false), a_mode (0), k_mode(1), d_pace (100000),
 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),
+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), noconstrain (false),
-h_min(-1), h_max(-1),	h_scale(-1),
-rho_min(0.0), rho_max(1.0),	rho_scale(-1),
+h_min(-1), h_max(-1), h_scale(-1),
+rho_min(0.0), rho_max(1.0), rho_scale(-1),
 logp_min(0.0), logp_max(0.0), logp_scale(-1),
 h_ngrid(10), rho_ngrid(10),
 s_min(0), s_max(300),
@@ -68,31 +59,22 @@ randseed(-1),
 window_cm(0), window_bp(0), window_ns(0), n_block(200),
 error(false),
 ni_subsample(0), n_cvt(1), n_vc(1), n_cat(0),
-time_total(0.0), time_G(0.0), time_eigen(0.0), time_UtX(0.0), time_UtZ(0.0), time_opt(0.0), time_Omega(0.0)
+time_total(0.0), time_G(0.0), time_eigen(0.0), time_UtX(0.0),
+time_UtZ(0.0), time_opt(0.0), time_Omega(0.0)
 {}
 
-
-//read files
-//obtain ns_total, ng_total, ns_test, ni_test
-void PARAM::ReadFiles (void)
-{
+// Read files: obtain ns_total, ng_total, ns_test, ni_test.
+void PARAM::ReadFiles (void) {
 	string file_str;
-	/*
-	//read continuous cat file
-	if (!file_mcatc.empty()) {
-	  if (ReadFile_mcatc (file_mcatc, mapRS2catc, n_cat)==false) {error=true;}
-	} else if (!file_catc.empty()) {
-	  if (ReadFile_catc (file_catc, mapRS2catc, n_cat)==false) {error=true;}
-	}
-	*/
-	//read cat file
+	
+	// Read cat file.
 	if (!file_mcat.empty()) {
 	  if (ReadFile_mcat (file_mcat, mapRS2cat, n_vc)==false) {error=true;}
 	} else if (!file_cat.empty()) {
 	  if (ReadFile_cat (file_cat, mapRS2cat, n_vc)==false) {error=true;}
 	}
 
-	//read snp weight files
+	// Read snp weight files.
 	if (!file_wcat.empty()) {
 	  if (ReadFile_wsnp (file_wcat, n_vc, mapRS2wcat)==false) {error=true;}
 	}
@@ -100,45 +82,60 @@ void PARAM::ReadFiles (void)
 	  if (ReadFile_wsnp (file_wsnp, mapRS2wsnp)==false) {error=true;}
 	}
 
-	//count number of kinship files
+	// Count number of kinship files.
 	if (!file_mk.empty()) {
 	  if (CountFileLines (file_mk, n_vc)==false) {error=true;}
 	}
 
-	//read snp set
+	// Read SNP set.
 	if (!file_snps.empty()) {
 		if (ReadFile_snps (file_snps, setSnps)==false) {error=true;}
 	} else {
 		setSnps.clear();
 	}
 
-	//for prediction
+	// For prediction.
 	if (!file_epm.empty()) {
-		if (ReadFile_est (file_epm, est_column, mapRS2est)==false) {error=true;}
-
+		if (ReadFile_est (file_epm, est_column, mapRS2est)==false) {
+		  error=true;
+		}
 		if (!file_bfile.empty()) {
 			file_str=file_bfile+".bim";
-			if (ReadFile_bim (file_str, snpInfo)==false) {error=true;}
-
+			if (ReadFile_bim (file_str, snpInfo)==false) {
+			  error=true;
+			}
 			file_str=file_bfile+".fam";
-			if (ReadFile_fam (file_str, indicator_pheno, pheno, mapID2num, p_column)==false) {error=true;}
+			if (ReadFile_fam (file_str, indicator_pheno, pheno,
+					  mapID2num, p_column)==false) {
+			  error=true;
+			}
 		}
 
 		if (!file_geno.empty()) {
-			if (ReadFile_pheno (file_pheno, indicator_pheno, pheno, p_column)==false) {error=true;}
+			if (ReadFile_pheno (file_pheno, indicator_pheno,
+					    pheno, p_column)==false) {
+			  error=true;
+			}
 
-			if (CountFileLines (file_geno, ns_total)==false) {error=true;}
+			if (CountFileLines (file_geno, ns_total)==false) {
+			  error=true;
+			}
 		}
 
 		if (!file_ebv.empty() ) {
-			if (ReadFile_column (file_ebv, indicator_bv, vec_bv, 1)==false) {error=true;}
+			if (ReadFile_column (file_ebv, indicator_bv,
+					     vec_bv, 1)==false) {
+			  error=true;
+			}
 		}
 
 		if (!file_log.empty() ) {
-			if (ReadFile_log (file_log, pheno_mean)==false) {error=true;}
+			if (ReadFile_log (file_log, pheno_mean)==false) {
+			  error=true;
+			}
 		}
 
-		//convert indicator_pheno to indicator_idv
+		// Convert indicator_pheno to indicator_idv.
 		int k=1;
 		for (size_t i=0; i<indicator_pheno.size(); i++) {
 			k=1;
@@ -153,10 +150,12 @@ void PARAM::ReadFiles (void)
 		return;
 	}
 
-	//read covariates before the genotype files
+	// Read covariates before the genotype files.
 	if (!file_cvt.empty() ) {
-		if (ReadFile_cvt (file_cvt, indicator_cvt, cvt, n_cvt)==false) {error=true;}
-
+		if (ReadFile_cvt (file_cvt, indicator_cvt,
+				  cvt, n_cvt)==false) {
+		  error=true;
+		}
 		if ((indicator_cvt).size()==0) {
 			n_cvt=1;
 		}
@@ -165,102 +164,135 @@ void PARAM::ReadFiles (void)
 	}
 
 	if (!file_gxe.empty() ) {
-	  if (ReadFile_column (file_gxe, indicator_gxe, gxe, 1)==false) {error=true;}
+	  if (ReadFile_column (file_gxe, indicator_gxe, gxe, 1)==false) {
+	    error=true;
+	  }
 	}
 	if (!file_weight.empty() ) {
-	  if (ReadFile_column (file_weight, indicator_weight, weight, 1)==false) {error=true;}
+	  if (ReadFile_column (file_weight, indicator_weight,
+			       weight, 1)==false) {
+	    error=true;
+	  }
 	}
 
 
-	// WJA added
-	//read genotype and phenotype file for bgen format
+	// WJA added.
+	// Read genotype and phenotype file for bgen format.
 	if (!file_oxford.empty()) {
 		file_str=file_oxford+".sample";
-		if (ReadFile_sample(file_str, indicator_pheno, pheno, p_column,indicator_cvt, cvt, n_cvt)==false) {error=true;}
+		if (ReadFile_sample(file_str, indicator_pheno, pheno, p_column,
+				    indicator_cvt, cvt, n_cvt)==false) {
+		  error=true;
+		}
 		if ((indicator_cvt).size()==0) {
 			n_cvt=1;
 		}
-	//	n_cvt=1;
 
-		//post-process covariates and phenotypes, obtain ni_test, save all useful covariates
+		// Post-process covariates and phenotypes, obtain
+		// ni_test, save all useful covariates.
 		ProcessCvtPhen();
 
-		//obtain covariate matrix
+		// Obtain covariate matrix.
 		gsl_matrix *W=gsl_matrix_alloc (ni_test, n_cvt);
 		CopyCvt (W);
 
 		file_str=file_oxford+".bgen";
-		if (ReadFile_bgen (file_str, setSnps, W, indicator_idv, indicator_snp, snpInfo, maf_level, miss_level, hwe_level, r2_level, ns_test)==false) {error=true;}
+		if (ReadFile_bgen (file_str, setSnps, W, indicator_idv,
+				   indicator_snp, snpInfo, maf_level,
+				   miss_level, hwe_level, r2_level,
+				   ns_test)==false) {
+		  error=true;
+		}
 		gsl_matrix_free(W);
 
 		ns_total=indicator_snp.size();
 	}
 
 
-	//read genotype and phenotype file for plink format
+	// Read genotype and phenotype file for PLINK format.
 	if (!file_bfile.empty()) {
 		file_str=file_bfile+".bim";
 		snpInfo.clear();
 		if (ReadFile_bim (file_str, snpInfo)==false) {error=true;}
 
-		//if both fam file and pheno files are used, use phenotypes inside the pheno file
+		// If both fam file and pheno files are used, use
+		// phenotypes inside the pheno file.
 		if (!file_pheno.empty()) {
-		  //phenotype file before genotype file
-		  if (ReadFile_pheno (file_pheno, indicator_pheno, pheno, p_column)==false) {error=true;}
+		  
+		  // Phenotype file before genotype file.
+		  if (ReadFile_pheno (file_pheno, indicator_pheno, pheno,
+				      p_column)==false) {error=true;}
 		} else {
 		  file_str=file_bfile+".fam";
-		  if (ReadFile_fam (file_str, indicator_pheno, pheno, mapID2num, p_column)==false) {error=true;}
+		  if (ReadFile_fam (file_str, indicator_pheno, pheno,
+				    mapID2num, p_column)==false) {error=true;}
 		}
 
-		//post-process covariates and phenotypes, obtain ni_test, save all useful covariates
+		// Post-process covariates and phenotypes, obtain
+		// ni_test, save all useful covariates.
 		ProcessCvtPhen();
 
-		//obtain covariate matrix
+		// Obtain covariate matrix.
 		gsl_matrix *W=gsl_matrix_alloc (ni_test, n_cvt);
 		CopyCvt (W);
 
 		file_str=file_bfile+".bed";
-		if (ReadFile_bed (file_str, setSnps, W, indicator_idv, indicator_snp, snpInfo, maf_level, miss_level, hwe_level, r2_level, ns_test)==false) {error=true;}
-
+		if (ReadFile_bed (file_str, setSnps, W, indicator_idv,
+				  indicator_snp, snpInfo, maf_level,
+				  miss_level, hwe_level, r2_level,
+				  ns_test) == false) {
+		  error=true;
+		}
 		gsl_matrix_free(W);
-
 		ns_total=indicator_snp.size();
 	}
 
-	//read genotype and phenotype file for bimbam format
+	// Read genotype and phenotype file for BIMBAM format.
 	if (!file_geno.empty()) {
-		//annotation file before genotype file
+	  
+	        // Annotation file before genotype file.
 		if (!file_anno.empty() ) {
-			if (ReadFile_anno (file_anno, mapRS2chr, mapRS2bp, mapRS2cM)==false) {error=true;}
+			if (ReadFile_anno (file_anno, mapRS2chr, mapRS2bp,
+					   mapRS2cM)==false) {
+			  error=true;
+			}
 		}
 
-		//phenotype file before genotype file
-		if (ReadFile_pheno (file_pheno, indicator_pheno, pheno, p_column)==false) {error=true;}
+		// Phenotype file before genotype file.
+		if (ReadFile_pheno (file_pheno, indicator_pheno, pheno,
+				    p_column) == false) {
+		  error=true;
+		}
 
-		//post-process covariates and phenotypes, obtain ni_test, save all useful covariates
+		// Post-process covariates and phenotypes, obtain
+		// ni_test, save all useful covariates.
 		ProcessCvtPhen();
 
-		//obtain covariate matrix
+		// Obtain covariate matrix.
 		gsl_matrix *W=gsl_matrix_alloc (ni_test, n_cvt);
 		CopyCvt (W);
 
-		if (ReadFile_geno (file_geno, setSnps, W, indicator_idv, indicator_snp, maf_level, miss_level, hwe_level, r2_level, mapRS2chr, mapRS2bp, mapRS2cM, snpInfo, ns_test)==false) {error=true;}
-
+		if (ReadFile_geno (file_geno, setSnps, W, indicator_idv,
+				   indicator_snp, maf_level, miss_level,
+				   hwe_level, r2_level, mapRS2chr, mapRS2bp,
+				   mapRS2cM, snpInfo, ns_test)==false) {
+		  error=true;
+		}
 		gsl_matrix_free(W);
-
 		ns_total=indicator_snp.size();
 	}
 
 
-	//read genotype file for multiple plink files
+	// Read genotype file for multiple PLINK files.
 	if (!file_mbfile.empty()) {
 	  igzstream infile (file_mbfile.c_str(), igzstream::in);
-	  if (!infile) {cout<<"error! fail to open mbfile file: "<<file_mbfile<<endl; return;}
+	  if (!infile) {
+	    cout<<"error! fail to open mbfile file: " << file_mbfile<<endl;
+	    return;
+	  }
 
 	  string file_name;
-
 	  size_t t=0, ns_test_tmp=0;
-
 	  gsl_matrix *W;
 	  while (!safeGetline(infile, file_name).eof()) {
 		file_str=file_name+".bim";
@@ -268,25 +300,40 @@ void PARAM::ReadFiles (void)
 		if (ReadFile_bim (file_str, snpInfo)==false) {error=true;}
 
 		if (t==0) {
-		  //if both fam file and pheno files are used, use phenotypes inside the pheno file
+		  
+		  // If both fam file and pheno files are used, use
+		  // phenotypes inside the pheno file.
 		  if (!file_pheno.empty()) {
-		    //phenotype file before genotype file
-		    if (ReadFile_pheno (file_pheno, indicator_pheno, pheno, p_column)==false) {error=true;}
+		    
+		    // Phenotype file before genotype file.
+		    if (ReadFile_pheno (file_pheno, indicator_pheno, pheno,
+					p_column)==false) {
+		      error=true;
+		    }
 		  } else {
 		    file_str=file_name+".fam";
-		    if (ReadFile_fam (file_str, indicator_pheno, pheno, mapID2num, p_column)==false) {error=true;}
+		    if (ReadFile_fam (file_str, indicator_pheno, pheno,
+				      mapID2num, p_column)==false) {
+		      error=true;
+		    }
 		  }
 
-		  //post-process covariates and phenotypes, obtain ni_test, save all useful covariates
+		  // Post-process covariates and phenotypes, obtain
+		  // ni_test, save all useful covariates.
 		  ProcessCvtPhen();
 
-		  //obtain covariate matrix
+		  // Obtain covariate matrix.
 		  W=gsl_matrix_alloc (ni_test, n_cvt);
 		  CopyCvt (W);
 		}
 
 		file_str=file_name+".bed";
-		if (ReadFile_bed (file_str, setSnps, W, indicator_idv, indicator_snp, snpInfo, maf_level, miss_level, hwe_level, r2_level, ns_test_tmp)==false) {error=true;}
+		if (ReadFile_bed (file_str, setSnps, W, indicator_idv,
+				  indicator_snp, snpInfo, maf_level,
+				  miss_level, hwe_level, r2_level,
+				  ns_test_tmp)==false) {
+		  error=true;
+		}
 		mindicator_snp.push_back(indicator_snp);
 		msnpInfo.push_back(snpInfo);
 		ns_test+=ns_test_tmp;
@@ -303,30 +350,46 @@ void PARAM::ReadFiles (void)
 
 
 
-	//read genotype and phenotype file for multiple bimbam files
+	// Read genotype and phenotype file for multiple BIMBAM files.
 	if (!file_mgeno.empty()) {
-	  //annotation file before genotype file
+	  
+	  // Annotation file before genotype file.
 	  if (!file_anno.empty() ) {
-	    if (ReadFile_anno (file_anno, mapRS2chr, mapRS2bp, mapRS2cM)==false) {error=true;}
+	    if (ReadFile_anno (file_anno, mapRS2chr, mapRS2bp,
+			       mapRS2cM)==false) {
+	      error=true;
+	    }
 	  }
 
-	  //phenotype file before genotype file
-	  if (ReadFile_pheno (file_pheno, indicator_pheno, pheno, p_column)==false) {error=true;}
+	  // Phenotype file before genotype file.
+	  if (ReadFile_pheno (file_pheno, indicator_pheno, pheno,
+			      p_column)==false) {
+	    error=true;
+	  }
 
-	  //post-process covariates and phenotypes, obtain ni_test, save all useful covariates
+	  // Post-process covariates and phenotypes, obtain ni_test,
+	  // save all useful covariates.
 	  ProcessCvtPhen();
 
-	  //obtain covariate matrix
+	  // Obtain covariate matrix.
 	  gsl_matrix *W=gsl_matrix_alloc (ni_test, n_cvt);
 	  CopyCvt (W);
 
 	  igzstream infile (file_mgeno.c_str(), igzstream::in);
-	  if (!infile) {cout<<"error! fail to open mgeno file: "<<file_mgeno<<endl; return;}
+	  if (!infile) {
+	    cout<<"error! fail to open mgeno file: "<<file_mgeno<<endl;
+	    return;
+	  }
 
 	  string file_name;
 	  size_t ns_test_tmp;
 	  while (!safeGetline(infile, file_name).eof()) {
-	    if (ReadFile_geno (file_name, setSnps, W, indicator_idv, indicator_snp, maf_level, miss_level, hwe_level, r2_level, mapRS2chr, mapRS2bp, mapRS2cM, snpInfo, ns_test_tmp)==false) {error=true;}
+	    if (ReadFile_geno (file_name, setSnps, W, indicator_idv,
+			       indicator_snp, maf_level, miss_level,
+			       hwe_level, r2_level, mapRS2chr, mapRS2bp,
+			       mapRS2cM, snpInfo, ns_test_tmp)==false) {
+	      error=true;
+	    }
 
 	    mindicator_snp.push_back(indicator_snp);
 	    msnpInfo.push_back(snpInfo);
@@ -343,9 +406,10 @@ void PARAM::ReadFiles (void)
 
 
 	if (!file_gene.empty()) {
-		if (ReadFile_pheno (file_pheno, indicator_pheno, pheno, p_column)==false) {error=true;}
+		if (ReadFile_pheno (file_pheno, indicator_pheno, pheno,
+				    p_column)==false) {error=true;}
 
-		//convert indicator_pheno to indicator_idv
+		// Convert indicator_pheno to indicator_idv.
 		int k=1;
 		for (size_t i=0; i<indicator_pheno.size(); i++) {
 			k=1;
@@ -355,57 +419,69 @@ void PARAM::ReadFiles (void)
 			indicator_idv.push_back(k);
 		}
 
-		//post-process covariates and phenotypes, obtain ni_test, save all useful covariates
+		// Post-process covariates and phenotypes, obtain
+		// ni_test, save all useful covariates.
 		ProcessCvtPhen();
 
-		//obtain covariate matrix
+		// Obtain covariate matrix.
 		gsl_matrix *W=gsl_matrix_alloc (ni_test, n_cvt);
 		CopyCvt (W);
 
-		if (ReadFile_gene (file_gene, vec_read, snpInfo, ng_total)==false) {error=true;}
+		if (ReadFile_gene (file_gene, vec_read, snpInfo,
+				   ng_total)==false) {
+		  error=true;
+		}
 	}
 
 
-	//read is after gene file
+	// Read is after gene file.
 	if (!file_read.empty() ) {
-		if (ReadFile_column (file_read, indicator_read, vec_read, 1)==false) {error=true;}
+		if (ReadFile_column (file_read, indicator_read,
+				     vec_read, 1)==false) {
+		  error=true;
+		}
 
 		ni_test=0;
-		for (vector<int>::size_type i=0; i<(indicator_idv).size(); ++i) {
+		for (vector<int>::size_type i=0;
+		     i<(indicator_idv).size();
+		     ++i) {
 			indicator_idv[i]*=indicator_read[i];
 			ni_test+=indicator_idv[i];
 		}
 
 		if (ni_test==0) {
-			error=true;
-			cout<<"error! number of analyzed individuals equals 0. "<<endl;
-			return;
+		  error=true;
+		  cout<<"error! number of analyzed individuals equals 0. "<<
+		    endl;
+		  return;
 		}
 	}
 
-	//for ridge prediction, read phenotype only
+	// For ridge prediction, read phenotype only.
 	if (file_geno.empty() && file_gene.empty() && !file_pheno.empty()) {
-		if (ReadFile_pheno (file_pheno, indicator_pheno, pheno, p_column)==false) {error=true;}
+		if (ReadFile_pheno (file_pheno, indicator_pheno, pheno,
+				    p_column)==false) {
+		  error=true;
+		}
 
-		//post-process covariates and phenotypes, obtain ni_test, save all useful covariates
+		// Post-process covariates and phenotypes, obtain
+		// ni_test, save all useful covariates.
 		ProcessCvtPhen();
 	}
-
 	return;
 }
 
-
-
-
-
-
 void PARAM::CheckParam (void)
 {
 	struct stat fileInfo;
 	string str;
 
-	//check parameters
-	if (k_mode!=1 && k_mode!=2) {cout<<"error! unknown kinship/relatedness input mode: "<<k_mode<<endl; error=true;}
+	// Check parameters.
+	if (k_mode!=1 && k_mode!=2) {
+	  cout<<"error! unknown kinship/relatedness input mode: "<<
+	    k_mode<<endl;
+	  error=true;
+	}
 	if (a_mode!=1 && a_mode!=2 && a_mode!=3 && a_mode!=4 && a_mode!=5 && a_mode!=11 && a_mode!=12 && a_mode!=13 && a_mode!=14 && a_mode!=15 && a_mode!=21 && a_mode!=22 && a_mode!=25 && a_mode!=26 && a_mode!=27 && a_mode!=28 && a_mode!=31 && a_mode!=41 && a_mode!=42 && a_mode!=43 && a_mode!=51 && a_mode!=52 && a_mode!=53 && a_mode!=54 && a_mode!=61 && a_mode!=62 && a_mode!=63 && a_mode!=66 && a_mode!=67 && a_mode!=71)
 	{cout<<"error! unknown analysis mode: "<<a_mode<<". make sure -gk or -eigen or -lmm or -bslmm -predict or -calccov is sepcified correctly."<<endl; error=true;}
 	if (miss_level>1) {cout<<"error! missing level needs to be between 0 and 1. current value = "<<miss_level<<endl; error=true;}
@@ -447,13 +523,10 @@ void PARAM::CheckParam (void)
 		}
 	}
 
-	//sort (p_column.begin(), p_column.end() );
 	n_ph=p_column.size();
 
-
-
-	//only lmm option (and one prediction option) can deal with multiple phenotypes
-	//and no gene expression files
+	// Only LMM option (and one prediction option) can deal with
+	// multiple phenotypes and no gene expression files.
 	if (n_ph>1 && a_mode!=1 && a_mode!=2 && a_mode!=3 && a_mode!=4 && a_mode!=43) {
 		cout<<"error! the current analysis mode "<<a_mode<<" can not deal with multiple phenotypes."<<endl; error=true;
 	}
diff --git a/src/param.h b/src/param.h
index 72b7d85..18d5e36 100644
--- a/src/param.h
+++ b/src/param.h
@@ -1,6 +1,6 @@
 /*
-	Genome-wide Efficient Mixed Model Association (GEMMA)
-    Copyright (C) 2011  Xiang Zhou
+    Genome-wide Efficient Mixed Model Association (GEMMA)
+    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 <http://www.gnu.org/licenses/>.
+    along with this program. If not, see <http://www.gnu.org/licenses/>.
 */
 
 #ifndef __PARAM_H__
@@ -27,8 +27,6 @@
 
 using namespace std;
 
-
-
 class SNPINFO {
 public:
 	string chr;
@@ -40,37 +38,36 @@ public:
 	size_t n_miss;
 	double missingness;
 	double maf;
-	size_t n_idv;//number of non-missing individuals
-	size_t n_nb;//number of neighbours on the right hand side
-	size_t file_position;//snp location on file
+	size_t n_idv;         // Number of non-missing individuals.
+	size_t n_nb;          // Number of neighbours on the right hand side.
+	size_t file_position; // SNP location in file.
 };
 
-//results for lmm
+// Results for LMM.
 class SUMSTAT {
 public:
-	double beta;				//REML estimator for beta
-	double se;				//SE for beta
-	double lambda_remle;		//REML estimator for lambda
-	double lambda_mle;		//MLE estimator for lambda
-	double p_wald;			//p value from a Wald test
-	double p_lrt;				//p value from a likelihood ratio test
-	double p_score;			//p value from a score test
+	double beta;		// REML estimator for beta.
+	double se;		// SE for beta.
+	double lambda_remle;	// REML estimator for lambda.
+	double lambda_mle;	// MLE estimator for lambda.
+	double p_wald;		// p value from a Wald test.
+	double p_lrt;		// p value from a likelihood ratio test.
+	double p_score;		// p value from a score test.
 };
 
-//results for mvlmm
+// Results for mvLMM.
 class MPHSUMSTAT {
 public:
-	vector<double> v_beta;	//REML estimator for beta
-	double p_wald;			//p value from a Wald test
-	double p_lrt;				//p value from a likelihood ratio test
-	double p_score;			//p value from a score test
-	vector<double> v_Vg;	//estimator for Vg, right half
-	vector<double> v_Ve;	//estimator for Ve, right half
-	vector<double> v_Vbeta;	//estimator for Vbeta, right half
+	vector<double> v_beta;	// REML estimator for beta.
+	double p_wald;		// p value from a Wald test.
+	double p_lrt;		// p value from a likelihood ratio test.
+	double p_score;		// p value from a score test.
+	vector<double> v_Vg;	// Estimator for Vg, right half.
+	vector<double> v_Ve;	// Estimator for Ve, right half.
+	vector<double> v_Vbeta;	// Estimator for Vbeta, right half.
 };
 
-
-//hyper-parameters for bslmm
+// Hyper-parameters for BSLMM.
 class HYPBSLMM {
 public:
 	double h;
@@ -78,15 +75,11 @@ public:
 	double rho;
 	double pge;
 	double logp;
-
 	size_t n_gamma;
 };
 
-
-//header class
-class HEADER
-{
-
+// Header class.
+class HEADER {
 public:
     size_t rs_col;
     size_t chr_col;
@@ -108,28 +101,26 @@ public:
     size_t var_col;
     size_t ws_col;
     size_t cor_col;
-    size_t coln;//number of columns
+    size_t coln; // Number of columns.
     set<size_t> catc_col;
     set<size_t> catd_col;
 };
 
-
-
 class PARAM {
 public:
-	// IO related parameters
+	// IO-related parameters.
 	bool mode_silence;
-	int a_mode;				//analysis mode, 1/2/3/4 for Frequentist tests
-	int k_mode;				//kinship read mode: 1: n by n matrix, 2: id/id/k_value;
-	vector<size_t> p_column;			//which phenotype column needs analysis
-	size_t d_pace;		//display pace
+	int a_mode;  // Analysis mode, 1/2/3/4 for Frequentist tests
+        int k_mode;  // Kinship read mode: 1: n by n matrix, 2: id/id/k_value;
+	vector<size_t> p_column; // Which phenotype column needs analysis.
+	size_t d_pace;	 	 // Display pace
 
 	string file_bfile, file_mbfile;
 	string file_geno, file_mgeno;
 	string file_pheno;
-	string file_anno;		//optional
-	string file_gxe;		//optional
-	string file_cvt;		//optional
+	string file_anno; // Optional.
+	string file_gxe;  // Optional.
+	string file_cvt;  // Optional.
 	string file_cat, file_mcat;
 	string file_catc, file_mcatc;
 	string file_var;
@@ -144,26 +135,23 @@ public:
 	string file_bf, file_hyp;
 	string path_out;
 
-
-	string file_epm;		//estimated parameter file
-	string file_ebv;		//estimated breeding value file
-	string file_log;		//log file containing mean estimate
-
-	string file_read;		//file containing total number of reads
-	string file_gene;		//gene expression file
-
-	string file_snps;		//file containing analyzed snps or genes
-// WJA Added
+	string file_epm;  // Estimated parameter file.
+	string file_ebv;  // Estimated breeding value file.
+	string file_log;  // Log file containing mean estimate.
+	string file_read; // File containing total number of reads.
+	string file_gene; // Gene expression file.
+	string file_snps; // File containing analyzed SNPs or genes.
+  
+        // WJA added.
 	string file_oxford;
 
-
-	// QC related parameters
+	// QC-related parameters.
 	double miss_level;
 	double maf_level;
 	double hwe_level;
 	double r2_level;
 
-	// LMM related parameters
+	// LMM-related parameters.
 	double l_min;
 	double l_max;
 	size_t n_region;
@@ -172,16 +160,18 @@ public:
 	double pve_null, pve_se_null, pve_total, se_pve_total;
 	double vg_remle_null, ve_remle_null, vg_mle_null, ve_mle_null;
 	vector<double> Vg_remle_null, Ve_remle_null, Vg_mle_null, Ve_mle_null;
-	vector<double> VVg_remle_null, VVe_remle_null, VVg_mle_null, VVe_mle_null;
-	vector<double> beta_remle_null, se_beta_remle_null, beta_mle_null, se_beta_mle_null;
+        vector<double> VVg_remle_null, VVe_remle_null, VVg_mle_null;
+        vector<double> VVe_mle_null;
+        vector<double> beta_remle_null, se_beta_remle_null, beta_mle_null;
+        vector<double> se_beta_mle_null;
 	double p_nr;
 	double em_prec, nr_prec;
 	size_t em_iter, nr_iter;
 	size_t crt;
-	double pheno_mean;		//phenotype mean from bslmm fitting or for prediction
+	double pheno_mean; // Phenotype mean from BSLMM fitting or prediction.
 
-	//for fitting multiple variance components
-	//the first three are of size n_vc, and the next two are of size n_vc+1
+	// For fitting multiple variance components.
+	// The first 3 are of size (n_vc), and the next 2 are of size n_vc+1.
 	bool noconstrain;
 	vector<double> v_traceG;
 	vector<double> v_pve;
@@ -194,98 +184,141 @@ public:
 	vector<double> v_beta;
 	vector<double> v_se_beta;
 
-	// BSLMM MCMC related parameters
-	double h_min, h_max, h_scale;			//priors for h
-	double rho_min, rho_max, rho_scale;		//priors for rho
-	double logp_min, logp_max, logp_scale;		//priors for log(pi)
+	// BSLMM/MCMC-related parameters.
+	double h_min, h_max, h_scale;		// Priors for h.
+	double rho_min, rho_max, rho_scale;	// Priors for rho.
+	double logp_min, logp_max, logp_scale; 	// Priors for log(pi).
 	size_t h_ngrid, rho_ngrid;
-	size_t s_min, s_max;			//minimum and maximum number of gammas
-	size_t w_step;					//number of warm up/burn in iterations
-	size_t s_step;					//number of sampling iterations
-	size_t r_pace;					//record pace
-	size_t w_pace;					//write pace
-	size_t n_accept;				//number of acceptance
-	size_t n_mh;					//number of MH steps within each iteration
-	double geo_mean;				//mean of the geometric distribution
+	size_t s_min, s_max;			// Min & max. number of gammas.
+	size_t w_step;				// # warm up/burn in iter.
+	size_t s_step;				// # sampling iterations.
+	size_t r_pace;				// Record pace.
+	size_t w_pace;				// Write pace.
+	size_t n_accept;			// Number of acceptance.
+	size_t n_mh;				// # MH steps in each iter.
+	double geo_mean;			// Mean of geometric dist.
 	long int randseed;
 	double trace_G;
 
 	HYPBSLMM cHyp_initial;
 
-	//VARCOV related parameters
+	// VARCOV-related parameters.
 	double window_cm;
 	size_t window_bp;
 	size_t window_ns;
 
-	//vc related parameters
+	// vc-related parameters.
 	size_t n_block;
 
-	// Summary statistics
+	// Summary statistics.
 	bool error;
-	size_t ni_total, ni_test, ni_cvt, ni_study, ni_ref;	//number of individuals
-	size_t np_obs, np_miss;		//number of observed and missing phenotypes
-	size_t ns_total, ns_test, ns_study, ns_ref;	//number of snps
-	size_t ng_total, ng_test;	//number of genes
-	size_t ni_control, ni_case;	//number of controls and number of cases
-	size_t ni_subsample;            //number of subsampled individuals
-	//size_t ni_total_ref, ns_total_ref, ns_pair;//max number of individuals, number of snps and number of snp pairs in the reference panel
-	size_t n_cvt;			//number of covariates
-	size_t n_cat;			//number of continuous categories
-	size_t n_ph;			//number of phenotypes
-	size_t n_vc;			//number of variance components (including the diagonal matrix)
-	double time_total;		//record total time
-	double time_G;			//time spent on reading files the second time and calculate K
-	double time_eigen;		//time spent on eigen-decomposition
-	double time_UtX;		//time spent on calculating UX and Uy
-	double time_UtZ;		//time spent on calculating UtZ, for probit BSLMM
-	double time_opt;		//time spent on optimization iterations/or mcmc
-	double time_Omega;		//time spent on calculating Omega
-	double time_hyp;		//time spent on sampling hyper-parameters, in PMM
-	double time_Proposal;  //time spend on constructing the proposal distribution (i.e. the initial lmm or lm analysis)
-
-	// Data
-	vector<vector<double> > pheno;			//a vector record all phenotypes, NA replaced with -9
-	vector<vector<double> > cvt;			//a vector record all covariates, NA replaced with -9
-	vector<double> gxe;			//a vector record all covariates, NA replaced with -9
-	vector<double> weight;                          //a vector record weights for the individuals, which is useful for animal breeding studies
-	vector<vector<int> > indicator_pheno;			//a matrix record when a phenotype is missing for an individual; 0 missing, 1 available
-	vector<int> indicator_idv;				//indicator for individuals (phenotypes), 0 missing, 1 available for analysis
-	vector<int> indicator_snp;				//sequence indicator for SNPs: 0 ignored because of (a) maf, (b) miss, (c) non-poly; 1 available for analysis
-	vector< vector<int> >  mindicator_snp;				//sequence indicator for SNPs: 0 ignored because of (a) maf, (b) miss, (c) non-poly; 1 available for analysis
-	vector<int> indicator_cvt;				//indicator for covariates, 0 missing, 1 available for analysis
-	vector<int> indicator_gxe;				//indicator for gxe, 0 missing, 1 available for analysis
-	vector<int> indicator_weight;                           //indicator for weight, 0 missing, 1 available for analysis
-
-	vector<int> indicator_bv;				//indicator for estimated breeding value file, 0 missing, 1 available for analysis
-	vector<int> indicator_read;				//indicator for read file, 0 missing, 1 available for analysis
-	vector<double> vec_read;				//total number of reads
-	vector<double> vec_bv;					//breeding values
+  
+        // Number of individuals.
+	size_t ni_total, ni_test, ni_cvt, ni_study, ni_ref;
+
+        // Number of observed and missing phenotypes.
+	size_t np_obs, np_miss;
+
+        // Number of SNPs.
+	size_t ns_total, ns_test, ns_study, ns_ref;
+  
+	size_t ng_total, ng_test;   // Number of genes.
+	size_t ni_control, ni_case; // Number of controls and number of cases.
+	size_t ni_subsample;        // Number of subsampled individuals.
+	size_t n_cvt;		    // Number of covariates.
+	size_t n_cat;		    // Number of continuous categories.
+	size_t n_ph;		    // Number of phenotypes.
+	size_t n_vc;		    // Number of variance components
+                                    // (including the diagonal matrix).
+	double time_total;	    // Record total time.
+	double time_G;	 	    // Time spent on reading files the
+                                    // second time and calculate K.
+	double time_eigen;	    // Time spent on eigen-decomposition.
+	double time_UtX;	    // Time spent on calculating UX and Uy.
+	double time_UtZ;	    // Time calculating UtZ for probit BSLMM.
+	double time_opt;	    // Time on optimization iterations/MCMC.
+	double time_Omega;   	    // Time spent on calculating Omega.
+	double time_hyp;	    // Time sampling hyperparameters in PMM.
+	double time_Proposal;       // Time spent on constructing the
+				    // proposal distribution (i.e. the
+				    // initial LMM or LM analysis).
+
+	// Data.
+        // Vector recording all phenotypes (NA replaced with -9).
+	vector<vector<double> > pheno;
+
+        // Vector recording all covariates (NA replaced with -9).
+	vector<vector<double> > cvt;
+
+        // Vector recording all covariates (NA replaced with -9).
+	vector<double> gxe;
+
+        // Vector recording weights for the individuals, which is
+        // useful for animal breeding studies.
+	vector<double> weight;
+
+        // Matrix recording when a phenotype is missing for an
+        // individual; 0 missing, 1 available.
+	vector<vector<int> > indicator_pheno;
+
+        // Indicator for individuals (phenotypes): 0 missing, 1
+        // available for analysis
+	vector<int> indicator_idv;
+
+        // Sequence indicator for SNPs: 0 ignored because of (a) maf,
+        // (b) miss, (c) non-poly; 1 available for analysis.
+	vector<int> indicator_snp;
+
+        // Sequence indicator for SNPs: 0 ignored because of (a) maf,
+        // (b) miss, (c) non-poly; 1 available for analysis.
+	vector< vector<int> >  mindicator_snp;
+
+        // Indicator for covariates: 0 missing, 1 available for
+        // analysis.
+	vector<int> indicator_cvt;
+
+        // Indicator for gxe: 0 missing, 1 available for analysis.
+	vector<int> indicator_gxe;
+
+        // Indicator for weight: 0 missing, 1 available for analysis.
+	vector<int> indicator_weight;
+
+        // Indicator for estimated breeding value file: 0 missing, 1
+        // available for analysis.
+	vector<int> indicator_bv;
+
+	// Indicator for read file: 0 missing, 1 available for analysis.
+	vector<int> indicator_read;
+	vector<double> vec_read;	// Total number of reads.
+	vector<double> vec_bv;   	// Breeding values.
 	vector<size_t> est_column;
 
-	map<string, int> mapID2num;		//map small ID number to number, from 0 to n-1
-	map<string, string> mapRS2chr;		//map rs# to chromosome location
-	map<string, long int> mapRS2bp;		//map rs# to base position
-	map<string, double> mapRS2cM;		//map rs# to cM
-	map<string, double> mapRS2est;			//map rs# to parameters
-	map<string, size_t> mapRS2cat;          //map rs# to category number
-	map<string, vector<double> > mapRS2catc;          //map rs# to continuous categories
-	map<string, double> mapRS2wsnp;          //map rs# to snp weights
-	map<string, vector<double> > mapRS2wcat;          //map rs# to snp cat weights
-
-	vector<SNPINFO> snpInfo;		//record SNP information
-	vector< vector<SNPINFO> > msnpInfo;		//record SNP information
-	set<string> setSnps;			//a set of snps for analysis
-
-	//constructor
+	map<string, int> mapID2num;	// Map small ID to number, 0 to n-1.
+	map<string, string> mapRS2chr; 	// Map rs# to chromosome location.
+	map<string, long int> mapRS2bp;	// Map rs# to base position.
+	map<string, double> mapRS2cM;	// Map rs# to cM.
+	map<string, double> mapRS2est;	// Map rs# to parameters.
+	map<string, size_t> mapRS2cat;  // Map rs# to category number.
+	map<string, vector<double> > mapRS2catc; // Map rs# to cont. cat's.
+	map<string, double> mapRS2wsnp;          // Map rs# to SNP weights.
+	map<string, vector<double> > mapRS2wcat; // Map rs# to SNP cat weights.
+
+	vector<SNPINFO> snpInfo;	 	 // Record SNP information.
+	vector< vector<SNPINFO> > msnpInfo;	 // Record SNP information.
+	set<string> setSnps;			 // Set of snps for analysis.
+
+	// Constructor.
 	PARAM();
 
-	//functions
+	// Functions.
 	void ReadFiles ();
 	void CheckParam ();
 	void CheckData ();
 	void PrintSummary ();
-	void ReadGenotypes (gsl_matrix *UtX, gsl_matrix *K, const bool calc_K);
-	void ReadGenotypes (vector<vector<unsigned char> > &Xt, gsl_matrix *K, const bool calc_K);
+	void ReadGenotypes (gsl_matrix *UtX, gsl_matrix *K,
+			    const bool calc_K);
+	void ReadGenotypes (vector<vector<unsigned char> > &Xt,
+			    gsl_matrix *K, const bool calc_K);
 	void CheckCvt ();
 	void CopyCvt (gsl_matrix *W);
 	void CopyA (size_t flag, gsl_matrix *A);
@@ -295,15 +328,27 @@ public:
 	void CopyCvtPhen (gsl_matrix *W, gsl_vector *y, size_t flag);
 	void CopyCvtPhen (gsl_matrix *W, gsl_matrix *Y, size_t flag);
 	void CalcKin (gsl_matrix *matrix_kin);
-	void CalcS (const map<string, double> &mapRS2wA, const map<string, double> &mapRS2wK, const gsl_matrix *W, gsl_matrix *A, gsl_matrix *K, gsl_matrix *S, gsl_matrix *Svar, gsl_vector *ns);
-	void WriteVector (const gsl_vector *q, const gsl_vector *s, const size_t n_total, const string suffix);
+	void CalcS (const map<string, double> &mapRS2wA,
+		    const map<string, double> &mapRS2wK,
+		    const gsl_matrix *W, gsl_matrix *A, gsl_matrix *K,
+		    gsl_matrix *S, gsl_matrix *Svar, gsl_vector *ns);
+	void WriteVector (const gsl_vector *q, const gsl_vector *s,
+			  const size_t n_total, const string suffix);
 	void WriteVar (const string suffix);
 	void WriteMatrix (const gsl_matrix *matrix_U, const string suffix);
 	void WriteVector (const gsl_vector *vector_D, const string suffix);
 	void CopyRead (gsl_vector *log_N);
-	void ObtainWeight (const set<string> &setSnps_beta, map<string, double> &mapRS2wK);
-	void UpdateWeight (const size_t pve_flag, const map<string, double> &mapRS2wK, const size_t ni_test, const gsl_vector *ns, map<string, double> &mapRS2wA);
-	void UpdateSNPnZ (const map<string, double> &mapRS2wA, const map<string, string> &mapRS2A1, const map<string, double> &mapRS2z, gsl_vector *w, gsl_vector *z, vector<size_t> &vec_cat);
+	void ObtainWeight (const set<string> &setSnps_beta, map<string,
+			   double> &mapRS2wK);
+	void UpdateWeight (const size_t pve_flag,
+			   const map<string,double> &mapRS2wK,
+			   const size_t ni_test, const gsl_vector *ns,
+			   map<string, double> &mapRS2wA);
+	void UpdateSNPnZ (const map<string, double> &mapRS2wA,
+			  const map<string, string> &mapRS2A1,
+			  const map<string, double> &mapRS2z,
+			  gsl_vector *w, gsl_vector *z,
+			  vector<size_t> &vec_cat);
 	void UpdateSNP (const map<string, double> &mapRS2wA);
 };