/*      This function computes the coverage probability attached to
  n pairs of bounds on the ecdf.
	The algorithm is Noe's recursion as described in Shorack and
  Wellner pp. 362-263.
	Specifically:

	find P( a[i] < X(n) <= b[i] for 1<=i<=n)
	    =P( Fa_i < U(n) <= Fb_i for 1<=i<=n)

  when the X(n) are order statistics from a sample of size n from F.


	Symbols are the same as in Shorack and Wellner, including
  subscripts.  Therefore some 0 elements of arrays are not used, etc.

        This version uses dynamic allocation of storage, so
  there is no preset sample size.  In rare circumstances it
  fails to free up all the space it allocates.  First of all
  if the problem is too large to handle, it doesn't free up
  what it has already managed to allocate.  Second if the input
  is absurd (as determined by noe_compute_cgh) the input is
  rejected without freeing up already allocated space.

        Coded by Art Owen, in 1987, with modifications in
  1993 and 1994.  The code has been checked in the following
  ways: - it matches a numerical example in Shorack and Wellner
        - it agrees closely with published values when tested
          on some Kolmogorov-Smirnov bands.  (The published values
	  are the widths of KS bands required to get 95% coverage
	  for a particular n.  These values are usually only
	  given to two or three significant figures.)
  I have not noticed any problems with this code for n<=1000.
  For larger n (eg 1800 or more) I have seen unexplained odd
  behavior.  This may be due to the calling program or it may
  be due to accumulation of roundoff eror within this function.

  */

extern double *get_vector(), **get_matrix();
extern int    *get_ivector();

#include <stdio.h>
#include <math.h>
#include <stdlib.h>

/*
void qsort(base, nel, width, compar)
void *base;
size_t nel, width;
int (*compar)();
*/

int dcompar( a, b )
double *a, *b;
{
printf("  a %g b %g\n",*a,*b);
if(  *a < *b  )return -1;
if(  *a > *b  )return  1;
return 0;
}


double ecdf_cover( n, a, b, F )

int     n;
double  *a, *b, (*F)();   /* a, b of length n+1 */

{
double  *c, *p, *coef, ans, *oQ, *nQ;
int     *g, *h;
int     i, k, m, twonp2;

if(  n < 1  )return(1.0); /* 0.0 would also be sensible! */

twonp2 = 2*n+2;
c = get_vector( twonp2 ); 
p = get_vector( twonp2 ); 
oQ = get_vector( twonp2 );
nQ = get_vector( twonp2 );
g = get_ivector( twonp2 ); 
h = get_ivector( twonp2 ); 
coef = get_vector( n+1 );

if(  !c || !p || !oQ || !nQ || !g || !h || !coef  ){
  fprintf(stderr,"Not enough memory to handle n=%d.\n",n);
  /* Note: if e.g.  (c && !nQ), we haven't free'd c */
  return(-1.0);
}

if(  !noe_compute_cgh(  n, a, b, c, g, h )  )return(0.0); /*Doesn't free*/
noe_compute_p(    n, c, F, p );


oQ[0] = 1.0;
for(  m=1; m<=2*n; m++  ){
    /*  Compute relevant part of column m */
    for(  i=h[m+1]-1; i<= g[m-1]; i++  ){
	coef[i] = 1.0;
	for(  k=i-1;  k>=h[m]-1;  k-- ){
	    coef[k] = coef[k+1] * p[m] * ((double) k+1)/((double) i-k);
	    }
	for(  k=i;  k>=h[m]-1;  k-- )
	  coef[k] *= oQ[k];

	nQ[i] = 0.0;
	for(  k=h[m]-1; k<=i; k++ ){
	    nQ[i] += coef[k];
	    }
      }
    nQ[g[m-1]+1] = 0.0;
    for(  i=h[m+1]-1; i<= g[m-1]+1; i++  )
      oQ[i] = nQ[i];
    }

ans = nQ[n];

free_vector(  coef, n+1 );
free_ivector( h, twonp2 ); 
free_ivector( g, twonp2 ); 
free_vector( oQ, twonp2 );
free_vector( nQ, twonp2 );
free_vector( p, twonp2 ); 
free_vector( c, twonp2 ); 

return ans;
}


noe_compute_p(    n, c, F, p  )

/*      Compute the F probabilities of all the intervals */

int     n;
double  *c, (*F)(), *p;

{
int     m;

p[1] = (*F)(c[1]) /* - (*F)(c[0]) */  ;   /* c[0] = -infinity  */
for(  m=2; m<=2*n; m++  ){
    p[m] = (*F)(c[m]) - (*F)(c[m-1]);
    }
p[2*n+1] = 1.0 - (*F)(c[2*n]);            /* c[2n+1] = +infinity */

/*
for(  m=1; m<=2*n+1; m++  ){
    printf(" %g %g %g \n",c[m-1],c[m],p[m]);
    }
*/



}



noe_compute_cgh(  n, a, b, c, g, h )

/*      Compute c,g,h for Noe's recursion */

int     n;
double  *a, *b, *c;
int     *g, *h;

{
int     ia, ib, m;
int     i, nprobs;

/*      Check that a and b are increasing and that a <= b*/

nprobs = (a[1]>b[1]);
for(  i=2; i<=n; i++ ){
    nprobs += (a[i-1] > a[i]);
    nprobs += (b[i-1] > b[i]);
    nprobs += (a[i  ] > b[i]);
    /*
    if(a[i-1] > a[i])printf("%g=a[%d]>a[%d]=%g\n",a[i-1],i-1,i,a[i]);
    if(b[i-1] > b[i])printf("%g=b[%d]>b[%d]=%g\n",b[i-1],i-1,i,b[i]);
    if(a[i  ] > b[i])printf("%g=a[%d]>b[%d]=%g\n",a[i  ],i  ,i,b[i]);
    */
}

if(  nprobs  ){
  /*
    fprintf(stderr,"noe_compute_cgh improper ordering, nprobs=%d\n",nprobs);
  */
    return(0);  /* Problems detected */
    }

/*      Merge a and b into c    */

g[0] = 0;
h[0] = 1;
h[1] = 1;
a[n+1] = b[n]+1;
b[n+1] = a[n]+1;

for(  ia=ib=1; ia+ib-1<=2*n;  ){
    m = ia+ib-1;
    if(  a[ia] <  b[ib]  ){
	c[m]   = a[ia];
	g[m]   = g[m-1]+1;
	h[m+1] = h[m];
	ia++;
	}
    else{
	c[m]   = b[ib];
	g[m]   = g[m-1];
	h[m+1] = h[m]+1;
	ib++;
	}
    }

/*
printf("             0 %d \n", g[0]);

for(  m=1; m<=2*n; m++ ){
    printf(" %g %d %d %d %d %d \n", c[m], m, g[m], h[m]-1,
	h[m+1]-1,g[m-1]);
   }
printf("            %d %d\n",2*n+1,h[2*n+1]-1);
*/

return(1); /* It worked */



}