Can you close that issue and open it in SDDP.jl instead? GitHub · Where software is built
Sorry, it’s done: Performance issue with value in JuMP when using SDDP · Issue #869 · odow/SDDP.jl · GitHub
1 Like
I just took a slight look at the procedure of JuMP.value. To be honest, it’s highly complex—no wonder it’s slow.
I would like to ask if JuMP could provide a pass-by-address method, besides the existing pass-by-value method?
I wrote one
import JuMP, Gurobi
const GRB_ENV = Gurobi.Env()
function unsafe_array_value!(X, x)
o = JuMP.unsafe_backend(JuMP.owner_model(x[1]))
for j = eachindex(X)
Gurobi.GRBgetdblattrelement(o, "X",
Gurobi.c_column(o, JuMP.index(x[j])),
view(X, j))
end
end
function unsafe_test()
m = JuMP.Model(() -> Gurobi.Optimizer(GRB_ENV))
JuMP.@variable(m, x[i = 1:3, j = 1:4] >= i + j)
JuMP.@objective(m, Min, sum(x))
JuMP.optimize!(m)
X = similar(x, Float64)
unsafe_array_value!(X, x)
X
end
unsafe_test()
For direct_model, do this
function scalar!(X, j, m, x)
o = m.moi_backend
Gurobi.GRBgetdblattrelement(o, "X", Gurobi.c_column(o, JuMP.index(x[j])), view(X, j))
end;
value!(X, m, x) = foreach(j -> scalar!(X, j, m, x), eachindex(X));
function solve_mst_and_up_value!(model, s, θ, β)
JuMP.optimize!(model)
JuMP.termination_status(model) == JuMP.OPTIMAL || error()
s.ub.x = JuMP.objective_bound(model)
value!(s.β, model, β)
value!(s.θ, model, θ)
end;
function shot!(ref; model = model, θ = θ, β = β, rEF = rEF)
s = rEF.x
@lock mst_lock solve_mst_and_up_value!(model, s, θ, β) # `s` gets updated/mutated here
s_tmp = ref.x
setfield!(ref, :x, s) # the ref gets updated here
setfield!(rEF, :x, s_tmp)
end;
This part is not correct.
You cannot connect the JuMP.index to an unsafe_backend.
Please read the docstring of unsafe_backend very carefully:
I have no plans to change how value works. I’m yet to see any evidence that it is a performance bottleneck in any real problem.
1 Like
What is this API for? (JuMP → src → variables.jl)
function value(var_value::Function, v::GenericVariableRef)
return var_value(v)
end
It’s weird:
julia> x
_[1]
julia> JuMP.value(x)
0.0
julia> JuMP.value(-, x) # What is the purpose of this??
-_[1]
Why not write -x directly?
For stuff like value(start_value, x), or for querying value(var_value, expr) where var_value is not value. vaule(-, x) is not the intended use-case. The weak assumption is really var_value(::AbstractVariableRef)::Real.
Originally added in Add `value` for constraints with a custom `var_value` function by dourouc05 · Pull Request #2317 · jump-dev/JuMP.jl · GitHub
I was a bit surprised that your Ref method has zero allocate. Why?
But my value! method is still 2x faster than (MOI’s method + @.) (see the Comments of the last 2 lines):
import MathOptInterface as MOI
import Gurobi
using BenchmarkTools
# Retrieve by MOI
const VP = MOI.VariablePrimal();
pass_by_value_scalar(o, x) = MOI.get(o, VP, x);
# Retrieve by my method
function scalar!(X, j, o, x)
Gurobi.GRBgetdblattrelement(o, "X",
Gurobi.c_column(o, x[j]),
view(X, j)
)
end;
value!(X, o, x) = foreach(j -> scalar!(X, j, o, x), eachindex(X));
# Test begins
o = Gurobi.Optimizer();
N = 9999;
x = MOI.add_variables(o, N); X = similar(x, Float64);
MOI.optimize!(o)
@btime value!($X, $o, $x) # 406.750 μs (0 allocations: 0 bytes)
@btime @. $X = pass_by_value_scalar($o, $x); # 870.432 μs (0 allocations: 0 bytes)
This is the best solution that I can come up with. I don’t know if anyone can figure out a faster one.
Another question, @odow:
what is the purpose of this wrapper?
# !!! danger
# We want to avoid being too specific in the type arguments to avoid method
# ambiguity. For example,
# get(::ModelLike, ::AbstractVariableAttribute, ::Vector{VariableIndex})
# would not allow us to define get(::SomeModel, ::AnyAttribute, ::Vector).
function get(model::ModelLike, attr::AnyAttribute, idxs::Vector)
if isempty(idxs)
return Vector{attribute_value_type(attr)}()
end
return get.(model, attr, idxs)
end
Why not let users themselves write get.?
The original intent was to allow solvers to provide an efficient way of querying multiple values in the same API call.We didn’t make much progress on this and most people use get. or value. because querying the value is essentially never a bottleneck in any real-world application.
Write nice, simple, high-level algorithms. Don’t worry about the low-level details and optimisations. You’ll find that most opportunities for improvements come from higher-level agorithmic choices. If your code is spending all its time querying a result value, you’re doing something wrong.
After gaining some experiences, I have to admit you’re right. So I won’t go back to use MOI because it will cost me lots of effort I can ill afford.
And I come up with this solution, where the master model’s solution is exposed to the outside world via solution_ref.x only. The master needs to keeps two initially allocated solution containers (=tuple).
Do you think this idea is advisable?
function retrieve!(solution_mutated, model) # this is allocate free (I hope so)
x = model[:x]
y = model[:y]
@. solution_mutated.x = value(x)
@. solution_mutated.y = value(y)
nothing
end
# build the master model
model = JuMP.Model()
@variable(model, x[1:3])
@variable(model, y[1:4])
@expression(model, e, sum(x) + sum(y))
@constraint(model, e >= 0)
@objective(model, Min, e)
# define these 3 containers
solution_exposed = (x = fill(0., 3), y = fill(0., 4))
solution_mutated = (x = fill(2., 3), y = fill(2., 4))
solution_ref = Ref(solution_exposed)
# do some iterations
@constraint(model, e >= 1) # add a new cut
optimize!(model)
retrieve!(solution_mutated, model) # get the updated solution _without_ vector allocating
solution_exposed, solution_mutated = solution_mutated, solution_exposed # __swap__
setfield!(solution_ref, :x, solution_exposed) # update the referrence (in an instant)
@constraint(model, e >= 2) # add a new cut
optimize!(model)
retrieve!(solution_mutated, model)
If you think this idea is correct, then probably this can serve as an answer to this topic.
Why so complicated? Until you benchmark that the allocations are a problem, go for the simplest possible thing.
get_solution(model) = (; x = value(model[:x]), y = value(model[:y]))
model = JuMP.Model()
@variable(model, x[1:3])
@variable(model, y[1:4])
@expression(model, e, sum(x) + sum(y))
@constraint(model, e >= 0)
@objective(model, Min, e)
solution_ref = Ref{Any}()
optimize!(model)
solution_ref[] = get_solution(model)
@constraint(model, e >= 1)
optimize!(model)
solution_ref[] = get_solution(model)
@constraint(model, e >= 2)
optimize!(model)
solution_ref[] = get_solution(model)
1 Like