/*
 * treetest.c
 *
 *
 * Part of TREE-PUZZLE 5.2 (July 2004)
 *
 * (c) 2003-2004 by Heiko A. Schmidt, Korbinian Strimmer, and Arndt von Haeseler
 * (c) 1999-2003 by Heiko A. Schmidt, Korbinian Strimmer,
 *                  M. Vingron, and Arndt von Haeseler
 * (c) 1995-1999 by Korbinian Strimmer and Arndt von Haeseler
 *
 * All parts of the source except where indicated are distributed under
 * the GNU public licence.  See http://www.opensource.org for details.
 *
 * ($Id$)
 *
 */

#include <treetest.h>

#ifndef STDOUT
#	define STDOUT stdout
#endif

/*****************************************************************************/
/* ELW/SR Test (Expected Likelihood Weights, Strimmer & Rambaut, 2002)       */
/*****************************************************************************/
void elw_test(ivector  Alias,
              dmatrix  allsites,
              int      numutrees,
              int      numsites,
              double   siglevel,
              int      numboots,
              ivector *elw_test_passed,
              dvector *elw_Support)
{

	/* 
	 * Determine posterior probabilties and support values
	 * for each hypothesis and store results in public arrays
	 * posterior, support etc which will automatically be
	 * created by this procedure.
	 *
	 * allsites    log-likelihoods of each pattern
	 * Alias       map of patterns to sites in sequence
	 * numboots    number of bootstraps
	 */	

	dvector deltaL;          /* delta of likelihoods (double[numH]) */
	dvector elw_support;     /* (double[numH]) */
	dvector posterior;       /* (double[numH]) */
	ivector likelihoodOrder; /* (int[numH]) */
	ivector supportOrder;    /* (int[numH]) */
	dvector rs;              /* (int[numH]) */

	double sum1;
	double sum;
	int numH;		/* number of hypotheses */
	int numSites;		/* number of sites */
	int s, p;
	int i, j, k;
	/* int j; */
	double maxL;
	double maxLogL;
	ivector sr_test_res; /* [numutrees]; */
	int elw_numboots = 1000;


	/* number of hypothesis  */
	numH = numutrees;
	
	/* allocate public arrays */
	deltaL = new_dvector(numH);
	elw_support = new_dvector(numH);
	posterior = new_dvector(numH);
	likelihoodOrder = new_ivector(numH);
	supportOrder = new_ivector(numH);	
	sr_test_res = new_ivector(numutrees);

	*elw_test_passed = sr_test_res;
	*elw_Support = elw_support;

	/* number of sites */
	numSites = numsites;
		

	/******************************************/
	/* Compute log-likelihoods, their order,  */
	/* their deltas and their posteriors      */
	/******************************************/

	/* initialize delta vector for each tree */
	for (j = 0; j < numSites; j++) {
		if (Alias == NULL) {
			p = j;
		} else {
			p = Alias[j];
		}

		for (k = 0; k < numH; k++) {
			deltaL[k] -= allsites[k][p];
		}
	}

	/* sort likelihoods -> sorted index list likelihoodOrder */
	/* insertion sort */
	{
		int tmp;
		int idx;
		for (k=0; k < numH; k++) likelihoodOrder[k] = k;
		for (j=0; j < numH-1; j++) {
			idx=j;
			for (k=j+1; k < numH; k++) {
				if (deltaL[likelihoodOrder[k]] < deltaL[likelihoodOrder[idx]]) idx=k;
			}
			tmp = likelihoodOrder[j];
			likelihoodOrder[j] = likelihoodOrder[idx];
			likelihoodOrder[idx] = tmp;
		}
	}


	/* Compute deltas */
	maxL= -deltaL[likelihoodOrder[0]];
	for (j = 0; j < numH; j++) {
		deltaL[j] = -(deltaL[j]+maxL);
	}

	/* compute posterior probabilities  */
	sum1 = 0.0;
	for (j = 0; j < numH; j++) {
		posterior[j] = exp(deltaL[j]);
		sum1 += posterior[j];
	}
	for (j = 0; j < numH; j++) {
		posterior[j] = posterior[j]/sum1;
	}

	/* reverse sign of delta L */
	for (j = 0; j < numH; j++) {
		deltaL[j] = -deltaL[j];
	}
	deltaL[likelihoodOrder[0]] = 0.0;


	/* Bootstrap/Resample data */
	/***************************/

	/* temporary memory  */
	rs = new_dvector(numH);
		
	/* MersenneTwisterFast mt = new MersenneTwisterFast(); */
	for (i = 0; i < elw_numboots; i++) {		/* num bootstraps */
		for (k = 0; k < numH; k++) {
			rs[k] = 0;
		}

		for (j = 0; j < numSites; j++) {	/* bootstrapping sites for ...*/
			/* int s = mt.nextInt(numSites); */
			s = randominteger(numSites);

			if (Alias == NULL) {
				p = s;
			} else {
				p = Alias[s];
			}

			for (k = 0; k < numH; k++) {	/* ...each utree <k> */
				/* rs[k] += pLogL[k][p]; */
				rs[k] += allsites[k][p];
			}
		}
			
		/* find ml hypothesis */
		/* maxLogL = findMax(rs); */
		{
			int best = 0; /* can be removed (HAS) */
			int idx;
			double max = rs[0];
			for (idx = 1; idx < numH; idx++) {
				if (rs[idx] > max) {
					best = idx; /* can be removed (HAS) */
					max = rs[idx];
				}
			}
			maxLogL = max;
		}
			
		/* compute log-likelihood difference */
		for (k = 0; k < numH; k++) {
			rs[k] = rs[k] - maxLogL;
		}

		/* compute posteriors and sum over resampled data set */
		sum = 0.0;
		for (k = 0; k < numH; k++) {
			rs[k] = exp(rs[k]);
			sum += rs[k];
		}
		for (k = 0; k < numH; k++) {
			elw_support[k] += rs[k]/sum;
		}
	}  
		
	/* compute support values */
	for (j = 0; j < numH; j++) {
		elw_support[j] = elw_support[j]/elw_numboots;
	}

	/* determine order of elw_support (smallest->largest) */
	/* HeapSort.sort(elw_support, supportOrder); */
	{
		int tmp;
		int idx;
		for (k=0; k < numH; k++) supportOrder[k] = k;
		for (j=0; j < numH-1; j++) {
			idx=j;
			for (k=j+1; k < numH; k++) {
				if (elw_support[supportOrder[k]] > elw_support[supportOrder[idx]]) idx=k;
			}
			tmp = supportOrder[j];
			supportOrder[j] = supportOrder[idx];
			supportOrder[idx] = tmp;
		}
	}

	sum = 0.0;
	for (k=0; k<numH; k++) {
		if(sum <= .95) {
			sum += elw_support[supportOrder[k]];
			sr_test_res[supportOrder[k]] = 1;     /* within confidence set */
		} else {
			/* sum += elw_support[supportOrder[k]]; */	/* shouldn't the sum increase after exclusion ??? (HAS) */
			sum += elw_support[supportOrder[k]];		/* included to code */
			sr_test_res[supportOrder[k]] = 0;     /* excluded from confidence set */
		}
	}

} /* elw_test */




