#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <malloc.h>
#include <time.h>

static void ReadInfo(int *pNumRow, int *pNumCol, int *pM, int *pN);
static void ReadTable (int NumRow, int NumCol, int **table, int *rowSum, int *colSum, int **s0, int *s0RowSum);
static double SSquareBar(int **a, int NumRow, int NumCol);
static void  InitialSetting (double *pSampleMean, double *pSampleVar, 
			     int *pNumBad, double *pPValue);
static void GenerateATable (int **t, int **s0, int NumRow, int NumCol, 
			   int *rowSum, int *colSum, int *s0RowSum, 
			   double *pCurrentSample, int *pBadSample);
static void UpdatePValue (int badSample, double currentSample, int **t, 
			  int NumRow, int NumCol, double sSquare0Bar, 
			  double *pSampleMean, 
		          double *pSampleVar, double *pPValue, int *pNumBad);
static void PrintOutput (FILE *outFile, int N, int numBad, double sampleMean, double sampleVar, 
			 double pValue);
static double  NChooseM(int n, int m);
static double Mean(double x[], int n);
static int RandomSample2(double weight[], int n);
static double Drafting(int w[], int start, int end, int n, int a[],
		       int s[], int NumRow);
static int ChooseOneInteger(double w[], int start, int end, int NumRow);
static double R(int k, double w[], int start, int end, int NumRow);
static double T(int i, double w[], int start, int end);
static double RandomReal(double low, double high);
static int RandomInteger(int low, int high);

main()
{
	int i, NumRow, NumCol, M, N, numEstimate, numSample, numBad, badSample;
	double sampleMean, currentSample, cvsquare1, sampleVar;
        double sSquare0Bar,  pValue;
	int **table, **t, **s0,  *s0RowSum, *rowSum, *colSum;
	FILE *outFile;

	ReadInfo(&NumRow, &NumCol, &M, &N);

	//allocate space for pointers
	table=malloc(NumRow*sizeof(int *));
	t=malloc(NumRow*sizeof(int *));
	s0=malloc(NumRow*sizeof(int *));
	rowSum=malloc(NumRow*sizeof(int));
	colSum=malloc(NumCol*sizeof(int));
	s0RowSum=malloc(NumRow*sizeof(int));
        for (i=0; i<NumRow; i++) {  
	  table[i]=malloc(NumCol*sizeof(int));
	  t[i]=malloc(NumCol*sizeof(int));
	  s0[i]=malloc(NumCol*sizeof(int));
        }	
	//end of allocation

	srand((int) time(NULL));

	ReadTable(NumRow, NumCol, table, rowSum, colSum, s0, s0RowSum);

	sSquare0Bar=SSquareBar(table, NumRow, NumCol);

	outFile=fopen("result.txt", "w");
	if (outFile != NULL) {	

	  for (numEstimate=0; numEstimate<M; numEstimate++) {
	    InitialSetting(&sampleMean, &sampleVar, &numBad, &pValue);

	    for (numSample=0; numSample<N; numSample++) {
	      GenerateATable(t, s0, NumRow, NumCol, rowSum, colSum, 
					 s0RowSum, &currentSample, &badSample);

	      UpdatePValue(badSample, currentSample, t, NumRow, NumCol, 
		       sSquare0Bar, &sampleMean, &sampleVar, &pValue, &numBad);
	    }	    
	    PrintOutput(outFile, N, numBad, sampleMean, sampleVar, pValue);
	  }
	}
	else { printf("File creating failed!\n");}
	fclose(outFile);
}

static void ReadInfo (int *pNumRow, int *pNumCol, int *pM, int *pN)
{
        int NumRow, NumCol, M, N;

	printf("\nPlease enter the number of rows in the table?\n");
	scanf("%d", &NumRow);
	printf("\nPlease enter the number of columns in the table?\n");
	scanf("%d", &NumCol);
	printf("\nHow many estimates do you need?\n");
	scanf("%d", &M);
	printf("\nHow many samples should each estimate be based on?\n");
	scanf("%d", &N);	

	*pNumRow=NumRow;
	*pNumCol=NumCol;
	*pM=M;
	*pN=N;
}

