#include <stdio.h>
#include <stdlib.h>
#include <malloc.h>
#include <time.h>
unsigned long*Table,*ColSum,*RowSum,LowBound,UpBound,*TempTable;
unsigned long CumColSum[2],CumRowSum[2];
unsigned long TTable[2][2],TColSum[2],TRowSum[2];
unsigned long Sum,temp,TSum,Size,Total;
int TCol,TRow,Ind_C,Ind_R,Col,Row;
double CHstats,Hstats,temp1;
double Mvol,Hvol,Dvol,CHstats,Hstats;
double Mvol_size,Hvol_size,Dvol_size;
double Weight;
unsigned long generate()
{
unsigned long int i;
unsigned long int intu;
intu=(65539*rand())%Total;
for(i=0;i<1;i++)
{
intu=(65539*intu)%Total;
}
return (unsigned long)((double)intu/Total*(UpBound-LowBound+1)+LowBound);
}
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;
}
//Evaluate the chi-square statistics
double chisq(unsigned long*table)
{
int i,j;
double chis2=0,tt;
for(i=0;i<Row;i++)
for(j=0;j<Col;j++)
{
tt=table[i*Col+j]-(double)RowSum[i]*ColSum[j]/(double)Sum;
chis2+=tt*tt*(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)
{
int i,j;
double ht=0,ht1=0,ht2=0;
for(i=0;i<Col;i++)
ht2+=ColSum[i]*ColSum[i];
for(i=0;i<Row;i++)
{
ht1+=RowSum[i]*RowSum[i];
for(j=0;j<Col;j++)
{
ht+=table[i*Col+j]*table[i*Col+j];
}
}
return ht-ht1*ht2/Sum/Sum;
}
//Build up the tree and Evaluate the statistics
void FillIn()
{
if(TCol<Col-1&&TRow<Row-1)
{
if(TRow==0)
{
CumColSum[0]=0;
CumColSum[1]=0;
}
else
{
for(Ind_R=0;Ind_R<TRow;Ind_R++)
{
if(Ind_R==0){CumColSum[0]=TempTable[TCol];}
else{CumColSum[0]+=TempTable[Ind_R*Col+TCol];}
temp=0;
for(Ind_C=0;Ind_C<=TCol-1;Ind_C++)
temp+=TempTable[Ind_R*Col+Ind_C];
if(Ind_R==0)CumColSum[1]=temp;
else CumColSum[1]+=temp;
if(Ind_R==0)TSum=RowSum[Ind_R];
else TSum+=RowSum[Ind_R];
}
CumColSum[1]=TSum-CumColSum[1];
}
for(Ind_C=TCol;Ind_C<Col;Ind_C++)
{
if(Ind_C==TCol) TColSum[1]=ColSum[Ind_C];
else TColSum[1]+=ColSum[Ind_C];
}
TSum=TColSum[1]-CumColSum[1];
if(TCol==0)
{
CumRowSum[0]=0;
}
else
{
for(Ind_C=0;Ind_C<TCol;Ind_C++)
{
if(Ind_C==0){CumRowSum[0]=TempTable[TRow*Col];}
else{CumRowSum[0]+=TempTable[TRow*Col+Ind_C];}
}
}
TColSum[0]=ColSum[TCol]-CumColSum[0];
TRowSum[0]=RowSum[TRow]-CumRowSum[0];
TColSum[1]=TSum-TColSum[0];
TRowSum[1]=TSum-TRowSum[0];
LowBound=Max(0,TColSum[0]-TRowSum[1]);
UpBound=Min(TColSum[0],TRowSum[0]);
Weight=Weight*(UpBound-LowBound+1);
TempTable[TRow*Col+TCol]=generate();
}
else if(TRow<Row-1)
{
temp=0;
for(Ind_C=0;Ind_C<TCol;Ind_C++)
{
temp+=TempTable[TRow*Col+Ind_C];
}
LowBound=UpBound=RowSum[TRow]-temp;
TempTable[TRow*Col+TCol]=LowBound;
}
else
{
temp=0;
for(Ind_R=0;Ind_R<TRow;Ind_R++)
{
temp+=TempTable[Ind_R*Col+TCol];
}
LowBound=UpBound=ColSum[TCol]-temp;
TempTable[TRow*Col+TCol]=LowBound;
if((TRow==Row-1)&&(TCol==Col-1))
{
temp1=chisq(&TempTable[0]);
/*M-Volumn Test*/
Mvol_size+=Weight;
if(temp1<CHstats)
Mvol+=Weight;
/*D-Volumn Test*/
if((Table[0]-ColSum[0]*RowSum[0]/(double)Sum)*(TempTable[0]-ColSum[0]*RowSum[0]/(double)Sum)>=0)
{
Dvol_size+=Weight;
if(temp1<CHstats)
Dvol+=Weight;
}
/*H-Volumn Test*/
temp1=H_T(&TempTable[0]);
if(temp1*Hstats>=0)
{
Hvol_size+=Weight;
if(sign(temp1)*temp1<sign(Hstats)*Hstats)
Hvol+=Weight;
}
}
}
}
void evaluate()
{
int i,j,k,Limit;
unsigned long l;
/*Evaluate the Chi-square statistic of observed table*/
CHstats=chisq(&Table[0]);
/*Evaluate the H statistics of observed table*/
Hstats=H_T(&Table[0]);
srand((unsigned)time(NULL));
Limit=Min(Row-1,Col-1);
for(l=0;l<Size;l++)
{
for(k=0;k<Limit;k++)
{
for(i=k;i<Col;i++)
{
if(i==0&&k==0)
{
TRowSum[0]=RowSum[0];
TColSum[0]=ColSum[0];
TRowSum[1]=Sum-RowSum[0];
TColSum[1]=Sum-ColSum[0];
LowBound=Max(0,TColSum[0]-TRowSum[1]);
UpBound=Min(TColSum[0],TRowSum[0]);
TRow=0;
TCol=0;
Weight=1;
FillIn();
}
else
{
TRow=k;
TCol=i;
FillIn();
}
}
for(j=k+1;j<Row;j++)
{
TRow=j;
TCol=k;
FillIn();
}
}
if(Row<Col)
{
TRow=Row-1;
for(i=Row-1;i<Col;i++)
{
TCol=i;
FillIn();
}
}
else
{
TCol=Col-1;
for(i=Col-1;i<Row;i++)
{
TRow=i;
FillIn();
}
}
}
Mvol=Mvol/Mvol_size;
Dvol=Dvol/Dvol_size;
Hvol=sign(Hstats)*Hvol/Hvol_size;
}
main()
{
FILE*in,*out;
char Input[256],Output[256];
int i,j,cell,length, flag;
unsigned long control;
/*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);
}
printf("\nHow many Monte Carlo samples?\n");
scanf("%ld",&Size);
fscanf(in,"%d %d %d",&Row,&Col, &cell);
for(i=0;i<Row;i++)
{
for(j=0;j<Col;j++)
{
fscanf(in,"%ld",&control);
}
}
length=0;
do
{
fscanf(in,"%ld",&control);
if(control==0)
length++;
}while(control==0);
length+=Row*Col;
rewind(in);
fprintf(out,"%s","Mvol Hvol\n");
Total=1;
for(i=0;i<31;i++)
Total=Total*2;
while(!feof(in))
{
Mvol=Hvol=Dvol=0;
Mvol_size=Hvol_size=Dvol_size=0;
Sum=0;
/*Enter the dimensions information of the matrix*/
fscanf(in,"%d %d %d",&Row,&Col,&cell);
Table=(unsigned long*)malloc(Row*Col*sizeof(unsigned long));
TempTable=(unsigned long*)malloc(Row*Col*sizeof(unsigned long));
ColSum=(unsigned long*)malloc(Col*sizeof(unsigned long));
RowSum=(unsigned long*)malloc(Row*sizeof(unsigned long));
/*Input the Table*/
for(i=0;i<2;i++)
{
CumRowSum[i]=0;
CumColSum[i]=0;
for(j=0;j<2;j++)
TTable[i][j]=0;
TColSum[i]=0;
TRowSum[i]=0;
}
for(i=0;i<Row;i++)
{
for(j=0;j<Col;j++)
{
TempTable[i*Col+j]=0;
fscanf(in,"%ld",&Table[i*Col+j]);
Sum+=Table[i*Col+j];
if(j==0)RowSum[i]=Table[i*Col+j];
else RowSum[i]+=Table[i*Col+j];
if(i==0)ColSum[j]=Table[i*Col+j];
else ColSum[j]+=Table[i*Col+j];
}
}
for(i=0;i<length-cell;i++)
fscanf(in,"%ld",&control);
flag=0;
for(i=0;i<Row;i++) {
if (RowSum[i] == 0) {
flag=1;
break;
}
}
if (flag == 0) {
for(j=0;j<Col;j++) {
if (ColSum[j] == 0) {
flag=1;
break;
}
}
}
if (flag == 1) {
printf("Error message: Not all row and column sums are
positive. Please remove rows and columns with zero sums.\n");
break;
}
evaluate();
fprintf(out,"%lf %lf\n",Mvol,Hvol);
free(Table);
free(TempTable);
free(ColSum);
free(RowSum);
fscanf(in,"\n");
}
fclose(out);
fclose(in);
}
