diff --git a/bench/runbench.jl b/bench/runbench.jl new file mode 100644 index 0000000..3d0f0d0 --- /dev/null +++ b/bench/runbench.jl @@ -0,0 +1,148 @@ +using SparseArrays +using LinearAlgebra +using Random +using Clp +using MathOptInterface +const MOI = MathOptInterface + +using TimerOutputs + +struct RandomLP + rows::Int + cols::Int + dens::Float64 +end + +function generate_moi_problem(model, At, b, c; + var_bounds = false, scalar = true) + cols, rows = size(At) + x = MOI.add_variables(model, cols) + A_cols = rowvals(At) + A_vals = nonzeros(At) + if var_bounds + for col in 1:cols + MOI.add_constraint(model, MOI.SingleVariable(x[col]), + MOI.LessThan(10.0)) + MOI.add_constraint(model, MOI.SingleVariable(x[col]), + MOI.GreaterThan(-10.0)) + end + end + if scalar + for row in 1:rows + MOI.add_constraint(model, MOI.ScalarAffineFunction( + [MOI.ScalarAffineTerm(A_vals[i], x[A_cols[i]]) for i in nzrange(At, row)], 0.0), + MOI.LessThan(b[row])) + end + else + for row in 1:rows + MOI.add_constraint(model, MOI.VectorAffineFunction( + [MOI.VectorAffineTerm(1, + MOI.ScalarAffineTerm(A_vals[i], x[A_cols[i]]) + ) for i in nzrange(At, row)], [-b[row]]), + MOI.Nonpositives(1)) + end + end + objective = MOI.ScalarAffineFunction( + [MOI.ScalarAffineTerm(c[i], x[i]) for i in findall(!iszero, c)], + 0.0) + MOI.set(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), objective) + MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE) + return x +end + +function random_data(seed, data) + + rows = data.rows + cols = data.cols + density = data.dens + + p_neg_element = 0.0 # 0.25 + + rng = Random.MersenneTwister(seed) + + f_A(r, n) = ifelse.(rand(r, n) .> p_neg_element, 1, -1) .* (15 .+ 30 .* rand(r, n)) + + At = sprand(rng, cols, rows, density, f_A) + b = 50 * rand(rng, rows) + + # not using signs now + sign = ifelse.(rand(rng, rows) .> 0.2, 'L', 'G') + + f_c(r, n) = 20 .* 2 .* (rand(r, n) .- 0.5) + c = sprand(rng, cols, 0.5, f_c) + + return At, b, c +end + +function bridged_cache_and_solver() + model = MOI.Bridges.full_bridge_optimizer(MOI.Utilities.CachingOptimizer( + MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()), + MOI.Utilities.MANUAL), Float64) + clp = Clp.Optimizer() + MOI.set(clp, MOI.Silent(), true) + return model, clp +end +function cache_and_solver() + model = MOI.Utilities.CachingOptimizer( + MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()), + MOI.Utilities.MANUAL) + clp = Clp.Optimizer() + MOI.set(clp, MOI.Silent(), true) + return model, clp +end +function bridged_cached_solver() + model = MOI.Bridges.full_bridge_optimizer(MOI.Utilities.CachingOptimizer( + MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()), + Clp.Optimizer()), Float64) + MOI.set(model, MOI.Silent(), true) + return model +end +function cached_solver() + model = MOI.Utilities.CachingOptimizer( + MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()), + Clp.Optimizer()) + MOI.set(model, MOI.Silent(), true) + return model +end + +function time_build_and_solve(to_build, to_solve, At, b, c, scalar = true) + @timeit "build" x = generate_moi_problem(to_build, At, b, c, scalar = scalar) + if to_build !== to_solve + @timeit "copy" MOI.copy_to(to_solve, to_build, copy_names = false) + end + @time @timeit "opt" MOI.optimize!(to_solve) + MOI.get(to_solve, MOI.ObjectiveValue()) + val = MOI.get(to_solve, MOI.SolveTime()) + println(val) +end + +function solve_clp(seed, data; time_limit_sec=Inf) + + reset_timer!() + + for i in 1:seed + + @timeit "data" At, b, c = random_data(i, data) + + bridged_cache, pure_solver = bridged_cache_and_solver() + @timeit "bc + s" time_build_and_solve(bridged_cache, pure_solver, At, b, c) + + cache, pure_solver2 = cache_and_solver() + @timeit "c + s" time_build_and_solve(cache, pure_solver2, At, b, c) + + full_solver = bridged_cached_solver() + @timeit "bcs" time_build_and_solve(full_solver, full_solver, At, b, c) + + full_solver = bridged_cached_solver() + @timeit "bcs + v" time_build_and_solve(full_solver, full_solver, At, b, c, false) + + cache_solver = cached_solver() + @timeit "cs" time_build_and_solve(cache_solver, cache_solver, At, b, c) + + end + + print_timer() + +end + +solve_clp(10, RandomLP(10000, 20000, 0.01); time_limit_sec=5) \ No newline at end of file diff --git a/src/MOI_wrapper/MOI_wrapper.jl b/src/MOI_wrapper/MOI_wrapper.jl index d229f84..9a08726 100644 --- a/src/MOI_wrapper/MOI_wrapper.jl +++ b/src/MOI_wrapper/MOI_wrapper.jl @@ -185,7 +185,7 @@ _add_bounds(lb, ::Vector{Float64}, i, s::MOI.GreaterThan{Float64}) = lb[i] = s.l _add_bounds(lb, ub, i, s::MOI.EqualTo{Float64}) = lb[i], ub[i] = s.value, s.value _add_bounds(lb, ub, i, s::MOI.Interval{Float64}) = lb[i], ub[i] = s.lower, s.upper -function _extract_bound_data(src, mapping, lb, ub, S) +function _extract_bound_data(src, mapping, lb, ub, ::Type{S}) where S for con_index in MOI.get( src, MOI.ListOfConstraintIndices{MOI.SingleVariable, S}() ) @@ -220,15 +220,33 @@ _bounds(s::MOI.LessThan{Float64}) = (-Inf, s.upper) _bounds(s::MOI.EqualTo{Float64}) = (s.value, s.value) _bounds(s::MOI.Interval{Float64}) = (s.lower, s.upper) -function _extract_row_data(src, mapping, lb, ub, I, J, V, S) +function add_sizehint!(vec, n) + len = length(vec) + return sizehint!(vec, len + n) +end + +function _extract_row_data(src, mapping, lb, ub, I, J, V, ::Type{S}) where S row = length(I) == 0 ? 1 : I[end] + 1 - for c_index in MOI.get( + list = MOI.get( src, MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{Float64}, S}() ) + add_sizehint!(lb, length(list)) + add_sizehint!(ub, length(list)) + n_terms = 0 + fs = Array{MOI.ScalarAffineFunction{Float64}}(undef, length(list)) + for (i,c_index) in enumerate(list) f = MOI.get(src, MOI.ConstraintFunction(), c_index) + fs[i] = f l, u = _bounds(MOI.get(src, MOI.ConstraintSet(), c_index)) push!(lb, l - f.constant) push!(ub, u - f.constant) + n_terms += length(f.terms) + end + add_sizehint!(I, n_terms) + add_sizehint!(J, n_terms) + add_sizehint!(V, n_terms) + for (i,c_index) in enumerate(list) + f = fs[i]#MOI.get(src, MOI.ConstraintFunction(), c_index) for term in f.terms push!(I, row) push!(J, Cint(mapping.varmap[term.variable_index].value)) @@ -242,12 +260,7 @@ function _extract_row_data(src, mapping, lb, ub, I, J, V, S) return end -function MOI.copy_to( - dest::Optimizer, - src::MOI.ModelLike; - copy_names::Bool = false -) - @assert MOI.is_empty(dest) +function test_data(src, dest) for (F, S) in MOI.get(src, MOI.ListOfConstraints()) if !MOI.supports_constraint(dest, F, S) throw(MOI.UnsupportedConstraint{F, S}("Clp.Optimizer does not support constraints of type $F-in-$S.")) @@ -257,19 +270,30 @@ function MOI.copy_to( if !MOI.supports(dest, MOI.ObjectiveFunction{fobj_type}()) throw(MOI.UnsupportedAttribute(MOI.ObjectiveFunction(fobj_type))) end +end + +function MOI.copy_to( + dest::Optimizer, + src::MOI.ModelLike; + copy_names::Bool = false +) + @assert MOI.is_empty(dest) + test_data(src, dest) + mapping = MOI.Utilities.IndexMap() N, c = _copy_to_columns(dest, src, mapping) cl, cu = fill(-Inf, N), fill(Inf, N) rl, ru, I, J, V = Float64[], Float64[], Cint[], Cint[], Float64[] - for S in ( - MOI.GreaterThan{Float64}, - MOI.LessThan{Float64}, - MOI.EqualTo{Float64}, - MOI.Interval{Float64}, - ) - _extract_bound_data(src, mapping, cl, cu, S) - _extract_row_data(src, mapping, rl, ru, I, J, V, S) - end + + _extract_bound_data(src, mapping, cl, cu, MOI.GreaterThan{Float64}) + _extract_row_data(src, mapping, rl, ru, I, J, V, MOI.GreaterThan{Float64}) + _extract_bound_data(src, mapping, cl, cu, MOI.LessThan{Float64}) + _extract_row_data(src, mapping, rl, ru, I, J, V, MOI.LessThan{Float64}) + _extract_bound_data(src, mapping, cl, cu, MOI.EqualTo{Float64}) + _extract_row_data(src, mapping, rl, ru, I, J, V, MOI.EqualTo{Float64}) + _extract_bound_data(src, mapping, cl, cu, MOI.Interval{Float64}) + _extract_row_data(src, mapping, rl, ru, I, J, V, MOI.Interval{Float64}) + M = Cint(length(rl)) A = SparseArrays.sparse(I, J, V, M, N) Clp_loadProblem(