### finding the bandwidth

kern_e = function(x){ 0.75 * (1 - x^2) * ifelse(abs(x) <= 1, 1, 0) }        # Epanechnikov
kern_u = function(x){ 0.5 * ifelse(abs(x) <= 1, 1, 0) }                     # Uniform
kern_t = function(x){ (1 - abs(x)) * ifelse(abs(x) <= 1, 1, 0) }            # Triangle
kern_q = function(x){ 15 / 16 * (1 - x^2)^2 * ifelse(abs(x) <= 1, 1, 0) }   # Quartic (Biweight)
# 
# 
# 
hzi = function(z, zi, k){
    d = sort( abs( z - zi ) )
    k = as.integer(k)
    return(d[k])
}
# 
# 
# 
fkizi = function(z, zi, hzi){
    a = 1 / (n - 1) / hzi
    b = sum( kern_u( (z - zi) / hzi ) ) - kern_u(0)
    return(a * b)
}
# 
# 
# 
cv = function(z, k){
    cv = rep(0,length(z))
    for(i in 1:length(z)){
        cv[i] = log( fkizi(z, z[i], hzi(z, z[i], k) ) )
    }
    return(mean(cv))
}
# 
# 
# 
h = function(z){
tmp = optimize(cv, c(1, length(z)), maximum = T, z = z)
k = floor(tmp$maximum)
k
q = rep(0,length(z))
for(i in 1:length(z)){
    q[i] = hzi(z, z[i], k)
}
tmp = list(h = mean(q), k = k)
return(tmp)
}

# ============================
# tricky distribution function
# 
HXYn = function(X, x, Y, y){
    Y = as.matrix(Y); if(ncol(Y) == 1){ Y = t(Y) }; y = as.matrix(y);
    y.tmp = Y
    for(qq in 1:nrow(Y)){
        y.tmp[qq,] = ifelse(Y[qq,] >= y[qq], 1, 0)
    }
    y.tmp1 = apply(y.tmp, MARGIN = 2, FUN = "mean")
    y.tmp2 = ifelse(y.tmp1 == 1, 1, 0)
    # 
    X = as.matrix(X);if(ncol(X) == 1){ X = t(X) }; x = as.matrix(x);
    x.tmp = X
    for(qq in 1:nrow(X)){
        x.tmp[qq,] = ifelse(X[qq,] <= x[qq], 1, 0)
    }
    x.tmp1 = apply(x.tmp, MARGIN = 2, FUN = "mean")
    x.tmp2 = ifelse(x.tmp1 == 1, 1, 0)
    # 
    xy.tmp = sum(y.tmp2 * x.tmp2) / ncol(Y)
    return(xy.tmp)
}
# 
# ======================
# finding Y or X
# 
YX = function(Y, y, minmax = min){
    Y = as.matrix(Y); if(ncol(Y) == 1){ Y = t(Y) }
    y = as.matrix(y); if(ncol(y) != 1){ y = t(y) }
    n = ncol(Y);
    yv = y %*% matrix(1, nrow = 1, ncol = n);
    ratio = Y / yv;
    apply(ratio, MARGIN = 2, FUN = "minmax");
}
# 
# 
# ===============================
# a-quantile efficiency measures 
# 
#  written by O. Badunenko June at July 4, 2007
# 
eff_a = function(X, Y, x0, y0, alpha = 0.95){
    # 
    X = as.matrix(X); if(ncol(X) == 1){ X = t(X) }
    Y = as.matrix(Y); if(ncol(Y) == 1){ Y = t(Y) }
    n = ncol(X);
    p = nrow(X);
    q = nrow(Y);
    # 
    # INPUT orientation
    # 
    y.tmp = Y
    for(j in 1:q){
        y.tmp[j,] = ifelse(Y[j,] >= y0[j], 1, 0)
    }
    y.tmp1 = apply(y.tmp, MARGIN = 2, FUN = "mean")
    # 
    X.a = YX(Y = X, y = x0, minmax = max);
    X.a = X.a[y.tmp1 == 1];
    # 
    aMy = (1 - alpha) * n * HXYn(X = X, x = rep(Inf,p), Y = Y, y = y0)
    if( (aMy - floor(aMy) ) > 0 ){aMy = floor(aMy) + 1}
    theta = sort(X.a)[aMy] # this is the score
    # 
    # OUTPUT orientation
    # 
    x.tmp = X
    for(j in 1:p){
        x.tmp[j,] = ifelse(X[j,] <= x0[j], 1, 0)
    }
    x.tmp1 = apply(x.tmp, MARGIN = 2, FUN = "mean")
    # 
    Y.a = YX(Y = Y, y = y0, minmax = min);
    Y.a = Y.a[x.tmp1 == 1];
    aNx = alpha * n * HXYn(X = X, x = x0, Y = Y, y = rep(0,q))
    if( (aNx - floor(aNx) ) > 0 ){aNx = floor(aNx) + 1}
    lambda = sort(Y.a)[aNx] # this is the score
    result = list(output = lambda, input = theta)
    return(result)
}

