/*
 * puzzle2.c
 *
 *
 * Part of TREE-PUZZLE 5.2 (February 2005)
 *
 * (c) 2003-2005 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$)
 *
 */

#ifdef HAVE_CONFIG_H
#  include <config.h>
#endif


#define EXTERN extern

#include "puzzle.h"
#include "ml.h"
#include <string.h>

#if PARALLEL
#	include "ppuzzle.h"
#	include "sched.h"
#endif /* PARALLEL */

/******************************************************************************/
/*** for debugging ***/
/******************************************************************************/

#define BININPUTDEBUG 1
#undef BININPUTDEBUG

int PPP2=0;
/* fprintf(stderr, "PPP2: %d (%s:%d)\n", PPP2++, __FILE__, __LINE__); */

/* fputid(stderr, 65); fprintf(stderr, " = 65 (%s:%d - 666)\n", __FILE__, __LINE__); */

void fprintfdbl2hex2(FILE *fp, double *ddd)
{
fprintf(fp,"%-16.9f= %02X%02x %02X%02x %02X%02x %02X%02x", *ddd,
	((unsigned char*)(ddd))[0],  ((unsigned char*)(ddd))[1],  
	((unsigned char*)(ddd))[2],  ((unsigned char*)(ddd))[3],  
	((unsigned char*)(ddd))[4],  ((unsigned char*)(ddd))[5],  
	((unsigned char*)(ddd))[6],  ((unsigned char*)(ddd))[7]);
}

void fprintfdbl2hex(FILE *fp, double ddd)
{
fprintf(fp,"%-16.9f= %02X%02x %02X%02x %02X%02x %02X%02x", ddd,
	((unsigned char*)(&ddd))[0],  ((unsigned char*)(&ddd))[1],  
	((unsigned char*)(&ddd))[2],  ((unsigned char*)(&ddd))[3],  
	((unsigned char*)(&ddd))[4],  ((unsigned char*)(&ddd))[5],  
	((unsigned char*)(&ddd))[6],  ((unsigned char*)(&ddd))[7]);
}

/******************************************************************************/
/* sequences                                                                  */
/******************************************************************************/

/* read ten characters of current line as identifier */
void readid(FILE *infp,        /* in:  file pointer to read from */
            int t,             /* in:  current taxon number */
            cmatrix  identif)  /* io:  taxon names (10 char w/o stop) */
{
	int i, j, flag, ci;

	for (i = 0; i < 10; i++) {
		ci = fgetc(infp);
		if (ci == EOF || !isprint(ci)) {
			fprintf(STDOUT, "\n\n\nUnable to proceed (no name for sequence %d)\n\n\n", t+1);
#			if PARALLEL
				PP_Finalize();
#			endif
   			tp_exit(1, NULL, FALSE, __FILE__, __LINE__, exit_wait_optn);
		}
		identif[t][i] = (char) ci;
	}	
	/* convert leading blanks in taxon name to underscores */
	flag = FALSE;
	for (i = 9; i > -1; i--) {
		if (flag == FALSE) {
			if (identif[t][i] != ' ') flag = TRUE; 
		} else {
			if (identif[t][i] == ' ') identif[t][i] = '_';
		}
	}
	/* check whether this name is already used */
	for (i = 0; i < t; i++) { /* compare with all other taxa */
		flag = TRUE; /* assume identity */
		for (j = 0; (j < 10) && (flag == TRUE); j++)
			if (identif[t][j] != identif[i][j])
				flag = FALSE;
		if (flag) {
			fprintf(STDOUT, "\n\n\nUnable to proceed (multiple occurrence of sequence name '");
			fputid(STDOUT, t);
			fprintf(STDOUT, "')\n\n\n");
#			if PARALLEL
				PP_Finalize();
#			endif
   			tp_exit(1, NULL, FALSE, __FILE__, __LINE__, exit_wait_optn);
		}
	}
} /* readid */

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


/* read next allowed character */
char readnextcharacter(FILE *ifp, int notu, int nsite)
{
	char c;

	/* ignore blanks and control characters except newline (UNIX,DOS) or CR (Mac,DOS) */
	do {
		if (fscanf(ifp, "%c", &c) != 1) {
			fprintf(STDOUT, "\n\n\nUnable to proceed (missing character at position %d in sequence '", nsite + 1);
			fputid(STDOUT, notu);
			fprintf(STDOUT, "')\n\n\n");
#			if PARALLEL
				PP_Finalize();
#			endif
   			tp_exit(1, NULL, FALSE, __FILE__, __LINE__, exit_wait_optn);
		}
	} while (c == ' ' || (iscntrl((int) c) && (c != '\n') && (c != '\r')));
	return c;
} /* readnextcharacter */

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


/* skip rest of the line */
void skiprestofline(FILE *ifp,    /* input file stream */
                    int   notu,   /* taxon number      - for error msg. */
                    int   nsite)  /* sequence position - for error msg. */
{
	int ci;

	/* read chars until the first newline or CR */
	do{
		ci = fgetc(ifp);
		if (ci == EOF) {
			fprintf(STDOUT, "Unable to proceed (missing newline at position %d in sequence '", nsite + 1);
			fputid(STDOUT, notu);
			fprintf(STDOUT, "')\n\n\n");
#			if PARALLEL
				PP_Finalize();
#			endif
   			tp_exit(1, NULL, FALSE, __FILE__, __LINE__, exit_wait_optn);
		}
	} while (((char) ci != '\n') && ((char) ci != '\r'));
} /* skiprestofline */

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


/* skip control characters and blanks */
void skipcntrl(FILE *ifp,    /* input file stream */
               int   notu,   /* taxon number      - for error msg. */
               int   nsite)  /* sequence position - for error msg. */
{
	int ci;

	/* read over all control characters and blanks */
	do {
		ci = fgetc(ifp);
		if (ci == EOF) {
			fprintf(STDOUT, "\n\n\nUnable to proceed (missing character at position %d in sequence '", nsite + 1);
			fputid(STDOUT, notu);
			fprintf(STDOUT, "')\n\n\n");
#			if PARALLEL
				PP_Finalize();
#			endif
   			tp_exit(1, NULL, FALSE, __FILE__, __LINE__, exit_wait_optn);
		}
	} while (iscntrl(ci) || (char) ci == ' ');
	/* go one character back */
	if (ungetc(ci, ifp) == EOF) {
		fprintf(STDOUT, "\n\n\nUnable to proceed (positioning error at position %d in sequence '", nsite + 1);
		fputid(STDOUT, notu);
		fprintf(STDOUT, "')\n\n\n");
#		if PARALLEL
			PP_Finalize();
#		endif
   		tp_exit(1, NULL, FALSE, __FILE__, __LINE__, exit_wait_optn);
	}
} /* skipcntrl */

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


/* read sequences of one data set */
void getseqs(FILE    *ifp,      /* in:  input file stream */
             int      Maxspc,   /* in:  number of taxa */
             int      Maxseqc,  /* in:  number of sites */
             cmatrix *seqch,    /* out: alignment matrix */
             cmatrix  identif)  /* io:  taxon names (10 char w/o stop) */
{
	int notu, nsite, endofline, linelength, i;
	char c;
	cmatrix seqchars;
	
	seqchars = new_cmatrix(Maxspc, Maxseqc);
	/* read all characters */
	nsite = 0; /* next site to be read */
	while (nsite < Maxseqc) {
		/* read first taxon */
		notu = 0;

		/* only Maxspc, Maxseqc read so far: */
		/* go to next true line */
		skiprestofline(ifp, notu, nsite); 

		skipcntrl(ifp, notu, nsite);

		if (nsite == 0) readid(ifp, notu, identif);
		endofline = FALSE;
		linelength = 0;		
		do {
			c = readnextcharacter(ifp, notu, nsite + linelength);
			if ((c == '\n') || (c == '\r')) endofline = TRUE;
			else if (c == '.') {
				fprintf(STDOUT, "\n\n\nUnable to proceed (invalid character '.' at position ");
				fprintf(STDOUT, "%d in first sequence)\n\n\n", nsite + linelength + 1);
#				if PARALLEL
					PP_Finalize();
#				endif
   				tp_exit(1, NULL, FALSE, __FILE__, __LINE__, exit_wait_optn);
			     } else if (nsite + linelength < Maxseqc) {
				       /* change to upper case */
				       seqchars[notu][nsite + linelength] = (char) toupper((int) c);
				       linelength++;
			            } else {
				       endofline = TRUE;
				       skiprestofline(ifp, notu, nsite + linelength);
			            }
		} while (!endofline);	
		if (linelength == 0) {
			fprintf(STDOUT, "\n\n\nUnable to proceed (line with length 0 at position %d in sequence '", nsite + 1);
			fputid(STDOUT, notu);
			fprintf(STDOUT, "')\n\n\n");
#			if PARALLEL
				PP_Finalize();
#			endif
   			tp_exit(1, NULL, FALSE, __FILE__, __LINE__, exit_wait_optn);
		}
		/* read other taxa */
		for (notu = 1; notu < Maxspc; notu++) {
			/* go to next true line */
			if (notu != 1) skiprestofline(ifp, notu, nsite);
			skipcntrl(ifp, notu, nsite);
			if (nsite == 0) readid(ifp, notu, identif);
			for (i = nsite; i < nsite + linelength; i++) {
				c = readnextcharacter(ifp, notu, i);
				if ((c == '\n') || (c == '\r')) { /* too short */
					fprintf(STDOUT, "\n\n\nUnable to proceed (line to short at position %d in sequence '", i + 1);
					fputid(STDOUT, notu);
					fprintf(STDOUT, "')\n\n\n");
#					if PARALLEL
						PP_Finalize();
#					endif
   					tp_exit(1, NULL, FALSE, __FILE__, __LINE__, exit_wait_optn);
				} else if (c == '.') {
					seqchars[notu][i] = seqchars[0][i];
				} else {
					/* change to upper case */
					seqchars[notu][i] = (char) toupper((int) c);
				}
			}
		}
		nsite = nsite + linelength;
	}
	*seqch =seqchars;

} /* getseqs */

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


/* initialize identifer arrays */
void initid(int      t,        /* in:  number of taxa */
            cmatrix *identif,  /* out: name array w/o end of string */
            cmatrix *namestr)  /* out: name array with end of string */
{
	int i, j;

	*identif = new_cmatrix(t, 10);
	*namestr = new_cmatrix(t, 11);
	for (i = 0; i < t; i++) {
		(*namestr)[i][0] = '\0';
		(*namestr)[i][9] = '\0';
		for (j = 0; j < 10; j++) {
			(*identif)[i][j] = ' ';
		} /* for j */
	} /* for i */
} /* initid */

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


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


/* copy undelimited identifer array to '\0' delimited identifer array */
void identif2namestr(int     num,      /* number of taxa         */
                     cmatrix Identif,  /* non-delimited names    */
                     cmatrix Namestr)  /* proper delimited names */
{
	int i, j;
	int laggingspc;

	for (i = 0; i < num; i++) {
		laggingspc = 1;
		for (j = 9; j >= 0; j--) {
			Namestr[i][j] = Identif[i][j];
			if (laggingspc == 1) {
				if (Identif[i][j] != ' ') {
					Namestr[i][j+1] = '\0';
					laggingspc = 0;
				} 
			} 

		} /* for j (letter 10-1) */
	} /* for i (species) */
} /* identif2namestr */

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


/* print identifier of specified taxon in full 10 char length */
void fputid10(FILE *ofp, int t)
{	
	int i;

	for (i = 0; i < 10; i++) fputc(Identif[t][i], ofp);
} /* fputid10 */

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


/* print identifier of specified taxon up to first space */
int fputid(FILE *ofp, int t)
{	
	int i;
	
	i = 0;
/*
	while (i < 10) {
 		if (Identif[t][i] != ' ') {
			fputc(Identif[t][i], ofp);
		}
		i++;
	}
*/
	while (i < 10 && Identif[t][i] != ' ') {
		fputc(Identif[t][i], ofp);
		i++;
	}
	return i;
} /* fputid */

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


