From e5d51f314214d780ef5fbe91b50d25fad8472d54 Mon Sep 17 00:00:00 2001
From: Peter Carbonetto
Date: Thu, 8 Jun 2017 22:26:31 -0500
Subject: Removed commented-out code from gemma.cpp.

---
 src/lmm.cpp | 1373 ++++++++++++++++++++++++++++++++---------------------------
 1 file changed, 734 insertions(+), 639 deletions(-)

(limited to 'src/lmm.cpp')

diff --git a/src/lmm.cpp b/src/lmm.cpp
index a707534..73a9232 100644
--- a/src/lmm.cpp
+++ b/src/lmm.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
@@ -16,8 +16,6 @@
     along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */
 
-
-
 #include <iostream>
 #include <fstream>
 #include <sstream>
@@ -25,6 +23,7 @@
 #include <iomanip>
 #include <cmath>
 #include <iostream>
+#include <assert.h>
 #include <stdio.h>
 #include <stdlib.h>
 #include <bitset>
@@ -34,8 +33,6 @@
 #include "gsl/gsl_matrix.h"
 #include "gsl/gsl_linalg.h"
 #include "gsl/gsl_blas.h"
-
-
 #include "gsl/gsl_cdf.h"
 #include "gsl/gsl_roots.h"
 #include "gsl/gsl_min.h"
@@ -45,22 +42,11 @@
 #include "eigenlib.h"
 #include "lapack.h"
 #include "gzstream.h"
-
-#ifdef FORCE_FLOAT
-#include "lmm_float.h"
-#else
 #include "lmm.h"
-#endif
-
 
 using namespace std;
 
