/**************************************************************************
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 reversed Berk-Jones statistic.


User specifies the following:
    N, the number of data points  (3 <= N <= 500)
    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_rev
		    F
		    PofR_n_rev

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

#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;

  double n_ones=0.0;
  lambda=get_vector(nvalues+1);

  /* These are upper limits for lambda[3].  They depend on ALPHA. */
  if (ALPHA < 0.00001)
    lambda[2]=1.0985;
  if (ALPHA < 0.002 && ALPHA >= 0.00001)
    lambda[2]=1.09;
  if (ALPHA < 0.4 && ALPHA >= 0.002)
    lambda[2]=1.0;
  if (ALPHA >=0.4)
    lambda[2]=0.35;

  int n; 
  for(n=3; n <= nvalues; n++) {
  
  double answer;
  double LL, UL;
  double alpha = ALPHA;

  LL = 1e-16;
  UL = lambda[n-1];
  /* There are certain ALPHA values where we need a tighter bound upper bound
     on lambda[n] than lambda[n-1] */  
  if (ALPHA < 0.015 && n==4)
    UL=lambda[n-1]+0.28;
  if (ALPHA < 0.005 && n==5)
    UL=lambda[n-1]+0.2;
  if (ALPHA < 0.0005 && n==6)
    UL=lambda[n-1]+0.15;
  if (ALPHA < 0.00001 && n==7)
    UL=lambda[n-1]+0.2;
  if (ALPHA < 0.000001 && n==8)
    UL=lambda[n-1]+0.2;
  if (ALPHA < 0.0000001 && n==9)
    UL=lambda[n-1]+0.2;

  answer=zbrent(PofR_n_rev, n_ones, alpha, n, data, LL, UL, 1e-20);
  /*checks to see if the number of 1's in the upper band must be increased */
  if (answer < 0) {
    n_ones++;
    answer=zbrent(PofR_n_rev, n_ones, 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);


  /* the case for the first order statistic */
  x=data[1];
  b[1]=zbrent(K_rev,lambda[n],x,n,data,1.0/n,1.0,1e-20);

  int z;   /* the number of ones in the upper confidence bound */  
  z= (int) n_ones;
  j=n-z;

 /* the cases for the other order statistics */
  for (i=2; i<j; i++) {
    
    x = data[i-1];
    LL= EmpDF(x,n,data);
    UL= 1.0;
    b[i]=zbrent(K_rev,lambda[n],x,n,data,LL,UL,1e-20);
    
}
  x=data[i-1];
  LL=EmpDF(x,n,data);
  UL=1.0;
  b[j]=zbrent(K_rev,lambda[n],x,n,data,LL,UL,1e-20);

  for(i=(j+1) ;i<=n; i++){
    b[i]=1.0;
  }
  
  for (i=1; i<=n; i++){
    a[i]=1-b[n-i+1];
    
  }
  /*
  FILE *out;
  out=fopen("lambda_r.dat","w");
  
  for(i=3; i <= n; i++) {
    fprintf(out, "%lf \n",lambda[i]);
  }
  
  if (fclose(out) != 0)
    printf("Error in closing file %s\n", "lambda_r.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;
}