#  written by O. Badunenko. Summer of '07

library(FEAR)

folder = "c:/temp/" # indicate your's

extreme.s03 = function(x = x, y = y, B = 200, fig = "eps", png.dim = 714){
n = ncol(x); p = nrow(x); q = nrow(y); pq = p + q
# 
# preparing a table for results of output and input oriented 
# "leave-one-out order-m and FDH" efficiency measures
# called "rez_o" and "rez_i", respectively
# 
rez_o = matrix(,nrow = 3*n, ncol = 7);
rez_i = matrix(,nrow = 3*n, ncol = 7);
name = rep(0,length(m));
for(qq in 1:length(m)){     
    name[qq] = paste("m =", m[qq])
}
colnames(rez_o) = c(name, "FDH")
rownames(rez_o) = sort( c( c(1:n),c(1:n),c(1:n) ) ); # can be changed for real names or id's
colnames(rez_i) = c(name, "FDH")
rownames(rez_i) = sort( c( c(1:n),c(1:n),c(1:n) ) ); # can be changed for real names or id's
# 
# filling in the table
# achtung: these are reciprocals of scores, calculated by ``orderm'' procedure
# 
time1 = proc.time()
for(i in 1:n){
    for(j in 1:length(m)){
#        tmp = orderm(XOBS = as.matrix(x[,i]), YOBS = as.matrix(y[,i]), 
#                     XREF = x[,-i], YREF = y[,-i],
#                     ORIENTATION = 3, M = m[j], NREP = 200)
        # observations
        x0 = as.matrix(x[,i]);
        y0 = as.matrix(y[,i]);
        X = as.matrix(x[,-i]); if(ncol(X) == 1){ X = t(X) }
        Y = as.matrix(y[,-i]); if(ncol(Y) == 1){ Y = t(Y) }
        tmp = eff_m(X = X, Y = Y, x0 = x0, y0 = y0, m = m[j], B = B)
        a = (i - 1)*3 + 1; b = a + 1; c = b + 1; d = length(m) + 1;
        rez_o[a,j] = tmp$output[1];
        rez_i[a,j] = tmp$input[1];
        rez_o[a,d] = 1/fdh(XOBS = x0, YOBS = y0, XREF = X, YREF = Y, ORIENTATION = 2)[1,1]
        rez_i[a,d] = 1/fdh(XOBS = x0, YOBS = y0, XREF = X, YREF = Y, ORIENTATION = 1)[1,1]
        rez_o[b,j] = tmp$output[2];
        rez_i[b,j] = tmp$input[2];
        rez_o[c,j] = tmp$output[3];
        rez_i[c,j] = tmp$input[3];
    }
}
time2 = proc.time()
time = (time2-time1)/60
time1 =  as.vector(time[1])
cat("   Time elapsed is", time1, "minutes","\n\n\n")
# rez_o[1:8,]
write.table(rez_o, file = paste("",folder,"rez_o.csv", sep = ""), sep = ",", row.names = F, col.names = T)
write.table(rez_i, file = paste("",folder,"rez_i.csv", sep = ""), sep = ",", row.names = F, col.names = T)
# 
# 
# 
# 
#   start a semi-automatic warning procedure
# 

# 
# table to determine out of frontier points
# 
out_o = matrix(,nrow = n, ncol = length(m))
for(i in 1:n){
    for(j in 1:length(m)){
        a = (i - 1)*3 + 1; b = a + 1;
        out_o[i,j] = rez_o[a,j] + 1.645 * rez_o[b,j]
    }
}
# 
out_i = matrix(,nrow = n, ncol = length(m))
for(i in 1:n){
    for(j in 1:length(m)){
        a = (i - 1)*3 + 1; b = a + 1;
        out_i[i,j] = rez_i[a,j] - 1.645 * rez_i[b,j]
    }
}
# 
# function count those greater than "a"
greater = function(x, a){
    tmp = ifelse(x >= a, 1, 0)
    return(sum(tmp, na.rm = T))
}
# 
# function count those smaller than "a"
smaller = function(x, a){
    tmp = ifelse(x <= a, 1, 0)
    return(sum(tmp, na.rm = T))
}
# 
# table with number of those outside frontier
# 
outside_o = matrix(nrow = length(alpha), ncol = length(m));
outside_i = matrix(nrow = length(alpha), ncol = length(m));
name_ = rep(0,length(alpha));
for(qq in 1:length(alpha)){     
    name_[qq] = paste("alpha =", alpha[qq])
}
colnames(outside_o) = name;
rownames(outside_o) = name_;
colnames(outside_i) = name;
rownames(outside_i) = name_;
# 
for(i in 1:length(alpha)){
    for(j in 1:length(m)){
        outside_o[i,j] = smaller( out_o[,j], 1 - alpha[i] ) / length(na.omit(out_o[,j]))
        outside_i[i,j] = greater( out_i[,j], 1 + alpha[i] ) / length(na.omit(out_i[,j]))
    }
}
outside_o
outside_i


# 
# plot for input orientation
# 
if(fig == "eps"){
    postscript(paste("",folder,"input.eps", sep = ""),
    horizontal = F, height = 8, width = 8, paper="special") 
} else {
    png(filename = paste("",folder,"input.png", sep = ""), width = png.dim, height = png.dim,
    units = "px", pointsize = 12, bg = "white", res = NA,
    restoreConsole = TRUE)
}
y_max = round(max(outside_i*1.05),d = 1)
y_min = round(min(outside_i*0.95),d = 1)
plot(0,0,type="n",axes = F,xlim=c(0,160),ylim=c(y_min,y_max),
xlab="Value of m",
ylab="Proportion of sample pints outside m-frontier",
main="Input Orientation")
xx = seq(0, 160, by = 20)
axis(side = 1, at = xx)
yy = seq(y_min, y_max, by = 0.1)
axis(side = 2, at = yy, las = 1)
abline(v = xx, col = "grey", lty = 2)
abline(h = yy, col = "grey", lty = 2)
for (i in 1:length(alpha)){
    lines(m,outside_i[i,], lty = 1, lwd = 2, col = i)
    print(i)}
legend(117,y_max,paste("alpha = ", alpha), bg = "white", col=1:length(alpha), lty = 1, lwd = 2)
dev.off()
# 
# 
# plot for input orientation
# 
if(fig == "eps"){
    postscript(paste("",folder,"output.eps", sep = ""),
    horizontal = F, height = 8, width = 8, paper="special") 
} else {
    png(filename = paste("",folder,"output.png", sep = ""), width = png.dim, height = png.dim,
    units = "px", pointsize = 12, bg = "white", res = NA,
    restoreConsole = TRUE)
}
y_max = round(max(outside_o*1.05),d = 1)
y_min = round(min(outside_o*0.95),d = 1)
plot(0,0,type="n",axes = F,xlim=c(0,160),ylim=c(y_min,y_max),
xlab="Value of m",
ylab="Proportion of sample points outside m-frontier",
main="Ouput Orientation")
xx = seq(0, 160, by = 20)
axis(side = 1, at = xx)
yy = seq(y_min, y_max, by = 0.1)
axis(side = 2, at = yy, las = 1)
abline(v = xx, col = "grey", lty = 2)
abline(h = yy, col = "grey", lty = 2)
for (i in 1:length(alpha)){
    lines(m,outside_o[i,], lty = 1, lwd = 2, col = i)
    print(i)}
legend(117,y_max,paste("alpha = ", alpha), bg = "white", col=1:length(alpha), lty = 1, lwd = 2)
dev.off()

tmp = list(outside_o = outside_o,outside_i = outside_i, minutes = round(time1, d = 2))
return(tmp)
}

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

# HOW TO USE

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

data(ccr)
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))

# folder with all stuff of project
folder = "c:/temp/";

# 
# tuning parameter
# 
m = c (10, 25, 50, 75, 100, 150);
alpha = c(0.2, 0.3, 0.4, 0.5);

extreme.s03(x = x, y = y, B = 2)
extreme.s03(x = x, y = y, B = 2, fig = "png", png.dim = 680)