https://github.com/javierbarbero/DataEnvelopmentAnalysis.jl
Tip revision: 7caeaf2c5a4476b7ff367ae2647a78a4d507b327 authored by Javier Barbero on 28 August 2024, 06:52:42 UTC
Bump version to 0.9.1
Bump version to 0.9.1
Tip revision: 7caeaf2
deaholder.jl
# This file contains functions for the Hölder DEA model
"""
AbstractHolderDEAModel
An abstract type representing a Holder DEA model.
"""
abstract type AbstractHolderDEAModel <: AbstractTechnicalDEAModel end
"""
HolderL1DEAModel
An data structure representing a Höler L1 DEA model.
"""
struct HolderL1DEAModel <: AbstractHolderDEAModel
n::Int64
m::Int64
s::Int64
orient::Symbol
rts::Symbol
isweighted::Bool
dmunames::Union{Vector{AbstractString},Nothing}
eff::Vector
effmin::Vector
slackX::Matrix
slackY::Matrix
lambda::SparseMatrixCSC{Float64, Int64}
Xtarget::Matrix
Ytarget::Matrix
end
"""
HolderL2DEAModel
An data structure representing a Höler L2 DEA model.
"""
struct HolderL2DEAModel <: AbstractHolderDEAModel
n::Int64
m::Int64
s::Int64
orient::Symbol
rts::Symbol
isweighted::Bool
dmunames::Union{Vector{AbstractString},Nothing}
eff::Vector
slackX::Matrix
slackY::Matrix
lambda::SparseMatrixCSC{Float64, Int64}
Xtarget::Matrix
Ytarget::Matrix
end
"""
HolderLInfDEAModel
An data structure representing a Höler LInf DEA model.
"""
struct HolderLInfDEAModel <: AbstractHolderDEAModel
n::Int64
m::Int64
s::Int64
orient::Symbol
rts::Symbol
isweighted::Bool
dmunames::Union{Vector{AbstractString},Nothing}
eff::Vector
slackX::Matrix
slackY::Matrix
lambda::SparseMatrixCSC{Float64, Int64}
Xtarget::Matrix
Ytarget::Matrix
end
"""
deaholder(X, Y; l)
Compute the Hölder distance function model using data envelopment analysis for inputs `X` and outputs `Y`,
using Hölder norm `l`.
# Hölder norm `l` specification
- `1`.
- `2`.
- `Inf`.
# Optional Arguments
- `weigt=false`: set to `true` for weighted (weakly) Hölder distance function.
- `orient=:Input`: chooses the radially oriented input mode. For the radially oriented output model choose `:Output`.
- `rts=:CRS`: chooses constant returns to scale. For variable returns to scale choose `:VRS`.
- `slack=true`: computes input and output slacks.
- `Xref=X`: Identifies the reference set of inputs against which the units are evaluated.
- `Yref=Y`: Identifies the reference set of outputs against which the units are evaluated.
- `names`: a vector of strings with the names of the decision making units.
"""
function deaholder(X::Union{Matrix,Vector}, Y::Union{Matrix,Vector};
l::Union{Int64,Float64}, weight::Bool = false,
orient::Symbol = :Graph, rts::Symbol = :CRS, slack = true,
Xref::Union{Matrix,Vector,Nothing} = nothing, Yref::Union{Matrix, Vector,Nothing} = nothing,
names::Union{Vector{<: AbstractString},Nothing} = nothing,
optimizer::Union{DEAOptimizer,Nothing} = nothing)::AbstractHolderDEAModel
# Check parameters
nx, m = size(X, 1), size(X, 2)
ny, s = size(Y, 1), size(Y, 2)
if Xref === nothing Xref = X end
if Yref === nothing Yref = Y end
nrefx, mref = size(Xref, 1), size(Xref, 2)
nrefy, sref = size(Yref, 1), size(Yref, 2)
if nx != ny
throw(DimensionMismatch("number of rows in X and Y ($nx, $ny) are not equal"));
end
if nrefx != nrefy
throw(DimensionMismatch("number of rows in Xref and Yref ($nrefx, $nrefy) are not equal"));
end
if m != mref
throw(DimensionMismatch("number of columns in X and Xref ($m, $mref) are not equal"));
end
if s != sref
throw(DimensionMismatch("number of columns in Y and Yref ($s, $sref) are not equal"));
end
if (orient != :Input) && (orient != :Output) && (orient != :Graph)
throw(ArgumentError("`orient` must be :Input or :Output"));
end
# Check l
if l == 1
return deaholderl1(X, Y, weight = weight, orient = orient, rts = rts, slack = slack, Xref = Xref, Yref = Yref, names = names, optimizer = optimizer)
elseif l == 2
return deaholderl2(X, Y, weight = weight, orient = orient, rts = rts, slack = slack, Xref = Xref, Yref = Yref, names = names, optimizer = optimizer)
elseif l == Inf
return deaholderlinf(X, Y, weight = weight, orient = orient, rts = rts, slack = slack, Xref = Xref, Yref = Yref, names = names, optimizer = optimizer)
else
throw(ArgumentError("`l` must be 1, 2 or Inf"));
end
end
function deaholderl1(X::Union{Matrix,Vector}, Y::Union{Matrix,Vector};
weight::Bool = false,
orient::Symbol = :Graph, rts::Symbol = :CRS, slack = true,
Xref::Union{Matrix,Vector,Nothing} = nothing, Yref::Union{Matrix, Vector,Nothing} = nothing,
names::Union{Vector{<: AbstractString},Nothing} = nothing,
optimizer::Union{DEAOptimizer,Nothing} = nothing)::HolderL1DEAModel
# Get parameters
nx, m = size(X, 1), size(X, 2)
ny, s = size(Y, 1), size(Y, 2)
if Xref === nothing Xref = X end
if Yref === nothing Yref = Y end
nrefx, mref = size(Xref, 1), size(Xref, 2)
nrefy, sref = size(Yref, 1), size(Yref, 2)
# Compute DDF efficiency with 1 for each input and output
n = nx
nref = nrefx
effXY = fill(Inf, (n,m+s))
lambdaXY = Array{SparseMatrixCSC{Float64, Int64}}(undef, m + s)
if slack == true
slackXms = Array{Matrix}(undef, m + s)
slackYms = Array{Matrix}(undef, m + s)
end
if orient == :Graph || orient == :Input
for i = 1:m
Gxi = zeros(n,m)
if weight
Gxi[:,i] .= X[:,i]
else
Gxi[:,i] .= 1
end
tempdea = deaddf(X, Y, Gx = Gxi, Gy = :Zeros, rts = rts, Xref = Xref, Yref = Yref, optimizer = optimizer)
effXY[:,i] = efficiency(tempdea)
lambdaXY[i] = tempdea.lambda
if slack == true
slackXms[i] = slacks(tempdea, :X)
slackYms[i] = slacks(tempdea, :Y)
end
end
end
if orient == :Graph || orient == :Output
for i = 1:s
Gyi = zeros(n,s)
if weight
Gyi[:,i] .= Y[:,i]
else
Gyi[:,i] .= 1
end
tempdea = deaddf(X, Y, Gx = :Zeros, Gy = Gyi, rts = rts, Xref = Xref, Yref = Yref, optimizer = optimizer)
effXY[:,i + m] = efficiency(tempdea)
lambdaXY[i + m] = tempdea.lambda
if slack == true
slackXms[i + m] = slacks(tempdea, :X)
slackYms[i + m] = slacks(tempdea, :Y)
end
end
end
# Obtain the minimum efficiency scores
effi = zeros(n)
effmin = fill(0, n)
lambda = spzeros(n, nref)
if slack == true
slackX = zeros(n, m)
slackY = zeros(n, s)
else
slackX = Array{Float64}(undef, 0, 0)
slackY = Array{Float64}(undef, 0, 0)
end
for i = 1:n
effi[i], effmin[i] = findmin(effXY[i,:])
lambda[i, :] = lambdaXY[effmin[i]][i,:]
if slack == true
slackX[i,:] = slackXms[effmin[i]][i,:]
slackY[i,:] = slackYms[effmin[i]][i,:]
end
end
# Get X and Y targets
Xtarget = convert.(Float64, X[:,:])
Ytarget = convert.(Float64, Y[:,:])
for i = 1:n
if effmin[i] <= m
if weight
Xtarget[i,effmin[i]] = X[i,effmin[i]] - effi[i] * X[i,effmin[i]]
else
Xtarget[i,effmin[i]] = X[i,effmin[i]] - effi[i]
end
else
if weight
Ytarget[i,effmin[i] - m] = Y[i,effmin[i] - m] + effi[i] * Y[i,effmin[i] - m]
else
Ytarget[i,effmin[i] - m] = Y[i,effmin[i] - m] + effi[i]
end
end
end
if slack == true
Xtarget = Xtarget - slackX
Ytarget = Ytarget + slackY
end
return HolderL1DEAModel(n, m, s, orient, rts, weight, names, effi, effmin, slackX, slackY, lambda, Xtarget, Ytarget)
end
function deaholderl2(X::Union{Matrix,Vector}, Y::Union{Matrix,Vector};
weight::Bool = false,
orient::Symbol = :Graph, rts::Symbol = :CRS, slack = true,
Xref::Union{Matrix,Vector,Nothing} = nothing, Yref::Union{Matrix, Vector,Nothing} = nothing,
names::Union{Vector{<: AbstractString},Nothing} = nothing,
optimizer::Union{DEAOptimizer,Nothing} = nothing)::HolderL2DEAModel
# Check parameters
nx, m = size(X, 1), size(X, 2)
ny, s = size(Y, 1), size(Y, 2)
if Xref === nothing Xref = X end
if Yref === nothing Yref = Y end
nrefx, mref = size(Xref, 1), size(Xref, 2)
nrefy, sref = size(Yref, 1), size(Yref, 2)
if rts != :VRS
throw(ArgumentError("`rts` must be :VRS if l = 2"));
end
# Default optimizer
if optimizer === nothing
throw(ArgumentError("An optimizer that supports SOS constraints must be specified for Hölder l = 2"));
end
noSOS = false
# Compute efficiency for each DMU
n = nx
nref = n
effi = zeros(n)
lambdaeff = spzeros(n, nref)
Xtarget = convert.(Float64, X[:,:])
Ytarget = convert.(Float64, Y[:,:])
slackX = zeros(n, m)
slackY = zeros(n, s)
for i=1:n
# Value of inputs and outputs to evaluate
x0 = X[i,:]
y0 = Y[i,:]
# Create the optimization model
deamodel = newdeamodel(optimizer)
if orient == :Input || orient == :Graph
@variable(deamodel, Xeff[1:m] >= 0)
set_start_value.(Xeff, x0)
end
if orient == :Output || orient == :Graph
@variable(deamodel, Yeff[1:s] >= 0)
set_start_value.(Yeff, y0)
end
if weight
if orient == :Graph
@objective(deamodel, Min, sum(((x0[j] - Xeff[j]) / x0[j])^2 for j in 1:m) + sum(((Yeff[j] - y0[j]) / y0[j])^2 for j in 1:s) )
elseif orient == :Input
@objective(deamodel, Min, sum(((x0[j] - Xeff[j]) / x0[j])^2 for j in 1:m) )
elseif orient == :Output
@objective(deamodel, Min, sum(((Yeff[j] - y0[j]) / y0[j])^2 for j in 1:s) )
end
else
if orient == :Graph
@objective(deamodel, Min, sum((x0[j] - Xeff[j])^2 for j in 1:m) + sum((Yeff[j] - y0[j])^2 for j in 1:s) )
elseif orient == :Input
@objective(deamodel, Min, sum((x0[j] - Xeff[j])^2 for j in 1:m) )
elseif orient == :Output
@objective(deamodel, Min, sum((Yeff[j] - y0[j])^2 for j in 1:s) )
end
end
@variable(deamodel, gammatau[1:2, 1:nref] >= 0)
@variable(deamodel, vsX[1:2,1:m] >= 0)
@variable(deamodel, musY[1:2,1:s] >= 0)
@variable(deamodel, psi)
set_start_value.(gammatau[1,i], 1)
if orient == :Input || orient == :Graph
@constraint(deamodel, [j in 1:m], sum(Xref[t,j] * gammatau[1,t] for t in 1:nref) == Xeff[j] - vsX[2,j])
elseif orient == :Output
@constraint(deamodel, [j in 1:m], sum(Xref[t,j] * gammatau[1,t] for t in 1:nref) == x0[j] - vsX[2,j])
end
if orient == :Output || orient == :Graph
@constraint(deamodel, [j in 1:s], sum(Yref[t,j] * gammatau[1,t] for t in 1:nref) == Yeff[j] + musY[2,j])
elseif orient == :Input
@constraint(deamodel, [j in 1:s], sum(Yref[t,j] * gammatau[1,t] for t in 1:nref) == y0[j] + musY[2,j])
end
@constraint(deamodel, sum(gammatau[1,:]) == 1)
if orient == :Graph
@constraint(deamodel, sum(vsX[1,t] for t in 1:m) + sum(musY[1,t] for t in 1:s) == 1)
elseif orient == :Input
@constraint(deamodel, sum(vsX[1,t] for t in 1:m) == 1)
elseif orient == :Output
@constraint(deamodel, sum(musY[1,t] for t in 1:s) == 1)
end
@constraint(deamodel, [j in 1:nref], sum(vsX[1,t] * Xref[j,t] for t in 1:m) - sum(musY[1,t] * Yref[j,t] for t in 1:s) + psi - gammatau[2,j] == 0 )
# SOS constraints
try
@constraint(deamodel, [j in 1:nref], gammatau[:,j] in SOS1())
catch
@constraint(deamodel, [j in 1:nref], gammatau[1,j] * gammatau[2,j] == 0)
noSOS = true
end
try
@constraint(deamodel, [j in 1:m], vsX[:,j] in SOS1())
catch
@constraint(deamodel, [j in 1:m], vsX[1,j] * vsX[2,j] == 0)
noSOS = true
end
try
@constraint(deamodel, [j in 1:s], musY[:,j] in SOS1())
catch
@constraint(deamodel, [j in 1:s], musY[1,j] * musY[2,j] == 0)
noSOS = true
end
# Optimize and return results
JuMP.optimize!(deamodel)
effi[i] = JuMP.objective_value(deamodel)
if effi[i] < 0
# If objective function is negative change to 0, to avoid error when calculatin the square root.
@warn ("DMU $i negative value in the objective function: $(effi[i]). Replaced to 0")
effi[i] = 0
end
effi[i] = sqrt(effi[i] )
lambdaeff[i,:] = JuMP.value.(gammatau[1,:])
if orient == :Input || orient == :Graph
Xtarget[i,:] = JuMP.value.(Xeff)
end
if orient == :Output || orient == :Graph
Ytarget[i,:] = JuMP.value.(Yeff)
end
slackX[i,:] = JuMP.value.(vsX[2,:])
slackY[i,:] = JuMP.value.(musY[2,:])
if (termination_status(deamodel) != MOI.OPTIMAL) && (termination_status(deamodel) != MOI.LOCALLY_SOLVED)
@warn ("DMU $i termination status: $(termination_status(deamodel)). Primal status: $(primal_status(deamodel)). Dual status: $(dual_status(deamodel))")
end
end
# Add slacks to X and Y targets
Xtarget = Xtarget - slackX
Ytarget = Ytarget + slackY
# Display warning if model not solved with SOS constraints
if noSOS
@warn ("Model solved with $(optimizer.optimizer) is innacuate. Use a solver that supports SOS1 constraints.")
end
return HolderL2DEAModel(n, m, s, orient, rts, weight, names, effi, slackX, slackY, lambdaeff, Xtarget, Ytarget)
end
function deaholderlinf(X::Union{Matrix,Vector}, Y::Union{Matrix,Vector};
weight::Bool = true,
orient::Symbol = :Graph, rts::Symbol = :CRS, slack = true,
Xref::Union{Matrix,Vector,Nothing} = nothing, Yref::Union{Matrix, Vector,Nothing} = nothing,
names::Union{Vector{<: AbstractString},Nothing} = nothing,
optimizer::Union{DEAOptimizer,Nothing} = nothing)::HolderLInfDEAModel
# Get parameters
nx, m = size(X, 1), size(X, 2)
ny, s = size(Y, 1), size(Y, 2)
n = nx
# Buil direction based on orientation
if orient == :Input
if weight
Gx = X
else
Gx = :Ones
end
Gy = :Zeros
elseif orient == :Output
Gx = :Zeros
if weight
Gy = Y
else
Gy = :Ones
end
elseif orient == :Graph
if weight
Gx = X
Gy = Y
else
Gx = :Ones
Gy = :Ones
end
end
# Get DDF technical efficiency
tempdea = deaddf(X, Y, Gx = Gx, Gy = Gy, rts = rts, slack = slack, Xref = Xref, Yref = Yref, optimizer = optimizer)
effi = efficiency(tempdea)
lambda = peersmatrix(tempdea)
slackX = slacks(tempdea, :X)
slackY = slacks(tempdea, :Y)
Xtarget = targets(tempdea, :X)
Ytarget = targets(tempdea, :Y)
return HolderLInfDEAModel(n, m, s, orient, rts, weight, names, effi, slackX, slackY, lambda, Xtarget, Ytarget)
end
function efficiency(model::HolderL1DEAModel, type::Symbol)::Vector
if (type == :min) return model.effmin end
throw(ArgumentError("$(typeof(model)) with orienation has no efficiency $(type)"));
end
function Base.show(io::IO, x::HolderL1DEAModel)
compact = get(io, :compact, false)
n = nobs(x)
m = ninputs(x)
s = noutputs(x)
isweighted = x.isweighted
dmunames = names(x)
eff = efficiency(x)
effmin = x.effmin
slackX = slacks(x, :X)
slackY = slacks(x, :Y)
hasslacks = ! isempty(slackX)
effmincol = Array{String}(undef,n)
for i = 1:n
if effmin[i] <= m
effmincol[i] = "X" * string(Int(effmin[i]))
else
effmincol[i] = "Y" * string(Int(effmin[i] - m))
end
end
if !compact
print(io, "Hölder L1 DEA Model \n")
print(io, "DMUs = ", n)
print(io, "; Inputs = ", m)
print(io, "; Outputs = ", s)
print(io, "\n")
print(io, "Orientation = ", string(x.orient))
print(io, "; Returns to Scale = ", string(x.rts))
print(io, "\n")
if isweighted
print(io, "Weighted (weakly) Hölder distance function \n")
end
if hasslacks == true
show(io, CoefTable(hcat(eff, effmincol, slackX, slackY), ["efficiency"; "minimum"; ["slackX$i" for i in 1:m ]; ["slackY$i" for i in 1:s ]], dmunames))
else
show(io, CoefTable(hcat(eff, effmincol), ["efficiency"; "minimum"], dmunames))
end
end
end
function Base.show(io::IO, x::HolderL2DEAModel)
compact = get(io, :compact, false)
n = nobs(x)
m = ninputs(x)
s = noutputs(x)
isweighted = x.isweighted
dmunames = names(x)
eff = efficiency(x)
slackX = slacks(x, :X)
slackY = slacks(x, :Y)
hasslacks = ! isempty(slackX)
if !compact
print(io, "Hölder L2 DEA Model \n")
print(io, "DMUs = ", n)
print(io, "; Inputs = ", m)
print(io, "; Outputs = ", s)
print(io, "\n")
print(io, "Orientation = ", string(x.orient))
print(io, "; Returns to Scale = ", string(x.rts))
print(io, "\n")
if isweighted
print(io, "Weighted (weakly) Hölder distance function \n")
end
show(io, CoefTable(hcat(eff, slackX, slackY), ["efficiency"; ["slackX$i" for i in 1:m ]; ["slackY$i" for i in 1:s ]], dmunames))
end
end
function Base.show(io::IO, x::HolderLInfDEAModel)
compact = get(io, :compact, false)
n = nobs(x)
m = ninputs(x)
s = noutputs(x)
isweighted = x.isweighted
dmunames = names(x)
eff = efficiency(x)
slackX = slacks(x, :X)
slackY = slacks(x, :Y)
hasslacks = ! isempty(slackX)
if !compact
print(io, "Hölder LInf DEA Model \n")
print(io, "DMUs = ", n)
print(io, "; Inputs = ", m)
print(io, "; Outputs = ", s)
print(io, "\n")
print(io, "Orientation = ", string(x.orient))
print(io, "; Returns to Scale = ", string(x.rts))
print(io, "\n")
if isweighted
print(io, "Weighted (weakly) Hölder distance function \n")
end
if hasslacks == true
show(io, CoefTable(hcat(eff, slackX, slackY), ["efficiency"; ["slackX$i" for i in 1:m ]; ["slackY$i" for i in 1:s ]], dmunames))
else
show(io, CoefTable(hcat(eff), ["efficiency"], dmunames))
end
end
end