#enter the following
hfile<-"../../../models/uvr1/sim";  #simulation file
kpars<-6;    #number of parameters
outf<-"out";  #output file



#do not write anything below
LDPAR<-250; 
ncrit<-5;
nline<-40;
pcrit<-c(0.001, 0.01, 0.05, 0.1, 0.5);


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();
};

#preliminaries
crit_qchisq(1-pcrit,kpars)/2;
ng_floor(niter/nline);

in1_in1[-(1:2)];
dim(in1)_c(npar+4,niter);     #reshape in1

#extract iter and pdat
iter_in1[1,1:(nline*ng)];
pdat_in1[4,1:(nline*ng)];   
gmax_max(pdat);
jter_iter[gmax==pdat][1];
cat(file=outf,append=T," Largest log data pdf: ", gmax,"\n");
cat(file=outf,append=T,"         at iteration: ", jter,"\n");

#get summary statistics by group
cat(file=outf,append=T,"                                   ");
cat(file=outf,append=T,"Fraction downhill from chisquare(",kpars,") at p=\n");
cat(file=outf,append=T," Group    Iterations   Max log data pdf");
for(i in 1:ncrit) 
    cat(file=outf,append=T,"   ",pcrit[i]);
sp_function(i)
{
  nover_rep(0,5);
  iter1_(i-1)*ng+1;
  iter2_i*ng;
  x_pdat[iter>=iter1 & iter <=iter2];
  amax_max(x);
  y_gmax-x;
  for (j in 1:ncrit) 
     nover[j]_sum(y>crit[j]);
  cat(file=outf,append=T,i,iter1, iter2, amax);
  for(j in 1:ncrit) 
     cat(file=outf,append=T,"   ",nover[j]/ng);
  cat(file=outf,append=T,"\n");
}
for(i in 1:nline) 
   sp(i);
   
#clean up
#remove(objects());
q();