/* read first line of sequence data set */
void getsizesites(FILE *ifp,      /* in: input file stream */
                  int  *Maxspc,   /* out: number of taxa */
                  int  *Maxseqc)  /* out: number of sites */
{
	if (fscanf(ifp, "%d", Maxspc) != 1) {
			fprintf(STDOUT, "\n\n\nUnable to proceed (missing number of sequences)\n\n\n");
#			if PARALLEL
				PP_Finalize();
#			endif
   			tp_exit(1, NULL, FALSE, __FILE__, __LINE__, exit_wait_optn);
	}
	if (fscanf(ifp, "%d", Maxseqc) != 1) {
			fprintf(STDOUT, "\n\n\nUnable to proceed (missing number of sites)\n\n\n");
#			if PARALLEL
				PP_Finalize();
#			endif
   			tp_exit(1, NULL, FALSE, __FILE__, __LINE__, exit_wait_optn);
	}
	
	if (*Maxspc < 3) {
		fprintf(STDOUT, "\n\n\nUnable to proceed (less than 3 sequences - no tree possible)\n\n\n");
#		if PARALLEL
			PP_Finalize();
#		endif
   		tp_exit(1, NULL, FALSE, __FILE__, __LINE__, exit_wait_optn);
	}
	seqnumcheck = SEQNUM_OK;
	if (*Maxspc < 4) {
		fprintf(STDOUT, "\n\n\nLess than 4 sequences: No quartet methods possible!!!\n\n\n");
		fprintf(STDOUT, "Parameter estimation and tree testing is still possible.\n\n\n");
		seqnumcheck = SEQNUM_TOOFEW;
	}
#if 0
	if (*Maxspc > 257) 
#endif

	if (sizeof(uli) <= 4) {
		/* fprintf(stderr, "ULI=4\n"); */
		if (*Maxspc > 257) {
			fprintf(STDOUT, "\n\n\nMore than 257 sequences: No quartet puzzling available!!!\n");
			fprintf(STDOUT, "because your compiler can only handle ULIs of 4byte size!\n");
			fprintf(STDOUT, "It usually works on 64bit computers!\n\n");
			fprintf(STDOUT, "Parameter estimation, likelihood mapping, and tree testing is still \npossible for larger sets.\n\n\n");
			seqnumcheck = SEQNUM_TOOMANY;
		}
	} else {
		/* fprintf(stderr, "ULI=8\n"); */
		if (*Maxspc > 65538) {
			fprintf(STDOUT, "\n\n\nMore than 65538 sequences: No quartet puzzling available!!!\n");
			fprintf(STDOUT, "Parameter estimation, likelihood mapping, and tree testing is still possible.\n\n\n");
			seqnumcheck = SEQNUM_TOOMANY;
		}
	}
	/* correct default analysis types */
	if ((seqnumcheck!=SEQNUM_OK) && (puzzlemode==QUARTPUZ))
		puzzlemode = (puzzlemode + 1) % NUMPUZZLEMODES;
	if ((seqnumcheck==SEQNUM_TOOFEW) && (typ_optn==LIKMAPING_OPTN))
		typ_optn = (typ_optn + 1) % NUMTYPES;
#if 0
	if (*Maxspc < 4) {
		fprintf(STDOUT, "\n\n\nUnable to proceed (less than 4 sequences)\n\n\n");
#		if PARALLEL
			PP_Finalize();
#		endif
   		tp_exit(1, NULL, FALSE, __FILE__, __LINE__, exit_wait_optn);
	}
	if (*Maxspc > 257) {
		fprintf(STDOUT, "\n\n\nUnable to proceed (more than 257 sequences)\n\n\n");
#		if PARALLEL
			PP_Finalize();
#		endif
   		tp_exit(1, NULL, FALSE, __FILE__, __LINE__, exit_wait_optn);
	}
#		endif
	if (*Maxseqc < 1) {
		fprintf(STDOUT, "\n\n\nUnable to proceed (no sequence sites)\n\n\n");
#		if PARALLEL
			PP_Finalize();
#		endif
   		tp_exit(1, NULL, FALSE, __FILE__, __LINE__, exit_wait_optn);
	}
} /* getsizesites */

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


/* read alignment from file */
void readsequencefile(FILE    *seqfp,     /* in:  sequence input file stream */
                      int     *Maxspc,    /* out: number of taxa             */
                      int     *Maxseqc,   /* out: number of sites            */
                      cmatrix *identif,   /* out: name array w/o end of str  */
                      cmatrix *namestr,   /* out: name array with end of str */
                      cmatrix *Seqchars)  /* out: alignment matrix           */
{
	getsizesites(seqfp, Maxspc, Maxseqc);
	initid(*Maxspc, identif, namestr);
	getseqs(seqfp, *Maxspc, *Maxseqc, Seqchars, *identif);
	identif2namestr(*Maxspc, *identif, *namestr);
} /* readsequencefile */

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


/* read subsets from file */
void readsubsetfile(FILE    *seqfp,              /* in:  sequence input file stream */
                    int      Maxspc,             /* in:  number of taxa             */
                    cmatrix  namestr,            /* in:  names taxa (seq file)      */
                    int     *Maxsubset,          /* out: number of subsets          */
                    imatrix *ss_setovlgraph,     /* out: size of overlap >= 3 between 2 subsets */
                    imatrix *ss_setoverlaps,     /* out: size of overlap between 2 subsets */
                    imatrix *ss_setovllist,      /* out: list with ovlerlapping subsets */
                    ivector *ss_setovllistsize,  /* out: size of list with ovlerlapping subsets */
                    imatrix *ss_matrix,          /* out: boolean list: taxon in set? */
                    imatrix *ss_list,            /* out: list of taxa in set */
                    ivector *ss_listsize)         /* out: size of list with taxa */
{
	int ss_maxspc;
	int n;
	int set1, set2;
	int tmp;

        cmatrix  ss_identarr;
        cmatrix  ss_namearr;
        cmatrix  ss_strmatr;

	/* read the subset file sizes */
	getsizesites(seqfp, &ss_maxspc, Maxsubset);

	/* check sizes */
	if (ss_maxspc != Maxspc) {
		fprintf(STDOUT, "ERROR: Number of taxa in subsetfile does not match!!! (%d != %d)\n", ss_maxspc, Maxspc);
#		if PARALLEL
			PP_SendDone();
			MPI_Finalize();
#		endif /* PARALLEL */
   		tp_exit(1, NULL, FALSE, __FILE__, __LINE__, exit_wait_optn);
	} else {
		fprintf(STDOUT, "   %d species, %d subsets or genes\n", ss_maxspc, *Maxsubset);
	}

	/* initialize name arrays for checking */
	initid(ss_maxspc, &ss_identarr, &ss_namearr);

	/* read the data */
	getseqs(seqfp, ss_maxspc, *Maxsubset, &ss_strmatr, ss_identarr);

	/* check the file names */
	identif2namestr(ss_maxspc, ss_identarr, ss_namearr);
	for(n = 0; n<Maxspc; n++) {
		if (0 != strcmp(ss_namearr[n], Namestr[n])) {
		/* if (ss_maxspc != Maxspc) { */
			fprintf(STDOUT, "ERROR: Taxon names or order in subset file do not match!!!\n       (%d \"%s\" != \"%s\")\n", n+1, namestr[n], ss_namearr[n]);
#			if PARALLEL
				PP_SendDone();
				MPI_Finalize();
#			endif /* PARALLEL */
   			tp_exit(1, NULL, FALSE, __FILE__, __LINE__, exit_wait_optn);
		} else {
			fprintf(STDOUT, "   %3d. %s\n", n+1, ss_namearr[n]);

		}
		
	}
	fprintf(STDOUT, "\n");

	free_cmatrix(ss_identarr);
	free_cmatrix(ss_namearr);

	*ss_setovlgraph = new_imatrix(*Maxsubset, *Maxsubset);
	*ss_setoverlaps = new_imatrix(*Maxsubset, *Maxsubset);

	*ss_setovllist     = new_imatrix(*Maxsubset, *Maxsubset);
	*ss_setovllistsize = new_ivector(*Maxsubset);

	*ss_matrix         = new_imatrix(*Maxsubset, Maxspc);

	*ss_list           = new_imatrix(*Maxsubset, Maxspc);
	*ss_listsize       = new_ivector(*Maxsubset);

	for (n=0; n<Maxspc; n++) {
		for (set1=0; set1<*Maxsubset; set1++) {
			switch (ss_strmatr[n][set1]) {
				case 'X':
				case '1':
                    			(*ss_matrix)[set1][n] = 1;
                    			(*ss_list)[set1][((*ss_listsize)[set1])++] = n;
					break;
				case '-':
				case '0':
                    			(*ss_matrix)[set1][n] = 0;
					break;
				default:
					fprintf(STDOUT, "ERROR: Unknown character in subset file!!! (\"%c\", taxon %d, site %d)\n", ss_strmatr[n][set1], n+1, set1+1);
#					if PARALLEL
						PP_SendDone();
						MPI_Finalize();
#					endif /* PARALLEL */
   					tp_exit(1, NULL, FALSE, __FILE__, __LINE__, exit_wait_optn);
					break;
			}
		}

		for (set1=1; set1<*Maxsubset; set1++) {
			for (set2=0; set2<set1; set2++) {
				if (((*ss_matrix)[set1][n]) && ((*ss_matrix)[set2][n])) 
					/* number of sequences in overlap */
                    			((*ss_setoverlaps)[set1][set2])++;
					
			}
		}
	} /* for Maxspc */

	for (set1=1; set1<*Maxsubset; set1++) {
		for (set2=0; set2<set1; set2++) {
			/* number of sequences in overlap */
                   	tmp = (*ss_setoverlaps)[set1][set2];
                   	(*ss_setoverlaps)[set2][set1] = tmp;
			if (tmp >= 3) {
				/* number of sequences (>=3) in overlap */
                   		(*ss_setovlgraph)[set1][set2] = tmp;
                   		(*ss_setovlgraph)[set2][set1] = tmp;
				/* list of overlapping sets */
                    		(*ss_setovllist)[set1][((*ss_setovllistsize)[set1])++] = set2;
                    		(*ss_setovllist)[set2][((*ss_setovllistsize)[set2])++] = set1;
			}
				
		}
	}
} /* readsubsetfile */

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


