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

#define MIN(x, y) ( (x)<(y) ? (x) : (y) )
#define MAX(x, y) ( (x)>(y) ? (x) : (y) )   

struct hapdata {
	int n;
	int n00;
	int n01;
	int n10;
	int n11;
	int n0q;
	int n1q;
	int nq0;
	int nq1;
	int nqq;
	double  d12;
	int ad ;
	} *hv ;	

double maxfactor, minfactor;
#define TOL  0.005

int npairs;
int *n  ;

char **fname ;
int numfiles ;
double cmax,cmin, cnum  ;
double conv, conlen;
double *****prob, *recrates , **ppoly ;
int nsitesp;

main(int argc, char *argv[] )
{
	int i,j,k, nmax,  n1, n2,  nsites,  howmany, sum, rsam  ;
	int n00, n01, n10, n11, n0q, n1q, nq0, nq1, nqq ;
	double pmle, findmax(int, double, double * ), c, cstart,
	   lnlikemshap( double,double, double, int,double *,double *****, double **), d12 , r, theta ;
	double minusone = -1.0 , alpha , lnlike , fac ;
	FILE *pout, *fopen(const char *, const char *);
	unsigned maxpair = 1000 ;
	double  ****getprobmat(char *, int *, int *, double **, double ** )  ;
	double cbest, convbest, bestlike, cl, fmin, fmax , rhosim, fsim ;
	double conlensim, conlenin ;
	int fnum, cnum, rhogrid, rhogeom, fgeom, argn, simflag, printflag ;
	double ffac, rhofac ;
	char adstr[5] ;
 
  if( argc < 7 ){
	printf("usage:  maxhap numfiles hxxrho ...  rhostart maxfactor minfactor fmin fmax fnum [-g] conlen [-p] [rho f conlen] <datafile \n"); exit(1);
	}
	
  simflag =  printflag = 0 ;
  rhogrid = rhogeom = fgeom = 0 ;
  cnum = 0 ;
  argn = 1 ;
  if( (argv[argn][0] < '1') || (argv[argn][0] > '9' ) ){ fprintf(stderr," First argument must be number of files.\n"); exit(1);}
  numfiles = atoi( argv[argn++]) ;
  fname = (char **)malloc( (size_t)numfiles*sizeof( char *) ) ;
  for( i=0; i< numfiles; i++) fname[i] = argv[argn++] ;
  if( argv[argn][1] == 'g' ) {    /*  grid of  rho values */
     if( argv[argn+1][0] == '-' ) { fprintf(stderr," command line error.\n"); exit(1); }
     cmin = atof( argv[++argn] ) ;
     if( argv[argn+1][0] == '-' ) { fprintf(stderr," command line error.\n"); exit(1); }
     cmax = atof( argv[++argn] ) ;
     if( argv[argn+1][0] == '-' ) { fprintf(stderr," command line error.\n"); exit(1); }
     cnum = atof( argv[++argn] ) ;
     rhogrid = 1 ;
     if( argv[++argn][1] == 'g'){argn++; rhogeom = 1 ; rhofac = pow(cmax/cmin, 1.0/(double)(cnum-1));}
     if( rhogeom && (cmin == 0.0 ) ){ fprintf(stderr,"with -g, rhomin must be > 0.0 .\n"); exit(1); }
  }
  else {
     if( argv[argn][0] == '-' ) { fprintf(stderr," command line error.\n"); exit(1); }
     cstart = atof( argv[argn++] ) ;
     if( cstart == 0.0 ) { fprintf(stderr," cstart must be > 0.\n"); exit(1);}
     if( argv[argn][0] == '-' ) { fprintf(stderr," command line error.\n"); exit(1); }
     maxfactor = atof( argv[argn++] ) ;
     if( argv[argn][0] == '-' ) { fprintf(stderr," command line error.\n"); exit(1); }
     minfactor = atof( argv[argn++] ) ;
  }
  if( argc <= argn) { fprintf(stderr,"not enough command line arguments.\n"); exit(1); }
  if( argv[argn][0] == '-' ) { fprintf(stderr," command line error.\n"); exit(1); }
  fmin = atof( argv[argn++] ) ;
  if( argc <= argn) { fprintf(stderr,"not enough command line arguments.\n"); exit(1); }
  if( argv[argn][0] == '-' ) { fprintf(stderr," command line error.\n"); exit(1); }
  fmax = atof( argv[argn++] ) ;
  if( fmax < fmin ) { fprintf(stderr,"fmax < fmin\n"); exit(0);}
  if( argc <= argn) { fprintf(stderr,"not enough command line arguments.\n"); exit(1); }
  if( argv[argn][0] == '-' ) { fprintf(stderr," command line error.\n"); exit(1); }
  fnum = atoi( argv[argn++] ) ;
  if( argc <= argn) { fprintf(stderr,"not enough command line arguments.\n"); exit(1); }
  if( argv[argn][1] == 'g'){argn++; fgeom = 1 ; ffac = pow(fmax/fmin, 1.0/(double)(fnum-1));}
  if( fgeom && (fmin == 0.0 ) ){ fprintf(stderr,"with -g, fmin must be > 0.0 .\n"); exit(1); }
  if( argc <= argn) { fprintf(stderr,"not enough command line arguments.\n"); exit(1); }
  if( argv[argn][0] == '-' ) { fprintf(stderr," command line error.\n"); exit(1); }
  conlen = atof( argv[argn++] ) ;
  conlenin = conlen ;
  if( argc > argn ){
     if( argv[argn][1] == 'p' ){
        argn++;
        printflag = 1 ;
     }
  }
  if( argc > argn ){
      simflag = 1;
      rhosim = atof( argv[argn++]);
      if( rhosim < 0.0 ){ fprintf(stderr,"rho must be >= 0.0 .\n"); exit(1); }
      if( argc <= argn) { fprintf(stderr,"not enough command line arguments.\n"); exit(1); }
      fsim = atof( argv[argn++]);
      if( fsim < 0.0 ){ fprintf(stderr,"f must be >= 0.0 .\n"); exit(1); }
      if( argc <= argn) { fprintf(stderr,"not enough command line arguments.\n"); exit(1); }
      conlensim = atof( argv[argn++] );
      if( conlensim < 0.0 ){ fprintf(stderr,"conlen must be >= 0.0 .\n"); exit(1); }
  }
  for( i=0; i<argc; i++) fprintf(stdout,"%s ",argv[i] ) ;
  fprintf(stdout,"\n");
/*
	fprintf(stdout,"# rho_guess: %lf  fmin: %lf  fmax: %lf fnum: %d conversion tract length(mean): %lf\n",
	  cstart, fmin, fmax, fnum, conlen);
	if( argc >(numfiles+8) ) fprintf(stdout,"# cmin: %lf cmax: %lf cnum: %d\n", cmin, cmax, cnum);
	  
*/
  prob = (double *****)malloc( (size_t)numfiles*sizeof( double ****) ) ;
  ppoly = (double **)malloc( (size_t)numfiles*sizeof( double * ) ) ;
  n = (int *)malloc( (size_t)numfiles*sizeof( int) ) ;
  for( i=0; i<numfiles; i++){
	prob[i] = getprobmat(fname[i], &(n[i]),  &nsitesp, &recrates , &(ppoly[i]) ) ; 
	if( (i>0) && ( n[i] > n[i-1] ) ) { fprintf(stderr,"likefiles in wrong order.\n"); exit(0);}
  }
  hv  = (struct hapdata *) malloc( maxpair*sizeof( struct hapdata ) ) ; 
  nmax = 0 ;
  for( i=0; i<numfiles; i++) if( n[i] >nmax ) nmax = n[i] ;

   while( 1) {
	npairs = 0 ;
	while( 	scanf(" %*d %*d %lf %s %d %d %d %d %d %d %d %d %d", &d12 , adstr, 
	  &n00, &n01,  &n10, &n11, &n0q, &n1q, &nq0, &nq1, &nqq ) == 11 ) {
	   if( (npairs+1)  >= maxpair ){
		maxpair = npairs + 10 ;
		hv = (struct hapdata *) realloc( hv , maxpair*sizeof( struct hapdata) ) ;
	   }
	   sum = n00+n01+n10+n11+n0q+n1q+nq0+nq1+nqq ; 
	   if( (sum - nqq ) > nmax ) { printf(" sample too big for prob file.\n"); exit(1); }
	   else {
		rsam = sum - nqq ; 
		i = numfiles -1 ;
		while( rsam > n[i] )  i-- ;
		nqq  = n[i] - rsam ;
		hv[npairs].n = n[i] ;
	   }
	    if( adstr[0] == 'u' ) hv[npairs].ad = 0 ;
	    else hv[npairs].ad = 1 ;
	    hv[npairs].n00 = n00 ;
	    hv[npairs].n01 = n01 ;
	    hv[npairs].n10 = n10 ;
	    hv[npairs].n11 = n11 ;
	    hv[npairs].n0q = n0q ;
	    hv[npairs].n1q = n1q ;
	    hv[npairs].nq0 = nq0 ;
	    hv[npairs].nq1 = nq1 ;
	    hv[npairs].nqq = nqq ;
	    hv[npairs].d12 = fabs(d12) ; 
	    npairs++;
         }
			
         if( npairs == 0 ) break;
         printf("# npairs: %d\n", npairs);
	 scanf("%*s");

	if( rhogrid ){
           c = cmin;
	   conv = fmin;
           bestlike = lnlikemshap( c,conv,conlen,nsitesp,recrates,prob,ppoly  )  ;
	   cbest = cmin;
	   convbest = fmin;
       	   for( j=0; j< cnum ; j++ ){
	      conv = fmin ;
              for( i=0; i < fnum; i++) {
                 lnlike = lnlikemshap( c,conv,conlen,nsitesp,recrates,prob,ppoly  )  ;
		 if( printflag) fprintf(stdout,"%lf %lf %lf\n", c, conv, lnlike);
		 if( lnlike > bestlike) { cbest = c; convbest = conv; bestlike = lnlike;}
                 if( (fnum > 1) && fgeom ) conv *= ffac ;
	         else if( fnum > 1) conv  = fmin + (i+1)*(fmax-fmin)/(double)(fnum-1) ;
	         else conv = fmin ;
              }    
              if( (cnum > 1) && rhogeom ) c *= rhofac;
	      else if( cnum > 1)  c = cmin + (j+1)*(cmax-cmin)/(double)(cnum-1) ;
	      else c = cmin ;
	   }
	   fprintf(stdout,"Mle: %lf %lf %lf\n", cbest, convbest, bestlike);
        }
        else{
	   if( cstart > 0.0) {
        	conv = fmin;
		pmle = findmax( npairs, cstart, &bestlike ) ;
		if( printflag ) fprintf(stdout,"mle-rho(given f= %lf): %lf (%lf)\n", conv, pmle, -1.0*bestlike);
        	cbest = pmle;
        	convbest = conv;
          	for( i=1; i < fnum; i++) {
           	   if( !fgeom) conv  = fmin + i*(fmax-fmin)/(double)(fnum-1) ;
		   else conv *= ffac ;
		   pmle = findmax( npairs, cstart, &lnlike ) ;
		   if( printflag) fprintf(stdout,"mle-rho(given f= %lf): %lf (%lf)\n", conv, pmle, -1.0*lnlike);
                   if(  lnlike < bestlike) {
                      convbest = conv; cbest = pmle; bestlike = lnlike; 
		   }
                } 
		fprintf(stdout,"Mle: %lf %lf %lf \n", cbest, convbest, -1.0*bestlike ) ;
	    }
	 }
         if( simflag ){
	    fprintf(stdout," %lf %lf %lf %lf\n", rhosim, fsim,conlensim, 
		lnlikemshap( rhosim,fsim, conlensim, nsitesp,recrates,prob,ppoly)  ) ;
 	 }
	 fprintf(stdout,"\n");
  }
}	
	

	double
