#include<stdio.h>
#include<math.h>
#define SEQLENGTH 5000
#define NUMGENOMES 50
#define DEBUG 0

unsigned _stklen = 2048576;

double mle(int N,int na,int nb ,int nab, double rho,double alpha)
{
	double d;
	int invar = ((double)N) * rho;
	double p;
	double r=1,a=1,b=1;
	int m;
	double dN,l;
	N = N-invar;
	dN = (double)N;
	na = na - invar;
	nb = nb - invar;
	nab = nab - invar;
	
	
	if(alpha == 0)
	{
		if(((double)(na*nb))/(N*nab)<1)
		{
			p = ((double)(na*nb))/(nab*nab);

		}
		else
		{
			p = ((double)(N*N))/(na*nb);
		}
		d = log(p);
	}
	else
	{
		if(((double)(na*nb))/(N*nab)<1)
		{
			p = 2*pow(dN/nab,1/alpha) - ( pow(dN/na,1/alpha) + pow(dN/nb,1/alpha) );
			d = p*alpha;

		}
		else
		{
			 pow(dN/na,1/alpha) + pow(dN/nb,1/alpha) - 2;
			 d = p*alpha;
		}
		
	}

	return d;
	
}



int main(int argc, char *argv[])
{
	int n,N;
	int  data[SEQLENGTH][NUMGENOMES];
	char names[NUMGENOMES][11];
	char inmsg[10];
	double matrix[NUMGENOMES][NUMGENOMES];
	double rho=0,gamma=0;
	char c;
	int na[100];
	int nab[NUMGENOMES][NUMGENOMES];
	FILE *f,*g;
	int i,j,i1,i2;
	int m,reps;

	if(argc!=3)
	{
		printf("usage: cghdist <infile> <outfile>\n");
		exit(1);
	}
	
	f = fopen(argv[1],"r");
	if(f==NULL)
	{
		printf("Couldn't open `%s'\n",argv[1]);
		exit(1);
	}
	g = fopen(argv[2],"r");
	if(g!=NULL)
	{
		while(1)
		{
			printf("%s already exists.  Overwrite? y/n?\n",argv[2]);
			c = getchar();
			if(c=='y')
			{
				break;
			}
			if(c=='n') exit(1);	
		}
	}
	fclose(g);
	
	printf("Rate variation parameter alpha? (0 for fixed rate)\n");
	scanf("%lf",&gamma);
	printf("Proportion of genes invariant?\n");
	scanf("%lf",&rho);
	printf("Number of datasets?\n");
	scanf("%d",&m);
	
	/* Read data */
	
	
	g = fopen(argv[2],"w");
	
	if(g==NULL)
	{
		printf("Couldn't open 'outfile'\n");
		exit(1);
	}
	
	for(reps = 0;reps<m;reps++)
	{
	
	fscanf(f,"%d",&n);
	fscanf(f,"%d",&N);
	
	
	printf("Dataset %d:%d sequences of %d genes\n",reps+1,n,N);
	
	for(i=0;i<n+1;i++)
	{
		na[i] = 0;
		for(j=0;j<n+1;j++) nab[i][j] = 0;
	}

	for(i=0;i<n;i++)
	{
		fscanf(f,"\n");
		for(j=0;j<10;j++) {c = getc(f);names[i][j]=c;}
		names[i][10]='\0';
		/*printf("Read %s\n",names[i]);*/
		j=0;
		while(j<N)
		{
			c = getc(f);
			if(c=='1')			
			{
				data[j][i] = 1;
				j++;
				na[i]++;
			}
			else if(c=='0')
			{
				data[j][i] = 0;
				j++;
			}
		}
	}
			
	for(j=0;j<N;j++) data[j][n]=1;
	na[n]=N;
	strcpy(names[n],"Reference ");
	n++;
				
	/* Make dm */
	fprintf(g,"%d",n);
	for(i1=0;i1<n;i1++)
	{
		fprintf(g,"\n%s",names[i1]);
		for(i2=0;i2<n;i2++)
		{
			for(j=0;j<N;j++)
			{
				
				if(data[j][i1]==1 && data[j][i2]==1) nab[i1][i2]++;
			}
			matrix[i1][i2] = mle(N,na[i1],na[i2],nab[i1][i2],rho,gamma);
			fprintf(g,"%g ",matrix[i1][i2]);
		}
	}
	fprintf(g,"\n");
	
	}
	
	fclose(f);
	fclose(g);
}