#include <math.h>
#include <imsl.h>
#include <stdlib.h>
#include <stdio.h>
#define LD 50000
#define LDPAR 300
#define LVAL 100

FILE *fp1;

extern double robeq(double f);
extern double phi(double f);
double ak;
main()
{
  char hfile[40], s[]=" ";
/*INT*/
  int nd, nburn, ifunc, i, niter, npars, nuse, nskip, ss, jter,
      m, mtot, id, sc=-1, inc=1;
  double wtlog, pri, pdat, par[LDPAR], wtmax, g[LD], wt[LD], dpar[LVAL],
         qup[LVAL], amu, asd, qdown[LVAL],  *f, beta[LVAL], aa;

/* read in and check user info */
  printf("%10sPrior density ratio class robustness\n\n",s);
  scanf("%d%d%d", &nd,&nburn,&ifunc);
  if (pinti("Number of d values", nd,1,LVAL)) return -1;
  if (plistr(stdin,nd,dpar)) return -1;
  printf("  Number of burn-in iterations:%10.d\n",nburn);
  printf("          Function of interest:%10.d\n",ifunc);

/*open the simulation file and read 1st line */

  ss=scanf("%s", &hfile);
  if (ss!=1) {printf("error reading file name\n"); return -1;}  
  fp1=fopen(hfile,"r");
  if (fp1==NULL) {printf("failure to open file %s for input\n", hfile);
                  return -1;}
  ss=fscanf(fp1,"%d %d", &niter, &npars);
  if (ss != 2) {printf("read error, first record of \ninput simulation ");
                printf("file %s", hfile); return -1;} 

/* error check */
 
  if (pinti("Burn-in iterations", nburn,0,niter)) return -1;
  if (pinti("Function of interest", ifunc,1,npars)) return -1;

/* skip the burn-in iterations */

  nuse=niter-nburn;
  nskip=(nuse-1)/LD + 1;
  for (i=1; i<= nburn; i++) {
       if (readrec(npars,&jter,&wtlog,&pri,&pdat,par)) return -1;
  }

/* pass through file to get weights and values */
  
  m=0;
  for (i=1; i<= nuse; i=i+nskip) {
     m=m+1;
     if (readrec(npars,&jter,&wtlog,&pri,&pdat,par)) return -1;
     if (i==1) wtmax=wtlog;
     if (wtlog > wtmax) wtmax=wtlog;
     wt[m-1]=wtlog;
     g[m-1]=par[ifunc-1];
  }

/* convert from log weights to weights  and find moments*/
  
  mtot=m;
  for (m=1; m <= mtot; m++) 
     wt[m-1]=exp(wt[m-1]-wtmax);
  
  rob1(mtot,wt,g,nd,dpar,qup,&amu,&asd);
  for (i=0; i<=mtot-1; i++)
     g[i]=-g[i];
  rob1(mtot,wt,g,nd,dpar,qdown,&aa,&asd);
  for (i=0; i<=nd-1; i++)
     qdown[i]=-qdown[i];

  for (i=0; i <= nd-1; i++) {
     ak=exp(dpar[i]);
     f=imsl_d_zeros_fcn(robeq,0);
     beta[i]=(ak-1)*phi(*f)/(ak+(1-ak)*imsl_d_normal_cdf(*f));
  }
  printf("%15sRobustness analysis of function%3d\n\n",s,ifunc);
  printf(" Posterior mean%15.6e   Stan dev%15.6e\n\n",amu,asd);
  printf("%24sExact",s);
  printf("            DeRobertis & Hartigan asymptotic\n");
  printf("%9sd    Lower bound    Upper bound",s);
  printf("    Lower bound    Upper bound\n");
  for (i=0; i <= nd-1; i++) 
     printf("%10.6f%15.6e%15.6e%15.6e%15.6e\n",dpar[i],qdown[i],qup[i],
            amu-asd*beta[i],amu+asd*beta[i]);
  fclose(fp1); 
  return 1;  
} /* main*/