/* Integrate3.c   Last Modified 12/01/99 */
/* Use I.S. to integrate a function      */
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>

void main(int argc, char *argv[ ])
{
  int numiter, tdf, n=5;
  double *theta;
  /*posterior mode*/
  double mu[5]={0.939373344604329,0.687756583850085,0.746314362940164,
		0.535995464487262,0.658348448932649};
  /*sigma=inv information*/
  double S[5][5]={{0.00323236484115814,8.92422821613215e-06,-0.00204938393705854,
		   -0.00219592090303282,-0.00186224855423242},
		  {8.92422821613215e-06,0.00121940565631849,9.63867549999951e-05,
		   -1.58367399748023e-05,1.18084750158351e-05},
		  {-0.00204938393705854,9.63867549999951e-05,0.0029697813334756,
		   0.00153095797404572,0.00127434923503003},
		  {-0.00219592090303282,-1.58367399748023e-05,0.00153095797404572,
		   0.00358790969204995,0.00127273333926394},
		  {-0.00186224855423242,1.18084750158351e-05,0.00127434923503003,
		   0.00127273333926394,0.00255465829939452}};
  double *C, *K, *Sinv;
  FILE *outfile;
  double isw, nc=0.0, det=1.0;
  double *mean;
  double **matalloc(int das,int dbs);
  void rant(double *thv, double *muv, double *Cm, int ds, int dfs);
  double func3(double *x);
  double tdens(double *xv, double *muv, double *Cinvm, int ds, 
	       int dfs, double dets);
  void dpotrf_(char *ul, int *nd, double *a, int *lda, int *info);
  void dpotri_(char *ul, int *nd, double *a, int *lda, int *info);
  double *vecalloc(int nr);
  int  j, k, info;
  char uplo[]="U";

  if (argc!=4){
    printf("\n");
    printf("Usage:  %s numiter tdf outfile\n",argv[0]);
  }
  else if ((!(numiter=atoi(argv[1])))||(numiter<1) ){
    printf("\n");
    printf("%s:  Error:  second command line argument, %i, must be an integer (>0)\n",
	   argv[0],numiter);
    printf("  giving the number of iterations.\n");
    printf("\n");
    exit(-1);
  }
  else if ((!(tdf=atoi(argv[2])))||(tdf<3)){
    printf("\n");
    printf("%s:  Error:  third command line argument must be an integer (>2)\n",argv[0]);
    printf("  giving the number of degrees of freedom for the t distribution.\n");
    printf("\n");
    exit(-1);
  }
  else if (!(outfile=fopen(argv[3],"w"))){
    printf("\n");
    printf("%s:  Error opening output file %s.\n",argv[0],argv[3]);
    printf("\n");
    exit(-1);
  }
  else{
    C=vecalloc(n*n);
    for (j=0;j<n;j++){
      for (k=0;k<n;k++){
	C[n*j+k]=S[j][k];
	/* printf("%6.4e ",C[n*j+k]); */
      }
      /* printf("\n"); */
    }
    /* Cholesky square root of cov matrix */
    dpotrf_(uplo,&n,&C[0],&n,&info);
    if (info==0){
      printf("Note:  Cholesky routine exited sucessfully\n");
    }
    else{
      printf("ERROR:  Cholesky routine exited with error %i\n",info);
      printf("  see man page for dpotrf\n");
      exit(-1);
    }
    Sinv=vecalloc(n*n);
    for (j=0;j<n;j++){
      for (k=0;k<n;k++){
	Sinv[n*j+k]=C[n*j+k];
	/*printf("%6.4e ",Sinv[n*j+k]); */
      }
      /*printf("\n");*/
    }      
    /* Inverse of Cov matrix using its Cholesky square root*/
    dpotri_(uplo,&n,&Sinv[0],&n,&info);
    if (info==0){
      printf("Note:  inversion routine exited sucessfully\n");
    }
    else{
      printf("ERROR:  inversion routine exited with error %i\n",info);
      printf("  see man page for dpotri\n");
      exit(-1);
    }
    K=vecalloc(n*n);
    for (j=0;j<n;j++){
      for (k=0;k<n;k++){
	K[n*j+k]=Sinv[n*j+k];
	/* printf("%6.4e ",K[n*j+k]);*/
      }
      /*printf("\n");*/
    }      
    /* Cholesky square root of inverse of cov matrix */
    dpotrf_(uplo,&n,&K[0],&n,&info);
    if (info==0){
      printf("Note:  Cholesky routine (call 2) exited sucessfully\n");
    }
    else{
      printf("ERROR:  Cholesky routine (call 2) exited with error %i\n",info);
      printf("  see man page for dpotrf\n");
      exit(-1);
    }
    for (j=0;j<n;j++) det=det*K[n*j+j];
    mean=vecalloc(n);
    theta=vecalloc(n);
    for (k=1;k<=numiter;k++){
      rant(theta,mu,C,n,tdf);
      isw=func3(theta)/tdens(theta,mu,K,n,tdf,det);
      nc=nc+isw;
      for (j=0;j<n;j++) mean[j]=mean[j]+isw*theta[j];
      if (k%1000==0) fprintf(outfile,"%6.4e\n",isw);
    }
    nc=nc/((double)numiter);
    for (j=0;j<n;j++) mean[j]=mean[j]/(nc*((double)numiter));
    fclose(outfile);
    printf("Normalizing Constant=%6.4e\n",nc);
    for (j=0;j<n;j++) printf("Mean[%i]=%6.4e\n",j+1,mean[j]);
  }
}