/*make sure you should let validLength be larger than minLenOfRepeat.*/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#include <math.h>
#include <ctype.h>

#define maxCloneLen 501
#define wordLen  18
#define genomelen 182625
#define maxReadNum 2558
#define groupNum  400
#define classNum  7

#define MBIG 1000000000
#define MSEED 16183398
#define MZ 0
#define FAC (1.0/MBIG)
static long inext,inextp;
static long ma[56];
static long iff=0;


typedef struct readlist
{
  long value;
  long positionInRead;
  struct readlist *next;
  } mylist;


typedef struct myhashElement
{
  int flag;
  double value;
  long  number;
  mylist  read;
  struct myhashElement *next;
  } mytable;





long factorial(long n)
{
   long i,j=1;

   for (i=2;i<=n;i++)
      j  *= i;
   return j;
}



void myEM(long *group,double *a,double *b,int *count)
{
   int i,j,k;
   double t1,t2,t3,t4,t5,c[classNum],d[classNum];

   if (*count>32000)
	   return;
   else
	   (*count)++;
   for (t1=0,i=1;i<groupNum;i++)
   		 t1 += group[i];
   for (t2=0,i=1;i<classNum;i++)
		 t2 += b[i]*exp(-a[i]);

   group[0] = (long)(t1/(1/t2 -1));

   
   for (i=0;i<classNum;i++)
   {
	   c[i] = b[i];
	   d[i] = a[i];
   }


   for (t3=0,i=0;i<classNum;i++)
   {
     b[i] = 0;
     for ( t4=t5=0,k=0;k<groupNum;k++)
	 {
		 for (t1=0,j=0;j<classNum;j++)
			 {
			     t2 = d[j]/d[i];
				 t1 += (c[j]*pow(t2,k)*exp(d[i]-d[j]));
			 }
   	     t5  = group[k]*c[i]/t1;
         b[i] += t5;
		 t4 += t5*k;
  	 }
	 a[i]   = t4/b[i];
	 t3  += b[i];
   }


   for (i=0;i<classNum;i++)
        b[i] /= t3 ;
 
   for (t1=0,i=0;i<classNum;i++)
	   t1 += ((d[i]-a[i])*(d[i]-a[i])+(b[i]-c[i])*(b[i]-c[i]));

   if ( sqrt(t1)>1e-3 )
	   myEM(group,a,b,count);
}



void myHist(mytable *data,long *group)
{
    mytable *p1;
	long i,j,m=genomelen-maxCloneLen+1;

   for ( i=0;i<m; i++)
   {
      p1 = &data[i];
	  while(p1!=NULL)
	  {
	    if ( p1->number!=0 )
		{
 		    if ( p1->number>groupNum-1 )
				j = groupNum-1;
			else
				j = p1->number;
			group[j]+=1;
		}
		p1 = p1->next;
	  }
  }
return;
}


double encode(char *str)
{
  double t1,t2;
  int i;

 for (t1=0,i=0;i<wordLen;i++) 
  {
    switch(str[i])
	{
	case 'A': t2 = 0;break;
	case 'C': t2 = 1;break;
	case 'G': t2 = 2;break;
	case 'T': t2=3;
	}
    t1 = t1*4 + t2;
  }

 return(t1);
}


void putIntoHashTable(double t1,mytable *head, long readnum,int i)
{
  int flag,l;
  long n,k=genomelen-maxCloneLen+1,L;
  mytable *p1,*p2,*p3;
  mylist *lp1,*lp2,*lp3;

  L = (long)(fmod(t1,k));
  p1 = &head[L];
  if ( p1->number==0 )
  {
	  p1->value=t1;
	  p1->read.value=readnum;
	  p1->read.positionInRead=i;
	  p1->read.next = NULL;
	  p1->number=1;
  }
  else
  {
	  flag =0;
	  while(p1!=NULL&&flag==0)
	  {
	    if ( p1->value!=t1 )
		{
			p3 = p1;
			p1= p1->next;
		}
		else
		{
			(p1->number)++;
			 lp1 = &(p1->read);
			 lp2 = lp1;
             n = lp1->value;
			 while(lp1!=NULL&&n<=readnum)
			 {
               lp2 = lp1;
               lp1 = lp1->next;
			   if (lp1!=NULL)
			      n = lp1->value;
			 }
			 
			 lp3 = (mylist *)(calloc)(1,sizeof(mylist));
			 lp3->value = readnum;
	         lp3->positionInRead=i;
			 lp3->next = NULL;
             if (lp1!=NULL)
			 {
                if (lp1==lp2)
				{
    				  n = lp1->value;
	                  l = lp1->positionInRead;
		    		  lp1->value = readnum;
	                  lp1->positionInRead= i;
					  lp3->next = lp1->next;
					  lp3->value = n;
		              lp3->positionInRead=l;
					  lp1->next = lp3;
				}
				else
				{
     				lp3->next = lp1;
	    			lp2->next = lp3;
				}
			}
			else
				 lp2->next = lp3;
 		flag =1 ;
		}
	  }

	  if ( flag == 0 )
	  {
	   p2 = (mytable *)(calloc)(1,sizeof(mytable));
	   if ( p2==NULL)
	   {
	    printf("\nWhat's wrong with the p2 in myhash");
		exit(0);
	   }
	   p2->flag = 1;
	   p3->next = p2;
	   p2->next= NULL;
	   p2->read.value=readnum;
	   p2->read.positionInRead=i;
	   p2->read.next = NULL;
	   p2->value = t1;
	   p2->number =1;
	  }
  } 

}