static void ReadTable (int NumRow, int NumCol, int **table, int *rowSum, int *colSum, int **s0, int *s0RowSum)
{
        int i, j, k, m, n, temp;
        int *maxRange, *minRange;
	FILE *inFile;
	char nameInFile[256];

        maxRange=malloc(NumRow*sizeof(int));
        minRange=malloc(NumRow*sizeof(int));

	printf("\nPlease enter the data file name:\n");
	scanf("%s", nameInFile);
	while((inFile=fopen(nameInFile, "r"))==NULL)
	  {
	    printf("\nCan't read file. Please enter the data file name:\n");
	    scanf("%s", nameInFile);
	  }

        for (i=0; i<NumRow; i++) {  
          for (j=0; j<NumCol; j++) {
            fscanf(inFile, "%d", &table[i][j]);
          }
        }	

	for (i=0; i<NumRow; i++) { 
	  rowSum[i]=0;
	  for (j=0; j<NumCol; j++) {
	    rowSum[i]=rowSum[i]+table[i][j];
	  }
	}

	for (j=0; j<NumCol; j++) { 
	  colSum[j]=0;
	  for (i=0; i<NumRow; i++) {
	    colSum[j]=colSum[j]+table[i][j];
	  }
	}

	//change order of columns
	for (m=0; m<NumCol; m++) {
	  k=m;
	  for (n=m+1; n<NumCol; n++) {
	    if (colSum[n]>colSum[k])
	      k=n;
	  }                   
	  temp=colSum[m];
	  colSum[m]=colSum[k];
	  colSum[k]=temp;

	  for (i=0; i<NumRow; i++) {
	    temp=table[i][m];
	    table[i][m]=table[i][k];	
	    table[i][k]=temp;
	  }
	}                      	

	for (i=0; i<NumRow; i++) { 
	  maxRange[i]=-1;
	  minRange[i]=NumRow+1;
	  for (j=0; j<NumCol; j++) {
	    if (table[i][j] == 1) {
	      if (colSum[j]>maxRange[i]) maxRange[i]=colSum[j];
	      if (colSum[j]<minRange[i]) minRange[i]=colSum[j];
	    }
	  }
	}

	for (i=0; i<NumRow; i++) { 
	  s0RowSum[i]=0;
	  for (j=0; j<NumCol; j++) {
	    if (colSum[j]<minRange[i] || colSum[j]>maxRange[i]) {
	      s0[i][j]=1;
       	      s0RowSum[i]++;
	    }
	    else 
	      s0[i][j]=0;
	  }
	}

	free(maxRange);
	free(minRange);
}

static double SSquareBar(int **a, int NumRow, int NumCol)
{

	int i,j,k;
	double sum;
	int s;

	sum=0;

	for (i=0; i<NumRow; i++) {
		for (j=0; j<NumRow; j++) {
			s=0;
			for (k=0; k<NumCol; k++) {
				s=s+a[i][k]*a[j][k];
			}
			if (j!=i) {
				sum=sum+s*s;
			}
		}
	}
	return (sum/NumRow/(NumRow-1));
}

static void  InitialSetting (double *pSampleMean, double *pSampleVar, 
			     int *pNumBad, double *pPValue)
{    
	    *pSampleMean=0;
	    *pSampleVar=0;
	    *pNumBad=0;
	    *pPValue=0;
}

