%enter the following
cin = 'in';      %control file
out = 'out';     %output log file


%do not write anything below
LD=50000;
LDPAR=300;
CON3=0.9375;
hk1='   Uniform'; hk2='Triangular'; hk3='  Biweight';
hr1='        '; hr2='Quantile';
s=' ';

%open the control file and log file
fin=fopen(cin);
if fin==-1; disp('failure to open file %s'); disp(cin); quit; end;    
fout=fopen(out,'w');
if fout==-1; disp('failure to open file'); disp(out); quit; end;

%read the control file
kern=fscanf(fin,'%d',1);     hfrac=fscanf(fin,'%f',1);
krange=fscanf(fin,'%d',1);   [a,r]=fscanf(fin,'%f',[1,2]);
npts=fscanf(fin,'%d',1);     nburn=fscanf(fin,'%d',1);
ifunc=fscanf(fin,'%d',1);    hfile=fscanf(fin,'%s',1);
outfile=fscanf(fin,'%s',1);

%open the input file and output file
fprintf(fout,' Preparation of posterior density for plotting\n');
in=fopen(hfile);
if in==-1; fprintf(fout, 'failure to open file %s', hfile); quit; end;
pout=fopen(outfile,'w');
if pout==-1; fprintf(fout,'failure to open file %s', outfile); quit; end;

%check user info
if (kern<1) | (kern>3);
   fprintf(fout,'**Kernel type %4d not in range 1 to 3\n',kern); quit;
end;
if (hfrac<=1.0e-12) | (hfrac>=1);
   fprintf(fout,'**Smoothing fraction parameter %16.7e not in range 1.0e-12 to 1.0\n',hfrac); quit;
end;
if (krange<1) | (krange>2);
   fprintf(fout,'**Range type %4d not in range 1 to 2\n',krange); quit;
end;
if (a(1)>a(2));
   fprintf(fout,'inverted range limits\n'); quit;
end;
if krange==2;
   if (a(1)<=0) | (a(1)>=a(2));
      fprintf(fout,'**Range limit %16.7e not in range 0 to %16.7e\n',a(1),a(2));
      quit;
   end;
   if (a(2)<=a(1)) | (a(2)>=1);
      fprintf(fout,'**Range limit %16.7e not in range %16.7e to 1.0\n',a(2),a(1));
      quit;
   end;
end;
if (npts<1) | (npts>10000000);
   fprintf(fout,'**Points produced %4d not in range 1 to 10000000\n',npts); quit;
end;

%write user info to log file
if kern==1; hkern=hk1; elseif kern==2; hkern=hk2; else hkern=hk3; end;
if krange==1 hrange=hr1; else hrange=hr2; end;
fprintf(fout,'                                  Kernel type:%10s\n',hkern);
fprintf(fout,' Window width fraction of interquartile range:%10.6f\n',hfrac);
fprintf(fout,'%20s%10s ordinate range: %15.5e%15.5e\n',s,hrange,a(1),a(2));
fprintf(fout,'  Number of density ordinates written to file:%10d\n',npts);
fprintf(fout,'                 Number of burn-in iterations:%10d\n',nburn);  
fprintf(fout,'                         Function of interest:%10d\n',ifunc);

%open simulation file and read
niter=fscanf(in,'%d',1);   errcheck(niter,fout,0);
npars=fscanf(in,'%d',1);    errcheck(npars,fout,0);

%error check
if (nburn<0) | (nburn>niter);
   fprintf(fout,'**Burn-in iterations %4d not in range 0 to niter',nburn,niter);
   quit;
end;
if (ifunc<1) | (ifunc>npars);
   fprintf(fout,'**Function of interest %4d not in range 1 to npars',ifunc,npars);
   quit;
end;

nuse=niter-nburn;
nskip=floor((nuse-1)/LD)+1;
for iburn=1:nburn;
   [jter, wtlog, pri, pdat, par]=readrec(npars,in,fout);
end;

%pass through file to get weights and values
m=0;
for iter=1:nskip:nuse;
   m=m+1;
   [jter, wtlog, pri, pdat, par]=readrec(npars,in,fout);
   if iter==1; wtmax=wtlog; end;
   if (wtlog>wtmax) wtmax=wtlog; end;
   wt(m)=wtlog;
   g(m)=par(ifunc);
end; %for iter

%convert from log weights to weigths
mtot=m;
wsum=0;
for m=1:mtot;
   wt(m)=exp(wt(m)-wtmax);
   wsum=wsum+wt(m);
end; %for m

%sort weights and values in ascending order of values; get range points

[g,nv1]=sort(g);
g(mtot+1)=0; %trick
b1=0.25*wsum;    b2=0.75*wsum;
r1=a(1);         r2=a(2);
if krange==2; r1=a(1)*wsum; r2=a(2)*wsum; end;
wrun=0;
for m=1:mtot;
   w(m)=wt(nv1(m));
   wrun1=wrun+w(m);
   if (wrun<b1) & (wrun1>=b1); q1=g(m); end;
   if (wrun<b2) & (wrun1>=b2); q2=g(m); end;
   if krange==2;
      if (wrun<r1) & (wrun1>=r1); r1=g(m); end;
      if (wrun<r2) & (wrun1>=r2); r2=g(m); end;
   end;
   wrun=wrun1;
   w(m)=w(m)/wsum;
end; %for m

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

h1=1/h;    mbase=1;     z=r1;   zinc=(r2-r1)/(npts-1);  ak=0.5;
for ipts=1:npts;
   zlower=z-h;    zupper=z+h;
   den=0;
   m=mbase;
   while (g(m)<=zlower);
      mbase=mbase+1;
      m=mbase;
   end; %while
   while (m<=mtot)&(g(m)<zupper);
      if kern==2; ak=1-abs(h1*(g(m)-z)); end;
      if kern==3; ak=CON3*(1-(h1*(g(m)-z))^2)^2; end;
      den=den+w(m)*ak;
      m=m+1;
   end;
   fprintf(pout,'%20.8e%20.8e\n',z,h1*den);
   z=z+zinc;
end; %for ipts
fclose(in);
fclose(fout);
fclose(pout);
quit;