#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 
outf_"out";# output log file name 
pout_"outgr"; # name of the numerical file for plotting 


#do not write anything below
LDPAR<-300; 
LD_50000;
CON3_0.9375;
hkern_c("   Uniform","Triangular","  Biweight");
hrange_c("          ","  Quantile"); 

#script starts here

#check user info

cat(file=outf,'\n'); cat(file=pout,'\n');
cat(file=outf,append=T, ' Preparation of posterior density for plotting\n');
if (kern<1 || kern >3) {
  cat(file=outf,append=T,'Kernel type ',kern, 'not in range 1 to 3\n');
  q();
}    
if (hfrac<=1.0e-12 || hfrac >=1) {
  cat(file=outf,append=T,'Smoothing fraction parameter ',hfrac);
  cat(file=outf,append=T, 'not in range 1.0e-12 to 1.0\n');
  q();
}    
if (krange<1 || krange >2) {
  cat(file=outf,append=T,'Range type ',krange, 'not in range 1 to 2\n');
  q();
}    
if (a1 > a2) {
   cat(file=outf,append=T, 'Inverted range limits\n');
   q();
}
if (krange ==2) {
   if (a1<=0 || a1 >=a2) {
     cat(file=outf,append=T,'Range limit ',a1);
     cat(file=outf,append=T, 'not in range 0 to ', a2,'\n');
     q();
   }    
   if (a2<=a1 || a2 >=1) {
     cat(file=outf,append=T,'Range limit ',a2);
     cat(file=outf,append=T, 'not in range ', a1, ' to 1\n');
     q();
   }    
}
if (npts<1 || npts >10000000) {
  cat(file=outf,append=T,'Points produced ',npts, 'not in range 1 to 10000000\n');
  q();
}    

#write user info to log file
cat(file=outf,append=T,'                                  Kernel type: ');
cat(file=outf,append=T, hkern[kern],'\n');
cat(file=outf,append=T,' Window width fraction of interquartile range: ');
cat(file=outf,append=T, hfrac,'\n');
cat(file=outf,append=T,'                    ');
cat(file=outf,append=T,hrange[krange],' ordinate range', a1, a2, '\n');
cat(file=outf,append=T,'  Number of density ordinates written to file: ');
cat(file=outf,append=T,npts,'\n');
cat(file=outf,append=T, '                Number of burn-in iterations: ');
cat(file=outf,append=T,nburn,'\n');
cat(file=outf,append=T, '                        Function of interest: ');
cat(file=outf,append=T,ifunc,'\n');
                                                                               


in1_scan(hfile);
niter<-in1[1];
npars<-in1[2];
if (nburn<0 || nburn >niter) {
  cat(file=outf,append=T,'Burn-in iterations ',nburn);
  cat(file=outf,append=T, 'not in range 0 to ', niter, '\n');
  q();
}    
if (ifunc<1 || ifunc >npars) {
  cat(file=outf,append=T,'Function of interest ',ifunc);
  cat(file=outf,append=T, 'not in range 1 to ', npars, '\n');
  q();
}    

if (length(in1) != (niter*(4+npars) + 2))
{
  cat(file=outf, 
      "unexpected end of simulation input file reached at line ",
      (niter*npars +2),"\n");
  q();
};



in1_in1[-(1:2)];
dim(in1)_c(npars+4,niter);     #reshape in1
nuse_niter-nburn;
nskip_floor((nuse-1)/LD+1);
#mtot_floor((nuse-1)/nskip+1);

#extract weights and values, order in ascending order of values
wt_in1[2,];   
g_in1[4+ifunc,]; 
iter_seq(nburn+1, niter, by=nskip);
g_g[iter];
wt_wt[iter];
nv_order(g,na.last=NA);
g_g[nv];
wt_wt[nv];

wtmax_max(wt);
wt_exp(wt-wtmax);
wsum_sum(wt);

#get range points
b1_0.25*wsum;     b2_0.75*wsum;
r1_a1;            r2_a2;
if (krange==2) {
   r1_a1*wsum;     r2_a2*wsum;
}
wrun1_cumsum(wt);
wrun_wrun1-wt[1];
q1_g[wrun<b1 & wrun1>=b1];
q2_g[wrun<b2 & wrun1>=b2];
if (krange==2) {
   r1_g[wrun<r1 & wrun1>=r1];
   r2_g[wrun<r2 & wrun1>=r2];
}
wt_wt/wsum;

#compute density at ordinates and write to file
h_hfrac*(q2-q1);
if (krange==2) {
   cat(file=outf,append=T,'\n                               ');
   cat(file=outf,append=T,'Ordinate range: ', r1, r2, '\n');
}
cat(file=outf,append=T,'                                 ');
cat(file=outf,append=T,' Window width: ',h,'\n');

h1_1/h;  zinc_(r2-r1)/(npts-1);
z_seq(r1,r2,by=zinc);
findden_function(zz){
   zlower_zz-h;
   zupper_zz+h;
   gz_g[g>zlower & g<zupper];
   wz_wt[g>zlower & g<zupper];
   if (kern==1) {
      h1*sum(wz*(0.5));
   }
   else if (kern==2) {
      h1*sum(wz*(1-abs(h1*(gz-zz))));
   }
   else if (kern==3) {
      h1*sum(wz*CON3*(1-(h1*(gz-zz))^2)^2);
   }
};
den_sapply(z,findden);
for (i in 1:length(z)) {
   cat(file=pout,append=T,z[i],'   ', den[i],'\n');
}


#clean up
#remove(objects());
q();