# eff_a(Y = y, X = x, x0 = x[,1], y0 = y[,1])
# 
# Symmetric kernels with compact support
# 
kern_e = function(x){ 0.75 * (1 - x^2) * ifelse(abs(x) <= 1, 1, 0) }        # Epanechnikov
kern_u = function(x){ 0.5 * ifelse(abs(x) <= 1, 1, 0) }                     # Uniform
kern_t = function(x){ (1 - abs(x)) * ifelse(abs(x) <= 1, 1, 0) }            # Triangle
kern_q = function(x){ 15 / 16 * (1 - x^2)^2 * ifelse(abs(x) <= 1, 1, 0) }   # Quartic (Biweight)
# 
# 

# 
# calculating Rxz
#

Rxz = function(X, x0, Z, z0){
    X = as.matrix(X); if(ncol(X) == 1){ X = t(X) }; x0 = as.matrix(x0);
    x.tmp = X
    for(qq in 1:nrow(X)){
        x.tmp[qq,] = ifelse(X[qq,] <= x0[qq], 1, 0)
    }
    x.tmp1 = apply(x.tmp, MARGIN = 2, FUN = "mean")
    a = ifelse(x.tmp1 == 1, 1, 0)
    # 
    hz = h(Z)$h
    b = kern_e( (z0 - Z) / hz )
    c = sum( a * b )
    return(c)
}

# 
# 
#i = 2
#Rxz(X = x, x0 = x[,i], Z = z, z0 = z[,i])
# 
# 
# calculating Qyz
#
Qyz = function(Y, y0, Z, z0){
    Y = as.matrix(Y); if(ncol(Y) == 1){ Y = t(Y) }; y0 = as.matrix(y0);
    y.tmp = Y
    for(qq in 1:nrow(Y)){
        y.tmp[qq,] = ifelse(Y[qq,] >= y0[qq], 1, 0)
    }
    y.tmp1 = apply(y.tmp, MARGIN = 2, FUN = "mean")
    a = ifelse(y.tmp1 == 1, 1, 0)
    # 
    b = kern_e( (z0 - Z) / hz )
    c = sum( a * b )
    return(c)
}
# 
# 
#i = 2
#Qyz(Y = y, y0 = y[,i], Z = z, z0 = z[,i])
# 
# ==============
# calculating lk
# 
# do not forget: it has to be divided by Qyz
# 
lk = function(z0, Z.y, k){
    sum( kern_e( (z0 - Z.y[1:k]) / hz ) )
}
lk_ = function(z0, Z.x, k){
    sum( kern_e( (z0 - Z.x[k:length(Z.x)]) / hz ) )
}
# 
# ===========================================================
# a-quantile efficiency measures given environmental variable
# 
#  written by O. Badunenko June at July 4, 2007
#

