#enter the following
nburn<-100;  #burn in iterations
nwrite<-1;   #1 if a simulation file is to be written, 0 otherwise
hfile<-"sim"; #input data file
outfile<-"mlout"; #output file name
simout<-"mlsim";  #name of the simulation output file if nwrite=1


repar<- function()
#  This routine is used to reparameterize the model, if desired.
#  (Reparameterization is often desirable to make the form of the posterior
#  closer to a multivariate normal.  See documentation in directory above.)
#  The code in this file is used for the univariate regression model;
#  editing is required for new models.
#  If parameters are not transformed, it is still necessary to copy the
#  appropriate information from par into theta and return nk.
#  Inputs:
#  pri       Logarithm of prior distribution at current parameter vector
#  par       Current record read from posterior simulation file;
#            parameter values are typically contained within this record.
#  npar      Number of elements in par
#  Outputs:
#  theta     Transformed parameter vector
#  nk        Number of elements in transformed parameter vector
{
   nk<<-(npar+1)/2;
   theta<<-par[1:nk,];
   theta[nk,]<<-log(par[npar,]);
   pri<<-pri+theta[nk,];
};

repar0<- function()
#  This procedure is typically used to read from file information that
#  may be required to interpret par in repar().  repar0() is entered once,
#  before repar().
{
};

lrange<-function(nk,theta)
#  This function indicates whether or not a parameter vector (as
#  reparameterized in repar) is in the support of the prior
#  Inputs:
#  nk         Number of parameters in the vector
#  theta      Parameter vector
#  Outputs:
#  lall       Set 1 if theta in Euclidean nk-space is in the
#             support of the prior; otherwise set 0
#  lrange     Set 1 if theta is in the support of the prior, else 0
{ lall<<-T; T };


#do not write anything below
cat(file=outfile,'\n');
LD<-100; LDK<-100; LSET<-9;
ntaper<-c(4,8,15);


#script starts here
if (nburn<0) nburn_0; 
lwrite_T;
if (nwrite==0) 
  lwrite_F;

in1_scan(hfile);
niter<-in1[1];
npar<-in1[2];
if (nburn>niter)
  {
     cat(file=outfile,append=T,'burn in requested but only ',nburn);
     cat(file=outfile,append=T,'avaliable', niter);
     q();
  }    
if (length(in1) != (niter*(4+npar) + 2))
{
  cat(file=outfile, 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


nrange_niter;
repar0();
nuse<-niter-nburn;
cat(file=outfile, append=T,'Marginal likelihood analysis \n');
cat(file=outfile, append=T,'Total iterations available: ',niter, '\n');
cat(file=outfile, append=T, '        Burn-in iterations: ',nburn,'\n');
cat(file=outfile, append=T, '       Analysis iterations: ',nuse,'\n');                   
   # first pass
#extract wtlog, pri, pdat, par
wtlog_in1[2,];   pri_in1[3,];   pdat_in1[4,];   par_in1[5:(4+npar),]; 
#skip the burn in iterations:
wtlog_wtlog[-(1:nburn)];    pri_pri[-(1:nburn)]; 
pdat_pdat[-(1:nburn)];      par_par[,-(1:nburn)];  

repar();

wt_exp(wtlog);              wtsum_sum(wt);
base_(sum(pri)+sum(pdat)+(niter-nburn)*nk/2)/nuse;

multwt_function(x) {x*wt}
tmu_(apply(apply(theta,1,multwt),2,sum))/wtsum;
tvar_theta%*%(apply(theta,1,multwt));
tr_function(x) { y <- x; dim(y) <- c(length(x), 1); y }
tvar_tvar/wtsum-tr(tmu)%*%tmu;
tpre_solve(tvar);

#prepare normal density constants
set1_LSET+1;
dk_nk;
adj_(seq(LSET,1,by=-1))^(-1)*set1;
cut_ qchisq(adj^(-1),dk)/2; 

lall_F;

gau<-function(k,amu,var) {
   a1<<-t(chol(var)); #a1*a1'=h 
   v1_a1%*%rnorm(k,0,1);
   theta1<<-amu+t(v1);
   theta2<<-amu-t(v1);
};
gaua<-function(k,amu,var) {
   v1_a1%*%rnorm(k,0,1);
   theta1<<-amu+t(v1);
   theta2<<-amu-t(v1);
};
gauk<- function(theta) {-(theta-tmu)%*%tpre%*%t(theta-tmu)/2};

ldummy_lrange(nk,theta);   
if (!lall)
{
   gau(nk,tmu,tvar);
   top<-rep(0,LSET); den<-rep(0,LSET);
   th1_rep(0,nk*nrange); dim(th1)_c(nk,nrange);
   th2_th1;
   for (irange in 1:nrange) {
      gaua(nk,tmu,tvar);
      th1[,irange]_t(theta1);
      th2[,irange]_t(theta2);
   }
   add_ifelse(apply(th1,2,lrange),0.5,0)+ifelse(apply(th2,2,lrange),0.5,0);  
   size_-apply(th1,2,gauk);
   ftd<- function(i) {
      dd<<-sum(ifelse(cut[i]>=size, 1,0));
      tt<<-sum(ifelse(cut[i]>=size, add,0));
   }
   for (i in 1:LSET) {
      ftd(i); top[i]_tt; den[i]_dd;
   }
   top_ifelse(top>0, top, -1);
   adj_adj*den/top;
   cat(file=outfile,append=T,"\nAdjustment factors for coverage:\n");
   cat(file=outfile,append=T,adj);
} #!lall

if (lwrite)
  cat(file=simout, nuse, LSET,'\n');

#second pass
size_-apply(theta,2,gauk);
ratio_exp(base-size-pri-pdat);

findg<- function(i) {
   g<<-ifelse(cut[i]>=size,adj[i]*ratio,0);
   tt<<-sum(ifelse(cut[i]>=size, wt*g,0));
   dd<<-sum(ifelse(cut[i]>=size, wt,0));
}
det_ function(x)prod(eigen(x)$values)
aa_base+nk*log(2*pi)/2-log(abs(det(tpre)))/2;

top<-rep(0,LSET);
den<-rep(0,LSET);

cat(file=outfile,append=T,"\nMarginal Likelihood (ML) Computations:\n");
cat(file=outfile,append=T,"\n Version  ");
cat(file=outfile,append=T,"Prob    ");
cat(file=outfile,append=T,"Weighted fractions included  ");
cat(file=outfile,append=T,"Log ML    ");
if (lwrite) {  #initialize
   gg_rep(0, LSET*nuse);
   dim(gg)_c(LSET,nuse);
};
sp_function(i)
{
  findg(i); top[i]<<-tt; den[i]<<-dd;
  if (lwrite) gg[i,]<<-g;
  cat(file=outfile,append=T, "\n",i); 
  cat(file=outfile,append=T,"        ",(LSET+1-i)/set1); 
  cat(file=outfile,append=T,"        ",den[i]/wtsum); 
  cat(file=outfile,append=T, "                  ",log(wtsum/top[i])+aa,"\n"); 
}
for(i in 1:LSET) 
   sp(i);

if (lwrite) 
   for (i in 1:nuse) {
      cat(file=simout,append=T,i+nburn,"     ",wtlog[i], " 0 0\n"); 
      cat(file=simout, append=T,gg[,i], "\n");
   };

cat(file=outfile,append=T,"Inverse ML entries written to file must be \n");
cat(file=outfile,append=T,"must be scaled by exp(", -aa, ") to recover ");
cat(file=outfile,append=T,"correct inverse ML");


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