/* TODO: move to treetest.c */

/*****************************************************************************/
/* SH Test (Shimodaira & Hasegawa, 1999)                                     */
/*****************************************************************************/
void sh_test(ivector  Alias,		/* site pattern translat array */
             dmatrix  allsites,		/* tree site log-likelihoods */
             int      numutrees,	/* number of trees */
             int      numsites,		/* number of sites */
             double   siglevel,		/* significance level to test */
             int      numboots,		/* number of bootstraps to draw */
             ivector *sh_test_passed,	/* has tree passed tests */
             dvector *sh_Pval)		/* and p-value */
{

	/*
 	 * Shimodaira-Hasegawa-Test (1999) to
 	 * compare a set of evolutionary hypotheses.
	 * Compare all given hypotheses to the best (ML) hypothesis
	 * and store results in public arrays delta, pval
	 * (which will automatically be created by this procedure).
	 *
	 * allsites    log-likelihoods of each pattern
	 * Alias       map of patterns to sites in sequence
	 * numboots    number of bootstraps
	 */	
		
	int bestH; /* number of maximum likelihood hypothesis */
	dvector sh_delta; /* log-LH difference to ml hypothesis */
	dvector sh_pval; /* corresponding p-value */
	dvector sh_logL;
	int sh_numboots = 1000;
	double maxLogL;
	double m;
	dmatrix bs;             /* = new int[numH,sh_numboots]; */
	int best;
	double colmax;
	int count;
	double sh_prob = .05;
	ivector sh_test_res; /* [numutrees]; */
	int numH;
	int numSites;
	int i, j, k;
	int p, s;


	/* number of hypothesis */
	numH = numutrees;

	/* allocate memory for results */
	sh_delta    = new_dvector(numH);
	sh_pval     = new_dvector(numH);
	sh_test_res = new_ivector(numutrees);

	*sh_Pval = sh_pval;
	*sh_test_passed = sh_test_res;

	/* number of sites */
	numSites = numsites;

	/* log likelihood of each hypothesis */
	sh_logL = new_dvector(numH);
	for (i = 0; i < numSites; i++) {
		if (Alias == NULL) {
			p = i;
		} else {
			p = Alias[i];
		}

		for (j = 0; j < numH; j++) {
			/* sh_logL[j] += pLogL[j][p]; */
			sh_logL[j] += allsites[j][p];
		}
	}

	/* find maximum-likelihood hypothesis */
	bestH = 0;
	maxLogL = sh_logL[0];
	for (i = 1; i < numH; i++) {
		if (sh_logL[i] > maxLogL) {
			bestH = i;
			maxLogL = sh_logL[i];
		}
	}

	/* compute log-likelihood differences to best hypothesis */
	for (i = 0; i < numH; i++) {
		sh_delta[i] = sh_logL[bestH]-sh_logL[i];
	}

	/* allocate temporary memory for resampling procedure */
	bs = new_dmatrix(numH,sh_numboots);

	/* Resample data */
	/* MersenneTwisterFast mt = new MersenneTwisterFast(); */
	for (i = 0; i < sh_numboots; i++) {
		for (j = 0; j < numSites; j++) {
			/* int s = mt.nextInt(numSites); */
			s = randominteger(numSites);

			if (Alias == NULL) {
				p = s;
			} else {
				p = Alias[s];
			}

			for (k = 0; k < numH; k++) {
				/* rs[k][i] += pLogL[k][p]; */
				bs[k][i] += allsites[k][p];
			}
		}
	}  

	/* center resampled log-likelihoods */
	for (i = 0; i < numH; i++) {

		/* double m = DiscreteStatistics.mean(rs[i]); */
		m = 0.0;
		for (j = 0; j < sh_numboots; j++) {
			m += bs[i][j];
		}
		m /= sh_numboots;


		for (j = 0; j < sh_numboots; j++) {
			bs[i][j] = bs[i][j] - m;
		}
	}

	/* compute resampled log-likelihood differences */
	for (i = 0; i < sh_numboots; i++) {
		/* double max = findMaxInColumn(rs, i); */
		{
			best = 0;
			colmax = bs[0][i];
			for (j = 1; j < numH; j++) {
				if (bs[j][i] > colmax) {
					best = j;
					colmax = bs[j][i];
				}
			}
		}

		for (j = 0; j < numH; j++) {
			bs[j][i] = colmax - bs[j][i];
		}
	}

	/* compute p-values for each hypothesis */
	for (i = 0; i < numH; i++) {
		count = 0;
		for (j = 0; j < sh_numboots; j++) {
			if (bs[i][j] >= sh_delta[i]) {
				count++;
			}
		}

		sh_pval[i] = (double) count/(double) sh_numboots;
		if (sh_pval[i] >= sh_prob) {
			sh_test_res[i] = 1;
		} else {
			sh_test_res[i] = 0;
		}
	}

	/* free memory */
	free_dmatrix(bs);
	free_dvector(sh_logL);

} /* sh_test */




