#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#define LD 50000
#define LDPAR 300
#define CON3 0.9375

FILE *fp1, *fp2;

extern double ab(double a);
extern double sq(double a);
main()
{
  char *hkern[3]={"  Uniform", "Triangular", "Biweight"};
  char *hrange[2]={"         ", "Quantile"};
  int kern, krange, nburn, npts, ifunc, ss, niter, npars, nuse, nskip,
      iburn, jter, iter, m, mtot, nv1[LD], mbase, ipts;
  double a1, a2, hfrac, wtlog, pri, pdat, par[LDPAR], wtmax, g[LD], wt[LD],
      b1, b2, r1, r2, wsum, wrun, wrun1, q1, q2, w[LD], h, h1, z, zinc, ak,
      zlower, zupper, den;
  char hfile[40], s[]=" ";


/* read in and check user info */
  printf(" Preparation of posterior density for plotting\n");
  ss=scanf("%d%lf%d%lf%lf%d%d%d", &kern,&hfrac,&krange,&a1,&a2,&npts,&nburn,&ifunc);
  if (ss!=8) { printf("read error in control file, line 1\n"); return -1;}
  if (pinti("Kernel type", kern,1,3)) return -1;
  if (pintr("Smoothing fraction parameter", hfrac,1.0e-12,1.0)) return -1;
  if (pinti("Range type", krange,1,2)) return -1;
  if (a1 > a2) {printf("inverted range limits\n"); return -1;}
  if (krange == 2) {
     if (pintr("Range limit", a1,0.0,a2)) return -1;
     if (pintr("Range limit", a2,a1,1.0)) return -1;
  }
  if (pinti("Points produced", npts,1,10000000)) return -1;
  printf("                                  Kernel type:%10s\n",hkern[kern-1]);
  printf(" Window width fraction of interquartile range:%10.6f\n",hfrac);
  printf("%20s%10s ordinate range: %15.5e%15.5e\n",s,hrange[krange-1],a1,a2);
  printf("  Number of density ordinates written to file:%10d\n",npts);
  printf("                 Number of burn-in iterations:%10.d\n",nburn);
  printf("                         Function of interest:%10.d\n",ifunc);

/*open the simulation file and read 1st line */
  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, &npars);
  if (ss != 2) {printf("read error, first record of \ninput simulation ");
                printf("file %s", hfile); return -1;} 

/* error check */
  if (pinti("Burn-in iterations", nburn,0,niter)) return -1;
  if (pinti("Function of interest", ifunc,1,npars)) return -1;

/* open output file */
  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;}        

  nuse=niter-nburn;
  nskip=(nuse-1)/LD + 1;
  for (iburn=1; iburn <= nburn; iburn++) {
       if (readrec(npars,&jter,&wtlog,&pri,&pdat,par)) return -1;
  }
  
/* pass through file to get weights and values */
  m=0;
  for (iter=1; iter <= nuse; iter=iter+nskip) {
     m=m+1;
     if (readrec(npars,&jter,&wtlog,&pri,&pdat,par)) return -1;
     if (iter==1) wtmax=wtlog;
     if (wtlog > wtmax) wtmax=wtlog;
     wt[m-1]=wtlog;
     g[m-1]=par[ifunc-1];
  }
 
/* convert from log weights to weights */
  mtot=m;
  wsum=0.0;
  for (m=1; m <= mtot; m++) {
     wt[m-1]=exp(wt[m-1]-wtmax);
     nv1[m-1]=m;
     wsum=wsum+wt[m-1];
  }

/* sort weights and values in ascending order of values; get range points */
  bubble(mtot,g,g,nv1);

  b1=0.25*wsum;
  b2=0.75*wsum;
  r1=a1;
  r2=a2;
  if (krange==2) {r1=a1*wsum; r2=a2*wsum;}
  wrun=0;
  for (m=1; m <= mtot; m++) {
     w[m-1]=wt[nv1[m-1]-1];
     wrun1=wrun+w[m-1];
     if ((wrun<b1)&&(wrun1>=b1)) q1=g[m-1];
     if ((wrun<b2)&&(wrun1>=b2)) q2=g[m-1];
     if (krange==2) {
        if ((wrun<r1)&&(wrun1>=r1)) r1=g[m-1];
        if ((wrun<r2)&&(wrun1>=r2)) r2=g[m-1];
     }
     wrun=wrun1;
     w[m-1]=w[m-1]/wsum;
  }

/* compute density at ordinates and write to file */
  h=hfrac*(q2-q1);
  if (krange==2)
     printf("\n%31s Ordinate range:%15.5e%15.5e\n", s,r1,r2);
  printf("%33s Window width:%15.5e\n",s,h); 

  h1=1/h;     mbase=1;      z=r1;       zinc=(r2-r1)/(npts-1);
  ak=0.5;
  for (ipts=1; ipts <= npts; ipts++) {
     zlower=z-h; zupper=z+h;
     den=0;
     m=mbase;
     while (g[m-1] <= zlower) {
        mbase=mbase+1;
        m=mbase;
     }
     while ((m<=mtot)&&(g[m-1] < zupper)) {
        if (kern==2) ak=1.0-ab(h1*(g[m-1]-z)); 
        if (kern==3) ak=CON3*sq(1-sq(h1*(g[m-1]-z))); 
        den=den+w[m-1]*ak;
        m=m+1;
     }
     fprintf(fp2,"%20.8e%20.8e\n",z,h1*den);
     z=z+zinc; 
  }
  fclose(fp1); fclose(fp2);
  return 1;  
} /* main*/