/* permute taxon order */
void permutetaxa_ss(int      Maxspc,             /* in:  number of taxa             */
                    ivector  permutation,        /* permuted taxon order            */
                    cmatrix  namestr,            /* in:  names taxa (seq file)      */
                    int      Maxsubset,          /* out: number of subsets          */
                    imatrix  ss_setovlgraph,     /* out: size of overlap >= 3 between 2 subsets */
                    imatrix  ss_setoverlaps,     /* out: size of overlap between 2 subsets */
                    imatrix  ss_setovllist,      /* out: list with ovlerlapping subsets */
                    ivector  ss_setovllistsize,  /* out: size of list with ovlerlapping subsets */
                    imatrix  ss_matrix,          /* out: boolean list: taxon in set? */
                    imatrix  ss_list,            /* out: list of taxa in set */
                    ivector  ss_listsize)         /* out: size of list with taxa */
{
	ivector candidatesetlist;
	int     candidatesetsize = 0;
	ivector candidatetaxalist;
	int     candidatetaxasize = 0;
	ivector taxonstatus;			/* 0 - not done; 1 - candidate; 2 - inserted */
	ivector setstatus;			/* 0 - not used; 1 - front; 2 - used*/
	int     permpos;
	int     k;
	int     tmprand;
	int     currset;
	int     currtaxon;
	candidatesetlist  = new_ivector(Maxsubset);
	candidatetaxalist = new_ivector(Maxspc);
	taxonstatus       = new_ivector(Maxspc);
	setstatus         = new_ivector(Maxsubset);

	/* first subset */
	currset = randominteger(Maxsubset);
	setstatus[currset] = 2; /* inserted */

	/* init subset front nodes */
	/*   insert sets with ovl(currset) >= 3 into candidate set list */
	for (k=0; k<ss_setovllistsize[currset]; k++) {
		candidatesetlist[candidatesetsize] = ss_setovllist[currset][k];
		setstatus[candidatesetlist[candidatesetsize++]] = 1;
	}
	/*   insert taxa from currset into candidate taxon list */
	for (k=0; k<ss_listsize[currset]; k++) {
		candidatetaxalist[candidatetaxasize] = ss_list[currset][k];
		taxonstatus[candidatetaxalist[candidatetaxasize++]] = 1;
	}

	/* start collectin permutation */
	for (permpos=0; permpos<Maxspc; permpos++) {
		/* if no taxa left, choose next set */
		while (candidatetaxasize == 0) {
			/* next set */
			tmprand = randominteger(candidatesetsize);
			currset = candidatesetlist[tmprand];
			/* move last set entry to gap */
			candidatesetlist[tmprand] = candidatesetlist[--candidatesetsize]; /* XXXX */
			/* candidatesetlist[tmprand] = candidatesetlist[candidatesetsize--]; */
			setstatus[currset] = 2; /* inserted */
			/* add new sets (ovl>=3) to set candidate list */
			for (k=0; k<ss_setovllistsize[currset]; k++) {
				if (setstatus[ss_setovllist[currset][k]] == 0) {
					candidatesetlist[candidatesetsize] = ss_setovllist[currset][k];
					setstatus[candidatesetlist[candidatesetsize++]] = 1;
				}
			}
			/* add new taxa to taxon candidate list */
			for (k=0; k<ss_listsize[currset]; k++) {
				if (taxonstatus[ss_list[currset][k]] == 0) {
					candidatetaxalist[candidatetaxasize] = ss_list[currset][k];
					taxonstatus[candidatetaxalist[candidatetaxasize++]] = 1;
				}
			}

		} /* if candidatetaxasize == 0 */

		tmprand = randominteger(candidatetaxasize);
		currtaxon = candidatetaxalist[tmprand];

		candidatetaxalist[tmprand] = candidatetaxalist[--candidatetaxasize];
		permutation[permpos] = currtaxon;
		taxonstatus[currtaxon] = 2;

	} /* one random taxon after the other */
	
	free_ivector(candidatesetlist);
	free_ivector(candidatetaxalist);
	free_ivector(taxonstatus);
	free_ivector(setstatus);
} /* permutetaxa_ss */

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


/* permute taxon order */
void permutetaxa_weighted_ss(int      Maxspc,             /* in:  number of taxa             */
                    ivector  permutation,        /* permuted taxon order            */
                    cmatrix  namestr,            /* in:  names taxa (seq file)      */
                    int      Maxsubset,          /* out: number of subsets          */
                    imatrix  ss_setovlgraph,     /* out: size of overlap >= 3 between 2 subsets */
                    imatrix  ss_setoverlaps,     /* out: size of overlap between 2 subsets */
                    imatrix  ss_setovllist,      /* out: list with ovlerlapping subsets */
                    ivector  ss_setovllistsize,  /* out: size of list with ovlerlapping subsets */
                    imatrix  ss_matrix,          /* out: boolean list: taxon in set? */
                    imatrix  ss_list,            /* out: list of taxa in set */
                    ivector  ss_listsize)         /* out: size of list with taxa */
{
	ivector candidatesetlist;		/* vector containing candidate sets possible to be added next */
	int     candidatesetsize = 0;		/* its size (# sets in front list) */
	ivector candidatetaxalist;		/* vector containing candidate taxa to be added next */
	int     candidatetaxasize = 0;		/* its size (# taxa in added front list) */
	ivector taxonstatus;			/* 0 - not done; 1 - candidate; 2 - inserted */
	ivector setstatus;			/* 0 - not used; 1 - front; 2 - used*/
	int     permpos;
	int     k;
	int     tmprand;
	int     currset;
	int     currtaxon;
	candidatesetlist  = new_ivector(Maxsubset);
	candidatetaxalist = new_ivector(Maxspc);
	taxonstatus       = new_ivector(Maxspc);
	setstatus         = new_ivector(Maxsubset);

	/* first subset */
/* INITIALIZE SIZE VECTOR (HAS ;-) */
	currset = randominteger(Maxsubset);
/* RANDOM WEIGHTED BY SUBSET SIZE (HAS ;-) */
	setstatus[currset] = 2; /* inserted */

	/* init subset front nodes */
	/*   insert sets with ovl(currset) >= 3 into candidate set list */
	for (k=0; k<ss_setovllistsize[currset]; k++) {
		candidatesetlist[candidatesetsize] = ss_setovllist[currset][k];
		setstatus[candidatesetlist[candidatesetsize++]] = 1;
/* INITIALIZE OVERLAP SIZE VECTOR (HAS ;-) */
/* MAX EDGE: VECTOR JUST CONTAINS THE HIGHEST OVERLAP OF x TO ANY SET USED (HAS ;-) */
/* SUM EDGE: VECTOR CONTAINS THE OVERLAP TO THE CURRENT TAXON SET  (HAS ;-) */
	}
	/*   insert taxa from currset into candidate taxon list */
	for (k=0; k<ss_listsize[currset]; k++) {
		candidatetaxalist[candidatetaxasize] = ss_list[currset][k];
		taxonstatus[candidatetaxalist[candidatetaxasize++]] = 1;
	}

	/* start collecting permutation */
	for (permpos=0; permpos<Maxspc; permpos++) {
		/* if no taxa left, choose next set */
		while (candidatetaxasize == 0) {
			/* next set */
/* INITIALIZE OVERLAP SIZE VECTOR (HAS ;-) */
			tmprand = randominteger(candidatesetsize);
/* RANDOM WEIGHTED BY OVERLAP SIZE (MAX EDGE, SUM EDGE) (HAS ;-) */
/* UPDATE OVERLAP SIZE VECTOR (HAS ;-) */
/* MAX EDGE: VECTOR UPDATED ONLY IF AN OVERLAP TO THE 'NEXT SET' IS LARGER THAN CURRENT (HAS ;-) */
/* SUM EDGE: VECTOR UPDATED IF SOME TAXA IN 'NEXT SET' INCREASE THE OVERLAP TO THE CURRENT TAXON SET  (HAS ;-) */
			currset = candidatesetlist[tmprand];
			/* move last set entry to gap */
			candidatesetlist[tmprand] = candidatesetlist[--candidatesetsize]; /* XXXX */
			/* candidatesetlist[tmprand] = candidatesetlist[candidatesetsize--]; */
			setstatus[currset] = 2; /* inserted */
			/* add new sets (ovl>=3) to set candidate list */
			for (k=0; k<ss_setovllistsize[currset]; k++) {
				if (setstatus[ss_setovllist[currset][k]] == 0) {
					candidatesetlist[candidatesetsize] = ss_setovllist[currset][k];
					setstatus[candidatesetlist[candidatesetsize++]] = 1;
				}
			}
			/* add new taxa to taxon candidate list */
			for (k=0; k<ss_listsize[currset]; k++) {
				if (taxonstatus[ss_list[currset][k]] == 0) {
					candidatetaxalist[candidatetaxasize] = ss_list[currset][k];
					taxonstatus[candidatetaxalist[candidatetaxasize++]] = 1;
				}
			}

		} /* if candidatetaxasize == 0 */

		tmprand = randominteger(candidatetaxasize);
		currtaxon = candidatetaxalist[tmprand];

		candidatetaxalist[tmprand] = candidatetaxalist[--candidatetaxasize];
		permutation[permpos] = currtaxon;
		taxonstatus[currtaxon] = 2;

	} /* one random taxon after the other */
	
	free_ivector(candidatesetlist);
	free_ivector(candidatetaxalist);
	free_ivector(taxonstatus);
	free_ivector(setstatus);
} /* permutetaxa_weighted_ss */

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


/* print subsets */
void fprintfss(FILE    *ofp,              /* in:  output file stream         */
             int      Maxspc,             /* in:  number of taxa             */
             int      Maxsubset,          /* out: number of subsets          */
             imatrix  ss_setovlgraph,     /* out: size of overlap >= 3 between 2 subsets */
             imatrix  ss_setoverlaps,     /* out: size of overlap between 2 subsets */
             imatrix  ss_setovllist,      /* out: list with ovlerlapping subsets */
             ivector  ss_setovllistsize,  /* out: size of list with ovlerlapping subsets */
             imatrix  ss_matrix,          /* out: boolean list: taxon in set? */
             imatrix  ss_list,            /* out: list of taxa in set */
             ivector  ss_listsize)         /* out: size of list with taxa */
{
	int     s, s1, s2, t;


	for (s=0; s<Maxsubset; s++) {
		fprintf(ofp, "   set %d (%d taxa)\n", s, ss_listsize[s]);
			fprintf(ofp, "      taxa:        ");
		for (t=0; t<ss_listsize[s]; t++) {
			fprintf(ofp, "%3d", ss_list[s][t]);
			if (t+1 < ss_listsize[s]) fprintf(ofp, ",");
			if (((t+1) % 15) == 0)
				fprintf(ofp, "\n                   ");
		}
		fprintf(ofp, "\n      overlap with:");
		for (t=0; t<ss_setovllistsize[s]; t++) {
			fprintf(ofp, "%3d", ss_setovllist[s][t]);
			if (t+1 < ss_setovllistsize[s]) fprintf(ofp, ",");
			if (((t+1) % 15) == 0)
				fprintf(ofp, "\n                   ");
		}
		fprintf(ofp, "\n\n");
	}
	fprintf(ofp, "   Overlap Graph:\n");
	for (s1=0; s1<Maxsubset; s1++) {
		fprintf(ofp, "   %3d :\t", s1);
		for (s2=0; s2<Maxsubset; s2++) {
			if (s1==s2) 
				fprintf(ofp, "(%3d) ", ss_listsize[s1]);
			else 
				fprintf(ofp, " %3d  ", ss_setoverlaps[s1][s2]);
		}
		fprintf(ofp, "\n");
	}
	fprintf(ofp, "\n\n");
	fflush(ofp);
} /* fprintfss */

/******************/
/* check connectedness of subsets */
void checkss(FILE    *ofp,              /* in:  output file stream         */
             int      Maxspc,             /* in:  number of taxa             */
             int      Maxsubset,          /* out: number of subsets          */
             imatrix  ss_setovlgraph,     /* out: size of overlap >= 3 between 2 subsets */
             imatrix  ss_setoverlaps,     /* out: size of overlap between 2 subsets */
             imatrix  ss_setovllist,      /* out: list with ovlerlapping subsets */
             ivector  ss_setovllistsize,  /* out: size of list with ovlerlapping subsets */
             imatrix  ss_matrix,          /* out: boolean list: taxon in set? */
             imatrix  ss_list,            /* out: list of taxa in set */
             ivector  ss_listsize)         /* out: size of list with taxa */
{
	int     s1, s2;
	int numdone;
        int *setarr;
        setarr = new_ivector(Maxsubset);

	s1 = 0;
	numdone=0;
	do {
		setarr[s1] = 2; /* 1=seen, 2=visited */
		for (s2=0; s2<Maxsubset; s2++) {
			if ((setarr[s2]==0) && (ss_setoverlaps[s1][s2] >= 3))
		  	setarr[s2]=1; 
		}

		numdone++;
		s1=0;
		while ((s1 < Maxsubset) && (setarr[s1]!=1)) s1++;

	} while (s1 != Maxsubset);

	if (Maxsubset != numdone) {
		fprintf(stderr, "Unable to proceed. More than one connected component! (%d != %d)\n", Maxsubset, numdone);
		fprintf(stderr, "visited: \n  ");
		for (s2=0; s2<Maxsubset; s2++) 
			if (setarr[s2]==2) fprintf(stderr, "%d, ", s2);
		fprintf(stderr, "\nseen: \n  ");
		for (s2=0; s2<Maxsubset; s2++) 
			if (setarr[s2]==1) fprintf(stderr, "%d, ", s2);
		fprintf(stderr, "\nlost: \n  ");
		for (s2=0; s2<Maxsubset; s2++) 
			if (setarr[s2]==0) fprintf(stderr, "%d, ", s2);
   		tp_exit(1, NULL, FALSE, __FILE__, __LINE__, exit_wait_optn);
	} 
} /* checkss */

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


