/*
 *  bw_icw.c by A. Pepelyshev  Copyright (C) 2011
 *  This version distributed under GPL (version 2 or later)
 */

/*
====         ICV bandwidth selection rule              ====
  Savchuk, O. Y., Hart, J. D., Sheather, S. J. (2010) 
  Indirect cross-validation for density estimation.
  J. Amer. Statist. Assoc. 105, no. 489, 415–423.
*/

#include <math.h>

#ifndef max
#  define max(a,b) ((a) > (b) ? (a) : (b))
#  define min(a,b) ((a) < (b) ? (a) : (b))
#endif

#define abso(a) ((a) > 0 ? (a) : -(a) )


#if !defined(PI)
#  define PI 3.14159265
#endif

#define DELMAX 1000
/* Avoid slow and possibly error-producing underflows by cutting off at
   plus/minus sqrt(DELMAX) std deviations */


#define kernL(x,a,s) ( (1+(a))*exp(-(x)*(x)*0.5)-(a)/(s)*exp(-0.5*(x)*(x)/((s)*(s))) )
#define kernLL(x,a,s) ( (1+(a))*(1+(a))/sqrt(2)*exp(-(x)*(x)*0.25) - 2*(a)*(1+(a))/sqrt(1+(s)*(s))*exp(-0.5*(x)*(x)/(1+(s)*(s))) + (a)*(a)/((s)*sqrt(2))*exp(-0.25*(x)*(x)/((s)*(s))) )

 
void
band_icv_bin(int *n, int *nb, double *d, int *x, double *h, double *u, double *a, double *s)
{
    int   i, nn = (*n), nbin = *nb;
    double delta, hh = (*h), sum, term;

    sum = 0.0;
    for (i = 0; i < nbin; i++) {
	delta = i * (*d) / hh;
	delta *= delta;
	if (delta >= DELMAX) break;
	term = kernLL(delta,*a,*s) - 2*kernL(delta,*a,*s);
	sum += term * x[i];
    }
    *u = kernLL(0,*a,*s) / (2 * nn * hh * sqrt(PI)) + sum / (nn * nn * hh * sqrt(PI));
}

void
band_icv_first_local_min(int *n, int *nb, double *d, int *x, double *hlow, double *hupp, double *a, double *s)
{
 int i=-1, Do=1;
 double h1=(*hlow), h2=(*hupp), h, y1=0, y2=0, y3=0;
 while( (i<100) && Do )
 {
  i++;
  h = h1+i*0.01*(h2-h1); 
  band_icv_bin(n, nb, d, x, (&h), (&y3), a, s);
  if(i>=2)
  {
   if( (y2<y3) && (y2<y1) )
   {
     *hlow = h-2*0.01*(h2-h1);
     *hupp = h;
     Do = 0;
   }
  }
  y1 = y2;
  y2 = y3;
 }
}

void
band_ddiff_bin(int *n, int *nb, double *d, double *x, int *cnt)
{
    int   i, j, ddiff, nn = *n;
    double xmin, xmax, rang, dd;

    for (i = 0; i < *nb; i++) cnt[i] = 0;
    xmin = xmax = x[0];
    for (i = 1; i < nn; i++) {
	xmin = min(xmin, x[i]);
	xmax = max(xmax, x[i]);
    }
    rang = (xmax - xmin) * 1.01;
    *d = dd = rang / (*nb);
    for (i = 1; i < nn; i++) {
	for (j = 0; j < i; j++) {
          ddiff = abso(x[i]-x[j]) / dd;
	    cnt[ddiff]++;
	}
    }
}