findmax(int npairs, double cstart, double *plnlike )
{
	double L0, L05, L100 , lnlikemshap( double,double, double, int, double *, double *****, double ** ) ;
	float golden(float, float, float, float (*f)(float),float, float *),cx,tol,cle  ;
	float ax, bx,  mlike( float) ;
	double cn, cnn, lcn, lcnn, c0 ;

	/* check for zero max, and 100 max , if not use 0.0 and 100.0 as brackets */
	L0 = lnlikemshap( 0.0,conv,conlen, nsitesp,recrates, prob, ppoly ) ;
	cn = cstart ;
	lcn  = lnlikemshap( cn ,conv,conlen, nsitesp, recrates,prob, ppoly ) ;
	c0 = 0.0 ;

if( lcn > L0 ) {
   while(  cn <= maxfactor*cstart ) {
	cnn = cn*1.6 ;
	lcnn  = lnlikemshap( cnn ,conv, conlen, nsitesp, recrates,prob, ppoly ) ;
	if( lcnn < lcn ) {
	   ax = c0; bx = cn; cx = cnn;
	   break;
	   }
	if( lcnn > lcn ) c0 = cn;
	cn = cnn;
	lcn = lcnn;
	}
   if( cn > maxfactor*cstart ) { fprintf(stdout, " Warning: maximum value of rho.\n");  *plnlike = -1.0*lcn; return( cn);}
}

else{
  while( cn >= minfactor*cstart ) {
	cnn = 0.25*cn ;
	lcnn  = lnlikemshap( cnn , conv, conlen, nsitesp, recrates,prob, ppoly ) ;
	if( lcnn > L0 ) {
	    ax = 0.0;
	    bx = cnn;
	    cx = cn;
	    break;
	    }
	cn = cnn;
	lcn = lcnn;
	}
  if( cn <= minfactor*cstart ) { *plnlike = -1.0*L0 ; return( 0.0 ) ;}
}
	tol = TOL ;
	*plnlike = golden(ax, bx, cx, mlike, tol, &cle) ; 
	return( (double)cle ) ;
} 


	float
