From 9aa9f76cb16a0e71fea19dd125761034aaca4fd4 Mon Sep 17 00:00:00 2001
From: xiangzhou
Date: Mon, 22 Sep 2014 11:12:46 -0400
Subject: version 0.95alpha

---
 gemma.cpp | 1864 -------------------------------------------------------------
 1 file changed, 1864 deletions(-)
 delete mode 100644 gemma.cpp

(limited to 'gemma.cpp')

diff --git a/gemma.cpp b/gemma.cpp
deleted file mode 100644
index b8693a8..0000000
--- a/gemma.cpp
+++ /dev/null
@@ -1,1864 +0,0 @@
-/*
-	Genome-wide Efficient Mixed Model Association (GEMMA)
-    Copyright (C) 2011  Xiang Zhou
-
-    This program is free software: you can redistribute it and/or modify
-    it under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    This program is distributed in the hope that it will be useful,
-    but WITHOUT ANY WARRANTY; without even the implied warranty of
-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-    GNU General Public License for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with this program.  If not, see <http://www.gnu.org/licenses/>.
-*/
-
-#include <iostream>
-#include <fstream>
-#include <string>
-#include <cstring>
-#include <sys/stat.h>
-#include <ctime>
-#include <cmath>
-
-#include "gsl/gsl_vector.h"
-#include "gsl/gsl_matrix.h"
-#include "gsl/gsl_linalg.h"
-#include "gsl/gsl_blas.h"
-#include "gsl/gsl_eigen.h"
-#include "gsl/gsl_cdf.h"
-
-#include "lapack.h"  //for functions EigenDecomp
-
-#ifdef FORCE_FLOAT
-#include "io_float.h"   //for function ReadFile_kin
-#include "gemma_float.h"
-#include "vc_float.h"
-#include "lm_float.h"  //for LM class
-#include "bslmm_float.h"  //for BSLMM class
-#include "lmm_float.h"  //for LMM class, and functions CalcLambda, CalcPve, CalcVgVe
-#include "mvlmm_float.h"  //for MVLMM class
-#include "prdt_float.h"	//for PRDT class
-#include "mathfunc_float.h"	//for a few functions
-#else
-#include "io.h"
-#include "gemma.h"
-#include "vc.h"
-#include "lm.h"
-#include "bslmm.h"
-#include "lmm.h"
-#include "mvlmm.h"
-#include "prdt.h"
-#include "mathfunc.h"
-#endif
-
-
-using namespace std;
-
-
-
-GEMMA::GEMMA(void):	
-version("0.95alpha"), date("08/08/2014"), year("2011")
-{}
-
-void GEMMA::PrintHeader (void)
-{
-	cout<<endl;
-	cout<<"*********************************************************"<<endl;
-	cout<<"  Genome-wide Efficient Mixed Model Association (GEMMA) "<<endl;
-	cout<<"  Version "<<version<<", "<<date<<"                              "<<endl;
-	cout<<"  Visit                                                 "<<endl;
-	cout<<"     http://stephenslab.uchicago.edu/software.html      "<<endl;
-	cout<<"     http://home.uchicago.edu/~xz7/software.html        "<<endl;
-	cout<<"  For Possible Updates                                  "<<endl;
-	cout<<"  (C) "<<year<<" Xiang Zhou                                   "<<endl;
-	cout<<"  GNU General Public License                            "<<endl;
-	cout<<"  For Help, Type ./gemma -h                             "<<endl;
-	cout<<"*********************************************************"<<endl;
-	cout<<endl;
-	
-	return;
-}
-
-
-void GEMMA::PrintLicense (void)
-{
-	cout<<endl;
-	cout<<"The Software Is Distributed Under GNU General Public License, But May Also Require The Following Notifications."<<endl;
-	cout<<endl;
-	
-	cout<<"Including Lapack Routines In The Software May Require The Following Notification:"<<endl;
-	cout<<"Copyright (c) 1992-2010 The University of Tennessee and The University of Tennessee Research Foundation.  All rights reserved."<<endl;
-	cout<<"Copyright (c) 2000-2010 The University of California Berkeley. All rights reserved."<<endl;
-	cout<<"Copyright (c) 2006-2010 The University of Colorado Denver.  All rights reserved."<<endl;	
-	cout<<endl;
-	
-	cout<<"$COPYRIGHT$"<<endl;
-	cout<<"Additional copyrights may follow"<<endl;
-	cout<<"$HEADER$"<<endl;
-	cout<<"Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:"<<endl;
-	cout<<"- Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer."<<endl;
-	cout<<"- Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer listed in this license in the documentation and/or other materials provided with the distribution."<<endl;
-	cout<<"- Neither the name of the copyright holders nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission."<<endl;
-	cout<<"The copyright holders provide no reassurances that the source code provided does not infringe any patent, copyright, or any other "
-		<<"intellectual property rights of third parties.  The copyright holders disclaim any liability to any recipient for claims brought against "
-		<<"recipient by any third party for infringement of that parties intellectual property rights. "<<endl;
-	cout<<"THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS \"AS IS\" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT "
-		<<"LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT "
-		<<"OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT "
-		<<"LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY "
-		<<"THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE "
-		<<"OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE."<<endl;
-	cout<<endl;
-	
-	
-	
-	return;
-}
-
-
-
-void GEMMA::PrintHelp(size_t option)
-{
-	if (option==0) {
-		cout<<endl; 
-		cout<<" GEMMA version "<<version<<", released on "<<date<<endl;
-		cout<<" implemented by Xiang Zhou"<<endl; 
-		cout<<endl;
-		cout<<" type ./gemma -h [num] for detailed helps"<<endl;
-		cout<<" options: " << endl;
-		cout<<" 1: quick guide"<<endl;
-		cout<<" 2: file I/O related"<<endl;
-		cout<<" 3: SNP QC"<<endl;
-		cout<<" 4: calculate relatedness matrix"<<endl;
-		cout<<" 5: perform eigen decomposition"<<endl;
-		cout<<" 6: perform variance component estiamtion"<<endl;
-		cout<<" 7: fit a linear model"<<endl;
-		cout<<" 8: fit a linear mixed model"<<endl;
-		cout<<" 9: fit a multivariate linear mixed model"<<endl;
-		cout<<" 10: fit a Bayesian sparse linear mixed model"<<endl;
-		cout<<" 11: obtain predicted values"<<endl;
-		cout<<" 12: note"<<endl;
-		cout<<endl;
-	}	
-	
-	if (option==1) {
-		cout<<" QUICK GUIDE" << endl;
-		cout<<" to generate a relatedness matrix: "<<endl;
-		cout<<"         ./gemma -bfile [prefix] -gk [num] -o [prefix]"<<endl;
-		cout<<"         ./gemma -g [filename] -p [filename] -gk [num] -o [prefix]"<<endl;
-		cout<<" to perform eigen decomposition of the relatedness matrix: "<<endl;
-		cout<<"         ./gemma -bfile [prefix] -k [filename] -eigen -o [prefix]"<<endl;
-		cout<<"         ./gemma -g [filename] -p [filename] -k [filename] -eigen -o [prefix]"<<endl;
-		cout<<" to estimate variance components: "<<endl;
-		cout<<"         ./gemma -bfile [prefix] -k [filename] -vc -o [prefix]"<<endl;
-		cout<<"         ./gemma -p [filename] -k [filename] -vc -o [prefix]"<<endl;
-		cout<<"         ./gemma -bfile [prefix] -mk [filename] -vc -o [prefix]"<<endl;
-		cout<<"         ./gemma -p [filename] -mk [filename] -vc -o [prefix]"<<endl;
-		cout<<" to fit a linear mixed model: "<<endl;
-		cout<<"         ./gemma -bfile [prefix] -k [filename] -lmm [num] -o [prefix]"<<endl;
-		cout<<"         ./gemma -g [filename] -p [filename] -a [filename] -k [filename] -lmm [num] -o [prefix]"<<endl;	
-		cout<<" to fit a multivariate linear mixed model: "<<endl;
-		cout<<"         ./gemma -bfile [prefix] -k [filename] -lmm [num] -n [num1] [num2] -o [prefix]"<<endl;
-		cout<<"         ./gemma -g [filename] -p [filename] -a [filename] -k [filename] -lmm [num] -n [num1] [num2] -o [prefix]"<<endl;	
-		cout<<" to fit a Bayesian sparse linear mixed model: "<<endl;
-		cout<<"         ./gemma -bfile [prefix] -bslmm [num] -o [prefix]"<<endl;
-		cout<<"         ./gemma -g [filename] -p [filename] -a [filename] -bslmm [num] -o [prefix]"<<endl;
-		cout<<" to obtain predicted values: "<<endl;
-		cout<<"         ./gemma -bfile [prefix] -epm [filename] -emu [filename] -ebv [filename] -k [filename] -predict [num] -o [prefix]"<<endl;
-		cout<<"         ./gemma -g [filename] -p [filename] -epm [filename] -emu [filename] -ebv [filename] -k [filename] -predict [num] -o [prefix]"<<endl;
-		cout<<endl;
-	}
-	
-	if (option==2) {
-		cout<<" FILE I/O RELATED OPTIONS" << endl;
-		cout<<" -bfile    [prefix]       "<<" specify input PLINK binary ped file prefix."<<endl;	
-		cout<<"          requires: *.fam, *.bim and *.bed files"<<endl;	
-		cout<<"          missing value: -9"<<endl;
-		cout<<" -g        [filename]     "<<" specify input BIMBAM mean genotype file name"<<endl;
-		cout<<"          format: rs#1, allele0, allele1, genotype for individual 1, genotype for individual 2, ..."<<endl;	
-		cout<<"                  rs#2, allele0, allele1, genotype for individual 1, genotype for individual 2, ..."<<endl;	
-		cout<<"                  ..."<<endl;	
-		cout<<"          missing value: NA"<<endl;	
-		cout<<" -p        [filename]     "<<" specify input BIMBAM phenotype file name"<<endl;
-		cout<<"          format: phenotype for individual 1"<<endl;	
-		cout<<"                  phenotype for individual 2"<<endl;	
-		cout<<"                  ..."<<endl;
-		cout<<"          missing value: NA"<<endl;	
-		cout<<" -a        [filename]     "<<" specify input BIMBAM SNP annotation file name (optional)"<<endl;	
-		cout<<"          format: rs#1, base_position, chr_number"<<endl;	
-		cout<<"                  rs#2, base_position, chr_number"<<endl;	
-		cout<<"                  ..."<<endl;
-		cout<<" -k        [filename]     "<<" specify input kinship/relatedness matrix file name"<<endl;	
-		cout<<" -mk       [filename]     "<<" specify input file which contains a list of kinship/relatedness matrices"<<endl;	
-		cout<<" -u        [filename]     "<<" specify input file containing the eigen vectors of the kinship/relatedness matrix"<<endl;	
-		cout<<" -d        [filename]     "<<" specify input file containing the eigen values of the kinship/relatedness matrix"<<endl;	
-		cout<<" -c        [filename]     "<<" specify input covariates file name (optional)"<<endl;	
-		cout<<"          format: covariate 1 for individual 1, ... , covariate c for individual 1"<<endl;	
-		cout<<"                  covariate 1 for individual 2, ... , covariate c for individual 2"<<endl;	
-		cout<<"                  ..."<<endl;
-		cout<<"          missing value: NA"<<endl;	
-		cout<<"          note: the intercept (a column of 1s) may need to be included"<<endl;
-		cout<<" -epm      [filename]     "<<" specify input estimated parameter file name"<<endl;
-		cout<<" -en [n1] [n2] [n3] [n4]  "<<" specify values for the input estimated parameter file (with a header)"<<endl;
-		cout<<"          options: n1: rs column number"<<endl;
-		cout<<"                   n2: estimated alpha column number (0 to ignore)"<<endl;
-		cout<<"                   n3: estimated beta column number (0 to ignore)"<<endl;
-		cout<<"                   n4: estimated gamma column number (0 to ignore)"<<endl;
-		cout<<"          default: 2 4 5 6 if -ebv is not specified; 2 0 5 6 if -ebv is specified"<<endl;
-		cout<<" -ebv      [filename]     "<<" specify input estimated random effect (breeding value) file name"<<endl;
-		cout<<"          format: value for individual 1"<<endl;	
-		cout<<"                  value for individual 2"<<endl;	
-		cout<<"                  ..."<<endl;
-		cout<<"          missing value: NA"<<endl;	
-		cout<<" -emu      [filename]     "<<" specify input log file name containing estimated mean"<<endl;
-		cout<<" -mu       [num]          "<<" specify input estimated mean value"<<endl;
-		cout<<" -gene     [filename]     "<<" specify input gene expression file name"<<endl;
-		cout<<"          format: header"<<endl;	
-		cout<<"                  gene1, count for individual 1, count for individual 2, ..."<<endl;	
-		cout<<"                  gene2, count for individual 1, count for individual 2, ..."<<endl;	
-		cout<<"                  ..."<<endl;
-		cout<<"          missing value: not allowed"<<endl;	
-		cout<<" -r        [filename]     "<<" specify input total read count file name"<<endl;
-		cout<<"          format: total read count for individual 1"<<endl;	
-		cout<<"                  total read count for individual 2"<<endl;	
-		cout<<"                  ..."<<endl;
-		cout<<"          missing value: NA"<<endl;	
-		cout<<" -snps     [filename]     "<<" specify input snps file name to only analyze a certain set of snps"<<endl;
-		cout<<"          format: rs#1"<<endl;	
-		cout<<"                  rs#2"<<endl;	
-		cout<<"                  ..."<<endl;
-		cout<<"          missing value: NA"<<endl;	
-		cout<<" -silence                 "<<" silent terminal display"<<endl;
-		cout<<" -km       [num]          "<<" specify input kinship/relatedness file type (default 1)."<<endl;
-		cout<<"          options: 1: \"n by n matrix\" format"<<endl;
-		cout<<"                   2: \"id  id  value\" format"<<endl;
-		cout<<" -n        [num]          "<<" specify phenotype column in the phenotype/*.fam file (optional; default 1)"<<endl;	
-		cout<<" -pace     [num]          "<<" specify terminal display update pace (default 100000 SNPs or 100000 iterations)."<<endl;
-		cout<<" -outdir   [path]         "<<" specify output directory path (default \"./output/\")"<<endl; 
-		cout<<" -o        [prefix]       "<<" specify output file prefix (default \"result\")"<<endl;  
-		cout<<"          output: prefix.cXX.txt or prefix.sXX.txt from kinship/relatedness matrix estimation"<<endl;	
-		cout<<"          output: prefix.assoc.txt and prefix.log.txt form association tests"<<endl;	
-		cout<<endl;
-	}
-	
-	if (option==3) {
-		cout<<" SNP QC OPTIONS" << endl;
-		cout<<" -miss     [num]          "<<" specify missingness threshold (default 0.05)" << endl; 
-		cout<<" -maf      [num]          "<<" specify minor allele frequency threshold (default 0.01)" << endl; 
-		cout<<" -hwe      [num]          "<<" specify HWE test p value threshold (default 0; no test)" << endl; 
-		cout<<" -r2       [num]          "<<" specify r-squared threshold (default 0.9999)" << endl; 
-		cout<<" -notsnp                  "<<" minor allele frequency cutoff is not used" << endl; 
-		cout<<endl;
-	}
-	
-	if (option==4) {
-		cout<<" RELATEDNESS MATRIX CALCULATION OPTIONS" << endl;
-		cout<<" -gk       [num]          "<<" specify which type of kinship/relatedness matrix to generate (default 1)" << endl; 
-		cout<<"          options: 1: centered XX^T/p"<<endl;
-		cout<<"                   2: standardized XX^T/p"<<endl;
-		cout<<"          note: non-polymorphic SNPs are excluded "<<endl;
-		cout<<endl;
-	}
-	
-	if (option==5) {
-		cout<<" EIGEN-DECOMPOSITION OPTIONS" << endl;
-		cout<<" -eigen                   "<<" specify to perform eigen decomposition of the loaded relatedness matrix" << endl; 
-		cout<<endl;
-	}
-
-	if (option==6) {
-		cout<<" VARIANCE COMPONENT ESTIMATION OPTIONS" << endl;
-		cout<<" -vc                      "<<" specify to perform variance component estimation for the loaded relatedness matrix/matrices" << endl; 
-		cout<<endl;
-	}
-	
-	if (option==7) {
-		cout<<" LINEAR MODEL OPTIONS" << endl;		
-		cout<<" -lm       [num]         "<<" specify analysis options (default 1)."<<endl;
-		cout<<"          options: 1: Wald test"<<endl;
-		cout<<"                   2: Likelihood ratio test"<<endl;
-		cout<<"                   3: Score test"<<endl;
-		cout<<"                   4: 1-3"<<endl;
-		cout<<endl;
-	}
-	
-	if (option==8) {
-		cout<<" LINEAR MIXED MODEL OPTIONS" << endl;		
-		cout<<" -lmm      [num]         "<<" specify analysis options (default 1)."<<endl;
-		cout<<"          options: 1: Wald test"<<endl;		
-		cout<<"                   2: Likelihood ratio test"<<endl;
-		cout<<"                   3: Score test"<<endl;
-		cout<<"                   4: 1-3"<<endl;
-		cout<<"                   5: Parameter estimation in the null model only"<<endl;
-		cout<<" -lmin     [num]          "<<" specify minimal value for lambda (default 1e-5)" << endl; 
-		cout<<" -lmax     [num]          "<<" specify maximum value for lambda (default 1e+5)" << endl; 
-		cout<<" -region   [num]          "<<" specify the number of regions used to evaluate lambda (default 10)" << endl; 
-		cout<<endl;
-	}
-	
-	if (option==9) {
-		cout<<" MULTIVARIATE LINEAR MIXED MODEL OPTIONS" << endl;
-		cout<<" -pnr				     "<<" specify the pvalue threshold to use the Newton-Raphson's method (default 0.001)"<<endl;
-		cout<<" -emi				     "<<" specify the maximum number of iterations for the PX-EM method in the null (default 10000)"<<endl;
-		cout<<" -nri				     "<<" specify the maximum number of iterations for the Newton-Raphson's method in the null (default 100)"<<endl;
-		cout<<" -emp				     "<<" specify the precision for the PX-EM method in the null (default 0.0001)"<<endl;
-		cout<<" -nrp				     "<<" specify the precision for the Newton-Raphson's method in the null (default 0.0001)"<<endl;
-		cout<<" -crt				     "<<" specify to output corrected pvalues for these pvalues that are below the -pnr threshold"<<endl;
-		cout<<endl;
-	}
-	
-	if (option==10) {
-		cout<<" MULTI-LOCUS ANALYSIS OPTIONS" << endl;
-		cout<<" -bslmm	  [num]			 "<<" specify analysis options (default 1)."<<endl;
-		cout<<"          options: 1: BSLMM"<<endl;	
-		cout<<"                   2: standard ridge regression/GBLUP (no mcmc)"<<endl;	
-		cout<<"                   3: probit BSLMM (requires 0/1 phenotypes)"<<endl;			
-		
-		cout<<"   MCMC OPTIONS" << endl;
-		cout<<"   Prior" << endl;	
-		cout<<" -hmin     [num]          "<<" specify minimum value for h (default 0)" << endl; 
-		cout<<" -hmax     [num]          "<<" specify maximum value for h (default 1)" << endl; 
-		cout<<" -rmin     [num]          "<<" specify minimum value for rho (default 0)" << endl; 
-		cout<<" -rmax     [num]          "<<" specify maximum value for rho (default 1)" << endl; 
-		cout<<" -pmin     [num]          "<<" specify minimum value for log10(pi) (default log10(1/p), where p is the number of analyzed SNPs )" << endl; 
-		cout<<" -pmax     [num]          "<<" specify maximum value for log10(pi) (default log10(1) )" << endl; 	
-		cout<<" -smin     [num]          "<<" specify minimum value for |gamma| (default 0)" << endl; 
-		cout<<" -smax     [num]          "<<" specify maximum value for |gamma| (default 300)" << endl; 
-		
-		cout<<"   Proposal" << endl;
-		cout<<" -gmean    [num]          "<<" specify the mean for the geometric distribution (default: 2000)" << endl; 
-		cout<<" -hscale   [num]          "<<" specify the step size scale for the proposal distribution of h (value between 0 and 1, default min(10/sqrt(n),1) )" << endl; 
-		cout<<" -rscale   [num]          "<<" specify the step size scale for the proposal distribution of rho (value between 0 and 1, default min(10/sqrt(n),1) )" << endl; 
-		cout<<" -pscale   [num]          "<<" specify the step size scale for the proposal distribution of log10(pi) (value between 0 and 1, default min(5/sqrt(n),1) )" << endl; 
-		
-		cout<<"   Others" << endl;
-		cout<<" -w        [num]          "<<" specify burn-in steps (default 100,000)" << endl; 
-		cout<<" -s        [num]          "<<" specify sampling steps (default 1,000,000)" << endl; 
-		cout<<" -rpace    [num]          "<<" specify recording pace, record one state in every [num] steps (default 10)" << endl; 	
-		cout<<" -wpace    [num]          "<<" specify writing pace, write values down in every [num] recorded steps (default 1000)" << endl; 	
-		cout<<" -seed     [num]          "<<" specify random seed (a random seed is generated by default)" << endl; 	
-		cout<<" -mh       [num]          "<<" specify number of MH steps in each iteration (default 10)" << endl; 
-		cout<<"          requires: 0/1 phenotypes and -bslmm 3 option"<<endl;	
-		cout<<endl;
-	}
-	
-	if (option==11) {
-		cout<<" PREDICTION OPTIONS" << endl;
-		cout<<" -predict  [num]			 "<<" specify prediction options (default 1)."<<endl;
-		cout<<"          options: 1: predict for individuals with missing phenotypes"<<endl;	
-		cout<<"                   2: predict for individuals with missing phenotypes, and convert the predicted values to probability scale. Use only for files fitted with -bslmm 3 option"<<endl;	
-		cout<<endl;
-	}
-	
-	if (option==12) {
-		cout<<" NOTE"<<endl;
-		cout<<" 1. Only individuals with non-missing phenotoypes and covariates will be analyzed."<<endl;
-		cout<<" 2. Missing genotoypes will be repalced with the mean genotype of that SNP."<<endl;
-		cout<<" 3. For lmm analysis, memory should be large enough to hold the relatedness matrix and to perform eigen decomposition."<<endl;
-		cout<<" 4. For multivariate lmm analysis, use a large -pnr for each snp will increase computation time dramatically."<<endl;
-		cout<<" 5. For bslmm analysis, in addition to 3, memory should be large enough to hold the whole genotype matrix."<<endl;
-		cout<<endl;
-	}
-	
-	return;
-}
-
-
-
-void GEMMA::Assign(int argc, char ** argv, PARAM &cPar)
-{
-	string str;
-	
-	for(int i = 1; i < argc; i++) {		
-		if (strcmp(argv[i], "-bfile")==0 || strcmp(argv[i], "--bfile")==0 || strcmp(argv[i], "-b")==0) {
-			if(argv[i+1] == NULL || argv[i+1][0] == '-') {continue;}
-			++i;
-			str.clear();
-			str.assign(argv[i]);
-			cPar.file_bfile=str;
-		}
-		else if (strcmp(argv[i], "-silence")==0) {
-			cPar.mode_silence=true;
-		}
-		else if (strcmp(argv[i], "-g")==0) {
-			if(argv[i+1] == NULL || argv[i+1][0] == '-') {continue;}
-			++i;
-			str.clear();
-			str.assign(argv[i]);
-			cPar.file_geno=str;
-		}
-		else if (strcmp(argv[i], "-p")==0) {
-			if(argv[i+1] == NULL || argv[i+1][0] == '-') {continue;}
-			++i;
-			str.clear();
-			str.assign(argv[i]);
-			cPar.file_pheno=str;
-		}
-		else if (strcmp(argv[i], "-a")==0) {
-			if(argv[i+1] == NULL || argv[i+1][0] == '-') {continue;}
-			++i;
-			str.clear();
-			str.assign(argv[i]);
-			cPar.file_anno=str;
-		}
-		else if (strcmp(argv[i], "-k")==0) {
-			if(argv[i+1] == NULL || argv[i+1][0] == '-') {continue;}
-			++i;
-			str.clear();
-			str.assign(argv[i]);
-			cPar.file_kin=str;
-		}
-		else if (strcmp(argv[i], "-mk")==0) {
-			if(argv[i+1] == NULL || argv[i+1][0] == '-') {continue;}
-			++i;
-			str.clear();
-			str.assign(argv[i]);
-			cPar.file_mk=str;
-		}
-		else if (strcmp(argv[i], "-u")==0) {
-			if(argv[i+1] == NULL || argv[i+1][0] == '-') {continue;}
-			++i;
-			str.clear();
-			str.assign(argv[i]);
-			cPar.file_ku=str;
-		}
-		else if (strcmp(argv[i], "-d")==0) {
-			if(argv[i+1] == NULL || argv[i+1][0] == '-') {continue;}
-			++i;
-			str.clear();
-			str.assign(argv[i]);
-			cPar.file_kd=str;
-		}
-		else if (strcmp(argv[i], "-c")==0) {
-			if(argv[i+1] == NULL || argv[i+1][0] == '-') {continue;}
-			++i;
-			str.clear();
-			str.assign(argv[i]);
-			cPar.file_cvt=str;
-		}
-		else if (strcmp(argv[i], "-epm")==0) {
-			if(argv[i+1] == NULL || argv[i+1][0] == '-') {continue;}
-			++i;
-			str.clear();
-			str.assign(argv[i]);
-			cPar.file_epm=str;
-		}
-		else if (strcmp(argv[i], "-en")==0) {			
-			while (argv[i+1] != NULL && argv[i+1][0] != '-') {
-				++i;
-				str.clear();
-				str.assign(argv[i]);
-				cPar.est_column.push_back(atoi(str.c_str()));
-			}
-		}
-		else if (strcmp(argv[i], "-ebv")==0) {
-			if(argv[i+1] == NULL || argv[i+1][0] == '-') {continue;}
-			++i;
-			str.clear();
-			str.assign(argv[i]);
-			cPar.file_ebv=str;
-		}
-		else if (strcmp(argv[i], "-emu")==0) {
-			if(argv[i+1] == NULL || argv[i+1][0] == '-') {continue;}
-			++i;
-			str.clear();
-			str.assign(argv[i]);
-			cPar.file_log=str;
-		}
-		else if (strcmp(argv[i], "-mu")==0) {
-			if(argv[i+1] == NULL) {continue;}
-			++i;
-			str.clear();
-			str.assign(argv[i]);
-			cPar.pheno_mean=atof(str.c_str());
-		}
-		else if (strcmp(argv[i], "-gene")==0) {
-			if(argv[i+1] == NULL || argv[i+1][0] == '-') {continue;}
-			++i;
-			str.clear();
-			str.assign(argv[i]);
-			cPar.file_gene=str;
-		}
-		else if (strcmp(argv[i], "-r")==0) {
-			if(argv[i+1] == NULL || argv[i+1][0] == '-') {continue;}
-			++i;
-			str.clear();
-			str.assign(argv[i]);
-			cPar.file_read=str;
-		}
-		else if (strcmp(argv[i], "-snps")==0) {
-			if(argv[i+1] == NULL || argv[i+1][0] == '-') {continue;}
-			++i;
-			str.clear();
-			str.assign(argv[i]);
-			cPar.file_snps=str;
-		}
-		else if (strcmp(argv[i], "-km")==0) {
-			if(argv[i+1] == NULL || argv[i+1][0] == '-') {continue;}
-			++i;
-			str.clear();
-			str.assign(argv[i]);
-			cPar.k_mode=atoi(str.c_str());
-		}		
-		else if (strcmp(argv[i], "-n")==0) {
-			(cPar.p_column).clear();
-			while (argv[i+1] != NULL && argv[i+1][0] != '-') {
-				++i;
-				str.clear();
-				str.assign(argv[i]);
-				(cPar.p_column).push_back(atoi(str.c_str()));
-			}
-		}
-		else if (strcmp(argv[i], "-pace")==0) {
-			if(argv[i+1] == NULL || argv[i+1][0] == '-') {continue;}
-			++i;
-			str.clear();
-			str.assign(argv[i]);
-			cPar.d_pace=atoi(str.c_str());
-		}
-		else if (strcmp(argv[i], "-outdir")==0) {
-			if(argv[i+1] == NULL || argv[i+1][0] == '-') {continue;}
-			++i;
-			str.clear();
-			str.assign(argv[i]);
-			cPar.path_out=str;
-		}
-		else if (strcmp(argv[i], "-o")==0) {
-			if(argv[i+1] == NULL || argv[i+1][0] == '-') {continue;}
-			++i;
-			str.clear();
-			str.assign(argv[i]);
-			cPar.file_out=str;
-		}		
-		else if (strcmp(argv[i], "-miss")==0) {
-			if(argv[i+1] == NULL || argv[i+1][0] == '-') {continue;}
-			++i;
-			str.clear();
-			str.assign(argv[i]);
-			cPar.miss_level=atof(str.c_str());
-		}
-		else if (strcmp(argv[i], "-maf")==0) {
-			if(argv[i+1] == NULL || argv[i+1][0] == '-') {continue;}
-			++i;
-			str.clear();
-			str.assign(argv[i]);
-			if (cPar.maf_level!=-1) {cPar.maf_level=atof(str.c_str());}
-		}
-		else if (strcmp(argv[i], "-hwe")==0) {
-			if(argv[i+1] == NULL || argv[i+1][0] == '-') {continue;}
-			++i;
-			str.clear();
-			str.assign(argv[i]);
-			cPar.hwe_level=atof(str.c_str());
-		}
-		else if (strcmp(argv[i], "-r2")==0) {
-			if(argv[i+1] == NULL || argv[i+1][0] == '-') {continue;}
-			++i;
-			str.clear();
-			str.assign(argv[i]);
-			cPar.r2_level=atof(str.c_str());
-		}
-		else if (strcmp(argv[i], "-notsnp")==0) {
-			cPar.maf_level=-1;
-		}
-		else if (strcmp(argv[i], "-gk")==0) {
-			if (cPar.a_mode!=0) {cPar.error=true; cout<<"error! only one of -gk -eigen -vc -lm -lmm -bslmm -predict options is allowed."<<endl; break;}
-			if(argv[i+1] == NULL || argv[i+1][0] == '-') {cPar.a_mode=21; continue;}
-			++i;
-			str.clear();
-			str.assign(argv[i]);
-			cPar.a_mode=20+atoi(str.c_str());
-		}	
-		else if (strcmp(argv[i], "-eigen")==0) {
-			if (cPar.a_mode!=0) {cPar.error=true; cout<<"error! only one of -gk -eigen -vc -lm -lmm -bslmm -predict options is allowed."<<endl; break;}
-			if(argv[i+1] == NULL || argv[i+1][0] == '-') {cPar.a_mode=31; continue;}
-			++i;
-			str.clear();
-			str.assign(argv[i]);
-			cPar.a_mode=30+atoi(str.c_str());
-		}	
-		else if (strcmp(argv[i], "-vc")==0) {
-			if (cPar.a_mode!=0) {cPar.error=true; cout<<"error! only one of -gk -eigen -vc -lm -lmm -bslmm -predict options is allowed."<<endl; break;}
-			if(argv[i+1] == NULL || argv[i+1][0] == '-') {cPar.a_mode=61; continue;}
-			++i;
-			str.clear();
-			str.assign(argv[i]);
-			cPar.a_mode=60+atoi(str.c_str());
-		}	
-		else if (strcmp(argv[i], "-lm")==0) {
-			if (cPar.a_mode!=0) {cPar.error=true; cout<<"error! only one of -gk -eigen -vc -lm -lmm -bslmm -predict options is allowed."<<endl; break;}
-			if(argv[i+1] == NULL || argv[i+1][0] == '-') {cPar.a_mode=51; continue;}
-			++i;
-			str.clear();
-			str.assign(argv[i]);
-			cPar.a_mode=50+atoi(str.c_str());
-		}
-		else if (strcmp(argv[i], "-fa")==0 || strcmp(argv[i], "-lmm")==0) {
-			if (cPar.a_mode!=0) {cPar.error=true; cout<<"error! only one of -gk -eigen -vc -lm -lmm -bslmm -predict options is allowed."<<endl; break;}
-			if(argv[i+1] == NULL || argv[i+1][0] == '-') {cPar.a_mode=1; continue;}
-			++i;
-			str.clear();
-			str.assign(argv[i]);
-			cPar.a_mode=atoi(str.c_str());
-		}
-		else if (strcmp(argv[i], "-lmin")==0) {
-			if(argv[i+1] == NULL || argv[i+1][0] == '-') {continue;}
-			++i;
-			str.clear();
-			str.assign(argv[i]);
-			cPar.l_min=atof(str.c_str());
-		}
-		else if (strcmp(argv[i], "-lmax")==0) {
-			if(argv[i+1] == NULL || argv[i+1][0] == '-') {continue;}
-			++i;
-			str.clear();
-			str.assign(argv[i]);
-			cPar.l_max=atof(str.c_str());
-		}
-		else if (strcmp(argv[i], "-region")==0) {
-			if(argv[i+1] == NULL || argv[i+1][0] == '-') {continue;}
-			++i;
-			str.clear();
-			str.assign(argv[i]);
-			cPar.n_region=atoi(str.c_str());
-		}
-		else if (strcmp(argv[i], "-pnr")==0) {
-			if(argv[i+1] == NULL || argv[i+1][0] == '-') {continue;}
-			++i;
-			str.clear();
-			str.assign(argv[i]);
-			cPar.p_nr=atof(str.c_str());
-		}
-		else if (strcmp(argv[i], "-emi")==0) {
-			if(argv[i+1] == NULL || argv[i+1][0] == '-') {continue;}
-			++i;
-			str.clear();
-			str.assign(argv[i]);
-			cPar.em_iter=atoi(str.c_str());
-		}
-		else if (strcmp(argv[i], "-nri")==0) {
-			if(argv[i+1] == NULL || argv[i+1][0] == '-') {continue;}
-			++i;
-			str.clear();
-			str.assign(argv[i]);
-			cPar.nr_iter=atoi(str.c_str());
-		}
-		else if (strcmp(argv[i], "-emp")==0) {
-			if(argv[i+1] == NULL || argv[i+1][0] == '-') {continue;}
-			++i;
-			str.clear();
-			str.assign(argv[i]);
-			cPar.em_prec=atof(str.c_str());
-		}
-		else if (strcmp(argv[i], "-nrp")==0) {
-			if(argv[i+1] == NULL || argv[i+1][0] == '-') {continue;}
-			++i;
-			str.clear();
-			str.assign(argv[i]);
-			cPar.nr_prec=atof(str.c_str());
-		}
-		else if (strcmp(argv[i], "-crt")==0) {
-			cPar.crt=1;
-		}
-		else if (strcmp(argv[i], "-bslmm")==0) {
-			if (cPar.a_mode!=0) {cPar.error=true; cout<<"error! only one of -gk -eigen -vc -lm -lmm -bslmm -predict options is allowed."<<endl; break;}
-			if(argv[i+1] == NULL || argv[i+1][0] == '-') {cPar.a_mode=11; continue;}
-			++i;
-			str.clear();
-			str.assign(argv[i]);
-			cPar.a_mode=10+atoi(str.c_str());
-		}
-		else if (strcmp(argv[i], "-hmin")==0) {
-			if(argv[i+1] == NULL || argv[i+1][0] == '-') {continue;}
-			++i;
-			str.clear();
-			str.assign(argv[i]);
-			cPar.h_min=atof(str.c_str());
-		}
-		else if (strcmp(argv[i], "-hmax")==0) {
-			if(argv[i+1] == NULL || argv[i+1][0] == '-') {continue;}
-			++i;
-			str.clear();
-			str.assign(argv[i]);
-			cPar.h_max=atof(str.c_str());
-		}
-		else if (strcmp(argv[i], "-rmin")==0) {
-			if(argv[i+1] == NULL || argv[i+1][0] == '-') {continue;}
-			++i;
-			str.clear();
-			str.assign(argv[i]);
-			cPar.rho_min=atof(str.c_str());
-		}
-		else if (strcmp(argv[i], "-rmax")==0) {
-			if(argv[i+1] == NULL || argv[i+1][0] == '-') {continue;}
-			++i;
-			str.clear();
-			str.assign(argv[i]);
-			cPar.rho_max=atof(str.c_str());
-		}
-		else if (strcmp(argv[i], "-pmin")==0) {
-			if(argv[i+1] == NULL) {continue;}
-			++i;
-			str.clear();
-			str.assign(argv[i]);
-			cPar.logp_min=atof(str.c_str())*log(10.0);
-		}
-		else if (strcmp(argv[i], "-pmax")==0) {
-			if(argv[i+1] == NULL) {continue;}
-			++i;
-			str.clear();
-			str.assign(argv[i]);
-			cPar.logp_max=atof(str.c_str())*log(10.0);
-		}
-		else if (strcmp(argv[i], "-smin")==0) {
-			if(argv[i+1] == NULL || argv[i+1][0] == '-') {continue;}
-			++i;
-			str.clear();
-			str.assign(argv[i]);
-			cPar.s_min=atoi(str.c_str());
-		}
-		else if (strcmp(argv[i], "-smax")==0) {
-			if(argv[i+1] == NULL || argv[i+1][0] == '-') {continue;}
-			++i;
-			str.clear();
-			str.assign(argv[i]);
-			cPar.s_max=atoi(str.c_str());
-		}
-		else if (strcmp(argv[i], "-gmean")==0) {
-			if(argv[i+1] == NULL || argv[i+1][0] == '-') {continue;}
-			++i;
-			str.clear();
-			str.assign(argv[i]);
-			cPar.geo_mean=atof(str.c_str());
-		}
-		else if (strcmp(argv[i], "-hscale")==0) {
-			if(argv[i+1] == NULL || argv[i+1][0] == '-') {continue;}
-			++i;
-			str.clear();
-			str.assign(argv[i]);
-			cPar.h_scale=atof(str.c_str());
-		}
-		else if (strcmp(argv[i], "-rscale")==0) {
-			if(argv[i+1] == NULL || argv[i+1][0] == '-') {continue;}
-			++i;
-			str.clear();
-			str.assign(argv[i]);
-			cPar.rho_scale=atof(str.c_str());
-		}
-		else if (strcmp(argv[i], "-pscale")==0) {
-			if(argv[i+1] == NULL || argv[i+1][0] == '-') {continue;}
-			++i;
-			str.clear();
-			str.assign(argv[i]);
-			cPar.logp_scale=atof(str.c_str())*log(10.0);
-		}
-		else if (strcmp(argv[i], "-w")==0) {
-			if(argv[i+1] == NULL || argv[i+1][0] == '-') {continue;}
-			++i;
-			str.clear();
-			str.assign(argv[i]);
-			cPar.w_step=atoi(str.c_str());
-		}
-		else if (strcmp(argv[i], "-s")==0) {
-			if(argv[i+1] == NULL || argv[i+1][0] == '-') {continue;}
-			++i;
-			str.clear();
-			str.assign(argv[i]);
-			cPar.s_step=atoi(str.c_str());
-		}
-		else if (strcmp(argv[i], "-rpace")==0) {
-			if(argv[i+1] == NULL || argv[i+1][0] == '-') {continue;}
-			++i;
-			str.clear();
-			str.assign(argv[i]);
-			cPar.r_pace=atoi(str.c_str());
-		}
-		else if (strcmp(argv[i], "-wpace")==0) {
-			if(argv[i+1] == NULL || argv[i+1][0] == '-') {continue;}
-			++i;
-			str.clear();
-			str.assign(argv[i]);
-			cPar.w_pace=atoi(str.c_str());
-		}
-		else if (strcmp(argv[i], "-seed")==0) {
-			if(argv[i+1] == NULL || argv[i+1][0] == '-') {continue;}
-			++i;
-			str.clear();
-			str.assign(argv[i]);
-			cPar.randseed=atol(str.c_str());
-		}
-		else if (strcmp(argv[i], "-mh")==0) {
-			if(argv[i+1] == NULL || argv[i+1][0] == '-') {continue;}
-			++i;
-			str.clear();
-			str.assign(argv[i]);
-			cPar.n_mh=atoi(str.c_str());
-		}
-		else if (strcmp(argv[i], "-predict")==0) {
-			if (cPar.a_mode!=0) {cPar.error=true; cout<<"error! only one of -gk -eigen -vc -lm -lmm -bslmm -predict options is allowed."<<endl; break;}
-			if(argv[i+1] == NULL || argv[i+1][0] == '-') {cPar.a_mode=41; continue;}
-			++i;
-			str.clear();
-			str.assign(argv[i]);
-			cPar.a_mode=40+atoi(str.c_str());
-		}
-		else {cout<<"error! unrecognized option: "<<argv[i]<<endl; cPar.error=true; continue;}
-	}
-	
-	//change prediction mode to 43, if the epm file is not provided
-	if (cPar.a_mode==41 && cPar.file_epm.empty()) {cPar.a_mode=43;}
-	
-	return;
-}
-
-
-
-void GEMMA::BatchRun (PARAM &cPar) 
-{
-	clock_t time_begin, time_start;
-	time_begin=clock();
-
-	//Read Files
-	cout<<"Reading Files ... "<<endl;
-	cPar.ReadFiles();
-	if (cPar.error==true) {cout<<"error! fail to read files. "<<endl; return;}
-	cPar.CheckData();
-	if (cPar.error==true) {cout<<"error! fail to check data. "<<endl; return;}
-	//Prediction for bslmm	
-	if (cPar.a_mode==41 || cPar.a_mode==42) {
-		gsl_vector *y_prdt;
-		
-		y_prdt=gsl_vector_alloc (cPar.ni_total-cPar.ni_test);
-
-		//set to zero
-		gsl_vector_set_zero (y_prdt);
-		
-		PRDT cPRDT;
-		cPRDT.CopyFromParam(cPar);
-		
-		//add breeding value if needed
-		if (!cPar.file_kin.empty() && !cPar.file_ebv.empty()) {
-			cout<<"Adding Breeding Values ... "<<endl;
-			
-			gsl_matrix *G=gsl_matrix_alloc (cPar.ni_total, cPar.ni_total);
-			gsl_vector *u_hat=gsl_vector_alloc (cPar.ni_test);
-			
-			//read kinship matrix and set u_hat
-			vector<int> indicator_all;
-			size_t c_bv=0;
-			for (size_t i=0; i<cPar.indicator_idv.size(); i++) {
-				indicator_all.push_back(1);
-				if (cPar.indicator_bv[i]==1) {gsl_vector_set(u_hat, c_bv, cPar.vec_bv[i]); c_bv++;}
-			}
-			
-			ReadFile_kin (cPar.file_kin, indicator_all, cPar.mapID2num, cPar.k_mode, cPar.error, G);
-			if (cPar.error==true) {cout<<"error! fail to read kinship/relatedness file. "<<endl; return;}
-			
-			//read u			
-			cPRDT.AddBV(G, u_hat, y_prdt);					
-			
-			gsl_matrix_free(G);
-			gsl_vector_free(u_hat);
-		}
-
-		//add beta
-		if (!cPar.file_bfile.empty()) {
-			cPRDT.AnalyzePlink (y_prdt);
-		}
-		else {
-			cPRDT.AnalyzeBimbam (y_prdt);
-		}
-		
-		//add mu
-		gsl_vector_add_constant(y_prdt, cPar.pheno_mean);
-		
-		//convert y to probability if needed
-		if (cPar.a_mode==42) {
-			double d;
-			for (size_t i=0; i<y_prdt->size; i++) {
-				d=gsl_vector_get(y_prdt, i);
-				d=gsl_cdf_gaussian_P(d, 1.0);
-				gsl_vector_set(y_prdt, i, d);
-			}
-		}
-			
-			
-		cPRDT.CopyToParam(cPar);
-		
-		cPRDT.WriteFiles(y_prdt);
-		
-		gsl_vector_free(y_prdt);
-	}
-	
-	
-	//Prediction with kinship matrix only; for one or more phenotypes
-	if (cPar.a_mode==43) {
-		//first, use individuals with full phenotypes to obtain estimates of Vg and Ve		
-		gsl_matrix *Y=gsl_matrix_alloc (cPar.ni_test, cPar.n_ph);
-		gsl_matrix *W=gsl_matrix_alloc (Y->size1, cPar.n_cvt);		
-		gsl_matrix *G=gsl_matrix_alloc (Y->size1, Y->size1);
-		gsl_matrix *U=gsl_matrix_alloc (Y->size1, Y->size1); 
-		gsl_matrix *UtW=gsl_matrix_alloc (Y->size1, W->size2);
-		gsl_matrix *UtY=gsl_matrix_alloc (Y->size1, Y->size2);
-		gsl_vector *eval=gsl_vector_alloc (Y->size1);
-		
-		gsl_matrix *Y_full=gsl_matrix_alloc (cPar.ni_cvt, cPar.n_ph);
-		gsl_matrix *W_full=gsl_matrix_alloc (Y_full->size1, cPar.n_cvt);
-		//set covariates matrix W and phenotype matrix Y
-		//an intercept should be included in W, 
-		cPar.CopyCvtPhen (W, Y, 0);
-		cPar.CopyCvtPhen (W_full, Y_full, 1);
-				
-		gsl_matrix *Y_hat=gsl_matrix_alloc (Y_full->size1, cPar.n_ph);		
-		gsl_matrix *G_full=gsl_matrix_alloc (Y_full->size1, Y_full->size1);		
-		gsl_matrix *H_full=gsl_matrix_alloc (Y_full->size1*Y_hat->size2, Y_full->size1*Y_hat->size2);
-				
-		//read relatedness matrix G, and matrix G_full
-		ReadFile_kin (cPar.file_kin, cPar.indicator_idv, cPar.mapID2num, cPar.k_mode, cPar.error, G);
-		if (cPar.error==true) {cout<<"error! fail to read kinship/relatedness file. "<<endl; return;}
-		ReadFile_kin (cPar.file_kin, cPar.indicator_cvt, cPar.mapID2num, cPar.k_mode, cPar.error, G_full);
-		if (cPar.error==true) {cout<<"error! fail to read kinship/relatedness file. "<<endl; return;}
-				
-		//center matrix G
-		CenterMatrix (G);
-		CenterMatrix (G_full);
-		
-		//eigen-decomposition and calculate trace_G
-		cout<<"Start Eigen-Decomposition..."<<endl;
-		time_start=clock();	
-		cPar.trace_G=EigenDecomp (G, U, eval, 0);
-		cPar.trace_G=0.0;
-		for (size_t i=0; i<eval->size; i++) {
-			if (gsl_vector_get (eval, i)<1e-10) {gsl_vector_set (eval, i, 0);}
-			cPar.trace_G+=gsl_vector_get (eval, i);
-		}
-		cPar.trace_G/=(double)eval->size;
-		cPar.time_eigen=(clock()-time_start)/(double(CLOCKS_PER_SEC)*60.0);	
-		
-		//calculate UtW and Uty
-		CalcUtX (U, W, UtW);
-		CalcUtX (U, Y, UtY);
-
-		//calculate variance component and beta estimates
-		//and then obtain predicted values
-		if (cPar.n_ph==1) {
-			gsl_vector *beta=gsl_vector_alloc (W->size2);
-			gsl_vector *se_beta=gsl_vector_alloc (W->size2);
-			
-			double lambda, logl, vg, ve;
-			gsl_vector_view UtY_col=gsl_matrix_column (UtY, 0);
-
-			//obtain estimates
-			CalcLambda ('R', eval, UtW, &UtY_col.vector, cPar.l_min, cPar.l_max, cPar.n_region, lambda, logl);
-			CalcLmmVgVeBeta (eval, UtW, &UtY_col.vector, lambda, vg, ve, beta, se_beta);
-
-			cout<<"REMLE estimate for vg in the null model = "<<vg<<endl;
-			cout<<"REMLE estimate for ve in the null model = "<<ve<<endl;
-			cPar.vg_remle_null=vg; cPar.ve_remle_null=ve;
-			
-			//obtain Y_hat from fixed effects
-			gsl_vector_view Yhat_col=gsl_matrix_column (Y_hat, 0);			
-			gsl_blas_dgemv (CblasNoTrans, 1.0, W_full, beta, 0.0, &Yhat_col.vector);
-			
-			//obtain H
-			gsl_matrix_set_identity (H_full);
-			gsl_matrix_scale (H_full, ve);
-			gsl_matrix_scale (G_full, vg);
-			gsl_matrix_add (H_full, G_full);
-			
-			//free matrices			
-			gsl_vector_free(beta);
-			gsl_vector_free(se_beta);
-		} else {			
-			gsl_matrix *Vg=gsl_matrix_alloc (cPar.n_ph, cPar.n_ph);
-			gsl_matrix *Ve=gsl_matrix_alloc (cPar.n_ph, cPar.n_ph);
-			gsl_matrix *B=gsl_matrix_alloc (cPar.n_ph, W->size2);
-			gsl_matrix *se_B=gsl_matrix_alloc (cPar.n_ph, W->size2);
-			
-			//obtain estimates
-			CalcMvLmmVgVeBeta (eval, UtW, UtY, cPar.em_iter, cPar.nr_iter, cPar.em_prec, cPar.nr_prec, cPar.l_min, cPar.l_max, cPar.n_region, Vg, Ve, B, se_B);
-			
-			cout<<"REMLE estimate for Vg in the null model: "<<endl;
-			for (size_t i=0; i<Vg->size1; i++) {
-				for (size_t j=0; j<=i; j++) {
-					cout<<gsl_matrix_get(Vg, i, j)<<"\t";
-				}
-				cout<<endl;
-			}
-			cout<<"REMLE estimate for Ve in the null model: "<<endl;
-			for (size_t i=0; i<Ve->size1; i++) {
-				for (size_t j=0; j<=i; j++) {
-					cout<<gsl_matrix_get(Ve, i, j)<<"\t";
-				}
-				cout<<endl;
-			}
-			cPar.Vg_remle_null.clear();
-			cPar.Ve_remle_null.clear();
-			for (size_t i=0; i<Vg->size1; i++) {
-				for (size_t j=i; j<Vg->size2; j++) {
-					cPar.Vg_remle_null.push_back(gsl_matrix_get (Vg, i, j) );
-					cPar.Ve_remle_null.push_back(gsl_matrix_get (Ve, i, j) );
-				}
-			}
-			
-			//obtain Y_hat from fixed effects
-			gsl_blas_dgemm (CblasNoTrans, CblasTrans, 1.0, W_full, B, 0.0, Y_hat);
-			
-			//obtain H
-			KroneckerSym(G_full, Vg, H_full);
-			for (size_t i=0; i<G_full->size1; i++) {
-				gsl_matrix_view H_sub=gsl_matrix_submatrix (H_full, i*Ve->size1, i*Ve->size2, Ve->size1, Ve->size2);
-				gsl_matrix_add (&H_sub.matrix, Ve);
-			}
-			
-			//free matrices					
-			gsl_matrix_free (Vg);
-			gsl_matrix_free (Ve);
-			gsl_matrix_free (B);
-			gsl_matrix_free (se_B);
-		}
-					
-		PRDT cPRDT;
-		
-		cPRDT.CopyFromParam(cPar);
-		
-		cout<<"Predicting Missing Phentypes ... "<<endl;
-		time_start=clock();	
-		cPRDT.MvnormPrdt(Y_hat, H_full, Y_full);
-		cPar.time_opt=(clock()-time_start)/(double(CLOCKS_PER_SEC)*60.0);	
-
-		cPRDT.WriteFiles(Y_full);
-		
-		gsl_matrix_free(Y);
-		gsl_matrix_free(W);		
-		gsl_matrix_free(G);
-		gsl_matrix_free(U); 
-		gsl_matrix_free(UtW);
-		gsl_matrix_free(UtY);
-		gsl_vector_free(eval);
-		
-		gsl_matrix_free(Y_full);
-		gsl_matrix_free(Y_hat);
-		gsl_matrix_free(W_full);
-		gsl_matrix_free(G_full);		
-		gsl_matrix_free(H_full);
-	}
-	
-	
-	//Generate Kinship matrix
-	if (cPar.a_mode==21 || cPar.a_mode==22) {  
-		cout<<"Calculating Relatedness Matrix ... "<<endl;
-		
-		gsl_matrix *G=gsl_matrix_alloc (cPar.ni_total, cPar.ni_total);
-		
-		time_start=clock();
-		cPar.CalcKin (G);
-		cPar.time_G=(clock()-time_start)/(double(CLOCKS_PER_SEC)*60.0);
-		if (cPar.error==true) {cout<<"error! fail to calculate relatedness matrix. "<<endl; return;}
-		
-		if (cPar.a_mode==21) {
-			cPar.WriteMatrix (G, "cXX");
-		} else {
-			cPar.WriteMatrix (G, "sXX");
-		}
-		
-		gsl_matrix_free (G);
-	}
-	
-	
-	//LM
-	if (cPar.a_mode==51 || cPar.a_mode==52 || cPar.a_mode==53 || cPar.a_mode==54) {  //Fit LM
-		gsl_matrix *Y=gsl_matrix_alloc (cPar.ni_test, cPar.n_ph);
-		gsl_matrix *W=gsl_matrix_alloc (Y->size1, cPar.n_cvt);	
-		
-		//set covariates matrix W and phenotype matrix Y		
-		//an intercept should be included in W, 
-		cPar.CopyCvtPhen (W, Y, 0);
-		
-		//Fit LM or mvLM
-		if (cPar.n_ph==1) {			
-			LM cLm;
-			cLm.CopyFromParam(cPar);
-			
-			gsl_vector_view Y_col=gsl_matrix_column (Y, 0);
-			
-			if (!cPar.file_gene.empty()) {		
-				cLm.AnalyzeGene (W, &Y_col.vector); //y is the predictor, not the phenotype
-			} else if (!cPar.file_bfile.empty()) {
-				cLm.AnalyzePlink (W, &Y_col.vector);
-			} else {
-				cLm.AnalyzeBimbam (W, &Y_col.vector);
-			}
-			
-			cLm.WriteFiles();
-			cLm.CopyToParam(cPar);
-		}
-		/*
-		else {			 
-			MVLM cMvlm;
-			cMvlm.CopyFromParam(cPar);			
-			
-			if (!cPar.file_bfile.empty()) {
-				cMvlm.AnalyzePlink (W, Y);
-			} else {
-				cMvlm.AnalyzeBimbam (W, Y);
-			}
-			
-			cMvlm.WriteFiles();
-			cMvlm.CopyToParam(cPar);
-		}
-		*/
-		//release all matrices and vectors
-		gsl_matrix_free (Y);
-		gsl_matrix_free (W);
-	} 
-
-
-	//VC estimation with one or multiple kinship matrices
-	//REML approach only
-	//if file_kin or file_ku/kd is provided, then a_mode is changed to 5 already, in param.cpp
-	//for one phenotype only; 
-	if (cPar.a_mode==61) {
-		gsl_matrix *Y=gsl_matrix_alloc (cPar.ni_test, cPar.n_ph);
-		gsl_matrix *W=gsl_matrix_alloc (Y->size1, cPar.n_cvt);
-		gsl_matrix *G=gsl_matrix_alloc (Y->size1, Y->size1*cPar.n_vc );
-
-		//set covariates matrix W and phenotype matrix Y		
-		//an intercept should be included in W, 
-		cPar.CopyCvtPhen (W, Y, 0);
-
-		//read kinship matrices
-		if (!(cPar.file_mk).empty()) {
-		  ReadFile_mk (cPar.file_mk, cPar.indicator_idv, cPar.mapID2num, cPar.k_mode, cPar.error, G);
-		  if (cPar.error==true) {cout<<"error! fail to read kinship/relatedness file. "<<endl; return;}
-	
-		  //center matrix G, and obtain v_traceG
-		  double d=0;
-		  (cPar.v_traceG).clear();
-		  for (size_t i=0; i<cPar.n_vc; i++) {
-		    gsl_matrix_view G_sub=gsl_matrix_submatrix (G, 0, i*G->size1, G->size1, G->size1);
-		    CenterMatrix (&G_sub.matrix);
-		    d=0;
-		    for (size_t j=0; j<G->size1; j++) {
-		      d+=gsl_matrix_get (&G_sub.matrix, j, j);
-		    }
-		    d/=(double)G->size1;
-		    (cPar.v_traceG).push_back(d);
-		  }
-		} else if (!(cPar.file_kin).empty()) {
-			ReadFile_kin (cPar.file_kin, cPar.indicator_idv, cPar.mapID2num, cPar.k_mode, cPar.error, G);
-			if (cPar.error==true) {cout<<"error! fail to read kinship/relatedness file. "<<endl; return;}
-						
-			//center matrix G
-			CenterMatrix (G);
-
-			(cPar.v_traceG).clear();
-			double d=0;
-			for (size_t j=0; j<G->size1; j++) {
-			  d+=gsl_matrix_get (G, j, j);
-			}
-			d/=(double)G->size1;
-			(cPar.v_traceG).push_back(d);
-		}
-			/*
-			//eigen-decomposition and calculate trace_G
-			cout<<"Start Eigen-Decomposition..."<<endl;
-			time_start=clock();	
-	
-			if (cPar.a_mode==31) {
-				cPar.trace_G=EigenDecomp (G, U, eval, 1);
-			} else {
-				cPar.trace_G=EigenDecomp (G, U, eval, 0);
-			}
-
-			cPar.trace_G=0.0;
-			for (size_t i=0; i<eval->size; i++) {
-				if (gsl_vector_get (eval, i)<1e-10) {gsl_vector_set (eval, i, 0);}
-				cPar.trace_G+=gsl_vector_get (eval, i);
-			}
-			cPar.trace_G/=(double)eval->size;
-
-			cPar.time_eigen=(clock()-time_start)/(double(CLOCKS_PER_SEC)*60.0);	
-		} else {
-			ReadFile_eigenU (cPar.file_ku, cPar.error, U);
-			if (cPar.error==true) {cout<<"error! fail to read the U file. "<<endl; return;}
-			
-			ReadFile_eigenD (cPar.file_kd, cPar.error, eval);			
-			if (cPar.error==true) {cout<<"error! fail to read the D file. "<<endl; return;}
-			
-			cPar.trace_G=0.0;
-			for (size_t i=0; i<eval->size; i++) {
-				if (gsl_vector_get(eval, i)<1e-10) {gsl_vector_set(eval, i, 0);}
-			  	cPar.trace_G+=gsl_vector_get(eval, i);
-			}
-			cPar.trace_G/=(double)eval->size;
-		}
-		*/
-		//fit multiple variance components
-		if (cPar.n_ph==1) {
-		  //		  if (cPar.n_vc==1) {
-		    /*
-		    //calculate UtW and Uty	
-		    CalcUtX (U, W, UtW);
-		    CalcUtX (U, Y, UtY);
-
-		    gsl_vector_view beta=gsl_matrix_row (B, 0);
-		    gsl_vector_view se_beta=gsl_matrix_row (se_B, 0);
-		    gsl_vector_view UtY_col=gsl_matrix_column (UtY, 0);
-
-		    CalcLambda ('L', eval, UtW, &UtY_col.vector, cPar.l_min, cPar.l_max, cPar.n_region, cPar.l_mle_null, cPar.logl_mle_H0);
-		    CalcLmmVgVeBeta (eval, UtW, &UtY_col.vector, cPar.l_mle_null, cPar.vg_mle_null, cPar.ve_mle_null, &beta.vector, &se_beta.vector);
-
-		    cPar.beta_mle_null.clear();
-		    cPar.se_beta_mle_null.clear();
-		    for (size_t i=0; i<B->size2; i++) {
-		      cPar.beta_mle_null.push_back(gsl_matrix_get(B, 0, i) );
-		      cPar.se_beta_mle_null.push_back(gsl_matrix_get(se_B, 0, i) );
-		    }
-
-		    CalcLambda ('R', eval, UtW, &UtY_col.vector, cPar.l_min, cPar.l_max, cPar.n_region, cPar.l_remle_null, cPar.logl_remle_H0);
-		    CalcLmmVgVeBeta (eval, UtW, &UtY_col.vector, cPar.l_remle_null, cPar.vg_remle_null, cPar.ve_remle_null, &beta.vector, &se_beta.vector);
-		    cPar.beta_remle_null.clear();
-		    cPar.se_beta_remle_null.clear();
-		    for (size_t i=0; i<B->size2; i++) {
-		      cPar.beta_remle_null.push_back(gsl_matrix_get(B, 0, i) );
-		      cPar.se_beta_remle_null.push_back(gsl_matrix_get(se_B, 0, i) );
-		    }
-				
-		    CalcPve (eval, UtW, &UtY_col.vector, cPar.l_remle_null, cPar.trace_G, cPar.pve_null, cPar.pve_se_null);
-		    cPar.PrintSummary();
-				
-		    //calculate and output residuals
-		    if (cPar.a_mode==5) {
-		      gsl_vector *Utu_hat=gsl_vector_alloc (Y->size1);
-		      gsl_vector *Ute_hat=gsl_vector_alloc (Y->size1);
-		      gsl_vector *u_hat=gsl_vector_alloc (Y->size1);
-		      gsl_vector *e_hat=gsl_vector_alloc (Y->size1);
-		      gsl_vector *y_hat=gsl_vector_alloc (Y->size1);
-					
-		      //obtain Utu and Ute
-		      gsl_vector_memcpy (y_hat, &UtY_col.vector);
-		      gsl_blas_dgemv (CblasNoTrans, -1.0, UtW, &beta.vector, 1.0, y_hat);
-		      
-		      double d, u, e;
-		      for (size_t i=0; i<eval->size; i++) {
-			d=gsl_vector_get (eval, i);
-			u=cPar.l_remle_null*d/(cPar.l_remle_null*d+1.0)*gsl_vector_get(y_hat, i);
-			e=1.0/(cPar.l_remle_null*d+1.0)*gsl_vector_get(y_hat, i);
-			gsl_vector_set (Utu_hat, i, u);
-			gsl_vector_set (Ute_hat, i, e);
-		      }
-					
-		      //obtain u and e
-		      gsl_blas_dgemv (CblasNoTrans, 1.0, U, Utu_hat, 0.0, u_hat);
-		      gsl_blas_dgemv (CblasNoTrans, 1.0, U, Ute_hat, 0.0, e_hat);
-		      
-		      //output residuals					
-		      cPar.WriteVector(u_hat, "residU");
-		      cPar.WriteVector(e_hat, "residE");
-		      
-		      gsl_vector_free(u_hat);
-		      gsl_vector_free(e_hat);
-		      gsl_vector_free(y_hat);
-		    }	
-*/	
-		  //		  } else {
-		    gsl_vector_view Y_col=gsl_matrix_column (Y, 0);
-		    VC cVc;
-		    cVc.CopyFromParam(cPar); 
-		    cVc.CalcVCreml (G, W, &Y_col.vector);			
-		    cVc.CopyToParam(cPar);
-
-		    //obtain pve from sigma2
-		    //obtain se_pve from se_sigma2
-		    
-		    //}
-		} 
-
-		
-	}
-	
-	
-	//LMM or mvLMM or Eigen-Decomposition
-	if (cPar.a_mode==1 || cPar.a_mode==2 || cPar.a_mode==3 || cPar.a_mode==4 || cPar.a_mode==5 || cPar.a_mode==31) {  //Fit LMM or mvLMM or eigen
-		gsl_matrix *Y=gsl_matrix_alloc (cPar.ni_test, cPar.n_ph);
-		gsl_matrix *W=gsl_matrix_alloc (Y->size1, cPar.n_cvt);
-		gsl_matrix *B=gsl_matrix_alloc (Y->size2, W->size2);	//B is a d by c matrix
-		gsl_matrix *se_B=gsl_matrix_alloc (Y->size2, W->size2);
-		gsl_matrix *G=gsl_matrix_alloc (Y->size1, Y->size1);
-		gsl_matrix *U=gsl_matrix_alloc (Y->size1, Y->size1); 
-		gsl_matrix *UtW=gsl_matrix_alloc (Y->size1, W->size2);
-		gsl_matrix *UtY=gsl_matrix_alloc (Y->size1, Y->size2);
-		gsl_vector *eval=gsl_vector_alloc (Y->size1);
-				
-		//set covariates matrix W and phenotype matrix Y		
-		//an intercept should be included in W, 
-		cPar.CopyCvtPhen (W, Y, 0);
-				
-		//read relatedness matrix G	
-		if (!(cPar.file_kin).empty()) {
-			ReadFile_kin (cPar.file_kin, cPar.indicator_idv, cPar.mapID2num, cPar.k_mode, cPar.error, G);
-			if (cPar.error==true) {cout<<"error! fail to read kinship/relatedness file. "<<endl; return;}
-						
-			//center matrix G
-			CenterMatrix (G);
-			
-			//eigen-decomposition and calculate trace_G
-			cout<<"Start Eigen-Decomposition..."<<endl;
-			time_start=clock();	
-	
-			if (cPar.a_mode==31) {
-				cPar.trace_G=EigenDecomp (G, U, eval, 1);
-			} else {
-				cPar.trace_G=EigenDecomp (G, U, eval, 0);
-			}
-
-			cPar.trace_G=0.0;
-			for (size_t i=0; i<eval->size; i++) {
-				if (gsl_vector_get (eval, i)<1e-10) {gsl_vector_set (eval, i, 0);}
-				cPar.trace_G+=gsl_vector_get (eval, i);
-			}
-			cPar.trace_G/=(double)eval->size;
-
-			cPar.time_eigen=(clock()-time_start)/(double(CLOCKS_PER_SEC)*60.0);	
-		} else {
-			ReadFile_eigenU (cPar.file_ku, cPar.error, U);
-			if (cPar.error==true) {cout<<"error! fail to read the U file. "<<endl; return;}
-			
-			ReadFile_eigenD (cPar.file_kd, cPar.error, eval);			
-			if (cPar.error==true) {cout<<"error! fail to read the D file. "<<endl; return;}
-			
-			cPar.trace_G=0.0;
-			for (size_t i=0; i<eval->size; i++) {
-				if (gsl_vector_get(eval, i)<1e-10) {gsl_vector_set(eval, i, 0);}
-			  	cPar.trace_G+=gsl_vector_get(eval, i);
-			}
-			cPar.trace_G/=(double)eval->size;
-		}
-		
-		if (cPar.a_mode==31) {
-			cPar.WriteMatrix(U, "eigenU");
-			cPar.WriteVector(eval, "eigenD");
-		} else {
-			//calculate UtW and Uty	
-			CalcUtX (U, W, UtW);
-			CalcUtX (U, Y, UtY);			
-
-			//calculate REMLE/MLE estimate and pve for univariate model
-			if (cPar.n_ph==1) {
-				gsl_vector_view beta=gsl_matrix_row (B, 0);
-				gsl_vector_view se_beta=gsl_matrix_row (se_B, 0);
-				gsl_vector_view UtY_col=gsl_matrix_column (UtY, 0);
-
-				CalcLambda ('L', eval, UtW, &UtY_col.vector, cPar.l_min, cPar.l_max, cPar.n_region, cPar.l_mle_null, cPar.logl_mle_H0);
-				CalcLmmVgVeBeta (eval, UtW, &UtY_col.vector, cPar.l_mle_null, cPar.vg_mle_null, cPar.ve_mle_null, &beta.vector, &se_beta.vector);
-
-				cPar.beta_mle_null.clear();
-				cPar.se_beta_mle_null.clear();
-				for (size_t i=0; i<B->size2; i++) {
-					cPar.beta_mle_null.push_back(gsl_matrix_get(B, 0, i) );
-					cPar.se_beta_mle_null.push_back(gsl_matrix_get(se_B, 0, i) );
-				}
-
-				CalcLambda ('R', eval, UtW, &UtY_col.vector, cPar.l_min, cPar.l_max, cPar.n_region, cPar.l_remle_null, cPar.logl_remle_H0);
-				CalcLmmVgVeBeta (eval, UtW, &UtY_col.vector, cPar.l_remle_null, cPar.vg_remle_null, cPar.ve_remle_null, &beta.vector, &se_beta.vector);
-				cPar.beta_remle_null.clear();
-				cPar.se_beta_remle_null.clear();
-				for (size_t i=0; i<B->size2; i++) {
-					cPar.beta_remle_null.push_back(gsl_matrix_get(B, 0, i) );
-					cPar.se_beta_remle_null.push_back(gsl_matrix_get(se_B, 0, i) );
-				}
-				
-				CalcPve (eval, UtW, &UtY_col.vector, cPar.l_remle_null, cPar.trace_G, cPar.pve_null, cPar.pve_se_null);
-				cPar.PrintSummary();
-				
-				//calculate and output residuals
-				if (cPar.a_mode==5) {
-					gsl_vector *Utu_hat=gsl_vector_alloc (Y->size1);
-					gsl_vector *Ute_hat=gsl_vector_alloc (Y->size1);
-					gsl_vector *u_hat=gsl_vector_alloc (Y->size1);
-					gsl_vector *e_hat=gsl_vector_alloc (Y->size1);
-					gsl_vector *y_hat=gsl_vector_alloc (Y->size1);
-					
-					//obtain Utu and Ute
-					gsl_vector_memcpy (y_hat, &UtY_col.vector);
-					gsl_blas_dgemv (CblasNoTrans, -1.0, UtW, &beta.vector, 1.0, y_hat);
-					
-					double d, u, e;
-					for (size_t i=0; i<eval->size; i++) {
-						d=gsl_vector_get (eval, i);
-						u=cPar.l_remle_null*d/(cPar.l_remle_null*d+1.0)*gsl_vector_get(y_hat, i);
-						e=1.0/(cPar.l_remle_null*d+1.0)*gsl_vector_get(y_hat, i);
-						gsl_vector_set (Utu_hat, i, u);
-						gsl_vector_set (Ute_hat, i, e);
-					}
-					
-					//obtain u and e
-					gsl_blas_dgemv (CblasNoTrans, 1.0, U, Utu_hat, 0.0, u_hat);
-					gsl_blas_dgemv (CblasNoTrans, 1.0, U, Ute_hat, 0.0, e_hat);
-					
-					//output residuals					
-					cPar.WriteVector(u_hat, "residU");
-					cPar.WriteVector(e_hat, "residE");
-					
-					gsl_vector_free(u_hat);
-					gsl_vector_free(e_hat);
-					gsl_vector_free(y_hat);
-				}							
-			} 
-			
-			//Fit LMM or mvLMM
-			if (cPar.a_mode==1 || cPar.a_mode==2 || cPar.a_mode==3 || cPar.a_mode==4) {
-				if (cPar.n_ph==1) {			
-					LMM cLmm;
-					cLmm.CopyFromParam(cPar);
-					
-					gsl_vector_view Y_col=gsl_matrix_column (Y, 0);
-					gsl_vector_view UtY_col=gsl_matrix_column (UtY, 0);
-					
-					if (!cPar.file_gene.empty()) {		
-						cLmm.AnalyzeGene (U, eval, UtW, &UtY_col.vector, W, &Y_col.vector); //y is the predictor, not the phenotype
-					} else if (!cPar.file_bfile.empty()) {
-						cLmm.AnalyzePlink (U, eval, UtW, &UtY_col.vector, W, &Y_col.vector);
-					} else {
-						cLmm.AnalyzeBimbam (U, eval, UtW, &UtY_col.vector, W, &Y_col.vector);
-					}	
-					
-					cLmm.WriteFiles();
-					cLmm.CopyToParam(cPar);
-				} else {			 
-					MVLMM cMvlmm;
-					cMvlmm.CopyFromParam(cPar);			
-					
-					if (!cPar.file_bfile.empty()) {
-						cMvlmm.AnalyzePlink (U, eval, UtW, UtY);
-					} else {
-						cMvlmm.AnalyzeBimbam (U, eval, UtW, UtY);
-					}
-					
-					cMvlmm.WriteFiles();
-					cMvlmm.CopyToParam(cPar);
-				}
-			}
-		}
-		
-				
-		//release all matrices and vectors
-		gsl_matrix_free (Y);
-		gsl_matrix_free (W);
-		gsl_matrix_free(B);
-		gsl_matrix_free(se_B);
-		gsl_matrix_free (G);	
-		gsl_matrix_free (U);
-		gsl_matrix_free (UtW);
-		gsl_matrix_free (UtY);
-		gsl_vector_free (eval);
-	} 
-	
-	
-	//BSLMM
-	if (cPar.a_mode==11 || cPar.a_mode==12 || cPar.a_mode==13) {
-		gsl_vector *y=gsl_vector_alloc (cPar.ni_test);
-		gsl_matrix *W=gsl_matrix_alloc (y->size, cPar.n_cvt);	
-		gsl_matrix *G=gsl_matrix_alloc (y->size, y->size);
-		gsl_matrix *UtX=gsl_matrix_alloc (y->size, cPar.ns_test);	
-		
-		//set covariates matrix W and phenotype vector y		
-		//an intercept should be included in W, 
-		cPar.CopyCvtPhen (W, y, 0);
-		
-		//center y, even for case/control data
-		cPar.pheno_mean=CenterVector(y);
-
-		//run bslmm if rho==1
-		if (cPar.rho_min==1 && cPar.rho_max==1) {
-		  //read genotypes X (not UtX)
-		  cPar.ReadGenotypes (UtX, G, false);
-
-		  //perform BSLMM analysis
-		  BSLMM cBslmm;
-		  cBslmm.CopyFromParam(cPar);
-		  time_start=clock();	
-		  cBslmm.MCMC(UtX, y);
-		  cPar.time_opt=(clock()-time_start)/(double(CLOCKS_PER_SEC)*60.0);
-		  cBslmm.CopyToParam(cPar);
-		  //else, if rho!=1
-		} else {
-		gsl_matrix *U=gsl_matrix_alloc (y->size, y->size); 
-		gsl_vector *eval=gsl_vector_alloc (y->size);
-		gsl_matrix *UtW=gsl_matrix_alloc (y->size, W->size2);
-		gsl_vector *Uty=gsl_vector_alloc (y->size);
-
-		
-		//read relatedness matrix G		
-		if (!(cPar.file_kin).empty()) {		
-			cPar.ReadGenotypes (UtX, G, false);
-			
-			//read relatedness matrix G
-			ReadFile_kin (cPar.file_kin, cPar.indicator_idv, cPar.mapID2num, cPar.k_mode, cPar.error, G);
-			if (cPar.error==true) {cout<<"error! fail to read kinship/relatedness file. "<<endl; return;}
-			
-			//center matrix G
-			CenterMatrix (G);
-		} else {
-			cPar.ReadGenotypes (UtX, G, true);
-		}
-		
-		//eigen-decomposition and calculate trace_G
-		cout<<"Start Eigen-Decomposition..."<<endl;
-		time_start=clock();
-		cPar.trace_G=EigenDecomp (G, U, eval, 0);
-		cPar.trace_G=0.0;
-		for (size_t i=0; i<eval->size; i++) {
-			if (gsl_vector_get (eval, i)<1e-10) {gsl_vector_set (eval, i, 0);}
-			cPar.trace_G+=gsl_vector_get (eval, i);
-		}
-		cPar.trace_G/=(double)eval->size;
-		cPar.time_eigen=(clock()-time_start)/(double(CLOCKS_PER_SEC)*60.0);			
-		
-		//calculate UtW and Uty		
-		CalcUtX (U, W, UtW);
-		CalcUtX (U, y, Uty);
-		
-		//calculate REMLE/MLE estimate and pve
-		CalcLambda ('L', eval, UtW, Uty, cPar.l_min, cPar.l_max, cPar.n_region, cPar.l_mle_null, cPar.logl_mle_H0);
-		CalcLambda ('R', eval, UtW, Uty, cPar.l_min, cPar.l_max, cPar.n_region, cPar.l_remle_null, cPar.logl_remle_H0);
-		CalcPve (eval, UtW, Uty, cPar.l_remle_null, cPar.trace_G, cPar.pve_null, cPar.pve_se_null);
-		
-		cPar.PrintSummary();
-				
-		//Creat and calcualte UtX, use a large memory
-		cout<<"Calculating UtX..."<<endl;
-		time_start=clock();							
-		CalcUtX (U, UtX);
-		cPar.time_UtX=(clock()-time_start)/(double(CLOCKS_PER_SEC)*60.0);
-		
-		//perform BSLMM analysis
-		BSLMM cBslmm;
-		cBslmm.CopyFromParam(cPar);
-		time_start=clock();	
-		if (cPar.a_mode==12) {  //ridge regression				
-			cBslmm.RidgeR(U, UtX, Uty, eval, cPar.l_remle_null);
-		} else {	//Run MCMC
-			cBslmm.MCMC(U, UtX, Uty, eval, y);
-		}
-		cPar.time_opt=(clock()-time_start)/(double(CLOCKS_PER_SEC)*60.0);
-		cBslmm.CopyToParam(cPar);
-		
-		//release all matrices and vectors
-		gsl_matrix_free (G);	
-		gsl_matrix_free (U);
-		gsl_matrix_free (UtW);
-		gsl_vector_free (eval);
-		gsl_vector_free (Uty);
-
-		}
-		gsl_matrix_free (W);
-		gsl_vector_free (y);
-		gsl_matrix_free (UtX);
-	} 
-	
-	
-		
-	cPar.time_total=(clock()-time_begin)/(double(CLOCKS_PER_SEC)*60.0);
-	
-	return;
-}
-
-
-
-
-void GEMMA::WriteLog (int argc, char ** argv, PARAM &cPar) 
-{
-	string file_str;
-	file_str=cPar.path_out+"/"+cPar.file_out;
-	file_str+=".log.txt";
-	
-	ofstream outfile (file_str.c_str(), ofstream::out);
-	if (!outfile) {cout<<"error writing log file: "<<file_str.c_str()<<endl; return;}
-	
-	outfile<<"##"<<endl;
-	outfile<<"## GEMMA Version = "<<version<<endl;
-	
-	outfile<<"##"<<endl;
-	outfile<<"## Command Line Input = ";
-	for(int i = 1; i < argc; i++) {	
-		outfile<<argv[i]<<" ";
-	}
-	outfile<<endl;
-
-	outfile<<"##"<<endl;
-	time_t  rawtime; 
-	time(&rawtime);
-	tm *ptm = localtime (&rawtime);
-
-	outfile<<"## Date = "<<asctime(ptm)<<endl;
-	  //ptm->tm_year<<":"<<ptm->tm_month<<":"<<ptm->tm_day":"<<ptm->tm_hour<<":"<<ptm->tm_min<<endl;
-	
-	outfile<<"##"<<endl;
-	outfile<<"## Summary Statistics:"<<endl;
-	outfile<<"## number of total individuals = "<<cPar.ni_total<<endl;	
-	if (cPar.a_mode==43) {
-		outfile<<"## number of analyzed individuals = "<<cPar.ni_cvt<<endl;
-		outfile<<"## number of individuals with full phenotypes = "<<cPar.ni_test<<endl;
-	} else {
-		outfile<<"## number of analyzed individuals = "<<cPar.ni_test<<endl;
-	}
-	outfile<<"## number of covariates = "<<cPar.n_cvt<<endl;
-	outfile<<"## number of phenotypes = "<<cPar.n_ph<<endl;
-	if (cPar.a_mode==43) {
-		outfile<<"## number of observed data = "<<cPar.np_obs<<endl;
-		outfile<<"## number of missing data = "<<cPar.np_miss<<endl;
-	}
-	if (cPar.a_mode==61) {
-		outfile<<"## number of variance components = "<<cPar.n_vc<<endl;
-	}
-		
-	if (!(cPar.file_gene).empty()) {
-		outfile<<"## number of total genes = "<<cPar.ng_total<<endl;
-		outfile<<"## number of analyzed genes = "<<cPar.ng_test<<endl;		
-	} else if (cPar.file_epm.empty()) {	
-		outfile<<"## number of total SNPs = "<<cPar.ns_total<<endl;	
-		outfile<<"## number of analyzed SNPs = "<<cPar.ns_test<<endl;
-	} else {
-		outfile<<"## number of analyzed SNPs = "<<cPar.ns_test<<endl;
-	}
-	
-	if (cPar.a_mode==13) {
-		outfile<<"## number of cases = "<<cPar.ni_case<<endl;
-		outfile<<"## number of controls = "<<cPar.ni_control<<endl;
-	}
-
-
-	if (cPar.a_mode==61) {
-	  //	        outfile<<"## REMLE log-likelihood in the null model = "<<cPar.logl_remle_H0<<endl;
-		if (cPar.n_ph==1) {
-		  outfile<<"## pve estimate in the null model = ";
-		  for (size_t i=0; i<cPar.v_pve.size(); i++) {
-		    outfile<<"  "<<cPar.v_pve[i];
-		  }
-		  outfile<<endl;
-
-		  outfile<<"## se(pve) in the null model = ";
-		  for (size_t i=0; i<cPar.v_se_pve.size(); i++) {
-		    outfile<<"  "<<cPar.v_se_pve[i];
-		  }
-		  outfile<<endl;
-
-		  outfile<<"## sigma2 estimate in the null model = ";
-		  for (size_t i=0; i<cPar.v_sigma2.size(); i++) {
-		    outfile<<"  "<<cPar.v_sigma2[i];
-		  }
-		  outfile<<endl;
-
-		  outfile<<"## se(sigma2) in the null model = ";
-		  for (size_t i=0; i<cPar.v_se_sigma2.size(); i++) {
-		    outfile<<"  "<<cPar.v_se_sigma2[i];
-		  }
-		  outfile<<endl;
-		  /*
-			outfile<<"## beta estimate in the null model = ";
-			for (size_t i=0; i<cPar.beta_remle_null.size(); i++) {
-				outfile<<"  "<<cPar.beta_remle_null[i];
-			}
-			outfile<<endl;
-			outfile<<"## se(beta) = ";
-			for (size_t i=0; i<cPar.se_beta_remle_null.size(); i++) {
-				outfile<<"  "<<cPar.se_beta_remle_null[i];
-			}
-			outfile<<endl;
-		  */
-		}
-	}
-	
-	if (cPar.a_mode==1 || cPar.a_mode==2 || cPar.a_mode==3 || cPar.a_mode==4 || cPar.a_mode==5 || cPar.a_mode==11 || cPar.a_mode==12 || cPar.a_mode==13) {
-		outfile<<"## REMLE log-likelihood in the null model = "<<cPar.logl_remle_H0<<endl;
-		outfile<<"## MLE log-likelihood in the null model = "<<cPar.logl_mle_H0<<endl;
-		if (cPar.n_ph==1) {
-			//outfile<<"## lambda REMLE estimate in the null (linear mixed) model = "<<cPar.l_remle_null<<endl;
-			//outfile<<"## lambda MLE estimate in the null (linear mixed) model = "<<cPar.l_mle_null<<endl;	
-			outfile<<"## pve estimate in the null model = "<<cPar.pve_null<<endl;
-			outfile<<"## se(pve) in the null model = "<<cPar.pve_se_null<<endl;	
-			outfile<<"## vg estimate in the null model = "<<cPar.vg_remle_null<<endl;
-			outfile<<"## ve estimate in the null model = "<<cPar.ve_remle_null<<endl;	
-			outfile<<"## beta estimate in the null model = ";
-			for (size_t i=0; i<cPar.beta_remle_null.size(); i++) {
-				outfile<<"  "<<cPar.beta_remle_null[i];
-			}
-			outfile<<endl;
-			outfile<<"## se(beta) = ";
-			for (size_t i=0; i<cPar.se_beta_remle_null.size(); i++) {
-				outfile<<"  "<<cPar.se_beta_remle_null[i];
-			}
-			outfile<<endl;
-			
-		} else {
-			size_t c;
-			outfile<<"## REMLE estimate for Vg in the null model: "<<endl;			
-			for (size_t i=0; i<cPar.n_ph; i++) {
-				for (size_t j=0; j<=i; j++) {
-					c=(2*cPar.n_ph-min(i,j)+1)*min(i,j)/2+max(i,j)-min(i,j);
-					outfile<<cPar.Vg_remle_null[c]<<"\t";
-				}
-				outfile<<endl;
-			}
-			outfile<<"## se(Vg): "<<endl;	
-			for (size_t i=0; i<cPar.n_ph; i++) {
-				for (size_t j=0; j<=i; j++) {
-					c=(2*cPar.n_ph-min(i,j)+1)*min(i,j)/2+max(i,j)-min(i,j);
-					outfile<<sqrt(cPar.VVg_remle_null[c])<<"\t";
-				}
-				outfile<<endl;
-			}
-			outfile<<"## REMLE estimate for Ve in the null model: "<<endl;	
-			for (size_t i=0; i<cPar.n_ph; i++) {
-				for (size_t j=0; j<=i; j++) {
-					c=(2*cPar.n_ph-min(i,j)+1)*min(i,j)/2+max(i,j)-min(i,j);
-					outfile<<cPar.Ve_remle_null[c]<<"\t";
-				}
-				outfile<<endl;
-			}
-			outfile<<"## se(Ve): "<<endl;	
-			for (size_t i=0; i<cPar.n_ph; i++) {
-				for (size_t j=0; j<=i; j++) {
-					c=(2*cPar.n_ph-min(i,j)+1)*min(i,j)/2+max(i,j)-min(i,j);
-					outfile<<sqrt(cPar.VVe_remle_null[c])<<"\t";
-				}
-				outfile<<endl;
-			}
-			
-			outfile<<"## MLE estimate for Vg in the null model: "<<endl;
-			for (size_t i=0; i<cPar.n_ph; i++) {
-				for (size_t j=0; j<cPar.n_ph; j++) {
-					c=(2*cPar.n_ph-min(i,j)+1)*min(i,j)/2+max(i,j)-min(i,j);
-					outfile<<cPar.Vg_mle_null[c]<<"\t";
-				}
-				outfile<<endl;
-			}
-			outfile<<"## se(Vg): "<<endl;	
-			for (size_t i=0; i<cPar.n_ph; i++) {
-				for (size_t j=0; j<=i; j++) {
-					c=(2*cPar.n_ph-min(i,j)+1)*min(i,j)/2+max(i,j)-min(i,j);
-					outfile<<sqrt(cPar.VVg_mle_null[c])<<"\t";
-				}
-				outfile<<endl;
-			}
-			outfile<<"## MLE estimate for Ve in the null model: "<<endl;	
-			for (size_t i=0; i<cPar.n_ph; i++) {
-				for (size_t j=0; j<cPar.n_ph; j++) {
-					c=(2*cPar.n_ph-min(i,j)+1)*min(i,j)/2+max(i,j)-min(i,j);
-					outfile<<cPar.Ve_mle_null[c]<<"\t";
-				}
-				outfile<<endl;
-			}
-			outfile<<"## se(Ve): "<<endl;	
-			for (size_t i=0; i<cPar.n_ph; i++) {
-				for (size_t j=0; j<=i; j++) {
-					c=(2*cPar.n_ph-min(i,j)+1)*min(i,j)/2+max(i,j)-min(i,j);
-					outfile<<sqrt(cPar.VVe_mle_null[c])<<"\t";
-				}
-				outfile<<endl;
-			}
-			outfile<<"## estimate for B (d by c) in the null model (columns correspond to the covariates provided in the file): "<<endl;
-			for (size_t i=0; i<cPar.n_ph; i++) {
-				for (size_t j=0; j<cPar.n_cvt; j++) {
-					c=i*cPar.n_cvt+j;
-					outfile<<cPar.beta_remle_null[c]<<"\t";
-				}
-				outfile<<endl;
-			}
-			outfile<<"## se(B): "<<endl;
-			for (size_t i=0; i<cPar.n_ph; i++) {
-				for (size_t j=0; j<cPar.n_cvt; j++) {
-					c=i*cPar.n_cvt+j;
-					outfile<<cPar.se_beta_remle_null[c]<<"\t";
-				}
-				outfile<<endl;
-			}
-		}
-	}
-	
-	/*
-	if (cPar.a_mode==1 || cPar.a_mode==2 || cPar.a_mode==3 || cPar.a_mode==4 || cPar.a_mode==11 || cPar.a_mode==12 || cPar.a_mode==13) {
-		if (cPar.n_ph==1) {
-			outfile<<"## REMLE vg estimate in the null model = "<<cPar.vg_remle_null<<endl;
-			outfile<<"## REMLE ve estimate in the null model = "<<cPar.ve_remle_null<<endl;	
-		} else {
-			size_t c;
-			outfile<<"## REMLE estimate for Vg in the null model: "<<endl;			
-			for (size_t i=0; i<cPar.n_ph; i++) {
-				for (size_t j=0; j<=i; j++) {
-					c=(2*cPar.n_ph-min(i,j)+1)*min(i,j)/2+max(i,j)-min(i,j);
-					outfile<<cPar.Vg_remle_null[c]<<"\t";
-				}
-				outfile<<endl;
-			}
-			outfile<<"## REMLE estimate for Ve in the null model: "<<endl;	
-			for (size_t i=0; i<cPar.n_ph; i++) {
-				for (size_t j=0; j<=i; j++) {
-					c=(2*cPar.n_ph-min(i,j)+1)*min(i,j)/2+max(i,j)-min(i,j);
-					outfile<<cPar.Ve_remle_null[c]<<"\t";
-				}
-				outfile<<endl;
-			}
-		}
-	}
-	 */
-	
-	
-	if (cPar.a_mode==11 || cPar.a_mode==12 || cPar.a_mode==13) {
-		outfile<<"## estimated mean = "<<cPar.pheno_mean<<endl;
-	}
-	
-	if (cPar.a_mode==11 || cPar.a_mode==13) {	
-		outfile<<"##"<<endl;
-		outfile<<"## MCMC related:"<<endl;	
-		outfile<<"## initial value of h = "<<cPar.cHyp_initial.h<<endl;
-		outfile<<"## initial value of rho = "<<cPar.cHyp_initial.rho<<endl;
-		outfile<<"## initial value of pi = "<<exp(cPar.cHyp_initial.logp)<<endl;
-		outfile<<"## initial value of |gamma| = "<<cPar.cHyp_initial.n_gamma<<endl;
-		outfile<<"## random seed = "<<cPar.randseed<<endl;
-		outfile<<"## acceptance ratio = "<<(double)cPar.n_accept/(double)((cPar.w_step+cPar.s_step)*cPar.n_mh)<<endl;
-	}
-	
-	outfile<<"##"<<endl;
-	outfile<<"## Computation Time:"<<endl;
-	outfile<<"## total computation time = "<<cPar.time_total<<" min "<<endl;
-	outfile<<"## computation time break down: "<<endl;
-	if (cPar.a_mode==21 || cPar.a_mode==22 || cPar.a_mode==11 || cPar.a_mode==13) {
-		outfile<<"##      time on calculating relatedness matrix = "<<cPar.time_G<<" min "<<endl;
-	}
-	if (cPar.a_mode==31) {
-		outfile<<"##      time on eigen-decomposition = "<<cPar.time_eigen<<" min "<<endl;
-	}
-	if (cPar.a_mode==1 || cPar.a_mode==2 || cPar.a_mode==3 || cPar.a_mode==4 || cPar.a_mode==5 || cPar.a_mode==11 || cPar.a_mode==12 || cPar.a_mode==13) {
-		outfile<<"##      time on eigen-decomposition = "<<cPar.time_eigen<<" min "<<endl;
-		outfile<<"##      time on calculating UtX = "<<cPar.time_UtX<<" min "<<endl;		
-	}
-	if ((cPar.a_mode>=1 && cPar.a_mode<=4) || (cPar.a_mode>=51 && cPar.a_mode<=54) ) {
-		outfile<<"##      time on optimization = "<<cPar.time_opt<<" min "<<endl;
-	}
-	if (cPar.a_mode==11 || cPar.a_mode==13) {
-		outfile<<"##      time on proposal = "<<cPar.time_Proposal<<" min "<<endl;
-		outfile<<"##      time on mcmc = "<<cPar.time_opt<<" min "<<endl;
-		outfile<<"##      time on Omega = "<<cPar.time_Omega<<" min "<<endl;
-	}
-	if (cPar.a_mode==41 || cPar.a_mode==42) {
-		outfile<<"##      time on eigen-decomposition = "<<cPar.time_eigen<<" min "<<endl;
-	}
-	if (cPar.a_mode==43) {
-		outfile<<"##      time on eigen-decomposition = "<<cPar.time_eigen<<" min "<<endl;
-		outfile<<"##      time on predicting phenotypes = "<<cPar.time_opt<<" min "<<endl;
-	}
-	outfile<<"##"<<endl;
-	
-	outfile.close();
-	outfile.clear();
-	return;
-}
-
-
-- 
cgit v1.2.3