void myhash(char *str,mytable *head, long readnum)
{
  int i,length;
  double t1,t2,m=pow(4,wordLen);
  long  k=genomelen-maxCloneLen+1;

  length= strlen(str);
  if (length<wordLen)
	     return;
  t1 = encode(str);
  putIntoHashTable(t1,head,readnum,0);

  for ( i=wordLen;i<length;i++)
  {
    switch(str[i])
	{
	case 'A': t2 = 0;break;
	case 'C': t2 = 1;break;
	case 'G': t2 = 2;break;
	case 'T': t2=3;
	}
   t1 = fmod(t1*4,m)+t2;
   putIntoHashTable(t1,head,readnum,i-wordLen+1);   
 }


}


double random1(long *idum)
{
    long mj, mk;
	long i, ii,k;

	if(*idum < 0 || iff == 0)	{ /* initialization */
		iff = 1;
		mj = MSEED - (*idum < 0 ? -*idum : *idum);
		mj %= MBIG;
		ma[55] = mj;
		mk = 1;
		for(i = 1; i < 54; i ++)	{
			ii = (21 + i) % 55;
			ma[ii] = mk;
			mk = mj - mk;
			if(mk < MZ)	mk += MBIG;
			mj = ma[ii];
		}
		for(k = 1; k <= 4; k ++)	{
			for(i = 1; i <= 55; i ++)	{
				ma[i] -= ma[1 + (i + 30) % 55];
				if(ma[i] < MZ)	ma[i] += MBIG;
			}
		}
		inext = 0;
		inextp = 31;
		*idum = 1;
	}
	if( ++inext == 56) inext = 1;
	if( ++inextp == 56) inextp = 1;
	mj = ma[inext] - ma[inextp];
	if(mj < MZ) mj += MBIG;
	ma[inext] = mj;
	return(mj * FAC);
}


long  ran_number(long n, long *idum)
/*generate random number from 0 to n-1 */
{
	double	t;
	long    p;

	t = random1(idum);
	p = ( long ) (n * t);
	return( p );
}


void GetOccurrenceNum(char *str,mytable *head,long *occurrenceNum)
{
  mytable *p1;
  int i,n=0;
  double t1,t2,m=pow(4,wordLen);
  long  k=genomelen-maxCloneLen+1,L,length;

  length = strlen(str);
  t1 = encode(str);    
  L = (long)(fmod(t1,k));
  p1 = &head[L];
  while(p1->value!=t1)
	  p1 = p1->next;
  occurrenceNum[n]= p1->number;
  n++;
  for ( i=wordLen;i<length;i++)
  {
    switch(str[i])
	{
	case 'A': t2 = 0;break;
	case 'C': t2 = 1;break;
	case 'G': t2 = 2;break;
	case 'T': t2=3;
	}
    t1 = fmod(t1*4,m)+t2;
    L = (long)(fmod(t1,(double)(k)));
    p1 = &head[L];
    while(p1->value!=t1)
	  p1 = p1->next;
    occurrenceNum[n]= p1->number;
	n++;
  }
}



long partition(long *data1,long i, long j)
{
	long x,y;
	long m,n;
   
	x = data1[i];
    data1[i] =data1[j];
	data1[j] = x;
	n = j+1;
	m = i-1;
	while (m<n)
	{
	   do
	   {  n--;}
	   while(data1[n]>x);
	   do
	   {m++;}
	   while(data1[m]<x);
	   if (m<n)
	   {
	      y = data1[m];
		  data1[m]=data1[n];
		  data1[n]=y;
	   }
	}
  return n;
}



void quickSort(long *data1,long i,long j)
{
   long k;

   if (i<j)
   {
        k = partition(data1,i,j);
        quickSort(data1,i,k);
        quickSort(data1,k+1,j);
   }
}


void mySort(long *data1,int * data2,char *str1,long number)
{
   long i,j,k;
   int l;
   for (i=0;i<number-1;i++)
   for (j=i+1;j<number;j++)
	   if ( strlen(&str1[data1[i]*(maxCloneLen+1)])+data2[i]>strlen(&str1[data1[j]*(maxCloneLen+1)])+data2[j] )
	   {
	     k = data1[i];
		 data1[i] = data1[j];
		 data1[j] = k;
		 l = data2[i];
		 data2[i] = data2[j];
		 data2[j] = l;
	     }
}


