/*
Created by YongYong Shi, 2003.11.28
*/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

#ifndef RAND_MAX
#include <limits.h>
#define RAND_MAX INT_MAX
#endif

#define error(s,t) error_func(__LINE__,__FILE__,s,t)

int _CONTROL[2][3] = {{172, 6, 13},
		      {67, 115, 10}};
int _CASE[2][3] = {{173, 6, 13},
		   {71, 110, 11}};
		   
error_func(unsigned l,char *f,char *s1,char *s2)
{
  fprintf(stderr,"Error on line %d of source file %s:\n %s%s\n\n",l,f,s1,s2);
  exit(1);
  return 0;
}

main(int argc,char *argv[])
{
  int i, k, seed;
  float fs0, fs1, fn0, fn1, x;
  unsigned N, n;
  unsigned AC=0, CA=0, AD=0, DA=0, BC=0, CB=0, BD=0, DB=0;
  float O11, O12, O21, O22, C1, C2, R1, R2, T;
  float c2A, c2B, c2C, c2D;
  FILE *fi, *fo;
  
  if (argc<3)
  {
    printf("\n\nUsage: sim <input file> <output file>\n");
    return 1;
  }
  else
  {
    fi=fopen(argv[1],"r");
    if (fi==NULL) { error("Could not open file: ",argv[1]); return 1; }
    printf("Reading %s\n",argv[1]);
    fscanf(fi,"%f\n", &fs0);
    fscanf(fi,"%f\n", &fs1);
    fscanf(fi,"%f\n", &fn0);
    fscanf(fi,"%f\n", &fn1);
    fscanf(fi,"%u\n", &N);
    fscanf(fi,"%d\n", &seed);
    fclose(fi);
    fo=fopen(argv[2],"w");
    if (fo==NULL) { error("Could not open file: ",argv[2]); return 1; }
  }
  srand(seed);
  
  //start simulating
  for(n=0;n<N;n++)
  {
    k=n/50000*50000-n;
    if(!k)printf("Step %u\n", n);
    
    int South[2][2]={{0,0},{0,0}};
    int North[2][2]={{0,0},{0,0}};
    int NSHC[2][2]={{0,0},{0,0}};
    int NNHC[2][2]={{0,0},{0,0}};
    //---------------------------------------------------------------------------
    for(i=0;i<_CONTROL[0][0];i++)  //allocate allele to Normals of (South|NSHC) 
    {
      x=(float)rand()/(float)RAND_MAX;
      if(x<fs0)
      {
      	South[0][0]++;
      	NSHC[0][0]++;
      }
      else
      {
      	South[0][1]++;
      	NSHC[0][1]++;
      }
      x=(float)rand()/(float)RAND_MAX;
      if(x<fs0)
      {
      	South[0][0]++;
      	NSHC[0][0]++;
      }
      else
      {
      	South[0][1]++;
      	NSHC[0][1]++;
      }
    }
    for(i=0;i<_CONTROL[0][1];i++)  //allocate allele to Normals of (South|NNHC)
    {
      x=(float)rand()/(float)RAND_MAX;
      if(x<fs0)
      {
      	South[0][0]++;
      	NNHC[0][0]++;
      }
      else
      {
      	South[0][1]++;
      	NNHC[0][1]++;
      }
      x=(float)rand()/(float)RAND_MAX;
      if(x<fs0)
      {
      	South[0][0]++;
      	NNHC[0][0]++;
      }
      else
      {
      	South[0][1]++;
      	NNHC[0][1]++;
      }
    }
    for(i=0;i<_CONTROL[0][2];i++)  //allocate allele to Normals of (South|Undet)
    {
      x=(float)rand()/(float)RAND_MAX;
      if(x<fs0)
      {
      	South[0][0]++;
      }
      else
      {
      	South[0][1]++;
      }
      x=(float)rand()/(float)RAND_MAX;
      if(x<fs0)
      {
      	South[0][0]++;
      }
      else
      {
      	South[0][1]++;
      }
    }
    for(i=0;i<_CONTROL[1][0];i++)  //allocate allele to Normals of (North|NSHC)
    {
      x=(float)rand()/(float)RAND_MAX;
      if(x<fn0)
      {
      	North[0][0]++;
      	NSHC[0][0]++;
      }
      else
      {
      	North[0][1]++;
      	NSHC[0][1]++;
      }
      x=(float)rand()/(float)RAND_MAX;
      if(x<fn0)
      {
      	North[0][0]++;
      	NSHC[0][0]++;
      }
      else
      {
      	North[0][1]++;
      	NSHC[0][1]++;
      }
    }
    for(i=0;i<_CONTROL[1][1];i++)  //allocate allele to Normals of (North|NNHC)
    {
      x=(float)rand()/(float)RAND_MAX;
      if(x<fn0)
      {
      	North[0][0]++;
      	NNHC[0][0]++;
      }
      else
      {
      	North[0][1]++;
      	NNHC[0][1]++;
      }
      x=(float)rand()/(float)RAND_MAX;
      if(x<fn0)
      {
      	North[0][0]++;
      	NNHC[0][0]++;
      }
      else
      {
      	North[0][1]++;
      	NNHC[0][1]++;
      }
    }
    for(i=0;i<_CONTROL[1][2];i++)  //allocate allele to Normals of (North|Undet)
    {
      x=(float)rand()/(float)RAND_MAX;
      if(x<fn0)
      {
      	North[0][0]++;
      }
      else
      {
      	North[0][1]++;
      }
      x=(float)rand()/(float)RAND_MAX;
      if(x<fn0)
      {
      	North[0][0]++;
      }
      else
      {
      	North[0][1]++;
      }
    }
    //----------------------------------------------------------------------------
    for(i=0;i<_CASE[0][0];i++)  //allocate allele to Cases of (South|NSHC) 
    {
      x=(float)rand()/(float)RAND_MAX;
      if(x<fs1)
      {
      	South[1][0]++;
      	NSHC[1][0]++;
      }
      else
      {
      	South[1][1]++;
      	NSHC[1][1]++;
      }
      x=(float)rand()/(float)RAND_MAX;
      if(x<fs1)
      {
      	South[1][0]++;
      	NSHC[1][0]++;
      }
      else
      {
      	South[1][1]++;
      	NSHC[1][1]++;
      }
    }
    for(i=0;i<_CASE[0][1];i++)  //allocate allele to Cases of (South|NNHC)
    {
      x=(float)rand()/(float)RAND_MAX;
      if(x<fs1)
      {
      	South[1][0]++;
      	NNHC[1][0]++;
      }
      else
      {
      	South[1][1]++;
      	NNHC[1][1]++;
      }
      x=(float)rand()/(float)RAND_MAX;
      if(x<fs1)
      {
      	South[1][0]++;
      	NNHC[1][0]++;
      }
      else
      {
      	South[1][1]++;
      	NNHC[1][1]++;
      }
    }
    for(i=0;i<_CASE[0][2];i++)  //allocate allele to Cases of (South|Undet)
    {
      x=(float)rand()/(float)RAND_MAX;
      if(x<fs1)
      {
      	South[1][0]++;
      }
      else
      {
      	South[1][1]++;
      }
      x=(float)rand()/(float)RAND_MAX;
      if(x<fs1)
      {
      	South[1][0]++;
      }
      else
      {
      	South[1][1]++;
      }
    }
    for(i=0;i<_CASE[1][0];i++)  //allocate allele to Cases of (North|NSHC)
    {
      x=(float)rand()/(float)RAND_MAX;
      if(x<fn1)
      {
      	North[1][0]++;
      	NSHC[1][0]++;
      }
      else
      {
      	North[1][1]++;
      	NSHC[1][1]++;
      }
      x=(float)rand()/(float)RAND_MAX;
      if(x<fn1)
      {
      	North[1][0]++;
      	NSHC[1][0]++;
      }
      else
      {
      	North[1][1]++;
      	NSHC[1][1]++;
      }
    }
    for(i=0;i<_CASE[1][1];i++)  //allocate allele to Cases of (North|NNHC)
    {
      x=(float)rand()/(float)RAND_MAX;
      if(x<fn1)
      {
      	North[1][0]++;
      	NNHC[1][0]++;
      }
      else
      {
      	North[1][1]++;
      	NNHC[1][1]++;
      }
      x=(float)rand()/(float)RAND_MAX;
      if(x<fn1)
      {
      	North[1][0]++;
      	NNHC[1][0]++;
      }
      else
      {
      	North[1][1]++;
      	NNHC[1][1]++;
      }
    }
    for(i=0;i<_CASE[1][2];i++)  //allocate allele to Cases of (North|Undet)
    {
      x=(float)rand()/(float)RAND_MAX;
      if(x<fn1)
      {
      	North[1][0]++;
      }
      else
      {
      	North[1][1]++;
      }
      x=(float)rand()/(float)RAND_MAX;
      if(x<fn1)
      {
      	North[1][0]++;
      }
      else
      {
      	North[1][1]++;
      }
    }
    //-------------------------------------------End of allocation
    
    O11=South[0][0];
    O12=South[0][1];
    O21=South[1][0];
    O22=South[1][1];
    C1=O11+O21;
    C2=O12+O22;
    R1=O11+O12;
    R2=O21+O22;
    T=C1+C2;
    c2A=(abs(O11*O22-O12*O21)-T/2)*(abs(O11*O22-O12*O21)-T/2)*T/R1/R2/C1/C2;
    
    O11=North[0][0];
    O12=North[0][1];
    O21=North[1][0];
    O22=North[1][1];
    C1=O11+O21;
    C2=O12+O22;
    R1=O11+O12;
    R2=O21+O22;
    T=C1+C2;
    c2B=(abs(O11*O22-O12*O21)-T/2)*(abs(O11*O22-O12*O21)-T/2)*T/R1/R2/C1/C2;
    
    O11=NSHC[0][0];
    O12=NSHC[0][1];
    O21=NSHC[1][0];
    O22=NSHC[1][1];
    C1=O11+O21;
    C2=O12+O22;
    R1=O11+O12;
    R2=O21+O22;
    T=C1+C2;
    c2C=(abs(O11*O22-O12*O21)-T/2)*(abs(O11*O22-O12*O21)-T/2)*T/R1/R2/C1/C2;
    
    O11=NNHC[0][0];
    O12=NNHC[0][1];
    O21=NNHC[1][0];
    O22=NNHC[1][1];
    C1=O11+O21;
    C2=O12+O22;
    R1=O11+O12;
    R2=O21+O22;
    T=C1+C2;
    c2D=(abs(O11*O22-O12*O21)-T/2)*(abs(O11*O22-O12*O21)-T/2)*T/R1/R2/C1/C2;
    
    if(c2A>c2C)AC++;
    else CA++;
    if(c2A>c2D)AD++;
    else DA++;
    if(c2B>c2C)BC++;
    else CB++;
    if(c2B>c2D)BD++;
    else DB++;
  }
  
  printf("A>C: %u ; C>A: %u\n", AC, CA);
  printf("A>D: %u ; D>A: %u\n", AD, DA);
  printf("B>C: %u ; C>B: %u\n", BC, CB);
  printf("B>D: %u ; D>B: %u\n", BD, DB);
  
  fprintf(fo, "A>C: %u ; C>A: %u\n", AC, CA);
  fprintf(fo, "A>D: %u ; D>A: %u\n", AD, DA);
  fprintf(fo, "B>C: %u ; C>B: %u\n", BC, CB);
  fprintf(fo, "B>D: %u ; D>B: %u\n", BD, DB);
  
  fclose(fo);
  
  return 0;
}