/* enter the following */
hfile="../../../models/uvr1/sim";  /* simulation file */
kpars=6; /* number of parameters*/
outfile="out"; /* name of the output file */



      /* do not write anything below */
let LDPAR=250; 
output file=^outfile reset; 
let pcrit={0.001 0.01 0.05 0.1 0.5};
let ncrit=5;
let nline=40;

load in[]=^hfile;
niter=in[1];
npar=in[2];
if not (rows(in) == (niter*(4+npar) + 2));
   print "
          Error in simulation input file line 1 
          or unexpected end of file reached at line " (niter*npar +2) "
          or some of the data is not numeric \n";
   end;
endif;
if ((niter<=0) or (npar <= 0));
   print "Illegal first record of input simulation file" hfile niter npar;
endif;
if (npar > LDPAR);
   print "input simulation file " hfile  npar " parameters exceed limit of " LDPAR;
endif;

/* preliminaries */
crit=((0.5)*cdfchii(1-pcrit,kpars+zeros(1,ncrit)));
ng=floor(niter/nline);


/* get largest log data pdf*/
in=in[3:rows(in)];
in=reshape(in, niter, npar+4);
in=selif(in, in[.,1].<=ng*nline);
gmax=maxc(in[.,4]);
jter=maxindc(in[.,4]);
format /len 14,6;
print " Largest log data pdf: " gmax;
format /rdn 14,0;
print "         at iteration:" jter;

/* get summary statistics by group */
print "                                   ";;
format /rdn 3,0;
print "Fraction downhill from chisquare (" kpars ") at p=";
print " Group    Iterations   Max log data pdf";;
format /rd 7,3; 
for i (1,ncrit,1);
   print pcrit[i];;
endfor;
for i (1,nline,1);
   nover=zeros(1,ncrit);
   iter1=(i-1)*ng+1;
   iter2=i*ng;
   x=selif(in,in[.,1].>=iter1 .and in[.,1] .<=iter2);
   amax=maxc(x[.,4]);
   y=gmax-x[.,4];
   for j (1,ncrit,1);
      nover[j]=sumc(y.>crit[j]);
   endfor;
   format /rdn 6,0; print i;;
   format /rdn 7,0; print iter1 "-" iter2;;
   format /ren 18,6; print  amax;;
   format /rd 7,3;
   for j (1,ncrit,1);
      print nover[j]/ng;;
   endfor;
   print;
endfor;


end;