void main()
{
	 mytable *head,*p1,*p2;
	 long occurrenceNum[maxCloneLen-wordLen+2],i,j,k1;
     double a[classNum],t1,t2,t3,t4,b[classNum],coverage;
	 char *str1,*str,*seq,ch;
     long m1,randSeed=15,m3,group[groupNum],readPosition[maxReadNum]; 
	 long length,readNum;
     FILE *fp1,*fp2;
     mylist *lp1,*lp2;
     int count;

	 head = (mytable *)(calloc)(genomelen-maxCloneLen+1,sizeof(mytable));
     str1 = (char *)(calloc)((maxCloneLen+1)*maxReadNum,sizeof(char));
     str = (char *)(calloc)(maxCloneLen+1,sizeof(char));
     seq = (char *)(calloc)(genomelen+1,sizeof(char));
	 if ( str1==NULL||head==NULL||str==NULL)
     {
       printf("\nwe can't find so much space for str");
       exit(0);
	 }

     occurrenceNum[maxCloneLen-wordLen+1] = 0;
	 m1 = genomelen - maxCloneLen +1;
     fp1 = fopen("test2.txt","aw");
	 fp2 = fopen("c:\\xiaomanFile\\NM\\repeatData\\H_NH0140H23\\p1seq.txt","r");
	 if ( fp1==NULL||fp2==NULL)
	 {
		 printf("\nThere is no enough space for output the files");
		 exit(0);
	 }

   i = 0;
   while (i<genomelen)
   {
       ch = fgetc(fp2);
	   if(ch=='\n')
	      ch = fgetc(fp2);
       seq[i]= ch;
       i++;
   }
   fclose(fp2);
   seq[genomelen]='\0';

   k1 = randSeed;
   fprintf(fp1,"\nresult by test2.c  randSeed=%ld\n ",k1);
   readNum = 1096;
   if (readNum>maxReadNum)
   {
	   printf("\nPlease increase maxReadNum");
	   exit(0);
   }

   for (i=0;i<readNum;i++)
	  readPosition[i] = ran_number(genomelen-500+1,&randSeed);

   for ( j=0;j<groupNum;j++)
       group[j] = 0;

  for ( m3=0;m3<genomelen-maxCloneLen+1;m3++ )
  {
      head[m3].number = 0 ;
      head[m3].flag = 1;
      head[m3].value = -100;
      head[m3].next = NULL ;
  }
	 
  for (t4=0,i=0;i<readNum;i++ )
  {
     for (j=0;j<500;j++)
		 str[j] = seq[readPosition[i]+j];
	 str[500] = '\0';
	 length = 500;
	 t4 += ((double)(length));
	 m3 = i*(maxCloneLen+1);
     strcpy(&str1[m3],str);
	 str1[m3+length]='\0';
	 myhash(str,head,i);
  }
  free(seq);
  coverage = t4/genomelen; 
  myHist(head,group);
  for (j=0;j<classNum;j++)
  {
    a[j] = 5*j+8;
    b[j] = 1.0/classNum;
   }
  
  count = 0;
  myEM(group,a,b,&count);
/*
   for (k3=0;k3<groupNum;k3++)
	 {
      if (k3%20==0)
	 	 fprintf(fp1,"\n");
       fprintf(fp1,"%ld ",group[k3]);
	 } 
*/
  for (t1=0,j=0;j<classNum;j++)
    	t1 += (b[j]*a[j]);
  for ( j=0;j<classNum;j++)
       fprintf(fp1,"\nclass%f, percent=%f",a[j],b[j]*a[j]/t1);
  for (t3=0,j=0;j<classNum;j++)
       t3 += (b[j]*a[j]/a[0]);
  for (t2=0,j=1;j<groupNum;j++)
		t2 += group[j];
 fprintf(fp1,"\ncoverage = %f, genome1= %f,genome2 = %f",coverage,t3*t2,t4/a[0]);
 fprintf(fp1,"\nthe relative error of estimation of coverage=%f",a[0]/coverage);

    for (m3=0;m3<genomelen-maxCloneLen+1;m3++ )
	{
         p1 = head[m3].next;
		 head[m3].next =NULL;
	     while (p1!=NULL)
		 {
	         p2 = p1->next;
		     lp1 = &p1->read;
		     lp1 = lp1->next;
		     p1->read.next=NULL;
		     while (lp1!=NULL)
			 {
		       lp2 = lp1->next;
			   free(lp1);
			   lp1 = lp2;
			 }
	         free(p1);
	         p1 = p2;
		 }
	}


  for ( m3=0;m3<genomelen-maxCloneLen+1;m3++ )
   {
      head[m3].number = 0 ;
      head[m3].flag = 1;
      head[m3].value = -100;
      head[m3].next = NULL ;
  }


free(str1);
free(str);
free(head);
fclose(fp1);
}