#enter the following
nd_3;        #  number of d values
nburn_100; # burn in iterations 
ifunc_3;   # function of interest 
dpar_c(0.690,1,2) # d values
hfile_"../../../models/uvr1/sim"; # posterior simulation file 
outf_"out";# output log file name 


#do not write anything below
LDPAR<-300; 
LD_50000;
LVAL_100;

#script starts here

#check user info

cat(file=outf,'\n');
cat(file=outf,append=T, '          Prior density ratio class robustness\n\n');
if (nd<1 || nd >LVAL) {
  cat(file=outf,append=T,'Number of d values ',nd, 'not in range 1 to ',LVAL, '\n');
  q();
}    
for (i in 1:nd) 
   if (dpar[i]<=0) {
     cat(file=outf,append=T,'Element ',i, ' of list', dpar[i]);
     cat(file=outf,append=T, ' is not positive.\n');
     q();
   }    
#write user info to log file
cat(file=outf,append=T, '                Number of burn-in iterations: ');
cat(file=outf,append=T,nburn,'\n');
cat(file=outf,append=T, '                        Function of interest: ');
cat(file=outf,append=T,ifunc,'\n');
                                                                               


in1_scan(hfile);
niter<-in1[1];
npars<-in1[2];
if (nburn<0 || nburn >niter) {
  cat(file=outf,append=T,'Burn-in iterations ',nburn);
  cat(file=outf,append=T, 'not in range 0 to ', niter, '\n');
  q();
}    
if (ifunc<1 || ifunc >npars) {
  cat(file=outf,append=T,'Function of interest ',ifunc);
  cat(file=outf,append=T, 'not in range 1 to ', npars, '\n');
  q();
}    

if (length(in1) != (niter*(4+npars) + 2))
{
  cat(file=outf, 
      "unexpected end of simulation input file reached at line ",
      (niter*npars +2),"\n");
  q();
};



in1_in1[-(1:2)];
dim(in1)_c(npars+4,niter);     #reshape in1
nuse_niter-nburn;
nskip_floor((nuse-1)/LD+1);
#mtot_floor((nuse-1)/nskip+1);

#extract weights and values, order in ascending order of values
wt_in1[2,];   
g_in1[4+ifunc,]; 
iter_seq(nburn+1, niter, by=nskip);
g_g[iter];
wt_wt[iter];
wtmax_max(wt);
wt_exp(wt-wtmax);
findq<-function(ak,d,t,g,wg1,w1)
{
   ql_(t+ak*(wg1-t))/(d+ak*(w1-d));
   gg_g[2:length(g)]
   gg[length(g)]_g[length(g)]-1;
   qq_ql[ql>=g & ql<=gg];
   qq[1]
};
rob1<-function(g,wt,nd,dpar)
{
  nv_order(g,na.last=NA); g_g[nv]; wt_wt[nv];
  wg1_wt*g;
  w1_sum(wt); 
  d_cumsum(wt); t_cumsum(wg1);
  wg2_sum(wg1*g);
  wg1_sum(wg1); 
  amu<<-wg1/w1; asd<<-sqrt(wg2/w1-amu^2);
  for (i in 1:nd)
     q[i]<<-findq(exp(dpar[i]),d,t,g,wg1,w1);
};
q_rep(0,nd);
rob1(g,wt,nd,dpar);
qup_q;
g_-g;
rob1(g,wt,nd,dpar);
qdown_-q;
amu_-amu; 

robeq_function(x)
{
  (ak-1)*(dnorm(x)+x*pnorm(x))-ak*x 
};

beta_rep(0,nd);
sparin_scan("table");
dim(sparin)_c(2,100);

for (i in 1:nd){
    if (dpar[i] <= sparin[1,1])
    {  
      f_dpar[i]*sparin[2,1]/sparin[1,1];
      beta[i]_f;                             
    };
    else if (dpar[i] <= sparin[1,100])
    {
       fu_sparin[2,sparin[1,]>=dpar[i]];
       du_sparin[1,sparin[1,]>=dpar[i]];
       fl_sparin[2,sparin[1,]<dpar[i]];
       dl_sparin[1,sparin[1,]<dpar[i]];
       ll_length(fl);
       f_fl[ll]+(dpar[i]-dl[ll])*(fu[1]-fl[ll])/(du[1]-dl[ll]);
       beta[i]_f;
    };
}

cat(file=outf,append=T,'               Robustness analysis of function ');
cat(file=outf,append=T,ifunc,'\n\n');
cat(file=outf,append=T,' Posterior mean    ',amu);
cat(file=outf,append=T,' Stan dev    ',asd,'\n\n');
cat(file=outf,append=T,'                       Exact');
cat(file=outf,append=T,'            DeRobertis & Hartigan asymptotic\n');
cat(file=outf,append=T,'          d    Lower bound    Upper bound');
cat(file=outf,append=T,'    Lower bound    Upper bound\n');
for (i in 1:nd){
   if (dpar[i] <= sparin[1,100])
   {
      cat(file=outf,append=T,dpar[i],"   ",qdown[i],"    ",qup[i]);
      cat(file=outf,append=T,"    ",amu-asd*beta[i],"    ",amu+asd*beta[i],'\n');
   }
   else
   {
      cat(file=outf,append=T,dpar[i],"   ",qdown[i],"    ",qup[i]);
      cat(file=outf,append=T, "  d value too large\n");
   }
}



#clean up
#remove(objects());
q();