#enter the following
hfile<-"sima";  #investigator's simulation file
outf<-"newpsmout";  #output file for error messages to be written 
simout<-"client1";  #client's simulation file (output)

#edit the function client below

client<-
function(npar, thetai, pri, nc, gc, priorc)
###########################################################
#  This subroutine accomplishes reweighting of posterior simulator output
#  and/or construction of new functions of interest.  The code given here
#  leaves the prior and functions of interest unchanged, and therefore
#  accomplishes nothing.  This code must be replaced by code specific to
#  the client's interests.
#  Inputs:
#  npar    Number of functions of interest in the investigator's posterior
#          simulation file.
#  thetai  Vector of dimension npar, containing the parameter values read from
#          the investigator's simulation file in the current record.
#  pri     Logarithm of the investigator's prior evaluated at thetai
#  Outputs:
#  nc      Number of functions of interest in the client's posterior simulation
#          file.
#  gc      Vector of dimension nc, containing the values of the client's
#          function of interest
#  priorc  Logarithm of the client's prior evaluated at thetai
#
###################################
{
   nc<<-npar;
   priorc<<-pri;
   gc<<-thetai;
};
###################################

#do not write anything below
LDPAR<-100; 

#script starts here

in1_scan(hfile);
niter<-in1[1];
npar<-in1[2];
cat(file=outf,'\n');
if (niter<=0 || npar <=0) {
  cat(file=outf,append=T,'illegal first record of input simulation file ',niter, npar);
  q();
}    
if (npar > LDPAR) {
  cat(file=outf,append=T,'input simulation file ',npar);
  cat(file=outf,append=T,' parameters exceed limit of ', LDPAR);
  q();
}    
if (length(in1) != (niter*(4+npar) + 2))
{
  cat(file=outf, append=T,
      "unexpected end of simulation input file reached at line ",
      (niter*npar +2),"\n");
  q();
};
in1_in1[-(1:2)];
dim(in1)_c(npar+4,niter);     #reshape in1

#extract wtlog, pri, pdat, par
wtlog_in1[2,];   pri_in1[3,];   pdat_in1[4,];   thetai_in1[5:(4+npar),]; 

client(npar, thetai, pri, nc, gc, priorc);
if (nc<1 || nc >LDPAR) {
  cat(file=outf,append=T,'Clients function of interest ',nc);
  cat(file=outf,append=T,' not in range 1 to ',LDPAR);
  q();
}    
wtlog_wtlog+priorc-pri;
pri_priorc;

cat(file=simout,niter, nc, "\n");
for (i in 1:niter) {
   cat(file=simout,append=T,i,"   ",wtlog[i], "   ", pri[i], "   ", pdat[i], "\n"); 
   cat(file=simout, append=T,gc[,i], "\n");
};


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