/**************************************************************************
Leah Jager
November 2003

This program allows interaction by the user to calculate the associated 
(1-alpha) confidence intervals associated with lambda(1-alpha) found by 
inverting the Berk-Jones statistic.


User specifies the following:
    N, the number of data points  (2 <= N <= 1000)
    outputFileName, the file where the confidence bands are to be written
    ALPHA, the confidence level  (1e-8 <= ALPHA < 1) 


Functions and their sources:

Owen in memory.c:  get_vector
                   get_ivector
		   free_vector
		   free_ivector

Owen in noe.c:  dcompar
                ecdf_cover
		noe_compute_p
                noe_compute_cgh

Numerical Recipes in C in vw-d-b-method.c:  zbrent
                                            nrerror

Me in functions.c:  EmpDF
                    K
		    F
		    PofR_n_BJ

***************************************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <malloc.h>
#include "memory.c"
#include "noe.c"
#include "vw-d-b-method.c"
#include "functions.c"

#define inputFileName "simdata.dat" /*simulated data used to find ci's */

/*User specifies the following: */
#define outputFileName "ci.dat"  /*the file where you want the ci's stored*/
#define ALPHA 0.05  /*your alpha value*/
#define N 20  /*the length of your data set*/

/**************************************************************************/
int main(void) {

  /* read in the data from a file  */
  
  FILE *in;
  int nvalues=N;
  int i,j;

  in=fopen(inputFileName,"r");
 
  double *data;
  data = get_vector(nvalues+1);
  if( !data )
    fprintf(stderr, "Not enough memory to handle n=%d", nvalues);
  data[0] = 0.0;
  data[nvalues+1] = 1.0;
  
  for(i=1; i <= nvalues; i++)
    fscanf(in, "%lf", &data[i]);

  if (fclose(in) != 0)
    printf("Error in closing file %s\n", inputFileName);
  

  double *lambda;
  lambda=get_vector(nvalues+1);

  /* These are upper limits for lambda[2].  They depend on ALPHA. */
  if (ALPHA >= 0.05)
    lambda[1]=3.7;
  if (ALPHA < 0.05)
    lambda[1]=15.0;

  int n; 
  for(n=2; n <= nvalues; n++) {
    
  double answer;
  double LL, UL;
  double alpha = ALPHA;
  
  UL = lambda[n-1];
  LL = 1e-16;
  /* There are certain ALPHA values where we need a tighter bound upper bound
     on lambda[n] than lambda[n-1] */
  if (ALPHA < 0.000001 && n==2)
    UL=15.0;
  if (ALPHA < 0.000001 && n==3)
    UL=10.0;
  
  answer=zbrent(PofR_n_BJ,1.0, alpha, n, data, LL, UL, 1e-20);

  lambda[n] = answer;
  }  
  
  /*calculates the ci's for N*/ 
  n=nvalues;

  double *a;
  double *b;
  double x;
  double LL, UL;
  a = get_vector(n+1);
  b = get_vector(n+1);
  if( !a || !b )
    fprintf(stderr, "Not enough memory to handle n=%d", n);

  b[1]=zbrent(K,lambda[n],data[0],n,data,1e-16,1.0-1e-16,1e-20);

  for (i=2; i<=n; i++) {
    x = data[i-1];

    LL= EmpDF(x,n,data);
    UL= 1.0-1e-16;
    
    b[i]=zbrent(K,lambda[n],x,n,data,LL,UL,1e-20);
  }

  for (i=1; i<=n; i++){
    a[i]=1-b[n-i+1];
   
  }

  /* prints out all lambda values up to lambda[n] */
  /* 
  FILE *out;
  out=fopen("lambda_bj.dat","w");
  
  for(i=2; i <= n; i++) {
    fprintf(out, "%lf \n",lambda[i]);
  }
  
  if (fclose(out) != 0)
    printf("Error in closing file %s\n", "lambda_bj.dat");
  */  

  /* writes confidence band to file and closes file */
  
  FILE *out;
  out=fopen(outputFileName,"w");
  
  fprintf(out, "a \t\t b \n");
  
  for(i=1; i <= n; i++) {
    fprintf(out, "%lf \t %lf \n", a[i],b[i]);
  }
  
  if (fclose(out) != 0)
    printf("Error in closing file %s\n", outputFileName);
  
  /*frees memory*/  
  
  free_vector(a, n+1);
  free_vector(b, n+1);
   
  free_vector(data, nvalues+1);
  free_vector(lambda,nvalues+1); 
  return 0;
}