/* guess data type: NUCLEOTIDE:0, AMINOACID:1, BINARY:2 */
int guessdatatype(cmatrix Seqchars,   /* alignment matrix (Maxspc x Maxseqc) */
                  int     Maxspc,     /* number of taxa */
                  int     Maxseqc)    /* number of sites */
{
	uli numnucs, numchars, numbins;
	int notu, nsite;
	char c;
	
	/* count A, C, G, T, U, N */
	numnucs = 0;
	numchars = 0;
	numbins = 0;
	for (notu = 0; notu < Maxspc; notu++)
		for (nsite = 0; nsite < Maxseqc; nsite++) {
			c = Seqchars[notu][nsite];
			if (c == 'A' || c == 'C' || c == 'G' ||
			    c == 'T' || c == 'U' || c == 'N') numnucs++;
			if (c != '-' && c != '?') numchars++;
			if (c == '0' || c == '1') numbins++;
		}
	if (numchars == 0) numchars = 1;
	/* more than 85 % frequency means nucleotide data */
	if ((double) numnucs / (double) numchars > 0.85) 
		return NUCLEOTIDE; /* 0 */
	else if ((double) numbins / (double) numchars > 0.2) 
		return BINARY; /* 2 */
	else 
		return AMINOACID; /* 1 */
} /* guessdatatype */

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


/* translate characters into format used by ML engine */
void translatedataset(int      maxspc, 
		      int      maxseqc, 
		      int     *maxsite, 
		      cmatrix  seqchars, 
		      cmatrix *seqchar, 
		      ivector *seqgapchar, 
		      ivector *seqotherchar)
{	
	int notu, sn, co;
	char c;
	cvector code;
	
	if (*seqgapchar   != NULL) free_ivector(*seqgapchar);
	if (*seqotherchar != NULL) free_ivector(*seqotherchar);
	Seqgapchar   = new_ivector(Maxspc);
	Seqotherchar = new_ivector(Maxspc);

	/* Seqgapchar   = (int *) calloc(maxspc, sizeof(int)); */
	/* Seqotherchar = (int *) calloc(maxspc, sizeof(int)); */


	/* determine maxsite - number of ML sites per taxon */
	if (data_optn == NUCLEOTIDE && SH_optn) {
		if (SHcodon)
			*maxsite = maxseqc / 3;
		else
			*maxsite = maxseqc / 2; /* assume doublets */
		
	} else
		*maxsite = maxseqc;

	if (data_optn == NUCLEOTIDE && (*maxsite % 3) == 0  && !SH_optn) {	
		if (codon_optn == 1 || codon_optn == 2 || codon_optn == 3)
			*maxsite = *maxsite / 3; /* only one of the three codon positions */
		if (codon_optn == 4)
			*maxsite = 2*(*maxsite / 3); /* 1st + 2nd codon positions */
	}
	
	/* allocate memory */
	if (*seqchar != NULL) free_cmatrix(*seqchar);
	*seqchar = new_cmatrix(maxspc, *maxsite);

	/* code length */
	if (data_optn == NUCLEOTIDE && SH_optn)
		code = new_cvector(2);
	else
		code = new_cvector(1);
	
	/* decode characters */
	if (data_optn == NUCLEOTIDE && SH_optn) { /* SH doublets */
		
		for (notu = 0; notu < maxspc; notu++) {
			for (sn = 0; sn < *maxsite; sn++) {
					for (co = 0; co < 2; co++) {
						if (SHcodon)
							c = seqchars[notu][sn*3 + co];
						else
							c = seqchars[notu][sn*2 + co];
						code[co] = c;
					}
				(*seqchar)[notu][sn] = code2int(code, &((*seqgapchar)[notu]), &((*seqotherchar)[notu]));
			}
		}
		
	} else if (!(data_optn == NUCLEOTIDE && (maxseqc % 3) == 0)) { /* use all */

		for (notu = 0; notu < maxspc; notu++) {
			for (sn = 0; sn < *maxsite; sn++) {
				code[0] = seqchars[notu][sn];
				(*seqchar)[notu][sn] = code2int(code, &((*seqgapchar)[notu]), &((*seqotherchar)[notu]));
			}
		}

	} else { /* codons */
		
		for (notu = 0; notu < maxspc; notu++) {
			for (sn = 0; sn < *maxsite; sn++) {
				if (codon_optn == 1 || codon_optn == 2 || codon_optn == 3)
					code[0] = seqchars[notu][sn*3+codon_optn-1];
				else if (codon_optn == 4) {
					if ((sn % 2) == 0)
						code[0] = seqchars[notu][(sn/2)*3];
					else
						code[0] = seqchars[notu][((sn-1)/2)*3+1];
				} else
					code[0] = seqchars[notu][sn];
				(*seqchar)[notu][sn] = code2int(code, &((*seqgapchar)[notu]), &((*seqotherchar)[notu]));
			}
		}
	
	}

	free_cvector(code);

} /* translatedataset */

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



/* estimate mean base frequencies from translated data set */
void estimatebasefreqs()
/* void estimatebasefreqs(int countsitefreqs) */
{
	int tpmradix, i, j;
	uli all, *gene, *sitefreqs;
int countsitefreqs=savesitefreqs_optn;
FILE *sffp;
	
	tpmradix = gettpmradix();
	
	if (Freqtpm != NULL) free_dvector(Freqtpm);
	Freqtpm = new_dvector(tpmradix);
	
	if (Basecomp != NULL) free_imatrix(Basecomp);
	Basecomp = new_imatrix(Maxspc, tpmradix);
	
	/* close sitefreqs file */
	if (countsitefreqs) {
		openfiletowrite(&sffp, SITEFREQS, "site frequencies", stdinput_fp);
		/* openfiletowrite(&sffp, "output.sitefreqs", "site frequencies", stdinput_fp); */
		/* output header line */
		fprintf(sffp, "site");
		for (i = 0; i < tpmradix; i++) 
			fprintf(sffp, "\t%s", int2code(i));
		fprintf(sffp, "\tother\n");
	}

	/* alloc vector to count overall frequencies */
	gene = (uli *) calloc((size_t) (tpmradix + 1), sizeof(uli));
	if (gene == NULL) maerror("gene in estimatebasefreqs");
	
	/* alloc vector to count site frequencies */
	sitefreqs = (uli *) calloc((size_t) (tpmradix + 1), sizeof(uli));
	if (sitefreqs == NULL) maerror("sitefreqs in estimatebasefreqs");
	
	for (i = 0; i < tpmradix + 1; i++) gene[i] = 0;
	for (i = 0; i < Maxspc; i++)
		for (j = 0; j < tpmradix; j++) Basecomp[i][j] = 0;
/* 
	for (i = 0; i < Maxspc; i++)
		for (j = 0; j < tpmradix; j++) Basecomp[i][j] = 0;
*/ 
	for (j = 0; j < Maxsite; j++) { /* site by site */

		/* reset sitefreq vector for current site */
		for (i = 0; i < tpmradix + 1; i++) sitefreqs[i] = 0;

		for (i = 0; i < Maxspc; i++) { /* species by species */
			/* count overall freqs */
			gene[ (int) Seqchar[i][j] ]++;
			/* count site freqs */
			sitefreqs[ (int) Seqchar[i][j] ]++;
			if (Seqchar[i][j] != tpmradix) { /* normal character, otherwise gap/wildcard */
				Basecomp[i][(int) Seqchar[i][j]]++;
			}
		}

		/* output sitefreqs */
		if (countsitefreqs) {
			fprintf(sffp, "%d", j+1);
			for (i = 0; i < tpmradix + 1; i++) 
				fprintf(sffp, "\t%.4f", (double)sitefreqs[i]/Maxspc);
			fprintf(sffp, "\n");
		}

	} /* end for all site */

	all = Maxspc * Maxsite - gene[tpmradix];
	if (all != 0) { /* normal case (i.e. not all characters are wildcards) */
		for (i = 0; i < tpmradix; i++)
			Freqtpm[i] = (double) gene[i] / (double) all;
	} else { /* pathological case with only wildcards in data set? -> equal freqs */
		for (i = 0; i < tpmradix; i++)
			Freqtpm[i] = 1.0 / (double) tpmradix;
	}
	
	/* close sitefreqs file */
	if (countsitefreqs) {
		fclose(sffp);
	}

	free(gene);
	
	Frequ_optn = TRUE;
} /* estimatebasefreqs */

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


#if 0
/* count base substitution frequencies from translated data set */
void count_tstv(ulimatrix *ts, ulimatrix *tv)
{
	int seq1, seq2, x, i, j;
	ulimatrix ts_counts, tv_counts, countmatr;
	
	if (data_optn == NUCLEOTIDE) { /* nucleotides */
		countmatr = new_ulimatrix(5,5);
		countmatr = new_ulimatrix(5,5);
		ts_counts = new_ulimatrix(Maxspc,Maxspc);
		tv_counts = new_ulimatrix(Maxspc,Maxspc);
		
		/* 0 = A */
		/* 1 = C */
		/* 2 = G */
		/* 3 = T/U */

		/* ts: 0-2 + 1-3 */
		/* tv: 0-1 + 0-3 + 2-1 + 2-3 */

		for (seq1 = 0; seq1 < Maxspc; seq1++) {
			for (seq2 = seq1+1; seq2 < Maxspc; seq2++) {
				for (i = 0; i < 5; i++) {
					for (j = 0; j < 5; j++) {
						countmatr[i][j]=0;
					}
				}

				for (x = 0; x < Maxsite; x++) {
					countmatr[(int) Seqchar[seq1][x]][(int) Seqchar[seq2][x]]++;
				}

				ts_counts[seq1][seq2] = countmatr[0][2] + countmatr[2][0] + 
							countmatr[1][3] + countmatr[3][1];
				ts_counts[seq2][seq1] = ts_counts[seq1][seq2];
				tv_counts[seq1][seq2] = countmatr[0][1] + countmatr[1][0] + 
							countmatr[0][3] + countmatr[3][0] + 
							countmatr[2][1] + countmatr[1][2] + 
							countmatr[2][3] + countmatr[3][2];
				tv_counts[seq2][seq1] = tv_counts[seq1][seq2];
			}
		}

		free_ulimatrix(countmatr);
		*ts = ts_counts;
		*tv = tv_counts;

	} 
	if (data_optn == AMINOACID) { /* amino acids */
		/*'A': 0 */
		/*'R': 1 */
		/*'N': 2 */
		/*'D': 3 */
		/*'C': 4 */
		/*'Q': 5 */
		/*'E': 6 */
		/*'G': 7 */
		/*'H': 8 */
		/*'I': 9 */
		/*'K': 11 */
		/*'L': 10 */
		/*'M': 12 */
		/*'F': 13 */
		/*'P': 14 */
		/*'S': 15 */
		/*'T': 16 */
		/*'W': 17 */
		/*'Y': 18 */
		/*'V': 19 */
		/*'-', 'X': 20 */

	}


} /* count_tstv */
#endif

