#  This program is free software; you can redistribute it and/or modify
#  it under the terms of the GNU General Public License as published by
#  the Free Software Foundation; either version 2 of the License, or
#  (at your option) any later version.
#
#  This program is distributed in the hope that it will be useful,
#  but WITHOUT ANY WARRANTY; without even the implied warranty of
#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#  GNU General Public License for more details.
#
#  A copy of the GNU General Public License is available at
#  http://www.r-project.org/Licenses/

### copyright (C) 2011 A. Pepelyshev
### 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.


bw.icv <- function(x, nb = 1000, lower = 0.05*h_os, upper = h_os)
{
  
    if((n <- length(x)) < 2) stop("need at least 2 data points")
    if(!is.numeric(x)) stop("invalid 'x'")
    if(n<=100)
    {
    a <- 25.2 
    s <- 1.39
    }
    else
    {
     z <- log10(n)
     a <- 10^(3.390-1.093*z+0.025*z^3-0.00004*z^6)
     s <- 10^(-0.58+0.386*z-0.012*z^2)
    }
    h_os <- 3*(35*n*2*sqrt(pi))^(-1/5)*sd(x)
    Rphi <- 1/sqrt(2) 
    RL  <- (1+a)^2/sqrt(2) - 2*a*(1+a)/sqrt(1+s*s) + a*a/(s*sqrt(2)) 
    muL <- 1+a-a*s^2
    C <- (Rphi/RL*muL^2)^(1/5)

    lower_hL <- lower/C
    upper_hL <- upper/C

    ficv <- function(h, x, n, d, a, s)
        .C("band_icv_bin",
           as.integer(n),
           as.integer(length(x)),
           as.double(d),
           x,
           as.double(h),
           u = double(1),
           as.double(a),
           as.double(s))$u

    storage.mode(x) <- "double"
    Z <- .C("band_ddiff_bin",
            as.integer(n),
            as.integer(nb),
            d = double(1),
            x,
            cnt = integer(nb))
    d <- Z$d; cnt <- as.integer(Z$cnt)

    Z <- .C("band_icv_first_local_min",
           as.integer(n),
           as.integer(length(cnt)),
           as.double(d),
           cnt,
           h1=as.double(lower_hL),
           h2=as.double(upper_hL),
           as.double(a),
           as.double(s))
    lower_hL <- Z$h1
    upper_hL <- Z$h2

    hL <- optimize(ficv, c(lower_hL, upper_hL), tol=0.05*(upper_hL-lower_hL),
                  x=cnt, n=n, d=d, a=a, s=s)$minimum
    h <- hL*C
    if(h>h_os)
    {
     h <- h_os;
    }
   if(h < 1.1*lower | h > upper-0.1*lower)
        warning("minimum occurred at one end of the range")
    h
}

#################################
#======================
bw.icv.slow=function(x)
{
n <- length(x)

    if(n<=100)
    {
    a <- 25.2 
    s <- 1.39
    }
    else
    {
     z <- log10(n)
     a <- 10^(3.390-1.093*z+0.025*z^3-0.00004*z^6)
     s <- 10^(-0.58+0.386*z-0.012*z^2)
    }

    h_os <- 3*(35*n*2*sqrt(pi))^(-1/5)*sd(x)
    Rphi <- 1/sqrt(2) 
    RL  <- (1+a)^2/sqrt(2) - 2*a*(1+a)/sqrt(1+s*s) + a*a/(s*sqrt(2)) 
    muL <- 1+a-a*s^2
    C <- (Rphi/RL*muL^2)^(1/5)

    upper_hL <- h_os/C

hval=0.05*upper_hL+c(1:100)*0.95*upper_hL/100; y=hval;
h1=hval[1];
h2=upper_hL;
i=0; Do=1;
while( i<length(hval) & Do )
{
 i=i+1;
 y[i]=fun_qk(hval[i],x,a,s);
 if(i>=3)
 {
  if(y[i-1]<y[i]&y[i-1]<y[i-2])
  {
     h1=hval[i-2];
     h2=hval[i];
     Do=0;
  }
 }
}

res <- optimize(band_icv_crit, interval=c(h1,h2), tol=0.05*(h2-h1), x, a, s)
hL <- res$minimum

 h <- C*hL
 if(h>h_os)
 {
    h <- h_os;
 }
 h
}

#===========================================
band_icv_crit=function(h, x, a, s)
{

n=length(x);
s2=sqrt(2);

D=matrix(0,n,n);
for( i in (1:n))
{
 D[,i]=(x-x[i])/(h);
}
M2h=kernLL( D ,a,s)/(h);
Mh=kernL( D ,a,s)/(h);

for (i in (1:n) )
{
    Mh[i,i]=0;
}
y=(sum(M2h))/n^2-2/(n*(n-1))*(sum(Mh));
return(y);
}

#============================
kernL=function(x,a,s)
{
y=(1+a)*exp(-x^2*0.5)-a/s*exp(-0.5*(x/s)^2);
return(y);
}

#===========================================
kernLL=function(x,a,s)
{
s2=sqrt(2);
ss=sqrt(1+s^2);
y=(1+a)^2/s2*exp(-x^2*0.5/2)- 2*a*(1+a)/ss*exp(-0.5*x^2/(1+s^2)) +a^2/s/s2*exp(-0.5*x^2/(2*s^2));
return(y);
}