static void GenerateATable (int **t, int **s0, int NumRow, int NumCol, 
			   int *rowSum, int *colSum,
			  int *s0RowSum, double *pCurrentSample, int *pBadSample)
{
  int j, k, m, badSample, count;
  int *currentRowSum, *currentS0RowSum, *index, *rowSegment, *leftoverCol;
  int *oneIndex;
  int currentColSum, numPosition, temp;
  double currentSample, prob;

  	currentRowSum=malloc(NumRow*sizeof(int));
	currentS0RowSum=malloc(NumRow*sizeof(int));
        index=malloc(NumRow*sizeof(int));
	rowSegment=malloc(NumRow*sizeof(int));
	leftoverCol=malloc(NumRow*sizeof(int));
	oneIndex=malloc(NumRow*sizeof(int));

	      for (j=0;j<NumRow; j++) {
		currentS0RowSum[j]=s0RowSum[j];
		currentRowSum[j]=rowSum[j];
		for (m=0; m<NumCol; m++) {
		  t[j][m]=0;
		}
	      }

	      currentSample=1.0;	      
	      badSample=0;

	      for (j=0; j<NumCol; j++) {
		if (colSum[j]>0 && badSample==0) {
		  count=0;
		  currentColSum=colSum[j];

		  for (m=0; m<NumRow; m++) {
		    numPosition=NumCol-j-currentS0RowSum[m];
		    if (numPosition<currentRowSum[m]) {
		      badSample=1;
		      break;
		    }
		    else if (s0[m][j]==1 || currentRowSum[m]==0)
		      t[m][j]=0;
		    else if (currentRowSum[m]>0 &&
			     currentRowSum[m]==numPosition) {
		      t[m][j]=1;
		      currentColSum--;
		    }
		    else {
		      index[count]=m;
		      rowSegment[count]=currentRowSum[m];
      		      leftoverCol[count]=numPosition;
		      count++;
		    }
		  }

		  if (currentColSum<0 || currentColSum>count) {
		    badSample=1;
		  }

		  if (badSample==1) break;

		  if (count>0) {
		    prob=Drafting(rowSegment,0,count-1,currentColSum,oneIndex, leftoverCol, NumRow);
		    currentSample=currentSample/prob;

		    for (k=0; k<currentColSum; k++) {
		      t[index[oneIndex[k]]][j]=1;
		    }
		  }

		  temp=0;
		  for (m=0; m<NumRow; m++) {
		    currentRowSum[m]=currentRowSum[m]-t[m][j];
		    currentS0RowSum[m]=currentS0RowSum[m]-s0[m][j];
		    temp=temp+t[m][j];
		  }
		  if (temp != colSum[j]) {
		    badSample=1;
		    break;
		  }

		}
		if (colSum[j]==0) {
		  for (m=0; m<NumRow; m++) {
		    currentS0RowSum[m]=currentS0RowSum[m]-s0[m][j];
		  }
		}

	      }

	      if (badSample==0) {
		for (m=0; m<NumRow; m++) {
		  if ((currentRowSum[m]>1) || (s0[m][j]==1 &&  currentRowSum[m]>0)) 
		    badSample=1;
		}
	      }

	      *pCurrentSample=currentSample;
	      *pBadSample=badSample;

	      free(currentRowSum);
	      free(currentS0RowSum);
              free(index);
	      free(rowSegment);
	      free(leftoverCol);
	      free(oneIndex);
}

static void UpdatePValue (int badSample, double currentSample, int **t, 
			  int NumRow, int NumCol, double sSquare0Bar, 
			  double *pSampleMean, 
		          double *pSampleVar, double *pPValue, int *pNumBad) 
{

	      if (badSample==0) {
		(*pSampleMean) = (*pSampleMean)+currentSample;
		(*pSampleVar)=(*pSampleVar)+currentSample*currentSample;

		if (SSquareBar(t, NumRow, NumCol)>=sSquare0Bar)	
		  (*pPValue)=(*pPValue)+currentSample;

	      }
	      else {
		(*pNumBad)++;
	      }
}	

static void PrintOutput (FILE *outFile, int N, int numBad, double sampleMean, double sampleVar, 
			 double pValue)
{	    
            double cvsquare1;

	    cvsquare1=(N-numBad)*sampleVar/sampleMean/sampleMean-1;	    
	    pValue=pValue/sampleMean;
	    sampleMean=sampleMean/N;

            fprintf(outFile, "%15g %15g %15g %10d\n", sampleMean, pValue, cvsquare1, numBad);
	    printf("%15g %15g %15g %10d\n", sampleMean, pValue, cvsquare1, numBad);
}

static double  NChooseM(int n, int m)
{
	int i;
	double value;

	value=1.0;

	if ((m<=0) || (n<=0) || (n<m)) return(value);
	for (i=1; i<=m; i++) {
		value=value*(n-i+1)/i;
	}
	return (value);

}

static double Mean(double x[], int n)
{
	int i;
	double sum=0;

	for(i=0; i<n; i++) sum = sum+x[i];
	return (sum/n);
}

static int RandomSample2(double weight[], int n)
{
	int low, up, mid;
	double randomNum;

	randomNum=RandomReal(0, weight[n-1]);
	low=0;
	up=n-1;
	mid=n/2;

	while (1) {
		if (randomNum>=weight[mid]) {
			low=mid+1;
			mid=(mid+up)/2;
			if (low==up) {
				mid=up;
				break;
			}
		}
		else if ((mid>0) && (randomNum<weight[mid-1])) {
			up=mid-1;
			mid=(mid+low)/2;
		}
		else break;
	}
	return (mid);
}

