#  ===========================
#  order-m efficiency measures
#  ===========================
#  written by L. Simar, July 8, 2001 in MATLAB
#  translated by O. Badunenko June 25, 2007 at 15:23
# 
#  IN
#  ------
#     X         : Matrix of input(s) (p x n)
#     Y         : Matrix of output(s) (q x n)
#     x0, y0    : coordinate of the reference point
#     m         : Order of the frontier
#     B         : Number of Monte-Carlo replications
# 
#  OUT
#  ------
#     B          : used B
#     m          : used m
#     input      : input-oriented efficiency score of ORDER-m, std, number of points in x, 
#                  with output level greater or equal to y0
#     output     : output-oriented efficiency score of ORDER-m, std, number of points in y, 
#                  with output smallerr or equal to x0
# 
# 
# 
eff_m = function(X, Y, x0, y0, m = 25, B = 200){
# 
# 
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(i in 1:q){
    y.tmp[i,] = ifelse(Y[i,] >= y0[i], 1, 0)
}
y.tmp1 = apply(y.tmp, MARGIN = 2, FUN = "mean")
# 
XM = X[,y.tmp1 == 1]; XM = as.matrix(XM); if(ncol(XM) == 1){ XM = t(XM) }
# 
# --------------------------------
# OUTPUT orientation
# 
x.tmp = X
for(i in 1:p){
    x.tmp[i,] = ifelse(X[i,] <= x0[i], 1, 0)
}
x.tmp1 = apply(x.tmp, MARGIN = 2, FUN = "mean")
# 
YM = Y[,x.tmp1 == 1]; YM = as.matrix(YM); if(ncol(YM) == 1){ YM = t(YM) }
# 
#   start Monte-Carlo loop
# 
thetab = matrix(nrow = 2, ncol = B);
YMb = matrix(nrow = q,ncol = m)
XMb = matrix(nrow = p,ncol = m)
# 
for(i in 1:B){
#    for(qq in 1:m){
#        a1 = sample(c(1:ncol(XM)), 1, replace = T)
#        a2 = sample(c(1:ncol(YM)), 1, replace = T)
#        XMb[,qq] = XM[,a1]
#        YMb[,qq] = YM[,a2]
#    }
    if(ncol(YM) != 0){
        sample_y = floor(ncol(YM) * runif(m) + 1);  YMb = YM[,sample_y];
        y0v = y0 %*% matrix(1, nrow = 1, ncol = m);
        ratioy0 = YMb / y0v;
        thetab[2,i] = max(apply(ratioy0, MARGIN = 2, FUN = "min"));
    }
    if(ncol(XM) != 0){
        sample_x = floor(ncol(XM) * runif(m) + 1);  XMb = XM[,sample_x];
        x0v = x0 %*% matrix(1, nrow = 1, ncol = m);
        ratiox0 = XMb / x0v;
        thetab[1,i] = min(apply(ratiox0, MARGIN = 2, FUN = "max"));
    }
}
# 
# end of Monte-Carlo loop
# 
if(ncol(XM) == 0){
    cat("\n\n")
    cat("   order-m input frontier not available, y0 is too large","\n\n\n");
    I_theta = Inf;
    I_theta_sd = 0;
    } else {
        I_theta = mean(thetab[1,]);
        I_theta_sd = sd(thetab[1,]) / sqrt(B);
      }
if(ncol(YM) == 0){
    cat("\n\n")
    cat("   order-m output frontier not available, x0 is too small","\n\n\n");
    O_theta = 0;
    O_theta_sd = 0;
    } else {
        O_theta = mean(thetab[2,]);
        O_theta_sd = sd(thetab[2,]) / sqrt(B);
      }
input  = c(I_theta, I_theta_sd, ncol(XM));
output = c(O_theta, O_theta_sd, ncol(YM));
# 
# 
rez = list(B = B, m = m, input = input, output = output);
rez
return(rez)
# end of the function eff_m
}
#*********************************************************************
# y = t(matrix(c(ccr$y1,ccr$y2,ccr$y3),nrow = 70,ncol = 3))
# x = t(matrix(c(ccr$x1,ccr$x2,ccr$x3,ccr$x4,ccr$x5),nrow = 70,ncol = 5))
# i = 48
# y0 = as.matrix(y[,i]); x0 = as.matrix(x[,i]); X = x[,-i]; Y = y[,-i];
# 
# eff_m(X = x[,-i], Y = y[,-i], x0 = as.matrix(x[,i]), y0 = as.matrix(y[,i]), m = 10, B = 200)