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

(* DO NOT MODIFY ANYTHING BELOW *)

Off[General::spell];
Off[General::spell1];
LDPAR=250;
nline=40;
ncrit=5;
pcrit={0.001, 0.01, 0.05, 0.1, 0.5};
<<readrec.m; <<for.m; 
<<Statistics`ContinuousDistributions`;
OpenRead::noopen="Cannot open input file";
Check[h=OpenRead[hfile], Quit[1]];
out=OpenWrite[outfile];                               



crit=Table[(0.5)*Quantile[ChiSquareDistribution[kpars],1-pcrit[[i]]],{i,ncrit}];
ng=Floor[niter/nline];



(* First pass *)
SetStreamPosition[hfile,0];
Check[niter=Read[h, Number],
      WriteString[out, "read error in simulation input file, line 1 \n"];
      Quit[1]];
If[niter==EndOfFile,
   WriteString[out,"unexpected end of simulation input file, line 1 \n"];
   Quit[1];
  ];
Check[npar=Read[h,Number],
      WriteString[out,"read error in simulation input file, line 1 \n"];
      Quit[1]];
If[npar==EndOfFile,
   WriteString[out,"unexpected end of simulation input file, line 1 \n"];
   Quit[1];
  ];    
If[niter <= 0 || npar <=0,
   WriteString[out,"illegal first record of simulation input file "];
   WriteString[out,hfile, niter, npar];
   Quit[1];
  ];    
If[npar > LDPAR,
   WriteString[out,"input simulation file "];
   WriteString[out,hfile,npar, " parameters exceed limit of ", LDPAR];
   Quit[1];
  ];    


For[miter=1, miter <= (nline*ng), miter++, 
    readrec[h,npar]; 
    If[miter==1, gmax=pdat; iter=jter; ];     
    If[pdat>gmax, gmax=pdat; iter=jter; ];     
   ];  (* for iter=1 to niter *)      

WriteString[out,"Largest log data pdf:",e[gmax,14,6],"\n"];
WriteString[out,"        at iteration:",iter,"\n"];

(* second pass *)
WriteString[out,"                                   "];
WriteString[out,"Fraction downhill from chisquare(",kpars,") at p=\n"]; 
WriteString[out," Group    Iterations   Max log data pdf "];
For[icrit=1, icrit <= ncrit, icrit++, 
    WriteString[out, PaddedForm[pcrit[[icrit]],{6,3}] ];
   ];
WriteString[out,"\n"];
SetStreamPosition[hfile,0];
niter=Read[h, Number];
npar=Read[h, Number];
For[line=1, line <= nline, line++, 
    nover=Table[0,{i,ncrit}];
    For[ig=1, ig <= ng, ig++, 
        readrec[h,npar]; 
        If[ig==1, iter1=jter; amax=pdat; ];
        If[ig==ng, iter2=jter;];
        If[pdat > amax, amax=pdat;];
        For[icrit=1, icrit <= ncrit, icrit++, 
            If[(gmax-pdat)>crit[[icrit]], nover[[icrit]]=nover[[icrit]]+1;];
           ];
       ];
     WriteString[out,PaddedForm[line,5], PaddedForm[iter1,6],"-",
                     PaddedForm[iter2,6],"       ",e[amax,18,6]," " ];
     nover=nover/ng //N;
     For[icrit=1, icrit <= ncrit, icrit++, 
         WriteString[out,PaddedForm[nover[[icrit]],{6,3}] ];
        ];
     WriteString[out,"\n"];
   ];

Quit[1];
Close[hfile];