#include <stdio.h>
#include <math.h>

#define dim 2

/*Return the minimum value*/

unsigned long max(long x, long y)
{
	if(x<y)return (unsigned long)y;
	else return (unsigned long)x;
}

/*Return the maximum value*/

unsigned long min(unsigned long x,unsigned long y)
{
	if(x<y)return x;
	else return y;
}

/*Chisquare statistics*/

double chisq(unsigned long*table,unsigned long*ColSum,unsigned long*RowSum,unsigned long sum)
{   
	int i,j;
	double chis2=0;
	for(i=0;i<dim;i++)
		for(j=0;j<dim;j++)
		{
	       chis2+=pow(table[i*dim+j]-(double)RowSum[i]*ColSum[j]/(double)sum,2)*(double)sum/(RowSum[i]*ColSum[j]);
		}
	return chis2;
}

/*Sign function*/

int sign(double x)
{
	if(x>0)return 1;
	else return -1;
}

/*Excess of homozygosity*/

double H_T(unsigned long*table,unsigned long*ColSum,unsigned long*RowSum,unsigned long sum)
{
	int i,j;
	double ht=0,ht1=0,ht2=0;
	for(i=0;i<dim;i++)
	{
        ht1+=RowSum[i]*RowSum[i];
		ht2+=ColSum[i]*ColSum[i];
		for(j=0;j<dim;j++)
		{
			ht+=table[i*dim+j]*table[i*dim+j];
		}
	}
	return ht-ht1*ht2/sum/sum;
}

main()
{
	unsigned long Table[dim][dim],ColSum[dim],RowSum[dim],LowBound,UpBound,n;
	unsigned long TTable[dim][dim],sum;
	double Mvol=0,Dvol=0,Hvol=0,CHstats,Hstats,temp;
	double Mvol_size=0,Dvol_size=0,Hvol_size=0;
	int i,j, flag;
	FILE*in,*out;
	char Input[256],Output[256];

	/*Input the Table*/

	printf("\nPlease enter the data file name (including the folders,i.e. c:\\user\\table.txt):\n");
	scanf("%s",Input);
	printf("Please enter the output file name (including the folders,i.e. c:\\user\\output.txt):\n");
	scanf("%s",Output);
	while((in=fopen(Input,"r"))==NULL)
	{
	    printf("\nPlease enter the data file name (including the folders,i.e. c:\\user\\table.txt):\n");
	    scanf("%s",Input);
	}
	while((out=fopen(Output,"w"))==NULL)
	{
        printf("Please enter the output file name (including the folders,i.e. c:\\user\\output.txt):\n");
	    scanf("%s",Output);
	}
	fprintf(out,"%s","Dvol     Mvol     Hvol\n");
	while(!feof(in))
	{		
		Mvol=Dvol=Hvol=0;
		Mvol_size=Dvol_size=Hvol_size=0;
		sum=0;

	    for(i=0;i<dim;i++)
		{
		   for(j=0;j<dim;j++)
		   {			
			fscanf(in,"%ld",&Table[i][j]);

			sum+=Table[i][j];
			if(j==0)RowSum[i]=Table[i][j];
			else RowSum[i]+=Table[i][j];
			if(i==0)ColSum[j]=Table[i][j];
			else ColSum[j]+=Table[i][j];
		   }		   
		}

	    flag=0;
	    for(i=0;i<dim;i++) {
	      if (RowSum[i] == 0) {
		flag=1;
		break;
	      }
	    }
	    if (flag == 0) {
	      for(j=0;j<dim;j++) {
		if (ColSum[j] == 0) {
		  flag=1;
		  break;
		}
	      }
	    }

	    if (flag == 1) {
	      printf("Error message: Now all row and column sums are positive. Please remove rows and columns with zero sums.\n");
	       break;
	    }          

	    LowBound=max(0,ColSum[0]-RowSum[1]);
	    UpBound=min(ColSum[0],RowSum[0]);

	 /*Evaluate the Chi-square statistic of observed table*/
	    CHstats=chisq(&Table[0][0],&ColSum[0],&RowSum[0],sum);

	/*Evaluate the H statistics of observed table*/
	    Hstats=H_T(&Table[0][0],&ColSum[0],&RowSum[0],sum);

	    for(n=LowBound;n<=UpBound;n++)
		{
		   for(i=0;i<dim;i++)
			for(j=0;j<dim;j++)
			{
				if(i==0&&j==0)TTable[i][j]=n;
				else if(i==0)TTable[i][j]=RowSum[i]-n;
				else if(i==1)TTable[i][j]=ColSum[j]-TTable[i-1][j];
			}

		temp=chisq(&TTable[0][0],&ColSum[0],&RowSum[0],sum);

		/*M-Volumn Test*/
            if(temp<CHstats)
			  Mvol+=1;

		/*D-Volumn Test*/
		    if((Table[0][0]-ColSum[0]*RowSum[0]/(double)sum)*(TTable[0][0]-ColSum[0]*RowSum[0]/(double)sum)>=0)
			{
			  Dvol_size+=1;
			  if(temp<CHstats)
				Dvol+=1;
			}

        /*H-Volumn Test*/
        temp=H_T(&TTable[0][0],&ColSum[0],&RowSum[0],sum);
		if(temp*Hstats>=0)
		{
			Hvol_size+=1;
			if(sign(temp)*temp<sign(Hstats)*Hstats)
				Hvol+=1;
		}
	}
	Mvol_size=UpBound-LowBound+1;
	Mvol=Mvol/Mvol_size;
	Dvol=Dvol/Dvol_size;
	Hvol=sign(Hstats)*Hvol/Hvol_size;
	fprintf(out,"%lf %lf %lf\n",Dvol,Mvol,Hvol);
	fscanf(in,"\n");
	}
	fclose(out);
	fclose(in);
}