#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#define LG 1000
#define NG 100


FILE *fp1, *fp2;

main()
{
  char hfile[40], s[]=" ", c[]="%";
/* INT */
  int ss, nburn, nf, jf, niter, nrun, err, nuse, ns, iburn, ig, is,
      jter, iter, nsave,
      m, mm, lag, am,
      ntaper[3] = {4, 8, 15};
/* DOUBLE */
  double cd[NG], cn[NG], cdn[NG], cdd[NG], cnn[NG], cvar[NG],
         rdn[NG], rdd[NG], rnn[NG], rnd[NG], g[LG],
         wtlog=0, pri, pdat, td, tn, tdn, tdd, tnn, tvar,
         gd, gn, gdn, gdd, gnn, gvar;

  double ad, an, eg, varg, sdg;
  double varnum, sdnum;
  double bard, barn, ann, and, adn, add, snn, snd, sdd, att;

  ss=scanf("%d", &nburn);
  if (ss!=1) { printf("read error in control file, line 1\n"); return -1;}
  if (nburn < 0) nburn = 0;
  ss=scanf("%d", &nsave);
  if (ss!=1) { printf("read error in control file, line 1\n"); return -1;}
  if (nsave) {
     ss=scanf("%s", &hfile);
     if (ss!=1) {printf("error reading file name\n"); return -1;}
     fp2=fopen(hfile,"w");
     if (fp2==NULL) {printf("failure to open file %s for output\n", hfile);
                     return -1;}
  }
  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;}
  nrun = 0;
  while (scanf("%d ", &jf) == 1) {
    err = 0;
                                         // start error checking
    if (jf <= 0) {
       printf("function of interest negative; skipped \n");
       err =  1;
    }
    rewind(fp1);                        //rewind the file hfile
    ss=fscanf(fp1,"%d %d", &niter, &nf);
    if (ss != 2) {printf("read error, first record of \ninput simulation ");
                  printf("file %s", hfile); return -1;} 
    if (jf > nf) {
       printf("function of interest %d greater than", jf);
       printf(" limit of %d ; skipped \n", nf);
       err =  1;
    }
    else if (nburn > niter) {
       printf(" %d burn in requested but only ", nburn); 
       printf(" %d avaliable \n", niter);
       return 1;
    }       
    if (nf > LG) {
       printf("file contains %d functions per record, ", nf);
       printf("exceeds storage capacity of %d \n", LG);
       return 1;
    }     
                                       //end of error checking
    if (!err) {
   
       nrun = nrun + 1;
       for (iburn=1; iburn <= nburn; iburn++) {
          if (readrec(nf,&jter,&wtlog,&pri,&pdat,g)) return 1;
       }

       nuse = niter - nburn;
       if (nuse < NG) {
          printf("minimum %d analysis iterations are required; ", NG);    
          printf(" %d avaliable \n", nuse);
          return 1;
       }
       ns = nuse / NG; nuse = ns * NG;
       if (nrun == 1) {
          header1(niter,nburn,nuse);
       }

/* form sufficiency statistics for all the subsequent work */
    
       td=0; tn=0; tdn=0; tdd=0; tnn=0; tvar=0;
       for (ig=0; ig <= NG-1; ig++) {                // possible array problem
          gd=0; gn=0; gdn=0; gdd=0; gnn=0; gvar=0;
          for (is=0; is <= ns-1; is++) {             // possible array problem
             if (readrec(nf,&jter,&wtlog,&pri,&pdat,g)) return 1;
             ad=0;
             if (wtlog > -700) ad = exp(wtlog);
             an=ad*g[jf-1];        
             gd=gd+ad; gn=gn+an; gdn=gdn+ad*an;
             gdd=gdd+ad*ad; gnn=gnn+an*an; gvar=gvar+an*g[jf-1]; 
          } /* for is */
          td=td+gd; tn=tn+gn; tdn=tdn+gdn;
          tdd=tdd+gdd; tnn=tnn+gnn; tvar=tvar+gvar;
   
          cn[ig]=gn/ns; cdn[ig]=gdn/ns; cdd[ig]=gdd/ns;
          cnn[ig]=gnn/ns; cvar[ig]=gvar/ns; cd[ig]=gd/ns;
       } /*for ig */

     /* form posterior moment and standard deviation, and print.*/
       eg=tn/td; varg=tvar/td - eg*eg;
       sdg= -1;
       if (varg > 0) sdg=sqrt(varg);
       if (nsave) 
          fprintf(fp2, "%4d%6d%16.7e%16.7e\n", jf, nuse, eg, sdg);
       printf("\n Function%4d", jf);
       printf(" Posterior mean:%16.7e", eg);
       printf(" Posterior st dev:%16.7e \n\n", sdg);
    /* numerical standard error assuming no serial correlation */
       varnum=(tnn-2*eg*tdn+tdd*eg*eg)/(td*td);
       sdnum= -1;
       if (varnum > 0) sdnum=sqrt(varnum);
       printf("%14s Accuracy of approximation:%9s Numerical",s,s);
       printf("%8s Relative \n %49s Standard%9s Numerical \n", s,s,s);
       printf(" Method-- %40s Error%12s Efficiency \n",s,s);
       printf("%8s Assuming no serial correlation:",s);
       printf(" %20.7e %15.6f \n", sdnum, (varg / (nuse*varnum)));
       if (nsave) fprintf(fp2,"%16.7e  ", sdnum);

    /* get autocovariance function of grouped means */

       barn=tn/nuse; bard=td/nuse;
       for (ig=1; ig <= NG; ig++) {                  
          cn[ig-1]=cn[ig-1]-barn;
          cd[ig-1]=cd[ig-1]-bard;
       }
       for (lag=0; lag <= NG-1; lag++) {
          ann=0; add=0; and=0; adn=0;
          for (ig=lag+1; ig <= NG; ig++) {
             ann=ann+cn[ig-1]*cn[ig-1-lag];
             add=add+cd[ig-1]*cd[ig-1-lag];
             and=and+cn[ig-1]*cd[ig-1-lag];
             adn=adn+cd[ig-1]*cn[ig-1-lag];
          }
          rnn[lag]=ann/NG; rdd[lag]=add/NG; rnd[lag]=and/NG; rdn[lag]=adn/NG;
        }    
     /* Numerical standard error with tapered autocovariance functions */
       for (mm=0; mm<=2; mm++) {
          m=ntaper[mm]; am=m;     
          snn=rnn[0];  sdd=rdd[0];  snd=rnd[0];
          for (lag=1; lag <= m-1; lag++) {
             att=1-(lag+0.0)/am;        // trick !!
             snn=snn+2*att*rnn[lag];
             sdd=sdd+2*att*rdd[lag];
             snd=snd+att*(rdn[lag]+rnd[lag]);
          }
          varnum=ns*nuse*(snn-2*eg*snd+sdd*eg*eg)/(td*td); 
          sdnum= -1;
          if (varnum > 0) sdnum = sqrt(varnum);
          printf("Autocovariance function tapered to %3d %s",m,c);
          printf(" %20.7e %15.6f\n", sdnum, varg/(nuse*varnum));  
          if (nsave) fprintf(fp2,"%16.7e  ", sdnum);
       } /* for */
       if (nsave) fprintf(fp2,"\n");
    }  /* if (!err) */
  }  /* while */
  
if (nrun == 0)
     printf("error on input for function number, program stops.\n"); 
     
fclose(fp1);
fclose(fp2);
return 1;  
} /* main*/