#include <math.h>
#define LD 100
#define LDK 100
#define LSET 9


FILE *fp1, *fp2; 
int npar;
double pri, par[250], v1[LD], a1[LD][LD], v2[LD], a2[LD][LD];
extern double gauk(int k, double *amu, double *h, int ldh, double *x);
extern double gaun(int k, double *h, int ldh);

main()
{

int i,j;
float d;
/* INT */
  int nburn, nwrite, ss, niter, nrange, nuse, 
      nk=0, nk2, ik, jk, jter, iter, inc=1, info,
      stat, which=2, iset, ldummy, lall, irange;
  char hfile[40], simfile[40], s[]=" "; 
/*DOUBLE */
  double tvar[LDK][LDK], tmu[LDK], wtsum, base=0, aa, wt, theta[LDK],
         wtlog=0, pdat, invwt, tpre[LDK][LDK], set1=LSET+1, dk,
         bound, p, q, adj[LSET], cut[LSET], top[LSET], den[LSET], size,
         ratio, g[LSET], theta1[LDK], theta2[LDK], add, dum=0; 

  init();
  scanf("%d %d", &nburn, &nwrite);
  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;}
  if (nwrite) {
     ss=scanf("%s", &hfile);
     if (ss!=1) {printf("error reading file name\n"); return -1;}
     fp2=fopen(hfile,"w");
     if (fp2==NULL) {printf("failure to open file %s for output\n", hfile); 
                     return -1;}
  }
  fscanf(fp1,"%d %d", &niter, &npar);
  if (nburn > niter) {
       printf("%d burn in requested but only %d avaliable \n", nburn, niter);
       return -1;
  }       
  nrange=niter; 
  repar0();
  nuse = niter - nburn;
  headerml(niter,nburn,nuse);
  wtsum=0;
  dset(npar, 0.0, tmu, 1);
  umoset(npar,npar,tvar,LDK);
  for (iter=1; iter <= niter; iter++) {
    if (readrec(npar,&jter,&wtlog,&pri,&pdat,&par)) return 1; 
    if (iter > nburn) {
       repar(theta, &nk);
       if (iter==nburn+1) nk2=nk/2;
       base = base+pri+pdat+nk2; 
       wt=exp(wtlog); 
       wtsum=wtsum+wt;
       for (ik=0; ik<=nk-1; ik++) {
          aa=wt*theta[ik]; tmu[ik]=tmu[ik]+aa;
          for (jk=ik; jk<=nk-1; jk++) 
             tvar[ik][jk]=tvar[ik][jk]+aa*theta[jk];
       }
    }
  }
  base=base/nuse;
  invwt=1/wtsum;
  dscal_(&nk,&invwt,tmu,&inc);
  for (ik=0; ik<=nk-1; ik++) 
     for (jk=ik; jk<=nk-1; jk++) 
        tvar[ik][jk]=tvar[ik][jk]/wtsum -tmu[ik]*tmu[jk];
  dfull(nk,tvar); 
  dpdck(tvar,LDK,nk,&info);
  if (info!=0) {printf("singular parameter vector variance\n"); return -1;} 
  dlinds(nk,tvar,LDK,tpre,LDK);
  dk=nk;
  for (iset=1; iset<=LSET; iset++){
     p=(LSET+1-iset)/set1;
     adj[iset-1]=1/p;
     q=1-p;
     cdfchi(&which,&p,&q,&cut[iset-1],&dk,&stat,&bound);
     cut[iset-1]=cut[iset-1]/2;
  }
  ldummy=lrange2(nk,theta,&lall);
/* Integrate multivariate over parameter space */
  if (!lall) {
     gau(nk,tmu,tpre,LDK,theta1,theta2);
     dset(LSET,0.0,den,1);
     dset(LSET,0.0,top,1);
     for (irange=1; irange<=nrange; irange++) {
        gaua(nk,tmu,tpre,LDK,theta1,theta2);
        size=-gauk(nk,tmu,tpre,LDK,theta1);
        add=0.0;
        if (lrange(nk,theta1,&lall)) add=0.5;
        if (lrange(nk,theta2,&lall)) add=add+0.5;
        for (iset=1; iset<=LSET; iset++) {
           if (size <= cut[iset]) {
              den[iset]=den[iset]+1;
              top[iset]=top[iset]+add;
           }
        } /* for iset */
     } /* for irange */
     for (iset=1; iset<=LSET; iset++) {
        if (top[iset]<0) top[iset]=-1;
        adj[iset]=adj[iset]*den[iset]/top[iset]; 
     } /* for iset */
  printf("\n Adjustment factors for coverage \n"); 
  for (iset=1; iset<=LSET; iset++) 
     printf("%16.7e", adj[iset-1]);
  printf("\n");
  } /* if !lall */
/* construct Gelfand-Dey posterior moment and print */
  rewind(fp1);
  fscanf(fp1,"%d %d", &niter, &npar);
  dset(LSET, 0.0, top, 1);
  dset(LSET, 0.0, den, 1);
  if (nwrite) fprintf(fp2, "%16d%16d\n", nuse, LSET);
  for (iter=1; iter <= niter; iter++) {
    if (readrec(npar,&jter,&wtlog,&pri,&pdat,&par)) return 1; 
    if (iter > nburn) {
       repar(theta, &nk);
       wt=exp(wtlog);
       size=-gauk(nk,tmu,tpre,LDK,theta);
       ratio=exp(base-size-pri-pdat);
/* loop for all LSET alternative values of p  */
       for (iset=0; iset<=LSET-1; iset++) {
          g[iset]=0;
          if (size <= cut[iset]) {
             g[iset]=adj[iset]*ratio;
             top[iset]=top[iset]+wt*g[iset];
             den[iset]=den[iset]+wt;
          }  /* if size >= cut[iset] */ 
       } /*for iset */
    if (nwrite) {
       fprintf(fp2,"%16d %16.7e %16.7e %16.7e\n", iter, wtlog,dum,dum);
       for (iset=0; iset<=LSET-1; iset++)
          fprintf(fp2,"%16.7e", g[iset]);
       fprintf(fp2,"\n");
    } /* if nwrite */
    } /*if iter > nburn */
  } /* for iter */
/* print marginal likelihood */
aa=base-gaun(nk,tpre,LDK);
printf("\n\n%10s Marginal Likelihood (ML) Computations\n\n",s);
printf("Version   Prob %8s Weighted Fractions Included %6s Log ML \n",s,s);
for (iset=0; iset<=LSET-1; iset++){
   printf("%5d %12.6f %23.6f ",iset+1, (LSET-iset)/set1, den[iset]/wtsum);
   printf("%24.7e\n",log(wtsum/top[iset])+aa);
}
printf("\n Inverse ML entries written to file must be scaled by \n");
printf("exp(%16.7e) to recover the correct inverse ML\n",-aa);
fclose(fp1);
fclose(fp2);
return 1;
} //main

/*
for (i=1; i<=LSET; i++) printf("%f  ", cut[i-1]);
printf("\n");
for (i=1; i<=nk; i++) { for (j=1; j<=nk; j++) printf("%f  ", tpre[i-1][j-1]);
printf("\n"); }
printf("%f\n",invwt);
for (i=1; i<=nk; i++) printf("%le\t", theta[i-1]);
printf("\n");
*/