%enter the following
cin = 'mlin';      %control file
out = 'mlout';     %output file


%do not write anything below
LD=100; LDK=100; LSET=9;
fin=fopen(cin);
if fin==-1; disp('failure to open file'); disp(cin); quit; end;
fout=fopen(out,'w');
if fout==-1; disp('failure to open file'); disp(out); quit; end;
nburn=fscanf(fin,'%d',1);
nwrite=fscanf(fin,'%d',1);
hfile=fscanf(fin,'%s',1);
if nwrite; 
   simfile=fscanf(fin,'%s',1); simout=fopen(simfile,'w');
   if simout==-1; disp('failure to open file'); disp(simfile); quit; end;
end;
if nburn<0; nburn=0; end;

in=fopen(hfile);
if in==-1; fprintf(fout, 'failure to open file %s', hfile); quit; end;
niter=fscanf(in,'%d',1);   errcheck(niter,fout,0);
npar=fscanf(in,'%d',1);    errcheck(npar,fout,0);
if (nburn>niter);
   fprintf(fout,'%3d burn in requested but only ',nburn);
   fprintf(fout,'%3d avaliable', niter);
   quit;
end;
nrange=niter; 
repar0;
nuse=niter-nburn;
header1(fout,niter,nburn,nuse); 

%first pass
wtsum=0; tmu=zeros(1,npar); tvar=zeros(npar); base=0; 
for iter=1:nburn; readrec(npar,in,fout); end;
for iter=nburn+1:niter;
   [jter, wtlog, pri, pdat, par]=readrec(npar,in,fout);
   [nk,theta, pri]=repar(par,npar,pri);
   base=base+pri+pdat+nk/2; wt=exp(wtlog); wtsum=wtsum+wt;
   for ik=1:nk;
      aa=wt*theta(ik); tmu(ik)=tmu(ik)+aa;
      for jk=ik:nk;
         tvar(ik,jk)=tvar(ik,jk)+aa*theta(jk);
      end;
   end; %for ik
end; %for iter


tmu=tmu(1:nk);
tvar=tvar(1:nk,1:nk);
base=base/nuse; tmu=tmu/wtsum; 
tvar=tvar+tvar'-diag(diag(tvar));
tvar=tvar/wtsum-tmu'*tmu; tpre=inv(tvar);

%prepare normal density constants
set1=LSET+1; global df; df=nk; adj=(LSET:-1:1).^(-1)*10; 
for i=1:LSET
   cut(i)=0.5*chiinv(adj(i)^(-1));
end;
global lall; lall=0; ldummy=lrange(nk,theta);
if ~lall; 
  gau(nk,tmu,tpre); top=zeros(1,LSET); den=zeros(1,LSET);
  for irange=1:nrange;
     [theta1, theta2]=gaua(nk,tmu,tpre);
     size=(0.5)*(theta1-tmu)*tpre*(theta1-tmu)'; add=0;
     if lrange(nk,theta1); add=.5; end;
     if lrange(nk,theta2); add=add+.5; end;     
     for iset=1:LSET;
        if size<=cut(iset);
           den(iset)=den(iset)+1; top(iset)=top(iset)+add;
        end;
     end;
  end; %for irange
  for iset=1:LSET;
     if top(iset)<=0; top(iset)=-1; end;
  end;
  fprintf(fout, ' Adjustment factors for coverage\n');
  fprintf(fout, '%16.7e',adj); fprintf(fout,'\n');
end; %if ~lall

%second pass
frewind(in); niter=fscanf(in,'%d',1); npar=fscanf(in,'%d',1);
top=zeros(1,LSET); den=zeros(1,LSET);
if nwrite; fprintf(simout, '%5d %16d\n', nuse, LSET); end;
wtsum=0;
for iter=1:nburn; readrec(npar,in,fout); end;
for iter=nburn+1:niter;
   [jter, wtlog, pri, pdat, par]=readrec(npar,in,fout);
   [nk,theta, pri]=repar(par,npar,pri);
   wt=exp(wtlog); wtsum=wtsum+wt; size=(0.5)*(theta-tmu)*tpre*(theta-tmu)';
   ratio=exp(base-size-pri-pdat);
   for iset=1:LSET;
      g(iset)=0;
      if size<=cut(iset); 
         g(iset)=adj(iset)*ratio; top(iset)=top(iset)+wt*g(iset);
         den(iset)=den(iset)+wt;
      end; %if
   end; %for
   if nwrite; writerec(simout,iter,wtlog,g); end;
end;
aa=base+(0.5)*nk*log(2*pi)-(0.5)*log(abs(det(tpre)));
%continue from here
fprintf(fout, '      Marginal Likelihood (ML) Computations\n'); 
fprintf(fout, 'Version    Prob     Weighted Fractions Included    Log Ml\n');
for iset=1:LSET;
   fprintf(fout, '%3d',iset);
   fprintf(fout, '%9.6f',((LSET+1-iset)/set1));
   fprintf(fout, '%25.6f',den(iset)/wtsum);
   fprintf(fout, '%25.7e \n',(log(wtsum/top(iset))+aa));
end;
fprintf(fout, 'Inverse entries written to file must be scaled by\n');
fprintf(fout, 'exp(');
fprintf(fout, '%10.7e',-aa);
fprintf(fout, ') to recover correct inverse ML');

fclose(fin);
fclose(fout);
fclose(simout);
quit;