diff --git a/src/PiecewiseLinearOpt.jl b/src/PiecewiseLinearOpt.jl index 6646c48..4327bec 100644 --- a/src/PiecewiseLinearOpt.jl +++ b/src/PiecewiseLinearOpt.jl @@ -6,10 +6,10 @@ module PiecewiseLinearOpt using JuMP -import MathOptInterface -const MOI = MathOptInterface -using LinearAlgebra -using Random + +import LinearAlgebra +import MathOptInterface as MOI +import Random export PWLFunction, UnivariatePWLFunction, BivariatePWLFunction, piecewiselinear diff --git a/src/jump.jl b/src/jump.jl index cd79b05..431c978 100644 --- a/src/jump.jl +++ b/src/jump.jl @@ -220,8 +220,8 @@ function sos2_encoding_constraints!(m, λ, y, h, B) n = length(λ)-1 for b in B JuMP.@constraints(m, begin - dot(b,h[1])*λ[1] + sum(min(dot(b,h[v]),dot(b,h[v-1]))*λ[v] for v in 2:n) + dot(b,h[n])*λ[n+1] ≤ dot(b,y) - dot(b,h[1])*λ[1] + sum(max(dot(b,h[v]),dot(b,h[v-1]))*λ[v] for v in 2:n) + dot(b,h[n])*λ[n+1] ≥ dot(b,y) + LinearAlgebra.dot(b,h[1])*λ[1] + sum(min(LinearAlgebra.dot(b,h[v]),LinearAlgebra.dot(b,h[v-1]))*λ[v] for v in 2:n) + LinearAlgebra.dot(b,h[n])*λ[n+1] ≤ LinearAlgebra.dot(b,y) + LinearAlgebra.dot(b,h[1])*λ[1] + sum(max(LinearAlgebra.dot(b,h[v]),LinearAlgebra.dot(b,h[v-1]))*λ[v] for v in 2:n) + LinearAlgebra.dot(b,h[n])*λ[n+1] ≥ LinearAlgebra.dot(b,y) end) end return nothing @@ -365,9 +365,9 @@ function compute_hyperplanes(C::Vector{Vector{T}}) where T <: Number if indices == [n-1] break end - if rank(d[indices,:]) == length(indices) && length(indices) <= k-1 + if LinearAlgebra.rank(d[indices,:]) == length(indices) && length(indices) <= k-1 if length(indices) == k-1 - nullsp = nullspace(d[indices,:]) + nullsp = LinearAlgebra.nullspace(d[indices,:]) @assert size(nullsp,2) == 1 v = vec(nullsp) push!(spanners, canonical!(v)) @@ -406,7 +406,7 @@ function compute_hyperplanes(C::Vector{Vector{T}}) where T <: Number end function canonical!(v::Vector{Float64}) - normalize!(v) + LinearAlgebra.normalize!(v) for j in 1:length(v) if abs(v[j]) < 1e-8 v[j] = 0 @@ -582,7 +582,7 @@ function piecewiselinear(m::JuMP.Model, x₁::VarOrAff, x₂::VarOrAff, pwl::Biv A = [p¹[1] p¹[2] 1 p²[1] p²[2] 1 p³[1] p³[2] 1] - @assert rank(A) == 3 + @assert LinearAlgebra.rank(A) == 3 b = [0, 0, 1] q = A \ b @assert isapprox(q[1]*p¹[1] + q[2]*p¹[2] + q[3], 0, atol=1e-4) diff --git a/src/types.jl b/src/types.jl index 5f86985..f516220 100644 --- a/src/types.jl +++ b/src/types.jl @@ -44,7 +44,7 @@ function BivariatePWLFunction(x, y, fz::Function; pattern=:BestFit, seed=hash((l m = length(x) n = length(y) - mt = MersenneTwister(seed) + mt = Random.MersenneTwister(seed) # run for each square on [x[i],x[i+1]] × [y[i],y[i+1]] for i in 1:length(x)-1, j in 1:length(y)-1 SWt, NWt, NEt, SEt = LinearIndices((m,n))[i,j], LinearIndices((m,n))[i,j+1], LinearIndices((m,n))[i+1,j+1], LinearIndices((m,n))[i+1,j] diff --git a/test/runtests.jl b/test/runtests.jl index 8d2b81a..cbfe74d 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -3,22 +3,51 @@ # Use of this source code is governed by an MIT-style license that can be found # in the LICENSE.md file or at https://opensource.org/licenses/MIT. -using Cbc +using PiecewiseLinearOpt using Test -using LinearAlgebra +import Cbc import JuMP -import MathOptInterface -const MOI = MathOptInterface - -using PiecewiseLinearOpt - -# TODO: Add :SOS2 to this list, once Cbc supports it -methods_1D = (:CC,:MC,:Logarithmic,:LogarithmicIB,:ZigZag,:ZigZagInteger,:GeneralizedCelaya,:SymmetricCelaya,:Incremental,:DisaggLogarithmic) -# TODO: Add :SOS2 to this list, once Cbc supports it -# TODO: Add :MC to this list, Cbc (but not Gurobi) gives a different answer below, only for :MC (maybe a bug in Cbc?) -methods_2D = (:CC,:Logarithmic,:LogarithmicIB,:ZigZag,:ZigZagInteger,:GeneralizedCelaya,:SymmetricCelaya,:DisaggLogarithmic) -patterns_2D = (:Upper,:Lower,:BestFit,:UnionJack,:K1,:Random) # :OptimalTriangleSelection and :Stencil not supported currently +import LinearAlgebra +import MathOptInterface as MOI + +methods_1D = ( + :CC, + :MC, + :Logarithmic, + :LogarithmicIB, + :ZigZag, + :ZigZagInteger, + :GeneralizedCelaya, + :SymmetricCelaya, + :Incremental, + :DisaggLogarithmic, + # :SOS2, not supported by Cbc +) + +methods_2D = ( + :CC, + :Logarithmic, + :LogarithmicIB, + :ZigZag, + :ZigZagInteger, + :GeneralizedCelaya, + :SymmetricCelaya, + :DisaggLogarithmic, + # :SOS2, not supported by Cbc + # TODO: Add :MC to this list, Cbc (but not Gurobi) gives a different answer + # below, only for :MC (maybe a bug in Cbc?) +) +patterns_2D = ( + :Upper, + :Lower, + :BestFit, + :UnionJack, + :K1, + :Random, + # :OptimalTriangleSelection not supported currently + # :Stencil +) optimizer=JuMP.optimizer_with_attributes(Cbc.Optimizer, MOI.Silent() => true) @@ -28,17 +57,12 @@ optimizer=JuMP.optimizer_with_attributes(Cbc.Optimizer, MOI.Silent() => true) JuMP.@variable(model, x) z = piecewiselinear(model, x, range(1,stop=2π, length=8), sin, method=method) JuMP.@objective(model, Max, z) - JuMP.optimize!(model) - @test JuMP.termination_status(model) == MOI.OPTIMAL @test JuMP.value(x) ≈ 1.75474 rtol=1e-4 @test JuMP.value(z) ≈ 0.98313 rtol=1e-4 - JuMP.@constraint(model, x ≤ 1.5z) - JuMP.optimize!(model) - @test JuMP.termination_status(model) == MOI.OPTIMAL @test JuMP.value(x) ≈ 1.36495 rtol=1e-4 @test JuMP.value(z) ≈ 0.90997 rtol=1e-4 @@ -55,20 +79,15 @@ end f = (x1,x2) -> 2*(x1-1/3)^2 + 3*(x2-4/7)^4 z = piecewiselinear(model, x[1], x[2], BivariatePWLFunction(d, d, f, pattern=pattern), method=method) JuMP.@objective(model, Min, z) - JuMP.optimize!(model) - @test JuMP.termination_status(model) == MOI.OPTIMAL @test JuMP.value(x[1]) ≈ 0.285714 rtol=1e-4 @test JuMP.value(x[2]) ≈ 0.571429 rtol=1e-4 @test JuMP.value(z) ≈ 0.004535 rtol=1e-3 @test JuMP.objective_value(model) ≈ 0.004535 rtol=1e-3 @test JuMP.objective_value(model) ≈ JuMP.value(z) rtol=1e-3 - JuMP.@constraint(model, x[1] ≥ 0.6) - JuMP.optimize!(model) - @test JuMP.termination_status(model) == MOI.OPTIMAL @test JuMP.value(x[1]) ≈ 0.6 rtol=1e-4 @test JuMP.value(x[2]) ≈ 0.571428 rtol=1e-4 @@ -77,7 +96,7 @@ end @test JuMP.objective_value(model) ≈ JuMP.value(z) rtol=1e-3 end end -# + # println("\nbivariate optimal IB scheme tests") # @testset "2D: optimal IB, UnionJack" begin # model = JuMP.Model(JuMP.with_optimizer(Cbc.Optimizer)) @@ -87,20 +106,15 @@ end # f = (x,y) -> 2*(x-1/3)^2 + 3*(y-4/7)^4 # z = piecewiselinear(model, x, y, BivariatePWLFunction(d, d, f, pattern=:UnionJack), method=:OptimalIB, subsolver=solver) # JuMP.@objective(model, Min, z) -# # JuMP.optimize!(model) -# # @test JuMP.termination_status(model) == MOI.OPTIMAL # @test JuMP.value(x) ≈ 0.5 rtol=1e-4 # @test JuMP.value(y) ≈ 0.5 rtol=1e-4 # @test JuMP.value(z) ≈ 0.055634 rtol=1e-3 # @test getobjectivevalue(model) ≈ 0.055634 rtol=1e-3 # @test getobjectivevalue(model) ≈ JuMP.value(z) rtol=1e-3 -# # JuMP.@constraint(model, x ≥ 0.6) -# # JuMP.optimize!(model) -# # @test JuMP.termination_status(model) == MOI.OPTIMAL # @test JuMP.value(x) ≈ 0.6 rtol=1e-4 # @test JuMP.value(y) ≈ 0.5 rtol=1e-4