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


%do not write anything below
LDPAR=250;
pcrit=[0.001 0.01 0.05 0.1 0.5];
ncrit=5;
nline=40;

%preliminaries
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;
hfile=fscanf(fin,'%s',1);
kpars=fscanf(fin,'%d',1);
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 (niter <=0) | (npar<=0);
   fprintf(fout, 'illegal first record of simulation input file  %s %5d %5d', hfile, niter, npar); 
   quit; 
end;
if (npar>LDPAR);
   fprintf(fout, 'input simulation file  %s %5d parameters exceed limit of %5d',hfile, npar, LDPAR); 
   quit; 
end;
ng=floor(niter/nline);
global df;
df=kpars;
for i=1:ncrit;
   crit(i)=0.5*chiinv(1-pcrit(i));
end;

%first pass
for miter=1:(nline*ng);
   [iter, wtlog, pri, pdat, g]=readrec(npar,in,fout);
   if miter==1;
      gmax=pdat; jter=iter;
   end;
   if (pdat>gmax);
      gmax=pdat; jter=iter;
   end;
end; %for iter
fprintf(fout,'Largest log data pdf:%14.6e\n',gmax);
fprintf(fout,'        at iteration:%14d\n',jter);

%second pass
s=' ';
fprintf(fout,'%35s Fraction downhill from chisquare(%3d ) at p=\n',s,kpars);
fprintf(fout,' Group    Iterations   Max log data pdf');
for icrit=1:ncrit;
   fprintf(fout,'%8.3f',pcrit(icrit));
end;
fprintf(fout,'\n');
fprintf(fout,'\n');
frewind(in);
niter=fscanf(in,'%d',1);   npar=fscanf(in,'%d',1);   
for line=1:nline;
   nover=zeros(1,ncrit);
   for ig=1:ng;
      [iter, wtlog, pri, pdat, g]=readrec(npar,in,fout);
      if ig==1; iter1=iter; end;
      if ig==ng; iter2=iter; end;
      if ig==1; amax=pdat; end;
      if (pdat > amax); amax=pdat; end;
      for icrit=1:ncrit;
         if ((gmax-pdat)>crit(icrit)); 
            nover(icrit)=nover(icrit)+1;
         end;
      end;
   end;
   fprintf(fout,'%6d%7d-%7d%18.6e',line,iter1,iter2,amax);
   for icrit=1:ncrit;
      fprintf(fout,'%8.3f',nover(icrit)/ng);
   end;
   fprintf(fout,'\n');
end;

fclose(in);
fclose(fout);
quit;