eff_a_z = function(X, Y, Z, x0, y0, z0, alpha = 0.95){
    X = as.matrix(X); if(ncol(X) == 1){ X = t(X) }
    Y = as.matrix(Y); if(ncol(Y) == 1){ Y = t(Y) }
    n = ncol(X);
    p = nrow(X);
    q = nrow(Y);
    # 
    # INPUT orientation
    # 
    y.tmp = Y
    #    y0 = y[,1]
    #    x0 = x[,1]
    #    z0 = z[,1]
    for(j in 1:q){
        y.tmp[j,] = ifelse(Y[j,] >= y0[j], 1, 0)
    }
    y.tmp1 = apply(y.tmp, MARGIN = 2, FUN = "mean")
    # 
    X.a = YX(Y = X, y = x0, minmax = max);
    X.a = X.a[y.tmp1 == 1];
    Z.y = Z[y.tmp1 == 1];
    X.a.Z.y = rbind(X.a,Z.y)[, order(X.a,Z.y)]
    M.y = n * HXYn(X = X, x = rep(Inf,p), Y = Y, y = y0)
    lk.s = rep(0,M.y)
    if(length(X.a.Z.y) == 2){ X.a.Z.y = as.matrix(X.a.Z.y) }
    for(tt in 1:M.y){
        lk.s[tt] = lk(z0 = z0, Z.y = X.a.Z.y[2,], k = tt)
    }
    lk.s = c(0,lk.s) / Qyz(Y = Y, y0 = y0, Z = Z, z0 = z0)
    # 
    beta = 1 - alpha
    a = rep(0,M.y)
    for(rr in 1:M.y){
        a[rr] = ifelse(beta >= lk.s[rr] & beta < lk.s[(rr + 1)], 1, 0)
    }
    theta = X.a.Z.y[1,a == 1] # this is the score
    # 
    # OUTPUT orientation
    # 
    x.tmp = X
    #    y0 = y[,1]
    #    x0 = x[,1]
    #    z0 = z[,1]
    for(j in 1:p){
        x.tmp[j,] = ifelse(X[j,] <= x0[j], 1, 0)
    }
    x.tmp1 = apply(x.tmp, MARGIN = 2, FUN = "mean")
    # 
    Y.a = YX(Y = Y, y = y0, minmax = min);
    Y.a = Y.a[x.tmp1 == 1];
    Z.x = Z[x.tmp1 == 1]
    Y.a.Z.x = rbind(Y.a,Z.x)[, order(Y.a,Z.x)]
    N.x = n * HXYn(X = X, x = x0, Y = Y, y = rep(0,q))
    lk.s = rep(0,N.x)
    if(length(Y.a.Z.x) == 2){ Y.a.Z.x = as.matrix(Y.a.Z.x) }
    for(tt in 1:N.x){
        lk.s[tt] = lk_(z0 = z0, Z.x = Y.a.Z.x[2,], k = tt)
    }
    lk.s = c(lk.s,0) / Rxz(X = X, x0 = x0, Z = Z, z0 = z0)
    # 
    beta = 1 - alpha
    a = rep(0,N.x)
    for(rr in 1:N.x){
        a[rr] = ifelse(beta >= lk.s[(rr + 1)] & beta < lk.s[rr], 1, 0)
    }
    lambda = Y.a.Z.x[1,a == 1] # this is the score
    result = list(output = lambda, input = theta)
    return(result)
}

#####################

# USAGE

#####################

library(FEAR)
data(ccr)
y = t(matrix(c(ccr$y2,ccr$y3),nrow = 70,ncol = 2))
x = t(matrix(c(ccr$x1,ccr$x2,ccr$x3,ccr$x4,ccr$x5),nrow = 70,ncol = 5))
z = t(matrix(c(ccr$y1),nrow = 70,ncol = 1))

input  = matrix(nrow = 2, ncol = 70);
output = matrix(nrow = 2, ncol = 70);

n = 70
hz = h(z)$h

for(i in 1:70){

y0 = as.matrix(y[,i]);
x0 = as.matrix(x[,i]);
z0 = z[,i];

tmp = eff_a_z(X = x, Y = y, Z = z, x0 = x0, y0 = y0, z0 = z0)
input[1,i]  = tmp$input
output[1,i] = tmp$output
tmp = eff_a(X = x, Y = y, x0 = x0, y0 = y0)
input[2,i]  = tmp$input
output[2,i] = tmp$output

}

# Plots

threshold = 60
y_m = max( (output[1,]/output[2,])[z < threshold] )
y_min = min( (output[1,]/output[2,])[z < threshold] )
x_min = min( c(z) )

plot(z[z < threshold], (output[1,]/output[2,])[z < threshold],
ylim = c((y_min / 1.01),(1.01*y_m)),
xlim = c((x_min / 1.01),(1.01*threshold )),
xlab = "Z",
ylab = expression( paste( lambda[alpha] (x, "y | Z") / lambda[alpha] (x,y) ) ), 
main = paste("Effect of 'Z' on Order-alpha frontier", sep = ""))
lines(ksmooth(z[z < threshold], (output[1,]/output[2,])[z < threshold], "normal", bandwidth = hz),
col=3, lwd = 5)