(* enter the following *)
kern=3;     (* type of kernel to use *)
hfrac=0.5;  (* window width as a fraction of interquartile range *)
krange=2;   (* range code for values of g(theta) *)
a1=0.001;
a2=0.999;   
(* if krange=1, absolute plotting is specified and ordinates from a1 to
   a2 will be plotted,
   if krange=2, quantile plotting is specified and ordinates from qa1 to
   qa2 will be plotted, where qa1 and qa2 are the a1'th and a2'th
   posterior quantiles of g(theta) respectively *)
npts=200;   (* number of ordinates to plot *)
nburn=100;  (* burn in iterations *)
ifunc=1;     (*function of interest *)
hfile="sim"; (* posterior simulation file *)
outfile="out"; (* output log file name *)  
pout="outgr"; (* name of the numerical file for plotting *)

(* DO NOT MODIFY ANYTHING BELOW *)

Off[General::spell];
Off[General::spell1];
LD=50000; LDPAR=300; CON3=0.9375;
Array[g,LD]; Array[wt,LD]; 
hkern={"  Uniform", "Triangular", "Biweight"};
hrange={"          ", "  Quantile"};
<<readrec.m; <<for.m; 
OpenRead::noopen="Cannot open input file";

(* open files *)
Check[h=OpenRead[hfile], Quit[1]];
out=OpenWrite[outfile];                               
outs=OpenWrite[pout];                               

(* check user info *)
WriteString [out, "Prepation of posterior density for plotting\n"];
If[kern<1 || kern>3,
          WriteString[out,"Kernel type ", kern, "not in range 1 to 3\n"];
          Quit[1];
  ];  
If[hfrac<=10^(-12) || hfrac>=1,
          WriteString[out,"Smoothing fraction parameter ", e[hfrac,16,7]];
          WriteString[out, "not in range 1.0e-12 to 1.0\n"];
          Quit[1];
  ];  
If[krange<1 || krange>2,
          WriteString[out,"Range type ", krange, "not in range 1 to 2\n"];
          Quit[1];
  ];  
If[a1>a2,
          WriteString[out,"Inverted range limits\n"];
          Quit[1];
  ];  
If[krange==2,
       If[a1<=0 || a1>=a2,
          WriteString[out,"Range limit ", e[a1,16,7]];
          WriteString[out," not in range 0 to ", e[a2,16,7],"\n"];
          Quit[1];
         ];
       If[a2<=a1 || a2>=1,
          WriteString[out,"Range limit ", e[a2,16,7]];
          WriteString[out," not in range ", e[a1,16,7],"to 1.0 \n"];
          Quit[1];
         ];
     ];
If[npts<1 || npts>10000000,
          WriteString[out,"Points produced ", npts];
          WriteString[out,"not in range 1 to 10000000\n"];
          Quit[1];
  ];  

(* Write log file *)
WriteString[out,"                                  Kernel type: ",hkern[[kern]],"\n"];
WriteString[out," Window width fraction of interquartile range: ", hfrac, "\n"];
WriteString[out,"                    "];
WriteString[out,hrange[[krange]], " ordinate range: ",e[a1,15,5], " ", e[a2,15,5],"\n"];
WriteString[out,"  Number of density ordinates written to file: ", npts, "\n"]; 
WriteString[out,"                 Number of burn-in iterations: ", nburn, "\n"];
WriteString[out,"                         Function of interest: ", ifunc,"\n"];


(* open the simulation file and read the first line *)
(*SetStreamPosition[hfile,0];*)
Check[niter=Read[h, Number],
      WriteString[out, "read error in simulation input file, line 1 \n"];
      Quit[1]];
If[niter==EndOfFile,
   WriteString[out,"unexpected end of simulation input file, line 1 \n"];
   Quit[1];
  ];
Check[npars=Read[h,Number],
      WriteString[out,"read error in simulation input file, line 1 \n"];
      Quit[1]];
If[npars==EndOfFile,
   WriteString[out,"unexpected end of simulation input file, line 1 \n"];
   Quit[1];
  ];    

(* error check *)
If[nburn>niter || nburn<0,
   WriteString[out,"Burn-in iterations ", nburn, "not in range 0 to ", niter, "\n"];
   Quit[1];
  ];
If[ifunc>npars || ifunc<1,
   WriteString[out,"Function of interest ", ifunc, "not in range 1 to ", npars, "\n"];
   Quit[1];
  ];

nuse=niter-nburn;
nskip=Floor[(nuse-1)/LD + 1];
For[iburn=1, iburn <= nburn, iburn++, 
    readrec[h,npars];
   ];  (* for iburn=1 to nburn *)      

(* pass through file to get weights and values *)
m=0;
For[iter=1, iter <= nuse, iter=iter+nskip, 
    m=m+1;
    readrec[h,npars];
    If[iter==1,wtmax=wtlog;];
    If[wtlog>wtmax,wtmax=wtlog;];
    wt[m]=wtlog;
    g[m]=par[[ifunc]];
   ];  (* for iter=1 to nuse *)      
mtot=m;

(* Convert from log weights to weights *)
wt=Table[wt[i],{i,mtot}];
wt=Exp[wt-wtmax] //N;
wsum=Sum[wt[[i]],{i,1,mtot}];

(* sort weights and values in ascending order of values *)
gsorted=Sort[Table[{g[i],wt[[i]]},{i,mtot}]];
g=Table[gsorted[[i,1]],{i,mtot}];
w=Table[gsorted[[i,2]],{i,mtot}];

(* get range points *)
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++,
    wrun1=wrun+w[[m]];
    If[(wrun<b1)&&(wrun1>=b1), q1=g[[m]];];
    If[(wrun<b2)&&(wrun1>=b2), q2=g[[m]];];
    If[krange==2,
       If[(wrun<r1)&&(wrun1>=r1), r1=g[[m]];];
       If[(wrun<r2)&&(wrun1>=r2), r2=g[[m]];];
      ];
    wrun=wrun1;
   ]; (* for m=1 to mtot *)
w=w/wsum;

(* compute density at ordinates and write to file *)
h=hfrac*(q2-q1);
If[krange==2,
   WriteString[out,"\n                               "];
   WriteString[out,"Ordinate range: ",e[r1,15,5], " ", e[r2,15,5],"\n"];
  ];
WriteString[out,"                                 "];
WriteString[out,"Window width: ",e[h,15,5],"\n"];

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]] <= zlower),
          mbase=mbase+1;
          m=mbase;
         ]; 
    While[(m<=mtot && g[[m]]<zupper),
          If[kern==2, ak=1-Abs[h1*(g[[m]]-z)]];
          If[kern==3, ak=(1-(h1*(g[[m]]-z))^2)^2*CON3];
          den=den+w[[m]]*ak; 
          m=m+1; 
         ];
    hh1=h1*den;
    WriteString[outs,e[z,20,8],"   ",e[hh1,20,8],"\n"];
    z=z+zinc;
   ]; (* for ipts *)

Close[hfile];
Close[outfile];
Close[pout];
Quit[1];