/******************/
/* estimate mean base frequencies from translated data set */
void output_colfeqs(ulimatrix *ts, ulimatrix *tv)
{
	int seq1, seq2, x, i, j;
	ulimatrix ts_counts, tv_counts, countmatr;
	
	if (data_optn == NUCLEOTIDE) { /* nucleotides */
		countmatr = new_ulimatrix(5,5);
		ts_counts = new_ulimatrix(Maxspc,Maxspc);
		tv_counts = new_ulimatrix(Maxspc,Maxspc);
		
		/* 0 = A */
		/* 1 = C */
		/* 2 = G */
		/* 3 = T/U */

		/* ts: 0-2 + 1-3 */
		/* tv: 0-1 + 0+3 + 2-1 + 2-3 */

		for (seq1 = 0; seq1 < Maxspc; seq1++) {
			for (seq2 = seq1+1; seq2 < Maxspc; seq2++) {
				for (i = 0; i < 5; i++) {
					for (j = 0; j < 5; j++) {
						countmatr[i][j]=0;
					}
				}

				for (x = 0; x < Maxsite; x++) {
					countmatr[(int) Seqchar[seq1][x]][(int) Seqchar[seq2][x]]++;
				}

				ts_counts[seq1][seq2] = countmatr[0][2] + countmatr[2][0] + 
							countmatr[1][3] + countmatr[3][1];
				ts_counts[seq2][seq1] = ts_counts[seq1][seq2];
				tv_counts[seq1][seq2] = countmatr[0][1] + countmatr[1][0] + 
							countmatr[0][3] + countmatr[3][0] + 
							countmatr[2][1] + countmatr[1][2] + 
							countmatr[2][3] + countmatr[3][2];
				tv_counts[seq2][seq1] = tv_counts[seq1][seq2];
			}
		}

		free_ulimatrix(countmatr);
		*ts = ts_counts;
		*tv = tv_counts;

	}


} /* output_colfeqs */

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


/* guess model of substitution */
void guessmodel()
{
	double c1, c2, c3, c4, c5, c6, c7;
	dvector f;
	dmatrix a;
	int i;

	Dayhf_optn = FALSE;
	Jtt_optn = TRUE;
	mtrev_optn = FALSE;
	cprev_optn = FALSE;
	blosum62_optn = FALSE;
	vtmv_optn = FALSE;
	wag_optn = FALSE;
	lg_optn = FALSE;
	/* poisson_optn = FALSE; */
	/* equalin_optn = FALSE; */
	/* mtmam_optn = FALSE; */
	/* rtrev_optn = FALSE; */
	TSparam = 2.0;
	YRparam = 1.0;
	optim_optn = TRUE;
	HKY_optn = TRUE;
	print_GTR_optn = TRUE;
	GTR_optn = FALSE;
	TN_optn = FALSE;
	
	if (data_optn == AMINOACID) { /* amino acids */
		
		/* chi2 fit to amino acid frequencies */
		
		f = new_dvector(20);
		a = new_dmatrix(20,20);
		/* chi2 distance Dayhoff */
		dyhfdata(a, f);
		c1 = 0;
		for (i = 0; i < 20; i++)
			c1 = c1 + (Freqtpm[i]-f[i])*(Freqtpm[i]-f[i]);
		/* chi2 distance JTT */
		jttdata(a, f);
		c2 = 0;
		for (i = 0; i < 20; i++)
			c2 = c2 + (Freqtpm[i]-f[i])*(Freqtpm[i]-f[i]);
		/* chi2 distance mtREV */
		mtrevdata(a, f);
		c3 = 0;
		for (i = 0; i < 20; i++)
			c3 = c3 + (Freqtpm[i]-f[i])*(Freqtpm[i]-f[i]);
		/* chi2 distance VT */
		vtmvdata(a, f);
		c4 = 0;
		for (i = 0; i < 20; i++)
			c4 = c4 + (Freqtpm[i]-f[i])*(Freqtpm[i]-f[i]);
		/* chi2 distance WAG */
		wagdata(a, f);
		c5 = 0;
		for (i = 0; i < 20; i++)
			c5 = c5 + (Freqtpm[i]-f[i])*(Freqtpm[i]-f[i]);
		/* chi2 distance cpREV */
		cprev45data(a, f);
		c6 = 0;
		for (i = 0; i < 20; i++)
			c6 = c6 + (Freqtpm[i]-f[i])*(Freqtpm[i]-f[i]);
		/* chi2 distance LG */
		lgdata(a, f);
		c7 = 0;
		for (i = 0; i < 20; i++)
			c7 = c7 + (Freqtpm[i]-f[i])*(Freqtpm[i]-f[i]);

		free_dvector(f);
		free_dmatrix(a);

#ifndef CPREV
		if ((c1 < c2) && (c1 < c3) && (c1 < c4) && (c1 < c5)) {
		        /* c1 -> Dayhoff */
			Dayhf_optn = TRUE;
			Jtt_optn = FALSE;
			mtrev_optn = FALSE;
			cprev_optn = FALSE;
			vtmv_optn = FALSE;
			wag_optn = FALSE;
			lg_optn = FALSE;
			fprintf(STDOUT, "(consists very likely of amino acids encoded on nuclear DNA)\n");
		} else {
			if ((c2 < c3) && (c2 < c4) && (c2 < c5)) {
				/* c2 -> JTT */
				Dayhf_optn = FALSE;
				Jtt_optn = TRUE;
				mtrev_optn = FALSE;
				cprev_optn = FALSE;
				vtmv_optn = FALSE;
				wag_optn = FALSE;
				lg_optn = FALSE;
				fprintf(STDOUT, "(consists very likely of amino acids encoded on nuclear DNA)\n");
			} else {
				if ((c3 < c4) && (c3 < c5)) {
					/* c3 -> mtREV */
					Dayhf_optn = FALSE;
					Jtt_optn = FALSE;
					mtrev_optn = TRUE;
					cprev_optn = FALSE;
					vtmv_optn = FALSE;
					wag_optn = FALSE;
					lg_optn = FALSE;
					fprintf(STDOUT, "(consists very likely of amino acids encoded on mtDNA)\n");
				} else {
					if ((c4 < c5)) {
						/* c4 -> VT */
						Dayhf_optn = FALSE;
						Jtt_optn = FALSE;
						mtrev_optn = FALSE;
						cprev_optn = FALSE;
						vtmv_optn = TRUE;
						wag_optn = FALSE;
						lg_optn = FALSE;
						fprintf(STDOUT, "(consists very likely of amino acids encoded on nuclear DNA)\n");
					} else {
						if ((c5 < c7)) {
							/* c5 -> WAG */
							Dayhf_optn = FALSE;
							Jtt_optn = FALSE;
							mtrev_optn = FALSE;
							cprev_optn = FALSE;
							vtmv_optn = FALSE;
							wag_optn = TRUE;
							lg_optn = FALSE;
							fprintf(STDOUT, "(consists very likely of amino acids encoded on nuclear DNA)\n");
						} else {
								/* missing c6, cpREV */
								/* c7 -> LG */
								Dayhf_optn = FALSE;
								Jtt_optn = FALSE;
								mtrev_optn = FALSE;
								cprev_optn = FALSE;
								vtmv_optn = FALSE;
								wag_optn = FALSE;
								lg_optn = TRUE;
								fprintf(STDOUT, "(consists very likely of amino acids encoded on nuclear DNA)\n");
						} /* if c5 else c7 */
					} /* if c4 else c5 */
				} /* if c3 else c4 */
			} /* if c2 */
		} /* if c1 */

#else /* CPREV */

		if ((c1 < c2) && (c1 < c3) && (c1 < c4) && (c1 < c5) && (c1 < c6)) {
		        /* c1 -> Dayhoff */
			Dayhf_optn = TRUE;
			Jtt_optn = FALSE;
			mtrev_optn = FALSE;
			cprev_optn = FALSE;
			vtmv_optn = FALSE;
			wag_optn = FALSE;
			lg_optn = FALSE;
			fprintf(STDOUT, "(consists very likely of amino acids encoded on nuclear DNA)\n");
		} else {
			if ((c2 < c3) && (c2 < c4) && (c2 < c5) && (c2 < c6)) {
				/* c2 -> JTT */
				Dayhf_optn = FALSE;
				Jtt_optn = TRUE;
				mtrev_optn = FALSE;
				cprev_optn = FALSE;
				vtmv_optn = FALSE;
				wag_optn = FALSE;
				lg_optn = FALSE;
				fprintf(STDOUT, "(consists very likely of amino acids encoded on nuclear DNA)\n");
			} else {
				if ((c3 < c4) && (c3 < c5) && (c3 < c6)) {
					/* c3 -> mtREV */
					Dayhf_optn = FALSE;
					Jtt_optn = FALSE;
					mtrev_optn = TRUE;
					cprev_optn = FALSE;
					vtmv_optn = FALSE;
					wag_optn = FALSE;
					lg_optn = FALSE;
					fprintf(STDOUT, "(consists very likely of amino acids encoded on mtDNA)\n");
				} else {
					if ((c4 < c5) && (c4 < c6)) {
						/* c4 -> VT */
						Dayhf_optn = FALSE;
						Jtt_optn = FALSE;
						mtrev_optn = FALSE;
						cprev_optn = FALSE;
						vtmv_optn = TRUE;
						wag_optn = FALSE;
						lg_optn = FALSE;
						fprintf(STDOUT, "(consists very likely of amino acids encoded on nuclear DNA)\n");
					} else {
						if (c5 < c6) {
							/* c5 -> WAG */
							Dayhf_optn = FALSE;
							Jtt_optn = FALSE;
							mtrev_optn = FALSE;
							cprev_optn = FALSE;
							vtmv_optn = FALSE;
							wag_optn = TRUE;
							lg_optn = FALSE;
							fprintf(STDOUT, "(consists very likely of amino acids encoded on nuclear DNA)\n");
						} else {
							if (c6 < c7) {
								/* c6 -> cpREV */
								Dayhf_optn = FALSE;
								Jtt_optn = FALSE;
								mtrev_optn = FALSE;
								cprev_optn = TRUE;
								vtmv_optn = FALSE;
								wag_optn = FALSE;
								lg_optn = FALSE;
								fprintf(STDOUT, "(consists very likely of amino acids encoded on cpDNA)\n");
							} else {
								/* if (c7) */
								/* c7 -> LG */
								Dayhf_optn = FALSE;
								Jtt_optn = FALSE;
								mtrev_optn = FALSE;
								cprev_optn = FALSE;
								vtmv_optn = FALSE;
								wag_optn = FALSE;
								lg_optn = TRUE;
								fprintf(STDOUT, "(consists very likely of amino acids encoded on nuclear DNA)\n");
							} /* if c6 else c7 */
						} /* if c5 else c6 */
					} /* if c4 else c5 */
				} /* if c3 else c4 */
			} /* if c2 */
		} /* if c1 */
#endif /* CPREV */

	} else if (data_optn == NUCLEOTIDE) {
		fprintf(STDOUT, "(consists very likely of nucleotides)\n");
	} else {
		fprintf(STDOUT, "(consists very likely of binary state data)\n");
	}
} /* guessmodel */


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

/* parts outsourced to treesort.c/treesort.h */
#if 0 /* TREESORT_H */
#endif /* TREESORT_H */

/* parts outsourced to consensus.c/consensus.h */
#if 0 /* CONSENSUS_H */
#endif /* CONSENSUS_H */


/******************************************************************************/
/* storing and evaluating quartet branching information                       */
/******************************************************************************/

/* general remarks:

	for a quartet with the taxa a, b, c, d there are
	three possible binary trees:
	
		1)  (a,b)-(c,d)	[001]
		2)  (a,c)-(b,d)	[010]
		3)  (a,d)-(b,c)	[100]
	
	For every quartet information about its branching structure is
	stored. With the functions  readquartet  and  writequartet
	this information can be accessed. For every quartet (a,b,c,d)
	with a < b < c < d (taxa) the branching information is encoded
	using 4 bits:
	
	value          8             4             2             1
	        +-------------+-------------+-------------+-------------+
	        |  not used   |   tree 3    |   tree 2    |   tree 1    |
	        +-------------+-------------+-------------+-------------+

	If the branching structure of the taxa corresponds to one of the
	three trees the corresponding bit is set. If the branching structure
	is unclear because two of the three trees have the same maximum
	likelihood value the corresponding two bits are set. If the branching
	structure is completely unknown all the bits are set (the highest
	bit is always cleared because it is not used).

*/