-
-
-
-
-void LMM::CopyFromParam (PARAM &cPar)
-{
+void LMM::CopyFromParam (PARAM &cPar) {
 	a_mode=cPar.a_mode;
 	d_pace=cPar.d_pace;
 
@@ -69,7 +55,8 @@ void LMM::CopyFromParam (PARAM &cPar)
 	file_out=cPar.file_out;
 	path_out=cPar.path_out;
 	file_gene=cPar.file_gene;
-	// WJA added
+	
+	// WJA added.
 	file_oxford=cPar.file_oxford;
 
 	l_min=cPar.l_min;
@@ -97,9 +84,7 @@ void LMM::CopyFromParam (PARAM &cPar)
 	return;
 }
 
-
-void LMM::CopyToParam (PARAM &cPar)
-{
+void LMM::CopyToParam (PARAM &cPar) {
 	cPar.time_UtX=time_UtX;
 	cPar.time_opt=time_opt;
 
@@ -108,83 +93,120 @@ void LMM::CopyToParam (PARAM &cPar)
 	return;
 }
 
-
-
-void LMM::WriteFiles ()
-{
+void LMM::WriteFiles () {
 	string file_str;
 	file_str=path_out+"/"+file_out;
 	file_str+=".assoc.txt";
 
 	ofstream outfile (file_str.c_str(), ofstream::out);
-	if (!outfile) {cout<<"error writing file: "<<file_str.c_str()<<endl; return;}
+	if (!outfile) {
+	  cout<<"error writing file: "<<file_str.c_str()<<endl;
+	  return;
+	}
 
 	if (!file_gene.empty()) {
 		outfile<<"geneID"<<"\t";
 
 		if (a_mode==1) {
-			outfile<<"beta"<<"\t"<<"se"<<"\t"<<"l_remle"<<"\t"<<"p_wald"<<endl;
+			outfile<<"beta"<<"\t"<<"se"<<"\t"<<"l_remle"<<
+			  "\t"<<"p_wald"<<endl;
 		} else if (a_mode==2) {
 			outfile<<"l_mle"<<"\t"<<"p_lrt"<<endl;
 		} else if (a_mode==3) {
 			outfile<<"beta"<<"\t"<<"se"<<"\t"<<"p_score"<<endl;
 		} else if (a_mode==4) {
-			outfile<<"beta"<<"\t"<<"se"<<"\t"<<"l_remle"<<"\t"<<"l_mle"<<"\t"<<"p_wald"<<"\t"<<"p_lrt"<<"\t"<<"p_score"<<endl;
+			outfile<<"beta"<<"\t"<<"se"<<"\t"<<"l_remle"<<
+			  "\t"<<"l_mle"<<"\t"<<"p_wald"<<"\t"<<"p_lrt"<<
+			  "\t"<<"p_score"<<endl;
 		} else {}
 
 		for (vector<SUMSTAT>::size_type t=0; t<sumStat.size(); ++t) {
 			outfile<<snpInfo[t].rs_number<<"\t";
 
 			if (a_mode==1) {
-				outfile<<scientific<<setprecision(6)<<sumStat[t].beta<<"\t"<<sumStat[t].se<<"\t"<<sumStat[t].lambda_remle<<"\t"<<sumStat[t].p_wald <<endl;
+				outfile<<scientific<<setprecision(6)<<
+				  sumStat[t].beta<<"\t"<<sumStat[t].se<<"\t"<<
+				  sumStat[t].lambda_remle<<"\t"<<
+				  sumStat[t].p_wald <<endl;
 			} else if (a_mode==2) {
-				outfile<<scientific<<setprecision(6)<<sumStat[t].lambda_mle<<"\t"<<sumStat[t].p_lrt<<endl;
+				outfile<<scientific<<setprecision(6)<<
+				  sumStat[t].lambda_mle<<"\t"<<
+				  sumStat[t].p_lrt<<endl;
 			} else if (a_mode==3) {
-				outfile<<scientific<<setprecision(6)<<sumStat[t].beta<<"\t"<<sumStat[t].se<<"\t"<<sumStat[t].p_score<<endl;
+				outfile<<scientific<<setprecision(6)<<
+				  sumStat[t].beta<<"\t"<<sumStat[t].se<<
+				  "\t"<<sumStat[t].p_score<<endl;
 			} else if (a_mode==4) {
-				outfile<<scientific<<setprecision(6)<<sumStat[t].beta<<"\t"<<sumStat[t].se<<"\t"<<sumStat[t].lambda_remle<<"\t"<<sumStat[t].lambda_mle<<"\t"<<sumStat[t].p_wald <<"\t"<<sumStat[t].p_lrt<<"\t"<<sumStat[t].p_score<<endl;
+				outfile<<scientific<<setprecision(6)<<
+				  sumStat[t].beta<<"\t"<<sumStat[t].se<<"\t"<<
+				  sumStat[t].lambda_remle<<"\t"<<
+				  sumStat[t].lambda_mle<<"\t"<<
+				  sumStat[t].p_wald <<"\t"<<
+				  sumStat[t].p_lrt<<"\t"<<
+				  sumStat[t].p_score<<endl;
 			} else {}
 		}
 	}  else {
-		outfile<<"chr"<<"\t"<<"rs"<<"\t"<<"ps"<<"\t"<<"n_miss"<<"\t"<<"allele1"<<"\t"<<"allele0"<<"\t"<<"af"<<"\t";
+		outfile<<"chr"<<"\t"<<"rs"<<"\t"<<"ps"<<"\t"<<"n_miss"<<"\t"
+		       <<"allele1"<<"\t"<<"allele0"<<"\t"<<"af"<<"\t";
 
 		if (a_mode==1) {
-			outfile<<"beta"<<"\t"<<"se"<<"\t"<<"l_remle"<<"\t"<<"p_wald"<<endl;
+			outfile<<"beta"<<"\t"<<"se"<<"\t"<<"l_remle"<<"\t"
+			       <<"p_wald"<<endl;
 		} else if (a_mode==2) {
 			outfile<<"l_mle"<<"\t"<<"p_lrt"<<endl;
 		} else if (a_mode==3) {
 			outfile<<"beta"<<"\t"<<"se"<<"\t"<<"p_score"<<endl;
 		} else if (a_mode==4) {
-			outfile<<"beta"<<"\t"<<"se"<<"\t"<<"l_remle"<<"\t"<<"l_mle"<<"\t"<<"p_wald"<<"\t"<<"p_lrt"<<"\t"<<"p_score"<<endl;
+			outfile<<"beta"<<"\t"<<"se"<<"\t"<<"l_remle"<<"\t"
+			       <<"l_mle"<<"\t"<<"p_wald"<<"\t"<<"p_lrt"<<
+			  "\t"<<"p_score"<<endl;
 		} else {}
 
 		size_t t=0;
 		for (size_t i=0; i<snpInfo.size(); ++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"<<snpInfo[i].a_minor<<"\t"<<snpInfo[i].a_major<<"\t"<<fixed<<setprecision(3)<<snpInfo[i].maf<<"\t";
+			outfile<<snpInfo[i].chr<<"\t"<<snpInfo[i].rs_number<<
+			  "\t"<<snpInfo[i].base_position<<"\t"<<
+			  snpInfo[i].n_miss<<"\t"<<snpInfo[i].a_minor<<"\t"<<
+			  snpInfo[i].a_major<<"\t"<<fixed<<setprecision(3)<<
+			  snpInfo[i].maf<<"\t";
 
 			if (a_mode==1) {
-				outfile<<scientific<<setprecision(6)<<sumStat[t].beta<<"\t"<<sumStat[t].se<<"\t"<<sumStat[t].lambda_remle<<"\t"<<sumStat[t].p_wald <<endl;
+				outfile<<scientific<<setprecision(6)<<
+				  sumStat[t].beta<<"\t"<<sumStat[t].se<<
+				  "\t"<<sumStat[t].lambda_remle<<"\t"<<
+				  sumStat[t].p_wald <<endl;
 			} else if (a_mode==2) {
-				outfile<<scientific<<setprecision(6)<<sumStat[t].lambda_mle<<"\t"<<sumStat[t].p_lrt<<endl;
+				outfile<<scientific<<setprecision(6)<<
+				  sumStat[t].lambda_mle<<"\t"<<
+				  sumStat[t].p_lrt<<endl;
 			} else if (a_mode==3) {
-				outfile<<scientific<<setprecision(6)<<sumStat[t].beta<<"\t"<<sumStat[t].se<<"\t"<<sumStat[t].p_score<<endl;
+				outfile<<scientific<<setprecision(6)<<
+				  sumStat[t].beta<<"\t"<<sumStat[t].se<<
+				  "\t"<<sumStat[t].p_score<<endl;
 			} else if (a_mode==4) {
-				outfile<<scientific<<setprecision(6)<<sumStat[t].beta<<"\t"<<sumStat[t].se<<"\t"<<sumStat[t].lambda_remle<<"\t"<<sumStat[t].lambda_mle<<"\t"<<sumStat[t].p_wald <<"\t"<<sumStat[t].p_lrt<<"\t"<<sumStat[t].p_score<<endl;
+				outfile<<scientific<<setprecision(6)<<
+				  sumStat[t].beta<<"\t"<<sumStat[t].se<<
+				  "\t"<<sumStat[t].lambda_remle<<"\t"<<
+				  sumStat[t].lambda_mle<<"\t"<<
+				  sumStat[t].p_wald <<"\t"<<
+				  sumStat[t].p_lrt<<"\t"<<
+				  sumStat[t].p_score<<endl;
 			} else {}
 			t++;
 		}
 	}
 
-
 	outfile.close();
 	outfile.clear();
 	return;
 }
 
-void CalcPab (const size_t n_cvt, const size_t e_mode, const gsl_vector *Hi_eval, const gsl_matrix *Uab, const gsl_vector *ab, gsl_matrix *Pab)
-{
+void CalcPab (const size_t n_cvt, const size_t e_mode,
+	      const gsl_vector *Hi_eval, const gsl_matrix *Uab,
+	      const gsl_vector *ab, gsl_matrix *Pab) {
 	size_t index_ab, index_aw, index_bw, index_ww;
 	double p_ab;
 	double ps_ab, ps_aw, ps_bw, ps_ww;
@@ -194,23 +216,26 @@ void CalcPab (const size_t n_cvt, const size_t e_mode, const gsl_vector *Hi_eval
 			for (size_t b=a; b<=n_cvt+2; ++b) {
 				index_ab=GetabIndex (a, b, n_cvt);
 				if (p==0) {
-					gsl_vector_const_view Uab_col=gsl_matrix_const_column (Uab, index_ab);
-					gsl_blas_ddot (Hi_eval, &Uab_col.vector, &p_ab);
-					if (e_mode!=0) {p_ab=gsl_vector_get (ab, index_ab)-p_ab;}
-					gsl_matrix_set (Pab, 0, index_ab, p_ab);
+				  gsl_vector_const_view Uab_col=
+				    gsl_matrix_const_column (Uab, index_ab);
+				  gsl_blas_ddot(Hi_eval,&Uab_col.vector,&p_ab);
+				  if (e_mode!=0) {
+				    p_ab=gsl_vector_get (ab, index_ab)-p_ab;
+				  }
+				  gsl_matrix_set (Pab, 0, index_ab, p_ab);
 				}
 				else {
-					index_aw=GetabIndex (a, p, n_cvt);
-					index_bw=GetabIndex (b, p, n_cvt);
-					index_ww=GetabIndex (p, p, n_cvt);
-
-					ps_ab=gsl_matrix_get (Pab, p-1, index_ab);
-					ps_aw=gsl_matrix_get (Pab, p-1, index_aw);
-					ps_bw=gsl_matrix_get (Pab, p-1, index_bw);
-					ps_ww=gsl_matrix_get (Pab, p-1, index_ww);
-
-					p_ab=ps_ab-ps_aw*ps_bw/ps_ww;
-					gsl_matrix_set (Pab, p, index_ab, p_ab);
+				  index_aw=GetabIndex (a, p, n_cvt);
+				  index_bw=GetabIndex (b, p, n_cvt);
+				  index_ww=GetabIndex (p, p, n_cvt);
+				  
+				  ps_ab=gsl_matrix_get (Pab, p-1, index_ab);
+				  ps_aw=gsl_matrix_get (Pab, p-1, index_aw);
+				  ps_bw=gsl_matrix_get (Pab, p-1, index_bw);
+				  ps_ww=gsl_matrix_get (Pab, p-1, index_ww);
+				  
+				  p_ab=ps_ab-ps_aw*ps_bw/ps_ww;
+				  gsl_matrix_set (Pab, p, index_ab, p_ab);
 				}
 			}
 		}
@@ -218,9 +243,9 @@ void CalcPab (const size_t n_cvt, const size_t e_mode, const gsl_vector *Hi_eval
 	return;
 }
 
-
-void CalcPPab (const size_t n_cvt, const size_t e_mode, const gsl_vector *HiHi_eval, const gsl_matrix *Uab, const gsl_vector *ab, const gsl_matrix *Pab, gsl_matrix *PPab)
-{
+void CalcPPab (const size_t n_cvt, const size_t e_mode,
+	       const gsl_vector *HiHi_eval, const gsl_matrix *Uab,
+	       const gsl_vector *ab, const gsl_matrix *Pab, gsl_matrix *PPab) {
 	size_t index_ab, index_aw, index_bw, index_ww;
 	double p2_ab;
 	double ps2_ab, ps_aw, ps_bw, ps_ww, ps2_aw, ps2_bw, ps2_ww;
@@ -230,28 +255,33 @@ void CalcPPab (const size_t n_cvt, const size_t e_mode, const gsl_vector *HiHi_e
 			for (size_t b=a; b<=n_cvt+2; ++b) {
 				index_ab=GetabIndex (a, b, n_cvt);
 				if (p==0) {
-					gsl_vector_const_view Uab_col=gsl_matrix_const_column (Uab, index_ab);
-					gsl_blas_ddot (HiHi_eval, &Uab_col.vector, &p2_ab);
-					if (e_mode!=0) {p2_ab=p2_ab-gsl_vector_get (ab, index_ab)+2.0*gsl_matrix_get (Pab, 0, index_ab);}
-					gsl_matrix_set (PPab, 0, index_ab, p2_ab);
+				  gsl_vector_const_view Uab_col=
+				    gsl_matrix_const_column (Uab, index_ab);
+				  gsl_blas_ddot (HiHi_eval, &Uab_col.vector,
+						 &p2_ab);
+				  if (e_mode!=0) {
+				    p2_ab=p2_ab-gsl_vector_get(ab,index_ab) +
+				      2.0*gsl_matrix_get (Pab, 0, index_ab);
+				  }
+				  gsl_matrix_set (PPab, 0, index_ab, p2_ab);
 				}
 				else {
-					index_aw=GetabIndex (a, p, n_cvt);
-					index_bw=GetabIndex (b, p, n_cvt);
-					index_ww=GetabIndex (p, p, n_cvt);
-
-					ps2_ab=gsl_matrix_get (PPab, p-1, index_ab);
-					ps_aw=gsl_matrix_get (Pab, p-1, index_aw);
-					ps_bw=gsl_matrix_get (Pab, p-1, index_bw);
-					ps_ww=gsl_matrix_get (Pab, p-1, index_ww);
-					ps2_aw=gsl_matrix_get (PPab, p-1, index_aw);
-					ps2_bw=gsl_matrix_get (PPab, p-1, index_bw);
-					ps2_ww=gsl_matrix_get (PPab, p-1, index_ww);
-
-					p2_ab=ps2_ab+ps_aw*ps_bw*ps2_ww/(ps_ww*ps_ww);
-					p2_ab-=(ps_aw*ps2_bw+ps_bw*ps2_aw)/ps_ww;
-					gsl_matrix_set (PPab, p, index_ab, p2_ab);
-
+				  index_aw=GetabIndex (a, p, n_cvt);
+				  index_bw=GetabIndex (b, p, n_cvt);
+				  index_ww=GetabIndex (p, p, n_cvt);
+				  
+				  ps2_ab=gsl_matrix_get (PPab, p-1, index_ab);
+				  ps_aw=gsl_matrix_get (Pab, p-1, index_aw);
+				  ps_bw=gsl_matrix_get (Pab, p-1, index_bw);
+				  ps_ww=gsl_matrix_get (Pab, p-1, index_ww);
+				  ps2_aw=gsl_matrix_get (PPab, p-1, index_aw);
+				  ps2_bw=gsl_matrix_get (PPab, p-1, index_bw);
+				  ps2_ww=gsl_matrix_get (PPab, p-1, index_ww);
+				  
+				  p2_ab=ps2_ab+ps_aw*ps_bw*
+				    ps2_ww/(ps_ww*ps_ww);
+				  p2_ab-=(ps_aw*ps2_bw+ps_bw*ps2_aw)/ps_ww;
+				  gsl_matrix_set (PPab, p, index_ab, p2_ab);
 				}
 			}
 		}
@@ -259,44 +289,56 @@ void CalcPPab (const size_t n_cvt, const size_t e_mode, const gsl_vector *HiHi_e
 	return;
 }
 
-
-void CalcPPPab (const size_t n_cvt, const size_t e_mode, const gsl_vector *HiHiHi_eval, const gsl_matrix *Uab, const gsl_vector *ab, const gsl_matrix *Pab, const gsl_matrix *PPab, gsl_matrix *PPPab)
-{
+void CalcPPPab (const size_t n_cvt, const size_t e_mode,
+		const gsl_vector *HiHiHi_eval, const gsl_matrix *Uab,
+		const gsl_vector *ab, const gsl_matrix *Pab,
+		const gsl_matrix *PPab, gsl_matrix *PPPab) {
 	size_t index_ab, index_aw, index_bw, index_ww;
 	double p3_ab;
-	double ps3_ab, ps_aw, ps_bw, ps_ww, ps2_aw, ps2_bw, ps2_ww, ps3_aw, ps3_bw, ps3_ww;
+	double ps3_ab, ps_aw, ps_bw, ps_ww, ps2_aw, ps2_bw, ps2_ww,
+	  ps3_aw, ps3_bw, ps3_ww;
 
 	for (size_t p=0; p<=n_cvt+1; ++p) {
 		for (size_t a=p+1; a<=n_cvt+2; ++a) {
 			for (size_t b=a; b<=n_cvt+2; ++b) {
 				index_ab=GetabIndex (a, b, n_cvt);
 				if (p==0) {
-					gsl_vector_const_view Uab_col=gsl_matrix_const_column (Uab, index_ab);
-					gsl_blas_ddot (HiHiHi_eval, &Uab_col.vector, &p3_ab);
-					if (e_mode!=0) {p3_ab=gsl_vector_get (ab, index_ab)-p3_ab+3.0*gsl_matrix_get (PPab, 0, index_ab)-3.0*gsl_matrix_get (Pab, 0, index_ab);}
-					gsl_matrix_set (PPPab, 0, index_ab, p3_ab);
+				  gsl_vector_const_view Uab_col=
+				    gsl_matrix_const_column (Uab, index_ab);
+				  gsl_blas_ddot (HiHiHi_eval, &Uab_col.vector,
+						 &p3_ab);
+				  if (e_mode!=0) {
+				    p3_ab=gsl_vector_get (ab, index_ab)-
+				      p3_ab+3.0*gsl_matrix_get(PPab,0,index_ab)
+				      -3.0*gsl_matrix_get (Pab, 0, index_ab);
+				  }
+				  gsl_matrix_set (PPPab, 0, index_ab, p3_ab);
 				}
 				else {
-					index_aw=GetabIndex (a, p, n_cvt);
-					index_bw=GetabIndex (b, p, n_cvt);
-					index_ww=GetabIndex (p, p, n_cvt);
-
-					ps3_ab=gsl_matrix_get (PPPab, p-1, index_ab);
-					ps_aw=gsl_matrix_get (Pab, p-1, index_aw);
-					ps_bw=gsl_matrix_get (Pab, p-1, index_bw);
-					ps_ww=gsl_matrix_get (Pab, p-1, index_ww);
-					ps2_aw=gsl_matrix_get (PPab, p-1, index_aw);
-					ps2_bw=gsl_matrix_get (PPab, p-1, index_bw);
-					ps2_ww=gsl_matrix_get (PPab, p-1, index_ww);
-					ps3_aw=gsl_matrix_get (PPPab, p-1, index_aw);
-					ps3_bw=gsl_matrix_get (PPPab, p-1, index_bw);
-					ps3_ww=gsl_matrix_get (PPPab, p-1, index_ww);
-
-					p3_ab=ps3_ab-ps_aw*ps_bw*ps2_ww*ps2_ww/(ps_ww*ps_ww*ps_ww);
-					p3_ab-=(ps_aw*ps3_bw+ps_bw*ps3_aw+ps2_aw*ps2_bw)/ps_ww;
-					p3_ab+=(ps_aw*ps2_bw*ps2_ww+ps_bw*ps2_aw*ps2_ww+ps_aw*ps_bw*ps3_ww)/(ps_ww*ps_ww);
-
-					gsl_matrix_set (PPPab, p, index_ab, p3_ab);
+				  index_aw=GetabIndex (a, p, n_cvt);
+				  index_bw=GetabIndex (b, p, n_cvt);
+				  index_ww=GetabIndex (p, p, n_cvt);
+				  
+				  ps3_ab=gsl_matrix_get (PPPab, p-1, index_ab);
+				  ps_aw=gsl_matrix_get (Pab, p-1, index_aw);
+				  ps_bw=gsl_matrix_get (Pab, p-1, index_bw);
+				  ps_ww=gsl_matrix_get (Pab, p-1, index_ww);
+				  ps2_aw=gsl_matrix_get (PPab, p-1, index_aw);
+				  ps2_bw=gsl_matrix_get (PPab, p-1, index_bw);
+				  ps2_ww=gsl_matrix_get (PPab, p-1, index_ww);
+				  ps3_aw=gsl_matrix_get (PPPab, p-1, index_aw);
+				  ps3_bw=gsl_matrix_get (PPPab, p-1, index_bw);
+				  ps3_ww=gsl_matrix_get (PPPab, p-1, index_ww);
+
+				  p3_ab=ps3_ab-ps_aw*ps_bw*ps2_ww*ps2_ww
+				    /(ps_ww*ps_ww*ps_ww);
+				  p3_ab-=(ps_aw*ps3_bw+ps_bw*ps3_aw +
+					  ps2_aw*ps2_bw)/ps_ww;
+				  p3_ab+=(ps_aw*ps2_bw*ps2_ww+ps_bw*
+					  ps2_aw*ps2_ww+ps_aw*ps_bw*ps3_ww)/
+				    (ps_ww*ps_ww);
+				  
+				  gsl_matrix_set (PPPab, p, index_ab, p3_ab);
 				}
 			}
 		}
@@ -304,10 +346,7 @@ void CalcPPPab (const size_t n_cvt, const size_t e_mode, const gsl_vector *HiHiH
 	return;
 }
 
-
-
-double LogL_f (double l, void *params)
-{
+double LogL_f (double l, void *params) {
 	FUNC_PARAM *p=(FUNC_PARAM *) params;
 	size_t n_cvt=p->n_cvt;
 	size_t ni_test=p->ni_test;
@@ -325,7 +364,11 @@ double LogL_f (double l, void *params)
 
 	gsl_vector_memcpy (v_temp, p->eval);
 	gsl_vector_scale (v_temp, l);
-	if (p->e_mode==0) {gsl_vector_set_all (Hi_eval, 1.0);} else {gsl_vector_memcpy (Hi_eval, v_temp);}
+	if (p->e_mode==0) {
+	  gsl_vector_set_all (Hi_eval, 1.0);
+	} else {
+	  gsl_vector_memcpy (Hi_eval, v_temp);
+	}
 	gsl_vector_add_constant (v_temp, 1.0);
 	gsl_vector_div (Hi_eval, v_temp);
 
@@ -348,13 +391,7 @@ double LogL_f (double l, void *params)
 	return f;
 }
 
-
-
-
-
-
-double LogL_dev1 (double l, void *params)
-{
+double LogL_dev1 (double l, void *params) {
 	FUNC_PARAM *p=(FUNC_PARAM *) params;
 	size_t n_cvt=p->n_cvt;
 	size_t ni_test=p->ni_test;
@@ -374,7 +411,11 @@ double LogL_dev1 (double l, void *params)
 
 	gsl_vector_memcpy (v_temp, p->eval);
 	gsl_vector_scale (v_temp, l);
-	if (p->e_mode==0) {gsl_vector_set_all (Hi_eval, 1.0);} else {gsl_vector_memcpy (Hi_eval, v_temp);}
+	if (p->e_mode==0) {
+	  gsl_vector_set_all (Hi_eval, 1.0);
+	} else {
+	  gsl_vector_memcpy (Hi_eval, v_temp);
+	}
 	gsl_vector_add_constant (v_temp, 1.0);
 	gsl_vector_div (Hi_eval, v_temp);
 
@@ -407,18 +448,18 @@ double LogL_dev1 (double l, void *params)
 	return dev1;
 }
 
-
-
-
-double LogL_dev2 (double l, void *params)
-{
+double LogL_dev2 (double l, void *params) {
 	FUNC_PARAM *p=(FUNC_PARAM *) params;
 	size_t n_cvt=p->n_cvt;
 	size_t ni_test=p->ni_test;
 	size_t n_index=(n_cvt+2+1)*(n_cvt+2)/2;
 
 	size_t nc_total;
-	if (p->calc_null==true) {nc_total=n_cvt;} else {nc_total=n_cvt+1;}
+	if (p->calc_null==true) {
+	  nc_total=n_cvt;
+	} else {
+	  nc_total=n_cvt+1;
+	}
 
 	double dev2=0.0, trace_Hi=0.0, trace_HiHi=0.0;
 	size_t index_yy;
@@ -433,7 +474,11 @@ double LogL_dev2 (double l, void *params)
 
 	gsl_vector_memcpy (v_temp, p->eval);
 	gsl_vector_scale (v_temp, l);
-	if (p->e_mode==0) {gsl_vector_set_all (Hi_eval, 1.0);} else {gsl_vector_memcpy (Hi_eval, v_temp);}
+	if (p->e_mode==0) {
+	  gsl_vector_set_all (Hi_eval, 1.0);
+	} else {
+	  gsl_vector_memcpy (Hi_eval, v_temp);
+	}
 	gsl_vector_add_constant (v_temp, 1.0);
 	gsl_vector_div (Hi_eval, v_temp);
 
@@ -453,7 +498,8 @@ double LogL_dev2 (double l, void *params)
 
 	CalcPab (n_cvt, p->e_mode, Hi_eval, p->Uab, p->ab, Pab);
 	CalcPPab (n_cvt, p->e_mode, HiHi_eval, p->Uab, p->ab, Pab, PPab);
-	CalcPPPab (n_cvt, p->e_mode, HiHiHi_eval, p->Uab, p->ab, Pab, PPab, PPPab);
+	CalcPPPab (n_cvt, p->e_mode, HiHiHi_eval, p->Uab, p->ab, Pab, PPab,
+		   PPPab);
 
 	double trace_HiKHiK=((double)ni_test+trace_HiHi-2*trace_Hi)/(l*l);
 
@@ -465,7 +511,8 @@ double LogL_dev2 (double l, void *params)
 	double yPKPy=(P_yy-PP_yy)/l;
 	double yPKPKPy=(P_yy+PPP_yy-2.0*PP_yy)/(l*l);
 
-	dev2=0.5*trace_HiKHiK-0.5*(double)ni_test*(2.0*yPKPKPy*P_yy-yPKPy*yPKPy)/(P_yy*P_yy);
+	dev2=0.5*trace_HiKHiK-0.5*(double)ni_test*
+	  (2.0*yPKPKPy*P_yy-yPKPy*yPKPy)/(P_yy*P_yy);
 
 	gsl_matrix_free (Pab);
 	gsl_matrix_free (PPab);
@@ -478,12 +525,7 @@ double LogL_dev2 (double l, void *params)
 	return dev2;
 }
 
-
-
-
-
-void LogL_dev12 (double l, void *params, double *dev1, double *dev2)
-{
+void LogL_dev12 (double l, void *params, double *dev1, double *dev2) {
 	FUNC_PARAM *p=(FUNC_PARAM *) params;
 	size_t n_cvt=p->n_cvt;
 	size_t ni_test=p->ni_test;
@@ -505,7 +547,11 @@ void LogL_dev12 (double l, void *params, double *dev1, double *dev2)
 
 	gsl_vector_memcpy (v_temp, p->eval);
 	gsl_vector_scale (v_temp, l);
-	if (p->e_mode==0) {gsl_vector_set_all (Hi_eval, 1.0);} else {gsl_vector_memcpy (Hi_eval, v_temp);}
+	if (p->e_mode==0) {
+	  gsl_vector_set_all (Hi_eval, 1.0);
+	} else {
+	  gsl_vector_memcpy (Hi_eval, v_temp);
+	}
 	gsl_vector_add_constant (v_temp, 1.0);
 	gsl_vector_div (Hi_eval, v_temp);
 
@@ -525,7 +571,8 @@ void LogL_dev12 (double l, void *params, double *dev1, double *dev2)
 
 	CalcPab (n_cvt, p->e_mode, Hi_eval, p->Uab, p->ab, Pab);
 	CalcPPab (n_cvt, p->e_mode, HiHi_eval, p->Uab, p->ab, Pab, PPab);
-	CalcPPPab (n_cvt, p->e_mode, HiHiHi_eval, p->Uab, p->ab, Pab, PPab, PPPab);
+	CalcPPPab (n_cvt, p->e_mode, HiHiHi_eval, p->Uab, p->ab, Pab, PPab,
+		   PPPab);
 
 	double trace_HiK=((double)ni_test-trace_Hi)/l;
 	double trace_HiKHiK=((double)ni_test+trace_HiHi-2*trace_Hi)/(l*l);
@@ -540,7 +587,8 @@ void LogL_dev12 (double l, void *params, double *dev1, double *dev2)
 	double yPKPKPy=(P_yy+PPP_yy-2.0*PP_yy)/(l*l);
 
 	*dev1=-0.5*trace_HiK+0.5*(double)ni_test*yPKPy/P_yy;
-	*dev2=0.5*trace_HiKHiK-0.5*(double)ni_test*(2.0*yPKPKPy*P_yy-yPKPy*yPKPy)/(P_yy*P_yy);
+	*dev2=0.5*trace_HiKHiK-0.5*(double)ni_test*
+	  (2.0*yPKPKPy*P_yy-yPKPy*yPKPy)/(P_yy*P_yy);
 
 	gsl_matrix_free (Pab);
 	gsl_matrix_free (PPab);
@@ -553,10 +601,7 @@ void LogL_dev12 (double l, void *params, double *dev1, double *dev2)
 	return;
 }
 
-
-
-double LogRL_f (double l, void *params)
-{
+double LogRL_f (double l, void *params) {
 	FUNC_PARAM *p=(FUNC_PARAM *) params;
 	size_t n_cvt=p->n_cvt;
 	size_t ni_test=p->ni_test;
@@ -564,7 +609,9 @@ double LogRL_f (double l, void *params)
 
 	double df;
 	size_t nc_total;
-	if (p->calc_null==true) {nc_total=n_cvt; df=(double)ni_test-(double)n_cvt; }
+	if (p->calc_null==true) {
+	  nc_total=n_cvt; df=(double)ni_test-(double)n_cvt;
+	}
 	else {nc_total=n_cvt+1; df=(double)ni_test-(double)n_cvt-1.0;}
 
 	double f=0.0, logdet_h=0.0, logdet_hiw=0.0, d;
@@ -577,7 +624,11 @@ double LogRL_f (double l, void *params)
 
 	gsl_vector_memcpy (v_temp, p->eval);
 	gsl_vector_scale (v_temp, l);
-	if (p->e_mode==0) {gsl_vector_set_all (Hi_eval, 1.0);} else {gsl_vector_memcpy (Hi_eval, v_temp);}
+	if (p->e_mode==0) {
+	  gsl_vector_set_all (Hi_eval, 1.0);
+	} else {
+	  gsl_vector_memcpy (Hi_eval, v_temp);
+	}
 	gsl_vector_add_constant (v_temp, 1.0);
 	gsl_vector_div (Hi_eval, v_temp);
 
@@ -590,7 +641,7 @@ double LogRL_f (double l, void *params)
 	gsl_vector_set_all (v_temp, 1.0);
 	CalcPab (n_cvt, p->e_mode, v_temp, p->Uab, p->ab, Iab);
 
-	//calculate |WHiW|-|WW|
+	// Calculate |WHiW|-|WW|.
 	logdet_hiw=0.0;
 	for (size_t i=0; i<nc_total; ++i) {
 		index_ww=GetabIndex (i+1, i+1, n_cvt);
@@ -612,10 +663,7 @@ double LogRL_f (double l, void *params)
 	return f;
 }
 
-
-
-double LogRL_dev1 (double l, void *params)
-{
+double LogRL_dev1 (double l, void *params) {
 	FUNC_PARAM *p=(FUNC_PARAM *) params;
 	size_t n_cvt=p->n_cvt;
 	size_t ni_test=p->ni_test;
@@ -623,8 +671,14 @@ double LogRL_dev1 (double l, void *params)
 
 	double df;
 	size_t nc_total;
-	if (p->calc_null==true) {nc_total=n_cvt; df=(double)ni_test-(double)n_cvt; }
-	else {nc_total=n_cvt+1; df=(double)ni_test-(double)n_cvt-1.0;}
+	if (p->calc_null==true) {
+	  nc_total=n_cvt;
+	  df=(double)ni_test-(double)n_cvt;
+	}
+	else {
+	  nc_total=n_cvt+1;
+	  df=(double)ni_test-(double)n_cvt-1.0;
+	}
 
 	double dev1=0.0, trace_Hi=0.0;
 	size_t index_ww;
@@ -637,7 +691,11 @@ double LogRL_dev1 (double l, void *params)
 
 	gsl_vector_memcpy (v_temp, p->eval);
 	gsl_vector_scale (v_temp, l);
-	if (p->e_mode==0) {gsl_vector_set_all (Hi_eval, 1.0);} else {gsl_vector_memcpy (Hi_eval, v_temp);}
+	if (p->e_mode==0) {
+	  gsl_vector_set_all (Hi_eval, 1.0);
+	} else {
+	  gsl_vector_memcpy (Hi_eval, v_temp);
+	}
 	gsl_vector_add_constant (v_temp, 1.0);
 	gsl_vector_div (Hi_eval, v_temp);
 
@@ -654,7 +712,7 @@ double LogRL_dev1 (double l, void *params)
 	CalcPab (n_cvt, p->e_mode, Hi_eval, p->Uab, p->ab, Pab);
 	CalcPPab (n_cvt, p->e_mode, HiHi_eval, p->Uab, p->ab, Pab, PPab);
 
-	//calculate tracePK and trace PKPK
+	// Calculate tracePK and trace PKPK.
 	double trace_P=trace_Hi;
 	double ps_ww, ps2_ww;
 	for (size_t i=0; i<nc_total; ++i) {
@@ -665,7 +723,7 @@ double LogRL_dev1 (double l, void *params)
 	}
 	double trace_PK=(df-trace_P)/l;
 
-	//calculate yPKPy, yPKPKPy
+	// Calculate yPKPy, yPKPKPy.
 	index_ww=GetabIndex (n_cvt+2, n_cvt+2, n_cvt);
 	double P_yy=gsl_matrix_get (Pab, nc_total, index_ww);
 	double PP_yy=gsl_matrix_get (PPab, nc_total, index_ww);
@@ -682,11 +740,7 @@ double LogRL_dev1 (double l, void *params)
 	return dev1;
 }
 
-
-
-
-double LogRL_dev2 (double l, void *params)
-{
+double LogRL_dev2 (double l, void *params) {
 	FUNC_PARAM *p=(FUNC_PARAM *) params;
 	size_t n_cvt=p->n_cvt;
 	size_t ni_test=p->ni_test;
@@ -694,8 +748,14 @@ double LogRL_dev2 (double l, void *params)
 
 	double df;
 	size_t nc_total;
-	if (p->calc_null==true) {nc_total=n_cvt; df=(double)ni_test-(double)n_cvt; }
-	else {nc_total=n_cvt+1; df=(double)ni_test-(double)n_cvt-1.0;}
+	if (p->calc_null==true) {
+	  nc_total=n_cvt;
+	  df=(double)ni_test-(double)n_cvt;
+	}
+	else {
+	  nc_total=n_cvt+1;
+	  df=(double)ni_test-(double)n_cvt-1.0;
+	}
 
 	double dev2=0.0, trace_Hi=0.0, trace_HiHi=0.0;
 	size_t index_ww;
@@ -710,7 +770,11 @@ double LogRL_dev2 (double l, void *params)
 
 	gsl_vector_memcpy (v_temp, p->eval);
 	gsl_vector_scale (v_temp, l);
-	if (p->e_mode==0) {gsl_vector_set_all (Hi_eval, 1.0);} else {gsl_vector_memcpy (Hi_eval, v_temp);}
+	if (p->e_mode==0) {
+	  gsl_vector_set_all (Hi_eval, 1.0);
+	} else {
+	  gsl_vector_memcpy (Hi_eval, v_temp);
+	}
 	gsl_vector_add_constant (v_temp, 1.0);
 	gsl_vector_div (Hi_eval, v_temp);
 
@@ -730,9 +794,10 @@ double LogRL_dev2 (double l, void *params)
 
 	CalcPab (n_cvt, p->e_mode, Hi_eval, p->Uab, p->ab, Pab);
 	CalcPPab (n_cvt, p->e_mode, HiHi_eval, p->Uab, p->ab, Pab, PPab);
-	CalcPPPab (n_cvt, p->e_mode, HiHiHi_eval, p->Uab, p->ab, Pab, PPab, PPPab);
+	CalcPPPab (n_cvt, p->e_mode, HiHiHi_eval, p->Uab, p->ab, Pab,
+		   PPab, PPPab);
 
-	//calculate tracePK and trace PKPK
+	// Calculate tracePK and trace PKPK.
 	double trace_P=trace_Hi, trace_PP=trace_HiHi;
 	double ps_ww, ps2_ww, ps3_ww;
 	for (size_t i=0; i<nc_total; ++i) {
@@ -745,7 +810,7 @@ double LogRL_dev2 (double l, void *params)
 	}
 	double trace_PKPK=(df+trace_PP-2.0*trace_P)/(l*l);
 
-	//calculate yPKPy, yPKPKPy
+	// Calculate yPKPy, yPKPKPy.
 	index_ww=GetabIndex (n_cvt+2, n_cvt+2, n_cvt);
 	double P_yy=gsl_matrix_get (Pab, nc_total, index_ww);
 	double PP_yy=gsl_matrix_get (PPab, nc_total, index_ww);
@@ -766,11 +831,7 @@ double LogRL_dev2 (double l, void *params)
 	return dev2;
 }
 
-
-
-
-void LogRL_dev12 (double l, void *params, double *dev1, double *dev2)
-{
+void LogRL_dev12 (double l, void *params, double *dev1, double *dev2) {
 	FUNC_PARAM *p=(FUNC_PARAM *) params;
 	size_t n_cvt=p->n_cvt;
 	size_t ni_test=p->ni_test;
@@ -778,8 +839,14 @@ void LogRL_dev12 (double l, void *params, double *dev1, double *dev2)
 
 	double df;
 	size_t nc_total;
-	if (p->calc_null==true) {nc_total=n_cvt; df=(double)ni_test-(double)n_cvt; }
-	else {nc_total=n_cvt+1; df=(double)ni_test-(double)n_cvt-1.0;}
+	if (p->calc_null==true) {
+	  nc_total=n_cvt;
+	  df=(double)ni_test-(double)n_cvt;
+	}
+	else {
+	  nc_total=n_cvt+1;
+	  df=(double)ni_test-(double)n_cvt-1.0;
+	}
 
 	double trace_Hi=0.0, trace_HiHi=0.0;
 	size_t index_ww;
@@ -794,7 +861,11 @@ void LogRL_dev12 (double l, void *params, double *dev1, double *dev2)
 
 	gsl_vector_memcpy (v_temp, p->eval);
 	gsl_vector_scale (v_temp, l);
-	if (p->e_mode==0) {gsl_vector_set_all (Hi_eval, 1.0);} else {gsl_vector_memcpy (Hi_eval, v_temp);}
+	if (p->e_mode==0) {
+	  gsl_vector_set_all (Hi_eval, 1.0);
+	} else {
+	  gsl_vector_memcpy (Hi_eval, v_temp);
+	}
 	gsl_vector_add_constant (v_temp, 1.0);
 	gsl_vector_div (Hi_eval, v_temp);
 
@@ -814,9 +885,10 @@ void LogRL_dev12 (double l, void *params, double *dev1, double *dev2)
 
 	CalcPab (n_cvt, p->e_mode, Hi_eval, p->Uab, p->ab, Pab);
 	CalcPPab (n_cvt, p->e_mode, HiHi_eval, p->Uab, p->ab, Pab, PPab);
-	CalcPPPab (n_cvt, p->e_mode, HiHiHi_eval, p->Uab, p->ab, Pab, PPab, PPPab);
+	CalcPPPab (n_cvt, p->e_mode, HiHiHi_eval, p->Uab, p->ab, Pab,
+		   PPab, PPPab);
 
-	//calculate tracePK and trace PKPK
+	// Calculate tracePK and trace PKPK.
 	double trace_P=trace_Hi, trace_PP=trace_HiHi;
 	double ps_ww, ps2_ww, ps3_ww;
 	for (size_t i=0; i<nc_total; ++i) {
@@ -830,7 +902,7 @@ void LogRL_dev12 (double l, void *params, double *dev1, double *dev2)
 	double trace_PK=(df-trace_P)/l;
 	double trace_PKPK=(df+trace_PP-2.0*trace_P)/(l*l);
 
-	//calculate yPKPy, yPKPKPy
+	// Calculate yPKPy, yPKPKPy.
 	index_ww=GetabIndex (n_cvt+2, n_cvt+2, n_cvt);
 	double P_yy=gsl_matrix_get (Pab, nc_total, index_ww);
 	double PP_yy=gsl_matrix_get (PPab, nc_total, index_ww);
@@ -839,7 +911,8 @@ void LogRL_dev12 (double l, void *params, double *dev1, double *dev2)
 	double yPKPKPy=(P_yy+PPP_yy-2.0*PP_yy)/(l*l);
 
 	*dev1=-0.5*trace_PK+0.5*df*yPKPy/P_yy;
-	*dev2=0.5*trace_PKPK-0.5*df*(2.0*yPKPKPy*P_yy-yPKPy*yPKPy)/(P_yy*P_yy);
+	*dev2=0.5*trace_PKPK-0.5*df*(2.0*yPKPKPy*P_yy-yPKPy*yPKPy)/
+	  (P_yy*P_yy);
 
 	gsl_matrix_free (Pab);
 	gsl_matrix_free (PPab);
@@ -849,18 +922,11 @@ void LogRL_dev12 (double l, void *params, double *dev1, double *dev2)
 	gsl_vector_free (HiHiHi_eval);
 	gsl_vector_free (v_temp);
 
-	return ;
+	return;
 }
 
-
-
-
-
-
-
-
-void LMM::CalcRLWald (const double &l, const FUNC_PARAM &params, double &beta, double &se, double &p_wald)
-{
+void LMM::CalcRLWald (const double &l, const FUNC_PARAM &params,
+		      double &beta, double &se, double &p_wald) {
 	size_t n_cvt=params.n_cvt;
 	size_t n_index=(n_cvt+2+1)*(n_cvt+2)/2;
 
@@ -872,7 +938,11 @@ void LMM::CalcRLWald (const double &l, const FUNC_PARAM &params, double &beta, d
 
 	gsl_vector_memcpy (v_temp, params.eval);
 	gsl_vector_scale (v_temp, l);
-	if (params.e_mode==0) {gsl_vector_set_all (Hi_eval, 1.0);} else {gsl_vector_memcpy (Hi_eval, v_temp);}
+	if (params.e_mode==0) {
+	  gsl_vector_set_all (Hi_eval, 1.0);
+	} else {
+	  gsl_vector_memcpy (Hi_eval, v_temp);
+	}
 	gsl_vector_add_constant (v_temp, 1.0);
 	gsl_vector_div (Hi_eval, v_temp);
 
@@ -890,17 +960,15 @@ void LMM::CalcRLWald (const double &l, const FUNC_PARAM &params, double &beta, d
 	double tau=(double)df/Px_yy;
 	se=sqrt(1.0/(tau*P_xx));
 	p_wald=gsl_cdf_fdist_Q ((P_yy-Px_yy)*tau, 1.0, df);
-//	p_wald=gsl_cdf_chisq_Q ((P_yy-Px_yy)*tau, 1);
 
 	gsl_matrix_free (Pab);
 	gsl_vector_free (Hi_eval);
 	gsl_vector_free (v_temp);
-	return ;
+	return;
 }
 
-
-void LMM::CalcRLScore (const double &l, const FUNC_PARAM &params, double &beta, double &se, double &p_score)
-{
+void LMM::CalcRLScore (const double &l, const FUNC_PARAM &params,
+		       double &beta, double &se, double &p_score) {
 	size_t n_cvt=params.n_cvt;
 	size_t n_index=(n_cvt+2+1)*(n_cvt+2)/2;
 
@@ -912,7 +980,11 @@ void LMM::CalcRLScore (const double &l, const FUNC_PARAM &params, double &beta,
 
 	gsl_vector_memcpy (v_temp, params.eval);
 	gsl_vector_scale (v_temp, l);
-	if (params.e_mode==0) {gsl_vector_set_all (Hi_eval, 1.0);} else {gsl_vector_memcpy (Hi_eval, v_temp);}
+	if (params.e_mode==0) {
+	  gsl_vector_set_all (Hi_eval, 1.0);
+	} else {
+	  gsl_vector_memcpy (Hi_eval, v_temp);
+	}
 	gsl_vector_add_constant (v_temp, 1.0);
 	gsl_vector_div (Hi_eval, v_temp);
 
@@ -930,24 +1002,16 @@ void LMM::CalcRLScore (const double &l, const FUNC_PARAM &params, double &beta,
 	double tau=(double)df/Px_yy;
 	se=sqrt(1.0/(tau*P_xx));
 
-	p_score=gsl_cdf_fdist_Q ((double)ni_test*P_xy*P_xy/(P_yy*P_xx), 1.0, df);
-//	p_score=gsl_cdf_chisq_Q ((double)ni_test*P_xy*P_xy/(P_yy*P_xx), 1);
+	p_score=gsl_cdf_fdist_Q ((double)ni_test*P_xy*P_xy/(P_yy*P_xx),
+				 1.0, df);
 
 	gsl_matrix_free (Pab);
 	gsl_vector_free (Hi_eval);
 	gsl_vector_free (v_temp);
-	return ;
+	return;
 }
 
-
-
-
-
-
-
-
-void CalcUab (const gsl_matrix *UtW, const gsl_vector *Uty, gsl_matrix *Uab)
-{
+void CalcUab (const gsl_matrix *UtW, const gsl_vector *Uty, gsl_matrix *Uab) {
 	size_t index_ab;
 	size_t n_cvt=UtW->size2;
 
@@ -958,20 +1022,26 @@ void CalcUab (const gsl_matrix *UtW, const gsl_vector *Uty, gsl_matrix *Uab)
 
 		if (a==n_cvt+2) {gsl_vector_memcpy (u_a, Uty);}
 		else {
-			gsl_vector_const_view UtW_col=gsl_matrix_const_column (UtW, a-1);
-			gsl_vector_memcpy (u_a, &UtW_col.vector);
+		  gsl_vector_const_view UtW_col=
+		    gsl_matrix_const_column (UtW, a-1);
+		  gsl_vector_memcpy (u_a, &UtW_col.vector);
 		}
 
 		for (size_t b=a; b>=1; --b) {
 			if (b==n_cvt+1) {continue;}
 
 			index_ab=GetabIndex (a, b, n_cvt);
-			gsl_vector_view Uab_col=gsl_matrix_column (Uab, index_ab);
+			gsl_vector_view Uab_col=
+			  gsl_matrix_column (Uab, index_ab);
 
-			if (b==n_cvt+2) {gsl_vector_memcpy (&Uab_col.vector, Uty);}
+			if (b==n_cvt+2) {
+			  gsl_vector_memcpy (&Uab_col.vector, Uty);
+			}
 			else {
-				gsl_vector_const_view UtW_col=gsl_matrix_const_column (UtW, b-1);
-				gsl_vector_memcpy (&Uab_col.vector, &UtW_col.vector);
+				gsl_vector_const_view UtW_col=
+				  gsl_matrix_const_column (UtW, b-1);
+				gsl_vector_memcpy (&Uab_col.vector,
+						   &UtW_col.vector);
 			}
 
 			gsl_vector_mul(&Uab_col.vector, u_a);
@@ -982,9 +1052,8 @@ void CalcUab (const gsl_matrix *UtW, const gsl_vector *Uty, gsl_matrix *Uab)
 	return;
 }
 
-
-void CalcUab (const gsl_matrix *UtW, const gsl_vector *Uty, const gsl_vector *Utx, gsl_matrix *Uab)
-{
+void CalcUab (const gsl_matrix *UtW, const gsl_vector *Uty,
+	      const gsl_vector *Utx, gsl_matrix *Uab) {
 	size_t index_ab;
 	size_t n_cvt=UtW->size2;
 
@@ -993,10 +1062,13 @@ void CalcUab (const gsl_matrix *UtW, const gsl_vector *Uty, const gsl_vector *Ut
 		gsl_vector_view Uab_col=gsl_matrix_column (Uab, index_ab);
 
 		if (b==n_cvt+2) {gsl_vector_memcpy (&Uab_col.vector, Uty);}
-		else if (b==n_cvt+1) {gsl_vector_memcpy (&Uab_col.vector, Utx);}
+		else if (b==n_cvt+1) {
+		  gsl_vector_memcpy (&Uab_col.vector, Utx);
+		}
 		else {
-			gsl_vector_const_view UtW_col=gsl_matrix_const_column (UtW, b-1);
-			gsl_vector_memcpy (&Uab_col.vector, &UtW_col.vector);
+		  gsl_vector_const_view UtW_col=
+		    gsl_matrix_const_column (UtW, b-1);
+		  gsl_vector_memcpy (&Uab_col.vector, &UtW_col.vector);
 		}
 
 		gsl_vector_mul(&Uab_col.vector, Utx);
@@ -1005,10 +1077,7 @@ void CalcUab (const gsl_matrix *UtW, const gsl_vector *Uty, const gsl_vector *Ut
 	return;
 }
 
-
-
-void Calcab (const gsl_matrix *W, const gsl_vector *y, gsl_vector *ab)
-{
+void Calcab (const gsl_matrix *W, const gsl_vector *y, gsl_vector *ab) {
 	size_t index_ab;
 	size_t n_cvt=W->size2;
 
@@ -1019,10 +1088,12 @@ void Calcab (const gsl_matrix *W, const gsl_vector *y, gsl_vector *ab)
 	for (size_t a=1; a<=n_cvt+2; ++a) {
 		if (a==n_cvt+1) {continue;}
 
-		if (a==n_cvt+2) {gsl_vector_memcpy (v_a, y);}
+		if (a==n_cvt+2) {
+		  gsl_vector_memcpy (v_a, y);
+		}
 		else {
-			gsl_vector_const_view W_col=gsl_matrix_const_column (W, a-1);
-			gsl_vector_memcpy (v_a, &W_col.vector);
+		  gsl_vector_const_view W_col=gsl_matrix_const_column (W, a-1);
+		  gsl_vector_memcpy (v_a, &W_col.vector);
 		}
 
 		for (size_t b=a; b>=1; --b) {
@@ -1030,10 +1101,13 @@ void Calcab (const gsl_matrix *W, const gsl_vector *y, gsl_vector *ab)
 
 			index_ab=GetabIndex (a, b, n_cvt);
 
-			if (b==n_cvt+2) {gsl_vector_memcpy (v_b, y);}
+			if (b==n_cvt+2) {
+			  gsl_vector_memcpy (v_b, y);
+			}
 			else {
-				gsl_vector_const_view W_col=gsl_matrix_const_column (W, b-1);
-				gsl_vector_memcpy (v_b, &W_col.vector);
+			  gsl_vector_const_view W_col=
+			    gsl_matrix_const_column (W, b-1);
+			  gsl_vector_memcpy (v_b, &W_col.vector);
 			}
 
 			gsl_blas_ddot (v_a, v_b, &d);
@@ -1046,9 +1120,8 @@ void Calcab (const gsl_matrix *W, const gsl_vector *y, gsl_vector *ab)
 	return;
 }
 
-
-void Calcab (const gsl_matrix *W, const gsl_vector *y, const gsl_vector *x, gsl_vector *ab)
-{
+void Calcab (const gsl_matrix *W, const gsl_vector *y, const gsl_vector *x,
+	     gsl_vector *ab) {
 	size_t index_ab;
 	size_t n_cvt=W->size2;
 
@@ -1061,8 +1134,8 @@ void Calcab (const gsl_matrix *W, const gsl_vector *y, const gsl_vector *x, gsl_
 		if (b==n_cvt+2) {gsl_vector_memcpy (v_b, y);}
 		else if (b==n_cvt+1) {gsl_vector_memcpy (v_b, x);}
 		else {
-			gsl_vector_const_view W_col=gsl_matrix_const_column (W, b-1);
-			gsl_vector_memcpy (v_b, &W_col.vector);
+		  gsl_vector_const_view W_col=gsl_matrix_const_column (W, b-1);
+		  gsl_vector_memcpy (v_b, &W_col.vector);
 		}
 
 		gsl_blas_ddot (x, v_b, &d);
@@ -1070,31 +1143,31 @@ void Calcab (const gsl_matrix *W, const gsl_vector *y, const gsl_vector *x, gsl_
 	}
 
 	gsl_vector_free (v_b);
-
 	return;
 }
 
-
-
-
-
-void LMM::AnalyzeGene (const gsl_matrix *U, const gsl_vector *eval, const gsl_matrix *UtW, const gsl_vector *Utx, const gsl_matrix *W, const gsl_vector *x)
-{
+void LMM::AnalyzeGene (const gsl_matrix *U, const gsl_vector *eval,
+		       const gsl_matrix *UtW, const gsl_vector *Utx,
+		       const gsl_matrix *W, const gsl_vector *x) {
 	igzstream infile (file_gene.c_str(), igzstream::in);
-	if (!infile) {cout<<"error reading gene expression file:"<<file_gene<<endl; return;}
+	if (!infile) {
+	  cout<<"error reading gene expression file:"<<file_gene<<endl;
+	  return;
+	}
 
 	clock_t time_start=clock();
 
 	string line;
 	char *ch_ptr;
 
-	double lambda_mle=0, lambda_remle=0, beta=0, se=0, p_wald=0, p_lrt=0, p_score=0;
+	double lambda_mle=0, lambda_remle=0, beta=0, se=0, p_wald=0;
+	double p_lrt=0, p_score=0;
 	double logl_H1=0.0, logl_H0=0.0, l_H0;
 	int c_phen;
-	string rs; //gene id
+	string rs; // Gene id.
 	double d;
 
-	//Calculate basic quantities
+	// Calculate basic quantities.
 	size_t n_index=(n_cvt+2+1)*(n_cvt+2)/2;
 
 	gsl_vector *y=gsl_vector_alloc (U->size1);
@@ -1102,12 +1175,14 @@ void LMM::AnalyzeGene (const gsl_matrix *U, const gsl_vector *eval, const gsl_ma
 	gsl_matrix *Uab=gsl_matrix_alloc (U->size2, n_index);
 	gsl_vector *ab=gsl_vector_alloc (n_index);
 
-	//header
+	// Header.
 	getline(infile, line);
 
 	for (size_t t=0; t<ng_total; t++) {
 		!safeGetline(infile, line).eof();
-		if (t%d_pace==0 || t==ng_total-1) {ProgressBar ("Performing Analysis ", t, ng_total-1);}
+		if (t%d_pace==0 || t==ng_total-1) {
+		  ProgressBar ("Performing Analysis ", t, ng_total-1);
+		}
 		ch_ptr=strtok ((char *)line.c_str(), " , \t");
 		rs=ch_ptr;
 
@@ -1126,7 +1201,7 @@ void LMM::AnalyzeGene (const gsl_matrix *U, const gsl_vector *eval, const gsl_ma
 		gsl_blas_dgemv (CblasTrans, 1.0, U, y, 0.0, Uty);
 		time_UtX+=(clock()-time_start)/(double(CLOCKS_PER_SEC)*60.0);
 
-		//calculate null
+		// Calculate null.
 		time_start=clock();
 
 		gsl_matrix_set_zero (Uab);
@@ -1135,32 +1210,36 @@ void LMM::AnalyzeGene (const gsl_matrix *U, const gsl_vector *eval, const gsl_ma
 		FUNC_PARAM param0={false, ni_test, n_cvt, eval, Uab, ab, 0};
 
 		if (a_mode==2 || a_mode==3 || a_mode==4) {
-			CalcLambda('L', param0, l_min, l_max, n_region, l_H0, logl_H0);
+		  CalcLambda('L', param0, l_min, l_max, n_region,
+			     l_H0, logl_H0);
 		}
 
-		//calculate alternative
+		// Calculate alternative.
 		CalcUab(UtW, Uty, Utx, Uab);
 		FUNC_PARAM param1={false, ni_test, n_cvt, eval, Uab, ab, 0};
 
-		//3 is before 1
+		//3 is before 1.
 		if (a_mode==3 || a_mode==4) {
 			CalcRLScore (l_H0, param1, beta, se, p_score);
 		}
 
 		if (a_mode==1 || a_mode==4) {
-			CalcLambda ('R', param1, l_min, l_max, n_region, lambda_remle, logl_H1);
-			CalcRLWald (lambda_remle, param1, beta, se, p_wald);
+		  CalcLambda ('R', param1, l_min, l_max, n_region,
+			      lambda_remle, logl_H1);
+		  CalcRLWald (lambda_remle, param1, beta, se, p_wald);
 		}
 
 		if (a_mode==2 || a_mode==4) {
-			CalcLambda ('L', param1, l_min, l_max, n_region, lambda_mle, logl_H1);
-			p_lrt=gsl_cdf_chisq_Q (2.0*(logl_H1-logl_H0), 1);
+		  CalcLambda ('L', param1, l_min, l_max, n_region,
+			      lambda_mle, logl_H1);
+		  p_lrt=gsl_cdf_chisq_Q (2.0*(logl_H1-logl_H0), 1);
 		}
 
 		time_opt+=(clock()-time_start)/(double(CLOCKS_PER_SEC)*60.0);
 
-		//store summary data
-		SUMSTAT SNPs={beta, se, lambda_remle, lambda_mle, p_wald, p_lrt, p_score};
+		// Store summary data.
+		SUMSTAT SNPs={beta, se, lambda_remle, lambda_mle,
+			      p_wald, p_lrt, p_score};
 		sumStat.push_back(SNPs);
     }
 	cout<<endl;
@@ -1176,27 +1255,27 @@ void LMM::AnalyzeGene (const gsl_matrix *U, const gsl_vector *eval, const gsl_ma
 	return;
 }
 
-
-
-
-
-void LMM::AnalyzeBimbam (const gsl_matrix *U, const gsl_vector *eval, const gsl_matrix *UtW, const gsl_vector *Uty, const gsl_matrix *W, const gsl_vector *y)
-{
+void LMM::AnalyzeBimbam (const gsl_matrix *U, const gsl_vector *eval,
+			 const gsl_matrix *UtW, const gsl_vector *Uty,
+			 const gsl_matrix *W, const gsl_vector *y) {
 	igzstream infile (file_geno.c_str(), igzstream::in);
-//	ifstream infile (file_geno.c_str(), ifstream::in);
-	if (!infile) {cout<<"error reading genotype file:"<<file_geno<<endl; return;}
+	if (!infile) {
+	  cout<<"error reading genotype file:"<<file_geno<<endl;
+	  return;
+	}
 
 	clock_t time_start=clock();
 
 	string line;
 	char *ch_ptr;
 
-	double lambda_mle=0, lambda_remle=0, beta=0, se=0, p_wald=0, p_lrt=0, p_score=0;
+	double lambda_mle=0, lambda_remle=0, beta=0, se=0, p_wald=0;
+	double p_lrt=0, p_score=0;
 	double logl_H1=0.0;
 	int n_miss, c_phen;
 	double geno, x_mean;
 
-	//Calculate basic quantities
+	// Calculate basic quantities.
 	size_t n_index=(n_cvt+2+1)*(n_cvt+2)/2;
 
 	gsl_vector *x=gsl_vector_alloc (U->size1);
@@ -1205,7 +1284,7 @@ void LMM::AnalyzeBimbam (const gsl_matrix *U, const gsl_vector *eval, const gsl_
 	gsl_matrix *Uab=gsl_matrix_alloc (U->size2, n_index);
 	gsl_vector *ab=gsl_vector_alloc (n_index);
 
-	//create a large matrix
+	// Create a large matrix.
 	size_t msize=10000;
 	gsl_matrix *Xlarge=gsl_matrix_alloc (U->size1, msize);
 	gsl_matrix *UtXlarge=gsl_matrix_alloc (U->size1, msize);
@@ -1213,10 +1292,6 @@ void LMM::AnalyzeBimbam (const gsl_matrix *U, const gsl_vector *eval, const gsl_
 
 	gsl_matrix_set_zero (Uab);
 	CalcUab (UtW, Uty, Uab);
-//	if (e_mode!=0) {
-//		gsl_vector_set_zero (ab);
-//		Calcab (W, y, ab);
-//	}
 
 	//start reading genotypes and analyze
 	size_t c=0, t_last=0;
@@ -1225,9 +1300,10 @@ void LMM::AnalyzeBimbam (const gsl_matrix *U, const gsl_vector *eval, const gsl_
 	  t_last++;
 	}
 	for (size_t t=0; t<indicator_snp.size(); ++t) {
-//		if (t>1) {break;}
 		!safeGetline(infile, line).eof();
-		if (t%d_pace==0 || t==(ns_total-1)) {ProgressBar ("Reading SNPs  ", t, ns_total-1);}
+		if (t%d_pace==0 || t==(ns_total-1)) {
+		  ProgressBar ("Reading SNPs  ", t, ns_total-1);
+		}
 		if (indicator_snp[t]==0) {continue;}
 
 		ch_ptr=strtok ((char *)line.c_str(), " , \t");
@@ -1240,7 +1316,9 @@ void LMM::AnalyzeBimbam (const gsl_matrix *U, const gsl_vector *eval, const gsl_
 			ch_ptr=strtok (NULL, " , \t");
 			if (indicator_idv[i]==0) {continue;}
 
-			if (strcmp(ch_ptr, "NA")==0) {gsl_vector_set(x_miss, c_phen, 0.0); n_miss++;}
+			if (strcmp(ch_ptr, "NA")==0) {
+			  gsl_vector_set(x_miss, c_phen, 0.0); n_miss++;
+			}
 			else {
 				geno=atof(ch_ptr);
 
@@ -1254,20 +1332,11 @@ void LMM::AnalyzeBimbam (const gsl_matrix *U, const gsl_vector *eval, const gsl_
 		x_mean/=(double)(ni_test-n_miss);
 
 		for (size_t i=0; i<ni_test; ++i) {
-			if (gsl_vector_get (x_miss, i)==0) {gsl_vector_set(x, i, x_mean);}
-			//geno=gsl_vector_get(x, i);
-			//if (x_mean>1) {
-			//	gsl_vector_set(x, i, 2-geno);
-			//}
+			if (gsl_vector_get (x_miss, i)==0) {
+			  gsl_vector_set(x, i, x_mean);
+			}
 		}
 
-		/*
-		//calculate statistics
-		time_start=clock();
-		gsl_blas_dgemv (CblasTrans, 1.0, U, x, 0.0, Utx);
-		time_UtX+=(clock()-time_start)/(double(CLOCKS_PER_SEC)*60.0);
-		*/
-
 		gsl_vector_view Xlarge_col=gsl_matrix_column (Xlarge, c%msize);
 		gsl_vector_memcpy (&Xlarge_col.vector, x);
 		c++;
@@ -1276,49 +1345,52 @@ void LMM::AnalyzeBimbam (const gsl_matrix *U, const gsl_vector *eval, const gsl_
 		  size_t l=0;
 		  if (c%msize==0) {l=msize;} else {l=c%msize;}
 
-		  gsl_matrix_view Xlarge_sub=gsl_matrix_submatrix(Xlarge, 0, 0, Xlarge->size1, l);
-		  gsl_matrix_view UtXlarge_sub=gsl_matrix_submatrix(UtXlarge, 0, 0, UtXlarge->size1, l);
+		  gsl_matrix_view Xlarge_sub=
+		    gsl_matrix_submatrix(Xlarge, 0, 0, Xlarge->size1, l);
+		  gsl_matrix_view UtXlarge_sub=
+		    gsl_matrix_submatrix(UtXlarge, 0, 0, UtXlarge->size1, l);
 
 		  time_start=clock();
-		  eigenlib_dgemm ("T", "N", 1.0, U, &Xlarge_sub.matrix, 0.0, &UtXlarge_sub.matrix);
+		  eigenlib_dgemm ("T", "N", 1.0, U, &Xlarge_sub.matrix,
+				  0.0, &UtXlarge_sub.matrix);
 		  time_UtX+=(clock()-time_start)/(double(CLOCKS_PER_SEC)*60.0);
 
 		  gsl_matrix_set_zero (Xlarge);
 
 		  for (size_t i=0; i<l; i++) {
-		    gsl_vector_view UtXlarge_col=gsl_matrix_column (UtXlarge, i);
+		    gsl_vector_view UtXlarge_col=
+		      gsl_matrix_column (UtXlarge, i);
 		    gsl_vector_memcpy (Utx, &UtXlarge_col.vector);
 
 		    CalcUab(UtW, Uty, Utx, Uab);
-		    //		if (e_mode!=0) {
-		    //			Calcab (W, y, x, ab);
-		    //		}
 
 		    time_start=clock();
-		    FUNC_PARAM param1={false, ni_test, n_cvt, eval, Uab, ab, 0};
+		    FUNC_PARAM param1=
+		      {false, ni_test, n_cvt, eval, Uab, ab, 0};
 
-		    //3 is before 1
+		    // 3 is before 1.
 		    if (a_mode==3 || a_mode==4) {
 		      CalcRLScore (l_mle_null, param1, beta, se, p_score);
 		    }
 
 		    if (a_mode==1 || a_mode==4) {
-		      CalcLambda ('R', param1, l_min, l_max, n_region, lambda_remle, logl_H1);
+		      CalcLambda ('R', param1, l_min, l_max, n_region,
+				  lambda_remle, logl_H1);
 		      CalcRLWald (lambda_remle, param1, beta, se, p_wald);
 		    }
 
 		    if (a_mode==2 || a_mode==4) {
-		      CalcLambda ('L', param1, l_min, l_max, n_region, lambda_mle, logl_H1);
+		      CalcLambda ('L', param1, l_min, l_max, n_region,
+				  lambda_mle, logl_H1);
 		      p_lrt=gsl_cdf_chisq_Q (2.0*(logl_H1-logl_mle_H0), 1);
 		    }
 
-		    //if (p_wald<0) {cout<<t<<"\t"<<i<<"\t"<<l<<endl;}
-		    //if (x_mean>1) {beta*=-1;}
+		    time_opt+=(clock()-time_start)/
+		      (double(CLOCKS_PER_SEC)*60.0);
 
-		    time_opt+=(clock()-time_start)/(double(CLOCKS_PER_SEC)*60.0);
-
-		    //store summary data
-		    SUMSTAT SNPs={beta, se, lambda_remle, lambda_mle, p_wald, p_lrt, p_score};
+		    // Store summary data.
+		    SUMSTAT SNPs={beta, se, lambda_remle, lambda_mle,
+				  p_wald, p_lrt, p_score};
 
 		    sumStat.push_back(SNPs);
 		  }
@@ -1341,14 +1413,9 @@ void LMM::AnalyzeBimbam (const gsl_matrix *U, const gsl_vector *eval, const gsl_
 	return;
 }
 
-
-
-
-
-
-
-void LMM::AnalyzePlink (const gsl_matrix *U, const gsl_vector *eval, const gsl_matrix *UtW, const gsl_vector *Uty, const gsl_matrix *W, const gsl_vector *y)
-{
+void LMM::AnalyzePlink (const gsl_matrix *U, const gsl_vector *eval,
+			const gsl_matrix *UtW, const gsl_vector *Uty,
+			const gsl_matrix *W, const gsl_vector *y) {
 	string file_bed=file_bfile+".bed";
 	ifstream infile (file_bed.c_str(), ios::binary);
 	if (!infile) {cout<<"error reading bed file:"<<file_bed<<endl; return;}
@@ -1358,12 +1425,13 @@ void LMM::AnalyzePlink (const gsl_matrix *U, const gsl_vector *eval, const gsl_m
 	char ch[1];
 	bitset<8> b;
 
-	double lambda_mle=0, lambda_remle=0, beta=0, se=0, p_wald=0, p_lrt=0, p_score=0;
+	double lambda_mle=0, lambda_remle=0, beta=0, se=0, p_wald=0;
+	double p_lrt=0, p_score=0;
 	double logl_H1=0.0;
 	int n_bit, n_miss, ci_total, ci_test;
 	double geno, x_mean;
 
-	//Calculate basic quantities
+	// Calculate basic quantities.
 	size_t n_index=(n_cvt+2+1)*(n_cvt+2)/2;
 
 	gsl_vector *x=gsl_vector_alloc (U->size1);
@@ -1371,7 +1439,7 @@ void LMM::AnalyzePlink (const gsl_matrix *U, const gsl_vector *eval, const gsl_m
 	gsl_matrix *Uab=gsl_matrix_alloc (U->size2, n_index);
 	gsl_vector *ab=gsl_vector_alloc (n_index);
 
-	//create a large matrix
+	// Create a large matrix.
 	size_t msize=10000;
 	gsl_matrix *Xlarge=gsl_matrix_alloc (U->size1, msize);
 	gsl_matrix *UtXlarge=gsl_matrix_alloc (U->size1, msize);
@@ -1379,16 +1447,12 @@ void LMM::AnalyzePlink (const gsl_matrix *U, const gsl_vector *eval, const gsl_m
 
 	gsl_matrix_set_zero (Uab);
 	CalcUab (UtW, Uty, Uab);
-//	if (e_mode!=0) {
-//		gsl_vector_set_zero (ab);
-//		Calcab (W, y, ab);
-//	}
 
-	//calculate n_bit and c, the number of bit for each snp
+	// Calculate n_bit and c, the number of bit for each SNP.
 	if (ni_total%4==0) {n_bit=ni_total/4;}
 	else {n_bit=ni_total/4+1; }
 
-	//print the first three majic numbers
+	// Print the first three magic numbers.
 	for (int i=0; i<3; ++i) {
 		infile.read(ch,1);
 		b=ch[0];
@@ -1400,31 +1464,44 @@ void LMM::AnalyzePlink (const gsl_matrix *U, const gsl_vector *eval, const gsl_m
 	  t_last++;
 	}
 	for (vector<SNPINFO>::size_type t=0; t<snpInfo.size(); ++t) {
-		if (t%d_pace==0 || t==snpInfo.size()-1) {ProgressBar ("Reading SNPs  ", t, snpInfo.size()-1);}
+		if (t%d_pace==0 || t==snpInfo.size()-1) {
+		  ProgressBar ("Reading SNPs  ", t, snpInfo.size()-1);
+		}
 		if (indicator_snp[t]==0) {continue;}
 
-		infile.seekg(t*n_bit+3);		//n_bit, and 3 is the number of magic numbers
+		// n_bit, and 3 is the number of magic numbers.
+		infile.seekg(t*n_bit+3);
 
-		//read genotypes
+		// Read genotypes.
 		x_mean=0.0;	n_miss=0; ci_total=0; ci_test=0;
 		for (int i=0; i<n_bit; ++i) {
 			infile.read(ch,1);
 			b=ch[0];
-			for (size_t j=0; j<4; ++j) {                //minor allele homozygous: 2.0; major: 0.0;
-				if ((i==(n_bit-1)) && ci_total==(int)ni_total) {break;}
-				if (indicator_idv[ci_total]==0) {ci_total++; continue;}
 
-				if (b[2*j]==0) {
-					if (b[2*j+1]==0) {gsl_vector_set(x, ci_test, 2); x_mean+=2.0; }
-					else {gsl_vector_set(x, ci_test, 1); x_mean+=1.0; }
-				}
-				else {
-					if (b[2*j+1]==1) {gsl_vector_set(x, ci_test, 0); }
-					else {gsl_vector_set(x, ci_test, -9); n_miss++; }
-				}
-
-				ci_total++;
-				ci_test++;
+			// Minor allele homozygous: 2.0; major: 0.0.
+			for (size_t j=0; j<4; ++j) {                
+			  if ((i==(n_bit-1)) && ci_total==(int)ni_total) {
+			    break;
+			  }
+			  if (indicator_idv[ci_total]==0) {
+			    ci_total++;
+			    continue;
+			  }
+			  
+			  if (b[2*j]==0) {
+			    if (b[2*j+1]==0) {
+			      gsl_vector_set(x, ci_test, 2);
+			      x_mean+=2.0;
+			    }
+			    else {gsl_vector_set(x, ci_test, 1); x_mean+=1.0; }
+			  }
+			  else {
+			    if (b[2*j+1]==1) {gsl_vector_set(x, ci_test, 0); }
+			    else {gsl_vector_set(x, ci_test, -9); n_miss++; }
+			  }
+			  
+			  ci_total++;
+			  ci_test++;
 			}
 		}
 
@@ -1432,19 +1509,12 @@ void LMM::AnalyzePlink (const gsl_matrix *U, const gsl_vector *eval, const gsl_m
 
 		for (size_t i=0; i<ni_test; ++i) {
 			geno=gsl_vector_get(x,i);
-			if (geno==-9) {gsl_vector_set(x, i, x_mean); geno=x_mean;}
-			//if (x_mean>1) {
-			//gsl_vector_set(x, i, 2-geno);
-			//}
+			if (geno==-9) {
+			  gsl_vector_set(x, i, x_mean);
+			  geno=x_mean;
+			}
 		}
 
-		/*
-		//calculate statistics
-		time_start=clock();
-		gsl_blas_dgemv (CblasTrans, 1.0, U, x, 0.0, Utx);
-		time_UtX+=(clock()-time_start)/(double(CLOCKS_PER_SEC)*60.0);
-		*/
-
 		gsl_vector_view Xlarge_col=gsl_matrix_column (Xlarge, c%msize);
 		gsl_vector_memcpy (&Xlarge_col.vector, x);
 		c++;
@@ -1453,52 +1523,56 @@ void LMM::AnalyzePlink (const gsl_matrix *U, const gsl_vector *eval, const gsl_m
 		  size_t l=0;
 		  if (c%msize==0) {l=msize;} else {l=c%msize;}
 
-		  gsl_matrix_view Xlarge_sub=gsl_matrix_submatrix(Xlarge, 0, 0, Xlarge->size1, l);
-		  gsl_matrix_view UtXlarge_sub=gsl_matrix_submatrix(UtXlarge, 0, 0, UtXlarge->size1, l);
+		  gsl_matrix_view Xlarge_sub=
+		    gsl_matrix_submatrix(Xlarge, 0, 0, Xlarge->size1, l);
+		  gsl_matrix_view UtXlarge_sub=
+		    gsl_matrix_submatrix(UtXlarge, 0, 0, UtXlarge->size1, l);
 
 		  time_start=clock();
-		  eigenlib_dgemm ("T", "N", 1.0, U, &Xlarge_sub.matrix, 0.0, &UtXlarge_sub.matrix);
+		  eigenlib_dgemm ("T", "N", 1.0, U, &Xlarge_sub.matrix,
+				  0.0, &UtXlarge_sub.matrix);
 		  time_UtX+=(clock()-time_start)/(double(CLOCKS_PER_SEC)*60.0);
 
 		  gsl_matrix_set_zero (Xlarge);
 
 		  for (size_t i=0; i<l; i++) {
-		    gsl_vector_view UtXlarge_col=gsl_matrix_column (UtXlarge, i);
+		    gsl_vector_view UtXlarge_col=
+		      gsl_matrix_column (UtXlarge, i);
 		    gsl_vector_memcpy (Utx, &UtXlarge_col.vector);
 
 		    CalcUab(UtW, Uty, Utx, Uab);
-		    //		if (e_mode!=0) {
-		    //			Calcab (W, y, x, ab);
-		    //		}
 
 		    time_start=clock();
-		    FUNC_PARAM param1={false, ni_test, n_cvt, eval, Uab, ab, 0};
+		    FUNC_PARAM param1={false, ni_test, n_cvt, eval,
+				       Uab, ab, 0};
 
-		    //3 is before 1, for beta
+		    // 3 is before 1, for beta.
 		    if (a_mode==3 || a_mode==4) {
 		      CalcRLScore (l_mle_null, param1, beta, se, p_score);
 		    }
 
 		    if (a_mode==1 || a_mode==4) {
-		      CalcLambda ('R', param1, l_min, l_max, n_region, lambda_remle, logl_H1);
+		      CalcLambda ('R', param1, l_min, l_max, n_region,
+				  lambda_remle, logl_H1);
 		      CalcRLWald (lambda_remle, param1, beta, se, p_wald);
 		    }
 
 		    if (a_mode==2 || a_mode==4) {
-		      CalcLambda ('L', param1, l_min, l_max, n_region, lambda_mle, logl_H1);
+		      CalcLambda ('L', param1, l_min, l_max, n_region,
+				  lambda_mle, logl_H1);
 		      p_lrt=gsl_cdf_chisq_Q (2.0*(logl_H1-logl_mle_H0), 1);
 		    }
 
-		    //if (x_mean>1) {beta*=-1;}
+		    time_opt+=(clock()-time_start)/
+		      (double(CLOCKS_PER_SEC)*60.0);
 
-		    time_opt+=(clock()-time_start)/(double(CLOCKS_PER_SEC)*60.0);
-
-		    //store summary data
-		    SUMSTAT SNPs={beta, se, lambda_remle, lambda_mle, p_wald, p_lrt, p_score};
+		    // Store summary data.
+		    SUMSTAT SNPs={beta, se, lambda_remle, lambda_mle,
+				  p_wald, p_lrt, p_score};
 		    sumStat.push_back(SNPs);
 		  }
 		}
-    }
+	}
 	cout<<endl;
 
 	gsl_vector_free (x);
@@ -1515,25 +1589,25 @@ void LMM::AnalyzePlink (const gsl_matrix *U, const gsl_vector *eval, const gsl_m
 	return;
 }
 
-
-// WJA added
-#include <assert.h>
-void LMM::Analyzebgen (const gsl_matrix *U, const gsl_vector *eval, const gsl_matrix *UtW, const gsl_vector *Uty, const gsl_matrix *W, const gsl_vector *y)
-{
+// WJA added.
+void LMM::Analyzebgen (const gsl_matrix *U, const gsl_vector *eval,
+		       const gsl_matrix *UtW, const gsl_vector *Uty,
+		       const gsl_matrix *W, const gsl_vector *y) {
 	string file_bgen=file_oxford+".bgen";
 	ifstream infile (file_bgen.c_str(), ios::binary);
-	if (!infile) {cout<<"error reading bgen file:"<<file_bgen<<endl; return;}
-
+	if (!infile) {
+	  cout<<"error reading bgen file:"<<file_bgen<<endl;
+	  return;
+	}
 
 	clock_t time_start=clock();
-
-
-	double lambda_mle=0, lambda_remle=0, beta=0, se=0, p_wald=0, p_lrt=0, p_score=0;
+	double lambda_mle=0, lambda_remle=0, beta=0, se=0, p_wald=0;
+	double p_lrt=0, p_score=0;
 	double logl_H1=0.0;
 	int n_miss, c_phen;
 	double geno, x_mean;
 
-	//Calculate basic quantities
+	// Calculate basic quantities.
 	size_t n_index=(n_cvt+2+1)*(n_cvt+2)/2;
 
 	gsl_vector *x=gsl_vector_alloc (U->size1);
@@ -1542,7 +1616,7 @@ void LMM::Analyzebgen (const gsl_matrix *U, const gsl_vector *eval, const gsl_ma
 	gsl_matrix *Uab=gsl_matrix_alloc (U->size2, n_index);
 	gsl_vector *ab=gsl_vector_alloc (n_index);
 
-	//create a large matrix
+	// Create a large matrix.
 	size_t msize=10000;
 	gsl_matrix *Xlarge=gsl_matrix_alloc (U->size1, msize);
 	gsl_matrix *UtXlarge=gsl_matrix_alloc (U->size1, msize);
@@ -1550,12 +1624,8 @@ void LMM::Analyzebgen (const gsl_matrix *U, const gsl_vector *eval, const gsl_ma
 
 	gsl_matrix_set_zero (Uab);
 	CalcUab (UtW, Uty, Uab);
-//	if (e_mode!=0) {
-//		gsl_vector_set_zero (ab);
-//		Calcab (W, y, ab);
-//	}
 
-	// read in header
+	// Read in header.
 	uint32_t bgen_snp_block_offset;
 	uint32_t bgen_header_length;
 	uint32_t bgen_nsamples;
@@ -1573,11 +1643,11 @@ void LMM::Analyzebgen (const gsl_matrix *U, const gsl_vector *eval, const gsl_ma
 	infile.read(reinterpret_cast<char*>(&bgen_flags),4);
 	bgen_snp_block_offset-=4;
 	bool CompressedSNPBlocks=bgen_flags&0x1;
-//	bool LongIds=bgen_flags&0x4;
 
 	infile.ignore(bgen_snp_block_offset);
 
-	double bgen_geno_prob_AA, bgen_geno_prob_AB, bgen_geno_prob_BB, bgen_geno_prob_non_miss;
+	double bgen_geno_prob_AA, bgen_geno_prob_AB, bgen_geno_prob_BB;
+	double bgen_geno_prob_non_miss;
 
 	uint32_t bgen_N;
 	uint16_t bgen_LS;
@@ -1593,11 +1663,10 @@ void LMM::Analyzebgen (const gsl_matrix *U, const gsl_vector *eval, const gsl_ma
 	string id;
 	string rs;
 	string chr;
-	std::cout<<"Warning: WJA hard coded SNP missingness threshold of 10%"<<std::endl;
+	std::cout << "Warning: WJA hard coded SNP missingness " <<
+	  "threshold of 10%"<<std::endl;
 
-
-
-	//start reading genotypes and analyze
+	// Start reading genotypes and analyze.
 	size_t c=0, t_last=0;
 	for (size_t t=0; t<indicator_snp.size(); ++t) {
 	  if (indicator_snp[t]==0) {continue;}
@@ -1605,11 +1674,12 @@ void LMM::Analyzebgen (const gsl_matrix *U, const gsl_vector *eval, const gsl_ma
 	}
 	for (size_t t=0; t<indicator_snp.size(); ++t)
 	{
-
-//		if (t>1) {break;}
-		if (t%d_pace==0 || t==(ns_total-1)) {ProgressBar ("Reading SNPs  ", t, ns_total-1);}
+		if (t%d_pace==0 || t==(ns_total-1)) {
+		  ProgressBar ("Reading SNPs  ", t, ns_total-1);
+		}
 		if (indicator_snp[t]==0) {continue;}
-		// read SNP header
+		
+		// Read SNP header.
 		id.clear();
 		rs.clear();
 		chr.clear();
@@ -1641,89 +1711,86 @@ void LMM::Analyzebgen (const gsl_matrix *U, const gsl_vector *eval, const gsl_ma
 		bgen_B_allele.resize(bgen_LB);
 		infile.read(&bgen_B_allele[0], bgen_LB);
 
-
-
-
 		uint16_t unzipped_data[3*bgen_N];
 
 		if (indicator_snp[t]==0) {
 			if(CompressedSNPBlocks)
-				infile.read(reinterpret_cast<char*>(&bgen_P),4);
+			  infile.read(reinterpret_cast<char*>(&bgen_P),4);
 			else
-				bgen_P=6*bgen_N;
+			  bgen_P=6*bgen_N;
 
 			infile.ignore(static_cast<size_t>(bgen_P));
 
 			continue;
 		}
 
-
-		if(CompressedSNPBlocks)
-		{
-
-
+		if(CompressedSNPBlocks) {
 			infile.read(reinterpret_cast<char*>(&bgen_P),4);
 			uint8_t zipped_data[bgen_P];
 
 			unzipped_data_size=6*bgen_N;
 
-			infile.read(reinterpret_cast<char*>(zipped_data),bgen_P);
+			infile.read(reinterpret_cast<char*>(zipped_data),
+				    bgen_P);
 
-			int result=uncompress(reinterpret_cast<Bytef*>(unzipped_data), reinterpret_cast<uLongf*>(&unzipped_data_size), reinterpret_cast<Bytef*>(zipped_data), static_cast<uLong> (bgen_P));
+			int result=
+			  uncompress(reinterpret_cast<Bytef*>(unzipped_data),
+			    reinterpret_cast<uLongf*>(&unzipped_data_size),
+			    reinterpret_cast<Bytef*>(zipped_data),
+			    static_cast<uLong> (bgen_P));
 			assert(result == Z_OK);
 
 		}
 		else
 		{
 
-			bgen_P=6*bgen_N;
-			infile.read(reinterpret_cast<char*>(unzipped_data),bgen_P);
+		  bgen_P=6*bgen_N;
+		  infile.read(reinterpret_cast<char*>(unzipped_data),bgen_P);
 		}
 
 		x_mean=0.0; c_phen=0; n_miss=0;
 		gsl_vector_set_zero(x_miss);
 		for (size_t i=0; i<bgen_N; ++i) {
-			if (indicator_idv[i]==0) {continue;}
-
-
-				bgen_geno_prob_AA=static_cast<double>(unzipped_data[i*3])/32768.0;
-				bgen_geno_prob_AB=static_cast<double>(unzipped_data[i*3+1])/32768.0;
-				bgen_geno_prob_BB=static_cast<double>(unzipped_data[i*3+2])/32768.0;
-				// WJA
-				bgen_geno_prob_non_miss=bgen_geno_prob_AA+bgen_geno_prob_AB+bgen_geno_prob_BB;
-				if (bgen_geno_prob_non_miss<0.9) {gsl_vector_set(x_miss, c_phen, 0.0); n_miss++;}
-				else {
-
-					bgen_geno_prob_AA/=bgen_geno_prob_non_miss;
-					bgen_geno_prob_AB/=bgen_geno_prob_non_miss;
-					bgen_geno_prob_BB/=bgen_geno_prob_non_miss;
-
-					geno=2.0*bgen_geno_prob_BB+bgen_geno_prob_AB;
-
-					gsl_vector_set(x, c_phen, geno);
-					gsl_vector_set(x_miss, c_phen, 1.0);
-					x_mean+=geno;
-			}
-			c_phen++;
+		  if (indicator_idv[i]==0) {continue;}
+		  
+		  bgen_geno_prob_AA=
+		    static_cast<double>(unzipped_data[i*3])/32768.0;
+		  bgen_geno_prob_AB=
+		    static_cast<double>(unzipped_data[i*3+1])/32768.0;
+		  bgen_geno_prob_BB=
+		    static_cast<double>(unzipped_data[i*3+2])/32768.0;
+		  
+		  // WJA.
+		  bgen_geno_prob_non_miss = bgen_geno_prob_AA +
+		    bgen_geno_prob_AB+bgen_geno_prob_BB;
+		  if (bgen_geno_prob_non_miss<0.9) {
+		    gsl_vector_set(x_miss, c_phen, 0.0);
+		    n_miss++;
+		  }
+		  else {
+		    
+		    bgen_geno_prob_AA/=bgen_geno_prob_non_miss;
+		    bgen_geno_prob_AB/=bgen_geno_prob_non_miss;
+		    bgen_geno_prob_BB/=bgen_geno_prob_non_miss;
+		    
+		    geno=2.0*bgen_geno_prob_BB+bgen_geno_prob_AB;
+		    
+		    gsl_vector_set(x, c_phen, geno);
+		    gsl_vector_set(x_miss, c_phen, 1.0);
+		    x_mean+=geno;
+		  }
+		  c_phen++;
 		}
 
 		x_mean/=static_cast<double>(ni_test-n_miss);
 
 		for (size_t i=0; i<ni_test; ++i) {
-			if (gsl_vector_get (x_miss, i)==0) {gsl_vector_set(x, i, x_mean);}
+			if (gsl_vector_get (x_miss, i)==0) {
+			  gsl_vector_set(x, i, x_mean);
+			}
 			geno=gsl_vector_get(x, i);
-			//if (x_mean>1) {
-			//gsl_vector_set(x, i, 2-geno);
-			//}
 		}
 
-		/*
-		//calculate statistics
-		time_start=clock();
-		gsl_blas_dgemv (CblasTrans, 1.0, U, x, 0.0, Utx);
-		time_UtX+=(clock()-time_start)/(double(CLOCKS_PER_SEC)*60.0);
-		*/
-
 		gsl_vector_view Xlarge_col=gsl_matrix_column (Xlarge, c%msize);
 		gsl_vector_memcpy (&Xlarge_col.vector, x);
 		c++;
@@ -1732,48 +1799,51 @@ void LMM::Analyzebgen (const gsl_matrix *U, const gsl_vector *eval, const gsl_ma
 		  size_t l=0;
 		  if (c%msize==0) {l=msize;} else {l=c%msize;}
 
-		  gsl_matrix_view Xlarge_sub=gsl_matrix_submatrix(Xlarge, 0, 0, Xlarge->size1, l);
-		  gsl_matrix_view UtXlarge_sub=gsl_matrix_submatrix(UtXlarge, 0, 0, UtXlarge->size1, l);
+		  gsl_matrix_view Xlarge_sub=
+		    gsl_matrix_submatrix(Xlarge, 0, 0, Xlarge->size1, l);
+		  gsl_matrix_view UtXlarge_sub=
+		    gsl_matrix_submatrix(UtXlarge, 0, 0, UtXlarge->size1, l);
 
 		  time_start=clock();
-		  eigenlib_dgemm ("T", "N", 1.0, U, &Xlarge_sub.matrix, 0.0, &UtXlarge_sub.matrix);
+		  eigenlib_dgemm ("T", "N", 1.0, U, &Xlarge_sub.matrix,
+				  0.0, &UtXlarge_sub.matrix);
 		  time_UtX+=(clock()-time_start)/(double(CLOCKS_PER_SEC)*60.0);
 
 		  gsl_matrix_set_zero (Xlarge);
 
 		  for (size_t i=0; i<l; i++) {
-		    gsl_vector_view UtXlarge_col=gsl_matrix_column (UtXlarge, i);
+		    gsl_vector_view UtXlarge_col=
+		      gsl_matrix_column (UtXlarge, i);
 		    gsl_vector_memcpy (Utx, &UtXlarge_col.vector);
 
 		    CalcUab(UtW, Uty, Utx, Uab);
-		    //		if (e_mode!=0) {
-		    //			Calcab (W, y, x, ab);
-		    //		}
 
 		    time_start=clock();
-		    FUNC_PARAM param1={false, ni_test, n_cvt, eval, Uab, ab, 0};
+		    FUNC_PARAM param1={false,ni_test,n_cvt,eval,Uab,ab,0};
 
-		    //3 is before 1
+		    // 3 is before 1.
 		    if (a_mode==3 || a_mode==4) {
 		      CalcRLScore (l_mle_null, param1, beta, se, p_score);
 		    }
 
 		    if (a_mode==1 || a_mode==4) {
-		      CalcLambda ('R', param1, l_min, l_max, n_region, lambda_remle, logl_H1);
+		      CalcLambda ('R', param1, l_min, l_max, n_region,
+				  lambda_remle, logl_H1);
 		      CalcRLWald (lambda_remle, param1, beta, se, p_wald);
 		    }
 
 		    if (a_mode==2 || a_mode==4) {
-		      CalcLambda ('L', param1, l_min, l_max, n_region, lambda_mle, logl_H1);
+		      CalcLambda ('L', param1, l_min, l_max, n_region,
+				  lambda_mle, logl_H1);
 		      p_lrt=gsl_cdf_chisq_Q (2.0*(logl_H1-logl_mle_H0), 1);
 		    }
 
-		    //if (x_mean>1) {beta*=-1;}
-
-		    time_opt+=(clock()-time_start)/(double(CLOCKS_PER_SEC)*60.0);
+		    time_opt+=(clock()-time_start)/
+		      (double(CLOCKS_PER_SEC)*60.0);
 
-		    //store summary data
-		    SUMSTAT SNPs={beta, se, lambda_remle, lambda_mle, p_wald, p_lrt, p_score};
+		    // Store summary data.
+		    SUMSTAT SNPs={beta, se, lambda_remle, lambda_mle,
+				  p_wald, p_lrt, p_score};
 		    sumStat.push_back(SNPs);
 		  }
 		}
@@ -1793,13 +1863,13 @@ void LMM::Analyzebgen (const gsl_matrix *U, const gsl_vector *eval, const gsl_ma
 	infile.clear();
 
 	return;
-
 }
 
-
-
-void MatrixCalcLR (const gsl_matrix *U, const gsl_matrix *UtX, const gsl_vector *Uty, const gsl_vector *K_eval, const double l_min, const double l_max, const size_t n_region, vector<pair<size_t, double> > &pos_loglr)
-{
+void MatrixCalcLR (const gsl_matrix *U, const gsl_matrix *UtX,
+		   const gsl_vector *Uty, const gsl_vector *K_eval,
+		   const double l_min, const double l_max,
+		   const size_t n_region,
+		   vector<pair<size_t, double> > &pos_loglr) {
 	double logl_H0, logl_H1, log_lr, lambda0, lambda1;
 
 	gsl_vector *w=gsl_vector_alloc (Uty->size);
@@ -1812,7 +1882,7 @@ void MatrixCalcLR (const gsl_matrix *U, const gsl_matrix *UtX, const gsl_vector
 	gsl_vector_view Utw_col=gsl_matrix_column (Utw, 0);
 	gsl_blas_dgemv (CblasTrans, 1.0, U, w, 0.0, &Utw_col.vector);
 
-	CalcUab (Utw, Uty, Uab) ;
+	CalcUab (Utw, Uty, Uab);
 	FUNC_PARAM param0={true, Uty->size, 1, K_eval, Uab, ab, 0};
 
 	CalcLambda('L', param0, l_min, l_max, n_region, lambda0, logl_H0);
@@ -1822,7 +1892,8 @@ void MatrixCalcLR (const gsl_matrix *U, const gsl_matrix *UtX, const gsl_vector
 		CalcUab(Utw, Uty, &UtX_col.vector, Uab);
 		FUNC_PARAM param1={false, UtX->size1, 1, K_eval, Uab, ab, 0};
 
-		CalcLambda ('L', param1, l_min, l_max, n_region, lambda1, logl_H1);
+		CalcLambda ('L', param1, l_min, l_max, n_region, lambda1,
+			    logl_H1);
 		log_lr=logl_H1-logl_H0;
 
 		pos_loglr.push_back(make_pair(i,log_lr) );
@@ -1836,17 +1907,21 @@ void MatrixCalcLR (const gsl_matrix *U, const gsl_matrix *UtX, const gsl_vector
 	return;
 }
 
-
-
-
-void CalcLambda (const char func_name, FUNC_PARAM &params, const double l_min, const double l_max, const size_t n_region, double &lambda, double &logf)
-{
-	if (func_name!='R' && func_name!='L' && func_name!='r' && func_name!='l') {cout<<"func_name only takes 'R' or 'L': 'R' for log-restricted likelihood, 'L' for log-likelihood."<<endl; return;}
+void CalcLambda (const char func_name, FUNC_PARAM &params,
+		 const double l_min, const double l_max,
+		 const size_t n_region, double &lambda, double &logf) {
+	if (func_name!='R' && func_name!='L' && func_name!='r' &&
+	    func_name!='l') {
+	  cout << "func_name only takes 'R' or 'L': 'R' for " <<
+	    "log-restricted likelihood, 'L' for log-likelihood." << endl;
+	  return;
+	}
 
 	vector<pair<double, double> > lambda_lh;
 
-	//evaluate first order derivates in different intervals
-	double lambda_l, lambda_h, lambda_interval=log(l_max/l_min)/(double)n_region;
+	// Evaluate first-order derivates in different intervals.
+	double lambda_l, lambda_h, lambda_interval=
+	  log(l_max/l_min)/(double)n_region;
 	double dev1_l, dev1_h, logf_l, logf_h;
 
 	for (size_t i=0; i<n_region; ++i) {
@@ -1867,7 +1942,7 @@ void CalcLambda (const char func_name, FUNC_PARAM &params, const double l_min, c
 		}
 	}
 
-	//if derivates do not change signs in any interval
+	// If derivates do not change signs in any interval.
 	if (lambda_lh.empty()) {
 		if (func_name=='R' || func_name=='r') {
 			logf_l=LogRL_f (l_min, &params);
@@ -1878,10 +1953,17 @@ void CalcLambda (const char func_name, FUNC_PARAM &params, const double l_min, c
 			logf_h=LogL_f (l_max, &params);
 		}
 
-		if (logf_l>=logf_h) {lambda=l_min; logf=logf_l;} else {lambda=l_max; logf=logf_h;}
+		if (logf_l>=logf_h) {
+		  lambda=l_min;
+		  logf=logf_l;
+		} else {
+		  lambda=l_max;
+		  logf=logf_h;
+		}
 	}
 	else {
-		//if derivates change signs
+	  
+		// If derivates change signs.
 		int status;
 		int iter=0, max_iter=100;
 		double l, l_temp;
@@ -1916,41 +1998,46 @@ void CalcLambda (const char func_name, FUNC_PARAM &params, const double l_min, c
 		s_fdf=gsl_root_fdfsolver_alloc(T_fdf);
 
 		for (vector<double>::size_type i=0; i<lambda_lh.size(); ++i) {
-			lambda_l=lambda_lh[i].first; lambda_h=lambda_lh[i].second;
-
-			gsl_root_fsolver_set (s_f, &F, lambda_l, lambda_h);
-
-			do {
-				iter++;
-				status=gsl_root_fsolver_iterate (s_f);
-				l=gsl_root_fsolver_root (s_f);
-				lambda_l=gsl_root_fsolver_x_lower (s_f);
-				lambda_h=gsl_root_fsolver_x_upper (s_f);
-				status=gsl_root_test_interval (lambda_l, lambda_h, 0, 1e-1);
-			}
-			while (status==GSL_CONTINUE && iter<max_iter);
-
-			iter=0;
-
-			gsl_root_fdfsolver_set (s_fdf, &FDF, l);
-
-			do {
-				iter++;
-				status=gsl_root_fdfsolver_iterate (s_fdf);
-				l_temp=l;
-				l=gsl_root_fdfsolver_root (s_fdf);
-				status=gsl_root_test_delta (l, l_temp, 0, 1e-5);
-			}
-			while (status==GSL_CONTINUE && iter<max_iter && l>l_min && l<l_max);
-
-			l=l_temp;
-			if (l<l_min) {l=l_min;}
-			if (l>l_max) {l=l_max;}
-			if (func_name=='R' || func_name=='r') {logf_l=LogRL_f (l, &params);} else {logf_l=LogL_f (l, &params);}
-
-			if (i==0) {logf=logf_l; lambda=l;}
-			else if (logf<logf_l) {logf=logf_l; lambda=l;}
-			else {}
+		  lambda_l=lambda_lh[i].first; lambda_h=lambda_lh[i].second;
+		  gsl_root_fsolver_set (s_f, &F, lambda_l, lambda_h);
+
+		  do {
+		    iter++;
+		    status=gsl_root_fsolver_iterate (s_f);
+		    l=gsl_root_fsolver_root (s_f);
+		    lambda_l=gsl_root_fsolver_x_lower (s_f);
+		    lambda_h=gsl_root_fsolver_x_upper (s_f);
+		    status=gsl_root_test_interval(lambda_l,lambda_h,0,1e-1);
+		  }
+		  while (status==GSL_CONTINUE && iter<max_iter);
+		  
+		  iter=0;
+		  
+		  gsl_root_fdfsolver_set (s_fdf, &FDF, l);
+		  
+		  do {
+		    iter++;
+		    status=gsl_root_fdfsolver_iterate (s_fdf);
+		    l_temp=l;
+		    l=gsl_root_fdfsolver_root (s_fdf);
+		    status=gsl_root_test_delta (l, l_temp, 0, 1e-5);
+		  }
+		  while (status==GSL_CONTINUE &&
+			 iter<max_iter &&
+			 l>l_min && l<l_max);
+
+		  l=l_temp;
+		  if (l<l_min) {l=l_min;}
+		  if (l>l_max) {l=l_max;}
+		  if (func_name=='R' || func_name=='r') {
+		    logf_l=LogRL_f (l, &params);
+		  } else {
+		    logf_l=LogL_f (l, &params);
+		  }
+		  
+		  if (i==0) {logf=logf_l; lambda=l;}
+		  else if (logf<logf_l) {logf=logf_l; lambda=l;}
+		  else {}
 		}
 		gsl_root_fsolver_free (s_f);
 		gsl_root_fdfsolver_free (s_fdf);
@@ -1971,14 +2058,17 @@ void CalcLambda (const char func_name, FUNC_PARAM &params, const double l_min, c
 	return;
 }
 
-
-
-
-
-//calculate lambda in the null model
-void CalcLambda (const char func_name, const gsl_vector *eval, const gsl_matrix *UtW, const gsl_vector *Uty, const double l_min, const double l_max, const size_t n_region, double &lambda, double &logl_H0)
-{
-	if (func_name!='R' && func_name!='L' && func_name!='r' && func_name!='l') {cout<<"func_name only takes 'R' or 'L': 'R' for log-restricted likelihood, 'L' for log-likelihood."<<endl; return;}
+// Calculate lambda in the null model.
+void CalcLambda (const char func_name, const gsl_vector *eval,
+		 const gsl_matrix *UtW, const gsl_vector *Uty,
+		 const double l_min, const double l_max,
+		 const size_t n_region, double &lambda, double &logl_H0) {
+	if (func_name!='R' && func_name!='L' && func_name!='r' &&
+	    func_name!='l') {
+	  cout<<"func_name only takes 'R' or 'L': 'R' for " <<
+	    "log-restricted likelihood, 'L' for log-likelihood." << endl;
+	  return;
+	}
 
 	size_t n_cvt=UtW->size2, ni_test=UtW->size1;
 	size_t n_index=(n_cvt+2+1)*(n_cvt+2)/2;
@@ -1988,10 +2078,6 @@ void CalcLambda (const char func_name, const gsl_vector *eval, const gsl_matrix
 
 	gsl_matrix_set_zero (Uab);
 	CalcUab (UtW, Uty, Uab);
-//	if (e_mode!=0) {
-//		gsl_vector_set_zero (ab);
-//		Calcab (W, y, ab);
-//	}
 
 	FUNC_PARAM param0={true, ni_test, n_cvt, eval, Uab, ab, 0};
 
@@ -2003,10 +2089,10 @@ void CalcLambda (const char func_name, const gsl_vector *eval, const gsl_matrix
 	return;
 }
 
-
-//obtain REMLE estimate for PVE using lambda_remle
-void CalcPve (const gsl_vector *eval, const gsl_matrix *UtW, const gsl_vector *Uty, const double lambda, const double trace_G, double &pve, double &pve_se)
-{
+// Obtain REMLE estimate for PVE using lambda_remle.
+void CalcPve (const gsl_vector *eval, const gsl_matrix *UtW,
+	      const gsl_vector *Uty, const double lambda,
+	      const double trace_G, double &pve, double &pve_se) {
 	size_t n_cvt=UtW->size2, ni_test=UtW->size1;
 	size_t n_index=(n_cvt+2+1)*(n_cvt+2)/2;
 
@@ -2015,10 +2101,6 @@ void CalcPve (const gsl_vector *eval, const gsl_matrix *UtW, const gsl_vector *U
 
 	gsl_matrix_set_zero (Uab);
 	CalcUab (UtW, Uty, Uab);
-	//	if (e_mode!=0) {
-	//		gsl_vector_set_zero (ab);
-	//		Calcab (W, y, ab);
-	//	}
 
 	FUNC_PARAM param0={true, ni_test, n_cvt, eval, Uab, ab, 0};
 
@@ -2032,11 +2114,13 @@ void CalcPve (const gsl_vector *eval, const gsl_matrix *UtW, const gsl_vector *U
 	return;
 }
 
-//obtain REML estimate for Vg and Ve using lambda_remle
-//obtain beta and se(beta) for coefficients
-//ab is not used when e_mode==0
-void CalcLmmVgVeBeta (const gsl_vector *eval, const gsl_matrix *UtW, const gsl_vector *Uty, const double lambda, double &vg, double &ve, gsl_vector *beta, gsl_vector *se_beta)
-{
+// Obtain REML estimate for Vg and Ve using lambda_remle.
+// Obtain beta and se(beta) for coefficients.
+// ab is not used when e_mode==0.
+void CalcLmmVgVeBeta (const gsl_vector *eval, const gsl_matrix *UtW,
+		      const gsl_vector *Uty, const double lambda,
+		      double &vg, double &ve, gsl_vector *beta,
+		      gsl_vector *se_beta) {
 	size_t n_cvt=UtW->size2, ni_test=UtW->size1;
 	size_t n_index=(n_cvt+2+1)*(n_cvt+2)/2;
 
@@ -2059,7 +2143,7 @@ void CalcLmmVgVeBeta (const gsl_vector *eval, const gsl_matrix *UtW, const gsl_v
 	gsl_vector_add_constant (v_temp, 1.0);
 	gsl_vector_div (Hi_eval, v_temp);
 
-	//calculate beta
+	// Calculate beta.
 	gsl_matrix_memcpy (HiW, UtW);
 	for (size_t i=0; i<UtW->size2; i++) {
 		gsl_vector_view HiW_col=gsl_matrix_column(HiW, i);
@@ -2074,7 +2158,7 @@ void CalcLmmVgVeBeta (const gsl_vector *eval, const gsl_matrix *UtW, const gsl_v
 	LUSolve (WHiW, pmt, WHiy, beta);
 	LUInvert (WHiW, pmt, Vbeta);
 
-	//calculate vg and ve
+	// Calculate vg and ve.
 	CalcPab (n_cvt, 0, Hi_eval, Uab, ab, Pab);
 
 	size_t index_yy=GetabIndex (n_cvt+2, n_cvt+2, n_cvt);
@@ -2083,12 +2167,12 @@ void CalcLmmVgVeBeta (const gsl_vector *eval, const gsl_matrix *UtW, const gsl_v
 	ve=P_yy/(double)(ni_test-n_cvt);
 	vg=ve*lambda;
 
-	//with ve, calculate se(beta)
+	// With ve, calculate se(beta).
 	gsl_matrix_scale(Vbeta, ve);
 
-	//obtain se_beta
+	// Obtain se_beta.
 	for (size_t i=0; i<Vbeta->size1; i++) {
-		gsl_vector_set (se_beta, i, sqrt(gsl_matrix_get(Vbeta, i, i) ) );
+		gsl_vector_set (se_beta, i, sqrt(gsl_matrix_get(Vbeta,i,i)));
 	}
 
 	gsl_matrix_free(Uab);
@@ -2105,29 +2189,28 @@ void CalcLmmVgVeBeta (const gsl_vector *eval, const gsl_matrix *UtW, const gsl_v
 	return;
 }
 
-
-
-
-
-
-
-void LMM::AnalyzeBimbamGXE (const gsl_matrix *U, const gsl_vector *eval, const gsl_matrix *UtW, const gsl_vector *Uty, const gsl_matrix *W, const gsl_vector *y, const gsl_vector *env)
-{
+void LMM::AnalyzeBimbamGXE (const gsl_matrix *U, const gsl_vector *eval,
+			    const gsl_matrix *UtW, const gsl_vector *Uty,
+			    const gsl_matrix *W, const gsl_vector *y,
+			    const gsl_vector *env) {
 	igzstream infile (file_geno.c_str(), igzstream::in);
-//	ifstream infile (file_geno.c_str(), ifstream::in);
-	if (!infile) {cout<<"error reading genotype file:"<<file_geno<<endl; return;}
+	if (!infile) {
+	  cout<<"error reading genotype file:"<<file_geno<<endl;
+	  return;
+	}
 
 	clock_t time_start=clock();
 
 	string line;
 	char *ch_ptr;
 
-	double lambda_mle=0, lambda_remle=0, beta=0, se=0, p_wald=0, p_lrt=0, p_score=0;
+	double lambda_mle=0, lambda_remle=0, beta=0, se=0, p_wald=0;
+	double p_lrt=0, p_score=0;
 	double logl_H1=0.0, logl_H0=0.0;
 	int n_miss, c_phen;
 	double geno, x_mean;
 
-	//Calculate basic quantities
+	// Calculate basic quantities.
 	size_t n_index=(n_cvt+2+2+1)*(n_cvt+2+2)/2;
 
 	gsl_vector *x=gsl_vector_alloc (U->size1);
@@ -2137,24 +2220,21 @@ void LMM::AnalyzeBimbamGXE (const gsl_matrix *U, const gsl_vector *eval, const g
 	gsl_vector *ab=gsl_vector_alloc (n_index);
 
 	gsl_matrix *UtW_expand=gsl_matrix_alloc (U->size1, UtW->size2+2);
-	gsl_matrix_view UtW_expand_mat=gsl_matrix_submatrix(UtW_expand, 0, 0, U->size1, UtW->size2);
+	gsl_matrix_view UtW_expand_mat=
+	  gsl_matrix_submatrix(UtW_expand, 0, 0, U->size1, UtW->size2);
 	gsl_matrix_memcpy (&UtW_expand_mat.matrix, UtW);
-	gsl_vector_view UtW_expand_env=gsl_matrix_column(UtW_expand, UtW->size2);
+	gsl_vector_view UtW_expand_env=
+	  gsl_matrix_column(UtW_expand, UtW->size2);
 	gsl_blas_dgemv (CblasTrans, 1.0, U, env, 0.0, &UtW_expand_env.vector);
-	gsl_vector_view UtW_expand_x=gsl_matrix_column(UtW_expand, UtW->size2+1);
-
-	//gsl_matrix_set_zero (Uab);
-	//	CalcUab (UtW, Uty, Uab);
-//	if (e_mode!=0) {
-//		gsl_vector_set_zero (ab);
-//		Calcab (W, y, ab);
-//	}
-
-	//start reading genotypes and analyze
+	gsl_vector_view UtW_expand_x=
+	  gsl_matrix_column(UtW_expand, UtW->size2+1);
+	
+	// Start reading genotypes and analyze.
 	for (size_t t=0; t<indicator_snp.size(); ++t) {
-//		if (t>1) {break;}
 		!safeGetline(infile, line).eof();
-		if (t%d_pace==0 || t==(ns_total-1)) {ProgressBar ("Reading SNPs  ", t, ns_total-1);}
+		if (t%d_pace==0 || t==(ns_total-1)) {
+		  ProgressBar ("Reading SNPs  ", t, ns_total-1);
+		}
 		if (indicator_snp[t]==0) {continue;}
 
 		ch_ptr=strtok ((char *)line.c_str(), " , \t");
@@ -2167,7 +2247,10 @@ void LMM::AnalyzeBimbamGXE (const gsl_matrix *U, const gsl_vector *eval, const g
 			ch_ptr=strtok (NULL, " , \t");
 			if (indicator_idv[i]==0) {continue;}
 
-			if (strcmp(ch_ptr, "NA")==0) {gsl_vector_set(x_miss, c_phen, 0.0); n_miss++;}
+			if (strcmp(ch_ptr, "NA")==0) {
+			  gsl_vector_set(x_miss, c_phen, 0.0);
+			  n_miss++;
+			}
 			else {
 				geno=atof(ch_ptr);
 
@@ -2181,17 +2264,19 @@ void LMM::AnalyzeBimbamGXE (const gsl_matrix *U, const gsl_vector *eval, const g
 		x_mean/=(double)(ni_test-n_miss);
 
 		for (size_t i=0; i<ni_test; ++i) {
-			if (gsl_vector_get (x_miss, i)==0) {gsl_vector_set(x, i, x_mean);}
+			if (gsl_vector_get (x_miss, i)==0) {
+			  gsl_vector_set(x, i, x_mean);
+			}
 			geno=gsl_vector_get(x, i);
 			if (x_mean>1) {
 				gsl_vector_set(x, i, 2-geno);
 			}
 		}
 
-
-		//calculate statistics
+		// Calculate statistics.
 		time_start=clock();
-		gsl_blas_dgemv (CblasTrans, 1.0, U, x, 0.0, &UtW_expand_x.vector);
+		gsl_blas_dgemv (CblasTrans, 1.0, U, x, 0.0,
+				&UtW_expand_x.vector);
 		gsl_vector_mul (x, env);
 		gsl_blas_dgemv (CblasTrans, 1.0, U, x, 0.0, Utx);
 		time_UtX+=(clock()-time_start)/(double(CLOCKS_PER_SEC)*60.0);
@@ -2201,29 +2286,29 @@ void LMM::AnalyzeBimbamGXE (const gsl_matrix *U, const gsl_vector *eval, const g
 
 		if (a_mode==2 || a_mode==4) {
 		  FUNC_PARAM param0={true, ni_test, n_cvt+2, eval, Uab, ab, 0};
-		  CalcLambda ('L', param0, l_min, l_max, n_region, lambda_mle, logl_H0);
+		  CalcLambda ('L', param0, l_min, l_max, n_region,
+			      lambda_mle, logl_H0);
 		}
 
 		CalcUab(UtW_expand, Uty, Utx, Uab);
-//		if (e_mode!=0) {
-//			Calcab (W, y, x, ab);
-//		}
 
 		time_start=clock();
 		FUNC_PARAM param1={false, ni_test, n_cvt+2, eval, Uab, ab, 0};
 
-		//3 is before 1
+		// 3 is before 1.
 		if (a_mode==3 || a_mode==4) {
 			CalcRLScore (l_mle_null, param1, beta, se, p_score);
 		}
 
 		if (a_mode==1 || a_mode==4) {
-			CalcLambda ('R', param1, l_min, l_max, n_region, lambda_remle, logl_H1);
+			CalcLambda ('R', param1, l_min, l_max, n_region,
+				    lambda_remle, logl_H1);
 			CalcRLWald (lambda_remle, param1, beta, se, p_wald);
 		}
 
 		if (a_mode==2 || a_mode==4) {
-			CalcLambda ('L', param1, l_min, l_max, n_region, lambda_mle, logl_H1);
+			CalcLambda ('L', param1, l_min, l_max, n_region,
+				    lambda_mle, logl_H1);
 			p_lrt=gsl_cdf_chisq_Q (2.0*(logl_H1-logl_H0), 1);
 		}
 
@@ -2231,10 +2316,11 @@ void LMM::AnalyzeBimbamGXE (const gsl_matrix *U, const gsl_vector *eval, const g
 
 		time_opt+=(clock()-time_start)/(double(CLOCKS_PER_SEC)*60.0);
 
-		//store summary data
-		SUMSTAT SNPs={beta, se, lambda_remle, lambda_mle, p_wald, p_lrt, p_score};
+		// Store summary data.
+		SUMSTAT SNPs={beta, se, lambda_remle, lambda_mle,
+			      p_wald, p_lrt, p_score};
 		sumStat.push_back(SNPs);
-    }
+	}
 	cout<<endl;
 
 	gsl_vector_free (x);
@@ -2251,14 +2337,10 @@ void LMM::AnalyzeBimbamGXE (const gsl_matrix *U, const gsl_vector *eval, const g
 	return;
 }
 
-
-
-
-
-
-
-void LMM::AnalyzePlinkGXE (const gsl_matrix *U, const gsl_vector *eval, const gsl_matrix *UtW, const gsl_vector *Uty, const gsl_matrix *W, const gsl_vector *y, const gsl_vector *env)
-{
+void LMM::AnalyzePlinkGXE (const gsl_matrix *U, const gsl_vector *eval,
+			   const gsl_matrix *UtW, const gsl_vector *Uty,
+			   const gsl_matrix *W, const gsl_vector *y,
+			   const gsl_vector *env) {
 	string file_bed=file_bfile+".bed";
 	ifstream infile (file_bed.c_str(), ios::binary);
 	if (!infile) {cout<<"error reading bed file:"<<file_bed<<endl; return;}
@@ -2268,12 +2350,13 @@ void LMM::AnalyzePlinkGXE (const gsl_matrix *U, const gsl_vector *eval, const gs
 	char ch[1];
 	bitset<8> b;
 
-	double lambda_mle=0, lambda_remle=0, beta=0, se=0, p_wald=0, p_lrt=0, p_score=0;
+	double lambda_mle=0, lambda_remle=0, beta=0, se=0, p_wald=0;
+	double p_lrt=0, p_score=0;
 	double logl_H1=0.0, logl_H0=0.0;
 	int n_bit, n_miss, ci_total, ci_test;
 	double geno, x_mean;
 
-	//Calculate basic quantities
+	// Calculate basic quantities.
 	size_t n_index=(n_cvt+2+2+1)*(n_cvt+2+2)/2;
 
 	gsl_vector *x=gsl_vector_alloc (U->size1);
@@ -2282,56 +2365,64 @@ void LMM::AnalyzePlinkGXE (const gsl_matrix *U, const gsl_vector *eval, const gs
 	gsl_vector *ab=gsl_vector_alloc (n_index);
 
 	gsl_matrix *UtW_expand=gsl_matrix_alloc (U->size1, UtW->size2+2);
-	gsl_matrix_view UtW_expand_mat=gsl_matrix_submatrix(UtW_expand, 0, 0, U->size1, UtW->size2);
+	gsl_matrix_view UtW_expand_mat=
+	  gsl_matrix_submatrix(UtW_expand, 0, 0, U->size1, UtW->size2);
 	gsl_matrix_memcpy (&UtW_expand_mat.matrix, UtW);
-	gsl_vector_view UtW_expand_env=gsl_matrix_column(UtW_expand, UtW->size2);
+	gsl_vector_view UtW_expand_env=
+	  gsl_matrix_column(UtW_expand, UtW->size2);
 	gsl_blas_dgemv (CblasTrans, 1.0, U, env, 0.0, &UtW_expand_env.vector);
-	gsl_vector_view UtW_expand_x=gsl_matrix_column(UtW_expand, UtW->size2+1);
-
-	//gsl_matrix_set_zero (Uab);
-	//CalcUab (UtW, Uty, Uab);
-//	if (e_mode!=0) {
-//		gsl_vector_set_zero (ab);
-//		Calcab (W, y, ab);
-//	}
+	gsl_vector_view UtW_expand_x=
+	  gsl_matrix_column(UtW_expand, UtW->size2+1);
 
-	//calculate n_bit and c, the number of bit for each snp
+	// Calculate n_bit and c, the number of bit for each SNP.
 	if (ni_total%4==0) {n_bit=ni_total/4;}
 	else {n_bit=ni_total/4+1; }
 
-	//print the first three majic numbers
+	// Print the first three magic numbers.
 	for (int i=0; i<3; ++i) {
 		infile.read(ch,1);
 		b=ch[0];
 	}
 
-
 	for (vector<SNPINFO>::size_type t=0; t<snpInfo.size(); ++t) {
-		if (t%d_pace==0 || t==snpInfo.size()-1) {ProgressBar ("Reading SNPs  ", t, snpInfo.size()-1);}
-		if (indicator_snp[t]==0) {continue;}
+	  if (t%d_pace==0 || t==snpInfo.size()-1) {
+	    ProgressBar ("Reading SNPs  ", t, snpInfo.size()-1);
+	  }
+	  if (indicator_snp[t]==0) {continue;}
 
-		infile.seekg(t*n_bit+3);		//n_bit, and 3 is the number of magic numbers
+	        // n_bit, and 3 is the number of magic numbers
+		infile.seekg(t*n_bit+3);
 
-		//read genotypes
+		// Read genotypes.
 		x_mean=0.0;	n_miss=0; ci_total=0; ci_test=0;
 		for (int i=0; i<n_bit; ++i) {
 			infile.read(ch,1);
 			b=ch[0];
-			for (size_t j=0; j<4; ++j) {                //minor allele homozygous: 2.0; major: 0.0;
-				if ((i==(n_bit-1)) && ci_total==(int)ni_total) {break;}
-				if (indicator_idv[ci_total]==0) {ci_total++; continue;}
-
-				if (b[2*j]==0) {
-					if (b[2*j+1]==0) {gsl_vector_set(x, ci_test, 2); x_mean+=2.0; }
-					else {gsl_vector_set(x, ci_test, 1); x_mean+=1.0; }
-				}
-				else {
-					if (b[2*j+1]==1) {gsl_vector_set(x, ci_test, 0); }
-					else {gsl_vector_set(x, ci_test, -9); n_miss++; }
-				}
 
-				ci_total++;
-				ci_test++;
+			// Minor allele homozygous: 2.0; major: 0.0.
+			for (size_t j=0; j<4; ++j) {                
+			  if ((i==(n_bit-1)) && ci_total==(int)ni_total) {
+			    break;
+			  }
+			  if (indicator_idv[ci_total]==0) {
+			    ci_total++;
+			    continue;
+			  }
+			  
+			  if (b[2*j]==0) {
+			    if (b[2*j+1]==0) {
+			      gsl_vector_set(x, ci_test, 2);
+			      x_mean+=2.0;
+			    }
+			    else {gsl_vector_set(x, ci_test, 1); x_mean+=1.0; }
+			  }
+			  else {
+			    if (b[2*j+1]==1) {gsl_vector_set(x, ci_test, 0); }
+			    else {gsl_vector_set(x, ci_test, -9); n_miss++; }
+			  }
+			  
+			  ci_total++;
+			  ci_test++;
 			}
 		}
 
@@ -2339,15 +2430,19 @@ void LMM::AnalyzePlinkGXE (const gsl_matrix *U, const gsl_vector *eval, const gs
 
 		for (size_t i=0; i<ni_test; ++i) {
 			geno=gsl_vector_get(x,i);
-			if (geno==-9) {gsl_vector_set(x, i, x_mean); geno=x_mean;}
+			if (geno==-9) {
+			  gsl_vector_set(x, i, x_mean);
+			  geno=x_mean;
+			}
 			if (x_mean>1) {
-				gsl_vector_set(x, i, 2-geno);
+			  gsl_vector_set(x, i, 2-geno);
 			}
 		}
 
-		//calculate statistics
+		// Calculate statistics.
 		time_start=clock();
-		gsl_blas_dgemv (CblasTrans, 1.0, U, x, 0.0, &UtW_expand_x.vector);
+		gsl_blas_dgemv (CblasTrans, 1.0, U, x, 0.0,
+				&UtW_expand_x.vector);
 		gsl_vector_mul (x, env);
 		gsl_blas_dgemv (CblasTrans, 1.0, U, x, 0.0, Utx);
 		time_UtX+=(clock()-time_start)/(double(CLOCKS_PER_SEC)*60.0);
@@ -2357,39 +2452,39 @@ void LMM::AnalyzePlinkGXE (const gsl_matrix *U, const gsl_vector *eval, const gs
 
 		if (a_mode==2 || a_mode==4) {
 		  FUNC_PARAM param0={true, ni_test, n_cvt+2, eval, Uab, ab, 0};
-		  CalcLambda ('L', param0, l_min, l_max, n_region, lambda_mle, logl_H0);
+		  CalcLambda ('L', param0, l_min, l_max, n_region,
+			      lambda_mle, logl_H0);
 		}
 
 		CalcUab(UtW_expand, Uty, Utx, Uab);
 
-//		if (e_mode!=0) {
-//			Calcab (W, y, x, ab);
-//		}
-
 		time_start=clock();
 		FUNC_PARAM param1={false, ni_test, n_cvt+2, eval, Uab, ab, 0};
 
-		//3 is before 1, for beta
+		// 3 is before 1, for beta.
 		if (a_mode==3 || a_mode==4) {
 			CalcRLScore (l_mle_null, param1, beta, se, p_score);
 		}
 
 		if (a_mode==1 || a_mode==4) {
-			CalcLambda ('R', param1, l_min, l_max, n_region, lambda_remle, logl_H1);
-			CalcRLWald (lambda_remle, param1, beta, se, p_wald);
+		  CalcLambda ('R', param1, l_min, l_max, n_region,
+			      lambda_remle, logl_H1);
+		  CalcRLWald (lambda_remle, param1, beta, se, p_wald);
 		}
 
 		if (a_mode==2 || a_mode==4) {
-			CalcLambda ('L', param1, l_min, l_max, n_region, lambda_mle, logl_H1);
-			p_lrt=gsl_cdf_chisq_Q (2.0*(logl_H1-logl_H0), 1);
+		  CalcLambda ('L', param1, l_min, l_max, n_region,
+			      lambda_mle, logl_H1);
+		  p_lrt=gsl_cdf_chisq_Q (2.0*(logl_H1-logl_H0), 1);
 		}
 
 		if (x_mean>1) {beta*=-1;}
 
 		time_opt+=(clock()-time_start)/(double(CLOCKS_PER_SEC)*60.0);
 
-		//store summary data
-		SUMSTAT SNPs={beta, se, lambda_remle, lambda_mle, p_wald, p_lrt, p_score};
+		// Store summary data.
+		SUMSTAT SNPs={beta, se, lambda_remle, lambda_mle, p_wald,
+			      p_lrt, p_score};
 		sumStat.push_back(SNPs);
     }
 	cout<<endl;
-- 
cgit v1.2.3