Removed FORCE_FLOAT from io.cpp.

1 files changed, 1161 insertions, 881 deletions

diff --git a/src/io.cpp b/src/io.cpp
index 40f1c16..3a4bc3c 100644
--- a/src/io.cpp
+++ b/src/io.cpp
@@ -27,8 +27,10 @@
 #include <set>
 #include <cstring>
 #include <cmath>
+#include <cstdint>
 #include <stdio.h>
 #include <stdlib.h>
+#include <assert.h>
 
 #include "gsl/gsl_vector.h"
 #include "gsl/gsl_matrix.h"
@@ -310,8 +312,10 @@ bool ReadFile_column (const string &file_pheno, vector<int> &indicator_idv,
 		if (strcmp(ch_ptr, "NA")==0) {
 		  indicator_idv.push_back(0);
 		  pheno.push_back(-9);
-		} // Pheno is different from pimass2.
+		} 
 		else {
+
+		  // Pheno is different from pimass2.
 		  p=atof(ch_ptr);
 		  indicator_idv.push_back(1);
 		  pheno.push_back(p);
@@ -487,15 +491,18 @@ bool ReadFile_bim (const string &file_bim, vector<SNPINFO> &snpInfo) {
 }
 
 // Read .fam file.
-bool ReadFile_fam (const string &file_fam, vector<vector<int> > &indicator_pheno, vector<vector<double> > &pheno, map<string, int> &mapID2num, const vector<size_t> &p_column)
-{
+bool ReadFile_fam (const string &file_fam,
+		   vector<vector<int> > &indicator_pheno,
+		   vector<vector<double> > &pheno,
+		   map<string, int> &mapID2num,
+		   const vector<size_t> &p_column) {
 	indicator_pheno.clear();
 	pheno.clear();
 	mapID2num.clear();
 
 	igzstream infile (file_fam.c_str(), igzstream::in);
-	//ifstream infile (file_fam.c_str(), ifstream::in);
-	if (!infile) {cout<<"error opening .fam file: "<<file_fam<<endl; return false;}
+	if (!infile) {
+	  cout<<"error opening .fam file: "<<file_fam<<endl; return false;}
 
 	string line;
 	char *ch_ptr;
@@ -528,12 +535,19 @@ bool ReadFile_fam (const string &file_fam, vector<vector<int> > &indicator_pheno
 		while (i<p_max ) {
 			if (mapP2c.count(i+1)!=0 ) {
 				if (strcmp(ch_ptr, "NA")==0) {
-					ind_pheno_row[mapP2c[i+1]]=0; pheno_row[mapP2c[i+1]]=-9;
+				  ind_pheno_row[mapP2c[i+1]]=0;
+				  pheno_row[mapP2c[i+1]]=-9;
 				} else {
-					p=atof(ch_ptr);
+				  p=atof(ch_ptr);
 
-					if (p==-9) {ind_pheno_row[mapP2c[i+1]]=0; pheno_row[mapP2c[i+1]]=-9;}
-					else {ind_pheno_row[mapP2c[i+1]]=1; pheno_row[mapP2c[i+1]]=p;}
+				  if (p==-9) {
+				    ind_pheno_row[mapP2c[i+1]]=0;
+				    pheno_row[mapP2c[i+1]]=-9;
+				  }
+				  else {
+				    ind_pheno_row[mapP2c[i+1]]=1;
+				    pheno_row[mapP2c[i+1]]=p;
+				  }
 				}
 			}
 			i++;
@@ -551,20 +565,26 @@ bool ReadFile_fam (const string &file_fam, vector<vector<int> > &indicator_pheno
 	return true;
 }
 
-
-
-
-
-
-//Read bimbam mean genotype file, the first time, to obtain #SNPs for analysis (ns_test) and total #SNP (ns_total)
-bool ReadFile_geno (const string &file_geno, const set<string> &setSnps, const gsl_matrix *W, vector<int> &indicator_idv, vector<int> &indicator_snp, const double &maf_level, const double &miss_level, const double &hwe_level, const double &r2_level, map<string, string> &mapRS2chr, map<string, long int> &mapRS2bp, map<string, double> &mapRS2cM, vector<SNPINFO> &snpInfo, size_t &ns_test)
-{
+// Read bimbam mean genotype file, the first time, to obtain #SNPs for
+// analysis (ns_test) and total #SNP (ns_total).
+bool ReadFile_geno (const string &file_geno, const set<string> &setSnps,
+		    const gsl_matrix *W, vector<int> &indicator_idv,
+		    vector<int> &indicator_snp, const double &maf_level,
+		    const double &miss_level, const double &hwe_level,
+		    const double &r2_level,
+		    map<string, string> &mapRS2chr,
+		    map<string, long int> &mapRS2bp,
+		    map<string, double> &mapRS2cM,
+		    vector<SNPINFO> &snpInfo,
+		    size_t &ns_test) {
 	indicator_snp.clear();
 	snpInfo.clear();
 
 	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 false;}
+	if (!infile) {
+	  cout<<"error reading genotype file:"<<file_geno<<endl;
+	  return false;
+	}
 
 	gsl_vector *genotype=gsl_vector_alloc (W->size1);
 	gsl_vector *genotype_miss=gsl_vector_alloc (W->size1);
@@ -575,7 +595,6 @@ bool ReadFile_geno (const string &file_geno, const set<string> &setSnps, const g
 	gsl_permutation * pmt=gsl_permutation_alloc (W->size2);
 
 	gsl_blas_dgemm(CblasTrans, CblasNoTrans, 1.0, W, W, 0.0, WtW);
-	//eigenlib_dgemm("T", "N", 1.0, W, W, 0.0, WtW);
 	int sig;
 	LUDecomp (WtW, pmt, &sig);
 	LUInvert (WtW, pmt, WtWi);
@@ -616,7 +635,8 @@ bool ReadFile_geno (const string &file_geno, const set<string> &setSnps, const g
 		major=ch_ptr;
 
 		if (setSnps.size()!=0 && setSnps.count(rs)==0) {
-		  SNPINFO sInfo={"-9", rs, -9, -9, minor, major, 0, -9, -9, 0, 0, file_pos};
+		  SNPINFO sInfo={"-9", rs, -9, -9, minor, major, 0, -9, -9,
+				 0, 0, file_pos};
 		  snpInfo.push_back(sInfo);
 		  indicator_snp.push_back(0);
 
@@ -634,7 +654,12 @@ bool ReadFile_geno (const string &file_geno, const set<string> &setSnps, const g
 			ch_ptr=strtok (NULL, " , \t");
 			if (indicator_idv[i]==0) {continue;}
 
-			if (strcmp(ch_ptr, "NA")==0) {gsl_vector_set (genotype_miss, c_idv, 1); n_miss++; c_idv++; continue;}
+			if (strcmp(ch_ptr, "NA")==0) {
+			  gsl_vector_set (genotype_miss, c_idv, 1);
+			  n_miss++;
+			  c_idv++;
+			  continue;
+			}
 
 			geno=atof(ch_ptr);
 			if (geno>=0 && geno<=0.5) {n_0++;}
@@ -643,8 +668,6 @@ bool ReadFile_geno (const string &file_geno, const set<string> &setSnps, const g
 
 			gsl_vector_set (genotype, c_idv, geno);
 
-//			if (geno<0) {n_miss++; continue;}
-
 			if (flag_poly==0) {geno_old=geno; flag_poly=2;}
 			if (flag_poly==2 && geno!=geno_old) {flag_poly=1;}
 
@@ -654,24 +677,38 @@ bool ReadFile_geno (const string &file_geno, const set<string> &setSnps, const g
 		}
 		maf/=2.0*(double)(ni_test-n_miss);
 
-		SNPINFO sInfo={chr, rs, cM, b_pos, minor, major, n_miss, (double)n_miss/(double)ni_test, maf, ni_test-n_miss, 0, file_pos};
+		SNPINFO sInfo={chr, rs, cM, b_pos, minor, major, n_miss,
+			       (double)n_miss/(double)ni_test, maf,
+			       ni_test-n_miss, 0, file_pos};
 		snpInfo.push_back(sInfo);
 		file_pos++;
 
-		if ( (double)n_miss/(double)ni_test > miss_level) {indicator_snp.push_back(0); continue;}
+		if ( (double)n_miss/(double)ni_test > miss_level) {
+		  indicator_snp.push_back(0);
+		  continue;
+		}
 
-		if ( (maf<maf_level || maf> (1.0-maf_level)) && maf_level!=-1 ) {indicator_snp.push_back(0); continue;}
+		if ((maf<maf_level || maf> (1.0-maf_level)) && maf_level!=-1) {
+		  indicator_snp.push_back(0);
+		  continue;
+		}
 
 		if (flag_poly!=1) {indicator_snp.push_back(0); continue;}
 
 		if (hwe_level!=0 && maf_level!=-1) {
-			if (CalcHWE(n_0, n_2, n_1)<hwe_level) {indicator_snp.push_back(0); continue;}
+			if (CalcHWE(n_0, n_2, n_1)<hwe_level) {
+			  indicator_snp.push_back(0);
+			  continue;
+			}
 		}
 
-		//filter SNP if it is correlated with W
-		//unless W has only one column, of 1s
+		// Filter SNP if it is correlated with W unless W has
+		// only one column, of 1s.
 		for (size_t i=0; i<genotype->size; ++i) {
-			if (gsl_vector_get (genotype_miss, i)==1) {geno=maf*2.0; gsl_vector_set (genotype, i, geno);}
+			if (gsl_vector_get (genotype_miss, i)==1) {
+			  geno=maf*2.0;
+			  gsl_vector_set (genotype, i, geno);
+			}
 		}
 
 		gsl_blas_dgemv (CblasTrans, 1.0, W, genotype, 0.0, Wtx);
@@ -679,7 +716,10 @@ bool ReadFile_geno (const string &file_geno, const set<string> &setSnps, const g
 		gsl_blas_ddot (genotype, genotype, &v_x);
 		gsl_blas_ddot (Wtx, WtWiWtx, &v_w);
 
-		if (W->size2!=1 && v_w/v_x >= r2_level) {indicator_snp.push_back(0); continue;}
+		if (W->size2!=1 && v_w/v_x >= r2_level) {
+		  indicator_snp.push_back(0);
+		  continue;
+		}
 
 		indicator_snp.push_back(1);
 		ns_test++;
@@ -699,19 +739,21 @@ bool ReadFile_geno (const string &file_geno, const set<string> &setSnps, const g
 	return true;
 }
 
-
-
-
-
-
-//Read bed file, the first time
-bool ReadFile_bed (const string &file_bed, const set<string> &setSnps, const gsl_matrix *W, vector<int> &indicator_idv, vector<int> &indicator_snp, vector<SNPINFO> &snpInfo, const double &maf_level, const double &miss_level, const double &hwe_level, const double &r2_level, size_t &ns_test)
-{
+// Read bed file, the first time.
+bool ReadFile_bed (const string &file_bed, const set<string> &setSnps,
+		   const gsl_matrix *W, vector<int> &indicator_idv,
+		   vector<int> &indicator_snp, vector<SNPINFO> &snpInfo,
+		   const double &maf_level, const double &miss_level,
+		   const double &hwe_level, const double &r2_level,
+		   size_t &ns_test) {
 	indicator_snp.clear();
 	size_t ns_total=snpInfo.size();
 
 	ifstream infile (file_bed.c_str(), ios::binary);
-	if (!infile) {cout<<"error reading bed file:"<<file_bed<<endl; return false;}
+	if (!infile) {
+	  cout<<"error reading bed file:"<<file_bed<<endl;
+	  return false;
+	}
 
 	gsl_vector *genotype=gsl_vector_alloc (W->size1);
 	gsl_vector *genotype_miss=gsl_vector_alloc (W->size1);
@@ -739,12 +781,12 @@ bool ReadFile_bed (const string &file_bed, const set<string> &setSnps, const gsl
 	}
 	ns_test=0;
 
-	//calculate n_bit and c, the number of bit for each snp
+	// Calculate n_bit and c, the number of bit for each snp.
 	size_t n_bit;
 	if (ni_total%4==0) {n_bit=ni_total/4;}
 	else {n_bit=ni_total/4+1;}
 
-	//ignore the first three majic numbers
+	// Ignore the first three magic numbers.
 	for (int i=0; i<3; ++i) {
 		infile.read(ch,1);
 		b=ch[0];
@@ -754,11 +796,14 @@ bool ReadFile_bed (const string &file_bed, const set<string> &setSnps, const gsl
 	size_t n_miss;
 	size_t n_0, n_1, n_2, c;
 
-	//start reading snps and doing association test
+	// Start reading snps and doing association test.
 	for (size_t t=0; t<ns_total; ++t) {
-	  infile.seekg(t*n_bit+3);		//n_bit, and 3 is the number of magic numbers
 
-		if (setSnps.size()!=0 && setSnps.count(snpInfo[t].rs_number)==0) {
+	  // n_bit, and 3 is the number of magic numbers.
+	  infile.seekg(t*n_bit+3); 
+
+		if (setSnps.size()!=0 &&
+		    setSnps.count(snpInfo[t].rs_number) == 0) {
 			snpInfo[t].n_miss=-9;
 			snpInfo[t].missingness=-9;
 			snpInfo[t].maf=-9;
@@ -767,24 +812,41 @@ bool ReadFile_bed (const string &file_bed, const set<string> &setSnps, const gsl
 			continue;
 		}
 
-		//read genotypes
+		// Read genotypes.
 		c=0; maf=0.0; n_miss=0; n_0=0; n_1=0; n_2=0;
 		c_idv=0; gsl_vector_set_zero (genotype_miss);
 		for (size_t 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;
+
+			// Minor allele homozygous: 2.0; major: 0.0;
+			for (size_t j=0; j<4; ++j) {
 				if ((i==(n_bit-1)) && c==ni_total) {break;}
 				if (indicator_idv[c]==0) {c++; continue;}
 				c++;
 
 				if (b[2*j]==0) {
-					if (b[2*j+1]==0) {gsl_vector_set(genotype, c_idv, 2.0); maf+=2.0; n_2++;}
-					else {gsl_vector_set(genotype, c_idv, 1.0); maf+=1.0; n_1++;}
+				  if (b[2*j+1]==0) {
+				    gsl_vector_set(genotype, c_idv, 2.0);
+				    maf+=2.0;
+				    n_2++;
+				  }
+				  else {
+				    gsl_vector_set(genotype, c_idv, 1.0);
+				    maf+=1.0;
+				    n_1++;
+				  }
 				}
 				else {
-					if (b[2*j+1]==1) {gsl_vector_set(genotype, c_idv, 0.0); maf+=0.0; n_0++;}
-					else {gsl_vector_set(genotype_miss, c_idv, 1); n_miss++; }
+				  if (b[2*j+1]==1) {
+				    gsl_vector_set(genotype, c_idv, 0.0);
+				    maf+=0.0;
+				    n_0++;
+				  }
+				  else {
+				    gsl_vector_set(genotype_miss, c_idv, 1);
+				    n_miss++;
+				  }
 				}
 				c_idv++;
 			}
@@ -798,20 +860,35 @@ bool ReadFile_bed (const string &file_bed, const set<string> &setSnps, const gsl
 		snpInfo[t].n_nb=0;
 		snpInfo[t].file_position=t;
 
-		if ( (double)n_miss/(double)ni_test > miss_level) {indicator_snp.push_back(0); continue;}
+		if ( (double)n_miss/(double)ni_test > miss_level) {
+		  indicator_snp.push_back(0);
+		  continue;
+		}
 
-		if ( (maf<maf_level || maf> (1.0-maf_level)) && maf_level!=-1 ) {indicator_snp.push_back(0); continue;}
+		if ((maf<maf_level || maf> (1.0-maf_level)) && maf_level!=-1) {
+		  indicator_snp.push_back(0);
+		  continue;
+		}
 
-		if ( (n_0+n_1)==0 || (n_1+n_2)==0 || (n_2+n_0)==0) {indicator_snp.push_back(0); continue;}
+		if ( (n_0+n_1)==0 || (n_1+n_2)==0 || (n_2+n_0)==0) {
+		  indicator_snp.push_back(0);
+		  continue;
+		}
 
 		if (hwe_level!=0 && maf_level!=-1) {
-			if (CalcHWE(n_0, n_2, n_1)<hwe_level) {indicator_snp.push_back(0); continue;}
+		  if (CalcHWE(n_0, n_2, n_1)<hwe_level) {
+		    indicator_snp.push_back(0);
+		    continue;
+		  }
 		}
 
-		//filter SNP if it is correlated with W
-		//unless W has only one column, of 1s
+		// Filter SNP if it is correlated with W unless W has
+		// only one column, of 1s.
 		for (size_t i=0; i<genotype->size; ++i) {
-			if (gsl_vector_get (genotype_miss, i)==1) {geno=maf*2.0; gsl_vector_set (genotype, i, geno);}
+			if (gsl_vector_get (genotype_miss, i)==1) {
+			  geno=maf*2.0;
+			  gsl_vector_set (genotype, i, geno);
+			}
 		}
 
 		gsl_blas_dgemv (CblasTrans, 1.0, W, genotype, 0.0, Wtx);
@@ -819,7 +896,10 @@ bool ReadFile_bed (const string &file_bed, const set<string> &setSnps, const gsl
 		gsl_blas_ddot (genotype, genotype, &v_x);
 		gsl_blas_ddot (Wtx, WtWiWtx, &v_w);
 
-		if (W->size2!=1 && v_w/v_x > r2_level) {indicator_snp.push_back(0); continue;}
+		if (W->size2!=1 && v_w/v_x > r2_level) {
+		  indicator_snp.push_back(0);
+		  continue;
+		}
 
 		indicator_snp.push_back(1);
 		ns_test++;
@@ -839,20 +919,14 @@ bool ReadFile_bed (const string &file_bed, const set<string> &setSnps, const gsl
 	return true;
 }
 
-
-
-
-
-//read the genotype for one SNP; remember to read empty lines
-//geno stores original genotypes without centering
-//missing values are replaced by mean
-bool Bimbam_ReadOneSNP (const size_t inc, const vector<int> &indicator_idv, igzstream &infile, gsl_vector *geno, double &geno_mean)
-{
+// Read the genotype for one SNP; remember to read empty lines.
+// Geno stores original genotypes without centering.
+// Missing values are replaced by mean.
+bool Bimbam_ReadOneSNP (const size_t inc, const vector<int> &indicator_idv,
+			igzstream &infile, gsl_vector *geno,
+			double &geno_mean) {
   size_t ni_total=indicator_idv.size();
 
-  //  if (infile.eof()) {infile.clear();}
-  //  infile.seekg(pos);
-
   string line;
   char *ch_ptr;
   bool flag=false;
@@ -896,16 +970,17 @@ bool Bimbam_ReadOneSNP (const size_t inc, const vector<int> &indicator_idv, igzs
   return flag;
 }
 
-
-//for plink, store SNPs as double too
-void Plink_ReadOneSNP (const int pos, const vector<int> &indicator_idv, ifstream &infile, gsl_vector *geno, double &geno_mean)
-{
+// For PLINK, store SNPs as double too.
+void Plink_ReadOneSNP (const int pos, const vector<int> &indicator_idv,
+		       ifstream &infile, gsl_vector *geno, double &geno_mean) {
   size_t ni_total=indicator_idv.size(), n_bit;
   if (ni_total%4==0) {n_bit=ni_total/4;}
   else {n_bit=ni_total/4+1;}
-  infile.seekg(pos*n_bit+3); //n_bit, and 3 is the number of magic numbers
 
-  //read genotypes
+  // n_bit, and 3 is the number of magic numbers.
+  infile.seekg(pos*n_bit+3); 
+
+  // Read genotypes.
   char ch[1];
   bitset<8> b;
 
@@ -916,7 +991,9 @@ void Plink_ReadOneSNP (const int pos, const vector<int> &indicator_idv, ifstream
   for (size_t 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;
+
+    // Minor allele homozygous: 2.0; major: 0.0.
+    for (size_t j=0; j<4; ++j) { 
       if ((i==(n_bit-1)) && c==ni_total) {break;}
       if (indicator_idv[c]==0) {c++; continue;}
       c++;
@@ -951,15 +1028,14 @@ void Plink_ReadOneSNP (const int pos, const vector<int> &indicator_idv, ifstream
   return;
 }
 
-
-
-
-
-void ReadFile_kin (const string &file_kin, vector<int> &indicator_idv, map<string, int> &mapID2num, const size_t k_mode, bool &error, gsl_matrix *G)
-{
+void ReadFile_kin (const string &file_kin, vector<int> &indicator_idv,
+		   map<string, int> &mapID2num, const size_t k_mode,
+		   bool &error, gsl_matrix *G) {
 	igzstream infile (file_kin.c_str(), igzstream::in);
-//	ifstream infile (file_kin.c_str(), ifstream::in);
-	if (!infile) {cout<<"error! fail to open kinship file: "<<file_kin<<endl; error=true; return;}
+	if (!infile) {
+	  cout<<"error! fail to open kinship file: "<<file_kin<<endl;
+	  error=true; return;
+	}
 
 	size_t ni_total=indicator_idv.size();
 
@@ -972,25 +1048,47 @@ void ReadFile_kin (const string &file_kin, vector<int> &indicator_idv, map<strin
 	if (k_mode==1) {
 		size_t i_test=0, i_total=0, j_test=0, j_total=0;
 		while (getline(infile, line)) {
-			if (i_total==ni_total) {cout<<"error! number of rows in the kinship file is larger than the number of phentypes."<<endl; error=true;}
+			if (i_total==ni_total) {
+			  cout<<"error! number of rows in the kinship "<<
+			    "file is larger than the number of phentypes."<<
+			    endl;
+			  error=true;
+			}
 
 			if (indicator_idv[i_total]==0) {i_total++; continue;}
 
 			j_total=0; j_test=0;
 			ch_ptr=strtok ((char *)line.c_str(), " , \t");
 			while (ch_ptr!=NULL) {
-				if (j_total==ni_total) {cout<<"error! number of columns in the kinship file is larger than the number of phentypes for row = "<<i_total<<endl; error=true;}
+				if (j_total==ni_total) {
+				  cout<<"error! number of columns in the "<<
+				    "kinship file is larger than the number"<<
+				    " of phentypes for row = "<<i_total<<endl;
+				  error=true;
+				}
 
 				d=atof(ch_ptr);
-				if (indicator_idv[j_total]==1) {gsl_matrix_set (G, i_test, j_test, d); j_test++;}
+				if (indicator_idv[j_total]==1) {
+				  gsl_matrix_set (G, i_test, j_test, d);
+				  j_test++;
+				}
 				j_total++;
 
 				ch_ptr=strtok (NULL, " , \t");
 			}
-			if (j_total!=ni_total) {cout<<"error! number of columns in the kinship file do not match the number of phentypes for row = "<<i_total<<endl; error=true;}
+			if (j_total!=ni_total) {
+			  cout<<"error! number of columns in the kinship "<<
+			    "file do not match the number of phentypes for "<<
+			    "row = "<<i_total<<endl;
+			  error=true;
+			}
 			i_total++; i_test++;
 		}
-		if (i_total!=ni_total) {cout<<"error! number of rows in the kinship file do not match the number of phentypes."<<endl; error=true;}
+		if (i_total!=ni_total) {
+		  cout<<"error! number of rows in the kinship file do "<<
+		    "not match the number of phentypes."<<endl;
+		  error=true;
+		}
 	}
 	else {
 		map<size_t, size_t> mapID2ID;
@@ -1010,14 +1108,24 @@ void ReadFile_kin (const string &file_kin, vector<int> &indicator_idv, map<strin
 			id2=ch_ptr;
 			ch_ptr=strtok (NULL, " , \t");
 			d=atof(ch_ptr);
-			if (mapID2num.count(id1)==0 || mapID2num.count(id2)==0) {continue;}
-			if (indicator_idv[mapID2num[id1]]==0 || indicator_idv[mapID2num[id2]]==0) {continue;}
+			if (mapID2num.count(id1)==0 ||
+			    mapID2num.count(id2)==0) {
+			  continue;
+			}
+			if (indicator_idv[mapID2num[id1]]==0 ||
+			    indicator_idv[mapID2num[id2]]==0) {
+			  continue;
+			}
 
 			n_id1=mapID2ID[mapID2num[id1]];
 			n_id2=mapID2ID[mapID2num[id2]];
 
 			Cov_d=gsl_matrix_get(G, n_id1, n_id2);
-			if (Cov_d!=0 && Cov_d!=d) {cout<<"error! redundant and unequal terms in the kinship file, for id1 = "<<id1<<" and id2 = "<<id2<<endl;}
+			if (Cov_d!=0 && Cov_d!=d) {
+			  cout<<"error! redundant and unequal terms in the "<<
+			    "kinship file, for id1 = "<<id1<<" and id2 = "<<
+			    id2<<endl;
+			}
 			else {
 				gsl_matrix_set(G, n_id1, n_id2, d);
 				gsl_matrix_set(G, n_id2, n_id1, d);
@@ -1031,19 +1139,24 @@ void ReadFile_kin (const string &file_kin, vector<int> &indicator_idv, map<strin
 	return;
 }
 
-
-void ReadFile_mk (const string &file_mk, vector<int> &indicator_idv, map<string, int> &mapID2num, const size_t k_mode, bool &error, gsl_matrix *G)
-{
+void ReadFile_mk (const string &file_mk, vector<int> &indicator_idv,
+		  map<string, int> &mapID2num, const size_t k_mode,
+		  bool &error, gsl_matrix *G) {
 	igzstream infile (file_mk.c_str(), igzstream::in);
-	if (!infile) {cout<<"error! fail to open file: "<<file_mk<<endl; error=true; return;}
+	if (!infile) {cout<<"error! fail to open file: "<<file_mk<<endl;
+	  error=true;
+	  return;
+	}
 
 	string file_kin, line;
 
 	size_t i=0;
 	while (getline(infile, line)) {
 	  file_kin=line.c_str();
-	  gsl_matrix_view G_sub=gsl_matrix_submatrix(G, 0, i*G->size1, G->size1, G->size1);
-	  ReadFile_kin (file_kin, indicator_idv, mapID2num, k_mode, error, &G_sub.matrix);
+	  gsl_matrix_view G_sub=gsl_matrix_submatrix(G, 0, i*G->size1,
+						     G->size1, G->size1);
+	  ReadFile_kin (file_kin, indicator_idv, mapID2num, k_mode,
+			error, &G_sub.matrix);
 	  i++;
 	}
 
@@ -1052,12 +1165,13 @@ void ReadFile_mk (const string &file_mk, vector<int> &indicator_idv, map<string,
 	return;
 }
 
-
-void ReadFile_eigenU (const string &file_ku, bool &error, gsl_matrix *U)
-{
+void ReadFile_eigenU (const string &file_ku, bool &error, gsl_matrix *U) {
 	igzstream infile (file_ku.c_str(), igzstream::in);
-//	ifstream infile (file_ku.c_str(), ifstream::in);
-	if (!infile) {cout<<"error! fail to open the U file: "<<file_ku<<endl; error=true; return;}
+	if (!infile) {
+	  cout<<"error! fail to open the U file: "<<file_ku<<endl;
+	  error=true;
+	  return;
+	}
 
 	size_t n_row=U->size1, n_col=U->size2, i_row=0, i_col=0;
 
@@ -1068,12 +1182,21 @@ void ReadFile_eigenU (const string &file_ku, bool &error, gsl_matrix *U)
 	double d;
 
 	while (getline(infile, line)) {
-		if (i_row==n_row) {cout<<"error! number of rows in the U file is larger than expected."<<endl; error=true;}
+		if (i_row==n_row) {
+		  cout<<"error! number of rows in the U file is larger "<<
+		    "than expected."<<endl;
+		  error=true;
+		}
 
 		i_col=0;
 		ch_ptr=strtok ((char *)line.c_str(), " , \t");
 		while (ch_ptr!=NULL) {
-			if (i_col==n_col) {cout<<"error! number of columns in the U file is larger than expected, for row = "<<i_row<<endl; error=true;}
+			if (i_col==n_col) {
+			  cout<<"error! number of columns in the U file "<<
+			    "is larger than expected, for row = "<<
+			    i_row<<endl;
+			  error=true;
+			}
 
 			d=atof(ch_ptr);
 			gsl_matrix_set (U, i_row, i_col, d);
@@ -1091,14 +1214,13 @@ void ReadFile_eigenU (const string &file_ku, bool &error, gsl_matrix *U)
 	return;
 }
 
-
-
-
-void ReadFile_eigenD (const string &file_kd, bool &error, gsl_vector *eval)
-{
+void ReadFile_eigenD (const string &file_kd, bool &error, gsl_vector *eval) {
 	igzstream infile (file_kd.c_str(), igzstream::in);
-//	ifstream infile (file_kd.c_str(), ifstream::in);
-	if (!infile) {cout<<"error! fail to open the D file: "<<file_kd<<endl; error=true; return;}
+	if (!infile) {
+	  cout<<"error! fail to open the D file: "<<file_kd<<endl;
+	  error=true;
+	  return;
+	}
 
 	size_t n_row=eval->size, i_row=0;
 
@@ -1109,13 +1231,21 @@ void ReadFile_eigenD (const string &file_kd, bool &error, gsl_vector *eval)
 	double d;
 
 	while (getline(infile, line)) {
-		if (i_row==n_row) {cout<<"error! number of rows in the D file is larger than expected."<<endl; error=true;}
+		if (i_row==n_row) {
+		  cout<<"error! number of rows in the D file is larger "<<
+		    "than expected."<<endl;
+		  error=true;
+		}
 
 		ch_ptr=strtok ((char *)line.c_str(), " , \t");
 		d=atof(ch_ptr);
 
 		ch_ptr=strtok (NULL, " , \t");
-		if (ch_ptr!=NULL) {cout<<"error! number of columns in the D file is larger than expected, for row = "<<i_row<<endl; error=true;}
+		if (ch_ptr!=NULL) {
+		  cout<<"error! number of columns in the D file is larger "<<
+		    "than expected, for row = "<<i_row<<endl;
+		  error=true;
+		}
 
 		gsl_vector_set (eval, i_row, d);
 
@@ -1128,14 +1258,15 @@ void ReadFile_eigenD (const string &file_kd, bool &error, gsl_vector *eval)
 	return;
 }
 
-
-
-//read bimbam mean genotype file and calculate kinship matrix
-bool BimbamKin (const string &file_geno, vector<int> &indicator_snp, const int k_mode, const int display_pace, gsl_matrix *matrix_kin)
-{
+// Read bimbam mean genotype file and calculate kinship matrix.
+bool BimbamKin (const string &file_geno, vector<int> &indicator_snp,
+		const int k_mode, const int display_pace,
+		gsl_matrix *matrix_kin) {
 	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 false;}
+	if (!infile) {
+	  cout<<"error reading genotype file:"<<file_geno<<endl;
+	  return false;
+	}
 
 	string line;
 	char *ch_ptr;
@@ -1147,7 +1278,7 @@ bool BimbamKin (const string &file_geno, vector<int> &indicator_snp, const int k
 	gsl_vector *geno=gsl_vector_alloc (ni_total);
 	gsl_vector *geno_miss=gsl_vector_alloc (ni_total);
 
-	//create a large matrix
+	// Create a large matrix.
 	size_t msize=10000;
 	gsl_matrix *Xlarge=gsl_matrix_alloc (ni_total, msize);
 	gsl_matrix_set_zero(Xlarge);
@@ -1155,7 +1286,9 @@ bool BimbamKin (const string &file_geno, vector<int> &indicator_snp, const int k
 	size_t ns_test=0;
 	for (size_t t=0; t<indicator_snp.size(); ++t) {
 		!safeGetline(infile, line).eof();
-		if (t%display_pace==0 || t==(indicator_snp.size()-1)) {ProgressBar ("Reading SNPs  ", t, indicator_snp.size()-1);}
+		if (t%display_pace==0 || t==(indicator_snp.size()-1)) {
+		  ProgressBar ("Reading SNPs  ", t, indicator_snp.size()-1);
+		}
 		if (indicator_snp[t]==0) {continue;}
 
 		ch_ptr=strtok ((char *)line.c_str(), " , \t");
@@ -1166,8 +1299,9 @@ bool BimbamKin (const string &file_geno, vector<int> &indicator_snp, const int k
 		gsl_vector_set_all(geno_miss, 0);
 		for (size_t i=0; i<ni_total; ++i) {
 			ch_ptr=strtok (NULL, " , \t");
-			if (strcmp(ch_ptr, "NA")==0) {gsl_vector_set(geno_miss, i, 0); n_miss++;}
-			else {
+			if (strcmp(ch_ptr, "NA")==0) {
+			  gsl_vector_set(geno_miss, i, 0); n_miss++;
+			} else {
 				d=atof(ch_ptr);
 				gsl_vector_set (geno, i, d);
 				gsl_vector_set (geno_miss, i, 1);
@@ -1180,36 +1314,27 @@ bool BimbamKin (const string &file_geno, vector<int> &indicator_snp, const int k
 		geno_var+=geno_mean*geno_mean*(double)n_miss;
 		geno_var/=(double)ni_total;
 		geno_var-=geno_mean*geno_mean;
-//		geno_var=geno_mean*(1-geno_mean*0.5);
 
 		for (size_t i=0; i<ni_total; ++i) {
-			if (gsl_vector_get (geno_miss, i)==0) {gsl_vector_set(geno, i, geno_mean);}
+			if (gsl_vector_get (geno_miss, i)==0) {
+			  gsl_vector_set(geno, i, geno_mean);
+			}
 		}
 
 		gsl_vector_add_constant (geno, -1.0*geno_mean);
 
-		/*
-		if (geno_var!=0) {
-		  if (k_mode==1) {
-		    gsl_blas_dsyr (CblasUpper, 1.0, geno, matrix_kin);
-		    //eigenlib_dsyr (1.0, geno, matrix_kin);
-		  } else if (k_mode==2) {
-		    gsl_blas_dsyr (CblasUpper, 1.0/geno_var, geno, matrix_kin);
-		    //eigenlib_dsyr (1.0/geno_var, geno, matrix_kin);
-		  } else {
-		    cout<<"Unknown kinship mode."<<endl;
-		  }
+		if (k_mode==2 && geno_var!=0) {
+		  gsl_vector_scale (geno, 1.0/sqrt(geno_var));
 		}
-		*/
-
-		if (k_mode==2 && geno_var!=0) {gsl_vector_scale (geno, 1.0/sqrt(geno_var));}
-		gsl_vector_view Xlarge_col=gsl_matrix_column (Xlarge, ns_test%msize);
+		gsl_vector_view Xlarge_col=
+		  gsl_matrix_column (Xlarge, ns_test%msize);
 		gsl_vector_memcpy (&Xlarge_col.vector, geno);
 
 		ns_test++;
 
 		if (ns_test%msize==0) {
-		  eigenlib_dgemm ("N", "T", 1.0, Xlarge, Xlarge, 1.0, matrix_kin);
+		  eigenlib_dgemm ("N", "T", 1.0, Xlarge, Xlarge, 1.0,
+				  matrix_kin);
 		  gsl_matrix_set_zero(Xlarge);
 		}
 	}
@@ -1238,16 +1363,14 @@ bool BimbamKin (const string &file_geno, vector<int> &indicator_snp, const int k
 	return true;
 }
 
-
-
-
-
-
-
-bool PlinkKin (const string &file_bed, vector<int> &indicator_snp, const int k_mode, const int display_pace, gsl_matrix *matrix_kin)
-{
+bool PlinkKin (const string &file_bed, vector<int> &indicator_snp,
+	       const int k_mode, const int display_pace,
+	       gsl_matrix *matrix_kin) {
 	ifstream infile (file_bed.c_str(), ios::binary);
-	if (!infile) {cout<<"error reading bed file:"<<file_bed<<endl; return false;}
+	if (!infile) {
+	  cout<<"error reading bed file:"<<file_bed<<endl;
+	  return false;
+	}
 
 	char ch[1];
 	bitset<8> b;
@@ -1261,12 +1384,12 @@ bool PlinkKin (const string &file_bed, vector<int> &indicator_snp, const int k_m
 	size_t ns_test=0;
 	int n_bit;
 
-	//create a large matrix
+	// Create a large matrix.
 	size_t msize=10000;
 	gsl_matrix *Xlarge=gsl_matrix_alloc (ni_total, msize);
 	gsl_matrix_set_zero(Xlarge);
 
-	//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; }
 
@@ -1277,26 +1400,46 @@ bool PlinkKin (const string &file_bed, vector<int> &indicator_snp, const int k_m
 	}
 
 	for (size_t t=0; t<indicator_snp.size(); ++t) {
-		if (t%display_pace==0 || t==(indicator_snp.size()-1)) {ProgressBar ("Reading SNPs  ", t, indicator_snp.size()-1);}
+		if (t%display_pace==0 || t==(indicator_snp.size()-1)) {
+		  ProgressBar ("Reading SNPs  ", t, indicator_snp.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.
 		geno_mean=0.0;	n_miss=0; ci_total=0; geno_var=0.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==ni_total) {break;}
+
+			// Minor allele homozygous: 2.0; major: 0.0.
+			for (size_t j=0; j<4; ++j) {                
+				if ((i==(n_bit-1)) && ci_total==ni_total) {
+				  break;
+				}
 
 				if (b[2*j]==0) {
-					if (b[2*j+1]==0) {gsl_vector_set(geno, ci_total, 2.0); geno_mean+=2.0; geno_var+=4.0; }
-					else {gsl_vector_set(geno, ci_total, 1.0); geno_mean+=1.0; geno_var+=1.0;}
+					if (b[2*j+1]==0) {
+					  gsl_vector_set(geno, ci_total, 2.0);
+					  geno_mean+=2.0;
+					  geno_var+=4.0;
+					}
+					else {
+					  gsl_vector_set(geno, ci_total, 1.0);
+					  geno_mean+=1.0;
+					  geno_var+=1.0;
+					}
 				}
 				else {
-					if (b[2*j+1]==1) {gsl_vector_set(geno, ci_total, 0.0); }
-					else {gsl_vector_set(geno, ci_total, -9.0); n_miss++; }
+					if (b[2*j+1]==1) {
+					  gsl_vector_set(geno,ci_total,0.0);
+					}
+					else {
+					  gsl_vector_set(geno,ci_total,-9.0);
+					  n_miss++;
+					}
 				}
 
 				ci_total++;
@@ -1307,7 +1450,6 @@ bool PlinkKin (const string &file_bed, vector<int> &indicator_snp, const int k_m
 		geno_var+=geno_mean*geno_mean*(double)n_miss;
 		geno_var/=(double)ni_total;
 		geno_var-=geno_mean*geno_mean;
-//		geno_var=geno_mean*(1-geno_mean*0.5);
 
 		for (size_t i=0; i<ni_total; ++i) {
 			d=gsl_vector_get(geno,i);
@@ -1316,22 +1458,17 @@ bool PlinkKin (const string &file_bed, vector<int> &indicator_snp, const int k_m
 
 		gsl_vector_add_constant (geno, -1.0*geno_mean);
 
-		/*
-		if (geno_var!=0) {
-			if (k_mode==1) {gsl_blas_dsyr (CblasUpper, 1.0, geno, matrix_kin);}
-			else if (k_mode==2) {gsl_blas_dsyr (CblasUpper, 1.0/geno_var, geno, matrix_kin);}
-			else {cout<<"Unknown kinship mode."<<endl;}
+		if (k_mode==2 && geno_var!=0) {
+		  gsl_vector_scale (geno, 1.0/sqrt(geno_var));
 		}
-		*/
-
-		if (k_mode==2 && geno_var!=0) {gsl_vector_scale (geno, 1.0/sqrt(geno_var));}
-		gsl_vector_view Xlarge_col=gsl_matrix_column (Xlarge, ns_test%msize);
+		gsl_vector_view Xlarge_col=
+		  gsl_matrix_column (Xlarge, ns_test%msize);
 		gsl_vector_memcpy (&Xlarge_col.vector, geno);
 
 		ns_test++;
 
 		if (ns_test%msize==0) {
-		  eigenlib_dgemm ("N", "T", 1.0, Xlarge, Xlarge, 1.0, matrix_kin);
+		  eigenlib_dgemm("N","T",1.0,Xlarge,Xlarge,1.0,matrix_kin);
 		  gsl_matrix_set_zero(Xlarge);
 		}
 	}
@@ -1360,16 +1497,16 @@ bool PlinkKin (const string &file_bed, vector<int> &indicator_snp, const int k_m
 	return true;
 }
 
-
-
-
-
-//Read bimbam mean genotype file, the second time, recode "mean" genotype and calculate K
-bool ReadFile_geno (const string &file_geno, vector<int> &indicator_idv, vector<int> &indicator_snp, gsl_matrix *UtX, gsl_matrix *K, const bool calc_K)
-{
+// Read bimbam mean genotype file, the second time, recode "mean"
+// genotype and calculate K.
+bool ReadFile_geno (const string &file_geno, vector<int> &indicator_idv,
+		    vector<int> &indicator_snp, gsl_matrix *UtX,
+		    gsl_matrix *K, const bool calc_K) {
 	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 false;}
+	if (!infile) {
+	  cout<<"error reading genotype file:"<<file_geno<<endl;
+	  return false;
+	}
 
 	string line;
 	char *ch_ptr;
@@ -1402,8 +1539,10 @@ bool ReadFile_geno (const string &file_geno, vector<int> &indicator_idv, vector<
 			ch_ptr=strtok (NULL, " , \t");
 			if (indicator_idv[j]==0) {continue;}
 
-			if (strcmp(ch_ptr, "NA")==0) {gsl_vector_set (genotype_miss, c_idv, 1); n_miss++;}
-			else {
+			if (strcmp(ch_ptr, "NA")==0) {
+			  gsl_vector_set (genotype_miss, c_idv, 1);
+			  n_miss++;
+			} else {
 				geno=atof(ch_ptr);
 				gsl_vector_set (genotype, c_idv, geno);
 				geno_mean+=geno;
@@ -1414,14 +1553,21 @@ bool ReadFile_geno (const string &file_geno, vector<int> &indicator_idv, vector<
 		geno_mean/=(double)(ni_test-n_miss);
 
 		for (size_t i=0; i<genotype->size; ++i) {
-			if (gsl_vector_get (genotype_miss, i)==1) {geno=0;}
-			else {geno=gsl_vector_get (genotype, i); geno-=geno_mean;}
+			if (gsl_vector_get (genotype_miss, i)==1) {
+			  geno=0;
+			}
+			else {
+			  geno=gsl_vector_get (genotype, i);
+			  geno-=geno_mean;
+			}
 
 			gsl_vector_set (genotype, i, geno);
 			gsl_matrix_set (UtX, i, c_snp, geno);
 		}
 
-		if (calc_K==true) {gsl_blas_dsyr (CblasUpper, 1.0, genotype, K);}
+		if (calc_K==true) {
+		  gsl_blas_dsyr (CblasUpper, 1.0, genotype, K);
+		}
 
 		c_snp++;
 	}
@@ -1446,14 +1592,18 @@ bool ReadFile_geno (const string &file_geno, vector<int> &indicator_idv, vector<
 	return true;
 }
 
-
-
-//compact version of the above function, using uchar instead of gsl_matrix
-bool ReadFile_geno (const string &file_geno, vector<int> &indicator_idv, vector<int> &indicator_snp, vector<vector<unsigned char> > &Xt, gsl_matrix *K, const bool calc_K, const size_t ni_test, const size_t ns_test)
-{
+// Compact version of the above function, using uchar instead of
+// gsl_matrix.
+bool ReadFile_geno (const string &file_geno, vector<int> &indicator_idv,
+		    vector<int> &indicator_snp,
+		    vector<vector<unsigned char> > &Xt,
+		    gsl_matrix *K, const bool calc_K, const size_t ni_test,
+		    const size_t ns_test) {
 	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 false;}
+	if (!infile) {
+	  cout<<"error reading genotype file:"<<file_geno<<endl;
+	  return false;
+	}
 
 	Xt.clear();
 	vector<unsigned char> Xt_row;
@@ -1490,7 +1640,10 @@ bool ReadFile_geno (const string &file_geno, vector<int> &indicator_idv, vector<
 			ch_ptr=strtok (NULL, " , \t");
 			if (indicator_idv[j]==0) {continue;}
 
-			if (strcmp(ch_ptr, "NA")==0) {gsl_vector_set (genotype_miss, c_idv, 1); n_miss++;} else {
+			if (strcmp(ch_ptr, "NA")==0) {
+			  gsl_vector_set (genotype_miss, c_idv, 1);
+			  n_miss++;
+			} else {
 				geno=atof(ch_ptr);
 				gsl_vector_set (genotype, c_idv, geno);
 				geno_mean+=geno;
@@ -1512,7 +1665,9 @@ bool ReadFile_geno (const string &file_geno, vector<int> &indicator_idv, vector<
 		}
 		Xt.push_back(Xt_row);
 
-		if (calc_K==true) {gsl_blas_dsyr (CblasUpper, 1.0, genotype, K);}
+		if (calc_K==true) {
+		  gsl_blas_dsyr (CblasUpper, 1.0, genotype, K);
+		}
 
 		c_snp++;
 	}
@@ -1537,14 +1692,16 @@ bool ReadFile_geno (const string &file_geno, vector<int> &indicator_idv, vector<
 	return true;
 }
 
-
-
-
-//Read bimbam mean genotype file, the second time, recode "mean" genotype and calculate K
-bool ReadFile_bed (const string &file_bed, vector<int> &indicator_idv, vector<int> &indicator_snp, gsl_matrix *UtX, gsl_matrix *K, const bool calc_K)
-{
+// Read bimbam mean genotype file, the second time, recode "mean"
+// genotype and calculate K.
+bool ReadFile_bed (const string &file_bed, vector<int> &indicator_idv,
+		   vector<int> &indicator_snp, gsl_matrix *UtX,
+		   gsl_matrix *K, const bool calc_K) {
 	ifstream infile (file_bed.c_str(), ios::binary);
-	if (!infile) {cout<<"error reading bed file:"<<file_bed<<endl; return false;}
+	if (!infile) {
+	  cout<<"error reading bed file:"<<file_bed<<endl;
+	  return false;
+	}
 
 	char ch[1];
 	bitset<8> b;
@@ -1558,7 +1715,7 @@ bool ReadFile_bed (const string &file_bed, vector<int> &indicator_idv, vector<in
 	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];
@@ -1572,28 +1729,44 @@ bool ReadFile_bed (const string &file_bed, vector<int> &indicator_idv, vector<in
 	size_t n_miss;
 	size_t c_idv=0, c_snp=0, c=0;
 
-	//start reading snps and doing association test
+	// Start reading snps and doing association test.
 	for (size_t t=0; t<ns_total; ++t) {
 		if (indicator_snp[t]==0) {continue;}
-		infile.seekg(t*n_bit+3);		//n_bit, and 3 is the number of magic numbers
 
-		//read genotypes
+		// n_bit, and 3 is the number of magic numbers.
+		infile.seekg(t*n_bit+3);		
+
+		// Read genotypes.
 		c_idv=0; geno_mean=0.0; n_miss=0; c=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;
+
+			// Minor allele homozygous: 2.0; major: 0.0.
+			for (size_t j=0; j<4; ++j) {
 			  if ((i==(n_bit-1)) && c==ni_total) {break;}
 				if (indicator_idv[c]==0) {c++; continue;}
 				c++;
 
 				if (b[2*j]==0) {
-					if (b[2*j+1]==0) {gsl_vector_set(genotype, c_idv, 2.0); geno_mean+=2.0;}
-					else {gsl_vector_set(genotype, c_idv, 1.0); geno_mean+=1.0;}
+				  if (b[2*j+1]==0) {
+				    gsl_vector_set(genotype, c_idv, 2.0);
+				    geno_mean+=2.0;
+				  }
+				  else {
+				    gsl_vector_set(genotype, c_idv, 1.0);
+				    geno_mean+=1.0;
+				  }
 				}
 				else {
-					if (b[2*j+1]==1) {gsl_vector_set(genotype, c_idv, 0.0); geno_mean+=0.0;}
-					else {gsl_vector_set(genotype, c_idv, -9.0); n_miss++;}
+				  if (b[2*j+1]==1) {
+				    gsl_vector_set(genotype, c_idv, 0.0);
+				    geno_mean+=0.0;
+				  }
+				  else {
+				    gsl_vector_set(genotype, c_idv, -9.0);
+				    n_miss++;
+				  }
 				}
 				c_idv++;
 			}
@@ -1610,7 +1783,9 @@ bool ReadFile_bed (const string &file_bed, vector<int> &indicator_idv, vector<in
 			gsl_matrix_set (UtX, i, c_snp, geno);
 		}
 
-		if (calc_K==true) {gsl_blas_dsyr (CblasUpper, 1.0, genotype, K);}
+		if (calc_K==true) {
+		  gsl_blas_dsyr (CblasUpper, 1.0, genotype, K);
+		}
 
 		c_snp++;
 	}
@@ -1633,14 +1808,17 @@ bool ReadFile_bed (const string &file_bed, vector<int> &indicator_idv, vector<in
 	return true;
 }
 
-
-
-
-//compact version of the above function, using uchar instead of gsl_matrix
-bool ReadFile_bed (const string &file_bed, vector<int> &indicator_idv, vector<int> &indicator_snp, vector<vector<unsigned char> > &Xt, gsl_matrix *K, const bool calc_K, const size_t ni_test, const size_t ns_test)
-{
+// Compact version of the above function, using uchar instead of gsl_matrix.
+bool ReadFile_bed (const string &file_bed, vector<int> &indicator_idv,
+		   vector<int> &indicator_snp,
+		   vector<vector<unsigned char> > &Xt, gsl_matrix *K,
+		   const bool calc_K, const size_t ni_test,
+		   const size_t ns_test) {
 	ifstream infile (file_bed.c_str(), ios::binary);
-	if (!infile) {cout<<"error reading bed file:"<<file_bed<<endl; return false;}
+	if (!infile) {
+	  cout<<"error reading bed file:"<<file_bed<<endl;
+	  return false;
+	}
 
 	Xt.clear();
 	vector<unsigned char> Xt_row;
@@ -1658,7 +1836,7 @@ bool ReadFile_bed (const string &file_bed, vector<int> &indicator_idv, vector<in
 	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];
@@ -1672,28 +1850,44 @@ bool ReadFile_bed (const string &file_bed, vector<int> &indicator_idv, vector<in
 	size_t n_miss;
 	size_t c_idv=0, c_snp=0, c=0;
 
-	//start reading snps and doing association test
+	// Start reading SNPs and doing association test.
 	for (size_t t=0; t<ns_total; ++t) {
 		if (indicator_snp[t]==0) {continue;}
-		infile.seekg(t*n_bit+3);		//n_bit, and 3 is the number of magic numbers
 
-		//read genotypes
+		// n_bit, and 3 is the number of magic numbers.
+		infile.seekg(t*n_bit+3);		
+
+		// Read genotypes.
 		c_idv=0; geno_mean=0.0; n_miss=0; c=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;
+
+			// Minor allele homozygous: 2.0; major: 0.0.
+			for (size_t j=0; j<4; ++j) {                
 			  if ((i==(n_bit-1)) && c==ni_total) {break;}
 				if (indicator_idv[c]==0) {c++; continue;}
 				c++;
 
 				if (b[2*j]==0) {
-					if (b[2*j+1]==0) {gsl_vector_set(genotype, c_idv, 2.0); geno_mean+=2.0;}
-					else {gsl_vector_set(genotype, c_idv, 1.0); geno_mean+=1.0;}
+				  if (b[2*j+1]==0) {
+				    gsl_vector_set(genotype, c_idv, 2.0);
+				    geno_mean+=2.0;
+				  }
+				  else {
+				    gsl_vector_set(genotype, c_idv, 1.0);
+				    geno_mean+=1.0;
+				  }
 				}
 				else {
-					if (b[2*j+1]==1) {gsl_vector_set(genotype, c_idv, 0.0); geno_mean+=0.0;}
-					else {gsl_vector_set(genotype, c_idv, -9.0); n_miss++;}
+				  if (b[2*j+1]==1) {
+				    gsl_vector_set(genotype, c_idv, 0.0);
+				    geno_mean+=0.0;
+				  }
+				  else {
+				    gsl_vector_set(genotype, c_idv, -9.0);
+				    n_miss++;
+				  }
 				}
 				c_idv++;
 			}
@@ -1713,7 +1907,9 @@ bool ReadFile_bed (const string &file_bed, vector<int> &indicator_idv, vector<in
 		}
 		Xt.push_back(Xt_row);
 
-		if (calc_K==true) {gsl_blas_dsyr (CblasUpper, 1.0, genotype, K);}
+		if (calc_K==true) {
+		  gsl_blas_dsyr (CblasUpper, 1.0, genotype, K);
+		}
 
 		c_snp++;
 	}
@@ -1736,18 +1932,15 @@ bool ReadFile_bed (const string &file_bed, vector<int> &indicator_idv, vector<in
 	return true;
 }
 
-
-
-
-
-
-
-bool ReadFile_est (const string &file_est, const vector<size_t> &est_column, map<string, double> &mapRS2est)
-{
+bool ReadFile_est (const string &file_est, const vector<size_t> &est_column,
+		   map<string, double> &mapRS2est) {
 	mapRS2est.clear();
 
 	ifstream infile (file_est.c_str(), ifstream::in);
-	if (!infile) {cout<<"error opening estimated parameter file: "<<file_est<<endl; return false;}
+	if (!infile) {
+	  cout<<"error opening estimated parameter file: "<<file_est<<endl;
+	  return false;
+	}
 
 	string line;
 	char *ch_ptr;
@@ -1755,7 +1948,7 @@ bool ReadFile_est (const string &file_est, const vector<size_t> &est_column, map
 	string rs;
 	double alpha, beta, gamma, d;
 
-	//header
+	// Header.
 	getline(infile, line);
 
 	size_t n=*max_element(est_column.begin(), est_column.end());
@@ -1778,7 +1971,9 @@ bool ReadFile_est (const string &file_est, const vector<size_t> &est_column, map
 			mapRS2est[rs]=d;
 		}
 		else {
-			cout<<"the same SNP occurs more than once in estimated parameter file: "<<rs<<endl; return false;
+		  cout << "the same SNP occurs more than once in estimated "<<
+		    "parameter file: "<<rs<<endl;
+		  return false;
 		}
 	}
 
@@ -1787,13 +1982,12 @@ bool ReadFile_est (const string &file_est, const vector<size_t> &est_column, map
 	return true;
 }
 
-
-
-bool CountFileLines (const string &file_input, size_t &n_lines)
-{
+bool CountFileLines (const string &file_input, size_t &n_lines) {
 	igzstream infile (file_input.c_str(), igzstream::in);
-	//ifstream infile (file_input.c_str(), ifstream::in);
-	if (!infile) {cout<<"error! fail to open file: "<<file_input<<endl; return false;}
+	if (!infile) {
+	  cout<<"error! fail to open file: "<<file_input<<endl;
+	  return false;
+	}
 
 	n_lines=count(istreambuf_iterator<char>(infile), istreambuf_iterator<char>(), '\n');
 	infile.seekg (0, ios::beg);
@@ -1801,16 +1995,17 @@ bool CountFileLines (const string &file_input, size_t &n_lines)
 	return true;
 }
 
-
-
-//Read gene expression file
-bool ReadFile_gene (const string &file_gene, vector<double> &vec_read, vector<SNPINFO> &snpInfo, size_t &ng_total)
-{
+// Read gene expression file.
+bool ReadFile_gene (const string &file_gene, vector<double> &vec_read,
+		    vector<SNPINFO> &snpInfo, size_t &ng_total) {
 	vec_read.clear();
 	ng_total=0;
 
 	igzstream infile (file_gene.c_str(), igzstream::in);
-	if (!infile) {cout<<"error! fail to open gene expression file: "<<file_gene<<endl; return false;}
+	if (!infile) {
+	  cout<<"error! fail to open gene expression file: "<<file_gene<<endl;
+	  return false;
+	}
 
 	string line;
 	char *ch_ptr;
@@ -1818,7 +2013,7 @@ bool ReadFile_gene (const string &file_gene, vector<double> &vec_read, vector<SN
 
 	size_t n_idv=0, t=0;
 
-	//header
+	// Header.
 	getline(infile, line);
 
 	while (getline(infile, line)) {
@@ -1841,9 +2036,13 @@ bool ReadFile_gene (const string &file_gene, vector<double> &vec_read, vector<SN
 			ch_ptr=strtok (NULL, " , \t");
 		}
 
-		if (t!=n_idv) {cout<<"error! number of columns doesn't match in row: "<<ng_total<<endl; return false;}
+		if (t!=n_idv) {
+		  cout<<"error! number of columns doesn't match in row: "<<
+		    ng_total<<endl;
+		  return false;
+		}
 
-		SNPINFO sInfo={"-9", rs, -9, -9, "-9", "-9", 0, -9, -9, 0, 0, 0};
+		SNPINFO sInfo={"-9",rs,-9,-9,"-9","-9",0,-9,-9,0,0,0};
 		snpInfo.push_back(sInfo);
 
 		ng_total++;
@@ -1855,28 +2054,28 @@ bool ReadFile_gene (const string &file_gene, vector<double> &vec_read, vector<SN
 	return true;
 }
 
-
-
-
-
-
-
 // WJA Added
-//Read Oxford sample file
-bool ReadFile_sample(const string &file_sample, vector<vector<int> > &indicator_pheno, vector<vector<double> > &pheno, const vector<size_t> &p_column, vector<int> &indicator_cvt, vector<vector<double> > &cvt, size_t &n_cvt)
-{
+// Read Oxford sample file.
+bool ReadFile_sample (const string &file_sample,
+		      vector<vector<int> > &indicator_pheno,
+		      vector<vector<double> > &pheno,
+		      const vector<size_t> &p_column,
+		      vector<int> &indicator_cvt,
+		      vector<vector<double> > &cvt, size_t &n_cvt) {
 	indicator_pheno.clear();
 	pheno.clear();
 	indicator_cvt.clear();
 
 	igzstream infile (file_sample.c_str(), igzstream::in);
 
-	if (!infile) {cout<<"error! fail to open sample file: "<<file_sample<<endl; return false;}
+	if (!infile) {
+	  cout<<"error! fail to open sample file: "<<file_sample<<endl;
+	  return false;
+	}
 
 	string line;
 	char *ch_ptr;
 
-
 	string id;
 	double p,d;
 
@@ -1888,8 +2087,6 @@ bool ReadFile_sample(const string &file_sample, vector<vector<int> > &indicator_
 	size_t num_p_in_file=0;
 	size_t num_cvt_in_file=0;
 
-//	size_t p_max=*max_element(p_column.begin(), p_column.end());
-
 	map<size_t, size_t> mapP2c;
 	for (size_t i=0; i<p_column.size(); i++) {
 		mapP2c[p_column[i]]=i;
@@ -1897,7 +2094,7 @@ bool ReadFile_sample(const string &file_sample, vector<vector<int> > &indicator_
 		ind_pheno_row.push_back(0);
 	}
 
-	// read header line1
+	// Read header line1.
 	if(!safeGetline(infile, line).eof()) {
 		ch_ptr=strtok((char *)line.c_str(), " \t");
 		if(strcmp(ch_ptr, "ID_1")!=0) {return false;}
@@ -1916,7 +2113,8 @@ bool ReadFile_sample(const string &file_sample, vector<vector<int> > &indicator_
 	vector<map<uint32_t, size_t> > cvt_factor_levels;
 
 	char col_type[num_cols];
-	// read header line2
+	
+	// Read header line2.
 	if(!safeGetline(infile, line).eof()) {
 		ch_ptr=strtok ((char *)line.c_str(), " \t");
 		if(strcmp(ch_ptr, "0")!=0) {return false;}
@@ -1927,13 +2125,17 @@ bool ReadFile_sample(const string &file_sample, vector<vector<int> > &indicator_
 		size_t it=0;
 		ch_ptr=strtok (NULL, " \t");
 		if(ch_ptr!=NULL)
-			while(ch_ptr!=NULL){
-				col_type[it++]=ch_ptr[0];
-				if(ch_ptr[0]=='D') {cvt_factor_levels.push_back(map<uint32_t, size_t>());num_cvt_in_file++;}
-				if(ch_ptr[0]=='C') {num_cvt_in_file++;}
-				if((ch_ptr[0]=='P')||(ch_ptr[0]=='B')) {num_p_in_file++;}
-				ch_ptr=strtok(NULL, " \t");
-			}
+		  while(ch_ptr!=NULL){
+		    col_type[it++]=ch_ptr[0];
+		    if(ch_ptr[0]=='D') {
+		      cvt_factor_levels.push_back(map<uint32_t,size_t>());
+		      num_cvt_in_file++;
+		    }
+		    if(ch_ptr[0]=='C') {num_cvt_in_file++;}
+		    if((ch_ptr[0]=='P')||(ch_ptr[0]=='B')) {
+		      num_p_in_file++;}
+		    ch_ptr=strtok(NULL, " \t");
+		  }
 
 	}
 
@@ -1943,7 +2145,6 @@ bool ReadFile_sample(const string &file_sample, vector<vector<int> > &indicator_
 
 		for(int it=0;it<3;it++){ch_ptr=strtok(NULL, " \t");}
 
-
 		size_t i=0;
 		size_t p_i=0;
 		size_t fac_cvt_i=0;
@@ -1952,42 +2153,62 @@ bool ReadFile_sample(const string &file_sample, vector<vector<int> > &indicator_
 
 			if((col_type[i]=='P')||(col_type[i]=='B'))
 			{
-				if (mapP2c.count(p_i+1)!=0) {
-					if (strcmp(ch_ptr, "NA")==0) {ind_pheno_row[mapP2c[p_i+1]]=0; pheno_row[mapP2c[p_i+1]]=-9;}
-					else {p=atof(ch_ptr); ind_pheno_row[mapP2c[p_i+1]]=1; pheno_row[mapP2c[p_i+1]]=p;}
-				}
-				p_i++;
+			  if (mapP2c.count(p_i+1)!=0) {
+			    if (strcmp(ch_ptr, "NA")==0) {
+			      ind_pheno_row[mapP2c[p_i+1]]=0;
+			      pheno_row[mapP2c[p_i+1]]=-9;
+			    }
+			    else {
+			      p=atof(ch_ptr);
+			      ind_pheno_row[mapP2c[p_i+1]]=1;
+			      pheno_row[mapP2c[p_i+1]]=p;
+			    }
+			  }
+			  p_i++;
 			}
 			if(col_type[i]=='D')
 			{
-				// NOTE THIS DOES NOT CHECK TO BE SURE LEVEL IS INTEGRAL i.e for atoi error
-				if (strcmp(ch_ptr, "NA")!=0) {uint32_t level=atoi(ch_ptr); if(cvt_factor_levels[fac_cvt_i].count(level) == 0) {cvt_factor_levels[fac_cvt_i][level]=cvt_factor_levels[fac_cvt_i].size();}}
-				fac_cvt_i++;
+			  
+			  // NOTE THIS DOES NOT CHECK TO BE SURE LEVEL
+			  // IS INTEGRAL i.e for atoi error.
+			  if (strcmp(ch_ptr, "NA")!=0) {
+			    uint32_t level=atoi(ch_ptr);
+			    if (cvt_factor_levels[fac_cvt_i].count(level)==0) {
+			      cvt_factor_levels[fac_cvt_i][level]=
+				cvt_factor_levels[fac_cvt_i].size();
+			    }
+			  }
+			  fac_cvt_i++;
 			}
 
 			ch_ptr=strtok (NULL, " \t");
 			i++;
 		}
 
-
 		indicator_pheno.push_back(ind_pheno_row);
 		pheno.push_back(pheno_row);
 
 	}
-	// close and reopen the file
+	
+	// Close and reopen the file.
  	infile.close();
  	infile.clear();
 
-	if(num_cvt_in_file>0)
-	{
+	if(num_cvt_in_file>0) {
 		igzstream infile2 (file_sample.c_str(), igzstream::in);
 
-		if (!infile2) {cout<<"error! fail to open sample file: "<<file_sample<<endl; return false;}
-		// skip header
+		if (!infile2) {
+		  cout<<"error! fail to open sample file: "<<
+		    file_sample<<endl;
+		  return false;
+		}
+		
+		// Skip header.
 		safeGetline(infile2, line);
 		safeGetline(infile2, line);
 
-		// pull in the covariates now we now the number of factor levels
+		// Pull in the covariates now we now the number of
+		// factor levels.
 		while (!safeGetline(infile2, line).eof()) {
 
 			vector<double> v_d; flag_na=0;
@@ -1995,38 +2216,51 @@ bool ReadFile_sample(const string &file_sample, vector<vector<int> > &indicator_
 
 			for(int it=0;it<3;it++){ch_ptr=strtok(NULL, " \t");}
 
-
 			size_t i=0;
 			size_t fac_cvt_i=0;
 			size_t num_fac_levels;
 			while (i<num_cols) {
-
-				if(col_type[i]=='C')
-				{
-					if (strcmp(ch_ptr, "NA")==0) {flag_na=1; d=-9;}
-					else {d=atof(ch_ptr);}
-
-					v_d.push_back(d);
+			  
+			  if(col_type[i]=='C') {
+			    if (strcmp(ch_ptr, "NA")==0) {flag_na=1; d=-9;}
+			    else {d=atof(ch_ptr);}
+			    
+			    v_d.push_back(d);
+			  }
+			  
+			  if(col_type[i]=='D') {
+			    
+			    // NOTE THIS DOES NOT CHECK TO BE SURE
+			    // LEVEL IS INTEGRAL i.e for atoi error.
+			    num_fac_levels=cvt_factor_levels[fac_cvt_i].size();
+			    if(num_fac_levels>1) {
+			      if (strcmp(ch_ptr, "NA")==0) {
+				flag_na=1;
+				for(size_t it=0;it<num_fac_levels-1; it++) {
+				  v_d.push_back(-9);
 				}
-
-
-				if(col_type[i]=='D')
-				{
-					// NOTE THIS DOES NOT CHECK TO BE SURE LEVEL IS INTEGRAL i.e for atoi error
-					num_fac_levels=cvt_factor_levels[fac_cvt_i].size();
-					if(num_fac_levels>1)
-					{
-						if (strcmp(ch_ptr, "NA")==0) {flag_na=1; for(size_t it=0;it<num_fac_levels-1; it++) {v_d.push_back(-9);}}
-						else {uint32_t level=atoi(ch_ptr); for(size_t it=0;it<num_fac_levels-1;it++) {cvt_factor_levels[fac_cvt_i][level]==it+1 ? v_d.push_back(1.0) : v_d.push_back(0.0); }}
-					}
-					fac_cvt_i++;
+			      }
+			      else {
+				uint32_t level=atoi(ch_ptr);
+				for(size_t it=0;it<num_fac_levels-1;it++) {
+				  cvt_factor_levels[fac_cvt_i][level]==it+1 ?
+				    v_d.push_back(1.0) :
+				    v_d.push_back(0.0);
 				}
-
-				ch_ptr=strtok (NULL, " \t");
-				i++;
+			      }
+			    }
+			    fac_cvt_i++;
+			  }
+			  
+			  ch_ptr=strtok (NULL, " \t");
+			  i++;
 			}
 
-			if (flag_na==0) {indicator_cvt.push_back(1);} else {indicator_cvt.push_back(0);}
+			if (flag_na==0) {
+			  indicator_cvt.push_back(1);
+			} else {
+			  indicator_cvt.push_back(0);
+			}
 			cvt.push_back(v_d);
 
 
@@ -2035,11 +2269,20 @@ bool ReadFile_sample(const string &file_sample, vector<vector<int> > &indicator_
 		if (indicator_cvt.empty()) {n_cvt=0;}
 		else {
 			flag_na=0;
-			for (vector<int>::size_type i=0; i<indicator_cvt.size(); ++i) {
+			for (vector<int>::size_type i=0;
+			     i<indicator_cvt.size();
+			     ++i) {
 				if (indicator_cvt[i]==0) {continue;}
 
-				if (flag_na==0) {flag_na=1; n_cvt=cvt[i].size();}
-					if (flag_na!=0 && n_cvt!=cvt[i].size()) {cout<<"error! number of covariates in row "<<i<<" do not match other rows."<<endl; return false;}
+				if (flag_na==0) {
+				  flag_na=1;
+				  n_cvt=cvt[i].size();
+				}
+				if (flag_na!=0 && n_cvt!=cvt[i].size()) {
+				  cout<<"error! number of covariates in row "<<
+				    i<<" do not match other rows."<<endl;
+				  return false;
+				}
 			}
 		}
 
@@ -2049,19 +2292,22 @@ bool ReadFile_sample(const string &file_sample, vector<vector<int> > &indicator_
  	return true;
 }
 
-
-
-// WJA Added
-//Read bgen file, the first time
-#include <cstdint>
-#include <assert.h>
-bool ReadFile_bgen(const string &file_bgen, const set<string> &setSnps, const gsl_matrix *W, vector<int> &indicator_idv, vector<int> &indicator_snp, vector<SNPINFO> &snpInfo, const double &maf_level, const double &miss_level, const double &hwe_level, const double &r2_level, size_t &ns_test)
-{
+// WJA Added.
+// Read bgen file, the first time.
+bool ReadFile_bgen(const string &file_bgen, const set<string> &setSnps,
+		   const gsl_matrix *W, vector<int> &indicator_idv,
+		   vector<int> &indicator_snp, vector<SNPINFO> &snpInfo,
+		   const double &maf_level, const double &miss_level,
+		   const double &hwe_level, const double &r2_level,
+		   size_t &ns_test) {
 
 	indicator_snp.clear();
 
 	ifstream infile (file_bgen.c_str(), ios::binary);
-	if (!infile) {cout<<"error reading bgen file:"<<file_bgen<<endl; return false;}
+	if (!infile) {
+	  cout<<"error reading bgen file:"<<file_bgen<<endl;
+	  return false;
+	}
 
 	gsl_vector *genotype=gsl_vector_alloc (W->size1);
 	gsl_vector *genotype_miss=gsl_vector_alloc (W->size1);
@@ -2075,8 +2321,8 @@ bool ReadFile_bgen(const string &file_bgen, const set<string> &setSnps, const gs
 	int sig;
 	LUDecomp (WtW, pmt, &sig);
 	LUInvert (WtW, pmt, WtWi);
-
-	// read in header
+	
+	// Read in header.
 	uint32_t bgen_snp_block_offset;
 	uint32_t bgen_header_length;
 	uint32_t bgen_nsamples;
@@ -2108,7 +2354,6 @@ bool ReadFile_bgen(const string &file_bgen, const set<string> &setSnps, const gs
 	string rs;
 	long int b_pos;
 	string chr;
-//	double cM;
 	string major;
 	string minor;
 	string id;
@@ -2116,17 +2361,19 @@ bool ReadFile_bgen(const string &file_bgen, const set<string> &setSnps, const gs
 	double v_x, v_w;
 	int c_idv=0;
 
-
 	double maf, geno, geno_old;
 	size_t n_miss;
 	size_t n_0, n_1, n_2;
 	int flag_poly;
 
-	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;
+	double bgen_geno_prob_BB, bgen_geno_prob_non_miss;
 
+	// Total number of samples in phenotype file.
+	size_t ni_total=indicator_idv.size();
 
-	size_t ni_total=indicator_idv.size();   // total number of samples in phenotype file
-	size_t ni_test=0;   // number of samples to use in test
+	// Number of samples to use in test.
+	size_t ni_test=0;   
 
 	uint32_t bgen_N;
 	uint16_t bgen_LS;
@@ -2141,13 +2388,9 @@ bool ReadFile_bgen(const string &file_bgen, const set<string> &setSnps, const gs
 	size_t unzipped_data_size;
 
 	for (size_t i=0; i<ni_total; ++i) {
-
-		ni_test+=indicator_idv[i];
+	  ni_test+=indicator_idv[i];
 	}
 
-
-
-//	ns_total=1;
 	for (size_t t=0; t<ns_total; ++t) {
 
 		id.clear();
@@ -2181,8 +2424,7 @@ bool ReadFile_bgen(const string &file_bgen, const set<string> &setSnps, const gs
 		bgen_B_allele.resize(bgen_LB);
 		infile.read(&bgen_B_allele[0], bgen_LB);
 
-
-		// should we switch according to MAF?
+		// Should we switch according to MAF?
 		minor=bgen_B_allele;
 		major=bgen_A_allele;
 		b_pos=static_cast<long int>(bgen_SNP_pos);
@@ -2196,16 +2438,15 @@ bool ReadFile_bgen(const string &file_bgen, const set<string> &setSnps, const gs
 			snpInfo.push_back(sInfo);
 			indicator_snp.push_back(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)
 		{
 			infile.read(reinterpret_cast<char*>(&bgen_P),4);
@@ -2213,36 +2454,48 @@ bool ReadFile_bgen(const string &file_bgen, const set<string> &setSnps, const gs
 
 			unzipped_data_size=6*bgen_N;
 
-			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));
+			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));
 			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);
 
 		}
 
-
 		maf=0; n_miss=0; flag_poly=0; geno_old=-9;
 		n_0=0; n_1=0; n_2=0;
 		c_idv=0;
 		gsl_vector_set_zero (genotype_miss);
 		for (size_t i=0; i<bgen_N; ++i) {
+		  
 			// CHECK this set correctly!
 			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;
-			bgen_geno_prob_non_miss=bgen_geno_prob_AA+bgen_geno_prob_AB+bgen_geno_prob_BB;
-
-			//CHECK 0.1 OK
-			if (bgen_geno_prob_non_miss<0.9) {gsl_vector_set (genotype_miss, c_idv, 1); n_miss++; c_idv++; 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;
+			bgen_geno_prob_non_miss=
+			  bgen_geno_prob_AA+bgen_geno_prob_AB+bgen_geno_prob_BB;
+
+			//CHECK 0.1 OK.
+			if (bgen_geno_prob_non_miss<0.9) {
+			  gsl_vector_set (genotype_miss, c_idv, 1);
+			  n_miss++;
+			  c_idv++;
+			  continue;
+			}
 
 			bgen_geno_prob_AA/=bgen_geno_prob_non_miss;
 			bgen_geno_prob_AB/=bgen_geno_prob_non_miss;
@@ -2255,7 +2508,7 @@ bool ReadFile_bgen(const string &file_bgen, const set<string> &setSnps, const gs
 
 			gsl_vector_set (genotype, c_idv, geno);
 
-			// CHECK WHAT THIS DOES
+			// CHECK WHAT THIS DOES.
 			if (flag_poly==0) {geno_old=geno; flag_poly=2;}
 			if (flag_poly==2 && geno!=geno_old) {flag_poly=1;}
 
@@ -2266,23 +2519,39 @@ bool ReadFile_bgen(const string &file_bgen, const set<string> &setSnps, const gs
 
 		maf/=2.0*static_cast<double>(ni_test-n_miss);
 
-		SNPINFO sInfo={chr, rs, -9, b_pos, minor, major, n_miss, (double)n_miss/(double)ni_test, maf};
+		SNPINFO sInfo={chr, rs, -9, b_pos, minor, major, n_miss,
+			       (double)n_miss/(double)ni_test, maf};
 		snpInfo.push_back(sInfo);
 
-		if ( (double)n_miss/(double)ni_test > miss_level) {indicator_snp.push_back(0); continue;}
+		if ( (double)n_miss/(double)ni_test > miss_level) {
+		  indicator_snp.push_back(0);
+		  continue;
+		}
 
-		if ( (maf<maf_level || maf> (1.0-maf_level)) && maf_level!=-1 ) {indicator_snp.push_back(0); continue;}
+		if ((maf<maf_level || maf> (1.0-maf_level)) && maf_level!=-1) {
+		  indicator_snp.push_back(0);
+		  continue;
+		}
 
-		if (flag_poly!=1) {indicator_snp.push_back(0); continue;}
+		if (flag_poly!=1) {
+		  indicator_snp.push_back(0);
+		  continue;
+		}
 
 		if (hwe_level!=0 && maf_level!=-1) {
-			if (CalcHWE(n_0, n_2, n_1)<hwe_level) {indicator_snp.push_back(0); continue;}
+		  if (CalcHWE(n_0, n_2, n_1)<hwe_level) {
+		    indicator_snp.push_back(0);
+		    continue;
+		  }
 		}
 
-		//filter SNP if it is correlated with W
-		//unless W has only one column, of 1s
+		// Filter SNP if it is correlated with W
+		// unless W has only one column, of 1s.
 		for (size_t i=0; i<genotype->size; ++i) {
-			if (gsl_vector_get (genotype_miss, i)==1) {geno=maf*2.0; gsl_vector_set (genotype, i, geno);}
+			if (gsl_vector_get (genotype_miss, i)==1) {
+			  geno=maf*2.0;
+			  gsl_vector_set (genotype, i, geno);
+			}
 		}
 
 		gsl_blas_dgemv (CblasTrans, 1.0, W, genotype, 0.0, Wtx);
@@ -2290,30 +2559,29 @@ bool ReadFile_bgen(const string &file_bgen, const set<string> &setSnps, const gs
 		gsl_blas_ddot (genotype, genotype, &v_x);
 		gsl_blas_ddot (Wtx, WtWiWtx, &v_w);
 
-		if (W->size2!=1 && v_w/v_x >= r2_level) {indicator_snp.push_back(0); continue;}
+		if (W->size2!=1 && v_w/v_x >= r2_level) {
+		  indicator_snp.push_back(0); continue;}
 
 		indicator_snp.push_back(1);
 		ns_test++;
 
 	}
 
-
-
-
 	return true;
-
 }
 
-
-//read oxford genotype file and calculate kinship matrix
-bool bgenKin (const string &file_oxford, vector<int> &indicator_snp, const int k_mode, const int display_pace, gsl_matrix *matrix_kin)
-{
+// Read oxford genotype file and calculate kinship matrix.
+bool bgenKin (const string &file_oxford, vector<int> &indicator_snp,
+	      const int k_mode, const int display_pace,
+	      gsl_matrix *matrix_kin) {
 	string file_bgen=file_oxford;
 	ifstream infile (file_bgen.c_str(), ios::binary);
-	if (!infile) {cout<<"error reading bgen file:"<<file_bgen<<endl; return false;}
-
+	if (!infile) {
+	  cout<<"error reading bgen file:"<<file_bgen<<endl;
+	  return false;
+	}
 
-	// read in header
+	// Read in header.
 	uint32_t bgen_snp_block_offset;
 	uint32_t bgen_header_length;
 	uint32_t bgen_nsamples;
@@ -2331,11 +2599,11 @@ bool bgenKin (const string &file_oxford, vector<int> &indicator_snp, const int k
 	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;
+	double bgen_geno_prob_BB, bgen_geno_prob_non_miss;
 
 	uint32_t bgen_N;
 	uint16_t bgen_LS;
@@ -2353,7 +2621,6 @@ bool bgenKin (const string &file_oxford, vector<int> &indicator_snp, const int k
 	string chr;
 	double genotype;
 
-
 	size_t n_miss;
 	double d, geno_mean, geno_var;
 
@@ -2364,7 +2631,9 @@ bool bgenKin (const string &file_oxford, vector<int> &indicator_snp, const int k
 	size_t ns_test=0;
 	for (size_t t=0; t<indicator_snp.size(); ++t) {
 
-		if (t%display_pace==0 || t==(indicator_snp.size()-1)) {ProgressBar ("Reading SNPs  ", t, indicator_snp.size()-1);}
+		if (t%display_pace==0 || t==(indicator_snp.size()-1)) {
+		  ProgressBar ("Reading SNPs  ", t, indicator_snp.size()-1);
+		}
 
 		id.clear();
 		rs.clear();
@@ -2396,74 +2665,76 @@ bool bgenKin (const string &file_oxford, vector<int> &indicator_snp, const int k
 		infile.read(reinterpret_cast<char*>(&bgen_LB),4);
 		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)
 		{
-
-
-			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);
-
-			int result=uncompress(reinterpret_cast<Bytef*>(unzipped_data), reinterpret_cast<uLongf*>(&unzipped_data_size), reinterpret_cast<Bytef*>(zipped_data), static_cast<uLong> (bgen_P));
+		  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);
+		  
+		  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);
 		}
 
-
-
 		geno_mean=0.0; n_miss=0; geno_var=0.0;
 		gsl_vector_set_all(geno_miss, 0);
 
 		for (size_t i=0; i<bgen_N; ++i) {
 
-
-				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(geno_miss, i, 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;
-
-					genotype=2.0*bgen_geno_prob_BB+bgen_geno_prob_AB;
-
-					gsl_vector_set(geno, i, genotype);
-					gsl_vector_set(geno_miss, i, 1.0);
-					geno_mean+=genotype;
-					geno_var+=genotype*genotype;
-				}
+		  
+		  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(geno_miss, i, 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;
+		    
+		    genotype=2.0*bgen_geno_prob_BB+bgen_geno_prob_AB;
+		    
+		    gsl_vector_set(geno, i, genotype);
+		    gsl_vector_set(geno_miss, i, 1.0);
+		    geno_mean+=genotype;
+		    geno_var+=genotype*genotype;
+		  }
 
 		}
 
@@ -2472,18 +2743,24 @@ bool bgenKin (const string &file_oxford, vector<int> &indicator_snp, const int k
 		geno_var+=geno_mean*geno_mean*(double)n_miss;
 		geno_var/=(double)ni_total;
 		geno_var-=geno_mean*geno_mean;
-//		geno_var=geno_mean*(1-geno_mean*0.5);
 
 		for (size_t i=0; i<ni_total; ++i) {
-			if (gsl_vector_get (geno_miss, i)==0) {gsl_vector_set(geno, i, geno_mean);}
+		  if (gsl_vector_get (geno_miss, i)==0) {
+		    gsl_vector_set(geno, i, geno_mean);
+		  }
 		}
 
 		gsl_vector_add_constant (geno, -1.0*geno_mean);
 
 		if (geno_var!=0) {
-			if (k_mode==1) {gsl_blas_dsyr (CblasUpper, 1.0, geno, matrix_kin);}
-			else if (k_mode==2) {gsl_blas_dsyr (CblasUpper, 1.0/geno_var, geno, matrix_kin);}
-			else {cout<<"Unknown kinship mode."<<endl;}
+		  if (k_mode==1) {
+		    gsl_blas_dsyr(CblasUpper,1.0,geno,matrix_kin);
+		  } else if (k_mode==2) {
+		    gsl_blas_dsyr(CblasUpper,1.0/geno_var,geno,matrix_kin);
+		  }
+		  else {
+		    cout<<"Unknown kinship mode."<<endl;
+		  }
 		}
 
 		ns_test++;
@@ -2508,42 +2785,23 @@ bool bgenKin (const string &file_oxford, vector<int> &indicator_snp, const int k
 	return true;
 }
 
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-//read header to determine which column contains which item
+// Read header to determine which column contains which item.
 bool ReadHeader_io (const string &line, HEADER &header)
 {
-  string rs_ptr[]={"rs","RS","snp","SNP","snps","SNPS","snpid","SNPID","rsid","RSID","MarkerName"};
+  string rs_ptr[]={"rs","RS","snp","SNP","snps","SNPS","snpid","SNPID",
+		   "rsid","RSID","MarkerName"};
   set<string> rs_set(rs_ptr, rs_ptr+11);
   string chr_ptr[]={"chr","CHR"};
   set<string> chr_set(chr_ptr, chr_ptr+2);
-  string pos_ptr[]={"ps","PS","pos","POS","base_position","BASE_POSITION", "bp", "BP"};
+  string pos_ptr[]={"ps","PS","pos","POS","base_position",
+		    "BASE_POSITION", "bp", "BP"};
   set<string> pos_set(pos_ptr, pos_ptr+8);
   string cm_ptr[]={"cm","CM"};
   set<string> cm_set(cm_ptr, cm_ptr+2);
   string a1_ptr[]={"a1","A1","allele1","ALLELE1","Allele1","INC_ALLELE"};
   set<string> a1_set(a1_ptr, a1_ptr+5);
-  string a0_ptr[]={"a0","A0","allele0","ALLELE0","Allele0","a2","A2","allele2","ALLELE2","Allele2","DEC_ALLELE"};
+  string a0_ptr[]={"a0","A0","allele0","ALLELE0","Allele0","a2","A2",
+		   "allele2","ALLELE2","Allele2","DEC_ALLELE"};
   set<string> a0_set(a0_ptr, a0_ptr+10);
 
   string z_ptr[]={"z","Z","z_score","Z_SCORE","zscore","ZSCORE"};
@@ -2568,7 +2826,10 @@ bool ReadHeader_io (const string &line, HEADER &header)
   string ncontrol_ptr[]={"ncontrol","NCONTROL","n_control","N_CONTROL"};
   set<string> ncontrol_set(ncontrol_ptr, ncontrol_ptr+4);
 
-  string af_ptr[]={"af","AF","maf","MAF","f","F","allele_freq","ALLELE_FREQ","allele_frequency","ALLELE_FREQUENCY","Freq.Allele1.HapMapCEU","FreqAllele1HapMapCEU", "Freq1.Hapmap"};
+  string af_ptr[]={"af","AF","maf","MAF","f","F","allele_freq",
+		   "ALLELE_FREQ","allele_frequency","ALLELE_FREQUENCY",
+		   "Freq.Allele1.HapMapCEU","FreqAllele1HapMapCEU",
+		   "Freq1.Hapmap"};
   set<string> af_set(af_ptr, af_ptr+13);
   string var_ptr[]={"var","VAR"};
   set<string> var_set(var_ptr, var_ptr+2);
@@ -2578,7 +2839,13 @@ bool ReadHeader_io (const string &line, HEADER &header)
   string cor_ptr[]={"cor","COR","r","R"};
   set<string> cor_set(cor_ptr, cor_ptr+4);
 
-  header.rs_col=0; header.chr_col=0; header.pos_col=0;  header.cm_col=0; header.a1_col=0; header.a0_col=0; header.z_col=0; header.beta_col=0; header.sebeta_col=0; header.chisq_col=0; header.p_col=0; header.n_col=0; header.nmis_col=0; header.nobs_col=0; header.ncase_col=0; header.ncontrol_col=0; header.af_col=0; header.var_col=0; header.ws_col=0; header.cor_col=0; header.coln=0;
+  header.rs_col=0; header.chr_col=0; header.pos_col=0;
+  header.cm_col=0; header.a1_col=0; header.a0_col=0; header.z_col=0;
+  header.beta_col=0; header.sebeta_col=0; header.chisq_col=0;
+  header.p_col=0; header.n_col=0; header.nmis_col=0;
+  header.nobs_col=0; header.ncase_col=0; header.ncontrol_col=0;
+  header.af_col=0; header.var_col=0; header.ws_col=0;
+  header.cor_col=0; header.coln=0;
 
   char *ch_ptr;
   string type;
@@ -2588,50 +2855,147 @@ bool ReadHeader_io (const string &line, HEADER &header)
   while (ch_ptr!=NULL) {
     type=ch_ptr;
     if (rs_set.count(type)!=0) {
-      if (header.rs_col==0) {header.rs_col=header.coln+1;} else {cout<<"error! more than two rs columns in the file."<<endl; n_error++;}
+      if (header.rs_col==0) {
+	header.rs_col=header.coln+1;
+      } else {
+	cout<<"error! more than two rs columns in the file."<<endl;
+	n_error++;
+      }
     } else if (chr_set.count(type)!=0) {
-      if (header.chr_col==0) {header.chr_col=header.coln+1;} else {cout<<"error! more than two chr columns in the file."<<endl; n_error++;}
+      if (header.chr_col==0) {
+	header.chr_col=header.coln+1;
+      } else {
+	cout<<"error! more than two chr columns in the file."<<endl;
+	n_error++;
+      }
     } else if (pos_set.count(type)!=0) {
-      if (header.pos_col==0) {header.pos_col=header.coln+1;} else {cout<<"error! more than two pos columns in the file."<<endl; n_error++;}
+      if (header.pos_col==0) {
+	header.pos_col=header.coln+1;
+      } else {
+	cout<<"error! more than two pos columns in the file."<<endl;
+	n_error++;
+      }
     } else if (cm_set.count(type)!=0) {
-      if (header.cm_col==0) {header.cm_col=header.coln+1;} else {cout<<"error! more than two cm columns in the file."<<endl; n_error++;}
+      if (header.cm_col==0) {
+	header.cm_col=header.coln+1;
+      } else {
+	cout<<"error! more than two cm columns in the file."<<endl;
+	n_error++;
+      }
     } else if (a1_set.count(type)!=0) {
-      if (header.a1_col==0) {header.a1_col=header.coln+1;} else {cout<<"error! more than two allele1 columns in the file."<<endl; n_error++;}
+      if (header.a1_col==0) {
+	header.a1_col=header.coln+1;
+      } else {
+	cout<<"error! more than two allele1 columns in the file."<<endl;
+	n_error++;
+      }
     } else if (a0_set.count(type)!=0) {
-      if (header.a0_col==0) {header.a0_col=header.coln+1;} else {cout<<"error! more than two allele0 columns in the file."<<endl; n_error++;}
+      if (header.a0_col==0) {
+	header.a0_col=header.coln+1;
+      } else {
+	cout<<"error! more than two allele0 columns in the file."<<endl;
+	n_error++;
+      }
     } else if (z_set.count(type)!=0) {
-      if (header.z_col==0) {header.z_col=header.coln+1;} else {cout<<"error! more than two z columns in the file."<<endl; n_error++;}
+      if (header.z_col==0) {
+	header.z_col=header.coln+1;
+      } else {
+	cout<<"error! more than two z columns in the file."<<endl;
+	n_error++;
+      }
     } else if (beta_set.count(type)!=0) {
-      if (header.beta_col==0) {header.beta_col=header.coln+1;} else {cout<<"error! more than two beta columns in the file."<<endl; n_error++;}
+      if (header.beta_col==0) {
+	header.beta_col=header.coln+1;
+      } else {
+	cout<<"error! more than two beta columns in the file."<<endl;
+	n_error++;
+      }
     } else if (sebeta_set.count(type)!=0) {
-      if (header.sebeta_col==0) {header.sebeta_col=header.coln+1;} else {cout<<"error! more than two se_beta columns in the file."<<endl; n_error++;}
+      if (header.sebeta_col==0) {
+	header.sebeta_col=header.coln+1;
+      } else {
+	cout<<"error! more than two se_beta columns in the file."<<endl;
+	n_error++;
+      }
     } else if (chisq_set.count(type)!=0) {
-      if (header.chisq_col==0) {header.chisq_col=header.coln+1;} else {cout<<"error! more than two z columns in the file."<<endl; n_error++;}
+      if (header.chisq_col==0) {
+	header.chisq_col=header.coln+1;
+      } else {
+	cout<<"error! more than two z columns in the file."<<endl;
+	n_error++;
+      }
     } else if (p_set.count(type)!=0) {
-      if (header.p_col==0) {header.p_col=header.coln+1;} else {cout<<"error! more than two p columns in the file."<<endl; n_error++;}
+      if (header.p_col==0) {
+	header.p_col=header.coln+1;
+      } else {
+	cout<<"error! more than two p columns in the file."<<endl;
+	n_error++;
+      }
     } else if (n_set.count(type)!=0) {
-      if (header.n_col==0) {header.n_col=header.coln+1;} else {cout<<"error! more than two n_total columns in the file."<<endl; n_error++;}
+      if (header.n_col==0) {
+	header.n_col=header.coln+1;
+      } else {
+	cout<<"error! more than two n_total columns in the file."<<endl;
+	n_
+	  error++;}
     } else if (nmis_set.count(type)!=0) {
-      if (header.nmis_col==0) {header.nmis_col=header.coln+1;} else {cout<<"error! more than two n_mis columns in the file."<<endl; n_error++;}
+      if (header.nmis_col==0) {header.nmis_col=header.coln+1;} else {
+	cout<<"error! more than two n_mis columns in the file."<<endl;
+	n_error++;
+      }
     } else if (nobs_set.count(type)!=0) {
-      if (header.nobs_col==0) {header.nobs_col=header.coln+1;} else {cout<<"error! more than two n_obs columns in the file."<<endl; n_error++;}
+      if (header.nobs_col==0) {
+	header.nobs_col=header.coln+1;
+      } else {
+	cout<<"error! more than two n_obs columns in the file."<<endl;
+	n_error++;
+      }
     } else if (ncase_set.count(type)!=0) {
-      if (header.ncase_col==0) {header.ncase_col=header.coln+1;} else {cout<<"error! more than two n_case columns in the file."<<endl; n_error++;}
+      if (header.ncase_col==0) {
+	header.ncase_col=header.coln+1;
+      } else {
+	cout<<"error! more than two n_case columns in the file."<<endl;
+	n_error++;
+      }
     } else if (ncontrol_set.count(type)!=0) {
-      if (header.ncontrol_col==0) {header.ncontrol_col=header.coln+1;} else {cout<<"error! more than two n_control columns in the file."<<endl; n_error++;}
+      if (header.ncontrol_col==0) {
+	header.ncontrol_col=header.coln+1;
+      } else {
+	cout<<"error! more than two n_control columns in the file."<<endl;
+	n_error++;
+      }
     } else if (ws_set.count(type)!=0) {
-      if (header.ws_col==0) {header.ws_col=header.coln+1;} else {cout<<"error! more than two window_size columns in the file."<<endl; n_error++;}
+      if (header.ws_col==0) {
+	header.ws_col=header.coln+1;
+      } else {
+	cout<<"error! more than two window_size columns in the file."<<endl;
+	n_error++;
+      }
     } else if (af_set.count(type)!=0) {
-      if (header.af_col==0) {header.af_col=header.coln+1;} else {cout<<"error! more than two af columns in the file."<<endl; n_error++;}
+      if (header.af_col==0) {
+	header.af_col=header.coln+1;
+      } else {
+	cout<<"error! more than two af columns in the file."<<endl;
+	n_error++;
+      }
     } else if (cor_set.count(type)!=0) {
-      if (header.cor_col==0) {header.cor_col=header.coln+1;} else {cout<<"error! more than two cor columns in the file."<<endl; n_error++;}
+      if (header.cor_col==0) {
+	header.cor_col=header.coln+1;
+      } else {
+	cout<<"error! more than two cor columns in the file."<<endl;
+	n_error++;
+      }
     } else {
       string str = ch_ptr;
       string cat = str.substr(str.size()-2, 2);
-      // continuous
+      
       if(cat == "_c" || cat =="_C"){
+
+        // continuous
 	header.catc_col.insert(header.coln+1);
-      } else { //discrete
+      } else {
+
+	// discrete
 	header.catd_col.insert(header.coln+1);
       }
     }
@@ -2640,7 +3004,10 @@ bool ReadHeader_io (const string &line, HEADER &header)
     header.coln++;
   }
 
-  if (header.cor_col!=0 && header.cor_col!=header.coln) {cout<<"error! the cor column should be the last column."<<endl; n_error++;}
+  if (header.cor_col!=0 && header.cor_col!=header.coln) {
+    cout<<"error! the cor column should be the last column."<<endl;
+    n_error++;
+  }
 
   if (header.rs_col==0) {
     if (header.chr_col!=0 && header.pos_col!=0) {
@@ -2650,34 +3017,38 @@ bool ReadHeader_io (const string &line, HEADER &header)
     }
   }
 
-  if (n_error==0) {return true;} else {return false;}
+  if (n_error==0) {
+    return true;
+  } else {
+    return false;
+  }
 }
 
-
-
-
-//read category file, record mapRS2in
-//the category file does not contain a null category
-//so if a snp has 0 for all categories, then it is not included in the analysis
-bool ReadFile_cat (const string &file_cat, map<string, size_t> &mapRS2cat, size_t &n_vc)
-{
+// Read category file, record mapRS2 in the category file does not
+// contain a null category so if a snp has 0 for all categories, then
+// it is not included in the analysis.
+bool ReadFile_cat (const string &file_cat, map<string, size_t> &mapRS2cat,
+		   size_t &n_vc) {
   mapRS2cat.clear();
 
   igzstream infile (file_cat.c_str(), igzstream::in);
-  if (!infile) {cout<<"error! fail to open category file: "<<file_cat<<endl; return false;}
+  if (!infile) {
+    cout<<"error! fail to open category file: "<<file_cat<<endl;
+    return false;
+  }
 
   string line;
   char *ch_ptr;
 
   string rs, chr, a1, a0, pos, cm;
-  size_t i_cat;// ns_vc=0;
+  size_t i_cat;
 
-  //read header
+  // Read header.
   HEADER header;
   !safeGetline(infile, line).eof();
   ReadHeader_io (line, header);
 
-  //use the header to count the number of categories
+  // Use the header to count the number of categories.
   n_vc=header.coln;
   if (header.rs_col!=0) {n_vc--;}
   if (header.chr_col!=0) {n_vc--;}
@@ -2686,7 +3057,7 @@ bool ReadFile_cat (const string &file_cat, map<string, size_t> &mapRS2cat, size_
   if (header.a1_col!=0) {n_vc--;}
   if (header.a0_col!=0) {n_vc--;}
 
-  //read the following lines to record mapRS2cat
+  // Read the following lines to record mapRS2cat.
   while (!safeGetline(infile, line).eof()) {
     ch_ptr=strtok ((char *)line.c_str(), " , \t");
 
@@ -2717,27 +3088,23 @@ bool ReadFile_cat (const string &file_cat, map<string, size_t> &mapRS2cat, size_
 
       ch_ptr=strtok (NULL, " , \t");
     }
-
-    //if (mapRS2cat.count(rs)==0) {mapRS2cat[rs]=n_vc+1; ns_vc++;}
   }
 
-  //if (ns_vc>0) {n_vc++;}
-
   infile.clear();
   infile.close();
 
   return true;
 }
 
-
-
-
-bool ReadFile_mcat (const string &file_mcat, map<string, size_t> &mapRS2cat, size_t &n_vc)
-{
+bool ReadFile_mcat (const string &file_mcat, map<string, size_t> &mapRS2cat,
+		    size_t &n_vc) {
   mapRS2cat.clear();
 
   igzstream infile (file_mcat.c_str(), igzstream::in);
-  if (!infile) {cout<<"error! fail to open mcategory file: "<<file_mcat<<endl; return false;}
+  if (!infile) {
+    cout<<"error! fail to open mcategory file: "<<file_mcat<<endl;
+    return false;
+  }
 
   string file_name;
   map<string, size_t> mapRS2cat_tmp;
@@ -2754,125 +3121,22 @@ bool ReadFile_mcat (const string &file_mcat, map<string, size_t> &mapRS2cat, siz
   return true;
 }
 
-
-
-
-
-/*
-//read the continuous category file, record mapR2catc
-bool ReadFile_catc (const string &file_cat, map<string, vector<double> > &mapRS2catc, size_t &n_cat)
-{
-  mapRS2catc.clear();
-
-  igzstream infile (file_cat.c_str(), igzstream::in);
-  if (!infile) {cout<<"error! fail to open category file: "<<file_cat<<endl; return false;}
-
-  string line;
-  char *ch_ptr;
-
-  string rs, chr, a1, a0, pos, cm;
-  size_t i_cat;// ns_vc=0;
-
-  //read header
-  HEADER header;
-  !safeGetline(infile, line).eof();
-  ReadHeader_io (line, header);
-
-  //use the header to count the number of categories
-  n_cat=header.coln;
-  if (header.rs_col!=0) {n_cat--;}
-  if (header.chr_col!=0) {n_cat--;}
-  if (header.pos_col!=0) {n_cat--;}
-  if (header.cm_col!=0) {n_cat--;}
-  if (header.a1_col!=0) {n_cat--;}
-  if (header.a0_col!=0) {n_cat--;}
-
-  //set up continous category
-  vector<double> catc;
-  for (size_t i=0; i<n_cat; i++) {
-    catc.push_back(0);
-  }
-
-  //read the following lines to record mapRS2cat
-  while (!safeGetline(infile, line).eof()) {
-    ch_ptr=strtok ((char *)line.c_str(), " , \t");
-
-    i_cat=0;
-    if (header.rs_col==0) {
-      rs=chr+":"+pos;
-    }
-
-    for (size_t i=0; i<header.coln; i++) {
-      if (header.rs_col!=0 && header.rs_col==i+1) {
-	rs=ch_ptr;
-      } else if (header.chr_col!=0 && header.chr_col==i+1) {
-	chr=ch_ptr;
-      } else if (header.pos_col!=0 && header.pos_col==i+1) {
-	pos=ch_ptr;
-      } else if (header.cm_col!=0 && header.cm_col==i+1) {
-	cm=ch_ptr;
-      } else if (header.a1_col!=0 && header.a1_col==i+1) {
-	a1=ch_ptr;
-      } else if (header.a0_col!=0 && header.a0_col==i+1) {
-	a0=ch_ptr;
-      } else {
-	catc[i_cat]=atof(ch_ptr);
-	i_cat++;
-      }
-
-      ch_ptr=strtok (NULL, " , \t");
-    }
-
-    if (mapRS2catc.count(rs)==0) {mapRS2catc[rs]=catc;}
-
-    //if (mapRS2cat.count(rs)==0) {mapRS2cat[rs]=n_vc+1; ns_vc++;}
-  }
-
-  //if (ns_vc>0) {n_vc++;}
-
-  infile.clear();
-  infile.close();
-
-  return true;
-}
-
-
-
-
-bool ReadFile_mcatc (const string &file_mcat, map<string, vector<double> > &mapRS2catc, size_t &n_cat)
-{
-  mapRS2catc.clear();
-
-  igzstream infile (file_mcat.c_str(), igzstream::in);
-  if (!infile) {cout<<"error! fail to open mcategory file: "<<file_mcat<<endl; return false;}
-
-  string file_name;
-  map<string, vector<double> > mapRS2catc_tmp;
-  size_t n_cat_tmp, t=0;
-
-  while (!safeGetline(infile, file_name).eof()) {
-    mapRS2catc_tmp.clear();
-    ReadFile_catc (file_name, mapRS2catc_tmp, n_cat_tmp);
-    mapRS2catc.insert(mapRS2catc_tmp.begin(), mapRS2catc_tmp.end());
-    if (t==0) {n_cat=n_cat_tmp;}
-    if (n_cat!=n_cat_tmp) {cout<<"number of category differs in different mcatc files."<<endl;;}
-
-    t++;
-  }
-
-  return true;
-}
-*/
-
-
-
-
-//read bimbam mean genotype file and calculate kinship matrix; this time, the kinship matrix is not centered, and can contain multiple K matrix
-bool BimbamKin (const string &file_geno, const int display_pace, const vector<int> &indicator_idv, const vector<int> &indicator_snp, const map<string, double> &mapRS2weight, const map<string, size_t> &mapRS2cat, const vector<SNPINFO> &snpInfo, const gsl_matrix *W, gsl_matrix *matrix_kin, gsl_vector *vector_ns)
-{
+// Read bimbam mean genotype file and calculate kinship matrix; this
+// time, the kinship matrix is not centered, and can contain multiple
+// K matrix.
+bool BimbamKin (const string &file_geno, const int display_pace,
+		const vector<int> &indicator_idv,
+		const vector<int> &indicator_snp,
+		const map<string, double> &mapRS2weight,
+		const map<string, size_t> &mapRS2cat,
+		const vector<SNPINFO> &snpInfo,
+		const gsl_matrix *W, gsl_matrix *matrix_kin,
+		gsl_vector *vector_ns) {
 	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 false;}
+	if (!infile) {
+	  cout<<"error reading genotype file:"<<file_geno<<endl;
+	  return false;
+	}
 
 	string line;
 	char *ch_ptr;
@@ -2902,7 +3166,7 @@ bool BimbamKin (const string &file_geno, const int display_pace, const vector<in
 	  ns_vec.push_back(0);
 	}
 
-	//create a large matrix
+	// Create a large matrix.
 	size_t msize=10000;
 	gsl_matrix *Xlarge=gsl_matrix_alloc (ni_test, msize*n_vc);
 	gsl_matrix_set_zero(Xlarge);
@@ -2910,14 +3174,16 @@ bool BimbamKin (const string &file_geno, const int display_pace, const vector<in
 	size_t ns_test=0;
 	for (size_t t=0; t<indicator_snp.size(); ++t) {
 		!safeGetline(infile, line).eof();
-		if (t%display_pace==0 || t==(indicator_snp.size()-1)) {ProgressBar ("Reading SNPs  ", t, indicator_snp.size()-1);}
+		if (t%display_pace==0 || t==(indicator_snp.size()-1)) {
+		  ProgressBar ("Reading SNPs  ", t, indicator_snp.size()-1);
+		}
 		if (indicator_snp[t]==0) {continue;}
 
 		ch_ptr=strtok ((char *)line.c_str(), " , \t");
 		ch_ptr=strtok (NULL, " , \t");
 		ch_ptr=strtok (NULL, " , \t");
 
-		rs=snpInfo[t].rs_number;//this line is new
+		rs=snpInfo[t].rs_number; // This line is new.
 
 		geno_mean=0.0; n_miss=0; geno_var=0.0;
 		gsl_vector_set_all(geno_miss, 0);
@@ -2926,13 +3192,15 @@ bool BimbamKin (const string &file_geno, const int display_pace, const vector<in
 		for (size_t i=0; i<indicator_idv.size(); ++i) {
 		  if (indicator_idv[i]==0) {continue;}
 			ch_ptr=strtok (NULL, " , \t");
-			if (strcmp(ch_ptr, "NA")==0) {gsl_vector_set(geno_miss, i, 0); n_miss++;}
+			if (strcmp(ch_ptr, "NA")==0) {
+			  gsl_vector_set(geno_miss, i, 0); n_miss++;
+			}
 			else {
-				d=atof(ch_ptr);
-				gsl_vector_set (geno, j, d);
-				gsl_vector_set (geno_miss, j, 1);
-				geno_mean+=d;
-				geno_var+=d*d;
+			  d=atof(ch_ptr);
+			  gsl_vector_set (geno, j, d);
+			  gsl_vector_set (geno_miss, j, 1);
+			  geno_mean+=d;
+			  geno_var+=d*d;
 			}
 			j++;
 		}
@@ -2941,10 +3209,11 @@ bool BimbamKin (const string &file_geno, const int display_pace, const vector<in
 		geno_var+=geno_mean*geno_mean*(double)n_miss;
 		geno_var/=(double)ni_test;
 		geno_var-=geno_mean*geno_mean;
-//		geno_var=geno_mean*(1-geno_mean*0.5);
 
 		for (size_t i=0; i<ni_test; ++i) {
-			if (gsl_vector_get (geno_miss, i)==0) {gsl_vector_set(geno, i, geno_mean);}
+		  if (gsl_vector_get (geno_miss, i)==0) {
+		    gsl_vector_set(geno, i, geno_mean);
+		  }
 		}
 
 		gsl_vector_add_constant (geno, -1.0*geno_mean);
@@ -2955,48 +3224,43 @@ bool BimbamKin (const string &file_geno, const int display_pace, const vector<in
 		gsl_blas_ddot (geno, geno, &geno_var);
 		geno_var/=(double)ni_test;
 
-		if (geno_var!=0 && (mapRS2weight.size()==0 || mapRS2weight.count(rs)!=0) ) {
+		if (geno_var!=0 && (mapRS2weight.size()==0 ||
+				    mapRS2weight.count(rs)!=0)) {
 		  if (mapRS2weight.size()==0) {
 		    d=1.0/geno_var;
 		  } else {
 		    d=mapRS2weight.at(rs)/geno_var;
 		  }
 
-		  /*
-		  if (n_vc==1 || mapRS2cat.size()==0 ) {
-		    gsl_blas_dsyr (CblasUpper, d, geno, matrix_kin);
-		    ns_vec[0]++;
-		  } else if (mapRS2cat.count(rs)!=0) {
-		      i_vc=mapRS2cat.at(rs);
-		      ns_vec[i_vc]++;
-		      gsl_matrix_view kin_sub=gsl_matrix_submatrix(matrix_kin, 0, ni_test*i_vc, ni_test, ni_test);
-		      gsl_blas_dsyr (CblasUpper, d, geno, &kin_sub.matrix);
-		      //eigenlib_dsyr (1.0, geno, matrix_kin);
-		  }
-		  */
-
 		  gsl_vector_scale (geno, sqrt(d));
 		  if (n_vc==1 || mapRS2cat.size()==0 ) {
-		    gsl_vector_view Xlarge_col=gsl_matrix_column (Xlarge, ns_vec[0]%msize);
+		    gsl_vector_view Xlarge_col=
+		      gsl_matrix_column(Xlarge,ns_vec[0]%msize);
 		    gsl_vector_memcpy (&Xlarge_col.vector, geno);
 		    ns_vec[0]++;
 
 		    if (ns_vec[0]%msize==0) {
-		      eigenlib_dgemm ("N", "T", 1.0, Xlarge, Xlarge, 1.0, matrix_kin);
+		      eigenlib_dgemm("N","T",1.0,Xlarge,Xlarge,1.0,matrix_kin);
 		      gsl_matrix_set_zero(Xlarge);
 		    }
 		  } else if (mapRS2cat.count(rs)!=0) {
 		    i_vc=mapRS2cat.at(rs);
 
-		    gsl_vector_view Xlarge_col=gsl_matrix_column (Xlarge, msize*i_vc+ns_vec[i_vc]%msize);
+		    gsl_vector_view Xlarge_col=
+		      gsl_matrix_column(Xlarge,msize*i_vc+ns_vec[i_vc]%msize);
 		    gsl_vector_memcpy (&Xlarge_col.vector, geno);
 
 		    ns_vec[i_vc]++;
 
 		    if (ns_vec[i_vc]%msize==0) {
-		      gsl_matrix_view X_sub=gsl_matrix_submatrix(Xlarge, 0, msize*i_vc, ni_test, msize);
-		      gsl_matrix_view kin_sub=gsl_matrix_submatrix(matrix_kin, 0, ni_test*i_vc, ni_test, ni_test);
-		      eigenlib_dgemm ("N", "T", 1.0, &X_sub.matrix, &X_sub.matrix, 1.0, &kin_sub.matrix);
+		      gsl_matrix_view X_sub=
+			gsl_matrix_submatrix(Xlarge,0,msize*i_vc,
+					     ni_test,msize);
+		      gsl_matrix_view kin_sub=
+			gsl_matrix_submatrix(matrix_kin,0,ni_test*i_vc,
+					     ni_test,ni_test);
+		      eigenlib_dgemm ("N", "T", 1.0, &X_sub.matrix,
+				      &X_sub.matrix, 1.0, &kin_sub.matrix);
 
 		      gsl_matrix_set_zero(&X_sub.matrix);
 		    }
@@ -3009,9 +3273,13 @@ bool BimbamKin (const string &file_geno, const int display_pace, const vector<in
 
 	for (size_t i_vc=0; i_vc<n_vc; i_vc++) {
 	  if (ns_vec[i_vc]%msize!=0) {
-	    gsl_matrix_view X_sub=gsl_matrix_submatrix(Xlarge, 0, msize*i_vc, ni_test, msize);
-	    gsl_matrix_view kin_sub=gsl_matrix_submatrix(matrix_kin, 0, ni_test*i_vc, ni_test, ni_test);
-	    eigenlib_dgemm ("N", "T", 1.0, &X_sub.matrix, &X_sub.matrix, 1.0, &kin_sub.matrix);
+	    gsl_matrix_view X_sub=
+	      gsl_matrix_submatrix(Xlarge, 0, msize*i_vc, ni_test, msize);
+	    gsl_matrix_view kin_sub=
+	      gsl_matrix_submatrix(matrix_kin, 0, ni_test*i_vc, ni_test,
+				   ni_test);
+	    eigenlib_dgemm ("N", "T", 1.0, &X_sub.matrix, &X_sub.matrix,
+			    1.0, &kin_sub.matrix);
 	  }
 	}
 
@@ -3047,16 +3315,19 @@ bool BimbamKin (const string &file_geno, const int display_pace, const vector<in
 	return true;
 }
 
-
-
-
-
-
-
-bool PlinkKin (const string &file_bed, const int display_pace, const vector<int> &indicator_idv, const vector<int> &indicator_snp, const map<string, double> &mapRS2weight, const map<string, size_t> &mapRS2cat, const vector<SNPINFO> &snpInfo, const gsl_matrix *W, gsl_matrix *matrix_kin, gsl_vector *vector_ns)
-{
+bool PlinkKin (const string &file_bed, const int display_pace,
+	       const vector<int> &indicator_idv,
+	       const vector<int> &indicator_snp,
+	       const map<string, double> &mapRS2weight,
+	       const map<string, size_t> &mapRS2cat,
+	       const vector<SNPINFO> &snpInfo,
+	       const gsl_matrix *W, gsl_matrix *matrix_kin,
+	       gsl_vector *vector_ns) {
 	ifstream infile (file_bed.c_str(), ios::binary);
-	if (!infile) {cout<<"error reading bed file:"<<file_bed<<endl; return false;}
+	if (!infile) {
+	  cout<<"error reading bed file:"<<file_bed<<endl;
+	  return false;
+	}
 
 	char ch[1];
 	bitset<8> b;
@@ -3089,58 +3360,68 @@ bool PlinkKin (const string &file_bed, const int display_pace, const vector<int>
 	  ns_vec.push_back(0);
 	}
 
-	//create a large matrix
+	// Create a large matrix.
 	size_t msize=10000;
 	gsl_matrix *Xlarge=gsl_matrix_alloc (ni_test, msize*n_vc);
 	gsl_matrix_set_zero(Xlarge);
 
-	//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 (size_t t=0; t<indicator_snp.size(); ++t) {
-		if (t%display_pace==0 || t==(indicator_snp.size()-1)) {ProgressBar ("Reading SNPs  ", t, indicator_snp.size()-1);}
+		if (t%display_pace==0 || t==(indicator_snp.size()-1)) {
+		  ProgressBar ("Reading SNPs  ", t, indicator_snp.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);
 
-		rs=snpInfo[t].rs_number;//this line is new
+		rs=snpInfo[t].rs_number; // This line is new.
 
-		//read genotypes
+		// Read genotypes.
 		geno_mean=0.0;	n_miss=0; ci_total=0; geno_var=0.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==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(geno, ci_test, 2.0); geno_mean+=2.0; geno_var+=4.0; }
-					else {gsl_vector_set(geno, ci_test, 1.0); geno_mean+=1.0; geno_var+=1.0;}
-				}
-				else {
-					if (b[2*j+1]==1) {gsl_vector_set(geno, ci_test, 0.0); }
-					else {gsl_vector_set(geno, ci_test, -9.0); n_miss++; }
-				}
-
-				ci_test++;
-				ci_total++;
-			}
+		  infile.read(ch,1);
+		  b=ch[0];
+
+		  // Minor allele homozygous: 2.0; major: 0.0;
+		  for (size_t j=0; j<4; ++j) {
+		    if ((i==(n_bit-1)) && ci_total==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(geno, ci_test, 2.0);
+			geno_mean+=2.0; geno_var+=4.0;
+		      }
+		      else {
+			gsl_vector_set(geno, ci_test, 1.0);
+			geno_mean+=1.0;
+			geno_var+=1.0;
+		      }
+		    }
+		    else {
+		      if (b[2*j+1]==1) {gsl_vector_set(geno, ci_test, 0.0); }
+		      else {gsl_vector_set(geno, ci_test, -9.0); n_miss++; }
+		    }
+		    
+		    ci_test++;
+		    ci_total++;
+		  }
 		}
 
-
 		geno_mean/=(double)(ni_test-n_miss);
 		geno_var+=geno_mean*geno_mean*(double)n_miss;
 		geno_var/=(double)ni_test;
 		geno_var-=geno_mean*geno_mean;
-//		geno_var=geno_mean*(1-geno_mean*0.5);
 
 		for (size_t i=0; i<ni_test; ++i) {
 			d=gsl_vector_get(geno,i);
@@ -3155,47 +3436,43 @@ bool PlinkKin (const string &file_bed, const int display_pace, const vector<int>
 		gsl_blas_ddot (geno, geno, &geno_var);
 		geno_var/=(double)ni_test;
 
-		if (geno_var!=0 && (mapRS2weight.size()==0 || mapRS2weight.count(rs)!=0) ) {
+		if (geno_var!=0 && (mapRS2weight.size()==0 ||
+				    mapRS2weight.count(rs)!=0)) {
 		  if (mapRS2weight.size()==0) {
 		    d=1.0/geno_var;
 		  } else {
 		    d=mapRS2weight.at(rs)/geno_var;
 		  }
 
-		  /*
-		  if (n_vc==1 || mapRS2cat.size()==0 ) {
-		    gsl_blas_dsyr (CblasUpper, d, geno, matrix_kin);
-		    ns_vec[0]++;
-		  } else if (mapRS2cat.count(rs)!=0) {
-		    i_vc=mapRS2cat.at(rs);
-		    ns_vec[i_vc]++;
-		    gsl_matrix_view kin_sub=gsl_matrix_submatrix(matrix_kin, 0, ni_test*i_vc, ni_test, ni_test);
-		    gsl_blas_dsyr (CblasUpper, d, geno, &kin_sub.matrix);
-		  }
-		  */
-
 		  gsl_vector_scale (geno, sqrt(d));
 		  if (n_vc==1 || mapRS2cat.size()==0 ) {
-		    gsl_vector_view Xlarge_col=gsl_matrix_column (Xlarge, ns_vec[0]%msize);
+		    gsl_vector_view Xlarge_col=
+		      gsl_matrix_column (Xlarge, ns_vec[0]%msize);
 		    gsl_vector_memcpy (&Xlarge_col.vector, geno);
 		    ns_vec[0]++;
 
 		    if (ns_vec[0]%msize==0) {
-		      eigenlib_dgemm ("N", "T", 1.0, Xlarge, Xlarge, 1.0, matrix_kin);
+		      eigenlib_dgemm("N","T",1.0,Xlarge,Xlarge,1.0,matrix_kin);
 		      gsl_matrix_set_zero(Xlarge);
 		    }
 		  } else if (mapRS2cat.count(rs)!=0) {
 		    i_vc=mapRS2cat.at(rs);
 
-		    gsl_vector_view Xlarge_col=gsl_matrix_column (Xlarge, msize*i_vc+ns_vec[i_vc]%msize);
+		    gsl_vector_view Xlarge_col=
+		      gsl_matrix_column(Xlarge,msize*i_vc+ns_vec[i_vc]%msize);
 		    gsl_vector_memcpy (&Xlarge_col.vector, geno);
 
 		    ns_vec[i_vc]++;
 
 		    if (ns_vec[i_vc]%msize==0) {
-		      gsl_matrix_view X_sub=gsl_matrix_submatrix(Xlarge, 0, msize*i_vc, ni_test, msize);
-		      gsl_matrix_view kin_sub=gsl_matrix_submatrix(matrix_kin, 0, ni_test*i_vc, ni_test, ni_test);
-		      eigenlib_dgemm ("N", "T", 1.0, &X_sub.matrix, &X_sub.matrix, 1.0, &kin_sub.matrix);
+		      gsl_matrix_view X_sub=
+			gsl_matrix_submatrix(Xlarge,0,msize*i_vc,ni_test,
+					     msize);
+		      gsl_matrix_view kin_sub=
+			gsl_matrix_submatrix(matrix_kin, 0, ni_test*i_vc,
+					     ni_test, ni_test);
+		      eigenlib_dgemm ("N", "T", 1.0, &X_sub.matrix,
+				      &X_sub.matrix, 1.0, &kin_sub.matrix);
 
 		      gsl_matrix_set_zero(&X_sub.matrix);
 		    }
@@ -3208,9 +3485,13 @@ bool PlinkKin (const string &file_bed, const int display_pace, const vector<int>
 
 	for (size_t i_vc=0; i_vc<n_vc; i_vc++) {
 	  if (ns_vec[i_vc]%msize!=0) {
-	    gsl_matrix_view X_sub=gsl_matrix_submatrix(Xlarge, 0, msize*i_vc, ni_test, msize);
-	    gsl_matrix_view kin_sub=gsl_matrix_submatrix(matrix_kin, 0, ni_test*i_vc, ni_test, ni_test);
-	    eigenlib_dgemm ("N", "T", 1.0, &X_sub.matrix, &X_sub.matrix, 1.0, &kin_sub.matrix);
+	    gsl_matrix_view X_sub=
+	      gsl_matrix_submatrix(Xlarge, 0, msize*i_vc, ni_test, msize);
+	    gsl_matrix_view kin_sub=
+	      gsl_matrix_submatrix(matrix_kin, 0, ni_test*i_vc,
+				   ni_test, ni_test);
+	    eigenlib_dgemm ("N", "T", 1.0, &X_sub.matrix, &X_sub.matrix,
+			    1.0, &kin_sub.matrix);
 	  }
 	}
 
@@ -3245,16 +3526,23 @@ bool PlinkKin (const string &file_bed, const int display_pace, const vector<int>
 	return true;
 }
 
-
-
-bool MFILEKin (const size_t mfile_mode, const string &file_mfile, const int display_pace, const vector<int> &indicator_idv, const vector<vector<int> > &mindicator_snp, const map<string, double> &mapRS2weight, const map<string, size_t> &mapRS2cat, const vector<vector<SNPINFO> > &msnpInfo, const gsl_matrix *W, gsl_matrix *matrix_kin, gsl_vector *vector_ns)
-{
+bool MFILEKin (const size_t mfile_mode, const string &file_mfile,
+	       const int display_pace, const vector<int> &indicator_idv,
+	       const vector<vector<int> > &mindicator_snp,
+	       const map<string, double> &mapRS2weight,
+	       const map<string, size_t> &mapRS2cat,
+	       const vector<vector<SNPINFO> > &msnpInfo,
+	       const gsl_matrix *W, gsl_matrix *matrix_kin,
+	       gsl_vector *vector_ns) {
   size_t n_vc=vector_ns->size, ni_test=matrix_kin->size1;
   gsl_matrix_set_zero(matrix_kin);
   gsl_vector_set_zero(vector_ns);
 
   igzstream infile (file_mfile.c_str(), igzstream::in);
-  if (!infile) {cout<<"error! fail to open mfile file: "<<file_mfile<<endl; return false;}
+  if (!infile) {
+    cout<<"error! fail to open mfile file: "<<file_mfile<<endl;
+    return false;
+  }
 
   string file_name;
 
@@ -3273,11 +3561,11 @@ bool MFILEKin (const size_t mfile_mode, const string &file_mfile, const int disp
     } else {
       BimbamKin (file_name, display_pace, indicator_idv, mindicator_snp[l], mapRS2weight, mapRS2cat, msnpInfo[l], W, kin_tmp, ns_tmp);
     }
-
-    //add ns
+    
+    // Add ns.
     gsl_vector_add(vector_ns, ns_tmp);
 
-    //add kin
+    // Add kin.
     for (size_t t=0; t<n_vc; t++) {
       for (size_t i=0; i<ni_test; ++i) {
 	for (size_t j=0; j<=i; ++j) {
@@ -3291,11 +3579,12 @@ bool MFILEKin (const size_t mfile_mode, const string &file_mfile, const int disp
     l++;
   }
 
-  //renormalize kin
+  // Renormalize kin.
   for (size_t t=0; t<n_vc; t++) {
     for (size_t i=0; i<ni_test; ++i) {
       for (size_t j=0; j<=i; ++j) {
-	d=gsl_matrix_get (matrix_kin, j, i+ni_test*t)/gsl_vector_get(vector_ns, t);
+	d=gsl_matrix_get (matrix_kin, j, i+ni_test*t)/
+	  gsl_vector_get(vector_ns, t);
 
 	gsl_matrix_set (matrix_kin, i, j+ni_test*t, d);
 	gsl_matrix_set (matrix_kin, j, i+ni_test*t, d);
@@ -3315,15 +3604,16 @@ bool MFILEKin (const size_t mfile_mode, const string &file_mfile, const int disp
 }
 
 
-
-
-//read var file, store mapRS2wsnp
-bool ReadFile_wsnp (const string &file_wsnp, map<string, double> &mapRS2weight)
-{
+// Read var file, store mapRS2wsnp.
+bool ReadFile_wsnp (const string &file_wsnp,
+		    map<string, double> &mapRS2weight) {
   mapRS2weight.clear();
 
   igzstream infile (file_wsnp.c_str(), igzstream::in);
-  if (!infile) {cout<<"error! fail to open snp weight file: "<<file_wsnp<<endl; return false;}
+  if (!infile) {
+    cout<<"error! fail to open snp weight file: "<<file_wsnp<<endl;
+    return false;
+  }
 
   char *ch_ptr;
   string line, rs;
@@ -3340,12 +3630,15 @@ bool ReadFile_wsnp (const string &file_wsnp, map<string, double> &mapRS2weight)
   return true;
 }
 
-bool ReadFile_wsnp (const string &file_wcat, const size_t n_vc, map<string, vector<double> > &mapRS2wvector)
-{
+bool ReadFile_wsnp (const string &file_wcat, const size_t n_vc,
+		    map<string, vector<double> > &mapRS2wvector) {
   mapRS2wvector.clear();
 
   igzstream infile (file_wcat.c_str(), igzstream::in);
-  if (!infile) {cout<<"error! fail to open snp weight file: "<<file_wcat<<endl; return false;}
+  if (!infile) {
+    cout<<"error! fail to open snp weight file: "<<file_wcat<<endl;
+    return false;
+  }
 
   char *ch_ptr;
   vector<double> weight;
@@ -3354,10 +3647,8 @@ bool ReadFile_wsnp (const string &file_wcat, const size_t n_vc, map<string, vect
   }
 
   string line, rs, chr, a1, a0, pos, cm;
-  //double af=0, var_x=0;
-  //size_t n_total=0, n_mis=0, n_obs=0, n_case=0, n_control=0;
-
-  //read header
+  
+  // Read header.
   HEADER header;
   !safeGetline(infile, line).eof();
   ReadHeader_io (line, header);
@@ -3366,7 +3657,6 @@ bool ReadFile_wsnp (const string &file_wcat, const size_t n_vc, map<string, vect
     if (isBlankLine(line)) {continue;}
     ch_ptr=strtok ((char *)line.c_str(), " , \t");
 
-    //n_total=0; n_mis=0; n_obs=0; n_case=0; n_control=0; n_case=0; af=0; var_x=0;
     size_t t=0;
     for (size_t i=0; i<header.coln; i++) {
       if (header.rs_col!=0 && header.rs_col==i+1) {rs=ch_ptr;}
@@ -3375,22 +3665,23 @@ bool ReadFile_wsnp (const string &file_wcat, const size_t n_vc, map<string, vect
       else if (header.cm_col!=0 && header.cm_col==i+1) {cm=ch_ptr; }
       else if (header.a1_col!=0 && header.a1_col==i+1) {a1=ch_ptr; }
       else if (header.a0_col!=0 && header.a0_col==i+1) {a0=ch_ptr; }
-      //else if (header.n_col!=0 && header.n_col==i+1) {n_total=atoi(ch_ptr); }
-      //else if (header.nmis_col!=0 && header.nmis_col==i+1) {n_mis=atoi(ch_ptr); }
-      //else if (header.nobs_col!=0 && header.nobs_col==i+1) {n_obs=atoi(ch_ptr); }
-      //else if (header.ncase_col!=0 && header.ncase_col==i+1) {n_case=atoi(ch_ptr); }
-      //else if (header.ncontrol_col!=0 && header.ncontrol_col==i+1) {n_control=atoi(ch_ptr); }
-      //else if (header.af_col!=0 && header.af_col==i+1) {af=atof(ch_ptr); }
-      //else if (header.var_col!=0 && header.var_col==i+1) {var_x=atof(ch_ptr); }
       else {
 	weight[t]=atof(ch_ptr); t++;
-	if (t>n_vc) {cout<<"error! Number of columns in the wcat file does not match that of cat file."; return false;}
+	if (t>n_vc) {
+	  cout<<"error! Number of columns in the wcat file does not "<<
+	    "match that of cat file.";
+	  return false;
+	}
       }
 
       ch_ptr=strtok (NULL, " , \t");
     }
 
-    if (t!=n_vc) {cout<<"error! Number of columns in the wcat file does not match that of cat file."; return false;}
+    if (t!=n_vc) {
+      cout<<"error! Number of columns in the wcat file does not "<<
+	"match that of cat file.";
+      return false;
+    }
 
     if (header.rs_col==0) {
       rs=chr+":"+pos;
@@ -3402,25 +3693,28 @@ bool ReadFile_wsnp (const string &file_wcat, const size_t n_vc, map<string, vect
   return true;
 }
 
-
-
-
-
-
-
-
-//read the beta file, save snp z scores in to z2_score, and save category into indicator_snp based on mapRS2var and set, and indicator_snp record the category number (from 1 to n_vc), and provide var if maf/var is not provided in the beta file
-//notice that indicator_snp contains ns_test snps, instead of ns_total snps
-//read the beta file for the second time, compute q, and Vq based on block jacknife
-//use the mapRS2var to select snps (and to ), calculate q
-//do a block-wise jacknife, and compute Vq
-void ReadFile_beta (const string &file_beta, const map<string, size_t> &mapRS2cat, const map<string, double> &mapRS2wA, vector<size_t> &vec_cat, vector<size_t> &vec_ni, vector<double> &vec_weight, vector<double> &vec_z2, size_t &ni_total, size_t &ns_total, size_t &ns_test)
-{
+// Read the beta file, save snp z scores in to z2_score, and save
+// category into indicator_snp based on mapRS2var and set, and
+// indicator_snp record the category number (from 1 to n_vc), and
+// provide var if maf/var is not provided in the beta file notice that
+// indicator_snp contains ns_test snps, instead of ns_total snps read
+// the beta file for the second time, compute q, and Vq based on block
+// jacknife use the mapRS2var to select snps (and to ), calculate q do
+// a block-wise jacknife, and compute Vq
+void ReadFile_beta (const string &file_beta,
+		    const map<string, size_t> &mapRS2cat,
+		    const map<string, double> &mapRS2wA,
+		    vector<size_t> &vec_cat, vector<size_t> &vec_ni,
+		    vector<double> &vec_weight, vector<double> &vec_z2,
+		    size_t &ni_total, size_t &ns_total, size_t &ns_test) {
   vec_cat.clear(); vec_ni.clear(); vec_weight.clear(); vec_z2.clear();
   ni_total=0; ns_total=0; ns_test=0;
 
   igzstream infile (file_beta.c_str(), igzstream::in);
-  if (!infile) {cout<<"error! fail to open beta file: "<<file_beta<<endl; return;}
+  if (!infile) {
+    cout<<"error! fail to open beta file: "<<file_beta<<endl;
+    return;
+  }
 
   string line;
   char *ch_ptr;
@@ -3430,27 +3724,25 @@ void ReadFile_beta (const string &file_beta, const map<string, size_t> &mapRS2ca
   double z=0, beta=0, se_beta=0, chisq=0, pvalue=0, zsquare=0, af=0, var_x=0;
   size_t n_total=0, n_mis=0, n_obs=0, n_case=0, n_control=0;
 
-  //read header
+  // Read header.
   HEADER header;
   !safeGetline(infile, line).eof();
   ReadHeader_io (line, header);
 
   if (header.n_col==0 ) {
-    if ( (header.nobs_col==0 && header.nmis_col==0) && (header.ncase_col==0 && header.ncontrol_col==0) ) {
+    if ( (header.nobs_col==0 && header.nmis_col==0) &&
+	 (header.ncase_col==0 && header.ncontrol_col==0) ) {
       cout<<"error! missing sample size in the beta file."<<endl;
     } else {
       cout<<"total sample size will be replaced by obs/mis sample size."<<endl;
     }
   }
 
-  if (header.z_col==0 && (header.beta_col==0 || header.sebeta_col==0) && header.chisq_col==0 && header.p_col==0) {
+  if (header.z_col==0 && (header.beta_col==0 || header.sebeta_col==0) &&
+      header.chisq_col==0 && header.p_col==0) {
     cout<<"error! missing z scores in the beta file."<<endl;
   }
-  /*
-  if (header.af_col==0 && header.var_col==0) {
-    cout<<"error! missing allele frequency in the beta file."<<endl;
-  }
-  */
+
   while (!safeGetline(infile, line).eof()) {
     if (isBlankLine(line)) {continue;}
     ch_ptr=strtok ((char *)line.c_str(), " , \t");
@@ -3467,7 +3759,9 @@ void ReadFile_beta (const string &file_beta, const map<string, size_t> &mapRS2ca
 
       if (header.z_col!=0 && header.z_col==i+1) {z=atof(ch_ptr);}
       if (header.beta_col!=0 && header.beta_col==i+1) {beta=atof(ch_ptr);}
-      if (header.sebeta_col!=0 && header.sebeta_col==i+1) {se_beta=atof(ch_ptr);}
+      if (header.sebeta_col!=0 && header.sebeta_col==i+1) {
+	se_beta=atof(ch_ptr);
+      }
       if (header.chisq_col!=0 && header.chisq_col==i+1) {chisq=atof(ch_ptr);}
       if (header.p_col!=0 && header.p_col==i+1) {pvalue=atof(ch_ptr);}
 
@@ -3475,8 +3769,9 @@ void ReadFile_beta (const string &file_beta, const map<string, size_t> &mapRS2ca
       if (header.nmis_col!=0 && header.nmis_col==i+1) {n_mis=atoi(ch_ptr);}
       if (header.nobs_col!=0 && header.nobs_col==i+1) {n_obs=atoi(ch_ptr);}
       if (header.ncase_col!=0 && header.ncase_col==i+1) {n_case=atoi(ch_ptr);}
-      if (header.ncontrol_col!=0 && header.ncontrol_col==i+1) {n_control=atoi(ch_ptr);}
-
+      if (header.ncontrol_col!=0 && header.ncontrol_col==i+1) {
+	n_control=atoi(ch_ptr);
+      }
       if (header.af_col!=0 && header.af_col==i+1) {af=atof(ch_ptr);}
       if (header.var_col!=0 && header.var_col==i+1) {var_x=atof(ch_ptr);}
 
@@ -3495,7 +3790,8 @@ void ReadFile_beta (const string &file_beta, const map<string, size_t> &mapRS2ca
       }
     }
 
-    //both z values and beta/se_beta have directions, while chisq/pvalue do not
+    // Both z values and beta/se_beta have directions, while
+    // chisq/pvalue do not.
     if (header.z_col!=0) {
       zsquare=z*z;
     } else if (header.beta_col!=0 && header.sebeta_col!=0) {
@@ -3507,13 +3803,14 @@ void ReadFile_beta (const string &file_beta, const map<string, size_t> &mapRS2ca
       zsquare=gsl_cdf_chisq_Qinv (pvalue, 1);
     } else {zsquare=0;}
 
-    //obtain var_x
+    // Obtain var_x.
     if (header.var_col==0 && header.af_col!=0) {
       var_x=2.0*af*(1.0-af);
     }
 
-    //if the snp is also present in cor file, then do calculations
-    if ( (mapRS2wA.size()==0 || mapRS2wA.count(rs)!=0) && (mapRS2cat.size()==0 || mapRS2cat.count(rs)!=0) && zsquare!=0) {
+    // If the SNP is also present in cor file, then do calculations.
+    if ( (mapRS2wA.size()==0 || mapRS2wA.count(rs)!=0) &&
+	 (mapRS2cat.size()==0 || mapRS2cat.count(rs)!=0) && zsquare!=0) {
       if (mapRS2cat.size()!=0) {
 	vec_cat.push_back(mapRS2cat.at(rs));
       } else {
@@ -3540,17 +3837,17 @@ void ReadFile_beta (const string &file_beta, const map<string, size_t> &mapRS2ca
   return;
 }
 
-
-
-
-
-
-void ReadFile_beta (const string &file_beta, const map<string, double> &mapRS2wA, map<string, string> &mapRS2A1, map<string, double> &mapRS2z)
-{
+void ReadFile_beta (const string &file_beta,
+		    const map<string, double> &mapRS2wA,
+		    map<string, string> &mapRS2A1,
+		    map<string, double> &mapRS2z) {
   mapRS2A1.clear(); mapRS2z.clear();
 
   igzstream infile (file_beta.c_str(), igzstream::in);
-  if (!infile) {cout<<"error! fail to open beta file: "<<file_beta<<endl; return;}
+  if (!infile) {
+    cout<<"error! fail to open beta file: "<<file_beta<<endl;
+    return;
+  }
 
   string line;
   char *ch_ptr;
@@ -3561,13 +3858,14 @@ void ReadFile_beta (const string &file_beta, const map<string, double> &mapRS2wA
   size_t n_total=0, n_mis=0, n_obs=0, n_case=0, n_control=0;
   size_t ni_total=0, ns_total=0, ns_test=0;
 
-  //read header
+  // Read header.
   HEADER header;
   !safeGetline(infile, line).eof();
   ReadHeader_io (line, header);
 
   if (header.n_col==0 ) {
-    if ( (header.nobs_col==0 && header.nmis_col==0) && (header.ncase_col==0 && header.ncontrol_col==0) ) {
+    if ((header.nobs_col==0 && header.nmis_col==0) &&
+	(header.ncase_col==0 && header.ncontrol_col==0)) {
       cout<<"error! missing sample size in the beta file."<<endl;
     } else {
       cout<<"total sample size will be replaced by obs/mis sample size."<<endl;
@@ -3577,11 +3875,7 @@ void ReadFile_beta (const string &file_beta, const map<string, double> &mapRS2wA
   if (header.z_col==0 && (header.beta_col==0 || header.sebeta_col==0)) {
     cout<<"error! missing z scores in the beta file."<<endl;
   }
-  /*
-  if (header.af_col==0 && header.var_col==0) {
-    cout<<"error! missing allele frequency in the beta file."<<endl;
-  }
-  */
+
   while (!safeGetline(infile, line).eof()) {
     if (isBlankLine(line)) {continue;}
     ch_ptr=strtok ((char *)line.c_str(), " , \t");
@@ -3598,7 +3892,9 @@ void ReadFile_beta (const string &file_beta, const map<string, double> &mapRS2wA
 
       if (header.z_col!=0 && header.z_col==i+1) {z=atof(ch_ptr);}
       if (header.beta_col!=0 && header.beta_col==i+1) {beta=atof(ch_ptr);}
-      if (header.sebeta_col!=0 && header.sebeta_col==i+1) {se_beta=atof(ch_ptr);}
+      if (header.sebeta_col!=0 && header.sebeta_col==i+1) {
+	se_beta=atof(ch_ptr);
+      }
       if (header.chisq_col!=0 && header.chisq_col==i+1) {chisq=atof(ch_ptr);}
       if (header.p_col!=0 && header.p_col==i+1) {pvalue=atof(ch_ptr);}
 
@@ -3606,7 +3902,9 @@ void ReadFile_beta (const string &file_beta, const map<string, double> &mapRS2wA
       if (header.nmis_col!=0 && header.nmis_col==i+1) {n_mis=atoi(ch_ptr);}
       if (header.nobs_col!=0 && header.nobs_col==i+1) {n_obs=atoi(ch_ptr);}
       if (header.ncase_col!=0 && header.ncase_col==i+1) {n_case=atoi(ch_ptr);}
-      if (header.ncontrol_col!=0 && header.ncontrol_col==i+1) {n_control=atoi(ch_ptr);}
+      if (header.ncontrol_col!=0 && header.ncontrol_col==i+1) {
+	n_control=atoi(ch_ptr);
+      }
 
       if (header.af_col!=0 && header.af_col==i+1) {af=atof(ch_ptr);}
       if (header.var_col!=0 && header.var_col==i+1) {var_x=atof(ch_ptr);}
@@ -3626,7 +3924,8 @@ void ReadFile_beta (const string &file_beta, const map<string, double> &mapRS2wA
       }
     }
 
-    //both z values and beta/se_beta have directions, while chisq/pvalue do not
+    // Both z values and beta/se_beta have directions, while
+    // chisq/pvalue do not.
     if (header.z_col!=0) {
       z=z;
     } else if (header.beta_col!=0 && header.sebeta_col!=0) {
@@ -3635,7 +3934,7 @@ void ReadFile_beta (const string &file_beta, const map<string, double> &mapRS2wA
       z=0;
     }
 
-    //if the snp is also present in cor file, then do calculations
+    // If the snp is also present in cor file, then do calculations.
     if ( (mapRS2wA.size()==0 || mapRS2wA.count(rs)!=0) ) {
       mapRS2z[rs]=z;
       mapRS2A1[rs]=a1;
@@ -3653,10 +3952,10 @@ void ReadFile_beta (const string &file_beta, const map<string, double> &mapRS2wA
   return;
 }
 
-
-
-void Calcq (const size_t n_block, const vector<size_t> &vec_cat, const vector<size_t> &vec_ni, const vector<double> &vec_weight, const vector<double> &vec_z2, gsl_matrix *Vq, gsl_vector *q, gsl_vector *s)
-{
+void Calcq (const size_t n_block, const vector<size_t> &vec_cat,
+	    const vector<size_t> &vec_ni, const vector<double> &vec_weight,
+	    const vector<double> &vec_z2, gsl_matrix *Vq, gsl_vector *q,
+	    gsl_vector *s) {
   gsl_matrix_set_zero (Vq);
   gsl_vector_set_zero (q);
   gsl_vector_set_zero (s);
@@ -3677,21 +3976,22 @@ void Calcq (const size_t n_block, const vector<size_t> &vec_cat, const vector<si
     mat_s.push_back(vec_s);
   }
 
-  //compute q and s
+  // Compute q and s.
   for (size_t i=0; i<vec_cat.size(); i++) {
-    //extract quantities
+    
+    // Extract quantities.
     cat=vec_cat[i];
     n_total=vec_ni[i];
     w=vec_weight[i];
     zsquare=vec_z2[i];
 
-    //compute q and s
+    // Compute q and s.
     vec_q[cat]+=(zsquare-1.0)*w/(double)n_total;
     vec_s[cat]+=w;
     n_snps[cat]++;
   }
 
-  //update q; vec_q is used again for computing Vq below
+  // Update q; vec_q is used again for computing Vq below.
   for (size_t i=0; i<q->size; i++) {
     if (vec_s[i]!=0) {
       gsl_vector_set(q, i, vec_q[i]/vec_s[i]);
@@ -3699,14 +3999,15 @@ void Calcq (const size_t n_block, const vector<size_t> &vec_cat, const vector<si
     gsl_vector_set(s, i, vec_s[i]);
   }
 
-  //compute Vq; divide SNPs in each category into evenly distributed blocks
+  // Compute Vq; divide SNPs in each category into evenly distributed
+  // blocks.
   size_t t=0, b=0, n_snp=0;
   double d, m, n;
   for (size_t l=0; l<q->size; l++) {
     n_snp=floor(n_snps[l]/n_block); t=0; b=0;
     if (n_snp==0) {continue;}
 
-    //initiate everything to zero
+    // Initiate everything to zero.
     for (size_t i=0; i<n_block; i++) {
       for (size_t j=0; j<q->size; j++) {
 	mat_q[i][j]=0;
@@ -3714,15 +4015,17 @@ void Calcq (const size_t n_block, const vector<size_t> &vec_cat, const vector<si
       }
     }
 
-    //record values
+    // Record values.
     for (size_t i=0; i<vec_cat.size(); i++) {
-      //extract quantities
+      
+      // Extract quantities.
       cat=vec_cat[i];
       n_total=vec_ni[i];
       w=vec_weight[i];
       zsquare=vec_z2[i];
 
-      //save quantities for computing Vq (which is not divided by n_total)
+      // Save quantities for computing Vq (which is not divided by
+      // n_total).
       mat_q[b][cat]+=(zsquare-1.0)*w;
       mat_s[b][cat]+=w;
 
@@ -3735,7 +4038,7 @@ void Calcq (const size_t n_block, const vector<size_t> &vec_cat, const vector<si
       }
     }
 
-    //center mat_q
+    // Center mat_q.
     for (size_t i=0; i<q->size; i++) {
       m=0; n=0;
       for (size_t k=0; k<n_block; k++) {
@@ -3755,7 +4058,7 @@ void Calcq (const size_t n_block, const vector<size_t> &vec_cat, const vector<si
       }
     }
 
-    //compute Vq for l'th row and l'th column only
+    // Compute Vq for l'th row and l'th column only.
     for (size_t i=0; i<q->size; i++) {
       d=0; n=0;
       for (size_t k=0; k<n_block; k++) {
@@ -3788,14 +4091,14 @@ void Calcq (const size_t n_block, const vector<size_t> &vec_cat, const vector<si
   return;
 }
 
-
-
-
-//read vector file
+// Read vector file.
 void ReadFile_vector (const string &file_vec, gsl_vector *vec)
 {
   igzstream infile (file_vec.c_str(), igzstream::in);
-  if (!infile) {cout<<"error! fail to open vector file: "<<file_vec<<endl; return;}
+  if (!infile) {
+    cout<<"error! fail to open vector file: "<<file_vec<<endl;
+    return;
+  }
 
   string line;
   char *ch_ptr;
@@ -3812,11 +4115,12 @@ void ReadFile_vector (const string &file_vec, gsl_vector *vec)
   return;
 }
 
-
-void ReadFile_matrix (const string &file_mat, gsl_matrix *mat)
-{
+void ReadFile_matrix (const string &file_mat, gsl_matrix *mat) {
   igzstream infile (file_mat.c_str(), igzstream::in);
-  if (!infile) {cout<<"error! fail to open matrix file: "<<file_mat<<endl; return;}
+  if (!infile) {
+    cout<<"error! fail to open matrix file: "<<file_mat<<endl;
+    return;
+  }
 
   string line;
   char *ch_ptr;
@@ -3836,10 +4140,13 @@ void ReadFile_matrix (const string &file_mat, gsl_matrix *mat)
   return;
 }
 
-void ReadFile_matrix (const string &file_mat, gsl_matrix *mat1, gsl_matrix *mat2)
-{
+void ReadFile_matrix (const string &file_mat, gsl_matrix *mat1,
+		      gsl_matrix *mat2) {
   igzstream infile (file_mat.c_str(), igzstream::in);
-  if (!infile) {cout<<"error! fail to open matrix file: "<<file_mat<<endl; return;}
+  if (!infile) {
+    cout<<"error! fail to open matrix file: "<<file_mat<<endl;
+    return;
+  }
 
   string line;
   char *ch_ptr;
@@ -3868,11 +4175,9 @@ void ReadFile_matrix (const string &file_mat, gsl_matrix *mat1, gsl_matrix *mat2
   return;
 }
 
-
-
-//read study file
-void ReadFile_study (const string &file_study, gsl_matrix *Vq_mat, gsl_vector *q_vec, gsl_vector *s_vec, size_t &ni)
-{
+// Read study file.
+void ReadFile_study (const string &file_study, gsl_matrix *Vq_mat,
+		     gsl_vector *q_vec, gsl_vector *s_vec, size_t &ni) {
   string Vqfile=file_study+".Vq.txt";
   string sfile=file_study+".size.txt";
   string qfile=file_study+".q.txt";
@@ -3895,19 +4200,16 @@ void ReadFile_study (const string &file_study, gsl_matrix *Vq_mat, gsl_vector *q
   return;
 }
 
-
-//read reference file
-void ReadFile_ref (const string &file_ref, gsl_matrix *S_mat, gsl_matrix *Svar_mat, gsl_vector *s_vec, size_t &ni)
-{
+// Read reference file.
+void ReadFile_ref (const string &file_ref, gsl_matrix *S_mat,
+		   gsl_matrix *Svar_mat, gsl_vector *s_vec, size_t &ni) {
   string sfile=file_ref+".size.txt";
   string Sfile=file_ref+".S.txt";
-  //string Vfile=file_ref+".V.txt";
 
   gsl_vector *s=gsl_vector_alloc (s_vec->size+1);
 
   ReadFile_vector(sfile, s);
   ReadFile_matrix(Sfile, S_mat, Svar_mat);
-  //ReadFile_matrix(Vfile, V_mat);
 
   double d;
   for (size_t i=0; i<s_vec->size; i++) {
@@ -3921,10 +4223,9 @@ void ReadFile_ref (const string &file_ref, gsl_matrix *S_mat, gsl_matrix *Svar_m
   return;
 }
 
-
-//read mstudy file
-void ReadFile_mstudy (const string &file_mstudy, gsl_matrix *Vq_mat, gsl_vector *q_vec, gsl_vector *s_vec, size_t &ni)
-{
+// Read mstudy file.
+void ReadFile_mstudy (const string &file_mstudy, gsl_matrix *Vq_mat,
+		      gsl_vector *q_vec, gsl_vector *s_vec, size_t &ni) {
   gsl_matrix_set_zero(Vq_mat);
   gsl_vector_set_zero(q_vec);
   gsl_vector_set_zero(s_vec);
@@ -3935,7 +4236,10 @@ void ReadFile_mstudy (const string &file_mstudy, gsl_matrix *Vq_mat, gsl_vector
   gsl_vector *s=gsl_vector_alloc (s_vec->size+1);
 
   igzstream infile (file_mstudy.c_str(), igzstream::in);
-  if (!infile) {cout<<"error! fail to open mstudy file: "<<file_mstudy<<endl; return;}
+  if (!infile) {
+    cout<<"error! fail to open mstudy file: "<<file_mstudy<<endl;
+    return;
+  }
 
   string file_name;
   double d1, d2, d;
@@ -3996,45 +4300,42 @@ void ReadFile_mstudy (const string &file_mstudy, gsl_matrix *Vq_mat, gsl_vector
   return;
 }
 
-//read reference file
-void ReadFile_mref (const string &file_mref, gsl_matrix *S_mat, gsl_matrix *Svar_mat, gsl_vector *s_vec, size_t &ni)
-{
+// Read reference file.
+void ReadFile_mref (const string &file_mref, gsl_matrix *S_mat,
+		    gsl_matrix *Svar_mat, gsl_vector *s_vec, size_t &ni) {
   gsl_matrix_set_zero(S_mat);
   gsl_matrix_set_zero(Svar_mat);
-  //  gsl_matrix_set_zero(V_mat);
   gsl_vector_set_zero(s_vec);
   ni=0;
 
-  //size_t n_vc=S_mat->size1;
   gsl_matrix *S_sub=gsl_matrix_alloc (S_mat->size1, S_mat->size2);
   gsl_matrix *Svar_sub=gsl_matrix_alloc (Svar_mat->size1, Svar_mat->size2);
-  //gsl_matrix *V_sub=gsl_matrix_alloc (V_mat->size1, V_mat->size2);
   gsl_vector *s=gsl_vector_alloc (s_vec->size+1);
 
   igzstream infile (file_mref.c_str(), igzstream::in);
-  if (!infile) {cout<<"error! fail to open mref file: "<<file_mref<<endl; return;}
+  if (!infile) {
+    cout<<"error! fail to open mref file: "<<file_mref<<endl;
+    return;
+  }
 
   string file_name;
   double d1, d2, d;
-  //size_t t_ij;
 
   while (!safeGetline(infile, file_name).eof()) {
     string sfile=file_name+".size.txt";
     string Sfile=file_name+".S.txt";
-    //string Vfile=file_name+".V.txt";
 
     ReadFile_vector(sfile, s);
     ReadFile_matrix(Sfile, S_sub, Svar_sub);
-    //ReadFile_matrix(Vfile, V_sub);
 
-    //update s_vec and ni
+    // Update s_vec and ni.
     for (size_t i=0; i<s_vec->size; i++) {
       d=gsl_vector_get (s, i)+gsl_vector_get (s_vec, i);
       gsl_vector_set (s_vec, i, d);
     }
     ni=max(ni, (size_t)gsl_vector_get (s, s_vec->size));
 
-    //update S and Svar from each file
+    // Update S and Svar from each file.
     for (size_t i=0; i<S_mat->size1; i++) {
       d1=gsl_vector_get(s, i);
       for (size_t j=0; j<S_mat->size2; j++) {
@@ -4049,30 +4350,9 @@ void ReadFile_mref (const string &file_mref, gsl_matrix *S_mat, gsl_matrix *Svar
 
     gsl_matrix_add (S_mat, S_sub);
     gsl_matrix_add (Svar_mat, Svar_sub);
-    /*
-    //update V from each file
-    for (size_t i=0; i<n_vc; i++) {
-      d1=gsl_vector_get(s, i);
-      for (size_t j=i; j<n_vc; j++) {
-	d2=gsl_vector_get(s, j);
-	t_ij=GetabIndex (i+1, j+1, n_vc-2);
-	for (size_t l=0; l<n_vc+1; l++) {
-	  if (l==n_vc) {d3=1;} else {d3=gsl_vector_get(s, l);}
-	  for (size_t m=0; m<n_vc+1; m++) {
-	    if (m==n_vc) {d4=1;} else {d4=gsl_vector_get(s, m);}
-
-	    d=gsl_matrix_get (V_sub, l, t_ij*(n_vc+1)+m)*d1*d2*d3*d4;
-	    gsl_matrix_set (V_sub, l, t_ij*(n_vc+1)+m, d);
-	  }
-	}
-      }
-    }
-
-    gsl_matrix_add (V_mat, V_sub);
-    */
   }
 
-  //final: update S and Svar
+  // Final: update S and Svar.
   for (size_t i=0; i<S_mat->size1; i++) {
     d1=gsl_vector_get(s_vec, i);
     if (d1==0) {continue;}