/*  This program computes the test statistic W on a data set:                      
    W=N5/N -(N3+N5)(N2+N5)/N^2  . Sites are divided into 5 types.
    N is the total number of sites with no gap or missing data.
    Needs
    - the number of sequences in each clade, or a additional header in the "phylip" datafile,
    specifying the clades names and the (abbreviated) names of sequences in those clades.
    - a Phylip Data file containing the sequences
    Result: W is written on the standard output
    Actually, as the program is now, it could take an invariable proportion of sites as
    an optional argument. It wouldn't be used.
*/


#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <ctype.h>
#include <time.h>
#include <math.h>

#ifdef __MWERKS__
#ifdef __DROP_CONSOLE__
#include "DropConsole.h"
#else
#include <console.h>
#endif
#endif


#include "global.h"

/*#define MAX_NAME_LEN 64   /* or #include "tree.h"  */
#define MAX_NAME_LEN 40 
#define MAX_NUM_SEQ 256

int printHelp,verbose;
char cladeFileName[256];
char treeFileName[256];
double proportionInvariable; /* argument of the program */
char **names;
char **sequences; 
int numseq[2];   /* number of sequences in clade 1 and clade 2. */
int numSequences, numSites;
int *clade1;     /* numbers of the clade1's sequences  */
int *clade2;
char ***cladeNames;

char cladeInfo; /* ='n' in case option -n, 'c' in case option -c and '0' in case no option (yet) */
FILE *clade_fv;
FILE *tree_fv;


static void PrintUsage()
{ 
	fprintf(stderr, "Usage: test-stat [-h][-q][-n # # | -c filename | -t filename] < datafile\n");
	fprintf(stderr, "  -h: Give this help message\n");
	fprintf(stderr, "  -v: verbose\n");
	fprintf(stderr, "  -n: #1 #2 = number of sites in the first and second clade.\n");
	fprintf(stderr, "      Sequences from the first clade are supposed to come first,\n");
	fprintf(stderr, "      and sequences from the second clade are supposed to come next.\n");
	fprintf(stderr, "  -c: the clade names and the corresponding sequence names\n");
	fprintf(stderr, "      (possibly abbreviated) are to be read in a file\n");
	fprintf(stderr, "  -t: the clades definition are to be read in a tree file,\n");
	fprintf(stderr, "      a clade is represented by a internal node. \n");
	fprintf(stderr, "      (not implemented yet) \n");
	fprintf(stderr, "  datafile: name of sequence data file \n");
}

void ReadParams(int argc, char **argv)
{
	char *P, ch;
	int i;
	
	cladeInfo='0';
	printHelp=0;
	verbose=0;

	for (i=1; i<argc; i++) {
		P=argv[i];

		/*		if (*P!='-') {
		  if (cladeFile) { 
		    fprintf(stderr, "Illegal command parameter: %s\n\n", argv[i]); 
		    PrintUsage(); 
		    exit(0); 
		  } 
		  strcpy(cladeFileName, argv[i]);  
		  cladeFile=1; 
		} 
		else  { */
		if (*P=='-' && toupper(*(P+1))=='C') {
		  if ((cladeInfo=='n')||(cladeInfo=='t')){
		    fprintf(stderr, "you can only use one of the options -n, -c and -t\n");
		    PrintUsage();
		    exit(0);
		  }
		  cladeInfo='c';
		  P=argv[i] + 2;
		  if (*P=='\0') {
		    i++;
		    P = argv[i];
		  }
		  strcpy(cladeFileName, P);
		} else {		  
		  while (*P) {
		    *P=toupper(*P); 
		    P++; 
		  }

		  P=argv[i];
		  P++;
		  ch=*P;
		  P++;
			switch (ch) {
				case 'H':
					printHelp=1;
				break;
				case 'V':
					verbose=1;
				break;
				case 'N':
				  if ((cladeInfo=='c')||(cladeInfo=='t')){
				    fprintf(stderr, "you can only use one of the options -n, -c and -t\n");
				    PrintUsage();
				    exit(0);
				  }
				  cladeInfo='n';
				  if (GetIntParams(argc, argv, &i, P, 2, numseq) || numseq[0] < 0 || numseq[1] < 0) {
				    fprintf(stderr, "Bad number of sequences: %d, %d\n\n", numseq[0], numseq[1]);
				    exit(0);
				  }
				break;
				default:
					fprintf(stderr, "Illegal command parameter: %s\n\n", argv[i]);
					PrintUsage();
					exit(0);
				break;
			}
		}
	}
}

void ReadFileParams()
{
/* reads the number of sequences and the number of sites from the data file */

	char ch, st[256];
	

	ch=fgetc(stdin);
	while (!feof(stdin) && isspace(ch)) 
		ch=fgetc(stdin);
		
	ungetc(ch, stdin);

	if (ch!='(' && isdigit(ch)) {
		fgets(st, 255, stdin);
		if ( sscanf( st, " %d %d", &numSequences, &numSites)!=2 ) {
			fprintf(stderr, "Unable to read parameters from standard input\n");
			exit(0);
		}
	}
}