static double Drafting(int w[], int start, int end, int n, int a[],
		       int s[], int NumRow)
{
        double *p, *newP;
	int i,j;
	double tmp;
	double prob;
	double *currentWeight;
	double r1, r2;
	double *weight;

	p=malloc(NumRow*sizeof(double));
	newP=malloc(NumRow*sizeof(double));
	currentWeight=malloc(NumRow*sizeof(double));
	weight=malloc(NumRow*sizeof(double));

	if (n==0) return (1);

	if (end-start+1==n) {
	  for (i=start; i<=end; i++) {
	    a[i-start]=i;
	  }
	  return(1);
	}

	for (i=start; i<=end; i++) {
		if (w[i]==0) {
			weight[i]=0;
		}
		else {
			weight[i]=exp(log(1.0*w[i]/(s[i]-w[i])));
		}
	}

	for (i=start; i<=end; i++) {
		tmp=weight[i];
		weight[i]=weight[end];
		p[i]=tmp*R(n-1, weight, start, end-1, NumRow);
		weight[i]=tmp;
		currentWeight[i]=weight[i];
	}

	a[0]=ChooseOneInteger(p, start, end, NumRow);
	currentWeight[a[0]]=0;
	j=0;

	while (n>(j+1)) {
		for (i=start; i<=end; i++) {
		  if (currentWeight[i]<=0.00000000001)

				newP[i]=0;
		  else if (abs(weight[i]-weight[a[j]])<0.0000000001) {

				currentWeight[i]=0;
				r1=R(n-j-2, currentWeight, start, end, NumRow);
				currentWeight[i]=weight[i];
				r2=R(n-j-1, currentWeight, start, end, NumRow);
				newP[i]=weight[i]*r1/(n-j-1)/r2;
			}
			else {
				newP[i]=(weight[a[j]]*p[i]-weight[i]*p[a[j]])/(n-j-1)/(weight[a[j]]-weight[i])/p[a[j]];
			}
		}
		j++;

		a[j]=ChooseOneInteger(newP, start, end, NumRow);
		currentWeight[a[j]]=0;

		for (i=start; i<=end; i++) {
			p[i]=newP[i];
		}
	}

	prob=1;
	for (i=0; i<n; i++) {
		prob=prob*weight[a[i]];
	}
	prob=prob/R(n, weight, start, end, NumRow);

	free(p);
	free(newP);
	free(currentWeight);
	free(weight);

	return(prob);
}

static int ChooseOneInteger(double weight[], int start, int end, int NumRow)
{
	int i;
	double *totalWeight, temp;

	totalWeight=malloc(NumRow*sizeof(double));

	totalWeight[start]=weight[start];
	for (i=start+1; i<=end; i++) {
		totalWeight[i]=totalWeight[i-1]+weight[i];
	}

	temp=RandomReal(0,totalWeight[end]);

	i=start;
	while (totalWeight[i]<=temp) {
		i++;
	}
	free(totalWeight);
	return (i);
}

static double R(int k, double w[], int start, int end, int NumRow)
{
	int i, j;
	double *t; 
	double *r, temp;

	t=malloc((NumRow+1)*sizeof(double));
	r=malloc(NumRow*sizeof(double));

	for (i=1; i<=k; i++) {
		t[i]=T(i, w, start, end);
	}

	r[0]=1;
	for (j=1; j<=k; j++) {
		r[j]=0;
		for (i=1; i<=j; i++) {
			if (i%2==0) 
				r[j]=r[j]-t[i]*r[j-i];
			else 
				r[j]=r[j]+t[i]*r[j-i];
		}
		r[j]=r[j]/j;
	}
	temp=r[k];
	free(t);
	free(r);
	return (temp);
}

static double T(int i, double w[], int start, int end) 
{
	double sum;
	int j, k;
	double tmp;

	sum=0;
	for (j=start; j<=end; j++) {
	  tmp=1.0;
	  for (k=0; k<i; k++) {
	    tmp=tmp*w[j];
	  }
	  sum = sum+tmp;
	}
	return (sum);
}

static int RandomInteger(int low, int high)
{
  int k;
  double d;

  d=(double) rand()/((double) RAND_MAX+1);
  k=(int) (d*(high-low+1));
  return (low+k);
}

static double RandomReal(double low, double high)
{
       double d;

       d=(double) rand()/((double) RAND_MAX+1);
       return (low+d*(high-low));
}