tkoolen / parametron.jl Goto Github PK
View Code? Open in Web Editor NEWEfficiently solving instances of a parameterized family of (possibly mixed-integer) linear/quadratic optimization problems in Julia
License: Other
Efficiently solving instances of a parameterized family of (possibly mixed-integer) linear/quadratic optimization problems in Julia
License: Other
e.g.:
julia> using SimpleQP
julia> using OSQP.MathOptInterfaceOSQP
julia> model = Model(OSQPOptimizer())
Model{Float64, OSQP.MathOptInterfaceOSQP.OSQPOptimizer}(…)
julia> solve!(model)
ERROR: Copy failed with status CopyOtherError. Message: Unsupported objective
Stacktrace:
[1] initialize!(::SimpleQP.Model{Float64,OSQP.MathOptInterfaceOSQP.OSQPOptimizer}) at /home/rdeits/locomotion/explorations/simpleqp-mpc/packages/v0.6/SimpleQP/src/model.jl:117
[2] solve!(::SimpleQP.Model{Float64,OSQP.MathOptInterfaceOSQP.OSQPOptimizer}) at /home/rdeits/locomotion/explorations/simpleqp-mpc/packages/v0.6/SimpleQP/src/model.jl:152
I posted my question here, but figured out this could be an issue with Parametron
.
I am using Julia version 1.3.1 (2019-12-30) on Mac OS, official binary downloaded from Julia website.
I added Parametron
using Julia's Pkg interface. It reported that the version of Parametron
installed was 0.8.0. However, on the github repo and in the Julia's registries, the current version is 0.9.1. No matter how I tried (removed the package then reinstalled it), the version is always 0.8.0. I have MathOptInterface v0.9.9.
What could be wrong? And how do I resolve the problem and get the latest version of Parametron?
Hey,
after skimming the documentation and the source code i could find no reference to generic non-linear optimization, is that on purpose, did i miss something or wasn't it required yet in your applications?
It looks like we're not quite handling literal_pow
correctly:
julia> using SimpleQP
julia> model = SimpleQP.MockModel()
SimpleQP.MockModel(SimpleQP.Parameter[], Base.RefValue{Int64}(1))
julia> p = Parameter(() -> 1.0, model)
Parameter{Float64, …}(…)
julia> e = @expression(p^2)
LazyExpression{FunctionWrapper{…}(LazyExpression{Base.#literal_pow, …}(…))}(…)
julia> e()
1.0
julia> SimpleQP.findallocs(e)
LazyExpression{FunctionWrapper{…}(LazyExpression{Base.#literal_pow, …}(…))}(…): 48 bytes
[1]: Base.#^
[2]: Parameter{Float64, …}(…): 0 bytes
[3]: DataType
It looks like constraint modification changed in jump-dev/MathOptInterface.jl#388 , and that broke SimpleQP's interface (I noticed this because the compatible version of Gurobi.jl needs MOI master). @tkoolen are you already working on this or is that something i can help with?
julia> using SimpleQP, Compat.LinearAlgebra
julia> dot(Variable[], Variable[])
ERROR: MethodError: no method matching zero(::Type{Variable})
Closest candidates are:
zero(::Type{LibGit2.GitHash}) at /Users/osx/buildbot/slave/package_osx64/build/usr/share/julia/stdlib/v0.7/LibGit2/src/oid.jl:210
zero(::Type{OldPkg.Resolve.VersionWeights.VWPreBuildItem}) at /Users/osx/buildbot/slave/package_osx64/build/usr/share/julia/stdlib/v0.7/OldPkg/src/resolve/versionweight.jl:80
zero(::Type{OldPkg.Resolve.VersionWeights.VWPreBuild}) at /Users/osx/buildbot/slave/package_osx64/build/usr/share/julia/stdlib/v0.7/OldPkg/src/resolve/versionweight.jl:116
...
Stacktrace:
[1] dot(::Array{Variable,1}, ::Array{Variable,1}) at /Users/osx/buildbot/slave/package_osx64/build/usr/share/julia/stdlib/v0.7/LinearAlgebra/src/generic.jl:719
[2] top-level scope at none:0
@objective
escapes its entire result, which assumes that setobjective!
and @expression
are defined in the caller's scope:
julia> module Foo
using SimpleQP: @objective, Model, Minimize
using OSQP.MathOptInterfaceOSQP: OSQPOptimizer
m = Model(OSQPOptimizer())
@objective(m, Minimize, 0)
end
WARNING: replacing module Foo
ERROR: UndefVarError: @expression not defined
Since all of the expressions are wrapped in a FunctionWrapper, findallocs
no longer prints out a nice tree. It looks like we can instead recurse into expr.f.obj[].args
instead of expr.args
to get the tree printing back.
Hi,
I am trying out Parametron v0.9.0 on Julia 1.1. Just running through the examples in the README.md and in the docs, and I can't seem to get the zero-allocations shown in any of the examples:
using OSQP
optimizer = OSQP.Optimizer()
# create a Parametron.Model, which holds problem information
using Parametron
using Random, LinearAlgebra
model = Model(optimizer)
# create decision variables and parameters
n = 8; m = 2
x = [Variable(model) for _ = 1 : n]
A = Parameter(rand!, zeros(n, n), model)
b = Parameter(rand!, zeros(n), model)
C = Parameter(rand!, zeros(m, n), model)
d = Parameter(zeros(m), model) do d
# do syntax makes it easy to create custom Parameters
rand!(d)
d .*= 2
end
# the @expression macro can be used to create 'lazy' expressions,
# which can be used in constraints or the objective function, and
# can be evaluated at a later time, automatically updating the
# Parameters in the process (if needed).
residual = @expression A * x - b
# set the objective function
@objective(model, Minimize, residual ⋅ residual)
# add the constraints. You could have multiple @constraint calls
# as well. ==, <=, and >= are supported.
@constraint(model, C * x == d)
solve!(model)
value.(model, x)
solve!(model)
using BenchmarkTools
@btime solve!($model)
gives
201.557 μs (12 allocations: 384 bytes)
Is there some change needed to reproduce the zero allocations?
Trying to constrain a scalar variable causes a few different errors:
optimizer = OSQP.MathOptInterfaceOSQP.OSQPOptimizer()
model = Model(optimizer)
x = Variable(model)
@constraint model x >= 0
solve!(model)
MethodError: no method matching length(::SimpleQP.Functions.AffineFunction{Int64})
Closest candidates are:
length(::SimpleVector) at essentials.jl:256
length(::Base.MethodList) at reflection.jl:558
length(::MethodTable) at reflection.jl:634
...
Stacktrace:
[1] add_nonnegative_constraint!(::SimpleQP.Model{Float64,OSQP.MathOptInterfaceOSQP.OSQPOptimizer}, ::SimpleQP.LazyExpression{FunctionWrappers.FunctionWrapper{SimpleQP.Functions.AffineFunction{Int64},Tuple{}},Tuple{}}) at /home/rdeits/locomotion/explorations/simple-qp-miqp/packages/v0.6/SimpleQP/src/model.jl:85
Trying to wrap the that scalar as a vector causes a different error:
@constraint model [x] >= 0
solve!(model)
MethodError: no method matching vecsubtract!(::Array{SimpleQP.Functions.AffineFunction{Int64},1}, ::Array{SimpleQP.Functions.Variable,1}, ::Int64)
Closest candidates are:
vecsubtract!(::AbstractArray{SimpleQP.Functions.AffineFunction{T},1}, ::AbstractArray{T,1} where T, ::AbstractArray{T,1} where T) where T at /home/rdeits/locomotion/explorations/simple-qp-miqp/packages/v0.6/SimpleQP/src/functions.jl:419
Stacktrace:
[1] macro expansion at /home/rdeits/locomotion/explorations/simple-qp-miqp/packages/v0.6/SimpleQP/src/lazyexpression.jl:20 [inlined]
[2] evalexpr at /home/rdeits/locomotion/explorations/simple-qp-miqp/packages/v0.6/SimpleQP/src/lazyexpression.jl:17 [inlined]
[3] LazyExpression at /home/rdeits/locomotion/explorations/simple-qp-miqp/packages/v0.6/SimpleQP/src/lazyexpression.jl:14 [inlined]
[4] wrap(::SimpleQP.LazyExpression{SimpleQP.Functions.#vecsubtract!,Tuple{Array{SimpleQP.Functions.AffineFunction{Int64},1},SimpleQP.LazyExpression{FunctionWrappers.FunctionWrapper{Array{SimpleQP.Functions.Variable,1},Tuple{}},Tuple{}},Int64}}) at /home/rdeits/locomotion/explorations/simple-qp-miqp/packages/v0.6/SimpleQP/src/lazyexpression.jl:146
and wrapping both x
and 0
produces a conversion error:
@constraint model [x] >= [0]
solve!(model)
MethodError: Cannot `convert` an object of type SimpleQP.Functions.AffineFunction{Int64} to an object of type SimpleQP.Functions.AffineFunction{Float64}
This may have arisen from a call to the constructor SimpleQP.Functions.AffineFunction{Float64}(...),
since type constructors fall back to convert methods.
Stacktrace:
[1] copy!(::IndexLinear, ::Array{SimpleQP.Functions.AffineFunction{Float64},1}, ::IndexLinear, ::Array{SimpleQP.Functions.AffineFunction{Int64},1}) at ./abstractarray.jl:656
[2] convert(::Type{Array{SimpleQP.Functions.AffineFunction{Float64},1}}, ::Array{SimpleQP.Functions.AffineFunction{Int64},1}) at ./array.jl:368
[3] macro expansion at /home/rdeits/locomotion/explorations/simple-qp-miqp/packages/v0.6/SimpleQP/src/lazyexpression.jl:20 [inlined]
[4] evalexpr at /home/rdeits/locomotion/explorations/simple-qp-miqp/packages/v0.6/SimpleQP/src/lazyexpression.jl:17 [inlined]
[5] LazyExpression at /home/rdeits/locomotion/explorations/simple-qp-miqp/packages/v0.6/SimpleQP/src/lazyexpression.jl:14 [inlined]
[6] optimize(::SimpleQP.LazyExpression{Base.#convert,Tuple{SimpleQP.LazyExpression{FunctionWrappers.FunctionWrapper{DataType,Tuple{}},Tuple{}},SimpleQP.LazyExpression{FunctionWrappers.FunctionWrapper{Array{SimpleQP.Functions.AffineFunction{Int64},1},Tuple{}},Tuple{}}}}, ::Type{T} where T, ::Type{Array{SimpleQP.Functions.AffineFunction{Int64},1}}) at /home/rdeits/locomotion/explorations/simple-qp-miqp/packages/v0.6/SimpleQP/src/lazyexpression.jl:135
[7] optimize_toplevel(::SimpleQP.LazyExpression{Base.#convert,Tuple{SimpleQP.LazyExpression{FunctionWrappers.FunctionWrapper{DataType,Tuple{}},Tuple{}},SimpleQP.LazyExpression{FunctionWrappers.FunctionWrapper{Array{SimpleQP.Functions.AffineFunction{Int64},1},Tuple{}},Tuple{}}}}) at /home/rdeits/locomotion/explorations/simple-qp-miqp/packages/v0.6/SimpleQP/src/lazyexpression.jl:55
[8] addconstraint!(::SimpleQP.Model{Float64,OSQP.MathOptInterfaceOSQP.OSQPOptimizer}, ::SimpleQP.LazyExpression{FunctionWrappers.FunctionWrapper{Array{SimpleQP.Functions.AffineFunction{Int64},1},Tuple{}},Tuple{}}, ::MathOptInterface.Nonnegatives) at /home/rdeits/locomotion/explorations/simple-qp-miqp/packages/v0.6/SimpleQP/src/model.jl:77
[9] add_nonnegative_constraint!(::SimpleQP.Model{Float64,OSQP.MathOptInterfaceOSQP.OSQPOptimizer}, ::SimpleQP.LazyExpression{FunctionWrappers.FunctionWrapper{Array{SimpleQP.Functions.AffineFunction{Int64},1},Tuple{}},Tuple{}}) at /home/rdeits/locomotion/explorations/simple-qp-miqp/packages/v0.6/SimpleQP/src/model.jl:85
FunctionWrappers is still broken on 1.0, and telling people to check out the tk/julia-0.7-quickfix
on my branch is not a real solution. I think to keep things moving for now we should just incorporate that version of FunctionWrappers into Parametron (it's just a single file). @rdeits, thoughts?
It's time 😉
I think "SimpleQP" is underselling how useful this package is (and not necessarily accurate if we decide to add MIPs in the future).
When defining constraints like @constraint(m, v >= 0)
for Vector
-valued v
on Julia 0.7 we get this deprecation warning
┌ Warning: `a::AbstractArray - b::Number` is deprecated, use `a .- b` instead.
│ caller = optimize_toplevel at lazyexpression.jl:59 [inlined]
└ @ Core ~/.julia/packages/SimpleQP/pHx4o/src/lazyexpression.jl:59
This comes from the minus expressions here but naively replacing those with broadcast
ed equivalents interferes with the current LazyExpression
optimizations. I feel like the answer is to either integrate base broadcasting into LazyExpression
s (a potentially daunting task) or to make a localized fix (e.g., maybe turn just this line into a broadcasted equivalent?); what do you think?
The mixed-integer tests are closer to working with the new LinQuadOptInterface tag, but GLPK master is still needed (probably because of jump-dev/GLPK.jl#71).
Currently, we only optimize when the vector has eltype <:Union{Variable, AffineFunction}
.
I've noticed that because everything inside an @constraint()
call gets LazyExpression-wrapped, there are some weird issues with referential transparency. For example, this fails:
x = Variable(model)
@constraint model [x] <= [0.]
with ArgumentError: Unhandled expression head: vect.
. But this works:
x = Variable(model)
lhs = [x]
rhs = [0.]
@constraint model lhs <= rhs
I realize this is a necessary result of the whole expression-wrapping stuff, but it would be nice to be able to work around it for things that don't actually need to be lazy-ified.
We could do that by supporting $
as a way to mark expressions which should just be evaluated rather than being wrapped. That would look like:
@constraint model $([x]) <= $([0.])
Here's a sample implementation of the macro that preserves any $
-interpolated expressions rather than wrapping them:
lazy_wrap(x) = esc(x)
function lazy_wrap(x::Expr)
if x.head == :line
return esc(x)
elseif x.head == :$ && length(x.args) == 1
return esc(x.args[1])
elseif x.head ∈ [:call, :block]
x.args .= lazy_wrap.(x.args)
if x.head == :call
return :(wrap(SimpleQP.optimize_toplevel(SimpleQP.LazyExpression($(x.args...)))))
else
return esc(x)
end
else
buf = IOBuffer()
dump(buf, x)
msg =
"""
Unhandled expression head: $(x.head). expr:
$(String(take!(buf)))
"""
return :(throw(ArgumentError($msg)))
end
end
macro expression(expr)
lazy_wrap(expr)
end
(I couldn't figure out how to do this with postwalk
cleanly, sadly)
I noticed when trying to build a large Parametron model that I was spending an unreasonable amount of time just constructing the initial model. From profiling, it looks like almost all of the time is spent inside the FunctionWrappersQuickFix.FunctionWrapper
constructor.
I was curious to see if the quick fix might be slowing things down, and it does appear that what we're doing is dramatically slower than the normal FunctionWrapper. For example:
julia> using FunctionWrappers, Parametron, BenchmarkTools
julia> function f1(x)
Parametron.FunctionWrappersQuickFix.FunctionWrapper{Int, Tuple{Int}}(y -> y + x)
end
f1 (generic function with 1 method)
julia> function f2(x)
FunctionWrappers.FunctionWrapper{Int, Tuple{Int}}(y -> y + x)
end
f2 (generic function with 1 method)
julia> @btime f1($1)
188.514 ns (2 allocations: 64 bytes)
Parametron.FunctionWrappersQuickFix.FunctionWrapper{Int64,Tuple{Int64}}(Ptr{Nothing} @0x00007f9498be3570, Ptr{Nothing} @0x00007f9493ed8010, Base.RefValue{getfield(Main, Symbol("##7#8")){Int64}}(getfield(Main, Symbol("##7#8")){Int64}(1)), getfield(Main, Symbol("##7#8")){Int64})
julia> @btime f2($1)
9.602 ns (2 allocations: 64 bytes)
FunctionWrappers.FunctionWrapper{Int64,Tuple{Int64}}(Ptr{Nothing} @0x00007f9498beada0, Ptr{Nothing} @0x00007f9493ed8020, Base.RefValue{getfield(Main, Symbol("##9#10")){Int64}}(getfield(Main, Symbol("##9#10")){Int64}(1)), getfield(Main, Symbol("##9#10")){Int64})
julia> versioninfo()
Julia Version 1.1.0
Commit 80516ca202 (2019-01-21 21:24 UTC)
Platform Info:
OS: Linux (x86_64-pc-linux-gnu)
CPU: Intel(R) Xeon(R) CPU E3-1535M v6 @ 3.10GHz
WORD_SIZE: 64
LIBM: libopenlibm
LLVM: libLLVM-6.0.1 (ORCJIT, skylake)
so FunctionWrappers.FunctionWrapper
is about 20x faster to construct for this example.
I'm not sure what exactly we can do about it, other than to keep making noise over at JuliaLang/julia#28669
I noticed a fair number of operations that don't yet have non-allocating optimizations, so i thought I'd put together a small example that shows the various issues:
function make_model()
optimizer = defaultoptimizer()
model = Model(optimizer)
v = [Variable(model) for _ in 1:2]
v0 = zeros(2)
Δt = 0.01
u = [Variable(model) for _ in 1:2]
H = Parameter(zeros(2, 2), model) do H
H[1, 1] = rand()
H[2, 2] = rand()
end
c = Parameter(zeros(2), model) do c
c[1] = rand()
c[2] = rand()
end
@constraint(model, H * (v - v0) == Δt * (u - c))
model
end
model = make_model()
SimpleQP.findallocs(model.zeroconstraints[1].expr)
which currently shows 3616 bytes.
I'll see if I can take a shot at fixing these.
This works:
model = Model(optimizer)
n = 2
x = [Variable(model) for _ = 1:n]
q = Parameter(zeros(n), model) do q
q[1] = rand()
q[2] = rand()
end
@objective(model, Minimize, transpose(x) * eye(2) * x + q' * x)
but reversing the order of operations in the objective like so:
@objective(model, Minimize, q' * x + transpose(x) * eye(2) * x)
gives:
MethodError: Cannot `convert` an object of type SimpleQP.Functions.AffineFunction{Float64} to an object of type SimpleQP.Functions.QuadraticFunction{Float64}
This may have arisen from a call to the constructor SimpleQP.Functions.QuadraticFunction{Float64}(...),
since type constructors fall back to convert methods.
Stacktrace:
[1] macro expansion at /home/rdeits/locomotion/explorations/simple-qp-miqp/packages/v0.6/SimpleQP/src/lazyexpression.jl:20 [inlined]
[2] evalexpr at /home/rdeits/locomotion/explorations/simple-qp-miqp/packages/v0.6/SimpleQP/src/lazyexpression.jl:17 [inlined]
[3] LazyExpression at /home/rdeits/locomotion/explorations/simple-qp-miqp/packages/v0.6/SimpleQP/src/lazyexpression.jl:14 [inlined]
[4] optimize(::SimpleQP.LazyExpression{Base.#+,Tuple{SimpleQP.LazyExpression{FunctionWrappers.FunctionWrapper{SimpleQP.Functions.AffineFunction{Float64},Tuple{}},Tuple{}},SimpleQP.LazyExpression{FunctionWrappers.FunctionWrapper{SimpleQP.Functions.QuadraticFunction{Float64},Tuple{}},Tuple{}}}}, ::DataType, ::DataType) at /home/rdeits/locomotion/explorations/simple-qp-miqp/packages/v0.6/SimpleQP/src/lazyexpression.jl:93
[5] optimize_toplevel(::SimpleQP.LazyExpression{Base.#+,Tuple{SimpleQP.LazyExpression{FunctionWrappers.FunctionWrapper{SimpleQP.Functions.AffineFunction{Float64},Tuple{}},Tuple{}},SimpleQP.LazyExpression{FunctionWrappers.FunctionWrapper{SimpleQP.Functions.QuadraticFunction{Float64},Tuple{}},Tuple{}}}}) at /home/rdeits/locomotion/explorations/simple-qp-miqp/packages/v0.6/SimpleQP/src/lazyexpression.jl:55
Creating a model with parameters whose initial values are NaN
seems to break the model, even if those parameters are changed to non-NaN before any call to solve!
. For example:
# Set up the model
m = Model(OSQP.Optimizer())
# Create an externally-controlled value. This is initially NaN, but
# we will set it to non-NaN later
outer_value = Ref(NaN)
# Create a parameter which just reads from our outer value
param = Parameter(m) do
outer_value[]
end
# Constrain that the variable x == param. At this point, outer_value[] == NaN
x = Variable(m)
@constraint m x == param
# Fix our outer_value so that it is non-NaN
outer_value[] = 1.0
# Solve. All the parameters are now non-NaN, but the solver still fails:
solve!(m)
objectivevalue(m)
gives an objective value of NaN
.
Moving the outer_value[] = 1.0
to before the call to @constraint
fixes the problem, but shouldn't be necessary as far as I know.
I'm seeing errors like the following frequently, and before diving in I'm curious to know if others have seen the same:
ERROR: LoadError: Gurobi.GurobiError(10003, "Element 1 of a double array is Nan.")
Stacktrace:
[1] add_constrs!(::Gurobi.Model, ::Array{Int32,1}, ::Array{Int32,1}, ::Array{Float64,1}, ::Array{Int8,1}, ::Array{Float64,1}) at /home/tim/.julia/dev/Gurobi/src/grb_constrs.jl:63
[2] add_constrs!(::Gurobi.Model, ::Array{Int64,1}, ::Array{Int64,1}, ::Array{Float64,1}, ::Array{Int8,1}, ::Array{Float64,1}) at /home/tim/.julia/dev/Gurobi/src/grb_constrs.jl:71
[3] add_linear_constraints! at /home/tim/.julia/dev/Gurobi/src/MOIWrapper.jl:131 [inlined]
[4] add_linear_constraint(::Gurobi.Optimizer, ::MathOptInterface.VectorAffineFunction{Float64}, ::Int8) at /home/tim/.julia/packages/LinQuadOptInterface/WqyOh/src/constraints/vectoraffine.jl:58
[5] addconstraint!(::Gurobi.Optimizer, ::MathOptInterface.VectorAffineFunction{Float64}, ::MathOptInterface.Nonnegatives) at /home/tim/.julia/packages/LinQuadOptInterface/WqyOh/src/constraints/vectoraffine.jl:12
[6] copyconstraints!(::Gurobi.Optimizer, ::Parametron.ParametronMOIModel{Float64}, ::Bool, ::MathOptInterface.Utilities.IndexMap, ::Type{MathOptInterface.VectorAffineFunction{Float64}}, ::Type{MathOptInterface.Nonnegatives}) at /home/tim/.julia/packages/MathOptInterface/QfYLx/src/Utilities/copy.jl:78
[7] defaultcopy!(::Gurobi.Optimizer, ::Parametron.ParametronMOIModel{Float64}, ::Bool) at /home/tim/.julia/packages/MathOptInterface/QfYLx/src/Utilities/copy.jl:111
[8] #copy!#71 at /home/tim/.julia/packages/LinQuadOptInterface/WqyOh/src/copy.jl:5 [inlined]
[9] copy! at /home/tim/.julia/packages/LinQuadOptInterface/WqyOh/src/copy.jl:4 [inlined]
[10] initialize!(::Parametron.Model{Float64,Gurobi.Optimizer}) at /home/tim/.julia/dev/Parametron/src/model.jl:114
[11] solve! at /home/tim/.julia/dev/Parametron/src/model.jl:149 [inlined]
[12] #fit_quadratic_lowerbound!#57(::Array{Float64,2}, ::Array{Float64,2}, ::Array{Float64,1}, ::Function, ::CoordinateSplittingPTrees.SymmetricArray{Float64,2,Array{Float64,2}}, ::getfield(CoordinateSplittingPTrees, Symbol("#updater#54")){typeof(value),Int64,Int64,DataStructures.PriorityQueue{Box{2,Float64,Float64,3,Float64,Interval{Float64}},Float64,Base.Order.ForwardOrdering},getfield(CoordinateSplittingPTrees, Symbol("#maybe_enqueue#53")){DataStructures.PriorityQueue{Box{2,Float64,Float64,3,Float64,Interval{Float64}},Float64,Base.Order.ForwardOrdering}},CoordinateSplittingPTrees.SymmetricArray{Float64,2,Array{Float64,2}},Array{Float64,1},Array{Float64,1},Array{Bool,1},Array{Float64,2},Array{Float64,2},Array{Float64,2},Array{Float64,1}}, ::Parametron.Model{Float64,Gurobi.Optimizer}, ::Array{Parametron.Functions.Variable,1}, ::Array{Parametron.Functions.Variable,1}, ::Box{2,Float64,Float64,3,Float64,Interval{Float64}}, ::Array{Float64,1}) at /home/tim/.julia/dev/CoordinateSplittingPTrees/src/cs2.jl:472
[13] (::getfield(CoordinateSplittingPTrees, Symbol("#kw##fit_quadratic_lowerbound!")))(::NamedTuple{(:X, :X2, :Δf),Tuple{Array{Float64,2},Array{Float64,2},Array{Float64,1}}}, ::typeof(CoordinateSplittingPTrees.fit_quadratic_lowerbound!), ::CoordinateSplittingPTrees.SymmetricArray{Float64,2,Array{Float64,2}}, ::Function, ::Parametron.Model{Float64,Gurobi.Optimizer}, ::Array{Parametron.Functions.Variable,1}, ::Array{Parametron.Functions.Variable,1}, ::Box{2,Float64,Float64,3,Float64,Interval{Float64}}, ::Array{Float64,1}) at ./none:0 (repeats 2 times)
[14] (::getfield(OptimizeIVQ, Symbol("#newton#7")){typeof(f),Box{2,Float64,Float64,3,Float64,Interval{Float64}},World{Float64,Float64,Float64},DataStructures.PriorityQueue{Box{2,Float64,Float64,3,Float64,Interval{Float64}},Float64,Base.Order.ForwardOrdering},DataStructures.PriorityQueue{Box{2,Float64,Float64,3,Float64,Interval{Float64}},Float64,Base.Order.ForwardOrdering},getfield(OptimizeIVQ, Symbol("#addqueues_iv#4")){typeof(f),getfield(OptimizeIVQ, Symbol("#addqueues#3")){DataStructures.PriorityQueue{Box{2,Float64,Float64,3,Float64,Interval{Float64}},Float64,Base.Order.ForwardOrdering},DataStructures.PriorityQueue{Box{2,Float64,Float64,3,Float64,Interval{Float64}},Float64,Base.Order.ForwardOrdering}}},getfield(CoordinateSplittingPTrees, Symbol("#updater#54")){typeof(value),Int64,Int64,DataStructures.PriorityQueue{Box{2,Float64,Float64,3,Float64,Interval{Float64}},Float64,Base.Order.ForwardOrdering},getfield(CoordinateSplittingPTrees, Symbol("#maybe_enqueue#53")){DataStructures.PriorityQueue{Box{2,Float64,Float64,3,Float64,Interval{Float64}},Float64,Base.Order.ForwardOrdering}},CoordinateSplittingPTrees.SymmetricArray{Float64,2,Array{Float64,2}},Array{Float64,1},Array{Float64,1},Array{Bool,1},Array{Float64,2},Array{Float64,2},Array{Float64,2},Array{Float64,1}},Parametron.Model{Float64,Gurobi.Optimizer},Array{Parametron.Functions.Variable,1},Array{Parametron.Functions.Variable,1},Array{Float64,2},Array{Float64,2},Array{Float64,1}})(::Float64) at /home/tim/.julia/dev/OptimizeIVQ/src/OptimizeIVQ.jl:78
[15] #minimize#1(::Int64, ::Float64, ::Int64, ::Function, ::typeof(f), ::Box{2,Float64,Float64,3,Float64,Interval{Float64}}) at /home/tim/.julia/dev/OptimizeIVQ/src/OptimizeIVQ.jl:111
[16] (::getfield(OptimizeIVQ, Symbol("#kw##minimize")))(::NamedTuple{(:nfit, :showprogress),Tuple{Int64,Int64}}, ::typeof(minimize), ::Function, ::Box{2,Float64,Float64,3,Float64,Interval{Float64}}) at ./none:0
[17] top-level scope at none:0
[18] include at ./boot.jl:317 [inlined]
[19] include_relative(::Module, ::String) at ./loading.jl:1041
[20] include at ./sysimg.jl:29 [inlined]
[21] includet(::String) at /home/tim/.julia/dev/Revise/src/Revise.jl:532
[22] top-level scope at none:0
in expression starting at /home/tim/data/interval_registration/ireg2.jl:104
I added some debugging output like this:
tim@diva:~/.julia/dev/Gurobi$ git diff
diff --git a/src/grb_attrs.jl b/src/grb_attrs.jl
index 2a2a591..69712be 100644
--- a/src/grb_attrs.jl
+++ b/src/grb_attrs.jl
@@ -304,6 +304,7 @@ function set_dblattrarray!(model::Model, name::String, start::Integer, len::Inte
ret = @grb_ccall(setdblattrarray, Cint,
(Ptr{Cvoid}, Ptr{UInt8}, Cint, Cint, Ptr{Float64}), model, name, start-1, len, fvec(values))
if ret != 0
+ @show name values
throw(GurobiError(model.env, ret))
end
nothing
diff --git a/src/grb_constrs.jl b/src/grb_constrs.jl
index c78febe..0bc79dd 100644
--- a/src/grb_constrs.jl
+++ b/src/grb_constrs.jl
@@ -60,6 +60,7 @@ function add_constrs!(model::Model, cbegins::IVec, inds::IVec, coeffs::FVec, rel
rel, rhs, C_NULL)
if ret != 0
+ @show cbegins inds coeffs rel rhs
throw(GurobiError(model.env, ret))
end
end
(and I'm pretty sure I've seen the same kind of error with GLPK), with the following result:
cbegins = Int32[1, 7, 13, 19, 25, 31, 37, 43, 49, 55, 61, 67, 73, 79, 85]
inds = Int32[1, 2, 3, 4, 5, 6, 1, 2, 3, 4, 5, 6, 1, 2, 3, 4, 5, 6, 1, 2, 3, 4, 5, 6, 1, 2, 3, 4, 5, 6, 1, 2, 3, 4, 5, 6, 1, 2, 3, 4, 5, 6, 1, 2, 3, 4, 5, 6, 1, 2, 3, 4, 5, 6, 1, 2, 3, 4, 5, 6, 1, 2, 3, 4, 5, 6, 1, 2, 3, 4, 5, 6, 1, 2, 3, 4, 5, 6, 1, 2, 3, 4, 5, 6, 1, 2, 3, 4, 5, 6]
coeffs = [6.91268e-310, 6.91268e-310, 6.91268e-310, 6.91266e-310, 6.91266e-310, 6.91266e-310, 6.91268e-310, 6.91268e-310, 6.91268e-310, 6.91266e-310, 6.91266e-310, 6.91266e-310, 6.91268e-310, 6.91268e-310, 6.91268e-310, 6.91266e-310, 6.91266e-310, 6.91266e-310, 6.91268e-310, 6.91268e-310, 6.91268e-310, 6.91266e-310, 6.91266e-310, 6.91266e-310, 6.91268e-310, 6.91268e-310, 6.91268e-310, 6.91266e-310, 6.91266e-310, 6.91266e-310, 6.91268e-310, 6.91268e-310, 6.91268e-310, 6.91266e-310, 6.91266e-310, 6.91266e-310, 6.91268e-310, 6.91268e-310, 6.91268e-310, 6.91266e-310, 6.91266e-310, 6.91266e-310, 6.91268e-310, 6.91268e-310, 6.91268e-310, 6.91266e-310, 6.91266e-310, 6.91266e-310, 6.91268e-310, 6.91268e-310, 6.91268e-310, 6.91266e-310, 6.91266e-310, 6.91266e-310, 6.91268e-310, 6.91268e-310, 6.91268e-310, 6.91266e-310, 6.91266e-310, 6.91266e-310, 6.91268e-310, 6.91268e-310, 6.91268e-310, 6.91266e-310, 6.91266e-310, 6.91266e-310, 6.91268e-310, 6.91268e-310, 6.91268e-310, 6.91266e-310, 6.91266e-310, 6.91266e-310, 6.91268e-310, 6.91268e-310, 6.91268e-310, 6.91266e-310, 6.91266e-310, 6.91268e-310, 6.91268e-310, 6.91268e-310, 6.91268e-310, 6.91266e-310, 6.91266e-310, 6.91266e-310, 6.91268e-310, 6.91268e-310, 6.91268e-310, 6.91266e-310, 6.91266e-310, 6.91266e-310]
rel = Int8[62, 62, 62, 62, 62, 62, 62, 62, 62, 62, 62, 62, 62, 62, 62]
rhs = [2.06668e-320, NaN, 9.88131e-324, 6.91268e-310, 6.91268e-310, 6.91266e-310, 6.91266e-310, 1.4822e-323, 9.88131e-324, 6.91268e-310, 6.91268e-310, -0.0, 6.91268e-310, 6.91266e-310, 6.91266e-310]
Looks like the C code is getting junk passed to it, perhaps just for the Float64
values.
I am using a slightly unconventional parameter-updating strategy (see #67 (comment)), but I've verified via a try/catch
block that when this error is thrown, the values in the Julia arrays that "back" the Parameter
s are perfectly normal-seeming, e.g.,
X = [0.0 0.0 0.0; 0.0 4.93827 0.0; 0.0 4.93827 -0.148148; -14.8148 4.93827 0.0; 0.0 -14.8148 0.0; 0.0 -14.8148 0.0; -14.8148 4.93827 -0.148148; 0.0 -14.8148 0.444444; 0.0 -14.8148 0.444444; 0.0 -14.8148 0.444444; -14.8148 -14.8148 0.444444; 0.0 0.0 0.444444; 0.0 -29.6296 0.444444; 0.0 4.93827 0.444444; -14.8148 -29.6296 0.444444]
X2 = [0.0 0.0 0.0; 0.0 12.1933 0.0; 0.0 12.1933 0.0109739; 109.739 12.1933 0.0; 0.0 109.739 0.0; 0.0 109.739 0.0; 109.739 12.1933 0.0109739; 0.0 109.739 0.0987654; 0.0 109.739 0.0987654; 0.0 109.739 0.0987654; 109.739 109.739 0.0987654; 0.0 0.0 0.0987654; 0.0 438.957 0.0987654; 0.0 12.1933 0.0987654; 109.739 438.957 0.0987654]
Δf = [0.0, 0.00402906, 0.0354306, 0.0369546, 0.0384715, 0.0384715, 0.0683561, 0.0830869, 0.0830869, 0.0830869, 0.0350268, 0.0943551, 0.0772324, 0.0983842, 0.0291723]
I tried running the example in the README. The Julia REPL crashed on the line
@objective(model, Minimize, residual ⋅ residual)
where the objective is set. Any fixes for this?
It looks like we're not currently handling :tuple
expressions:
julia> using Parametron
julia> @expression(reshape(A, (2, 2)))
ERROR: ArgumentError: Unhandled expression head: tuple
Original expression
Expr
head: Symbol call
args: Array{Any}((3,))
1: Symbol reshape
2: Symbol A
3: Expr
head: Symbol tuple
args: Array{Any}((2,))
1: Int64 2
2: Int64 2
Preprocessed expression
Expr
head: Symbol call
args: Array{Any}((3,))
1: Symbol reshape
2: Symbol A
3: Expr
head: Symbol tuple
args: Array{Any}((2,))
1: Int64 2
2: Int64 2
Current expression:
Expr
head: Symbol tuple
args: Array{Any}((2,))
1: Int64 2
2: Int64 2
We can probably play the same trick we do with vcat
, etc. and just preprocess :tuple
expressions into a call to Base.tuple(x.args...)
.
Currently, constrained optimization doesn't work with the upcoming Gurobi MOI interface. For example:
optimizer = Gurobi.GurobiOptimizer()
model = Model(optimizer)
x = Variable(model)
lhs = [x]
rhs = [0.]
@constraint model lhs <= rhs
solve!(model)
gives:
ArgumentError: ModelLike of type Gurobi.GurobiOptimizer does not support setting the attribute `ConstraintFunction`
Stacktrace:
[1] hidden_set!(::Gurobi.GurobiOptimizer, ::MathOptInterface.ConstraintFunction, ::MathOptInterface.ConstraintIndex{MathOptInterface.VectorAffineFunction{Float64},MathOptInterface.Nonpositives}, ::MathOptInterface.VectorAffineFunction{Float64}) at /home/rdeits/locomotion/explorations/simple-qp-miqp/packages/v0.6/MathOptInterface/src/attributes.jl:588
[2] update!(::SimpleQP.Constraint{Float64,MathOptInterface.Nonpositives}, ::SimpleQP.Model{Float64,Gurobi.GurobiOptimizer}) at /home/rdeits/locomotion/explorations/simple-qp-miqp/packages/v0.6/SimpleQP/src/model.jl:125
[3] update!(::SimpleQP.Model{Float64,Gurobi.GurobiOptimizer}) at /home/rdeits/locomotion/explorations/simple-qp-miqp/packages/v0.6/SimpleQP/src/model.jl:136
[4] solve!(::SimpleQP.Model{Float64,Gurobi.GurobiOptimizer}) at /home/rdeits/locomotion/explorations/simple-qp-miqp/packages/v0.6/SimpleQP/src/model.jl:149
(using SimpleQP master, Gurobi.jl with jump-dev/Gurobi.jl#125 , MOI master, and LQOI with JuliaOpt/LinQuadOptInterface.jl#29 ).
It looks like LInQuadOptInterface defines the available canset
/canmodify
methods, and it only supports things like scalar variable change rather than replacing the constraint function. Gurobi seems to have the relevant low-level APIs to change variable coefficients: https://www.gurobi.com/documentation/8.0/refman/c_grbchgcoeffs.html so I think this should at least be possible. For example, we could probably implement additional canset
/set!
methods in Gurobi.jl to supplement the ones implemented in LQOI. Or maybe we could update LQOI to allow setting an affine constraint function by using the available modify!
methods.
I'm trying to integrate the contact model from MomentumBasedControl.jl into my controller, and since I'm explicitly reasoning about the u
s, i need to extract the torques from my contact wrenches. That means doing something like @expression torque(J, wrench)
where J
is a Paramter{GeometricJacobian}
and wrench
is a Wrench{Variable}
. However, torque(J(), wrench)
is currently type-unstable because promote_type(Float64, Variable)
is Any
. Here's a non-parameterized example:
body = findbody(mechanism, "foot")
path_to_root = path(state.mechanism, body, root_body(state.mechanism))
J = geometric_jacobian(state, path_to_root)
w = Wrench(root_frame(mechanism),
SVector(Variable(1), Variable(2), Variable(3)),
SVector(Variable(1), Variable(2), Variable(3)))
torque(J, w)
which returns:
2-element Array{Any,1}:
-0.0 * x1 + -0.0 * x2 + -0.0 * x3 + -0.0 * x1 + -0.0 * x2 + -1.0 * x3 + 0.0
-0.0 * x1 + -0.0 * x2 + -0.0 * x3 + -0.0 * x1 + -0.0 * x2 + 1.0 * x3 + 0.0
Hi,
I want to use parametron but when I test it for a simple problem, my parameters are not updated when I resolve it.
using Gurobi
using Parametron
using Random, LinearAlgebra
const model=Model(Gurobi.Optimizer(OutputFlag=0));
n = 8; m = 2
x = [Variable(model) for _ = 1 : n]
A_val=ones(n,n);
A= Parameter(model, val=A_val);
b_val=zeros(n);
b= Parameter(model, val=b_val);
C_val=ones(m,n);
C= Parameter(model, val=C_val);
d_val=zeros(m);
d= Parameter(model, val=d_val);
residual = @expression A * x - b
@objective(model, Minimize, residual ⋅ b)
@constraint(model, C * x == d)
@time solve!(model) # curren
println("x: ", value.(model, x));
objvalue=objectivevalue(model);
println("Objective Value: ", objvalue, "(",terminationstatus(model),")");
A_val=randn(n,n);
A= Parameter(model, val=A_val);
b_val=randn(n)+randn(n);
b= Parameter(model, val=b_val);
C_val=randn(m,n)+randn(m,n);
C= Parameter(model, val=C_val);
d_val=rand(m)+rand(m);
d= Parameter(model, val=d_val);
@time solve!(model)
println("x: ", value.(model, x));
objvalue=objectivevalue(model);
println("Objective Value: ", objvalue, "(",terminationstatus(model),")");
Here is the result
5.465271 seconds (16.84 M allocations: 836.260 MiB, 8.30% gc time)
x: [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
Objective Value: 0.0(OPTIMAL)
0.000167 seconds (77 allocations: 8.906 KiB)
x: [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
Objective Value: 0.0(OPTIMAL)
We already have an optimization for transpose(x) * A * x
expressions, but it doesn't match x' * A * x
.
(just continuing our conversation from the QPControl bug)
What about a function defer(x)
which would do:
function defer(x, model)
Parameter(() -> x, model)
end
which would be useful for wrapping mutable objects that should be evaluated "later" rather than eagerly. I realize it's only a one-liner, but it would also remove the need for a let
block, so:
let x = x
Parameter(() -> x, model)
end
becomes:
defer(x, model)
We're not doing anything platform-specific, but it might be nice to test this on a 32-bit system just in case we accidentally assume Int === Int64
and break some of the optimizations. OSQP.jl has an appveyor build, so it should work OK.
Accessing fields inside an @expression
causes allocation, as we end up transforming x.y
into getfield(x, :y)
, which can't be inferred (at least in 0.6). For example:
frame = CartesianFrame3D()
model = SimpleQP.MockModel()
p = Parameter(model) do
Point3D(frame, 0., 0., 0.)
end
expr = @expression p.v
SimpleQP.findallocs(expr)
gives:
LazyExpression{FunctionWrapper{…}(LazyExpression{Core.#getfield, …}(…))}(…): 112 bytes
[1]: Parameter{RigidBodyDynamics.Spatial.Point3D{StaticArrays.SArray{Tuple{3},Float64,1,3}}, …}(…): 48 bytes
[2]: Symbol
We probably need to special-case Core.getfield
.
I just noticed that in @expression
and in @constraint
sparse matrices are not accepted.
For example, A=spzeros(n,m)
. If I use
A_expr = @expression A
or @constraint(model, A * X ==0)
It gives me this error:
ERROR: LoadError: MethodError: no method matching one(::Type{Any})
Closest candidates are:
one(::Type{Union{Missing, T}}) where T at missing.jl:87
one(::Missing) at missing.jl:83
one(::BitArray{2}) at bitarray.jl:400
...
Stacktrace:
[1] one(::Type{Any}) at ./missing.jl:87
[2] *(::SparseMatrixCSC{Float32,Int64}, ::Array{Variable,1})at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.1/SparseArrays/src/linalg.jl:53
Originally posted by @ssadat in #107 (comment)
A declarative, efficient, and flexible JavaScript library for building user interfaces.
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google ❤️ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.