/* TODO: move to treetest.c */

/* print tree statistics */
void printtreestats(FILE *ofp,
			dvector ulkl,
			dvector ulklc,
			int numutrees,
			dmatrix  allsites,	/* tree site log-likelihoods */
			dmatrix  allsitesc,	/* tree site log-likelihoods */
			ivector  Alias,		/* link site -> corresponding site pattern */
			int Maxsite,
			int Numptrn,
			ivector Weight,
			int compclock)
{
	int i, j;		/* counter variables */
	int besttree;		/* best tree */

	/* for the KH Test (Kishino & Hasegawa, 1989) */
	double bestlkl;		/* best likelihood */
	double difflkl;		/* best likelihood difference */
	double difflklps;	/* best likelihood difference per site */
	double temp;
	double sum;

	/* variables for the SH Test (Shimodaira & Hasegawa, 1999) */
	dvector sh_pval;		/* corresponding p-value */
	double sh_prob = .05;		/* significance level */
	int sh_numboots = 1000;		/* number of bootstrap samples */
	ivector sh_test_res;		/* result vector  [numutrees]; */

	/* variables for the ELW/SR Test (Expected Likelihood Weights) */
	/* (Strimmer & Rambaut, 2002) */

	dvector elw_support;         /* = new double[numH; */
	ivector elw_test_res;		/* result vector  [numutrees]; */
	int elw_numboots = 1000;	/* number of bootstrap samples */
	double elw_prob = .05;		/* significance level */

	/* variables for the one-sided KH Test using SH */
	dvector kh1_pval_tmp;		/* temp p-value arry for pairwise SH */
	ivector kh1_test_res_tmp;	/* temp result arry for pairwise SH */
	dvector kh1_pval;		/* p-values */
	ivector kh1_test_res;		/* result vector  [numutrees]; */
	dvector kh1_allsites_tmp[2];	/* temp arry for pairwise SH test */
	int kh1_numboots = 1000;	/* number of bootstrap samples */
	double kh1_prob = .05;		/* significance level */

	/* find best tree */
	besttree = 0;
	bestlkl = ulkl[0];
	for (i = 1; i < numutrees; i++) {
		if (ulkl[i] > bestlkl) {
			besttree = i;
			bestlkl = ulkl[i];
		}
	}

	/* one sided KH  = pairwise SH test between tree and besttree */
	fprintf(STDOUT, "Performing single sided KH test.\n");
	fflush(STDOUT);
	kh1_pval     = new_dvector(numutrees);
	kh1_test_res = new_ivector(numutrees);

	kh1_allsites_tmp[0] = allsites[besttree]; /* set best tree */
	for (i = 0; i < numutrees; i++) {
		if (i == besttree) { /* if best tree -> no test */
			kh1_pval[i] = 1.0;
			kh1_test_res[i] = 1;
		} else { /* other wise test pairwise SH */
			if ((ulkl[besttree]-ulkl[i]) < MINKHDIFF) {
				/* if approx equal to best tree -> no test */
				kh1_pval[i] = 1.0;
				kh1_test_res[i] = 1;
			} else { /* other wise test pairwise SH */
				kh1_allsites_tmp[1] = allsites[i]; /* set site log-lh */
				sh_test(Alias, kh1_allsites_tmp, 2, Maxsite, kh1_prob, kh1_numboots, &kh1_test_res_tmp, &kh1_pval_tmp); /* pairwise SH */
				kh1_pval[i] = kh1_pval_tmp[1]; /* store p-value */
				kh1_test_res[i] = kh1_test_res_tmp[1]; /* save result */
			}
		}
	}
	free_ivector(kh1_test_res_tmp); /* moved out of loop */
	free_dvector(kh1_pval_tmp);

	/* ELW */
	fprintf(STDOUT, "Performing ELW test.\n");
	fflush(STDOUT);
	elw_test(Alias, allsites, numutrees, Maxsite, elw_prob, elw_numboots, &elw_test_res, &elw_support);

	/* SH */
	fprintf(STDOUT, "Performing SH test.\n");
	fflush(STDOUT);
	sh_test(Alias, allsites, numutrees, Maxsite, sh_prob, sh_numboots, &sh_test_res, &sh_pval);
	
/*****************************************************************************/
/* two-sided KH Test (Kishino & Hasegawa, 1989)                              */
/* and output                                                                */
/*****************************************************************************/
	
	fprintf(ofp, "\n\nCOMPARISON OF USER TREES (NO CLOCK)\n\n");
#	ifdef KHTWOSIDED
		fprintf(ofp, "Tree   log L   difference    S.E.   Sig. worse   p-1sKH     p-SH       c-ELW   \n");
		fprintf(ofp, "-------------------------------------------------------------------------------\n");
#	else
			fprintf(ofp, "Tree   log L   difference    S.E.      p-1sKH     p-SH       c-ELW      2sKH\n");
			fprintf(ofp, "-------------------------------------------------------------------------------\n");
#	endif
	for (i = 0; i < numutrees; i++) {
		difflkl = ulkl[besttree]-ulkl[i];
		if (difflkl < MINKHDIFF) { 
			/* if lkls (approx) equal to best tree, no difference */
			difflkl = 0.0;
			fprintf(ofp, "%2d %10.2f %8.2f ", i+1, ulkl[i], difflkl);
		} else {
			fprintf(ofp, "%2d %10.2f %8.2f ", i+1, ulkl[i], difflkl);
		}
		/* fprintf(ofp, "%2d %10.2f %8.2f ", i+1, ulkl[i], difflkl); */
		if (i == besttree) {
#			ifdef KHTWOSIDED
				fprintf(ofp, " <-------------- best ");
#			else
				fprintf(ofp, " <---- best ");
				/* fprintf(ofp, " <---- best "); */
#			endif

			if (kh1_test_res[i] == 1)
				fprintf(ofp, "   %6.4f +", kh1_pval[i]);
			else
				fprintf(ofp, "   %6.4f -", kh1_pval[i]);

			if (sh_test_res[i] == 1)
				fprintf(ofp, "   %6.4f +", sh_pval[i]);
			else
				fprintf(ofp, "   %6.4f -", sh_pval[i]);

			if (elw_test_res[i] == 1)
				fprintf(ofp, "   %6.4f +", elw_support[i]);
			else
				fprintf(ofp, "   %6.4f -", elw_support[i]);

#			ifdef KHTWOSIDED
#			else
					fprintf(ofp, "    best");
#			endif
		} else {
			/* compute variance of Log L differences over sites */
#ifndef USE_WINDOWS
			difflklps = difflkl/(double)Maxsite;
#else
			difflklps = difflkl/(double)alimaxsite;
#endif
			sum = 0.0;
			for (j = 0; j < Numptrn; j++) {
				temp = allsites[besttree][j] - allsites[i][j] - difflklps;
				sum += temp*temp*Weight[j];
			}
#ifndef USE_WINDOWS
			sum = sqrt( fabs(sum/(Maxsite-1.0)*Maxsite) );
#else
			sum = sqrt(fabs(sum/(alimaxsite-1.0)*alimaxsite));
#endif
#			ifdef KHTWOSIDED
			fprintf(ofp, "%11.2f      ", sum);
				if (difflkl > 1.96*sum)
					fprintf(ofp, "yes  ");
				else
					fprintf(ofp, "no   ");
#			else
				fprintf(ofp, "%11.4f ", sum);
				/* fprintf(ofp, "%11.2f ", sum); */
#			endif

			if (kh1_test_res[i] == 1)
				fprintf(ofp, "   %6.4f +", kh1_pval[i]);
			else
				fprintf(ofp, "   %6.4f -", kh1_pval[i]);
			if (sh_test_res[i] == 1)
				fprintf(ofp, "   %6.4f +", sh_pval[i]);
			else
				fprintf(ofp, "   %6.4f -", sh_pval[i]);

			if (elw_test_res[i] == 1)
				fprintf(ofp, "   %6.4f +", elw_support[i]);
			else
				fprintf(ofp, "   %6.4f -", elw_support[i]);

#			ifdef KHTWOSIDED
#			else
#if 0
				if (difflkl > 1.96*sum)
					fprintf(ofp, "       - (diff=%.4f 1.96*sum=%.4e)", difflkl, 1.96*sum);
				else
					fprintf(ofp, "       + (diff=%.4f 1.96*sum=%.4e)", difflkl, 1.96*sum);
#endif
				if (difflkl > 1.96*sum)
					fprintf(ofp, "       -");
				else
					fprintf(ofp, "       +");
#			endif
		}
		fprintf(ofp, "\n");
	}
	
	fprintf(ofp, "\nThe columns show the results and p-values of the following tests:\n");
	fprintf(ofp, "1sKH - one sided KH test based on pairwise SH tests (Shimodaira-Hasegawa\n");
	fprintf(ofp, "       2000, Goldman et al., 2001, Kishino-Hasegawa 1989)\n");
	fprintf(ofp, "SH   - Shimodaira-Hasegawa test (2000)\n");
	fprintf(ofp, "ELW  - Expected Likelihood Weight (Strimmer-Rambaut 2002)\n");
	fprintf(ofp, "2sKH - two sided Kishino-Hasegawa test (1989)\n");
	fprintf(ofp, "\n");
	fprintf(ofp, "Plus signs denote the confidence sets. Minus signs denote significant\n");
	fprintf(ofp, "exclusion. All tests used 5%% significance level. 1sKH, SH, and ELW\n");
	fprintf(ofp, "performed 1000 resamplings using the RELL method.\n");
	fprintf(ofp, "1sKH and 2sKH are correct to the 2nd position after the the decimal\n");
	fprintf(ofp, "point of the log-likelihoods.\n");
	
	if (compclock) {
	
		/* find best tree */
		besttree = 0;
		bestlkl = ulklc[0];
		for (i = 1; i < numutrees; i++)
			if (ulklc[i] > bestlkl) {
				besttree = i;
				bestlkl = ulklc[i];
			}
	
		/* one sided KH */
		fprintf(STDOUT, "Performing single sided KH test (clock).\n");
		fflush(STDOUT);
		kh1_pval     = new_dvector(numutrees);
		kh1_test_res = new_ivector(numutrees);
	
		kh1_allsites_tmp[0] = allsitesc[besttree];
		for (i = 0; i < numutrees; i++) {
			if (i != besttree) {
				if ((ulklc[besttree]-ulklc[i]) < MINKHDIFF) {
					/* if approx equal to best tree -> no test */
					kh1_pval[i] = 1.0;
					kh1_test_res[i] = 1;
				} else { /* other wise test pairwise SH */
					kh1_allsites_tmp[1] = allsitesc[i];
					sh_test(Alias, kh1_allsites_tmp, 2, Maxsite, kh1_prob, kh1_numboots, &kh1_test_res_tmp, &kh1_pval_tmp);
					kh1_pval[i] = kh1_pval_tmp[1];
					kh1_test_res[i] = kh1_test_res_tmp[1];
					free_ivector(kh1_test_res_tmp);
					free_dvector(kh1_pval_tmp);
				}
			} else {
				kh1_pval[i] = 1.0;
				kh1_test_res[i] = 1;
			}
		}

		/* ELW */
		fprintf(STDOUT, "Performing ELW test (clock).\n");
		fflush(STDOUT);
		elw_test(Alias, allsitesc, numutrees, Maxsite, elw_prob, elw_numboots, &elw_test_res, &elw_support);

		/* SH */
		fprintf(STDOUT, "Performing SH test (clock).\n");
		fflush(STDOUT);
		sh_test(Alias, allsitesc, numutrees, Maxsite, sh_prob, sh_numboots, &sh_test_res, &sh_pval);

	
		fprintf(ofp, "\n\nCOMPARISON OF USER TREES (WITH CLOCK)\n\n");
#		ifdef KHTWOSIDED
			fprintf(ofp, "Tree   log L   difference    S.E.   Sig. worse   p-1sKH     p-SH       c-ELW   \n");
			fprintf(ofp, "-------------------------------------------------------------------------------\n");
#		else
			fprintf(ofp, "Tree   log L   difference    S.E.      p-1sKH     p-SH       c-ELW      2sKH\n");
			fprintf(ofp, "-------------------------------------------------------------------------------\n");
#		endif
		for (i = 0; i < numutrees; i++) {
			difflkl = ulklc[besttree]-ulklc[i];
			if (difflkl < MINKHDIFF) { 
				/* if lkls (approx) equal to best tree, no difference */
				difflkl = 0.0;
				fprintf(ofp, "%2d %10.2f %8.2f ", i+1, ulklc[i], difflkl);
			} else {
				fprintf(ofp, "%2d %10.2f %8.2f ", i+1, ulklc[i], difflkl);
			}
			if (i == besttree) {
#		ifdef KHTWOSIDED
				fprintf(ofp, " <-------------- best ");
#		else
				fprintf(ofp, " <---- best ");
#		endif

				if (kh1_test_res[i] == 1)
					fprintf(ofp, "   %6.4f +", kh1_pval[i]);
				else
					fprintf(ofp, "   %6.4f -", kh1_pval[i]);
				if (sh_test_res[i] == 1)
					fprintf(ofp, "   %6.4f +", sh_pval[i]);
				else
					fprintf(ofp, "   %6.4f -", sh_pval[i]);
	
				if (elw_test_res[i] == 1)
					fprintf(ofp, "   %6.4f +", elw_support[i]);
				else
					fprintf(ofp, "   %6.4f -", elw_support[i]);

#				ifdef KHTWOSIDED
#				else
						fprintf(ofp, "    best");
#				endif
			} else {
				/* compute variance of Log L differences over sites */
#ifndef USE_WINDOWS
				difflklps = difflkl/(double)Maxsite;
#else
				difflklps = difflkl/(double)alimaxsite;
#endif
				sum = 0.0;
				for (j = 0; j < Numptrn; j++) {
					temp = allsitesc[besttree][j] - allsitesc[i][j] - difflklps;
					sum += temp*temp*Weight[j];
				}
#ifndef USE_WINDOWS
				sum = sqrt(fabs(sum/(Maxsite-1.0)*Maxsite));
#else
				sum = sqrt(fabs(sum/(alimaxsite-1.0)*alimaxsite));
#endif
#				ifdef KHTWOSIDED
					fprintf(ofp, "%11.2f      ", sum);
					if (difflkl > 1.96*sum)
						fprintf(ofp, "yes  ");
					else
						fprintf(ofp, "no   ");
#				else
					fprintf(ofp, "%11.2f ", sum);
#				endif

				if (kh1_test_res[i] == 1)
					fprintf(ofp, "   %6.4f +", kh1_pval[i]);
				else
					fprintf(ofp, "   %6.4f -", kh1_pval[i]);
				if (sh_test_res[i] == 1)
					fprintf(ofp, "   %6.4f +", sh_pval[i]);
				else
					fprintf(ofp, "   %6.4f -", sh_pval[i]);
	
				if (elw_test_res[i] == 1)
					fprintf(ofp, "   %6.4f +", elw_support[i]);
				else
					fprintf(ofp, "   %6.4f -", elw_support[i]);

#				ifdef KHTWOSIDED
#				else
					if (difflkl > 1.96*sum)
						fprintf(ofp, "       -");
					else
						fprintf(ofp, "       +");
#				endif
			}
			fprintf(ofp, "\n");
		}
		fprintf(ofp, "\nThe columns show the results and p-values of the following tests:\n");
		fprintf(ofp, "1sKH - one sided KH test based on pairwise SH tests (Shimodaira-Hasegawa\n");
		fprintf(ofp, "       2000, Goldman et al., 2001, Kishino-Hasegawa 1989)\n");
		fprintf(ofp, "SH   - Shimodaira-Hasegawa test (2000)\n");
		fprintf(ofp, "ELW  - Expected Likelihood Weight (Strimmer-Rambaut 2002)\n");
		fprintf(ofp, "2sKH - two sided Kishino-Hasegawa test (1989)\n");
		fprintf(ofp, "\n");
		fprintf(ofp, "Plus signs denote the confidence sets. Minus signs denote significant\n");
		fprintf(ofp, "exclusion. All tests used 5%% significance level. 1sKH, SH, and ELW\n");
		fprintf(ofp, "performed 1000 resamplings using the RELL method.\n");
		fprintf(ofp, "1sKH and 2sKH are correct to the 2nd position after the the decimal\n");
		fprintf(ofp, "point of the log-likelihoods.\n");
#if 0
		fprintf(ofp, "\nColmn 5 gives the results of the (old) two-sided Kishino-Hasegawa test \n");
		fprintf(ofp, "following Kishino and Hasegawa (1989). It tests whether a likelihood is \n");
		fprintf(ofp, "significantly worse than the best one, marked with 'best'. This test should \n");
		fprintf(ofp, "only be used for data not having been determined with the data tested on.\n");
		fprintf(ofp, "1sKH is the one-sided KH test (Goldman et al., 2001). It is applicable\n");
		fprintf(ofp, "to test whether likelihoods are worse than the lielihood of the ML tree.\n");
		fprintf(ofp, "SH tests for the best trees (Shimodaira and Hasegawa, 2000).\n");
		fprintf(ofp, "Note that KH, 1sKH, and SH assume to have the 'true' topologies among the tested.\n");
		fprintf(ofp, "ELW (Expected likelihood weights) seems to work even without this restriction\n");
		fprintf(ofp, "(Strimmer and Rambaut, 2002). Still plausible trees should be among the tested.\n");
		fprintf(ofp, "For 1sKH, SH, and ELW plus signs '+' mark the topologies belonging to the\n");
		fprintf(ofp, "confidence sets the numbers give the p-values (1sKH and SH) or the confidence\n");
		fprintf(ofp, "weight (ELW).\n");
		fprintf(ofp, "All tests used 5%% significance level. 1sKH, SH, and ELW used 1000 resamplings.\n");
#endif
	}
} /* printtreestats */



/********************************************************/