/* 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 write anything below */
let LDPAR=300; 
let LD=50000;
let CON3=0.9375;
hk1="   Uniform";
hk2="Triangular";
hk3="  Biweight";
hr1="          ";
hr2="  Quantile";

output file=^outfile reset; 

/* check user info */
print " Preparation of posterior density for plotting\n";
if ((kern<1) or (kern > 3));
   format /rdn 4,0; print "Kernel type " kern "not in range 1 to 3\n"; end;
endif;
if ((hfrac<=(1.0e-12)) or (hfrac >= 1));
   format /len 16,7; 
   print "Smoothing fraction parameter " hfrac "not in range 1.0e-12 to 1\n"; 
   end;
endif;
if ((krange<1) or (krange > 2));
   format /rdn 4,0; print "Range type " krange "not in range 1 to 2\n"; end;
endif;
if (a1 > a2);
   print "Inverted range limits\n"; end;
endif;
if krange==2;
   if ((a1<=0) or (a1 >= a2));
      format /len 16,7; print "Range limit " a1 "not in range 0 to " a2 "\n";
      end;
   endif;
   if ((a2<=a1) or (a2 >= 1));
      format /len 16,7; print "Range limit " a2 "not in range " a1 " to 1.0\n";
      end;
   endif;
endif;
if ((npts<1) or (npts > 10000000));
   format /rdn 4,0; 
   print "Points produced " npts "not in range 1 to 10000000"; end;
endif;

/* write user info to log file */
if kern==1; hkern=hk1; elseif kern==2; hkern=hk2; else; hkern=hk3; endif;
if krange==1; hrange=hr1; else; hrange=hr2; endif;
print "                                  Kernel type: " hkern;
format /rd 10,6;
print " Window width fraction of interquartile range: " hfrac;
print "                    ";;
format /len 15,5;
print hrange " ordinate range: " a1 a2;
format /rdn 10,0;
print "  Number of density ordinates written to file: " npts;
print "                 Number of burn-in iterations: " nburn;
print "                         Function of interest: " ifunc;    


load in[]=^hfile;
niter=in[1];
npars=in[2];
if not (rows(in) == (niter*(4+npars) + 2));
   print "
          Error in simulation input file line 1 
          or unexpected end of file reached at line " (niter*npars +2) "
          or some of the data is not numeric \n";
   end;
endif;
if ((nburn<0) or (nburn > niter));
   format /rdn 4,0; 
   print "Burn-in iterations " nburn "not in range 1 to " niter "\n"; end;
endif;
if ((ifunc<1) or (ifunc > npars));
   format /rdn 4,0; 
   print "Function of interest" ifunc "not in range 1 to " npars "\n"; end;
endif;

nuse=niter-nburn;
nskip=floor((nuse-1)/LD+1);
in=in[3:rows(in)];
in=reshape(in, niter, npars+4);
jter=in[.,1];
mtot=floor((nuse-1)/nskip+1);
gw=zeros(mtot,2);
/*gw=selif(in,jter.>nburn .and nskip.*floor((jter+1)./nskip).==(jter+1));*/
gw[.,2]=selif(in[.,2],jter.>nburn .and nskip.*floor((jter+1)./nskip).==(jter+1));
gw[.,1]=selif(in[.,(4+ifunc)],jter.>nburn .and nskip.*floor((jter+1)./nskip).==(jter+1));
/* sort weights and values */
gw=sortc(gw,1);
wtmax=maxc(gw[.,2]);
gw[.,2]=exp(gw[.,2]-wtmax);
wsum=sumc(gw[.,2]);

/* get range points */
b1=0.25*wsum;         b2=0.75*wsum;
r1=a1;                r2=a2;
if krange==2; r1=a1*wsum; r2=a2*wsum; endif;
wrun1=cumsumc(gw[.,2]);
wrun=wrun1-gw[1,2];
q1=selif(gw[.,1], wrun .< b1 .and wrun1 .>= b1);
q2=selif(gw[.,1], wrun .< b2 .and wrun1 .>= b2);
if (krange==2);
   r1=selif(gw[.,1], wrun .< r1 .and wrun1 .>= r1);
   r2=selif(gw[.,1], wrun .< r2 .and wrun1 .>= r2);
endif;
gw[.,2]=gw[.,2]/wsum;

/* compute density at ordinates and write to file */
h=hfrac*(q2-q1);
format /len 15,5;
if (krange==2);
   print "\n                               Ordinate range: " r1 r2;
endif;
print "                                 Window width: " h;

output file=^pout reset; 
h1=1/h;   z=r1;   zinc=(r2-r1)/(npts-1);
format /len 20,8;
for ipts (1,npts,1);
   zlower=z-h;   zupper=z+h;
   /* put all the g(theta)'s that fall between in zlower and zupper in gz */
   gz=selif(gw[.,1], gw[.,1] .> zlower .and gw[.,1] .< zupper);
   wz=selif(gw[.,2], gw[.,1] .> zlower .and gw[.,1] .< zupper);
   if (kern==1);
      den=sumc(wz*0.5);
   elseif (kern==2);
      den=sumc( wz.*(1-abs(h1*(gz-z))) );
   elseif (kern==3);
      den=sumc( wz.*(CON3*(1-(h1*(gz-z))^2)^2) );
   endif;
   print z h1*den;
   z=z+zinc;
endfor;
   
end;