#include <math.h>
#define LDPAR 250


FILE *fp1;

main()
{

/* INT */
 
  int ss, niter, npar, kpars, stat, i, ng, iter, jter, nover[5], j,
      ig, iter1, iter2;
  int ncrit=5, nline=40, which=2;
  char hfile[40], s[]=" "; 
/*DOUBLE */
  double dkpars, crit[5], bound, p, q, dg,
         g[LDPAR], pri, wtlog, pdat, gmax, amax;
double tt;
  double pcrit[5]={0.001, 0.01, 0.05, 0.1, 0.5};
/* 
  double thetai[LDPAR],
  int ss, niter, npar, jter, iter, nc;
  char hfile[40], simfile[40], 
*/

  ss=scanf("%s", &hfile);
  if (ss!=1) {printf("error reading file name\n"); return -1;}
  fp1=fopen(hfile,"r");
  if (fp1==NULL) {printf("failure to open file %s for input\n", hfile); 
                  return -1;}

  ss=fscanf(fp1,"%d %d", &niter, &npar);
  if (ss != 2) {printf("read error, first record of \ninput simulation ");
                printf("file %s", hfile); return -1;}
  if ((niter <= 0) || (npar<=0)) {
     printf("illegal first record of  input simulation");  
     printf("file %s %16d%16d", hfile, niter, npar); return -1;
  }
  if (npar > LDPAR) {
     printf("input simulation file %s \n ", hfile);  
     printf("%16d parameters exceed limit of %16d \n", npar, LDPAR); return -1;
  }
  scanf("%d",&kpars);
  dkpars=kpars; 
  for (i=1; i <= ncrit; i++) {
     q=pcrit[i-1];
     p=1-q;
     cdfchi(&which,&p,&q,&crit[i-1],&dkpars,&stat,&bound);
     crit[i-1]=crit[i-1]/2;
  }
  ng=niter/nline;
  dg=ng;

/* first pass */
  for (i=1; i<= nline*ng; i++) {
    if (readrec(npar,&iter,&wtlog,&pri,&pdat,&g)) return 1; 
    if (i==1) {
       jter=iter; gmax=pdat;
    }
    if (pdat>gmax){
       jter=iter; gmax=pdat;
    }
  }
  printf(" Largest log data pdf: %14.6e \n", gmax);
  printf("         at iteration: %14d \n", jter);

/* second pass to get summary statistics by group */
  printf("%35s Fraction downhill from chisquare(%3d) at p=\n",s,kpars);
  printf(" Group    Iterations   Max log data pdf");
  for (i=1; i<=ncrit; i++)
     printf("   %5.3f",pcrit[i-1]);
  printf("\n");

  rewind(fp1);
  fscanf(fp1,"%d %d", &niter, &npar);
  for (i=1; i<=nline; i++) {
     for (j=1; j<=ncrit; j++)
        nover[j-1]=0;
     for (ig=1; ig<=ng; ig++) {
        if (readrec(npar,&iter,&wtlog,&pri,&pdat,&g)) return 1; 
        if (ig==1) iter1=iter;
        if (ig==ng) iter2=iter;
        if (ig==1) amax=pdat;
        if (pdat>amax) amax=pdat;
        for (j=1; j<=ncrit; j++)
           if ((gmax-pdat)>crit[j-1]) nover[j-1]=nover[j-1]+1;
     }
     printf("%6d%7d-%7d%18.6e",i,iter1,iter2,amax);
     for (j=1; j<=ncrit; j++)
        printf("   %5.3f",nover[j-1]/dg);
     printf("\n");
  }
fclose(fp1);
return 1;
} //main