/* allocate memory for quartets */
unsigned char *callocquartets(int taxa)
{
	uli nc, numch;
	unsigned char *qinfo;
	
	/* compute number of quartets */
	Numquartets = (uli) taxa*(taxa-1)*(taxa-2)*(taxa-3)/24;
	if (Numquartets % 2 == 0) { /* even number */
		numch = Numquartets/2;
	} else { /* odd number */
		numch = (Numquartets + 1)/2;
	}
	/* allocate memory */
	qinfo = (unsigned char *) calloc((size_t) numch, sizeof(unsigned char));
	if (qinfo == NULL) maerror("quartetinfo in callocquartets");
	for (nc = 0; nc < numch; nc++) qinfo[nc] = 0;
	return(qinfo);
} /* callocquartets */

/* free quartet memory */
void freequartets()
{	
	free(quartetinfo);
} /* freequartets */

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

/* read quartet info - a < b < c < d */
unsigned char readquartet(int a, int b, int c, int d)
{
	uli qnum;

	if (! ((a < b) && (b < c) && (c < d))) {
		fprintf(stderr, "\n\n\nHALT: PLEASE REPORT ERROR HS2 TO DEVELOPERS! (%d,%d,%d,%d)\n", a,b,c,d);
#		if PARALLEL
			PP_Finalize();
#		endif
   		tp_exit(999, NULL, FALSE, __FILE__, __LINE__, exit_wait_optn);

	}
	qnum = (uli) a
			+ (uli) b*(b-1)/2
			+ (uli) c*(c-1)*(c-2)/6
			+ (uli) d*(d-1)*(d-2)*(d-3)/24;
	if (qnum % 2 == 0) { /* even number */
		/* bits 0 to 3 */
		return (quartetinfo[qnum/2] & (unsigned char) 15);
	} else { /* odd number */
		/* bits 4 to 7 */
		return ((quartetinfo[(qnum-1)/2] & (unsigned char) 240)>>4);
	}
} /* readquartet */

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

/* write quartet info - a < b < c < d, 0 <= info <= 15 */
void writequartet(int a, int b, int c, int d, unsigned char info)
{
	uli qnum;

	qnum = (uli) a
			+ (uli) b*(b-1)/2
			+ (uli) c*(c-1)*(c-2)/6
			+ (uli) d*(d-1)*(d-2)*(d-3)/24;
	if (qnum % 2 == 0) { /* even number */
		/* bits 0 to 3 */
		quartetinfo[qnum/2] =
			((quartetinfo[qnum/2] & (unsigned char) 240) |
			(info & (unsigned char) 15));
	} else { /* odd number */
		/* bits 4 to 7 */
		quartetinfo[(qnum-1)/2] =
			((quartetinfo[(qnum-1)/2] & (unsigned char) 15) |
			((info & (unsigned char) 15)<<4));
	}
} /* writequartet */

/**************/
/* prototypes */
int openfiletowrite(FILE **, char[], char[], FILE *);
void closefile(FILE *);

/* sorts three doubles in descending order */
void sort3doubles(dvector num, ivector order)
{
	if (num[0] > num[1]) {
		if(num[2] > num[0]) {
			order[0] = 2;
			order[1] = 0;
			order[2] = 1;		
		} else if (num[2] < num[1]) {
			order[0] = 0;
			order[1] = 1;
			order[2] = 2;		
		} else {
			order[0] = 0;
			order[1] = 2;
			order[2] = 1;		
		}
	} else {
		if(num[2] > num[1]) {
			order[0] = 2;
			order[1] = 1;
			order[2] = 0;		
		} else if (num[2] < num[0]) {
			order[0] = 1;
			order[1] = 0;
			order[2] = 2;		
		} else {
			order[0] = 1;
			order[1] = 2;
			order[2] = 0;		
		}
	}
} /* sort3doubles */

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

/* compute Bayesian weights from log-lkls d1, d2, d3 */
unsigned char loglkl2weight(int    a,
                            int    b,
                            int    c,
                            int    i,
                            double d1,
                            double d2,
                            double d3,
                            int    usebestq)
{
	double onethird;
	unsigned char treebits[3];
	double templog;
	unsigned char tmpweight;
	double temp;
	double temp1, temp2, temp3;
	unsigned char discreteweight[3];

	double tttqweight[3];
	double tttsqdiff[3];
	int tttsqorder[3];
	int tttqworder[3];

	tttqweight[0] = d1;
	tttqweight[1] = d2;
	tttqweight[2] = d3;

	onethird = 1.0/3.0;
	treebits[0] = (unsigned char) 1;
	treebits[1] = (unsigned char) 2;
	treebits[2] = (unsigned char) 4;
	
	/* sort in descending order */
	sort3doubles(tttqweight, tttqworder);

	if (usebestq) {
		tttsqorder[2] = 2;
		discreteweight[tttsqorder[2]] = treebits[tttqworder[0]];
		if (tttqweight[tttqworder[0]] == tttqweight[tttqworder[1]]) {
			discreteweight[tttsqorder[2]] = discreteweight[tttsqorder[2]] || treebits[tttqworder[1]];
			if (tttqweight[tttqworder[1]] == tttqweight[tttqworder[2]]) {
				discreteweight[tttsqorder[2]] = discreteweight[tttsqorder[2]] || treebits[tttqworder[2]];
				discreteweight[tttsqorder[2]] = 7;
			} 
		}
	} else {

		/* compute Bayesian weights */
		templog = tttqweight[tttqworder[1]]-tttqweight[tttqworder[0]];
		if(templog < -TP_MAX_EXP_DIFF)	/* possible, since 1.0+exp(>36) == 1.0 */
			tttqweight[tttqworder[1]] = 0.0;
		else
			tttqweight[tttqworder[1]] = exp(templog);

	 	templog = tttqweight[tttqworder[2]]-tttqweight[tttqworder[0]];
		if(templog < -TP_MAX_EXP_DIFF)	/* possible, since 1.0+exp(>36) == 1.0 */
			tttqweight[tttqworder[2]] = 0.0;
		else
			tttqweight[tttqworder[2]] = exp(templog);

		tttqweight[tttqworder[0]] = 1.0;

		temp = tttqweight[0] + tttqweight[1] + tttqweight[2];

		tttqweight[0] = tttqweight[0]/temp;
		tttqweight[1] = tttqweight[1]/temp;
		tttqweight[2] = tttqweight[2]/temp;

		/* square deviations */
		temp1 = 1.0 - tttqweight[tttqworder[0]];
		tttsqdiff[0] = temp1 * temp1 +
		tttqweight[tttqworder[1]] * tttqweight[tttqworder[1]] +
		tttqweight[tttqworder[2]] * tttqweight[tttqworder[2]];
		discreteweight[0] = treebits[tttqworder[0]];
 
		temp1 = 0.5 - tttqweight[tttqworder[0]];
		temp2 = 0.5 - tttqweight[tttqworder[1]];
		tttsqdiff[1] = temp1 * temp1 + temp2 * temp2 +
		tttqweight[tttqworder[2]] * tttqweight[tttqworder[2]];
		discreteweight[1] = treebits[tttqworder[0]] + treebits[tttqworder[1]];

		temp1 = onethird - tttqweight[tttqworder[0]];
		temp2 = onethird - tttqweight[tttqworder[1]];
		temp3 = onethird - tttqweight[tttqworder[2]];
		tttsqdiff[2] = temp1 * temp1 + temp2 * temp2 + temp3 * temp3;
		discreteweight[2] = (unsigned char) 7;

		/* sort in descending order */
		sort3doubles(tttsqdiff, tttsqorder);
	}

	tmpweight = discreteweight[tttsqorder[2]];
	return(tmpweight);
} /* loglkl2weight */


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