mlike( float fc )
{
	double lnlikemshap( double, double, double,int, double *, double *****, double ** ) ;

  return( (float)( -1.0*lnlikemshap( (double) fc,conv,conlen, nsitesp, recrates, prob, ppoly ) ) ) ;
}


#define R 0.61803399
#define C (1.0-R)
#define SHFT2(a,b,c) (a)=(b);(b)=(c);
#define SHFT3(a,b,c,d) (a)=(b);(b)=(c);(c)=(d);

float golden(float ax,float bx,float cx,float (*f)(float) ,float tol,float *xmin)
{
        float f1,f2,x0,x1,x2,x3;

        x0=ax;
        x3=cx;
        if (fabs(cx-bx) > fabs(bx-ax)) {
                x1=bx;
                x2=bx+C*(cx-bx);
        } else {
                x2=bx;
                x1=bx-C*(bx-ax);
        }
        f1=(*f)(x1);
        f2=(*f)(x2);
        while ( ( fabs(x3-x0) > tol*(fabs(x1)+fabs(x2) ))&&((x0+x1+x2+x3)>0.00001)) {

                if (f2 < f1) {
                        SHFT3(x0,x1,x2,R*x1+C*x3)
                        SHFT2(f1,f2,(*f)(x2))
                } else {
                        SHFT3(x3,x2,x1,R*x2+C*x0)
                        SHFT2(f2,f1,(*f)(x1))
                }
        }
        if (f1 < f2) {
                *xmin=x1;
                return f1;
        } else {
                *xmin=x2;
                return f2;
        }
}
#undef C
#undef R
#undef SHFT2
#undef SHFT3

