From f3df6447b345c6b4dee79d9996696978520344bb Mon Sep 17 00:00:00 2001
From: Peter Carbonetto
Date: Mon, 19 Jun 2017 13:49:30 -0500
Subject: Removed FORCE_FLOAT from param.cpp.

---
 src/bslmm.cpp |   30 +-
 src/io.cpp    |   26 +-
 src/param.cpp | 1124 ++++++++++++++++++++++++++++++++++-----------------------
 3 files changed, 699 insertions(+), 481 deletions(-)

diff --git a/src/bslmm.cpp b/src/bslmm.cpp
index 92762e2..563b743 100644
--- a/src/bslmm.cpp
+++ b/src/bslmm.cpp
@@ -1730,10 +1730,14 @@ void BSLMM::MCMC (const gsl_matrix *X, const gsl_vector *y) {
 					gsl_vector_view Xtznew_sub=gsl_vector_subvector(Xtz_new, 0, rank_new.size());
 					gsl_vector_view betanew_sub=gsl_vector_subvector(beta_new, 0, rank_new.size());
 
-					gsl_matrix_memcpy(&Xold_sub.matrix, &Xnew_sub.matrix);
-					gsl_matrix_memcpy(&XtXold_sub.matrix, &XtXnew_sub.matrix);
-					gsl_vector_memcpy(&Xtzold_sub.vector, &Xtznew_sub.vector);
-					gsl_vector_memcpy(&betaold_sub.vector, &betanew_sub.vector);
+					gsl_matrix_memcpy(&Xold_sub.matrix,
+							  &Xnew_sub.matrix);
+					gsl_matrix_memcpy(&XtXold_sub.matrix,
+							  &XtXnew_sub.matrix);
+					gsl_vector_memcpy(&Xtzold_sub.vector,
+							  &Xtznew_sub.vector);
+					gsl_vector_memcpy(&betaold_sub.vector,
+							  &betanew_sub.vector);
 				}
 			} else {
 			  cHyp_new=cHyp_old;
@@ -1777,12 +1781,18 @@ void BSLMM::MCMC (const gsl_matrix *X, const gsl_vector *y) {
 					}
 				}
 
-				gsl_matrix_set(Result_hyp,w_col,0,cHyp_old.h);
-				gsl_matrix_set(Result_hyp,w_col,1,cHyp_old.pve);
-				gsl_matrix_set(Result_hyp,w_col,2,cHyp_old.rho);
-				gsl_matrix_set(Result_hyp,w_col,3,cHyp_old.pge);
-				gsl_matrix_set(Result_hyp,w_col,4,cHyp_old.logp);
-				gsl_matrix_set(Result_hyp,w_col,5,cHyp_old.n_gamma);
+				gsl_matrix_set(Result_hyp,w_col,0,
+					       cHyp_old.h);
+				gsl_matrix_set(Result_hyp,w_col,1,
+					       cHyp_old.pve);
+				gsl_matrix_set(Result_hyp,w_col,2,
+					       cHyp_old.rho);
+				gsl_matrix_set(Result_hyp,w_col,3,
+					       cHyp_old.pge);
+				gsl_matrix_set(Result_hyp,w_col,4,
+					       cHyp_old.logp);
+				gsl_matrix_set(Result_hyp,w_col,5,
+					       cHyp_old.n_gamma);
 				
 				for (size_t i=0; i<cHyp_old.n_gamma; ++i) {
 					pos=mapRank2pos[rank_old[i]]+1;
diff --git a/src/io.cpp b/src/io.cpp
index e8f00ed..40f1c16 100644
--- a/src/io.cpp
+++ b/src/io.cpp
@@ -4089,30 +4089,8 @@ void ReadFile_mref (const string &file_mref, gsl_matrix *S_mat, gsl_matrix *Svar
       if (i!=j) {gsl_matrix_set(Svar_mat, j, i, d);}
     }
   }
-  /*
-  //final: update V
-  for (size_t i=0; i<n_vc; i++) {
-    d1=gsl_vector_get(s_vec, i);
-    if (d1==0) {continue;}
-    for (size_t j=i; j<n_vc; j++) {
-      d2=gsl_vector_get(s_vec, j);
-      if (d2==0) {continue;}
-      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_vec, l);}
-	  if (d3==0) {continue;}
-	  for (size_t m=0; m<n_vc+1; m++) {
-	    if (m==n_vc) {d4=1;} else {d4=gsl_vector_get(s_vec, m);}
-	    if (d4==0) {continue;}
-
-	    d=gsl_matrix_get (V_mat, l, t_ij*(n_vc+1)+m)/(d1*d2*d3*d4);
-	    gsl_matrix_set (V_mat, l, t_ij*(n_vc+1)+m, d);
-	  }
-	}
-      }
-    }
-  */
-  //free matrices
+  
+  // Free matrices.
   gsl_matrix_free(S_sub);
   gsl_matrix_free(Svar_sub);
   gsl_vector_free(s);
diff --git a/src/param.cpp b/src/param.cpp
index a7be178..a56f5eb 100644
--- a/src/param.cpp
+++ b/src/param.cpp
@@ -175,7 +175,6 @@ void PARAM::ReadFiles (void) {
 	  }
 	}
 