/* checks out all possible quartets */
void computeallquartets()
{	
#ifdef BININPUTDEBUG
	unsigned long TTT=1;
#endif /* BININPUTDEBUG */

	double onethird;
	uli nq;
	unsigned char treebits[3];
#	if ! PARALLEL
		double brlens[3];
		unsigned char tmpweight;
		FILE *lhfp;
		int a, b, c, i;
#	endif

		onethird = 1.0/3.0;
		treebits[0] = (unsigned char) 1;
		treebits[1] = (unsigned char) 2;
		treebits[2] = (unsigned char) 4;
	
	if (show_optn) { /* list all unresolved quartets */
		openfiletowrite(&unresfp, UNRESOLVED, "unresolved quartet trees", stdinput_fp);
		fprintf(unresfp, "List of all completely unresolved quartets:\n\n");
	}

	nq = 0;
	badqs = 0;
	
	/* start timer - percentage of completed quartets */
	time(&time0);
	time1 = time0;
	mflag = 0;
	
#	if PARALLEL
	{
		schedtype sched; 
		int flag;
		MPI_Status stat;

		int dest = 1;
		uli qaddr  =0;
		uli qamount=0;
		int qblocksent = 0;
		int apr;
		uli Sumquartets;
		uli sq, noq;

		uli tmpfullresqs;
		uli tmppartresqs;
		uli tmpunresqs;
		uli tmpmissingqs;

		Sumquartets = numquarts(Maxspc);
		initsched(&sched, Sumquartets, PP_NumProcs-1, 4);
		qamount=SCHEDALG_ML_STEP(&sched);

		while (qamount > 0) {
			if (PP_emptyslave()) {
				PP_RecvQuartBlock(0, &sq, &noq, quartetinfo, &tmpfullresqs, &tmppartresqs, &tmpunresqs, &tmpmissingqs, &usebestq_optn, &apr);
				qblocksent -= noq;
			}
			dest = PP_getslave();
			PP_SendDoQuartBlock(dest, qaddr, qamount, usebestq_optn, (approxqp_optn ? QUART_APPROX : QUART_EXACT));
			qblocksent += qamount;
			qaddr += qamount;
			qamount=SCHEDALG_ML_STEP(&sched);

			MPI_Iprobe(MPI_ANY_SOURCE, PP_QUARTBLOCKSPECS, PP_Comm, &flag, &stat);
			while (flag) {
				PP_RecvQuartBlock(0, &sq, &noq, quartetinfo, &tmpfullresqs, &tmppartresqs, &tmpunresqs, &tmpmissingqs, &usebestq_optn, &apr);
				qblocksent -= noq;
				MPI_Iprobe(MPI_ANY_SOURCE, PP_QUARTBLOCKSPECS, PP_Comm, &flag, &stat);
			}
			checktime(&time0, &time1, &time2, qaddr - qblocksent, Sumquartets, &mflag);
		}
		while (qblocksent > 0) {
			PP_RecvQuartBlock(0, &sq, &noq, quartetinfo, &tmpfullresqs, &tmppartresqs, &tmpunresqs, &tmpmissingqs, &usebestq_optn, &apr);
			qblocksent -= noq;
			checktime(&time0, &time1, &time2, qaddr - qblocksent, Sumquartets, &mflag);
		}
	}
#	else /* ! PARALLEL */

		addtimes(GENERAL, &tarr);
		if (savequartlh_optn) {
			openfiletowrite(&lhfp, ALLQUARTLH, "all quartet likelihoods", stdinput_fp);
			if (saveqlhbin_optn) {
				if (!saveqlhblen_optn) writetpqfheader(Maxspc, lhfp, 3);
				else                   writetpqfheader(Maxspc, lhfp, 5);
#ifdef BININPUTDEBUG
				{
					double ttt;
					((unsigned char*)(&ttt))[0]=(unsigned char)0xaf;
					((unsigned char*)(&ttt))[1]=(unsigned char)0xfe;
					((unsigned char*)(&ttt))[2]=(unsigned char)0xaf;
					((unsigned char*)(&ttt))[3]=(unsigned char)0xfe;
					((unsigned char*)(&ttt))[4]=(unsigned char)0xab;
					((unsigned char*)(&ttt))[5]=(unsigned char)0xcc;
					((unsigned char*)(&ttt))[6]=(unsigned char)0xab;
					((unsigned char*)(&ttt))[7]=(unsigned char)0xcc;
					fprintf(stdout,"%ld: TTT: ", (long int)-1);
					fprintfdbl2hex(stdout, ttt);
					fprintf(stdout,"\n");
				}
#endif /* BININPUTDEBUG */
			} else {
				if (!saveqlhblen_optn) writetpqfheader(Maxspc, lhfp, 4);
				else                   writetpqfheader(Maxspc, lhfp, 6);
			}
		}

		for (i = 3; i < Maxspc; i++) 
			for (c = 2; c < i; c++) 
				for (b = 1; b < c; b++)
					for (a = 0; a < b; a++) {
						nq++;

						/* generate message every >TIMECHECK_INTERVAL seconds */
						/* check timer */
						checktime(&time0, &time1, &time2, nq, Numquartets, &mflag);

						/* maximum likelihood values */

						/* exact or approximate maximum likelihood values */
						compute_quartlklhds(a,b,c,i,&qweight[0],&qweight[1],&qweight[2], &brlens[0], &brlens[1], &brlens[2], (approxqp_optn ? QUART_APPROX : QUART_EXACT));

						if (savequartlh_optn) {
							if (saveqlhbin_optn) {
								if (saveqlhblen_optn) {
									/* write quartet likelihoods in binary format */
									fwrite(qweight, sizeof(double), 3, lhfp);
									/* write quartets' inner branch lengths in binary format */
									fwrite(brlens, sizeof(double), 3, lhfp);
								} else {
									/* write quartet likelihoods in binary format */
									fwrite(qweight, sizeof(double), 3, lhfp);
#ifdef BININPUTDEBUG
fprintf(stdout,"%ld: TTT: ", TTT++);
/* fprintf(stdout, "%-16.9f= ", qweight[0]); */
fprintfdbl2hex(stdout, qweight[0]);
fprintf(stdout,", ");
/* fprintf(stdout, "%-16.9f= ", qweight[1]); */
fprintfdbl2hex(stdout, qweight[1]);
fprintf(stdout,", ");
/* fprintf(stdout, "%-16.9f= ", qweight[2]); */
fprintfdbl2hex(stdout, qweight[2]);
fprintf(stdout,"\n");
#endif /* BININPUTDEBUG */
								}
							} else {
								if (saveqlhblen_optn) {
									/* write quartet likelihoods + quartets' inner branch lengths in ASCII format */
									fprintf(lhfp, "(%d,%d,%d,%d)\t%f\t%f\t%f\t%f\t%f\t%f\n", a, b, c, i, qweight[0], qweight[1], qweight[2], brlens[0], brlens[1], brlens[2]); 
								} else {
									/* write quartet likelihoods in ASCII format */
									fprintf(lhfp, "(%d,%d,%d,%d)\t%f\t%f\t%f\n", a, b, c, i, qweight[0], qweight[1], qweight[2]); 
	
								}
							}
						}

						tmpweight = loglkl2weight(a, b, c, i, qweight[0], qweight[1], qweight[2], usebestq_optn);
/* tmpweight = loglkl2weight(a, b, c, i, d1, d2, d3, usebestq_optn); */

						/* determine best discrete weight */
						/* writequartet(a, b, c, i, discreteweight[sqorder[2]]); */
						writequartet(a, b, c, i, tmpweight);


						/* compute sums of topologies per taxon, step by step */
						++(qinfomatr[8][a]);
						++(qinfomatr[8][b]);
						++(qinfomatr[8][c]);
						++(qinfomatr[8][i]);
						++(qinfomatr[tmpweight][a]);
						++(qinfomatr[tmpweight][b]);
						++(qinfomatr[tmpweight][c]);
						++(qinfomatr[tmpweight][i]);

						if ((tmpweight <= 2) || (tmpweight == 4)) {
							fullresqs++;
						} else {
							if (tmpweight == 7) {
								unresqs++;
							} else {
								if (tmpweight == 0) {
									missingqs++;
								} else {
									partresqs++;
								}
							}
						}

						/* counting completely unresolved quartets */
						if (tmpweight == 7) {
							badqs++;
							badtaxon[a]++;
							badtaxon[b]++;
							badtaxon[c]++;
							badtaxon[i]++;
							if (show_optn) {
								fputid10(unresfp, a);
								fprintf(unresfp, "  ");
								fputid10(unresfp, b);
								fprintf(unresfp, "  ");
								fputid10(unresfp, c);
								fprintf(unresfp, "  ");
								fputid(unresfp, i);
								fprintf(unresfp, "\n");
							}
						}
						addtimes(QUARTETS, &tarr);
		}
		if (savequartlh_optn) {
			if (saveqlhbin_optn) {
				double ttt;
				((unsigned char*)(&ttt))[0]=(unsigned char)0xaf;
				((unsigned char*)(&ttt))[1]=(unsigned char)0xfe;
				((unsigned char*)(&ttt))[2]=(unsigned char)0xaf;
				((unsigned char*)(&ttt))[3]=(unsigned char)0xfe;
				((unsigned char*)(&ttt))[4]=(unsigned char)0xab;
				((unsigned char*)(&ttt))[5]=(unsigned char)0xcc;
				((unsigned char*)(&ttt))[6]=(unsigned char)0xab;
				((unsigned char*)(&ttt))[7]=(unsigned char)0xcc;
  				fwrite(&ttt, sizeof(double), 1, lhfp);
#ifdef BININPUTDEBUG
				fprintf(stdout,"%ld: TTT: ", (long int)-2);
				fprintfdbl2hex(stdout, ttt);
				fprintf(stdout,"\n");
#endif /* BININPUTDEBUG */
  			}

			closefile(lhfp);
		}
		if (show_optn)
			closefile(unresfp);
		if (mflag == 1)
			fprintf(STDOUT, "\n");
#	endif /* PARALLEL */

} /* computeallquartets */


/* un-trueID-ed (HAS) */
/* check the branching structure between the leaves (not the taxa!)
   A, B, C, and I (A, B, C, I don't need to be ordered). As result,
   the two leaves that are closer related to each other than to leaf I
   are found in chooseX and chooseY while neighb keeps the one related
   to I, if the quartet is resolved. 
   If the branching structure is partly resolved, chooseX contains the 
   leave non-neighbor to I, while chooseY and neighb are set randomly 
   to the two remaining leaves. In the unresolved case all (chooseX,
   chooseY, neighb) are chosen randomly from the possible taxa */