void AllocateMemory()
{
	int i,j;
	
	names=(char **)AllocMem(sizeof(char *)*numSequences, "names", "AllocateMemory", 0);
	sequences=(char **)AllocMem(sizeof(char *)*numSequences, "sequences", "AllocateMemory", 0);

	for (i=0; i<numSequences; i++) {
		names[i]=(char *)AllocMem(sizeof(char)*(MAX_NAME_LEN+1), "names[]", "AllocateMemory", 0);
		sequences[i]=(char *)AllocMem(sizeof(char)*numSites, "sequences[]", "AllocateMemory", 0);
	}
	clade1=(int *)AllocMem(sizeof(int)*numSequences, "clade1", "AllocateMemory",0);
	clade2=(int *)AllocMem(sizeof(int)*numSequences, "clade2", "AllocateMemory",0);
	cladeNames =(char ***)AllocMem(sizeof(char **)*2, "cladeNames", "AllocateMemory", 0);
	for (i=0; i<2; i++){
	  cladeNames[i] =(char **)AllocMem(sizeof(char *)*(MAX_NUM_SEQ+1), "cladeNames[]", "AllocateMemory", 0);
	  for (j=0; j<numSequences+1; j++) {
		cladeNames[i][j]=(char *)AllocMem(sizeof(char)*(MAX_NAME_LEN+1), "cladeNames[][]", "AllocateMemory", 0);
	  }
	}
}

void ReadFile()
{
/* reads the names of the taxa and the sequences. */ 
	int n, b, i;
	char ch;
		
	n=0;
	do {
		ch=fgetc(stdin);
		while ( !feof(stdin) && isspace(ch)) 
			ch=fgetc(stdin);
			
		if ( feof(stdin) ) {
			fprintf(stderr, "Unexpected end of file on standard input\n"); 
			exit(0);
		}
	
		i=0;
		while ( i<MAX_NAME_LEN && !feof(stdin) && !isspace(ch) ) {
			names[n][i]=ch;
			ch=fgetc(stdin);
			i++;
		}
		names[n][i]='\0';
		if (i==0) {
			fprintf(stderr, "Name missing for species %d\n", n+1);
			exit(0);
		}
		while (!feof(stdin) && isspace(ch))
			ch=fgetc(stdin);
		
		if ( feof(stdin) ) {
			fprintf(stderr, "Unexpected end of file on standard input\n");
			exit(0);
		}
		
		b=0;
		while ( !feof(stdin) && b<numSites) {
			if ( !isspace(ch) ) {
				sequences[n][b]=ch;
				b++;
			}
			ch=toupper(fgetc(stdin));
		}
		
		if ( b<numSites ) {
			fprintf(stderr, "Unexpected end of file on standard input\n");
			exit(0);
		}
		
		if (verbose) fprintf(stderr, "%d: %s, bases read: %d\n", n+1, names[n], b);
		n++;
		
		if ( n<numSequences && feof(stdin) ) {
			fprintf(stderr, "Too few sequences in input file\n");
			exit(0);
		}
	} while ( n<numSequences );



	/* it should also read how many taxa are in clade1, how many in clade 2 : */
	/* initialize numseq[0] and [1], and initialize the sizes of *clade1 and *clade2 */
	/* It should read whose taxa are in clade1 and whose are in clade 2 : enter it in *clade1 and 2 */
 
}

