Nothing
# $Id: slack.R 241 2022-03-16 14:00:23Z X052717 $
# Calculate slack at the efficient points.
# Hvis data ikke er transponeret bliver x og y transponeret og alle
# beregninger laves transponeres. Det betyder at resultat skal
# transponeres inden afslut og retur.
slack <- function(X, Y, e, XREF=NULL, YREF=NULL, FRONT.IDX=NULL,
LP=FALSE, CONTROL=NULL)
{
if ( !methods::is(e,"Farrell") ) {
stop("In call of slack: argument 'e' must be of class 'Farrell'")
}
RTS <- e$RTS
if ( RTS == "fdh0" ) RTS <- "fdh"
if (LP) print(paste("slack: RTS =",RTS),quote=F)
# Hvis data er en data.frame saa tjek om det er numerisk data og lav
# dem i saa fald om til en matrix
X <- tjek_data(X)
Y <- tjek_data(Y)
if ( missing(XREF) || is.null(XREF) ) XREF <- X
if ( missing(YREF) || is.null(YREF) ) YREF <- Y
XREF <- tjek_data(XREF)
YREF <- tjek_data(YREF)
if ( e$TRANSPOSE ) {
if (LP) print("X and Y are transposed")
X <- t(X)
Y <- t(Y)
XREF <- t(XREF)
YREF <- t(YREF)
}
okr <- dim(XREF)[1]
if (LP) cat("okr =",okr,"\n")
if ( FALSE && length(FRONT.IDX) == 0 ) {
# Kun units med eff==1 bruges i referenceteknologi
FRONT.IDX <- which( abs(e$eff -1) < 1e-5 )
}
if ( length(FRONT.IDX) > 0 ) {
if ( !is.vector(FRONT.IDX) )
stop("FRONT.IDX is not a vector in method 'slack'")
XREF <- XREF[FRONT.IDX,, drop=FALSE]
YREF <- YREF[FRONT.IDX,, drop=FALSE]
if (LP) cat("FRONT.IDX =",FRONT.IDX,"\n")
}
K = dim(X)[1] # number of units, firms, DMUs
Kr = dim(XREF)[1] # number of units in reference set
m = dim(X)[2] # number of inputs
n = dim(Y)[2] # number of outputs
if (LP) cat("In slack: m n K Kr = ",m,n,K,Kr,"\n")
if ( m != dim(XREF)[2] )
stop("Number of inputs must be the same in X and XREF")
if ( n != dim(YREF)[2] )
stop("Number of outputs must be the same in Y and YREF")
if ( K != dim(Y)[1] )
stop("Number of units must be the same in X and Y")
if ( Kr != dim(YREF)[1] )
stop("Number of units must be the same in XREF and YREF")
# For at undgaa afrundingsproblemer i forbindelse med slacks skal
# in- og ouput skaleres til at vaere i omegnen af 1
mmm <- (colMeans(X))
nnn <- (colMeans(Y))
if ( min(mmm) < 1e-4 || max(mmm) > 1e4 ||
min(nnn) < 1e-4 || max(nnn) > 1e4 ) {
SKALERING <- TRUE
X <- X / matrix(mmm, nrow=K, ncol=m, byrow=TRUE)
XREF <- XREF / matrix(mmm, nrow=Kr, ncol=m, byrow=TRUE)
Y <- Y / matrix(nnn, nrow=K, ncol=n, byrow=TRUE)
YREF <- YREF / matrix(nnn, nrow=Kr, ncol=n, byrow=TRUE)
} else {
SKALERING <- FALSE
}
if ( e$ORIENTATION=="graph" ) {
## stop(paste0("The use of the method 'slack' does not work for orientation 'graph'\n",
## " Use option SLACK in method 'dea'"))
lps <- make.lp(1, 1)
lpcontr <- lp.control(lps)
tol <- (lpcontr$epsilon["epsint"])
# delete.lp(lps)
rm(lps, lpcontr)
objval <- eff <- e$eff
sx <- X*eff - e$lambda %*% XREF
sy <- -Y/eff + e$lambda %*% YREF
sx[abs(sx) < tol] <- 0
sy[abs(sy) < tol] <- 0
sum <- rowSums(sx) + rowSums(sy)
sum[abs(sum) < tol] <- 0
colnames(sx) <- paste("sx",1:m,sep="")
colnames(sy) <- paste("sy",1:n,sep="")
if ( FALSE && e$TRANSPOSE ) {
sx <- t(sx)
sy <- t(sy)
}
if ( SKALERING ) {
sx <- sx * matrix(mmm, nrow=K, ncol=m, byrow=TRUE)
sy <- sy * matrix(nnn, nrow=K, ncol=n, byrow=TRUE)
}
oe <- list(eff=eff, slack=sum>tol, sum=sum, objval=objval, sx=sx, sy=sy,
lambda=lambda, RTS=e$RTS, ORIENTATION=e$ORIENTATION,
TRANSPOSE=e$TRANSPOSE)
class(oe) <- "slack"
return(oe)
} ## if (e$ORIENTATION=="graph")
if ( RTS != "crs" && RTS != "add" ) {
rlamb <- 2
} else {
rlamb <- 0
}
if ( is.null(e$objval) )
eff <- e$eff
else
eff <- e$objval
if ( is.null(e$direct) ) {
if ( e$ORIENTATION == "in" ) {
E <- eff
FF <- rep(1,K) # Naturligt at bruge F, men F er FALSE i R
} else if ( e$ORIENTATION == "out" ) {
FF <- eff
E <- rep(1,K)
} else if ( e$ORIENTATION == "graph" ) {
# Nedenstaaende virker ikke altid, derfor bruges beregningen
# under "if ( e$ORIENTATION=="graph" )"
stop("We should never end here")
lps <- make.lp(1, 1)
lpcontr <- lp.control(lps)
tol <- lpcontr$epsilon["epsint"]
# delete.lp(lps)
# Afrunding af eff. kan goere, at det beregnede punkt paa frontier
# faktisk ligger uden for frontier. Derfor adderes 'tol' til eff.
# tallet for at sikre den bliver inde i teknologimaengden.
# tol <- 0
E <- eff + tol
FF <- 1/eff - tol
rm(lps, lpcontr, tol)
} else {
stop("Unknown orientation in slack: ", e$ORIENTATION)
}
} else {
mmd <- switch(e$ORIENTATION, "in"=m, "out"=n, "in-out"=m+n)
ob <- matrix(e$objval,nrow=K, ncol=mmd)
if ( is(e$direct, "matrix") && dim(e$direct)[1] > 1 ) {
dir <- e$direct
} else {
dir <- matrix(e$direct,nrow=K, ncol=mmd, byrow=TRUE)
}
if ( e$ORIENTATION=="in" ) {
dirRhs <- cbind(X - ob*dir, -Y)
} else if ( e$ORIENTATION=="out" ) {
dirRhs <- cbind(X, -Y - ob*dir)
} else if ( e$ORIENTATION=="in-out" ) {
dirRhs <- cbind(X -ob[,1:m,drop=FALSE]*dir[,1:m,drop=FALSE],
-Y -ob[,(m+1):(m+n),drop=FALSE]*dir[,(m+1):(m+n),drop=FALSE])
} else {
warning("Illegal ORIENTATION for argument DIRECT in slacks",
immediate. = TRUE)
}
} ## if ( is.null(e$direct) )
# Initialiser LP objekt
lps <- make.lp(m+n +rlamb, m+n+Kr)
lp.control(lps,
scaling=c("range", "equilibrate", "integers") # default scalering er 'geometric'
) # og den giver ikke altid tilfredsstillende resultat;
# curtisreid virker i mange tilfaelde slet ikke
name.lp(lps, paste("DEA-slack",RTS,",",e$ORIENTATION,"orientated",sep="-"))
# saet raekker i matrix med restriktioner
for ( h in 1:m )
set.row(lps,h, XREF[,h], (m+n+1):(m+n+Kr) )
for ( h in 1:n)
set.row(lps,m+h, -YREF[,h], (m+n+1):(m+n+Kr) )
for ( h in 1:(m+n) )
set.mat(lps,h,h,1)
# restriktioner paa lambda
if ( RTS != "crs" && RTS != "add" ) {
set.row(lps, m+n+1, c(rep(0,m+n),rep(-1,Kr)))
set.row(lps, m+n+2, c(rep(0,m+n),rep( 1,Kr)))
}
set.constr.type(lps, c(rep("=",m+n), rep("<=",rlamb)))
set.objfn(lps, rep(1, m+n), 1:(m+n))
lp.control(lps, sense="max")
if ( RTS == "fdh" ) {
set.type(lps,(m+n+1):(m+n+Kr),"binary")
delete.constraint(lps, m+n+1)
set.rhs(lps,1, m+n+1)
rlamb <- rlamb -1
} else if ( RTS == "vrs" ) {
set.rhs(lps, c(-1,1), (m+n+1):(m+n+2))
} else if ( RTS == "drs" ) {
set.rhs(lps, -1, m+n+1)
delete.constraint(lps, m+n+2)
rlamb <- rlamb -1
set.constr.type(lps, ">=", m+n+1)
} else if ( RTS == "irs" ) {
set.rhs(lps, 1, m+n+2)
delete.constraint(lps, m+n+1)
rlamb <- rlamb -1
set.constr.type(lps, ">=", m+n+1)
} else if ( RTS == "irs2" ) {
set.rhs(lps, 1, m+n+2)
delete.constraint(lps, m+n+1)
rlamb <- rlamb -1
set.constr.type(lps, ">=", m+n+1)
set.semicont(lps, 2:(1+Kr))
set.bounds(lps, lower=rep(1,Kr), columns=(m+n+1):(m+n+Kr))
} else if ( RTS == "add" ) {
set.type(lps, (m+n+1):(m+n+Kr),"integer")
} else if ( RTS == "fdh+" ) {
param <- e$param
low <- param["low"]
high <- param["high"]
set.rhs(lps, c(-low,high), (m+n+1):(m+n+2))
add.SOS(lps,"lambda", 1,1, (m+n+1):(m+n+Kr), rep(1, Kr))
}
if (LP) {
print("Right hand side not sat yet")
print(lps)
}
objval <- rep(NA,K) # vector with sums of slacks
sx <- matrix(NA,K,m)
sy <- matrix(NA,K,n)
lambda <- matrix(NA,K,Kr)
if (!missing(CONTROL)) set_control(lps, CONTROL)
for ( k in 1:K ) { # beregn for hver enhed
if (LP) cat("\n---\nFirm",k,"\n")
# her er der problem med afrunding, sammen med venstrsiden er
# det ikke sikkert at restriktionen holder naar der er
# sket afrunding i mellemregningerner.
if ( is.null(e$direct) ) {
rhs <- c(E[k] * X[k,], -FF[k] * Y[k,])
} else {
rhs <- dirRhs[k,]
}
set.rhs(lps, rhs, 1:(m+n))
if (LP) {
print(paste("Effektivitet er ", E[k]))
print(paste("Hoejresiden for firm",k))
if ( is.null(e$direct) ) {
print(E[k] * X[k,])
print( -FF[k] * Y[k,])
} else {
print(rhs)
}
print(lps)
}
if ( !is.null(LP) && k %in% LP ) {
write.lp(lps, paste(name.lp(lps),k,".mps",sep=""),
type="mps",use.names=TRUE)
}
set.basis(lps, default=TRUE)
status <- solve(lps)
if (LP) print(paste0("status = ", status))
if ( status != 0 ) {
if ( status == 2 || status == 3 ) {
if (LP) print(paste("Firm",k,"not in the technology set"), quote=F)
if (LP) print(paste("Status =",status))
objval[k] <- 0
sol <- rep(0,m+n+Kr)
sol[m+n+k] <- 1
} else {
print(paste("Error in solving for firm",k,": Status =",status),
quote=F)
objval[k] <- NA
sol <- NA
}
} else {
objval[k] <- get.objective(lps)
sol <- get.variables(lps)
}
sx[k,] <- sol[1:m]
sy[k,] <- sol[(m+1):(m+n)]
lambda[k,] <- sol[(m+n+1):(m+n+Kr)]
if (LP) {
print(paste("Obj.value =",get.objective(lps)))
cat(paste("Solutions for",k,": "))
print(sol)
}
} # loop for each firm; "for (k in 1:K)"
# Drop soejler i lambda der alene er nuller, dvs lambda soejler skal
# alene vaere reference firms
colnames(sx) <- paste("sx",1:m,sep="")
colnames(sy) <- paste("sy",1:n,sep="")
# For at undgaa at regneunoejagtigheder giver ikke-nul slack bliver
# slack sat til nul hvis de er mindre end den relative noejagitghed
# i input og output, noejagtighed der er i beregningerne. Det
# betyder at sum af slacks bliver sum mens objval er sum plus
# regneunoejagtighed som giver en forskel hvis X er meget stoerre
# eller mindre end 1.
lpcontr <- lp.control(lps)
eps <- lpcontr$epsilon["epsint"]
sx[abs(sx) < eps ] <- 0
sy[abs(sy) < eps ] <- 0
if ( SKALERING ) {
sx <- sx * matrix(mmm, nrow=K, ncol=m, byrow=TRUE)
sy <- sy * matrix(nnn, nrow=K, ncol=n, byrow=TRUE)
}
sum <- rowSums(sx) + rowSums(sy)
lambda[abs(lambda-1) < eps] <- 1 # taet ved 1
lambda[abs(lambda) < eps] <- 0 # taet ved 0
if ( length(FRONT.IDX)>0 ) {
colnames(lambda) <- paste("L",(1:okr)[FRONT.IDX],sep="")
} else {
colnames(lambda) <- paste("L",1:Kr,sep="")
}
if (FALSE && LP) {
print("Done with loop\nColumn names for lambda:")
print(colnames(lambda))
print(lambda)
}
if ( e$TRANSPOSE ) {
sx <- t(sx)
sy <- t(sy)
lambda <- t(lambda)
}
if ( FALSE && LP ) {
print("Faerdig med slack: lambda")
print(colnames(lambda))
print(lambda)
print("slack:")
print(paste("laengden af objval:",length(objval)))
print(objval==0)
# print(!slackZero)
}
oe <- list(eff=eff, slack=sum>eps, sum=sum, objval=objval, sx=sx, sy=sy,
lambda=lambda, RTS=e$RTS, ORIENTATION=e$ORIENTATION,
TRANSPOSE=e$TRANSPOSE)
class(oe) <- "slack"
return(oe)
} ## slack()
print.slack <- function(x, digits=4, ...) {
sum <- x$sum
print(sum, digits=digits, ...)
invisible(sum)
} ## print.slack
summary.slack <- function(object, digits=4, ...) {
eps <- 1e-7
cat("Efficiency and slacks:\n")
# if ( object$ORIENTATION != "out" )
# a <- cbind(E=round(object$eff, digits=2),
# "sum"=object$sum, object$sx, object$sy)
# else
# a <- cbind(F=round(object$eff,2),
# "sum"=object$sum, object$sx, object$sy)
# print(a,...)
cat("Number of firms with efficiency==1 and positive slacks:",
sum(abs(object$eff-1) < eps & object$slack , na.rm=TRUE),"\n" )
if ( dim(object$sx)[2]==1 ) {
SX <- object$sx > eps
} else {
SX <- rowSums(object$sx, na.rm=TRUE) > eps
}
if ( dim(object$sy)[2]==1 ) {
SY <- object$sy > eps
} else {
SY <- rowSums(object$sy, na.rm=TRUE) > eps
}
cat("Number of firms with:\n")
cat(" only x slacks: ", sum(SX & !SY, na.rm=TRUE), "\n")
cat(" only y slacks: ", sum(!SX & SY, na.rm=TRUE), "\n")
cat(" x and y slacks:", sum(SY & SX, na.rm=TRUE), "\n")
cat("\n x slacks: ", sum(SX, na.rm=TRUE), "\n")
cat( " y slacks: ", sum(SY, na.rm=TRUE), "\n")
cat( "all slacks: ", sum(SX | SY, na.rm=TRUE), "\n")
if (sum(is.na(SX) | is.na(SY))) cat("Number of firms Not Available for slacks: ",
sum(is.na(SX)|is.na(SY)),"\n")
# invisible(a)
if ( FALSE && dim(lambda(object))[2] <= 10 ) {
cat("Weights (lambda):\n")
x <- lambda(object)
xx <- round(x, digits)
# xx <- format(unclass(x), digits=digits)
if (any(ina <- is.na(x)))
xx[ina] <- ""
if ( any(i0 <- !ina & abs(x) < 1e-9) )
xx[i0] <- sub("0", ".", xx[i0])
print(xx, quote=FALSE, rigth=TRUE, ...)
invisible(xx)
}
} ## summary.slack()
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.