-
 	// WJA added.
 	// Read genotype and phenotype file for bgen format.
 	if (!file_oxford.empty()) {
@@ -281,7 +280,6 @@ void PARAM::ReadFiles (void) {
 		ns_total=indicator_snp.size();
 	}
 
-
 	// Read genotype file for multiple PLINK files.
 	if (!file_mbfile.empty()) {
 	  igzstream infile (file_mbfile.c_str(), igzstream::in);
@@ -515,37 +513,94 @@ void PARAM::CheckParam (void) {
 	  error=true;
 	}
 	if (h_max<h_min) {
-	  cout<<"error! maximum h value must be larger than the minimal value. current values = "<<h_max<<" and "<<h_min<<endl;
+	  cout<<"error! maximum h value must be larger than the minimal "<<
+	    "value. current values = "<<h_max<<" and "<<h_min<<endl;
+	  error=true;
+	}
+	if (s_max<s_min) {
+	  cout<<"error! maximum s value must be larger than the minimal "<<
+	    "value. current values = "<<s_max<<" and "<<s_min<<endl;
+	  error=true;
+	}
+	if (rho_max<rho_min) {
+	  cout<<"error! maximum rho value must be larger than the"<<
+	    "minimal value. current values = "<<rho_max<<" and "<<
+	    rho_min<<endl;
+	  error=true;
+	}
+	if (logp_max<logp_min) {
+	  cout<<"error! maximum logp value must be larger than the "<<
+	    "minimal value. current values = "<<logp_max/log(10)<<
+	    " and "<<logp_min/log(10)<<endl;
 	  error=true;
 	}
-	if (s_max<s_min) {cout<<"error! maximum s value must be larger than the minimal value. current values = "<<s_max<<" and "<<s_min<<endl; error=true;}
-	if (rho_max<rho_min) {cout<<"error! maximum rho value must be larger than the minimal value. current values = "<<rho_max<<" and "<<rho_min<<endl; error=true;}
-	if (logp_max<logp_min) {cout<<"error! maximum logp value must be larger than the minimal value. current values = "<<logp_max/log(10)<<" and "<<logp_min/log(10)<<endl; error=true;}
 
-	if (h_max>1) {cout<<"error! h values must be bewtween 0 and 1. current values = "<<h_max<<" and "<<h_min<<endl; error=true;}
-	if (rho_max>1) {cout<<"error! rho values must be between 0 and 1. current values = "<<rho_max<<" and "<<rho_min<<endl; error=true;}
-	if (logp_max>0) {cout<<"error! maximum logp value must be smaller than 0. current values = "<<logp_max/log(10)<<" and "<<logp_min/log(10)<<endl; error=true;}
-	if (l_max<l_min) {cout<<"error! maximum lambda value must be larger than the minimal value. current values = "<<l_max<<" and "<<l_min<<endl; error=true;}
+	if (h_max>1) {
+	  cout<<"error! h values must be bewtween 0 and 1. current "<<
+	    "values = "<<h_max<<" and "<<h_min<<endl;
+	  error=true;
+	}
+	if (rho_max>1) {
+	  cout<<"error! rho values must be between 0 and 1. current "<<
+	    "values = "<<rho_max<<" and "<<rho_min<<endl;
+	  error=true;
+	}
+	if (logp_max>0) {
+	  cout<<"error! maximum logp value must be smaller than 0. "<<
+	    "current values = "<<logp_max/log(10)<<" and "<<
+	    logp_min/log(10)<<endl;
+	  error=true;
+	}
+	if (l_max<l_min) {
+	  cout<<"error! maximum lambda value must be larger than the "<<
+	    "minimal value. current values = "<<l_max<<" and "<<l_min<<endl;
+	  error=true;
+	}
 
-	if (h_scale>1.0) {cout<<"error! hscale value must be between 0 and 1. current value = "<<h_scale<<endl; error=true;}
-	if (rho_scale>1.0) {cout<<"error! rscale value must be between 0 and 1. current value = "<<rho_scale<<endl; error=true;}
-	if (logp_scale>1.0) {cout<<"error! pscale value must be between 0 and 1. current value = "<<logp_scale<<endl; error=true;}
+	if (h_scale>1.0) {
+	  cout<<"error! hscale value must be between 0 and 1. "<<
+	    "current value = "<<h_scale<<endl;
+	  error=true;
+	}
+	if (rho_scale>1.0) {
+	  cout<<"error! rscale value must be between 0 and 1. "<<
+	    "current value = "<<rho_scale<<endl;
+	  error=true;
+	}
+	if (logp_scale>1.0) {
+	  cout<<"error! pscale value must be between 0 and 1. "<<
+	    "current value = "<<logp_scale<<endl;
+	  error=true;
+	}
 
-	if (rho_max==1 && rho_min==1 && a_mode==12) {cout<<"error! ridge regression does not support a rho parameter. current values = "<<rho_max<<" and "<<rho_min<<endl; error=true;}
+	if (rho_max==1 && rho_min==1 && a_mode==12) {
+	  cout<<"error! ridge regression does not support a rho "<<
+	    "parameter. current values = "<<rho_max<<" and "<<rho_min<<endl;
+	  error=true;
+	}
 
-	if (window_cm<0) {cout<<"error! windowcm values must be non-negative. current values = "<<window_cm<<endl; error=true;}
+	if (window_cm<0) {
+	  cout<<"error! windowcm values must be non-negative. "<<
+	    "current values = "<<window_cm<<endl;
+	  error=true;
+	}
 
 	if (window_cm==0 && window_bp==0 && window_ns==0) {
 	  window_bp=1000000;
 	}
 
-	//check p_column, and (no need to) sort p_column into ascending order
+	// Check p_column, and (no need to) sort p_column into
+	// ascending order.
 	if (p_column.size()==0) {
 		p_column.push_back(1);
 	} else {
 		for (size_t i=0; i<p_column.size(); i++) {
 			for (size_t j=0; j<i; j++) {
-				if (p_column[i]==p_column[j]) {cout<<"error! identical phenotype columns: "<<p_column[i]<<endl; error=true;}
+				if (p_column[i]==p_column[j]) {
+				  cout<<"error! identical phenotype "<<
+				    "columns: "<<p_column[i]<<endl;
+				  error=
+				    true;}
 			}
 		}
 	}
@@ -554,15 +609,22 @@ void PARAM::CheckParam (void) {
 
 	// Only LMM option (and one prediction option) can deal with
 	// multiple phenotypes and no gene expression files.
-	if (n_ph>1 && a_mode!=1 && a_mode!=2 && a_mode!=3 && a_mode!=4 && a_mode!=43) {
-		cout<<"error! the current analysis mode "<<a_mode<<" can not deal with multiple phenotypes."<<endl; error=true;
+	if (n_ph>1 && a_mode!=1 && a_mode!=2 && a_mode!=3 && a_mode!=4 &&
+	    a_mode!=43) {
+		cout<<"error! the current analysis mode "<<a_mode<<
+		  " can not deal with multiple phenotypes."<<endl;
+		error=true;
 	}
 	if (n_ph>1 && !file_gene.empty() ) {
-		cout<<"error! multiple phenotype analysis option not allowed with gene expression files. "<<endl; error=true;
+	  cout<<"error! multiple phenotype analysis option not "<<
+	    "allowed with gene expression files. "<<endl;
+	  error=true;
 	}
 
 	if (p_nr>1) {
-		cout<<"error! pnr value must be between 0 and 1. current value = "<<p_nr<<endl; error=true;
+	  cout<<"error! pnr value must be between 0 and 1. current value = "<<
+	    p_nr<<endl;
+	  error=true;
 	}
 
 	//check est_column
@@ -580,86 +642,162 @@ void PARAM::CheckParam (void) {
 		}
 	}
 
-	if (est_column.size()!=4) {cout<<"error! -en not followed by four numbers. current number = "<<est_column.size()<<endl; error=true;}
-	if (est_column[0]==0) {cout<<"error! -en rs column can not be zero. current number = "<<est_column.size()<<endl; error=true;}
+	if (est_column.size()!=4) {
+	  cout<<"error! -en not followed by four numbers. current number = "<<
+	    est_column.size()<<endl;
+	  error=true;
+	}
+	if (est_column[0]==0) {
+	  cout<<"error! -en rs column can not be zero. current number = "<<
+	    est_column.size()<<endl;
+	  error=true;
+	}
 
-	//check if files are compatible with each other, and if files exist
+	// Check if files are compatible with each other, and if files exist.
 	if (!file_bfile.empty()) {
 		str=file_bfile+".bim";
-		if (stat(str.c_str(),&fileInfo)==-1) {cout<<"error! fail to open .bim file: "<<str<<endl; error=true;}
+		if (stat(str.c_str(),&fileInfo)==-1) {
+		  cout<<"error! fail to open .bim file: "<<str<<endl;
+		  error=true;
+		}
 		str=file_bfile+".bed";
-		if (stat(str.c_str(),&fileInfo)==-1) {cout<<"error! fail to open .bed file: "<<str<<endl; error=true;}
+		if (stat(str.c_str(),&fileInfo)==-1) {
+		  cout<<"error! fail to open .bed file: "<<str<<endl;
+		  error=true;
+		}
 		str=file_bfile+".fam";
-		if (stat(str.c_str(),&fileInfo)==-1) {cout<<"error! fail to open .fam file: "<<str<<endl; error=true;}
+		if (stat(str.c_str(),&fileInfo)==-1) {
+		  cout<<"error! fail to open .fam file: "<<str<<endl;
+		  error=true;
+		}
 	}
 
 	if (!file_oxford.empty()) {
 		str=file_oxford+".bgen";
-		if (stat(str.c_str(),&fileInfo)==-1) {cout<<"error! fail to open .bgen file: "<<str<<endl; error=true;}
+		if (stat(str.c_str(),&fileInfo)==-1) {
+		  cout<<"error! fail to open .bgen file: "<<str<<endl;
+		  error=true;
+		}
 		str=file_oxford+".sample";
-		if (stat(str.c_str(),&fileInfo)==-1) {cout<<"error! fail to open .sample file: "<<str<<endl; error=true;}
+		if (stat(str.c_str(),&fileInfo)==-1) {
+		  cout<<"error! fail to open .sample file: "<<str<<endl;
+		  error=true;
+		}
 	}
 
 	if ((!file_geno.empty() || !file_gene.empty()) ) {
 		str=file_pheno;
-		if (stat(str.c_str(),&fileInfo)==-1) {cout<<"error! fail to open phenotype file: "<<str<<endl; error=true;}
+		if (stat(str.c_str(),&fileInfo)==-1) {
+		  cout<<"error! fail to open phenotype file: "<<str<<endl;
+		  error=true;
+		}
 	}
 
 	str=file_geno;
-	if (!str.empty() && stat(str.c_str(),&fileInfo)==-1 ) {cout<<"error! fail to open mean genotype file: "<<str<<endl; error=true;}
+	if (!str.empty() && stat(str.c_str(),&fileInfo)==-1 ) {
+	  cout<<"error! fail to open mean genotype file: "<<str<<endl;
+	  error=true;
+	}
 
 	str=file_gene;
-	if (!str.empty() && stat(str.c_str(),&fileInfo)==-1 ) {cout<<"error! fail to open gene expression file: "<<str<<endl; error=true;}
+	if (!str.empty() && stat(str.c_str(),&fileInfo)==-1 ) {
+	  cout<<"error! fail to open gene expression file: "<<str<<endl;
+	  error=true;
+	}
 
 	str=file_cat;
-	if (!str.empty() && stat(str.c_str(),&fileInfo)==-1 ) {cout<<"error! fail to open category file: "<<str<<endl; error=true;}
+	if (!str.empty() && stat(str.c_str(),&fileInfo)==-1 ) {
+	  cout<<"error! fail to open category file: "<<str<<endl;
+	  error=true;
+	}
 
 	str=file_mcat;
-	if (!str.empty() && stat(str.c_str(),&fileInfo)==-1 ) {cout<<"error! fail to open mcategory file: "<<str<<endl; error=true;}
+	if (!str.empty() && stat(str.c_str(),&fileInfo)==-1 ) {
+	  cout<<"error! fail to open mcategory file: "<<str<<endl;
+	  error=true;
+	}
 
 	str=file_beta;
-	if (!str.empty() && stat(str.c_str(),&fileInfo)==-1 ) {cout<<"error! fail to open beta file: "<<str<<endl; error=true;}
+	if (!str.empty() && stat(str.c_str(),&fileInfo)==-1 ) {
+	  cout<<"error! fail to open beta file: "<<str<<endl;
+	  error=true;
+	}
 
 	str=file_cor;
-	if (!str.empty() && stat(str.c_str(),&fileInfo)==-1 ) {cout<<"error! fail to open correlation file: "<<str<<endl; error=true;}
+	if (!str.empty() && stat(str.c_str(),&fileInfo)==-1 ) {
+	  cout<<"error! fail to open correlation file: "<<str<<endl;
+	  error=true;
+	}
 
 	if (!file_study.empty()) {
 	  str=file_study+".Vq.txt";
-		if (stat(str.c_str(),&fileInfo)==-1) {cout<<"error! fail to open .Vq.txt file: "<<str<<endl; error=true;}
+		if (stat(str.c_str(),&fileInfo)==-1) {
+		  cout<<"error! fail to open .Vq.txt file: "<<str<<endl;
+		  error=true;
+		}
 		str=file_study+".q.txt";
-		if (stat(str.c_str(),&fileInfo)==-1) {cout<<"error! fail to open .q.txt file: "<<str<<endl; error=true;}
+		if (stat(str.c_str(),&fileInfo)==-1) {
+		  cout<<"error! fail to open .q.txt file: "<<str<<endl;
+		  error=true;
+		}
 		str=file_study+".size.txt";
-		if (stat(str.c_str(),&fileInfo)==-1) {cout<<"error! fail to open .size.txt file: "<<str<<endl; error=true;}
+		if (stat(str.c_str(),&fileInfo)==-1) {
+		  cout<<"error! fail to open .size.txt file: "<<str<<endl;
+		  error=true;
+		}
 	}
 
 	if (!file_ref.empty()) {
 		str=file_ref+".S.txt";
-		if (stat(str.c_str(),&fileInfo)==-1) {cout<<"error! fail to open .S.txt file: "<<str<<endl; error=true;}
+		if (stat(str.c_str(),&fileInfo)==-1) {
+		  cout<<"error! fail to open .S.txt file: "<<str<<endl;
+		  error=true;
+		}
 		str=file_ref+".size.txt";
-		if (stat(str.c_str(),&fileInfo)==-1) {cout<<"error! fail to open .size.txt file: "<<str<<endl; error=true;}
+		if (stat(str.c_str(),&fileInfo)==-1) {
+		  cout<<"error! fail to open .size.txt file: "<<str<<endl;
+		  error=true;
+		}
 	}
 
 	str=file_mstudy;
-	if (!str.empty() && stat(str.c_str(),&fileInfo)==-1 ) {cout<<"error! fail to open mstudy file: "<<str<<endl; error=true;}
+	if (!str.empty() && stat(str.c_str(),&fileInfo)==-1 ) {
+	  cout<<"error! fail to open mstudy file: "<<str<<endl;
+	  error=true;
+	}
 
 	str=file_mref;
-	if (!str.empty() && stat(str.c_str(),&fileInfo)==-1 ) {cout<<"error! fail to open mref file: "<<str<<endl; error=true;}
+	if (!str.empty() && stat(str.c_str(),&fileInfo)==-1 ) {
+	  cout<<"error! fail to open mref file: "<<str<<endl;
+	  error=true;
+	}
 
 	str=file_mgeno;
-	if (!str.empty() && stat(str.c_str(),&fileInfo)==-1 ) {cout<<"error! fail to open mgeno file: "<<str<<endl; error=true;}
+	if (!str.empty() && stat(str.c_str(),&fileInfo)==-1 ) {
+	  cout<<"error! fail to open mgeno file: "<<str<<endl;
+	  error=true;
+	}
 
 	str=file_mbfile;
-	if (!str.empty() && stat(str.c_str(),&fileInfo)==-1 ) {cout<<"error! fail to open mbfile file: "<<str<<endl; error=true;}
+	if (!str.empty() && stat(str.c_str(),&fileInfo)==-1 ) {
+	  cout<<"error! fail to open mbfile file: "<<str<<endl;
+	  error=true;
+	}
 
 	size_t flag=0;
 	if (!file_bfile.empty()) {flag++;}
 	if (!file_geno.empty()) {flag++;}
 	if (!file_gene.empty()) {flag++;}
-	// WJA added
+	
+	// WJA added.
 	if (!file_oxford.empty()) {flag++;}
 
-	if (flag!=1 && a_mode!=15 && a_mode!=27 && a_mode!=28 && a_mode!=43 && a_mode!=5 && a_mode!=61 && a_mode!=62 && a_mode!=63 && a_mode!=66 && a_mode!=67) {
-		cout<<"error! either plink binary files, or bimbam mean genotype files, or gene expression files are required."<<endl; error=true;
+	if (flag!=1 && a_mode!=15 && a_mode!=27 && a_mode!=28 &&
+	    a_mode!=43 && a_mode!=5 && a_mode!=61 && a_mode!=62 &&
+	    a_mode!=63 && a_mode!=66 && a_mode!=67) {
+	  cout<<"error! either plink binary files, or bimbam mean"<<
+	    "genotype files, or gene expression files are required."<<endl;
+	  error=true;
 	}
 
 	if (file_pheno.empty() && (a_mode==43 || a_mode==5) ) {
@@ -668,28 +806,30 @@ void PARAM::CheckParam (void) {
 
 	if (a_mode==61 || a_mode==62) {
 	  if (!file_beta.empty()) {
-	    if ( file_mbfile.empty() && file_bfile.empty() && file_mgeno.empty() && file_geno.empty() && file_mref.empty() && file_ref.empty() ) {
-	      cout<<"error! missing genotype file or ref/mref file."<<endl; error=true;
+	    if ( file_mbfile.empty() && file_bfile.empty() &&
+		 file_mgeno.empty() && file_geno.empty() &&
+		 file_mref.empty() && file_ref.empty() ) {
+	      cout<<"error! missing genotype file or ref/mref file."<<endl;
+	      error=true;
 	    }
 	  } else if (!file_pheno.empty()) {
-	    if (file_kin.empty() && (file_ku.empty()||file_kd.empty()) && file_mk.empty() ) {
+	    if (file_kin.empty() && (file_ku.empty()||file_kd.empty()) &&
+		file_mk.empty() ) {
 	      cout<<"error! missing relatedness file. "<<endl;  error=true;
 	    }
-	    /*
-	  } else if (!file_cor.empty()) {
-	    if (file_beta.empty() ) {
-	      cout<<"error! missing cor file."<<endl; error=true;
-	    }
-	    */
-	  } else if ( (file_mstudy.empty() && file_study.empty()) || (file_mref.empty() && file_ref.empty() )  ) {
-	      cout<<"error! either beta file, or phenotype files or study/ref mstudy/mref files are required."<<endl; error=true;
+	  } else if ( (file_mstudy.empty() && file_study.empty()) ||
+		      (file_mref.empty() && file_ref.empty() )  ) {
+	    cout<<"error! either beta file, or phenotype files or "<<
+	      "study/ref mstudy/mref files are required."<<endl;
+	    error=true;
 	  }
 	}
 
 
 	if (a_mode==63) {
-	  if (file_kin.empty() && (file_ku.empty()||file_kd.empty()) && file_mk.empty() ) {
-	    cout<<"error! missing relatedness file. "<<endl;  error=true;
+	  if (file_kin.empty() && (file_ku.empty()||file_kd.empty()) &&
+	      file_mk.empty() ) {
+	    cout<<"error! missing relatedness file. "<<endl; error=true;
 	  }
 	  if ( file_pheno.empty() ) {
 	    cout<<"error! missing phenotype file."<<endl; error=true;
@@ -697,111 +837,197 @@ void PARAM::CheckParam (void) {
 	}
 
 	if (a_mode==66 || a_mode==67) {
-	  if (file_beta.empty() || ( file_mbfile.empty() && file_bfile.empty() && file_mgeno.empty() && file_geno.empty()) ) {
-	    cout<<"error! missing beta file or genotype file."<<endl; error=true;
+	  if (file_beta.empty() ||
+	      (file_mbfile.empty() && file_bfile.empty() &&
+	       file_mgeno.empty() && file_geno.empty()) ) {
+	    cout<<"error! missing beta file or genotype file."<<endl;
+	    error=true;
 	  }
 	}
 
 
-	if (!file_epm.empty() && file_bfile.empty() && file_geno.empty() ) {cout<<"error! estimated parameter file also requires genotype file."<<endl; error=true;}
-	if (!file_ebv.empty() && file_kin.empty()) {cout<<"error! estimated breeding value file also requires relatedness file."<<endl; error=true;}
+	if (!file_epm.empty() && file_bfile.empty() && file_geno.empty()) {
+	  cout<<"error! estimated parameter file also requires genotype "<<
+	    "file."<<endl;
+	  error=true;
+	}
+	if (!file_ebv.empty() && file_kin.empty()) {
+	  cout<<"error! estimated breeding value file also requires "<<
+	    "relatedness file."<<endl;
+	  error=true;
+	}
 
-	if (!file_log.empty() && pheno_mean!=0) {cout<<"error! either log file or mu value can be provide."<<endl; error=true;}
+	if (!file_log.empty() && pheno_mean!=0) {
+	  cout<<"error! either log file or mu value can be provide."<<endl;
+	  error=true;
+	}
 
 	str=file_snps;
-	if (!str.empty() && stat(str.c_str(),&fileInfo)==-1 ) {cout<<"error! fail to open snps file: "<<str<<endl; error=true;}
+	if (!str.empty() && stat(str.c_str(),&fileInfo)==-1 ) {
+	  cout<<"error! fail to open snps file: "<<str<<endl;
+	  error=true;
+	}
 
 	str=file_log;
-	if (!str.empty() && stat(str.c_str(),&fileInfo)==-1 ) {cout<<"error! fail to open log file: "<<str<<endl; error=true;}
+	if (!str.empty() && stat(str.c_str(),&fileInfo)==-1 ) {
+	  cout<<"error! fail to open log file: "<<str<<endl;
+	  error=true;
+	}
 
 	str=file_anno;
-	if (!str.empty() && stat(str.c_str(),&fileInfo)==-1 ) {cout<<"error! fail to open annotation file: "<<str<<endl; error=true;}
+	if (!str.empty() && stat(str.c_str(),&fileInfo)==-1 ) {
+	  cout<<"error! fail to open annotation file: "<<str<<endl;
+	  error=true;
+	}
 
 	str=file_kin;
-	if (!str.empty() && stat(str.c_str(),&fileInfo)==-1 ) {cout<<"error! fail to open relatedness matrix file: "<<str<<endl; error=true;}
+	if (!str.empty() && stat(str.c_str(),&fileInfo)==-1 ) {
+	  cout<<"error! fail to open relatedness matrix file: "<<str<<endl;
+	  error=true;
+	}
 
 	str=file_mk;
-	if (!str.empty() && stat(str.c_str(),&fileInfo)==-1 ) {cout<<"error! fail to open relatedness matrix file: "<<str<<endl; error=true;}
+	if (!str.empty() && stat(str.c_str(),&fileInfo)==-1 ) {
+	  cout<<"error! fail to open relatedness matrix file: "<<str<<endl;
+	  error=true;
+	}
 
 	str=file_cvt;
-	if (!str.empty() && stat(str.c_str(),&fileInfo)==-1 ) {cout<<"error! fail to open covariates file: "<<str<<endl; error=true;}
+	if (!str.empty() && stat(str.c_str(),&fileInfo)==-1 ) {
+	  cout<<"error! fail to open covariates file: "<<str<<endl;
+	  error=true;
+	}
 
 	str=file_gxe;
-	if (!str.empty() && stat(str.c_str(),&fileInfo)==-1 ) {cout<<"error! fail to open environmental covariate file: "<<str<<endl; error=true;}
+	if (!str.empty() && stat(str.c_str(),&fileInfo)==-1 ) {
+	  cout<<"error! fail to open environmental covariate file: "<<
+	    str<<endl;
+	  error=true;
+	}
 
 	str=file_weight;
-	if (!str.empty() && stat(str.c_str(),&fileInfo)==-1 ) {cout<<"error! fail to open the residual weight file: "<<str<<endl; error=true;}
+	if (!str.empty() && stat(str.c_str(),&fileInfo)==-1 ) {
+	  cout<<"error! fail to open the residual weight file: "<<str<<endl;
+	  error=true;
+	}
 
 	str=file_epm;
-	if (!str.empty() && stat(str.c_str(),&fileInfo)==-1 ) {cout<<"error! fail to open estimated parameter file: "<<str<<endl; error=true;}
+	if (!str.empty() && stat(str.c_str(),&fileInfo)==-1 ) {
+	  cout<<"error! fail to open estimated parameter file: "<<str<<endl;
+	  error=true;
+	}
 
 	str=file_ebv;
-	if (!str.empty() && stat(str.c_str(),&fileInfo)==-1 ) {cout<<"error! fail to open estimated breeding value file: "<<str<<endl; error=true;}
+	if (!str.empty() && stat(str.c_str(),&fileInfo)==-1 ) {
+	  cout<<"error! fail to open estimated breeding value file: "<<
+	    str<<endl;
+	  error=true;
+	}
 
 	str=file_read;
-	if (!str.empty() && stat(str.c_str(),&fileInfo)==-1 ) {cout<<"error! fail to open total read file: "<<str<<endl; error=true;}
+	if (!str.empty() && stat(str.c_str(),&fileInfo)==-1 ) {
+	  cout<<"error! fail to open total read file: "<<str<<endl;
+	  error=true;
+	}
 
-	//check if files are compatible with analysis mode
-	if (k_mode==2 && !file_geno.empty() ) {cout<<"error! use \"-km 1\" when using bimbam mean genotype file. "<<endl; error=true;}
+	// Check if files are compatible with analysis mode.
+	if (k_mode==2 && !file_geno.empty() ) {
+	  cout<<"error! use \"-km 1\" when using bimbam mean genotype "<<
+	    "file. "<<endl;
+	  error=true;
+	}
 
-	if ((a_mode==1 || a_mode==2 || a_mode==3 || a_mode==4 || a_mode==5 || a_mode==31) && (file_kin.empty() && (file_ku.empty()||file_kd.empty())) )  {cout<<"error! missing relatedness file. "<<endl;  error=true;}
+	if ((a_mode==1 || a_mode==2 || a_mode==3 || a_mode==4 ||
+	     a_mode==5 || a_mode==31) &&
+	    (file_kin.empty() && (file_ku.empty()||file_kd.empty()))) {
+	  cout<<"error! missing relatedness file. "<<endl;
+	  error=true;
+	}
 
-	if ((a_mode==43) && file_kin.empty())  {cout<<"error! missing relatedness file. -predict option requires -k option to provide a relatedness file."<<endl;  error=true;}
+	if ((a_mode==43) && file_kin.empty()) {
+	  cout<<"error! missing relatedness file. -predict option requires "<<
+	    "-k option to provide a relatedness file."<<endl;
+	  error=true;
+	}
 
-	if ((a_mode==11 || a_mode==12 || a_mode==13 || a_mode==14 || a_mode==16) && !file_cvt.empty() ) {cout<<"error! -bslmm option does not support covariates files."<<endl; error=true;}
+	if ((a_mode==11 || a_mode==12 || a_mode==13 || a_mode==14 ||
+	     a_mode==16) && !file_cvt.empty()) {
+	  cout<<"error! -bslmm option does not support covariates files."<<
+	    endl;
+	  error=true;
+	}
 
 	if (a_mode==41 || a_mode==42) {
-		if (!file_cvt.empty() ) {cout<<"error! -predict option does not support covariates files."<<endl; error=true;}
-		if (file_epm.empty() ) {cout<<"error! -predict option requires estimated parameter files."<<endl; error=true;}
+		if (!file_cvt.empty() ) {
+		  cout<<"error! -predict option does not support "<<
+		    "covariates files."<<endl;
+		  error=true;
+		}
+		if (file_epm.empty() ) {
+		  cout<<"error! -predict option requires estimated "<<
+		    "parameter files."<<endl;
+		  error=true;
+		}
 	}
 
 	if (file_beta.empty() && (a_mode==27 || a_mode==28) ) {
-		cout<<"error! beta effects file is required."<<endl; error=true;
+		cout<<"error! beta effects file is required."<<endl;
+		error=true;
 	}
 
 	return;
 }
 
-
-
-
-
 void PARAM::CheckData (void) {
-  if(file_oxford.empty())	// WJA NOTE: I added this condition so that covariates can be added through sample, probably not exactly what is wanted
+
+  // WJA NOTE: I added this condition so that covariates can be added
+  // through sample, probably not exactly what is wanted.
+  if(file_oxford.empty())	
 	{
- 	if ((file_cvt).empty() || (indicator_cvt).size()==0) {
- 		n_cvt=1;
- 	}
+	  if ((file_cvt).empty() || (indicator_cvt).size()==0) {
+	    n_cvt=1;
+	  }
 	}
 
   if ( (a_mode==66 || a_mode==67) && (v_pve.size()!=n_vc))  {
-	    cout<<"error! the number of pve estimates does not equal to the number of categories in the cat file:"<<v_pve.size()<<" "<<n_vc<<endl; error=true;
-	}
-
-	if ( (indicator_cvt).size()!=0 && (indicator_cvt).size()!=(indicator_idv).size()) {
-		error=true;
-		cout<<"error! number of rows in the covariates file do not match the number of individuals. "<<endl;
-		return;
-	}
-	if ( (indicator_gxe).size()!=0 && (indicator_gxe).size()!=(indicator_idv).size()) {
-		error=true;
-		cout<<"error! number of rows in the gxe file do not match the number of individuals. "<<endl;
-		return;
-	}
-	if ( (indicator_weight).size()!=0 && (indicator_weight).size()!=(indicator_idv).size()) {
-		error=true;
-		cout<<"error! number of rows in the weight file do not match the number of individuals. "<<endl;
-		return;
-	}
-
-	if ( (indicator_read).size()!=0 && (indicator_read).size()!=(indicator_idv).size()) {
-		error=true;
-		cout<<"error! number of rows in the total read file do not match the number of individuals. "<<endl;
-		return;
-	}
+    cout<<"error! the number of pve estimates does not equal to "<<
+      "the number of categories in the cat file:"<<v_pve.size()<<" "<<
+      n_vc<<endl;
+    error=true;
+  }
 
-	//calculate ni_total and ni_test, and set indicator_idv to 0 whenever indicator_cvt=0
-	//and calculate np_obs and np_miss
+  if ( (indicator_cvt).size()!=0 &&
+       (indicator_cvt).size()!=(indicator_idv).size()) {
+    error=true;
+    cout << "error! number of rows in the covariates file do not "<<
+      "match the number of individuals. "<<endl;
+    return;
+  }
+  if ( (indicator_gxe).size()!=0 && (indicator_gxe).size() !=
+       (indicator_idv).size()) {
+    error=true;
+    cout<<"error! number of rows in the gxe file do not match the number "<<
+      "of individuals. "<<endl;
+    return;
+  }
+  if ( (indicator_weight).size()!=0 &&
+       (indicator_weight).size()!=(indicator_idv).size()) {
+    error=true;
+    cout<<"error! number of rows in the weight file do not match "<<
+      "the number of individuals. "<<endl;
+    return;
+  }
+  
+  if ( (indicator_read).size()!=0 &&
+       (indicator_read).size()!=(indicator_idv).size()) {
+    error=true;
+    cout<<"error! number of rows in the total read file do not "<<
+      "match the number of individuals. "<<endl;
+    return;
+  }
+  
+	// Calculate ni_total and ni_test, and set indicator_idv to 0
+	// whenever indicator_cvt=0, and calculate np_obs and np_miss.
 	ni_total=(indicator_idv).size();
 
 	ni_test=0;
@@ -839,24 +1065,8 @@ void PARAM::CheckData (void) {
 		}
 	}
 
-	/*
-	if ((indicator_cvt).size()!=0) {
-		ni_test=0;
-		for (vector<int>::size_type i=0; i<(indicator_idv).size(); ++i) {
-			indicator_idv[i]*=indicator_cvt[i];
-			ni_test+=indicator_idv[i];
-		}
-	}
-
-	if ((indicator_read).size()!=0) {
-		ni_test=0;
-		for (vector<int>::size_type i=0; i<(indicator_idv).size(); ++i) {
-			indicator_idv[i]*=indicator_read[i];
-			ni_test+=indicator_idv[i];
-		}
-	}
-	*/
-	if (ni_test==0 && file_cor.empty() && file_mstudy.empty() && file_study.empty() && file_beta.empty() && file_bf.empty() ) {
+	if (ni_test==0 && file_cor.empty() && file_mstudy.empty() &&
+	    file_study.empty() && file_beta.empty() && file_bf.empty() ) {
 		error=true;
 		cout<<"error! number of analyzed individuals equals 0. "<<endl;
 		return;
@@ -865,23 +1075,27 @@ void PARAM::CheckData (void) {
 	if (a_mode==43) {
 		if (ni_cvt==ni_test) {
 			error=true;
-			cout<<"error! no individual has missing phenotypes."<<endl;
+			cout<<"error! no individual has missing "<<
+			  "phenotypes."<<endl;
 			return;
 		}
 		if ((np_obs+np_miss)!=(ni_cvt*n_ph)) {
 			error=true;
-			//cout<<ni_cvt<<"\t"<<ni_test<<"\t"<<ni_total<<"\t"<<np_obs<<"\t"<<np_miss<<"\t"<<indicator_cvt.size()<<endl;
-			cout<<"error! number of phenotypes do not match the summation of missing and observed phenotypes."<<endl;
+			cout<<"error! number of phenotypes do not match the "<<
+			  "summation of missing and observed phenotypes."<<
+			  endl;
 			return;
 		}
 	}
 
-	//output some information
-	if (file_cor.empty() && file_mstudy.empty() && file_study.empty() && a_mode!=15 && a_mode!=27 && a_mode!=28) {
+	// Output some information.
+	if (file_cor.empty() && file_mstudy.empty() && file_study.empty() &&
+	    a_mode!=15 && a_mode!=27 && a_mode!=28) {
 	  cout<<"## number of total individuals = "<<ni_total<<endl;
 	  if (a_mode==43) {
 	    cout<<"## number of analyzed individuals = "<<ni_cvt<<endl;
-	    cout<<"## number of individuals with full phenotypes = "<<ni_test<<endl;
+	    cout<<"## number of individuals with full phenotypes = "<<
+	      ni_test<<endl;
 	  } else {
 	    cout<<"## number of analyzed individuals = "<<ni_test<<endl;
 	  }
@@ -899,12 +1113,12 @@ void PARAM::CheckData (void) {
 	  } else {}
 	}
 
-	//set d_pace to 1000 for gene expression
+	// Set d_pace to 1000 for gene expression.
 	if (!file_gene.empty() && d_pace==100000) {
 		d_pace=1000;
 	}
 
-	//for case-control studies, count #cases and #controls
+	// For case-control studies, count # cases and # controls.
 	int flag_cc=0;
 	if (a_mode==13) {
 		ni_case=0;
@@ -920,53 +1134,88 @@ void PARAM::CheckData (void) {
 		cout<<"## number of controls = "<<ni_control<<endl;
 	}
 
-	if (flag_cc==1) {cout<<"Unexpected non-binary phenotypes for case/control analysis. Use default (BSLMM) analysis instead."<<endl; a_mode=11;}
+	if (flag_cc==1) {cout<<"Unexpected non-binary phenotypes for "<<
+	    "case/control analysis. Use default (BSLMM) analysis instead."<<
+	    endl;
+	  a_mode=11;
+	}
 
-	//set parameters for BSLMM
-	//and check for predict
+	// Set parameters for BSLMM and check for predict.
 	if (a_mode==11 || a_mode==12 || a_mode==13 || a_mode==14) {
-		if (a_mode==11) {n_mh=1;}
-		if (logp_min==0) {logp_min=-1.0*log((double)ns_test);}
+		if (a_mode==11) {
+		  n_mh=1;
+		}
+		if (logp_min==0) {
+		  logp_min=-1.0*log((double)ns_test);
+		}
 
-		if (h_scale==-1) {h_scale=min(1.0, 10.0/sqrt((double)ni_test) );}
-		if (rho_scale==-1) {rho_scale=min(1.0, 10.0/sqrt((double)ni_test) );}
-		if (logp_scale==-1) {logp_scale=min(1.0, 5.0/sqrt((double)ni_test) );}
+		if (h_scale==-1) {
+		  h_scale=min(1.0, 10.0/sqrt((double)ni_test) );
+		}
+		if (rho_scale==-1) {
+		  rho_scale=min(1.0, 10.0/sqrt((double)ni_test) );
+		}
+		if (logp_scale==-1) {
+		  logp_scale=min(1.0, 5.0/sqrt((double)ni_test) );
+		}
 
 		if (h_min==-1) {h_min=0.0;}
 		if (h_max==-1) {h_max=1.0;}
 
-		if (s_max>ns_test) {s_max=ns_test; cout<<"s_max is re-set to the number of analyzed SNPs."<<endl;}
-		if (s_max<s_min) {cout<<"error! maximum s value must be larger than the minimal value. current values = "<<s_max<<" and "<<s_min<<endl; error=true;}
+		if (s_max>ns_test) {
+		  s_max=ns_test;
+		  cout<<"s_max is re-set to the number of analyzed SNPs."<<
+		    endl;
+		}
+		if (s_max<s_min) {
+		  cout<<"error! maximum s value must be larger than the "<<
+		    "minimal value. current values = "<<s_max<<" and "<<
+		    s_min<<endl;
+		  error=true;
+		}
 	} else if (a_mode==41 || a_mode==42) {
 		if (indicator_bv.size()!=0) {
-			if (indicator_idv.size()!=indicator_bv.size()) {
-				cout<<"error! number of rows in the phenotype file does not match that in the estimated breeding value file: "<<indicator_idv.size()<<"\t"<<indicator_bv.size()<<endl;
-				error=true;
-			} else {
-				size_t flag_bv=0;
-				for (size_t i=0; i<(indicator_bv).size(); ++i) {
-					if (indicator_idv[i]!=indicator_bv[i]) {flag_bv++;}
-				}
-				if (flag_bv!=0) {
-					cout<<"error! individuals with missing value in the phenotype file does not match that in the estimated breeding value file: "<<flag_bv<<endl;
-					error=true;
-				}
-			}
+		  if (indicator_idv.size()!=indicator_bv.size()) {
+		    cout<<"error! number of rows in the "<<
+		      "phenotype file does not match that in the "<<
+		      "estimated breeding value file: "<<
+		      indicator_idv.size()<<"\t"<<indicator_bv.size()<<
+		      endl;
+		    error=true;
+		  } else {
+		    size_t flag_bv=0;
+		    for (size_t i=0; i<(indicator_bv).size(); ++i) {
+		      if (indicator_idv[i]!=indicator_bv[i]) {flag_bv++;}
+		    }
+		    if (flag_bv!=0) {
+		      cout<<"error! individuals with missing value in the "<<
+			"phenotype file does not match that in the "<<
+			"estimated breeding value file: "<<flag_bv<<endl;
+		      error=true;
+		    }
+		  }
 		}
 	}
 
-	if (a_mode==62 && !file_beta.empty() && mapRS2wcat.size()==0) {cout<<"vc analysis with beta files requires -wcat file."<<endl; error=true;}
-	if (a_mode==67 && mapRS2wcat.size()==0) {cout<<"ci analysis with beta files requires -wcat file."<<endl; error=true;}
+	if (a_mode==62 && !file_beta.empty() && mapRS2wcat.size()==0) {
+	  cout<<"vc analysis with beta files requires -wcat file."<<endl;
+	  error=true;
+	}
+	if (a_mode==67 && mapRS2wcat.size()==0) {
+	  cout<<"ci analysis with beta files requires -wcat file."<<endl;
+	  error=true;
+	}
 
-	//file_mk needs to contain more than one line
-	if (n_vc==1 && !file_mk.empty()) {cout<<"error! -mk file should contain more than one line."<<endl; error=true;}
+	// File_mk needs to contain more than one line.
+	if (n_vc==1 && !file_mk.empty()) {
+	  cout<<"error! -mk file should contain more than one line."<<endl;
+	  error=true;
+	}
 
 	return;
 }
 
-
-void PARAM::PrintSummary ()
-{
+void PARAM::PrintSummary () {
 	if (n_ph==1) {
 		cout<<"pve estimate ="<<pve_null<<endl;
 		cout<<"se(pve) ="<<pve_se_null<<endl;
@@ -976,39 +1225,46 @@ void PARAM::PrintSummary ()
 	return;
 }
 
-
-
 void PARAM::ReadGenotypes (gsl_matrix *UtX, gsl_matrix *K, const bool calc_K) {
 	string file_str;
 
 	if (!file_bfile.empty()) {
 		file_str=file_bfile+".bed";
-		if (ReadFile_bed (file_str, indicator_idv, indicator_snp, UtX, K, calc_K)==false) {error=true;}
+		if (ReadFile_bed (file_str, indicator_idv, indicator_snp,
+				  UtX, K, calc_K)==false) {
+		  error=true;
+		}
 	}
 	else {
-		if (ReadFile_geno (file_geno, indicator_idv, indicator_snp, UtX, K, calc_K)==false) {error=true;}
+		if (ReadFile_geno (file_geno, indicator_idv, indicator_snp,
+				   UtX, K, calc_K)==false) {
+		  error=true;
+		}
 	}
 
 	return;
 }
 
-
-void PARAM::ReadGenotypes (vector<vector<unsigned char> > &Xt, gsl_matrix *K, const bool calc_K) {
+void PARAM::ReadGenotypes (vector<vector<unsigned char> > &Xt, gsl_matrix *K,
+			   const bool calc_K) {
 	string file_str;
 
 	if (!file_bfile.empty()) {
 		file_str=file_bfile+".bed";
-		if (ReadFile_bed (file_str, indicator_idv, indicator_snp, Xt, K, calc_K, ni_test, ns_test)==false) {error=true;}
+		if (ReadFile_bed (file_str, indicator_idv, indicator_snp,
+				  Xt, K, calc_K, ni_test, ns_test)==false) {
+		  error=true;
+		}
 	} else {
-		if (ReadFile_geno (file_geno, indicator_idv, indicator_snp, Xt, K, calc_K, ni_test, ns_test)==false) {error=true;}
+		if (ReadFile_geno (file_geno, indicator_idv, indicator_snp,
+				   Xt, K, calc_K, ni_test, ns_test)==false) {
+		  error=true;
+		}
 	}
 
 	return;
 }
 
-
-
-
 void PARAM::CalcKin (gsl_matrix *matrix_kin)  {
 	string file_str;
 
@@ -1016,24 +1272,33 @@ void PARAM::CalcKin (gsl_matrix *matrix_kin)  {
 
 	if (!file_bfile.empty() ) {
 		file_str=file_bfile+".bed";
-		if (PlinkKin (file_str, indicator_snp, a_mode-20, d_pace, matrix_kin)==false) {error=true;}
+		if (PlinkKin (file_str, indicator_snp, a_mode-20, d_pace,
+			      matrix_kin)==false) {
+		  error=true;
+		}
 	}
 	else if (!file_oxford.empty() ) {
 		file_str=file_oxford+".bgen";
-		if (bgenKin (file_str, indicator_snp, a_mode-20, d_pace, matrix_kin)==false) {error=true;}
+		if (bgenKin (file_str, indicator_snp, a_mode-20, d_pace,
+			     matrix_kin)==false) {
+		  error=true;
+		}
  	}
 	else {
 		file_str=file_geno;
-		if (BimbamKin (file_str, indicator_snp, a_mode-20, d_pace, matrix_kin)==false) {error=true;}
+		if (BimbamKin (file_str, indicator_snp, a_mode-20, d_pace,
+			       matrix_kin)==false) {
+		  error=true;
+		}
 	}
 
 	return;
 }
 
-
-
-//from an existing n by nd A and K matrices, compute the d by d S matrix (which is not necessary symmetric)
-void compAKtoS (const gsl_matrix *A, const gsl_matrix *K, const size_t n_cvt, gsl_matrix *S) {
+// From an existing n by nd A and K matrices, compute the d by d S
+// matrix (which is not necessary symmetric).
+void compAKtoS (const gsl_matrix *A, const gsl_matrix *K, const size_t n_cvt,
+		gsl_matrix *S) {
   size_t n_vc=S->size1, ni_test=A->size1;
   double di, dj, tr_AK, sum_A, sum_K, s_A, s_K, sum_AK, tr_A, tr_K, d;
 
@@ -1070,17 +1335,16 @@ void compAKtoS (const gsl_matrix *A, const gsl_matrix *K, const size_t n_cvt, gs
     }
   }
 
-  //eigenlib_invert(Si);
-  //cout<<tr_KiKj<<" "<<s_KiKj<<" "<<sum_Ki<<" "<<sum_Kj<<" "<<si<<" "<<sj<<" "<<d*1000000<<endl;
   return;
 }
 
-
-
-//copied from lmm.cpp; is used in the following function compKtoV
-//map a number 1-(n_cvt+2) to an index between 0 and [(n_c+2)^2+(n_c+2)]/2-1
+// Copied from lmm.cpp; is used in the following function compKtoV
+// map a number 1-(n_cvt+2) to an index between 0 and [(n_c+2)^2+(n_c+2)]/2-1
 size_t GetabIndex (const size_t a, const size_t b, const size_t n_cvt) {
-	if (a>n_cvt+2 || b>n_cvt+2 || a<=0 || b<=0) {cout<<"error in GetabIndex."<<endl; return 0;}
+	if (a>n_cvt+2 || b>n_cvt+2 || a<=0 || b<=0) {
+	  cout<<"error in GetabIndex."<<endl;
+	  return 0;
+	}
 	size_t index;
 	size_t l, h;
 	if (b>a) {l=a; h=b;} else {l=b; h=a;}
@@ -1091,8 +1355,10 @@ size_t GetabIndex (const size_t a, const size_t b, const size_t n_cvt) {
 	return index;
 }
 
-//from an existing n by nd (centered) G matrix, compute the d+1 by d*(d-1)/2*(d+1) Q matrix
-//where inside i'th d+1 by d+1 matrix, each element is tr(KiKlKjKm)-r*tr(KmKiKl)-r*tr(KlKjKm)+r^2*tr(KlKm), where r=n/(n-1)
+// From an existing n by nd (centered) G matrix, compute the d+1 by
+// d*(d-1)/2*(d+1) Q matrix where inside i'th d+1 by d+1 matrix, each
+// element is tr(KiKlKjKm)-r*tr(KmKiKl)-r*tr(KlKjKm)+r^2*tr(KlKm),
+// where r=n/(n-1)
 void compKtoV (const gsl_matrix *G, gsl_matrix *V) {
   size_t n_vc=G->size2/G->size1, ni_test=G->size1;
 
@@ -1103,27 +1369,24 @@ void compKtoV (const gsl_matrix *G, gsl_matrix *V) {
   double d, tr, r=(double)ni_test/(double)(ni_test-1);
   size_t t, t_il, t_jm, t_lm, t_im, t_jl, t_ij;
 
-  //compute KiKj for all pairs of i and j (not including the identity matrix)
+  // Compute KiKj for all pairs of i and j (not including the identity
+  // matrix).
   t=0;
   for (size_t i=0; i<n_vc; i++) {
-    gsl_matrix_const_view Ki=gsl_matrix_const_submatrix(G, 0, i*ni_test, ni_test, ni_test);
+    gsl_matrix_const_view Ki=
+      gsl_matrix_const_submatrix(G, 0, i*ni_test, ni_test, ni_test);
     for (size_t j=i; j<n_vc; j++) {
-      gsl_matrix_const_view Kj=gsl_matrix_const_submatrix(G, 0, j*ni_test, ni_test, ni_test);
-      gsl_matrix_view KiKj_sub=gsl_matrix_submatrix (KiKj, 0, t*ni_test, ni_test, ni_test);
-      eigenlib_dgemm ("N", "N", 1.0, &Ki.matrix, &Kj.matrix, 0.0, &KiKj_sub.matrix);
+      gsl_matrix_const_view Kj=
+	gsl_matrix_const_submatrix(G, 0, j*ni_test, ni_test, ni_test);
+      gsl_matrix_view KiKj_sub=
+	gsl_matrix_submatrix (KiKj, 0, t*ni_test, ni_test, ni_test);
+      eigenlib_dgemm ("N", "N", 1.0, &Ki.matrix, &Kj.matrix, 0.0,
+		      &KiKj_sub.matrix);
       t++;
     }
   }
-  /*
-  for (size_t i=0; i<5; i++) {
-    for (size_t j=0; j<5; j++) {
-      cout<<gsl_matrix_get (G, i, j)<<" ";
-    }
-    cout<<endl;
-  }
-  */
 
-  //compute trKi, trKiKj
+  // Compute trKi, trKiKj.
   t=0;
   for (size_t i=0; i<n_vc; i++) {
     for (size_t j=i; j<n_vc; j++) {
@@ -1143,7 +1406,7 @@ void compKtoV (const gsl_matrix *G, gsl_matrix *V) {
     gsl_vector_set (trKi, i, tr);
   }
 
-  //compute V
+  // Compute V.
   for (size_t i=0; i<n_vc; i++) {
     for (size_t j=i; j<n_vc; j++) {
       t_ij=GetabIndex (i+1, j+1, n_vc-2);
@@ -1153,16 +1416,21 @@ void compKtoV (const gsl_matrix *G, gsl_matrix *V) {
 	    t_il=GetabIndex (i+1, l+1, n_vc-2);
 	    t_jm=GetabIndex (j+1, m+1, n_vc-2);
 	    t_lm=GetabIndex (l+1, m+1, n_vc-2);
-	    //cout<<ni_test<<" "<<r<<t_ij<<" "<<t_il<<" "<<t_jl<<" "<<endl;
 	    tr=0;
 	    for (size_t k=0; k<ni_test; k++) {
-	      gsl_vector_const_view KiKl_row=gsl_matrix_const_subrow (KiKj, k, t_il*ni_test, ni_test);
-	      gsl_vector_const_view KiKl_col=gsl_matrix_const_column (KiKj, t_il*ni_test+k);
-	      gsl_vector_const_view KjKm_row=gsl_matrix_const_subrow (KiKj, k, t_jm*ni_test, ni_test);
-	      gsl_vector_const_view KjKm_col=gsl_matrix_const_column (KiKj, t_jm*ni_test+k);
-
-	      gsl_vector_const_view Kl_row=gsl_matrix_const_subrow (G, k, l*ni_test, ni_test);
-	      gsl_vector_const_view Km_row=gsl_matrix_const_subrow (G, k, m*ni_test, ni_test);
+	      gsl_vector_const_view KiKl_row=
+		gsl_matrix_const_subrow (KiKj, k, t_il*ni_test, ni_test);
+	      gsl_vector_const_view KiKl_col=
+		gsl_matrix_const_column (KiKj, t_il*ni_test+k);
+	      gsl_vector_const_view KjKm_row=
+		gsl_matrix_const_subrow (KiKj, k, t_jm*ni_test, ni_test);
+	      gsl_vector_const_view KjKm_col=
+		gsl_matrix_const_column (KiKj, t_jm*ni_test+k);
+
+	      gsl_vector_const_view Kl_row=
+		gsl_matrix_const_subrow (G, k, l*ni_test, ni_test);
+	      gsl_vector_const_view Km_row=
+		gsl_matrix_const_subrow (G, k, m*ni_test, ni_test);
 
 	      if (i<=l && j<=m) {
 		gsl_blas_ddot (&KiKl_row.vector, &KjKm_col.vector, &d);
@@ -1201,9 +1469,12 @@ void compKtoV (const gsl_matrix *G, gsl_matrix *V) {
 	    t_jl=GetabIndex (j+1, l+1, n_vc-2);
 	    tr=0;
 	    for (size_t k=0; k<ni_test; k++) {
-	      gsl_vector_const_view KiKl_row=gsl_matrix_const_subrow (KiKj, k, t_il*ni_test, ni_test);
-	      gsl_vector_const_view KiKl_col=gsl_matrix_const_column (KiKj, t_il*ni_test+k);
-	      gsl_vector_const_view Kj_row=gsl_matrix_const_subrow (G, k, j*ni_test, ni_test);
+	      gsl_vector_const_view KiKl_row=
+		gsl_matrix_const_subrow (KiKj, k, t_il*ni_test, ni_test);
+	      gsl_vector_const_view KiKl_col=
+		gsl_matrix_const_column (KiKj, t_il*ni_test+k);
+	      gsl_vector_const_view Kj_row=
+		gsl_matrix_const_subrow (G, k, j*ni_test, ni_test);
 
 	      if (i<=l) {
 		gsl_blas_ddot (&KiKl_row.vector, &Kj_row.vector, &d);
@@ -1213,15 +1484,19 @@ void compKtoV (const gsl_matrix *G, gsl_matrix *V) {
 		tr+=d;
 	      }
 	    }
-	    tr+=-r*gsl_vector_get (trKiKj, t_il)-r*gsl_vector_get (trKiKj, t_jl)+r*r*gsl_vector_get (trKi, l);
+	    tr+=-r*gsl_vector_get (trKiKj, t_il) -
+	      r*gsl_vector_get (trKiKj, t_jl)+r*r*gsl_vector_get (trKi, l);
 	  } else if (l==n_vc && m!=n_vc) {
 	    t_jm=GetabIndex (j+1, m+1, n_vc-2);
 	    t_im=GetabIndex (i+1, m+1, n_vc-2);
 	    tr=0;
 	    for (size_t k=0; k<ni_test; k++) {
-	      gsl_vector_const_view KjKm_row=gsl_matrix_const_subrow (KiKj, k, t_jm*ni_test, ni_test);
-	      gsl_vector_const_view KjKm_col=gsl_matrix_const_column (KiKj, t_jm*ni_test+k);
-	      gsl_vector_const_view Ki_row=gsl_matrix_const_subrow (G, k, i*ni_test, ni_test);
+	      gsl_vector_const_view KjKm_row=
+		gsl_matrix_const_subrow (KiKj, k, t_jm*ni_test, ni_test);
+	      gsl_vector_const_view KjKm_col=
+		gsl_matrix_const_column (KiKj, t_jm*ni_test+k);
+	      gsl_vector_const_view Ki_row=
+		gsl_matrix_const_subrow (G, k, i*ni_test, ni_test);
 
 	      if (j<=m) {
 		gsl_blas_ddot (&KjKm_row.vector, &Ki_row.vector, &d);
@@ -1231,9 +1506,12 @@ void compKtoV (const gsl_matrix *G, gsl_matrix *V) {
 		tr+=d;
 	      }
 	    }
-	    tr+=-r*gsl_vector_get (trKiKj, t_im)-r*gsl_vector_get (trKiKj, t_jm)+r*r*gsl_vector_get (trKi, m);
+	    tr+=-r*gsl_vector_get (trKiKj, t_im) -
+	      r*gsl_vector_get (trKiKj, t_jm)+r*r*gsl_vector_get (trKi, m);
 	  } else {
-	    tr=gsl_vector_get (trKiKj, t_ij)-r*gsl_vector_get (trKi, i)-r*gsl_vector_get (trKi, j)+r*r*(double)(ni_test-1);
+	    tr=gsl_vector_get (trKiKj, t_ij) -
+	      r*gsl_vector_get (trKi, i) -
+	      r*gsl_vector_get (trKi, j)+r*r*(double)(ni_test-1);
 	  }
 
 	  gsl_matrix_set (V, l, t_ij*(n_vc+1)+m, tr);
@@ -1251,10 +1529,9 @@ void compKtoV (const gsl_matrix *G, gsl_matrix *V) {
   return;
 }
 
-
-
-//perform Jacknife sampling for variance of S
-void JackknifeAKtoS (const gsl_matrix *W, const gsl_matrix *A, const gsl_matrix *K, gsl_matrix *S, gsl_matrix *Svar) {
+// Perform Jacknife sampling for variance of S.
+void JackknifeAKtoS (const gsl_matrix *W, const gsl_matrix *A,
+		     const gsl_matrix *K, gsl_matrix *S, gsl_matrix *Svar) {
   size_t n_vc=Svar->size1, ni_test=A->size1, n_cvt=W->size2;
 
   vector<vector<vector<double> > > trAK, sumAK;
@@ -1262,10 +1539,7 @@ void JackknifeAKtoS (const gsl_matrix *W, const gsl_matrix *A, const gsl_matrix
   vector<double> vec_tmp;
   double di, dj, d, m, v;
 
-  //gsl_matrix *Stmp=gsl_matrix_alloc (n_vc, ni_test*n_vc);
-  //gsl_matrix *Stmp_sub=gsl_matrix_alloc (n_vc, n_vc);
-
-  //initialize and set all elements to zero
+  // Initialize and set all elements to zero.
   for (size_t i=0; i<ni_test; i++) {
     vec_tmp.push_back(0);
   }
@@ -1284,7 +1558,7 @@ void JackknifeAKtoS (const gsl_matrix *W, const gsl_matrix *A, const gsl_matrix
     sumAK.push_back(sumK);
   }
 
-  //run jackknife
+  // Run jackknife.
   for (size_t i=0; i<n_vc; i++) {
     for (size_t l=0; l<ni_test; l++) {
       for (size_t k=0; k<ni_test; k++) {
@@ -1344,8 +1618,6 @@ void JackknifeAKtoS (const gsl_matrix *W, const gsl_matrix *A, const gsl_matrix
 	  d/=(trA[i][t]-sumA[i][t])*(trK[j][t]-sumK[j][t]);
 	  d-=1/(double)(ni_test-n_cvt-1);
 	}
-	//gsl_matrix_set (Stmp, i, t*n_vc+j, d);
-	//gsl_matrix_set (Stmp, j, t*n_vc+i, d);
 	m+=d; v+=d*d;
       }
       m/=(double)ni_test;
@@ -1361,76 +1633,73 @@ void JackknifeAKtoS (const gsl_matrix *W, const gsl_matrix *A, const gsl_matrix
     }
   }
 
-  /*
-  for (size_t t=0; t<ni_test; t++) {
-    gsl_matrix_view Stmp_view=gsl_matrix_submatrix(Stmp, 0, t*n_vc, n_vc, n_vc);
-    gsl_matrix_memcpy (Stmp_sub, &Stmp_view.matrix);
-    eigenlib_invert(Stmp_sub);
-    gsl_matrix_memcpy (&Stmp_view.matrix, Stmp_sub);
-  }
-
-  for (size_t i=0; i<n_vc; i++) {
-    for (size_t j=i; j<n_vc; j++) {
-      m=0; v=0;
-      for (size_t t=0; t<ni_test; t++) {
-	d=gsl_matrix_get (Stmp, i, t*n_vc+j);
-	m+=d;
-	v+=d*d;
-      }
-      m/=(double)ni_test;
-      v/=(double)ni_test;
-      v-=m*m;
-      v*=(double)(ni_test-1);
-
-      gsl_matrix_set (Svar, i, j, v);
-      d=gsl_matrix_get (Si, i, j);
-      d=(double)ni_test*d-(double)(ni_test-1)*m;
-      gsl_matrix_set (Si, i, j, d);
-      if (i!=j) {gsl_matrix_set (Svar, j, i, v); gsl_matrix_set (Si, j, i, d);}
-    }
-  }
-
-  gsl_matrix_free (Stmp);
-  */
   return;
 }
 
-
-
-//compute the d by d S matrix with its d by d variance matrix of Svar, and the d+1 by d(d+1) matrix of Q for V(q)
-void PARAM::CalcS (const map<string, double> &mapRS2wA, const map<string, double> &mapRS2wK, const gsl_matrix *W, gsl_matrix *A, gsl_matrix *K, gsl_matrix *S, gsl_matrix *Svar, gsl_vector *ns)  {
+// Compute the d by d S matrix with its d by d variance matrix of
+// Svar, and the d+1 by d(d+1) matrix of Q for V(q).
+void PARAM::CalcS (const map<string, double> &mapRS2wA,
+		   const map<string, double> &mapRS2wK,
+		   const gsl_matrix *W, gsl_matrix *A,
+		   gsl_matrix *K, gsl_matrix *S,
+		   gsl_matrix *Svar, gsl_vector *ns)  {
   string file_str;
 
   gsl_matrix_set_zero (S);
   gsl_matrix_set_zero (Svar);
   gsl_vector_set_zero (ns);
 
-  //compute the kinship matrix G for multiple categories; these matrices are not centered, for convienence of Jacknife sampling
+  // Compute the kinship matrix G for multiple categories; these
+  // matrices are not centered, for convienence of Jacknife sampling.
   if (!file_bfile.empty() ) {
     file_str=file_bfile+".bed";
     if (mapRS2wA.size()==0) {
-      if (PlinkKin (file_str, d_pace, indicator_idv, indicator_snp, mapRS2wK, mapRS2cat, snpInfo, W, K, ns)==false) {error=true;}
+      if (PlinkKin (file_str, d_pace, indicator_idv, indicator_snp, mapRS2wK,
+		    mapRS2cat, snpInfo, W, K, ns)==false) {
+	error=true;
+      }
     } else {
-      if (PlinkKin (file_str, d_pace, indicator_idv, indicator_snp, mapRS2wA, mapRS2cat, snpInfo, W, A, ns)==false) {error=true;}
+      if (PlinkKin (file_str, d_pace, indicator_idv, indicator_snp, mapRS2wA,
+		    mapRS2cat, snpInfo, W, A, ns)==false) {
+	error=true;
+      }
     }
   } else if (!file_geno.empty()) {
     file_str=file_geno;
     if (mapRS2wA.size()==0) {
-      if (BimbamKin (file_str, d_pace, indicator_idv, indicator_snp, mapRS2wK, mapRS2cat, snpInfo, W, K, ns)==false) {error=true;}
+      if (BimbamKin (file_str, d_pace, indicator_idv, indicator_snp,
+		     mapRS2wK, mapRS2cat, snpInfo, W, K, ns)==false) {
+	error=true;
+      }
     } else {
-      if (BimbamKin (file_str, d_pace, indicator_idv, indicator_snp, mapRS2wA, mapRS2cat, snpInfo, W, A, ns)==false) {error=true;}
+      if (BimbamKin (file_str, d_pace, indicator_idv, indicator_snp,
+		     mapRS2wA, mapRS2cat, snpInfo, W, A, ns)==false) {
+	error=true;
+      }
     }
   } else if (!file_mbfile.empty() ){
     if (mapRS2wA.size()==0) {
-      if (MFILEKin (1, file_mbfile, d_pace, indicator_idv, mindicator_snp, mapRS2wK, mapRS2cat, msnpInfo, W, K, ns)==false) {error=true;}
+      if (MFILEKin (1, file_mbfile, d_pace, indicator_idv, mindicator_snp,
+		    mapRS2wK, mapRS2cat, msnpInfo, W, K, ns)==false) {
+	error=true;
+      }
     } else {
-      if (MFILEKin (1, file_mbfile, d_pace, indicator_idv, mindicator_snp, mapRS2wA, mapRS2cat, msnpInfo, W, A, ns)==false) {error=true;}
+      if (MFILEKin (1, file_mbfile, d_pace, indicator_idv, mindicator_snp,
+		    mapRS2wA, mapRS2cat, msnpInfo, W, A, ns)==false) {
+	error=true;
+      }
     }
   } else if (!file_mgeno.empty()) {
     if (mapRS2wA.size()==0) {
-      if (MFILEKin (0, file_mgeno, d_pace, indicator_idv, mindicator_snp, mapRS2wK, mapRS2cat, msnpInfo, W, K, ns)==false) {error=true;}
+      if (MFILEKin (0, file_mgeno, d_pace, indicator_idv, mindicator_snp,
+		    mapRS2wK, mapRS2cat, msnpInfo, W, K, ns)==false) {
+	error=true;
+      }
     } else {
-      if (MFILEKin (0, file_mgeno, d_pace, indicator_idv, mindicator_snp, mapRS2wA, mapRS2cat, msnpInfo, W, A, ns)==false) {error=true;}
+      if (MFILEKin (0, file_mgeno, d_pace, indicator_idv, mindicator_snp,
+		    mapRS2wA, mapRS2cat, msnpInfo, W, A, ns)==false) {
+	error=true;
+      }
     }
   }
 
@@ -1438,33 +1707,28 @@ void PARAM::CalcS (const map<string, double> &mapRS2wA, const map<string, double
     gsl_matrix_memcpy (A, K);
   }
 
-  //center and scale every kinship matrix inside G
+  // Center and scale every kinship matrix inside G.
   for (size_t i=0; i<n_vc; i++) {
-    gsl_matrix_view Ksub=gsl_matrix_submatrix(K, 0, i*ni_test, ni_test, ni_test);
+    gsl_matrix_view Ksub=gsl_matrix_submatrix(K,0,i*ni_test,ni_test,ni_test);
     CenterMatrix(&Ksub.matrix);
     ScaleMatrix(&Ksub.matrix);
 
-    gsl_matrix_view Asub=gsl_matrix_submatrix(A, 0, i*ni_test, ni_test, ni_test);
+    gsl_matrix_view Asub=gsl_matrix_submatrix(A,0,i*ni_test,ni_test,ni_test);
     CenterMatrix(&Asub.matrix);
     ScaleMatrix(&Asub.matrix);
   }
 
-  //based on G, compute S
+  // Cased on G, compute S.
   compAKtoS (A, K, W->size2, S);
 
-  //compute Svar and update S with Jacknife
+  // Compute Svar and update S with Jacknife.
   JackknifeAKtoS (W, A, K, S, Svar);
 
-  //based on G, compute a matrix Q that can be used to calculate the variance of q
-  //compKtoV (G, V);
-
   return;
 }
 
-
-
-void PARAM::WriteVector (const gsl_vector *q, const gsl_vector *s, const size_t n_total, const string suffix)
-{
+void PARAM::WriteVector (const gsl_vector *q, const gsl_vector *s,
+			 const size_t n_total, const string suffix) {
 	string file_str;
 	file_str=path_out+"/"+file_out;
 	file_str+=".";
@@ -1472,7 +1736,10 @@ void PARAM::WriteVector (const gsl_vector *q, const gsl_vector *s, const size_t
 	file_str+=".txt";
 
 	ofstream outfile (file_str.c_str(), ofstream::out);
-	if (!outfile) {cout<<"error writing file: "<<file_str.c_str()<<endl; return;}
+	if (!outfile) {
+	  cout<<"error writing file: "<<file_str.c_str()<<endl;
+	  return;
+	}
 
 	outfile.precision(10);
 
@@ -1491,10 +1758,7 @@ void PARAM::WriteVector (const gsl_vector *q, const gsl_vector *s, const size_t
 	return;
 }
 
-
-
-void PARAM::WriteVar (const string suffix)
-{
+void PARAM::WriteVar (const string suffix) {
   string file_str, rs;
 	file_str=path_out+"/"+file_out;
 	file_str+=".";
@@ -1502,7 +1766,10 @@ void PARAM::WriteVar (const string suffix)
 	file_str+=".txt.gz";
 
 	ogzstream outfile (file_str.c_str(), ogzstream::out);
-	if (!outfile) {cout<<"error writing file: "<<file_str.c_str()<<endl; return;}
+	if (!outfile) {
+	  cout<<"error writing file: "<<file_str.c_str()<<endl;
+	  return;
+	}
 
 	outfile.precision(10);
 
@@ -1528,11 +1795,7 @@ void PARAM::WriteVar (const string suffix)
 	return;
 }
 
-
-
-
-void PARAM::WriteMatrix (const gsl_matrix *matrix_U, const string suffix)
-{
+void PARAM::WriteMatrix (const gsl_matrix *matrix_U, const string suffix) {
 	string file_str;
 	file_str=path_out+"/"+file_out;
 	file_str+=".";
@@ -1540,7 +1803,10 @@ void PARAM::WriteMatrix (const gsl_matrix *matrix_U, const string suffix)
 	file_str+=".txt";
 
 	ofstream outfile (file_str.c_str(), ofstream::out);
-	if (!outfile) {cout<<"error writing file: "<<file_str.c_str()<<endl; return;}
+	if (!outfile) {
+	  cout<<"error writing file: "<<file_str.c_str()<<endl;
+	  return;
+	}
 
 	outfile.precision(10);
 
@@ -1556,9 +1822,7 @@ void PARAM::WriteMatrix (const gsl_matrix *matrix_U, const string suffix)
 	return;
 }
 
-
-void PARAM::WriteVector (const gsl_vector *vector_D, const string suffix)
-{
+void PARAM::WriteVector (const gsl_vector *vector_D, const string suffix) {
 	string file_str;
 	file_str=path_out+"/"+file_out;
 	file_str+=".";
@@ -1566,7 +1830,10 @@ void PARAM::WriteVector (const gsl_vector *vector_D, const string suffix)
 	file_str+=".txt";
 
 	ofstream outfile (file_str.c_str(), ofstream::out);
-	if (!outfile) {cout<<"error writing file: "<<file_str.c_str()<<endl; return;}
+	if (!outfile) {
+	  cout<<"error writing file: "<<file_str.c_str()<<endl;
+	  return;
+	}
 
 	outfile.precision(10);
 
@@ -1579,9 +1846,7 @@ void PARAM::WriteVector (const gsl_vector *vector_D, const string suffix)
 	return;
 }
 
-
-void PARAM::CheckCvt ()
-{
+void PARAM::CheckCvt () {
 	if (indicator_cvt.size()==0) {return;}
 
 	size_t ci_test=0;
@@ -1600,22 +1865,25 @@ void PARAM::CheckCvt ()
 	double v_min, v_max;
 	set<size_t> set_remove;
 
-	//check if any columns is an intercept
+	// Check if any columns is an intercept.
 	for (size_t i=0; i<W->size2; i++) {
 		gsl_vector_view w_col=gsl_matrix_column (W, i);
 		gsl_vector_minmax (&w_col.vector, &v_min, &v_max);
 		if (v_min==v_max) {flag_ipt=1; set_remove.insert (i);}
 	}
 
-	//add an intecept term if needed
+	// Add an intecept term if needed.
 	if (n_cvt==set_remove.size()) {
 		indicator_cvt.clear();
 		n_cvt=1;
 	} else if (flag_ipt==0) {
-		cout<<"no intecept term is found in the cvt file. a column of 1s is added."<<endl;
+	  cout<<"no intecept term is found in the cvt file. "<<
+	    "a column of 1s is added."<<endl;
 		for (vector<int>::size_type i=0; i<indicator_idv.size(); ++i) {
-			if (indicator_idv[i]==0 || indicator_cvt[i]==0) {continue;}
-			cvt[i].push_back(1.0);
+		  if (indicator_idv[i]==0 || indicator_cvt[i]==0) {
+		    continue;
+		  }
+		  cvt[i].push_back(1.0);
 		}
 
 		n_cvt++;
@@ -1626,11 +1894,10 @@ void PARAM::CheckCvt ()
 	return;
 }
 
-
-//post-process phentoypes, covariates
-void PARAM::ProcessCvtPhen ()
-{
-	//convert indicator_pheno to indicator_idv
+// Post-process phentoypes and covariates.
+void PARAM::ProcessCvtPhen () {
+  
+ 	// Convert indicator_pheno to indicator_idv.
 	int k=1;
 	indicator_idv.clear();
 	for (size_t i=0; i<indicator_pheno.size(); i++) {
@@ -1641,42 +1908,49 @@ void PARAM::ProcessCvtPhen ()
 		indicator_idv.push_back(k);
 	}
 
-	//remove individuals with missing covariates
+	// Remove individuals with missing covariates.
 	if ((indicator_cvt).size()!=0) {
-		for (vector<int>::size_type i=0; i<(indicator_idv).size(); ++i) {
+		for (vector<int>::size_type i=0;
+		     i<(indicator_idv).size();
+		     ++i) {
 			indicator_idv[i]*=indicator_cvt[i];
 		}
 	}
 
-	//remove individuals with missing gxe variables
+	// Remove individuals with missing gxe variables.
 	if ((indicator_gxe).size()!=0) {
-		for (vector<int>::size_type i=0; i<(indicator_idv).size(); ++i) {
+		for (vector<int>::size_type i=0;
+		     i<(indicator_idv).size();
+		     ++i) {
 			indicator_idv[i]*=indicator_gxe[i];
 		}
 	}
 
-	//remove individuals with missing residual weights
+	// Remove individuals with missing residual weights.
 	if ((indicator_weight).size()!=0) {
-		for (vector<int>::size_type i=0; i<(indicator_idv).size(); ++i) {
+		for (vector<int>::size_type i=0;
+		     i<(indicator_idv).size();
+		     ++i) {
 			indicator_idv[i]*=indicator_weight[i];
 		}
 	}
 
-	//obtain ni_test
+	// Obtain ni_test.
 	ni_test=0;
 	for (vector<int>::size_type i=0; i<(indicator_idv).size(); ++i) {
 	    if (indicator_idv[i]==0) {continue;}
 		ni_test++;
 	}
 
-
-
-	//if subsample number is set, perform a random sub-sampling to determine the subsampled ids
+	// If subsample number is set, perform a random sub-sampling
+	// to determine the subsampled ids.
 	if (ni_subsample!=0) {
 	  if (ni_test<ni_subsample) {
-	    cout<<"error! number of subsamples is less than number of analyzed individuals. "<<endl;
+	    cout<<"error! number of subsamples is less than number of"<<
+	      "analyzed individuals. "<<endl;
 	  } else {
-	    //set up random environment
+	    
+	    // Set up random environment.
 	    gsl_rng_env_setup();
 	    gsl_rng *gsl_r;
 	    const gsl_rng_type * gslType;
@@ -1686,12 +1960,13 @@ void PARAM::ProcessCvtPhen ()
 	      time (&rawtime);
 	      tm * ptm = gmtime (&rawtime);
 
-	      randseed = (unsigned) (ptm->tm_hour%24*3600+ptm->tm_min*60+ptm->tm_sec);
+	      randseed = (unsigned)
+		(ptm->tm_hour%24*3600+ptm->tm_min*60+ptm->tm_sec);
 	    }
 	    gsl_r = gsl_rng_alloc(gslType);
 	    gsl_rng_set(gsl_r, randseed);
 
-	    //from ni_test, sub-sample ni_subsample
+	    // From ni_test, sub-sample ni_subsample.
 	    vector<size_t> a, b;
 	    for (size_t i=0; i<ni_subsample; i++) {
               a.push_back(0);
@@ -1700,9 +1975,10 @@ void PARAM::ProcessCvtPhen ()
 	      b.push_back(i);
 	    }
 
-	    gsl_ran_choose (gsl_r, static_cast<void*>(&a[0]), ni_subsample, static_cast<void*>(&b[0]), ni_test, sizeof (size_t) );
+	    gsl_ran_choose (gsl_r, static_cast<void*>(&a[0]), ni_subsample,
+			    static_cast<void*>(&b[0]),ni_test,sizeof (size_t));
 
-	    //re-set indicator_idv and ni_test
+	    // Re-set indicator_idv and ni_test.
 	    int j=0;
 	    for (vector<int>::size_type i=0; i<(indicator_idv).size(); ++i) {
 	      if (indicator_idv[i]==0) {continue;}
@@ -1715,25 +1991,27 @@ void PARAM::ProcessCvtPhen ()
 	  }
 	}
 
-	//check ni_test
+	// Check ni_test.
 	if (ni_test==0 && a_mode!=15) {
 		error=true;
 		cout<<"error! number of analyzed individuals equals 0. "<<endl;
 		return;
 	}
 
-	//check covariates to see if they are correlated with each other, and to see if the intercept term is included
-	//after getting ni_test
-	//add or remove covariates
+	// Check covariates to see if they are correlated with each
+	// other, and to see if the intercept term is included.
+	// After getting ni_test.
+	// Add or remove covariates.
 	if (indicator_cvt.size()!=0) {
 		CheckCvt();
 	} else {
 		vector<double> cvt_row;
 		cvt_row.push_back(1);
 
-		for (vector<int>::size_type i=0; i<(indicator_idv).size(); ++i) {
+		for (vector<int>::size_type i=0;
+		     i<(indicator_idv).size();
+		     ++i) {
 			indicator_cvt.push_back(1);
-
 			cvt.push_back(cvt_row);
 		}
 	}
@@ -1741,11 +2019,7 @@ void PARAM::ProcessCvtPhen ()
 	return;
 }
 
-
-
-
-void PARAM::CopyCvt (gsl_matrix *W)
-{
+void PARAM::CopyCvt (gsl_matrix *W) {
 	size_t ci_test=0;
 
 	for (vector<int>::size_type i=0; i<indicator_idv.size(); ++i) {
@@ -1759,57 +2033,7 @@ void PARAM::CopyCvt (gsl_matrix *W)
 	return;
 }
 
-/*
-//if there is a column contains all 1's, then flag==1; otherwise flag=0
-void PARAM::CopyA (size_t flag, gsl_matrix *A)
-{
-	size_t ci_test=0;
-	string rs;
-	vector<size_t> flag_vec;
-	vector<double> catc;
-
-	for (size_t j=0; j<n_cat; j++) {
-	  flag_vec.push_back(0);
-	}
-
-	for (vector<int>::size_type i=0; i<indicator_snp.size(); ++i) {
-		if (indicator_snp[i]==0) {continue;}
-		rs=snpInfo[i].rs_number;
-
-		if (mapRS2catc.count(rs)==0) {
-		  for (size_t j=0; j<n_cat; j++) {
-		    gsl_matrix_set (A, ci_test, j, 0);
-		  }
-		} else {
-		  for (size_t j=0; j<n_cat; j++) {
-		    gsl_matrix_set (A, ci_test, j, mapRS2catc[rs][j]);
-		  }
-		}
-
-		if (ci_test==0) {
-		  for (size_t j=0; j<n_cat; j++) {
-		    catc.push_back(mapRS2catc[rs][j]);
-		  }
-		} else {
-		  for (size_t j=0; j<n_cat; j++) {
-		    if (catc[j]==mapRS2catc[rs][j]) {flag_vec[j]++;};
-		  }
-		}
-
-		ci_test++;
-	}
-
-	flag=0;
-	for (size_t j=0; j<n_cat; j++) {
-	  if (flag_vec[j]==0) {flag++;}
-	}
-
-	return;
-}
-*/
-
-void PARAM::CopyGxe (gsl_vector *env)
-{
+void PARAM::CopyGxe (gsl_vector *env) {
 	size_t ci_test=0;
 
 	for (vector<int>::size_type i=0; i<indicator_idv.size(); ++i) {
@@ -1821,8 +2045,7 @@ void PARAM::CopyGxe (gsl_vector *env)
 	return;
 }
 
-void PARAM::CopyWeight (gsl_vector *w)
-{
+void PARAM::CopyWeight (gsl_vector *w) {
 	size_t ci_test=0;
 
 	for (vector<int>::size_type i=0; i<indicator_idv.size(); ++i) {
@@ -1834,11 +2057,9 @@ void PARAM::CopyWeight (gsl_vector *w)
 	return;
 }
 
-
-//if flag=0, then use indicator_idv to load W and Y
-//else, use indicator_cvt to load them
-void PARAM::CopyCvtPhen (gsl_matrix *W, gsl_vector *y, size_t flag)
-{
+// If flag=0, then use indicator_idv to load W and Y;
+// else, use indicator_cvt to load them.
+void PARAM::CopyCvtPhen (gsl_matrix *W, gsl_vector *y, size_t flag) {
 	size_t ci_test=0;
 
 	for (vector<int>::size_type i=0; i<indicator_idv.size(); ++i) {
@@ -1859,10 +2080,9 @@ void PARAM::CopyCvtPhen (gsl_matrix *W, gsl_vector *y, size_t flag)
 	return;
 }
 
-//if flag=0, then use indicator_idv to load W and Y
-//else, use indicator_cvt to load them
-void PARAM::CopyCvtPhen (gsl_matrix *W, gsl_matrix *Y, size_t flag)
-{
+// If flag=0, then use indicator_idv to load W and Y;
+// else, use indicator_cvt to load them.
+void PARAM::CopyCvtPhen (gsl_matrix *W, gsl_matrix *Y, size_t flag) {
 	size_t ci_test=0;
 
 	for (vector<int>::size_type i=0; i<indicator_idv.size(); ++i) {
@@ -1885,12 +2105,7 @@ void PARAM::CopyCvtPhen (gsl_matrix *W, gsl_matrix *Y, size_t flag)
 	return;
 }
 
-
-
-
-
-void PARAM::CopyRead (gsl_vector *log_N)
-{
+void PARAM::CopyRead (gsl_vector *log_N) {
 	size_t ci_test=0;
 
 	for (vector<int>::size_type i=0; i<indicator_idv.size(); ++i) {
@@ -1902,10 +2117,8 @@ void PARAM::CopyRead (gsl_vector *log_N)
 	return;
 }
 
-
-
-void PARAM::ObtainWeight (const set<string> &setSnps_beta, map<string, double> &mapRS2wK)
-{
+void PARAM::ObtainWeight (const set<string> &setSnps_beta,
+			  map<string, double> &mapRS2wK) {
   mapRS2wK.clear();
 
   vector<double> wsum, wcount;
@@ -1921,7 +2134,10 @@ void PARAM::ObtainWeight (const set<string> &setSnps_beta, map<string, double> &
       if (indicator_snp[i]==0) {continue;}
 
       rs=snpInfo[i].rs_number;
-      if ( (setSnps_beta.size()==0 || setSnps_beta.count(rs)!=0) && (mapRS2wsnp.size()==0 || mapRS2wsnp.count(rs)!=0) && (mapRS2wcat.size()==0 || mapRS2wcat.count(rs)!=0) && (mapRS2cat.size()==0 || mapRS2cat.count(rs)!=0) ) {
+      if ( (setSnps_beta.size()==0 || setSnps_beta.count(rs)!=0) &&
+	   (mapRS2wsnp.size()==0 || mapRS2wsnp.count(rs)!=0) &&
+	   (mapRS2wcat.size()==0 || mapRS2wcat.count(rs)!=0) &&
+	   (mapRS2cat.size()==0 || mapRS2cat.count(rs)!=0) ) {
 	if (mapRS2wsnp.size()!=0) {
 	  mapRS2wK[rs]=mapRS2wsnp[rs];
 	  if (mapRS2cat.size()==0) {
@@ -1945,7 +2161,10 @@ void PARAM::ObtainWeight (const set<string> &setSnps_beta, map<string, double> &
 	if (indicator_snp[i]==0) {continue;}
 
 	rs=snpInfo[i].rs_number;
-	if ( (setSnps_beta.size()==0 || setSnps_beta.count(rs)!=0) && (mapRS2wsnp.size()==0 || mapRS2wsnp.count(rs)!=0) && (mapRS2wcat.size()==0 || mapRS2wcat.count(rs)!=0) && (mapRS2cat.size()==0 || mapRS2cat.count(rs)!=0) ) {
+	if ((setSnps_beta.size()==0 || setSnps_beta.count(rs)!=0) &&
+	    (mapRS2wsnp.size()==0 || mapRS2wsnp.count(rs)!=0) &&
+	    (mapRS2wcat.size()==0 || mapRS2wcat.count(rs)!=0) &&
+	    (mapRS2cat.size()==0 || mapRS2cat.count(rs)!=0) ) {
 	  if (mapRS2wsnp.size()!=0) {
 	    mapRS2wK[rs]=mapRS2wsnp[rs];
 	    if (mapRS2cat.size()==0) {
@@ -1967,7 +2186,9 @@ void PARAM::ObtainWeight (const set<string> &setSnps_beta, map<string, double> &
       wsum[i]/=wcount[i];
     }
 
-    for (map<string, double>::iterator it=mapRS2wK.begin(); it!=mapRS2wK.end(); ++it) {
+    for (map<string, double>::iterator it=mapRS2wK.begin();
+	 it!=mapRS2wK.end();
+	 ++it) {
       if (mapRS2cat.size()==0) {
 	it->second/=wsum[0];
       } else {
@@ -1978,10 +2199,12 @@ void PARAM::ObtainWeight (const set<string> &setSnps_beta, map<string, double> &
   return;
 }
 
-
-//pve_flag=0 then do not change pve; pve_flag==1, then change pve to 0 if pve < 0 and pve to 1 if pve > 1
-void PARAM::UpdateWeight (const size_t pve_flag, const map<string, double> &mapRS2wK, const size_t ni_test, const gsl_vector *ns, map<string, double> &mapRS2wA)
-{
+// If pve_flag=0 then do not change pve; pve_flag==1, then change pve
+// to 0 if pve < 0 and pve to 1 if pve > 1.
+void PARAM::UpdateWeight (const size_t pve_flag,
+			  const map<string, double> &mapRS2wK,
+			  const size_t ni_test, const gsl_vector *ns,
+			  map<string, double> &mapRS2wA) {
   double d;
   vector<double> wsum, wcount;
 
@@ -1990,7 +2213,9 @@ void PARAM::UpdateWeight (const size_t pve_flag, const map<string, double> &mapR
     wcount.push_back(0.0);
   }
 
-  for (map<string, double>::const_iterator it=mapRS2wK.begin(); it!=mapRS2wK.end(); ++it) {
+  for (map<string, double>::const_iterator it=mapRS2wK.begin();
+       it!=mapRS2wK.end();
+       ++it) {
     d=1;
     for (size_t i=0; i<n_vc; i++) {
       if (v_pve[i]>=1 && pve_flag==1) {
@@ -1998,7 +2223,8 @@ void PARAM::UpdateWeight (const size_t pve_flag, const map<string, double> &mapR
       } else if (v_pve[i]<=0 && pve_flag==1) {
 	d+=0;
       } else {
-	d+=(double)ni_test/gsl_vector_get(ns, i)*mapRS2wcat[it->first][i]*v_pve[i];
+	d+=(double)ni_test/gsl_vector_get(ns, i)*
+	  mapRS2wcat[it->first][i]*v_pve[i];
       }
     }
     mapRS2wA[it->first]=1/(d*d);
@@ -2016,7 +2242,9 @@ void PARAM::UpdateWeight (const size_t pve_flag, const map<string, double> &mapR
     wsum[i]/=wcount[i];
   }
 
-  for (map<string, double>::iterator it=mapRS2wA.begin(); it!=mapRS2wA.end(); ++it) {
+  for (map<string, double>::iterator it=mapRS2wA.begin();
+       it!=mapRS2wA.end();
+       ++it) {
     if (mapRS2cat.size()==0) {
       it->second/=wsum[0];
     } else {
@@ -2026,9 +2254,13 @@ void PARAM::UpdateWeight (const size_t pve_flag, const map<string, double> &mapR
   return;
 }
 
-// this function updates indicator_snp, and save z-scores and other values into vectors
-void PARAM::UpdateSNPnZ (const map<string, double> &mapRS2wA, const map<string, string> &mapRS2A1, const map<string, double> &mapRS2z, gsl_vector *w, gsl_vector *z, vector<size_t> &vec_cat)
-{
+// This function updates indicator_snp, and save z-scores and other
+// values into vectors.
+void PARAM::UpdateSNPnZ (const map<string, double> &mapRS2wA,
+			 const map<string, string> &mapRS2A1,
+			 const map<string, double> &mapRS2z,
+			 gsl_vector *w, gsl_vector *z,
+			 vector<size_t> &vec_cat) {
   gsl_vector_set_zero (w);
   gsl_vector_set_zero (z);
   vec_cat.clear();
@@ -2086,11 +2318,9 @@ void PARAM::UpdateSNPnZ (const map<string, double> &mapRS2wA, const map<string,
   return;
 }
 
-
-
-// this function updates indicator_snp, and save z-scores and other values into vectors
-void PARAM::UpdateSNP (const map<string, double> &mapRS2wA)
-{
+// This function updates indicator_snp, and save z-scores and other
+// values into vectors.
+void PARAM::UpdateSNP (const map<string, double> &mapRS2wA) {
   string rs;
   if (msnpInfo.size()==0) {
     for (size_t i=0; i<snpInfo.size(); i++) {
-- 
cgit v1.2.3