(* enter the following *)
nd=3;   (* number of d values *)
nburn=100;  (* burn in iterations *)
ifunc=3;     (*function of interest *)
dpar={0.69, 1, 2};    (* d values *)
hfile="../../../models/uvr1/sim"; (* posterior simulation file *)
outfile="out"; (* output log file name *)  

(* DO NOT MODIFY ANYTHING BELOW *)

Off[General::spell];
Off[General::spell1];
LD=50000; LDPAR=300; LVAL=100;
q=Table[0,{i,1,nd}];
<<Statistics`ContinuousDistributions`;
<<readrec.m; <<for.m; <<rob1.m; <<robeq.m;
OpenRead::noopen="Cannot open input file";

(* open files *)
Check[h=OpenRead[hfile], Quit[1]];
out=OpenWrite[outfile];                               

(* check user info *)
WriteString [out, "          Prior density ratio class robustness\n\n"];
If[nd<1 || nd>LVAL,
          WriteString[out,"Number of d values ", nd];
          WriteString[out,"not in range 1 to ", LVAL, "\n"];
          Quit[1];
  ];  
For[id=1, id <= nd, id++, 
    If[dpar[[id]]<=0,
       WriteString[out,"Element ", id, " of list ", dpar[[id]]];
       WriteString[out," is not positive.\n"];
       Quit[1];
      ];
   ];  
WriteString[out,"                 Number of burn-in iterations: ", PaddedForm[nburn,10], "\n"];
WriteString[out,"                         Function of interest: ", PaddedForm[ifunc,10],"\n"];
                                                                               


(* open the simulation file and read the first line *)
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[npars=Read[h,Number],
      WriteString[out,"read error in simulation input file, line 1 \n"];
      Quit[1]];
If[npars==EndOfFile,
   WriteString[out,"unexpected end of simulation input file, line 1 \n"];
   Quit[1];
  ];    

(* error check *)
If[nburn>niter || nburn<0,
   WriteString[out,"Burn-in iterations ", nburn, "not in range 0 to ", niter, "\n"];
   Quit[1];
  ];
If[ifunc>npars || ifunc<1,
   WriteString[out,"Function of interest ", ifunc, "not in range 1 to ", npars, "\n"];
   Quit[1];
  ];

nuse=niter-nburn;
nskip=Floor[(nuse-1)/LD + 1];
For[iburn=1, iburn <= nburn, iburn++, 
    readrec[h,npars];
   ];  (* for iburn=1 to nburn *)      

(* pass through file to get weights and values *)
m=0;
For[iter=1, iter <= nuse, iter=iter+nskip, 
    m=m+1;
    readrec[h,npars];
    If[iter==1,wtmax=wtlog;];
    If[wtlog>wtmax,wtmax=wtlog;];
    wt[m]=wtlog;
    g[m]=par[[ifunc]];
   ];  (* for iter=1 to nuse *)      
mtot=m;

(* Convert from log weights to weights *)
wt=Table[wt[i],{i,mtot}];
g=Table[g[i],{i,mtot}];
wt=Exp[wt-wtmax] //N;
rob1[mtot, wt, g, nd, dpar];
qup=Table[q[[i]],{i,1,nd}];
g=g*(-1);
rob1[mtot, wt, g, nd, dpar];
qdown=Table[-q[[i]],{i,1,nd}];
amu=amu*(-1);

ndist=NormalDistribution[0,1];
beta=Table[0,{i,1,nd}];
For[id=1, id<=nd, id++, 
    ak=Exp[dpar[[id]]];
    ff=FindRoot[robeq[x]==0,{x,0}] //N;
    f = x/.ff; 
    beta[[id]]=(ak-1)*phi[f]/(ak+(1-ak)*CDF[ndist,f]) //N;
   ]; 
WriteString[out,"               Robustness analysis of function ",PaddedForm[ifunc,3],"\n\n"];
WriteString[out," Posterior mean ",e[amu,15,6],"   Stan dev ",e[asd,15,6],"\n\n"];
WriteString[out,"                        Exact"];
WriteString[out,"            DeRobertis & Hartigan asymptotic\n"];
WriteString[out,"         d    Lower bound    Upper bound"];
WriteString[out,"    Lower bound    Upper bound\n"];
For[id=1, id<=nd, id++, 
    WriteString[out,PaddedForm[dpar[[id]],{10,6}],"    ",e[qdown[[id]],15,6]];
    WriteString[out,"    ",e[qup[[id]],15,6],"   ",e[amu-asd*beta[[id]],15,6]];
    WriteString[out,"    ",e[amu+asd*beta[[id]],15,6],"\n"];
   ]; 
Close[hfile];
Close[outfile];
Quit[1];