unsigned char checkquartet(
                  int  A,          /* quartet taxon ID                 */
                  int  B,          /* dito.                            */
                  int  C,          /* dito.                            */
                  int  I,          /* dito., to be inserted            */
                  int *chooseX,    /* (chooseX,chooseY | neighb,I)     */
                  int *chooseY,    /* chooseX+Y are non-neighbors of I */
                  int *neighb,     /* neighb is neighbors of I         */
                  int *status,     /* status of the quartet: 0=missing, 1=unresolved, 2=partly, 3=resolved */
                  int rootquartsonly_optn) /* use only root quartets   */
{
	int i, j;              /* counter */
	int a;                 /* temp variable for sorting */
	int leaf[5];           /* sorting array */
	int ipos;              /* position of leaf I in sorted array */
	unsigned char qresult; /* quartet topology */
	unsigned char qinfo;   /* quartet topology */
	int notunique = FALSE; /* non-unique topology */
	int rtemp;


        if (rootquartsonly_optn) { /* use only root quartets   */
		if (! ((A==outgroup) || (B==outgroup) || /* root not in quartet? */
                       (C==outgroup) || (I==outgroup))) {
			*status = QUARTET_MISSING;
			qinfo = 0;

			/* random setup */
			rtemp = randominteger(3);
			switch (rtemp) {
				case 0:
					*neighb  = A;
					*chooseX = B;
					*chooseY = C;
					break;
				case 1:
					*neighb  = B;
					*chooseX = C;
					*chooseY = A;
					break;
				case 2:
					*neighb  = C;
					*chooseX = A;
					*chooseY = B;
					break;
			}
			return(qinfo);
		} /* end if root not in quartet? */
	} /* end if rootquartsonly */

	/* The relationship between leaves and taxa is defined by trueID */
	leaf [1] = A; /* taxon ID of leaf 1 */
	leaf [2] = B; /* taxon ID of leaf 2 */
	leaf [3] = C; /* taxon ID of leaf 3 */
	leaf [4] = I; /* taxon ID of leaf 4 */

	/* sort for taxa */
	/* Source: Numerical Recipes (PIKSR2.C) */
	for (j = 2; j <= 4; j++) {
		a = leaf[j];
		i = j-1;
		while (i > 0 && leaf[i] > a) {
			leaf[i+1] = leaf[i];
			i--;
		}
		leaf[i+1] = a;
	}

	/* where is leaf I ? */
	ipos = 1;
	while (leaf[ipos] != I) ipos++;

	/* look up sequence quartet */
	qresult = readquartet(leaf[1], leaf[2], leaf[3], leaf[4]);
	qinfo = qresult;

	switch (qresult) {

		/*** one single branching structure ***/

		/* 001  (12|34) */
		case 1:		if (ipos == 1 || ipos == 2) {
					*chooseX = leaf[3];
					*chooseY = leaf[4];
					*neighb = (ipos==1 ? leaf[2] : leaf[1]);
				} else {
					*chooseX = leaf[1];
					*chooseY = leaf[2];
					*neighb = (ipos==3 ? leaf[4] : leaf[3]);
				}
				notunique = FALSE;
				*status = QUARTET_FULLY;
				break;

		/* 010  (13|24) */
		case 2:		if (ipos == 1 || ipos == 3) {
					*chooseX = leaf[2];
					*chooseY = leaf[4];
					*neighb = (ipos==1 ? leaf[3] : leaf[1]);
				} else {
					*chooseX = leaf[1];
					*chooseY = leaf[3];
					*neighb = (ipos==2 ? leaf[4] : leaf[2]);
				}
				notunique = FALSE;
				*status = QUARTET_FULLY;
				break;

		/* 100  (23|14) */
		case 4:		if (ipos == 1 || ipos == 4) {
					*chooseX = leaf[2];
					*chooseY = leaf[3];
					*neighb = (ipos==1 ? leaf[4] : leaf[1]);
				} else {
					*chooseX = leaf[1];
					*chooseY = leaf[4];
					*neighb = (ipos==2 ? leaf[3] : leaf[2]);
				}
				notunique = FALSE;
				*status = QUARTET_FULLY;
				break;

		/*** two possible branching structures ***/		

		/* 011  (12|34) + (13|24) */
		case 3:		if (randominteger(2)) qresult = 1;
				else qresult = 2;
				notunique = TRUE;
				*status = QUARTET_PARTLY;
				if (ipos == 1 || ipos == 4) { 
					/* I (=1 or 4) is neighbor to 2 or 3 */
					*chooseX=(ipos==1 ? leaf[4] : leaf[1]);
					if (randominteger(2)) {
						*neighb  = leaf[2]; 
						*chooseY = leaf[3];
					} else {
						*neighb  = leaf[3]; 
						*chooseY = leaf[2];
					}
				} else { 
					/* I (=2 or 3) is neighbor to 1 or 4 */
					*chooseX=(ipos==2 ? leaf[3] : leaf[2]);
					if (randominteger(2)) {
						*neighb  = leaf[1]; 
						*chooseY = leaf[4];
					} else {
						*neighb  = leaf[4]; 
						*chooseY = leaf[1];
					}
				}
				break;

		/* 101  (12|34) + (23|14) */
		case 5:		if (randominteger(2)) qresult = 1;
				else qresult = 4;
				notunique = TRUE;
				*status = QUARTET_PARTLY;
				if (ipos == 1 || ipos == 3) { 
					/* I (=1 or 3) is neighbor to 2 or 4 */
					*chooseX=(ipos==1 ? leaf[3] : leaf[1]);
					if (randominteger(2)) {
						*neighb  = leaf[2]; 
						*chooseY = leaf[4];
					} else {
						*neighb  = leaf[4]; 
						*chooseY = leaf[2];
					}
				} else { 
					/* I (=2 or 4) is neighbor to 1 or 3 */
					*chooseX=(ipos==2 ? leaf[4] : leaf[2]);
					if (randominteger(2)) {
						*neighb  = leaf[1]; 
						*chooseY = leaf[3];
					} else {
						*neighb  = leaf[3]; 
						*chooseY = leaf[1];
					}
				}
				break;

		/* 110  (13|24) + (23|14) */
		case 6:		if (randominteger(2)) qresult = 2;
				else qresult = 4;
				notunique = TRUE;
				*status = QUARTET_PARTLY;

				if (ipos == 3 || ipos == 4) { 
					/* I (=3 or 4) is neighbor to 1 or 2 */
					*chooseX=(ipos==3 ? leaf[4] : leaf[3]);
					if (randominteger(2)) {
						*neighb  = leaf[1]; 
						*chooseY = leaf[2];
					} else {
						*neighb  = leaf[2]; 
						*chooseY = leaf[1];
					}
				} else { 
					/* I (=1 or 2) is neighbor to 3 or 4 */
					*chooseX=(ipos==1 ? leaf[2] : leaf[1]);
					if (randominteger(2)) {
						*neighb  = leaf[3]; 
						*chooseY = leaf[4];
					} else {
						*neighb  = leaf[4]; 
						*chooseY = leaf[3];
					}
				}
				break;

		/*** three possible branching structures ***/

		/* 111  (12|34) + (13|24) + (23|14) */
		case 7:		
				rtemp = randominteger(3);
				switch (rtemp) {
					case 0:
						*neighb  = A;
						*chooseX = B;
						*chooseY = C;
						break;
					case 1:
						*neighb  = B;
						*chooseX = C;
						*chooseY = A;
						break;
					case 2:
						*neighb  = C;
						*chooseX = A;
						*chooseY = B;
						break;
				}
				notunique = TRUE;
				*status = QUARTET_UNRES;
				break;

		/* 0000, 1000, 1111 */  /* empty quartet: missing data -> OK, otherwise ERROR */
		case 0:		
		case 8:		
		case 15:		
				*status = QUARTET_MISSING;
				if (readsubset_optn) {
					rtemp = randominteger(3);
					switch (rtemp) {
						case 0:
							*neighb  = A;
							*chooseX = B;
							*chooseY = C;
							break;
						case 1:
							*neighb  = B;
							*chooseX = C;
							*chooseY = A;
							break;
						case 2:
							*neighb  = C;
							*chooseX = A;
							*chooseY = B;
							break;
					}
					/* *chooseX = -1; */
					/* *chooseY = -1; */
				} else {
#					if PARALLEL
						fprintf(STDOUT, "\n\n\n(%2d)HALT: PLEASE REPORT ERROR K1-PARALLEL TO DEVELOPERS : empty (%d,%d,%d,%d) = %ld (missing quartets are not expected here)\n\n\n", 
							PP_Myid, leaf[1], leaf[2], leaf[3], leaf[4],
							quart2num(leaf[1], leaf[2], leaf[3], leaf[4]));
							PP_Finalize();
   							tp_exit(999, NULL, FALSE, __FILE__, __LINE__, exit_wait_optn);
#					else
						fprintf(STDOUT, "\n\n\nHALT: PLEASE REPORT ERROR K1 TO DEVELOPERS: empty (%d,%d,%d,%d) = %ld (missing quartets are not expected here)\n\n\n", 
							leaf[1], leaf[2], leaf[3], leaf[4],
							quart2num(leaf[1], leaf[2], leaf[3], leaf[4]));
   							tp_exit(999, NULL, FALSE, __FILE__, __LINE__, exit_wait_optn);
#					endif
				}
				notunique = FALSE;
				break;

		default:	/* Program error [checkquartet] */
#				if PARALLEL
					fprintf(STDOUT, "\n\n\n(%2d)HALT: PLEASE REPORT ERROR K2-PARALLEL TO DEVELOPERS (%d,%d,%d,%d) = %ld (%d)\n\n\n", 
						PP_Myid, leaf[1], leaf[2], leaf[3], leaf[4],
						quart2num(leaf[1], leaf[2], leaf[3], leaf[4]), qresult);
						PP_Finalize();
   						tp_exit(999, NULL, FALSE, __FILE__, __LINE__, exit_wait_optn);
#				else
					fprintf(STDOUT, "\n\n\nHALT: PLEASE REPORT ERROR K2 TO DEVELOPERS (%d,%d,%d,%d) = %ld (%d)\n\n\n", 
						leaf[1], leaf[2], leaf[3], leaf[4],
						quart2num(leaf[1], leaf[2], leaf[3], leaf[4]), qresult);
   						tp_exit(999, NULL, FALSE, __FILE__, __LINE__, exit_wait_optn);
#				endif

	} /* switch (qresult) */

	return(qinfo);
} /* checkquartet */


/* un-trueID-ed (HAS) */
/* check the branching structure between the leaves (not the taxa!)
   A, B, C, and I (A, B, C, I don't need to be ordered). As a result,
   the two leaves that are closer related to each other than to leaf I
   are found in chooseX and chooseY. If the branching structure is
   not uniquely defined, chooseX and chooseY are chosen randomly
   from the possible taxa */
unsigned char checkquartet_random(
                  int  A,          /* quartet taxon ID                 */
                  int  B,          /* dito.                            */
                  int  C,          /* dito.                            */
                  int  I,          /* dito., to be inserted            */
                  int *chooseX,    /* (chooseX,chooseY | neighb,I)     */
                  int *chooseY,    /* chooseX+Y are non-neighbors of I */
                  int *neighb,     /* neighb is neighbors of I         */
         int       rootquartsonly_optn) /* in: use root quartets only  */
{
	int i, j;              /* counter */
	int a;                 /* temp variable for sorting */
	int leaf[5];           /* sorting array */
	int ipos;              /* position of leaf I in sorted array */
	unsigned char qresult; /* quartet topology */
	unsigned char qinfo;   /* quartet topology */
	int notunique = FALSE; /* non-unique topology */

	/* The relationship between leaves and taxa is defined by trueID */
	leaf [1] = A; /* taxon ID of leaf 1 */
	leaf [2] = B; /* taxon ID of leaf 2 */
	leaf [3] = C; /* taxon ID of leaf 3 */
	leaf [4] = I; /* taxon ID of leaf 4 */

	/* sort for taxa */
	/* Source: Numerical Recipes (PIKSR2.C) */
	for (j = 2; j <= 4; j++) {
		a = leaf[j];
		i = j-1;
		while (i > 0 && leaf[i] > a) {
			leaf[i+1] = leaf[i];
			i--;
		}
		leaf[i+1] = a;
	}

	/* where is leaf I ? */
	ipos = 1;
	while (leaf[ipos] != I) ipos++;

	/* look at sequence quartet */
	qresult = readquartet(leaf[1], leaf[2], leaf[3], leaf[4]);
	qinfo = qresult;

	/* chooseX and chooseY */
 	/* qresult is not unique, a covered unique qresult is chosen */
	/* and the while loop is circled again (optimizable?) (HAS)  */
	do {
		switch (qresult) {

			/* one single branching structure */

			/* 001 */
			case 1:		if (ipos == 1 || ipos == 2) {
						*chooseX = leaf[3];
						*chooseY = leaf[4];
						*neighb = (ipos==1 ? leaf[2] : leaf[1]);
					} else {
						*chooseX = leaf[1];
						*chooseY = leaf[2];
						*neighb = (ipos==3 ? leaf[4] : leaf[3]);
					}
					notunique = FALSE;
					break;

			/* 010 */
			case 2:		if (ipos == 1 || ipos == 3) {
						*chooseX = leaf[2];
						*chooseY = leaf[4];
						*neighb = (ipos==1 ? leaf[3] : leaf[1]);
					} else {
						*chooseX = leaf[1];
						*chooseY = leaf[3];
						*neighb = (ipos==2 ? leaf[4] : leaf[2]);
					}
					notunique = FALSE;
					break;

			/* 100 */
			case 4:		if (ipos == 1 || ipos == 4) {
						*chooseX = leaf[2];
						*chooseY = leaf[3];
						*neighb = (ipos==1 ? leaf[4] : leaf[1]);
					} else {
						*chooseX = leaf[1];
						*chooseY = leaf[4];
						*neighb = (ipos==2 ? leaf[3] : leaf[2]);
					}
					notunique = FALSE;
					break;

			/* two possible branching structures */		

			/* 011 */
			case 3:		if (randominteger(2)) qresult = 1;
					else qresult = 2;
					notunique = TRUE;
					break;

			/* 101 */
			case 5:		if (randominteger(2)) qresult = 1;
					else qresult = 4;
					notunique = TRUE;
					break;

			/* 110 */
			case 6:		if (randominteger(2)) qresult = 2;
					else qresult = 4;
					notunique = TRUE;
					break;

			/* three possible branching structures */

			/* 111 */
			case 7:		qresult = (1 << randominteger(3)); /* 1, 2, or 4 */
					notunique = TRUE;
					break;

			/* 000 */  /* empty quartet: missing data -> OK, otherwise ERROR */
			case 0:		if (readsubset_optn) {
						*chooseX = -1;
						*chooseY = -1;
					} else {
#						if PARALLEL
							fprintf(STDOUT, "\n\n\n(%2d)HALT: PLEASE REPORT ERROR K1-PARALLEL TO DEVELOPERS : empty (%d,%d,%d,%d) = %ld\n\n\n", 
								PP_Myid, leaf[1], leaf[2], leaf[3], leaf[4],
								quart2num(leaf[1], leaf[2], leaf[3], leaf[4]));
								PP_Finalize();
   								tp_exit(999, NULL, FALSE, __FILE__, __LINE__, exit_wait_optn);
#						else
							fprintf(STDOUT, "\n\n\nHALT: PLEASE REPORT ERROR K1 TO DEVELOPERS: empty (%d,%d,%d,%d) = %ld\n\n\n", 
								leaf[1], leaf[2], leaf[3], leaf[4],
								quart2num(leaf[1], leaf[2], leaf[3], leaf[4]));
   								tp_exit(999, NULL, FALSE, __FILE__, __LINE__, exit_wait_optn);
#						endif
					}
					notunique = FALSE;
					break;

			default:	/* Program error [checkquartet] */
#					if PARALLEL
						fprintf(STDOUT, "\n\n\n(%2d)HALT: PLEASE REPORT ERROR K2-PARALLEL TO DEVELOPERS (%d,%d,%d,%d) = %ld (%d)\n\n\n", 
							PP_Myid, leaf[1], leaf[2], leaf[3], leaf[4],
							quart2num(leaf[1], leaf[2], leaf[3], leaf[4]), qresult);
							PP_Finalize();
   							tp_exit(999, NULL, FALSE, __FILE__, __LINE__, exit_wait_optn);
#					else
						fprintf(STDOUT, "\n\n\nHALT: PLEASE REPORT ERROR K2 TO DEVELOPERS (%d,%d,%d,%d) = %ld (%d)\n\n\n", 
							leaf[1], leaf[2], leaf[3], leaf[4],
							quart2num(leaf[1], leaf[2], leaf[3], leaf[4]), qresult);
   							tp_exit(999, NULL, FALSE, __FILE__, __LINE__, exit_wait_optn);
#					endif
						
		} /* switch (qresult) */
	} while (notunique);

	return(qinfo);
} /* checkquartet_random */