int ReadCladeFile()
{
  /* the sequence names are known at this point. 
     this function reads a file with the following format:
     clade1name name1 name2 ...
     clade2name name7 name8 ...
     then each sequence whose name starts with name1 or name 2 or... is put in clade 1
     and each sequence whose name starts with name7 or name 8 or... is put in clade 2
     More precisely clades'sizes numseq[0] and numseq[1] are computed
     and the arrays clade1 and clade2 are filled with the numbers of the corresponding sequences.
  */
	int n, cl, i,s,i0,i1,diff;
	int numSeqClade[2];
	char ch;
	char *p;
	int alreadyfound,found;
		
	if ( (clade_fv=fopen(cladeFileName, "rt"))==NULL ) {
	  fprintf(stderr, "Error opening clade file: '%s'\n", cladeFileName);
	  exit(0);
	}

	ch=fgetc(clade_fv);
	for (cl=0; cl<2; cl++){
	  n=0;
	  i=0;
	  while ( !feof(clade_fv) && isspace(ch) ) { ch=fgetc(clade_fv);}
	  while ( !feof(clade_fv) && ((ch)!='\n') ) {
	    if (isspace(ch)) {
	      cladeNames[cl][n][i]='\0';
	      while ( !feof(clade_fv) && isspace(ch) && ((ch)!='\n')) {  ch=fgetc(clade_fv);}
	      if (!(isspace(ch))){
		n++;
		i=0;
	      }
	    }
	    if (!(isspace(ch))){
	      cladeNames[cl][n][i]=ch;
	      ch=fgetc(clade_fv);
	      i++;
	    }
	  }
	  if (n==0){
	    fprintf(stderr, "missing sequence names in clade file. I won't go further!\n");
	    exit(0);
	  }
	  numSeqClade[cl]=n;
	}
	/* now the names in each clade are known 
	 we still need to get the numbers of the sequences in each clade */

	i0=i1=0;
	for (s=0; s<numSequences; s++){
	  alreadyfound=0;
	  for (cl=0; cl<2; cl++){
	    found=0;
	    n=1;
	    while ((n<= numSeqClade[cl]) && (found==0)){
	      diff=0;
	      i=0;
	      while ((cladeNames[cl][n][i])&& !(diff)){
		if (names[s][i] != cladeNames[cl][n][i] ) diff=1;
		i++;
	      }
	      if (diff==0){/* check that names[s] is not a substring of cladeNames[cl][n]*/
		if (names[s][i] == cladeNames[cl][n][i]){ found=1;}
	      }
	      n++;
	    }
	    if (found){
	      if (cl==0){ 
		clade1[i0]=s;
		i0++;
		alreadyfound=1;
	      }
	      else{ 
		clade2[i1]=s;
		i1++;
		if (alreadyfound){
		  fprintf(stderr, "taxon %s cannot belong to both clades. I won't go further!\n",names[s]);
		  exit(0);
		}
	      }
	    }
	  }
	}
	numseq[0]=i0;
	numseq[1]=i1;
	for (cl=0; cl<2; cl++){
	  if (numseq[cl]==0){
	    fprintf(stderr, "no sequence in clade %d. I won't go further!\n",cl+1);
	    exit(0);
	  }
	}
	if (verbose) { for (cl=0; cl<2; cl++){
	  fprintf(stderr, "name of clade %d:\t",cl+1);
	  for (p=cladeNames[cl][0]; (*p != 0); p++) fprintf(stderr, "%c",*p);
	  fprintf(stderr, "\n");
	  fprintf(stderr, "there are %d sequences in this clade, whose numbers are:\n",numseq[cl]);
	  for (i=0; i<numseq[cl]; i++){
	    fprintf(stderr, "%d\t",(((cl==0)?clade1[i]:clade2[i])+1) );
	  }
	  fprintf(stderr, "\n");
	}}

}

		

int main(int argc, char **argv)
{
	int i, j, k;
	char ch;
	char **P;
	int flag0, flag1;
	int gap;

	int N3,N4,N5;    /* number of sites of type 3,4,5  */
	int N;        /* number of sites with no gap(s)  */
	double W;        /* Test statistic */


#ifdef __MWERKS__
#ifndef __DROP_CONSOLE__
	argc = ccommand(&argv);
#endif
#endif

	ReadParams(argc, argv);

	if (printHelp) {
		PrintUsage();
		exit(0);
	}


	ReadFileParams();
	AllocateMemory();	
	ReadFile();

	switch (cladeInfo){
	case '0':
	  fprintf(stderr, "No information available about clades\n");
	  exit(0);
	  break;
	case 'c':
	  ReadCladeFile();
	  break;
	case 'n':
	  for (i=0; i<numseq[0]; i++) clade1[i]=i;
	  for (i=0; i<numseq[1]; i++) clade2[i]=numseq[0]+i;
	  break;
	}

	if ( numseq[0]+numseq[1] > numSequences ) {
	  fprintf(stderr, "Bad number of sequences, compared to the clades' sizes\n");
	  exit(0);
	}


	N3=N4=N5=0;
	N=0;

	P=sequences;
	for (i=0; i<numSites; i++) {

	  gap=0;
	  for (j=0; j<numseq[0]; j++){        /* see how many gaps there are in this site */
	    if((*P[clade1[j]])=='-') gap++;
	  }
	  for (j=0; j<numseq[1]; j++){ 
	    if(*P[clade2[j]]=='-') gap++;
	  }
	  if ((verbose) && (gap!=0)) fprintf(stderr, "number of gaps for site %d: %d\n",i+1,gap);
	  if (gap==0){
	    N++;

	    flag0=0;
	    j=1;
	    do {                                /*see if the character varies across clade 1 */
	      if(*P[clade1[j-1]]!=*P[clade1[j]]) flag0=1;
	      j++;
	    }
	    while(j<numseq[0] && flag0==0);

	    flag1=0;
	    j=1;
	    do {                                /*see if the character varies across clade 2 */
	      if(*P[clade2[j-1]]!=*P[clade2[j]]) flag1=1;
	      j++;
	    }
	    while(j<numseq[1] && flag1==0);

	    if (flag0==1 && flag1==0) N3++;
	    if (flag0==0 && flag1==1) N4++;
	    if (flag0==1 && flag1==1) N5++;
	  }

	  for (j=0; j<numSequences; j++) P[j]++;  /* we go to the following site */
	}

	W=0.0;
	W += N5*N-(N3+N5)*(N4+N5);
	W /= (N*N);
	if (verbose){
	  fprintf(stderr, "\nnumber of sites with no gap: %d \nsites with gap(s) are ignored.\n",N);
	  fprintf(stderr, "number of sites of type 3: %d , type 4: %d  , type 5: %d\n",N3,N4,N5);
	}
	fprintf(stdout,"%2.8e\n", W);
	return(0);
}