diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 3743a4b..65db331 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -15,7 +15,7 @@ jobs: strategy: fail-fast: false matrix: - version: ['1.6', '1'] + version: ['1.10', '1'] os: [ubuntu-latest, macOS-latest, windows-latest] arch: [x64] include: diff --git a/.gitignore b/.gitignore index b1fb063..9071367 100644 --- a/.gitignore +++ b/.gitignore @@ -3,3 +3,5 @@ *.jl.mem Manifest.toml *.DS_Store +.DS_Store +.vscode/settings.json diff --git a/Project.toml b/Project.toml index 50d9b8c..07c3da3 100644 --- a/Project.toml +++ b/Project.toml @@ -6,18 +6,16 @@ version = "0.4.2" [deps] JuMP = "4076af6c-e467-56ae-b986-b466b2749572" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" -MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" [compat] -Cbc = "1" -JuMP = "0.23, 1" -MathOptInterface = "1" -julia = "1.6" +HiGHS = "1" +JuMP = "1" +julia = "1.10" [extras] -Cbc = "9961bab8-2fa3-5c5a-9d89-47fab24efd76" +HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["Cbc", "Test"] +test = ["Test", "HiGHS"] diff --git a/README.md b/README.md index abfab63..f4f6c19 100644 --- a/README.md +++ b/README.md @@ -33,37 +33,21 @@ Pkg.add("PiecewiseLinearOpt") ## Use with JuMP Current support is limited to modeling the graph of a continuous piecewise -linear function, either univariate or bivariate, with the goal of adding support -for the epigraphs of lower semicontinuous piecewise linear functions. +linear function, with a primary focus on univariate or bivariate functions. +There are also methods for more general multivariate problems. ### Univariate -Consider a piecewise linear function `f`. The function is described a domain `d`, -which is a set of breakpoints between pieces, and the function value `fd` at -those breakpoints: +Consider a piecewise linear function described by a domain `d`, +which is a set of breakpoints between pieces, and the function value at +those breakpoints given by the function `f` at those points: ```julia -julia> f(x) = sin(x) -f (generic function with 1 method) - julia> d = 0:0.5:2pi 0.0:0.5:6.0 -julia> fd = f.(d) -13-element Vector{Float64}: - 0.0 - 0.479425538604203 - 0.8414709848078965 - 0.9974949866040544 - 0.9092974268256817 - 0.5984721441039564 - 0.1411200080598672 - -0.35078322768961984 - -0.7568024953079282 - -0.977530117665097 - -0.9589242746631385 - -0.7055403255703919 - -0.27941549819892586 +julia> f(x) = sin(x) +f (generic function with 1 method) ``` To represent this function in a JuMP model, do: @@ -72,13 +56,14 @@ To represent this function in a JuMP model, do: using JuMP, PiecewiseLinearOpt model = Model() @variable(model, x) -z = PiecewiseLinearOpt.piecewiselinear(model, x, d, fd; method = :CC) +z = PiecewiseLinearOpt.piecewiselinear(model, x, d, f; method = Logarithmic()) @objective(model, Min, z) # minimize f(x) ``` ### Bivariate -Consider piecewise linear approximation for the function $f(x, y) = exp(x + y)$: +Consider a piecewise linear approximation for the function $f(x, y) = exp(x + y)$ +on a triangular grid with a best fit pattern: ```julia using JuMP, PiecewiseLinearOpt @@ -92,38 +77,34 @@ z = PiecewiseLinearOpt.piecewiselinear( 0:0.1:1, 0:0.1:1, (u, v) -> exp(u + v); - method = :DisaggLogarithmic, + method = SixStencil(), + pattern = :BestFit ) @objective(model, Min, z) ``` ## Methods +The following formualations are available in the package and is provided through the +`method` argument: + +Supported multivariate formulations: +* `ConvexCombination()` +* `DisaggregatedLogarithmic()` +* `MultipleChoice()`: Limited support as it currently needs an explicit formulations with hyperplanes + Supported univariate formulations: +* `Incremental()` +* `Logarithmic()` +* `LogarithmicIndependentBranching()` +* `NativeSOS2()` +* `ZigZagBinary()` +* `ZigZagInteger()` + +The following bivariate formulations are available and can be combined with most univariate +formulations to impose two axis-aligned SOS2 constraints. See the associated paper for more details. +* `K1(sos2_method)`: requires a K1 grid triangulation +* `UnionJack(sos2_method)`: requires a UnionJack grid triangulation +* `SixStencil(sos2_method)`: requires a grid triangulation +* `NineStencil(sos2_method)`: requires a grid triangulation -* Convex combination (`:CC`) -* Multiple choice (`:MC`) -* Native SOS2 branching (`:SOS2`) -* Incremental (`:Incremental`) -* Logarithmic (`:Logarithmic`; default) -* Disaggregated Logarithmic (`:DisaggLogarithmic`) -* Binary zig-zag (`:ZigZag`) -* General integer zig-zag (`:ZigZagInteger`) - -Supported bivariate formulations for entire constraint: - -* Convex combination (`:CC`) -* Multiple choice (`:MC`) -* Disaggregated Logarithmic (`:DisaggLogarithmic`) - -Also, you can use any univariate formulation for bivariate functions as well. -They will be used to impose two axis-aligned SOS2 constraints, along with the -"6-stencil" formulation for the triangle selection portion of the constraint. -See the associated paper for more details. In particular, the following are also -acceptable bivariate formulation choices: - -* Native SOS2 branching (`:SOS2`) -* Incremental (`:Incremental`) -* Logarithmic (`:Logarithmic`) -* Binary zig-zag (`:ZigZag`) -* General integer zig-zag (`:ZigZagInteger`) diff --git a/src/PiecewiseLinearOpt.jl b/src/PiecewiseLinearOpt.jl index 4327bec..71a5f8d 100644 --- a/src/PiecewiseLinearOpt.jl +++ b/src/PiecewiseLinearOpt.jl @@ -8,12 +8,68 @@ module PiecewiseLinearOpt using JuMP import LinearAlgebra -import MathOptInterface as MOI import Random export PWLFunction, UnivariatePWLFunction, BivariatePWLFunction, piecewiselinear include("types.jl") -include("jump.jl") + +mutable struct PWLData + counter::Int + PWLData() = new(0) +end + +function initPWL!(m::JuMP.Model) + if !haskey(m.ext, :PWL) + m.ext[:PWL] = PWLData() + end + return nothing +end + +const VarOrAff = Union{JuMP.VariableRef,JuMP.AffExpr} + +include("methods/util.jl") + +export Incremental, + LogarithmicEmbedding, + LogarithmicIndependentBranching, + NativeSOS2, + ZigZagBinary, + ZigZagInteger +include("methods/univariate/incremental.jl") + +include("methods/univariate/logarithmic_embedding.jl") +include("methods/univariate/logarithmic_independent_branching.jl") +include("methods/univariate/native_sos2.jl") +include("methods/univariate/zig_zag_binary.jl") +include("methods/univariate/zig_zag_integer.jl") +# ConvexCombination has an SOS2 formulation, so defer this until after the +# multivariate formulations are defined +include("methods/univariate/sos2_formulation_base.jl") + +# Consider the colloqial "log" to refer to the embedding formulation +const Logarithmic = LogarithmicEmbedding +export Logarithmic + +export K1, + NineStencil, + OptimalIndependentBranching, + OptimalTriangleSelection, + SixStencil, + UnionJack +include("methods/bivariate/k1.jl") +include("methods/bivariate/nine_stencil.jl") +include("methods/bivariate/optimal_independent_branching.jl") +include("methods/bivariate/optimal_triangle_selection.jl") +include("methods/bivariate/six_stencil.jl") +include("methods/bivariate/union_jack.jl") +include("methods/bivariate/common.jl") + +export ConvexCombination, DisaggregatedLogarithmic, MultipleChoice +include("methods/multivariate/convex_combination.jl") +include("methods/multivariate/disaggregated_logarithmic.jl") +include("methods/multivariate/multiple_choice.jl") + +include("pwlinear.jl") end # module diff --git a/src/jump.jl b/src/jump.jl deleted file mode 100644 index 70f1393..0000000 --- a/src/jump.jl +++ /dev/null @@ -1,1270 +0,0 @@ -# Copyright (c) 2016: Joey Huchette and contributors -# -# 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. - -# TODO: choose method based on problem size -defaultmethod() = :Logarithmic - -mutable struct PWLData - counter::Int - PWLData() = new(0) -end - -function initPWL!(m::JuMP.Model) - if !haskey(m.ext, :PWL) - m.ext[:PWL] = PWLData() - end - return -end - -const VarOrAff = Union{JuMP.VariableRef,JuMP.AffExpr} - -function piecewiselinear( - m::JuMP.Model, - x::VarOrAff, - d, - f::Function; - method = defaultmethod(), -) - initPWL!(m) - fd = [f(xx) for xx in d] - return piecewiselinear(m, x, d, fd; method = method) -end - -function piecewiselinear( - m::JuMP.Model, - x::VarOrAff, - d, - fd; - method = defaultmethod(), -) - return piecewiselinear(m, x, UnivariatePWLFunction(d, fd); method = method) -end - -function piecewiselinear( - m::JuMP.Model, - x::VarOrAff, - pwl::UnivariatePWLFunction; - method = defaultmethod(), -) - initPWL!(m) - counter = m.ext[:PWL].counter - counter += 1 - m.ext[:PWL].counter = counter - d = [_x[1] for _x in pwl.x] - fd = pwl.z - n = length(d) - - if n != length(fd) - error( - "You provided a different number of breakpoints ($n) and function values at the breakpoints ($(length(fd)))", - ) - end - if n == 0 - error( - "I don't know how to handle a piecewise linear function with no breakpoints", - ) - end - z = JuMP.@variable( - m, - lower_bound = minimum(fd), - upper_bound = maximum(fd), - base_name = "z_$counter" - ) - if n == 1 - JuMP.@constraint(m, x == d[1]) - JuMP.@constraint(m, z == fd[1]) - return z - end - if method == :Incremental - δ = JuMP.@variable( - m, - [1:n], - lower_bound = 0, - upper_bound = 1, - base_name = "δ_$counter" - ) - y = JuMP.@variable(m, [1:n-1], Bin, base_name = "y_$counter") - JuMP.@constraint( - m, - x == d[1] + sum(δ[i] * (d[i+1] - d[i]) for i in 1:n-1) - ) - JuMP.@constraint( - m, - z == fd[1] + sum(δ[i] * (fd[i+1] - fd[i]) for i in 1:n-1) - ) - for i in 1:n-1 - JuMP.@constraint(m, δ[i+1] ≤ y[i]) - JuMP.@constraint(m, y[i] ≤ δ[i]) - end - elseif method == :MC - x̂ = JuMP.@variable(m, [1:n-1], base_name = "x̂_$counter") - ẑ = JuMP.@variable(m, [1:n-1], base_name = "ẑ_$counter") - y = JuMP.@variable(m, [1:n-1], Bin, base_name = "y_$counter") - JuMP.@constraint(m, sum(y) == 1) - JuMP.@constraint(m, sum(x̂) == x) - JuMP.@constraint(m, sum(ẑ) == z) - Δ = [(fd[i+1] - fd[i]) / (d[i+1] - d[i]) for i in 1:n-1] - for i in 1:n-1 - JuMP.@constraints( - m, - begin - x̂[i] ≥ d[i] * y[i] - x̂[i] ≤ d[i+1] * y[i] - ẑ[i] == fd[i] * y[i] + Δ[i] * (x̂[i] - d[i] * y[i]) - end - ) - end - elseif method in (:DisaggLogarithmic, :DLog) - γ = JuMP.@variable( - m, - [i = 1:n, j = max(1, i - 1):min(n - 1, i)], - lower_bound = 0, - upper_bound = 1 - ) - JuMP.@constraint(m, sum(γ) == 1) - JuMP.@constraint( - m, - γ[1, 1] * d[1] + - sum((γ[i, i-1] + γ[i, i]) * d[i] for i in 2:n-1) + - γ[n, n-1] * d[n] == x - ) - JuMP.@constraint( - m, - γ[1, 1] * fd[1] + - sum((γ[i, i-1] + γ[i, i]) * fd[i] for i in 2:n-1) + - γ[n, n-1] * fd[n] == z - ) - r = ceil(Int, log2(n - 1)) - H = reflected_gray_codes(r) - y = JuMP.@variable(m, [1:r], Bin) - for j in 1:r - JuMP.@constraint( - m, - sum((γ[i, i] + γ[i+1, i]) * H[i][j] for i in 1:(n-1)) == y[j] - ) - end - else - # V-formulation methods - λ = JuMP.@variable( - m, - [1:n], - lower_bound = 0, - upper_bound = 1, - base_name = "λ_$counter" - ) - JuMP.@constraint(m, sum(λ) == 1) - JuMP.@constraint(m, sum(λ[i] * d[i] for i in 1:n) == x) - JuMP.@constraint(m, sum(λ[i] * fd[i] for i in 1:n) == z) - if method in (:Logarithmic, :Log) - sos2_logarithmic_formulation!(m, λ) - elseif method in (:LogarithmicIB, :LogIB) - sos2_logarithmic_IB_formulation!(m, λ) - elseif method == :CC - sos2_cc_formulation!(m, λ) - elseif method in (:ZigZag, :ZZB) - sos2_zigzag_formulation!(m, λ) - elseif method in (:ZigZagInteger, :ZZI) - sos2_zigzag_general_integer_formulation!(m, λ) - elseif method == :GeneralizedCelaya - sos2_generalized_celaya_formulation!(m, λ) - elseif method == :SymmetricCelaya - sos2_symmetric_celaya_formulation!(m, λ) - elseif method == :SOS2 - JuMP.@constraint(m, λ in MOI.SOS2([i for i in 1:n])) - else - error("Unrecognized method $method") - end - end - return z -end - -function sos2_cc_formulation!(m::JuMP.Model, λ) - counter = m.ext[:PWL].counter - n = length(λ) - y = JuMP.@variable(m, [1:n-1], Bin, base_name = "y_$counter") - JuMP.@constraint(m, sum(y) == 1) - JuMP.@constraint(m, λ[1] ≤ y[1]) - for i in 2:n-1 - JuMP.@constraint(m, λ[i] ≤ y[i-1] + y[i]) - end - JuMP.@constraint(m, λ[n] ≤ y[n-1]) - return -end - -# function sos2_mc_formulation!(m::JuMP.Model, λ) # not currently used -# counter = m.ext[:PWL].counter -# n = length(λ) -# γ = JuMP.@variable(m, [1:n-1, 1:n], base_name="γ_$counter") -# y = JuMP.@variable(m, [1:n-1], Bin, base_name="y_$counter") -# JuMP.@constraint(m, sum(y) == 1) -# JuMP.@constraint(m, sum(γ[i,:] for i in 1:n-1) .== λ) -# for i in 1:n-1 -# JuMP.@constraint(m, γ[i,i] + γ[i,i+1] ≥ y[i]) -# end -# return -# end - -function sos2_logarithmic_formulation!(m::JuMP.Model, λ) - counter = m.ext[:PWL].counter - n = length(λ) - 1 - k = ceil(Int, log2(n)) - y = JuMP.@variable(m, [1:k], Bin, base_name = "y_$counter") - sos2_encoding_constraints!( - m, - λ, - y, - reflected_gray_codes(k), - unit_vector_hyperplanes(k), - ) - return -end - -function sos2_logarithmic_IB_formulation!(m::JuMP.Model, λ) # IB: independent branching - counter = m.ext[:PWL].counter - n = length(λ) - 1 - k = ceil(Int, log2(n)) - y = JuMP.@variable(m, [1:k], Bin, base_name = "y_$counter") - _H = reflected_gray_codes(k) - d = length(_H) - H = Dict(i => _H[i] for i in 1:d) - H[0] = H[1] - H[d+1] = H[d] - for j in 1:k - JuMP.@constraints( - m, - begin - sum(λ[i] for i in 1:(n+1) if H[i-1][j] == H[i][j] == 1) ≤ y[j] - sum(λ[i] for i in 1:(n+1) if H[i-1][j] == H[i][j] == 0) ≤ - 1 - y[j] - end - ) - end - return -end - -function sos2_zigzag_formulation!(m::JuMP.Model, λ) - counter = m.ext[:PWL].counter - n = length(λ) - 1 - k = ceil(Int, log2(n)) - y = JuMP.@variable(m, [1:k], Bin, base_name = "y_$counter") - sos2_encoding_constraints!(m, λ, y, zigzag_codes(k), zigzag_hyperplanes(k)) - return -end - -function sos2_zigzag_general_integer_formulation!(m::JuMP.Model, λ) - counter = m.ext[:PWL].counter - n = length(λ) - 1 - k = ceil(Int, log2(n)) - # TODO: tighter upper_bounds - y = JuMP.@variable( - m, - [i = 1:k], - Int, - lower_bound = 0, - upper_bound = 2^(k - i), - base_name = "y_$counter" - ) - sos2_encoding_constraints!( - m, - λ, - y, - integer_zigzag_codes(k), - unit_vector_hyperplanes(k), - ) - return -end - -function sos2_generalized_celaya_formulation!(m::JuMP.Model, λ) - counter = m.ext[:PWL].counter - n = length(λ) - 1 - k = ceil(Int, log2(n)) - codes = generalized_celaya_codes(k) - lb = [minimum(t[i] for t in codes) for i in 1:k] - ub = [maximum(t[i] for t in codes) for i in 1:k] - y = JuMP.@variable( - m, - [i = 1:k], - Int, - lower_bound = lb[i], - upper_bound = ub[i], - base_name = "y_$counter" - ) - sos2_encoding_constraints!( - m, - λ, - y, - codes, - generalized_celaya_hyperplanes(k), - ) - return -end - -function sos2_symmetric_celaya_formulation!(m::JuMP.Model, λ) - counter = m.ext[:PWL].counter - n = length(λ) - 1 - k = ceil(Int, log2(n)) - codes = symmetric_celaya_codes(k) - lb = [minimum(t[i] for t in codes) for i in 1:k] - ub = [maximum(t[i] for t in codes) for i in 1:k] - y = JuMP.@variable( - m, - [i = 1:k], - Int, - lower_bound = lb[i], - upper_bound = ub[i], - base_name = "y_$counter" - ) - sos2_encoding_constraints!(m, λ, y, codes, symmetric_celaya_hyperplanes(k)) - return -end - -function sos2_encoding_constraints!(m, λ, y, h, B) - n = length(λ) - 1 - for b in B - JuMP.@constraints( - m, - begin - 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 -end - -function reflected_gray_codes(k::Int) - if k == 0 - return Vector{Int}[] - elseif k == 1 - return [[0], [1]] - else - codes′ = reflected_gray_codes(k - 1) - return vcat( - [vcat(code, 0) for code in codes′], - [vcat(code, 1) for code in reverse(codes′)], - ) - end -end - -function zigzag_codes(k::Int) - if k <= 0 - error("Invalid code length $k") - elseif k == 1 - return [[0], [1]] - else - codes′ = zigzag_codes(k - 1) - return vcat( - [vcat(code, 0) for code in codes′], - [vcat(code, 1) for code in codes′], - ) - end -end - -function integer_zigzag_codes(k::Int) - if k <= 0 - error("Invalid code length $k") - elseif k == 1 - return [[0], [1]] - else - codes′ = integer_zigzag_codes(k - 1) - offset = [2^(j - 2) for j in k:-1:2] - return vcat( - [vcat(code, 0) for code in codes′], - [vcat(code .+ offset, 1) for code in codes′], - ) - end -end - -function generalized_celaya_codes(k::Int) - if k <= 0 - error("Invalid code length $k") - elseif k == 1 - return [[0], [1]] - elseif k == 2 - return [[0, 0], [1, 1], [1, 0], [0, 1]] - else - codes′ = generalized_celaya_codes(k - 1) - n = length(codes′) - hp = Int(n / 2) - firstcodes = [codes′[i] for i in 1:hp] - secondcodes = [codes′[i] for i in (hp+1):n] - return vcat( - [vcat(codes′[i], 0) for i in 1:n], - [vcat(codes′[i], 1) for i in 1:hp], - [vcat(codes′[i], -1) for i in (hp+1):n], - ) - end -end - -function symmetric_celaya_codes(k::Int) - if k <= 0 - error("Invalid code length $k") - elseif k == 1 - return [[0], [1]] - else - codes′ = generalized_celaya_codes(k - 1) - n = length(codes′) - hp = Int(n / 2) - firstcodes = [codes′[i] for i in 1:hp] - secondcodes = [codes′[i] for i in (hp+1):n] - return vcat( - [vcat(codes′[i], 0) for i in 1:n], - [vcat(codes′[i], 1) for i in 1:hp], - [vcat(codes′[i], -1) for i in (hp+1):n], - ) - end -end - -function zigzag_hyperplanes(k::Int) - hps = Vector{Int}[] - for i in 1:k - hp = zeros(Int, k) - hp[i] = 1 - for j in (i+1):k - hp[j] = 2^(j - i - 1) - end - push!(hps, hp) - end - return hps -end - -function unit_vector_hyperplanes(k::Int) - hps = Vector{Int}[] - for i in 1:k - hp = zeros(Int, k) - hp[i] = 1 - push!(hps, hp) - end - return hps -end - -function generalized_celaya_hyperplanes(k::Int) - return compute_hyperplanes(generalized_celaya_codes(k)) -end - -function symmetric_celaya_hyperplanes(k::Int) - return compute_hyperplanes(symmetric_celaya_codes(k)) -end - -function compute_hyperplanes(C::Vector{Vector{T}}) where {T<:Number} - n = length(C) - k = length(C[1]) - d = zeros(n - 1, k) - for i in 1:n-1 - d[i, :] = C[i+1] - C[i] - end - if k <= 1 - error("Cannot process codes of length $k") - elseif k == 2 - @assert n == 4 - spanners = Vector{Float64}[] - for i in 1:n-1 - v = canonical!(d[i, :]) - v = [v[2], -v[1]] - push!(spanners, v) - end - indices = [1] - approx12 = isapprox(spanners[1], spanners[2]) - approx13 = isapprox(spanners[1], spanners[3]) - approx23 = isapprox(spanners[2], spanners[3]) - if !approx12 - push!(indices, 2) - end - if !approx13 && !(!approx12 && approx23) - push!(indices, 3) - end - return spanners[indices] - end - indices = [1, 2] - spanners = Vector{Float64}[] - while !isempty(indices) - if indices == [n - 1] - break - end - if LinearAlgebra.rank(d[indices, :]) == length(indices) && - length(indices) <= k - 1 - if length(indices) == k - 1 - nullsp = LinearAlgebra.nullspace(d[indices, :]) - @assert size(nullsp, 2) == 1 - v = vec(nullsp) - push!(spanners, canonical!(v)) - end - if indices[end] != n - 1 - push!(indices, indices[end] + 1) - else - pop!(indices) - indices[end] += 1 - end - else - if indices[end] != n - 1 - indices[end] += 1 - else - pop!(indices) - indices[end] += 1 - end - end - end - keepers = [1] - for i in 2:length(spanners) - alreadyin = false - for j in keepers - if isapprox(spanners[i], spanners[j]) - alreadyin = true - break - end - end - if !alreadyin - push!(keepers, i) - end - end - return spanners[keepers] -end - -function canonical!(v::Vector{Float64}) - LinearAlgebra.normalize!(v) - for j in 1:length(v) - if abs(v[j]) < 1e-8 - v[j] = 0 - end - end - sgn = sign(v[findfirst([x != 0 for x in v])]) - for j in 1:length(v) - if abs(v[j]) < 1e-8 - v[j] = 0 - else - v[j] *= sgn - end - end - return v -end - -function optimal_IB_scheme!(m::JuMP.Model, λ, pwl, subsolver) # IB: independent branching - m.ext[:OptimalIB] = Int[] - if !haskey(m.ext, :OptimalIBCache) - m.ext[:OptimalIBCache] = Dict() - end - T = pwl.T - n = maximum(maximum(t) for t in T) - J = 1:n - E = fill(false, n, n) - for t in T - for i in t, j in t - E[i, j] = true - end - end - if haskey(m.ext[:OptimalIBCache], pwl.T) - xx, yy = m.ext[:OptimalIBCache][pwl.T] - t = size(xx, 1) - else - t = ceil(Int, log2(length(T))) - xx = Array{Float64,2}(undef, 0, 0) - yy = Array{Float64,2}(undef, 0, 0) - while true - model = JuMP.Model(; solver = subsolver) - JuMP.@variable(model, x[1:t, 1:n], Bin) - JuMP.@variable(model, y[1:t, 1:n], Bin) - JuMP.@variable(model, z[1:t, 1:n, 1:n], Bin) - for j in 1:t - for r in J, s in J - if r >= s - continue - end - JuMP.@constraints( - model, - begin - z[j, r, s] <= x[j, r] + x[j, s] - z[j, r, s] <= x[j, r] + y[j, r] - z[j, r, s] <= x[j, s] + y[j, s] - z[j, r, s] <= y[j, r] + y[j, s] - z[j, r, s] >= x[j, r] + y[j, s] - 1 - z[j, r, s] >= x[j, s] + y[j, r] - 1 - end - ) - end - for r in J - JuMP.@constraint(model, x[j, r] + y[j, r] <= 1) - end - end - for r in J, s in J - if r >= s - continue - end - if E[r, s] - JuMP.@constraint(model, sum(z[j, r, s] for j in 1:t) == 0) - else - JuMP.@constraint(model, sum(z[j, r, s] for j in 1:t) >= 1) - end - end - JuMP.@objective(model, Min, sum(x) + sum(y)) - stat = JuMP.solve(model) - xx = JuMP.getvalue(x) - yy = JuMP.getvalue(y) - if any(isnan, xx) || any(isnan, yy) - t += 1 - else - break - end - end - m.ext[:OptimalIBCache][pwl.T] = (xx, yy) - end - y = JuMP.@variable(m, [1:t], Bin) - dˣ = [_x[1] for _x in pwl.x] - dʸ = [_x[2] for _x in pwl.x] - uˣ, uʸ = unique(dˣ), unique(dʸ) - @assert issorted(uˣ) - @assert issorted(uʸ) - nˣ, nʸ = length(uˣ), length(uʸ) - ˣtoⁱ = Dict(uˣ[i] => i for i in 1:nˣ) - ʸtoʲ = Dict(uʸ[i] => i for i in 1:nʸ) - for i in 1:t - JuMP.@constraints( - m, - begin - sum( - λ[ˣtoⁱ[pwl.x[j][1]], ʸtoʲ[pwl.x[j][2]]] for - j in J if xx[i, j] == 1 - ) ≤ y[i] - sum( - λ[ˣtoⁱ[pwl.x[j][1]], ʸtoʲ[pwl.x[j][2]]] for - j in J if yy[i, j] == 1 - ) ≤ 1 - y[i] - end - ) - end - push!(m.ext[:OptimalIB], t) - return -end - -function piecewiselinear( - m::JuMP.Model, - x::VarOrAff, - y::VarOrAff, - dˣ, - dʸ, - f::Function; - method = defaultmethod(), -) - return piecewiselinear( - m, - x, - y, - BivariatePWLFunction(dˣ, dʸ, f); - method = method, - ) -end - -function piecewiselinear( - m::JuMP.Model, - x₁::VarOrAff, - x₂::VarOrAff, - pwl::BivariatePWLFunction; - method = defaultmethod(), - subsolver = nothing, -) - if (method == :OptimalIB) && (subsolver === nothing) - error( - "No MIP solver provided to construct optimal IB scheme. Pass a solver object to the piecewiselinear function, e.g. piecewiselinear(m, x₁, x₂, bivariatefunc, method=:OptimalIB, subsolver=GurobiSolver())", - ) - end - initPWL!(m) - counter = m.ext[:PWL].counter - counter += 1 - m.ext[:PWL].counter = counter - dˣ = [_x[1] for _x in pwl.x] - dʸ = [_x[2] for _x in pwl.x] - uˣ, uʸ = unique(dˣ), unique(dʸ) - @assert issorted(uˣ) - @assert issorted(uʸ) - T = pwl.T - nˣ, nʸ = length(uˣ), length(uʸ) - if nˣ == 0 || nʸ == 0 - error( - "I don't know how to handle a piecewise linear function with zero breakpoints", - ) - elseif nˣ == 1 - @assert length(dʸ) == nʸ == length(pwl.z) - return piecewiselinear( - m, - x₂, - UnivariatePWLFunction(dʸ, [pwl.z[i] for i in 1:nʸ]), - ) - elseif nʸ == 1 - @assert length(dˣ) == nˣ == length(pwl.z) - return piecewiselinear( - m, - x₁, - UnivariatePWLFunction(dˣ, [pwl.z[i] for i in 1:nˣ]), - ) - end - ˣtoⁱ = Dict(uˣ[i] => i for i in 1:nˣ) - ʸtoʲ = Dict(uʸ[i] => i for i in 1:nʸ) - fd = Array{Float64}(undef, nˣ, nʸ) - for (v, fv) in zip(pwl.x, pwl.z) - # i is the linear index into pwl.x...really want (i,j) pair - fd[ˣtoⁱ[v[1]], ʸtoʲ[v[2]]] = fv - end - z = JuMP.@variable( - m, - lower_bound = minimum(fd), - upper_bound = maximum(fd), - base_name = "z_$counter" - ) - if method == :MC - x̂₁ = JuMP.@variable(m, [T], base_name = "x̂₁_$counter") - x̂₂ = JuMP.@variable(m, [T], base_name = "x̂₂_$counter") - ẑ = JuMP.@variable(m, [T], base_name = "ẑ_$counter") - y = JuMP.@variable(m, [T], Bin, base_name = "y_$counter") - JuMP.@constraint(m, sum(y) == 1) - JuMP.@constraint(m, sum(x̂₁) == x₁) - JuMP.@constraint(m, sum(x̂₂) == x₂) - JuMP.@constraint(m, sum(ẑ) == z) - for t in T - @assert length(t) == 3 - r¹, r², r³ = pwl.x[t[1]], pwl.x[t[2]], pwl.x[t[3]] - fz¹, fz², fz³ = pwl.z[t[1]], pwl.z[t[2]], pwl.z[t[3]] - for P in ([1, 2, 3], [2, 3, 1], [3, 1, 2]) - p¹, p², p³ = [r¹, r², r³][P] - # p¹, p², p³ = pwl.x[t[P[1]]], pwl.x[t[P[2]]], pwl.x[t[P[3]]] - - A = [ - p¹[1] p¹[2] 1 - p²[1] p²[2] 1 - p³[1] p³[2] 1 - ] - @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, - ) - @assert isapprox( - q[1] * p²[1] + q[2] * p²[2] + q[3], - 0, - atol = 1e-4, - ) - @assert isapprox( - q[1] * p³[1] + q[2] * p³[2] + q[3], - 1, - atol = 1e-4, - ) - JuMP.@constraint( - m, - q[1] * x̂₁[t] + q[2] * x̂₂[t] + q[3] * y[t] ≥ 0 - ) - end - A = [ - r¹[1] r¹[2] 1 - r²[1] r²[2] 1 - r³[1] r³[2] 1 - ] - b = [fz¹, fz², fz³] - q = A \ b - @assert isapprox( - q[1] * r¹[1] + q[2] * r¹[2] + q[3], - fz¹, - atol = 1e-4, - ) - @assert isapprox( - q[1] * r²[1] + q[2] * r²[2] + q[3], - fz², - atol = 1e-4, - ) - @assert isapprox( - q[1] * r³[1] + q[2] * r³[2] + q[3], - fz³, - atol = 1e-4, - ) - JuMP.@constraint( - m, - ẑ[t] == q[1] * x̂₁[t] + q[2] * x̂₂[t] + q[3] * y[t] - ) - end - elseif method in (:DisaggLogarithmic, :DLog) - T = pwl.T - X = pwl.x - Z = pwl.z - n = length(X) - γ = JuMP.@variable( - m, - [t = T, v = t], - lower_bound = 0, - upper_bound = 1, - base_name = "γ_$counter" - ) - Tv = Dict(v => Any[] for v in 1:n) - for t in T, v in t - push!(Tv[v], t) - end - JuMP.@constraint(m, sum(γ) == 1) - JuMP.@constraint( - m, - sum(sum(γ[t, i] for t in Tv[i]) * X[i][1] for i in 1:n) == x₁ - ) - JuMP.@constraint( - m, - sum(sum(γ[t, i] for t in Tv[i]) * X[i][2] for i in 1:n) == x₂ - ) - JuMP.@constraint( - m, - sum(sum(γ[t, i] for t in Tv[i]) * Z[i] for i in 1:n) == z - ) - r = ceil(Int, log2(length(T))) - H = reflected_gray_codes(r) - y = JuMP.@variable(m, [1:r], Bin, base_name = "y_$counter") - for j in 1:r - JuMP.@constraint( - m, - sum( - sum(γ[T[i], v] for v in T[i]) * H[i][j] for i in 1:length(T) - ) == y[j] - ) - end - else - # V-formulation methods - λ = JuMP.@variable( - m, - [1:nˣ, 1:nʸ], - lower_bound = 0, - upper_bound = 1, - base_name = "λ_$counter" - ) - JuMP.@constraint(m, sum(λ) == 1) - JuMP.@constraint(m, sum(λ[i, j] * uˣ[i] for i in 1:nˣ, j in 1:nʸ) == x₁) - JuMP.@constraint(m, sum(λ[i, j] * uʸ[j] for i in 1:nˣ, j in 1:nʸ) == x₂) - JuMP.@constraint( - m, - sum(λ[i, j] * fd[i, j] for i in 1:nˣ, j in 1:nʸ) == z - ) - if method == :CC - T = pwl.T - y = JuMP.@variable(m, [T], Bin, base_name = "y_$counter") - JuMP.@constraint(m, sum(y) == 1) - Ts = Dict((i, j) => Vector{Int}[] for i in 1:nˣ, j in 1:nʸ) - for t in T, ind in t - rˣ, rʸ = pwl.x[ind] - i, j = ˣtoⁱ[rˣ], ʸtoʲ[rʸ] - push!(Ts[(i, j)], t) - end - for i in 1:nˣ, j in 1:nʸ - JuMP.@constraint(m, λ[i, j] ≤ sum(y[t] for t in Ts[(i, j)])) - end - elseif method == :Incremental - error( - "Incremental formulation for bivariate functions is not currently implemented.", - ) - # TODO: implement algorithm of Geissler et al. (2012) - elseif method == :OptimalIB - error( - "Optimal IB formulation has not yet been updated for JuMP v0.19.", - ) - optimal_IB_scheme!(m, λ, pwl, subsolver) - else - # formulations with SOS2 along each dimension - Tx = [sum(λ[tˣ, tʸ] for tˣ in 1:nˣ) for tʸ in 1:nʸ] - Ty = [sum(λ[tˣ, tʸ] for tʸ in 1:nʸ) for tˣ in 1:nˣ] - - if method in (:Logarithmic, :Log) - sos2_logarithmic_formulation!(m, Tx) - sos2_logarithmic_formulation!(m, Ty) - elseif method in (:LogarithmicIB, :LogIB) - sos2_logarithmic_IB_formulation!(m, Tx) - sos2_logarithmic_IB_formulation!(m, Ty) - elseif method in (:ZigZag, :ZZB) - sos2_zigzag_formulation!(m, Tx) - sos2_zigzag_formulation!(m, Ty) - elseif method in (:ZigZagInteger, :ZZI) - sos2_zigzag_general_integer_formulation!(m, Tx) - sos2_zigzag_general_integer_formulation!(m, Ty) - elseif method == :GeneralizedCelaya - sos2_generalized_celaya_formulation!(m, Tx) - sos2_generalized_celaya_formulation!(m, Ty) - elseif method == :SymmetricCelaya - sos2_symmetric_celaya_formulation!(m, Tx) - sos2_symmetric_celaya_formulation!(m, Ty) - elseif method == :SOS2 - γˣ = JuMP.@variable( - m, - [1:nˣ], - lower_bound = 0, - upper_bound = 1, - base_name = "γˣ_$counter" - ) - γʸ = JuMP.@variable( - m, - [1:nʸ], - lower_bound = 0, - upper_bound = 1, - base_name = "γʸ_$counter" - ) - JuMP.@constraint( - m, - [tˣ in 1:nˣ], - γˣ[tˣ] == sum(λ[tˣ, tʸ] for tʸ in 1:nʸ) - ) - JuMP.@constraint( - m, - [tʸ in 1:nʸ], - γʸ[tʸ] == sum(λ[tˣ, tʸ] for tˣ in 1:nˣ) - ) - JuMP.@constraint(m, γˣ in MOI.SOS2([k for k in 1:nˣ])) - JuMP.@constraint(m, γʸ in MOI.SOS2([k for k in 1:nʸ])) - else - error("Unrecognized method $method") - end - - pattern = pwl.meta[:structure] - if pattern == :UnionJack - numT = 0 - minˣ, minʸ = uˣ[1], uʸ[1] - # find the index of the bottom-left point - idx = findfirst(pwl.x) do w - return w == (minˣ, minʸ) - end - for t in pwl.T - if idx in t - numT += 1 - end - end - @assert 1 <= numT <= 2 - - w = JuMP.@variable(m, binary = true, base_name = "w_$counter") - # if numT==1, then bottom-left is contained in only one point, and so needs separating; otherwise numT==2, and need to offset by one - JuMP.@constraints( - m, - begin - sum(λ[tx, ty] for tx in 1:2:nˣ, ty in numT:2:nʸ) ≤ w - sum(λ[tx, ty] for tx in 2:2:nˣ, ty in (3-numT):2:nʸ) ≤ - 1 - w - end - ) - elseif pattern == :K1 - # assumption is that the triangulation is uniform, as defined in types.jl - w = JuMP.@variable(m, [1:2], Bin, base_name = "w_$counter") - JuMP.@constraints( - m, - begin - sum( - λ[tx, ty] for tx in 1:nˣ, ty in 1:nʸ if - mod(tx, 2) == mod(ty, 2) && (tx + ty) in 2:4:(nˣ+nʸ) - ) ≤ w[1] - sum( - λ[tx, ty] for tx in 1:nˣ, ty in 1:nʸ if - mod(tx, 2) == mod(ty, 2) && (tx + ty) in 4:4:(nˣ+nʸ) - ) ≤ 1 - w[1] - sum( - λ[tx, ty] for tx in 1:nˣ, ty in 1:nʸ if - mod(tx, 2) != mod(ty, 2) && (tx + ty) in 3:4:(nˣ+nʸ) - ) ≤ w[2] - sum( - λ[tx, ty] for tx in 1:nˣ, ty in 1:nʸ if - mod(tx, 2) != mod(ty, 2) && (tx + ty) in 5:4:(nˣ+nʸ) - ) ≤ 1 - w[2] - end - ) - elseif pattern == :OptimalTriangleSelection - error( - "Optimal triangle selection formulation has not yet been updated for JuMP v0.19.", - ) - m.ext[:OptimalTriSelect] = Int[] - - if !haskey(m.ext, :OptimalTriSelectCache) - m.ext[:OptimalTriSelectCache] = Dict() - end - - J = [(i, j) for i in 1:nˣ, j in 1:nʸ] - E = Set{Tuple{Tuple{Int,Int},Tuple{Int,Int}}}() - for t in T - @assert length(t) == 3 - IJ = [(ˣtoⁱ[pwl.x[i][1]], ʸtoʲ[pwl.x[i][2]]) for i in t] - im = minimum(ij[1] for ij in IJ) - iM = maximum(ij[1] for ij in IJ) - jm = minimum(ij[2] for ij in IJ) - jM = maximum(ij[2] for ij in IJ) - @assert im < iM - @assert im < iM - if ((im, jm) in IJ) && ((iM, jM) in IJ) - push!(E, ((im, jM), (iM, jm))) - elseif ((im, jM) in IJ) && ((iM, jm) in IJ) - push!(E, ((im, jm), (iM, jM))) - else - error() - end - end - - if haskey(m.ext[:OptimalTriSelectCache], E) - xx, yy = m.ext[:OptimalTriSelectCache][E] - t = JuMP.size(xx, 1) - else - if subsolver === nothing - error( - "No MIP solver provided to construct optimal triangle selection. Pass a solver object to the piecewiselinear function, e.g. piecewiselinear(m, x₁, x₂, bivariatefunc, method=:Logarithmic, subsolver=GurobiSolver())", - ) - end - t = 1 - xx, yy = Array(Float64, 0, 0), Array(Float64, 0, 0) - while true - subm = JuMP.Model(; solver = subsolver) - JuMP.@variable(subm, xˢᵘᵇ[1:t, J], Bin) - JuMP.@variable(subm, yˢᵘᵇ[1:t, J], Bin) - JuMP.@variable(subm, zˢᵘᵇ[1:t, J, J], Bin) - - for j in 1:t - for r in J, s in J - # lexicographic ordering on points on grid - if r[1] > s[1] || (r[1] == s[1] && r[2] >= s[2]) - continue - end - JuMP.@constraints( - subm, - begin - zˢᵘᵇ[j, r, s] <= - xˢᵘᵇ[j, r] + xˢᵘᵇ[j, s] - zˢᵘᵇ[j, r, s] <= - xˢᵘᵇ[j, r] + yˢᵘᵇ[j, r] - zˢᵘᵇ[j, r, s] <= - xˢᵘᵇ[j, s] + yˢᵘᵇ[j, s] - zˢᵘᵇ[j, r, s] <= - yˢᵘᵇ[j, r] + yˢᵘᵇ[j, s] - zˢᵘᵇ[j, r, s] >= - xˢᵘᵇ[j, r] + yˢᵘᵇ[j, s] - 1 - zˢᵘᵇ[j, r, s] >= - xˢᵘᵇ[j, s] + yˢᵘᵇ[j, r] - 1 - end - ) - end - for r in J - JuMP.@constraint( - subm, - xˢᵘᵇ[j, r] + yˢᵘᵇ[j, r] <= 1 - ) - end - end - - for r in J, s in J - # lexicographic ordering on points on grid - (r[1] > s[1] || (r[1] == s[1] && r[2] >= s[2])) && - continue - if (r, s) in E - JuMP.@constraint( - subm, - sum(zˢᵘᵇ[j, r, s] for j in 1:t) >= 1 - ) - elseif max(abs(r[1] - s[1]), abs(r[2] - s[2])) == 1 - JuMP.@constraint( - subm, - sum(zˢᵘᵇ[j, r, s] for j in 1:t) == 0 - ) - end - end - - JuMP.@objective(subm, Min, sum(xˢᵘᵇ) + sum(yˢᵘᵇ)) - stat = JuMP.solve(subm) - if any(isnan, subm.colVal) - t += 1 - else - xx = JuMP.getvalue(xˢᵘᵇ) - yy = JuMP.getvalue(yˢᵘᵇ) - m.ext[:OptimalTriSelectCache][E] = (xx, yy) - break - end - end - end - y = JuMP.@variable( - m, - [1:t], - Bin, - base_name = "Δselect_$counter" - ) - - for i in 1:t - JuMP.@constraints( - m, - begin - sum(λ[v[1], v[2]] for v in J if xx[i, v] ≈ 1) ≤ - y[i] - sum(λ[v[1], v[2]] for v in J if yy[i, v] ≈ 1) ≤ - 1 - y[i] - end - ) - end - push!(m.ext[:OptimalTriSelect], t) - elseif pattern in (:Stencil, :Stencil9) - w = JuMP.@variable(m, [1:3, 1:3], Bin, base_name = "w_$counter") - for oˣ in 1:3, oʸ in 1:3 - innoT = fill(true, nˣ, nʸ) - for (i, j, k) in pwl.T - xⁱ, xʲ, xᵏ = pwl.x[i], pwl.x[j], pwl.x[k] - iiˣ, iiʸ = ˣtoⁱ[xⁱ[1]], ʸtoʲ[xⁱ[2]] - jjˣ, jjʸ = ˣtoⁱ[xʲ[1]], ʸtoʲ[xʲ[2]] - kkˣ, kkʸ = ˣtoⁱ[xᵏ[1]], ʸtoʲ[xᵏ[2]] - # check to see if one of the points in the triangle falls on the grid - if (mod1(iiˣ, 3) == oˣ && mod1(iiʸ, 3) == oʸ) || - (mod1(jjˣ, 3) == oˣ && mod1(jjʸ, 3) == oʸ) || - (mod1(kkˣ, 3) == oˣ && mod1(kkʸ, 3) == oʸ) - innoT[iiˣ, iiʸ] = false - innoT[jjˣ, jjʸ] = false - innoT[kkˣ, kkʸ] = false - end - end - JuMP.@constraints( - m, - begin - sum(λ[i, j] for i in oˣ:3:nˣ, j in oʸ:3:nʸ) ≤ - 1 - w[oˣ, oʸ] - sum( - λ[i, j] for i in 1:nˣ, j in 1:nʸ if innoT[i, j] - ) ≤ w[oˣ, oʸ] - end - ) - end - else - @assert pattern in (:Upper, :Lower, :BestFit, :Random) - # Eⁿᵉ[i,j] = true means that we must cover the edge {(i,j),(i+1,j+1)} - Eⁿᵉ = fill(false, nˣ - 1, nʸ - 1) - for (i, j, k) in pwl.T - xⁱ, xʲ, xᵏ = pwl.x[i], pwl.x[j], pwl.x[k] - iiˣ, iiʸ = ˣtoⁱ[xⁱ[1]], ʸtoʲ[xⁱ[2]] - jjˣ, jjʸ = ˣtoⁱ[xʲ[1]], ʸtoʲ[xʲ[2]] - kkˣ, kkʸ = ˣtoⁱ[xᵏ[1]], ʸtoʲ[xᵏ[2]] - IJ = [(iiˣ, iiʸ), (jjˣ, jjʸ), (kkˣ, kkʸ)] - im = min(iiˣ, jjˣ, kkˣ) - iM = max(iiˣ, jjˣ, kkˣ) - jm = min(iiʸ, jjʸ, kkʸ) - jM = max(iiʸ, jjʸ, kkʸ) - if ((im, jM) in IJ) && ((iM, jm) in IJ) - Eⁿᵉ[im, jm] = true - else - @assert (im, jm) in IJ && (iM, jM) in IJ - end - end - - # diagonal lines running from SW to NE. Grouped with an offset of 3. - wⁿᵉ = JuMP.@variable(m, [0:2], Bin, base_name = "wⁿᵉ_$counter") - for o in 0:2 - Aᵒ = Set{Tuple{Int,Int}}() - Bᵒ = Set{Tuple{Int,Int}}() - for offˣ in o:3:(nˣ-2) - SWinA = true # whether we put the SW corner of the next triangle to cover in set A - for i in (1+offˣ):(nˣ-1) - j = i - offˣ - if !(1 ≤ i ≤ nˣ - 1) - continue - end - if !(1 ≤ j ≤ nʸ - 1) - continue # should never happen - end - if Eⁿᵉ[i, j] # if we need to cover the edge... - if SWinA # figure out which set we need to put it in; this depends on previous triangle in our current line - push!(Aᵒ, (i, j)) - push!(Bᵒ, (i + 1, j + 1)) - else - push!(Aᵒ, (i + 1, j + 1)) - push!(Bᵒ, (i, j)) - end - SWinA = !SWinA - end - end - end - for offʸ in (3-o):3:(nʸ-1) - SWinA = true - for j in (offʸ+1):(nʸ-1) - i = j - offʸ - if !(1 ≤ i ≤ nˣ - 1) - continue - end - if Eⁿᵉ[i, j] - if SWinA - push!(Aᵒ, (i, j)) - push!(Bᵒ, (i + 1, j + 1)) - else - push!(Aᵒ, (i + 1, j + 1)) - push!(Bᵒ, (i, j)) - end - SWinA = !SWinA - end - end - end - JuMP.@constraints( - m, - begin - sum(λ[i, j] for (i, j) in Aᵒ) ≤ wⁿᵉ[o] - sum(λ[i, j] for (i, j) in Bᵒ) ≤ 1 - wⁿᵉ[o] - end - ) - end - - wˢᵉ = JuMP.@variable(m, [0:2], Bin, base_name = "wˢᵉ_$counter") - for o in 0:2 - Aᵒ = Set{Tuple{Int,Int}}() - Bᵒ = Set{Tuple{Int,Int}}() - for offˣ in o:3:(nˣ-2) - SEinA = true - # for i in (1+offˣ):-1:1 - # j = offˣ - i + 2 - for j in 1:(nʸ-1) - i = nˣ - j - offˣ - if !(1 ≤ i ≤ nˣ - 1) - continue - end - if !Eⁿᵉ[i, j] - if SEinA - push!(Aᵒ, (i + 1, j)) - push!(Bᵒ, (i, j + 1)) - else - push!(Aᵒ, (i, j + 1)) - push!(Bᵒ, (i + 1, j)) - end - SEinA = !SEinA - end - end - end - for offʸ in (3-o):3:(nʸ-1) - SEinA = true - for j in (offʸ+1):(nʸ-1) - i = nˣ - j + offʸ - if !(1 ≤ i ≤ nˣ - 1) - continue - end - if !Eⁿᵉ[i, j] - if SEinA - push!(Aᵒ, (i + 1, j)) - push!(Bᵒ, (i, j + 1)) - else - push!(Aᵒ, (i, j + 1)) - push!(Bᵒ, (i + 1, j)) - end - SEinA = !SEinA - end - end - end - JuMP.@constraints( - m, - begin - sum(λ[i, j] for (i, j) in Aᵒ) ≤ wˢᵉ[o] - sum(λ[i, j] for (i, j) in Bᵒ) ≤ 1 - wˢᵉ[o] - end - ) - end - end - end - end - return z -end diff --git a/src/methods/bivariate/common.jl b/src/methods/bivariate/common.jl new file mode 100644 index 0000000..fd01de7 --- /dev/null +++ b/src/methods/bivariate/common.jl @@ -0,0 +1,107 @@ +# Copyright (c) 2016: Joey Huchette and contributors +# +# 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. + +const BivariateSOS2Method = + Union{K1,OptimalTriangleSelection,NineStencil,SixStencil,UnionJack} + +function formulate_pwl!( + model::JuMP.Model, + input_vars::NTuple{2,VarOrAff}, + output_vars::NTuple{F,VarOrAff}, + pwl::PWLFunctionPointRep{2,F}, + method::BivariateSOS2Method, + direction::DIRECTION, +) where {F} + initPWL!(model) + counter = model.ext[:PWL].counter + counter += 1 + model.ext[:PWL].counter = counter + + grid = _continuous_gridpoints_or_die(pwl) + + xs, ys = grid.input_vals, grid.output_vals + u_1 = unique(x[1] for x in xs) + u_2 = unique(x[2] for x in xs) + @assert issorted(u_1) + @assert issorted(u_2) + + n_1, n_2 = size(xs) + + x_1_to_i = Dict(u_1[i] => i for i in 1:n_1) + x_2_to_j = Dict(u_2[j] => j for j in 1:n_2) + + λ = JuMP.@variable( + model, + [1:n_1, 1:n_2], + lower_bound = 0, + upper_bound = 1, + base_name = "λ_$counter" + ) + JuMP.@constraint(model, sum(λ) == 1) + JuMP.@constraint( + model, + sum(λ[i, j] * u_1[i] for i in 1:n_1, j in 1:n_2) == input_vars[1] + ) + JuMP.@constraint( + model, + sum(λ[i, j] * u_2[j] for i in 1:n_1, j in 1:n_2) == input_vars[2] + ) + for k in 1:F + rhs = sum(λ[i, j] * ys[i, j][k] for i in 1:n_1, j in 1:n_2) + _constrain_output_var(model, output_vars[k], rhs, direction) + end + + # formulations with SOS2 along each dimension + T_x = [sum(λ[i, j] for i in 1:n_1) for j in 1:n_2] + T_y = [sum(λ[i, j] for j in 1:n_2) for i in 1:n_1] + + formulate_sos2!(model, T_x, axis_method(method)) + formulate_sos2!(model, T_y, axis_method(method)) + + canonical_input_segments = _canonicalize_triangulation(pwl, grid) + # triangle_direction[i,j] = true means that triangle is of the form + # ----- + # | \ | + # ----- + # If false, it is of the form + # ----- + # | / | + # ----- + triangle_direction = fill(false, (n_1 - 1, n_2 - 1)) + for input_seg in canonical_input_segments + # TODO: Remove; assertion redundant with _check_triangulation above. + @assert length(input_seg) == 3 + i_min, i_max = + extrema([input_seg[1][1], input_seg[2][1], input_seg[3][1]]) + j_min, j_max = + extrema([input_seg[1][2], input_seg[2][2], input_seg[3][2]]) + if ((i_min, j_max) in input_seg) && ((i_max, j_min) in input_seg) + triangle_direction[i_min, j_min] = true + else + @assert (i_min, j_min) in input_seg + @assert (i_max, j_max) in input_seg + end + end + + return formulate_triangle_selection!( + model, + λ, + triangle_direction, + method, + pwl.structure, + ) +end + +function formulate_triangle_selection!( + model::JuMP.Model, + λ::Matrix{JuMP.VariableRef}, + triangle_direction::Matrix{Bool}, + method::BivariateSOS2Method, + structure::GridTriangulation, +) + return error( + "The triangulation structure $structure is not suppported for method $method", + ) +end diff --git a/src/methods/bivariate/k1.jl b/src/methods/bivariate/k1.jl new file mode 100644 index 0000000..14a62a7 --- /dev/null +++ b/src/methods/bivariate/k1.jl @@ -0,0 +1,50 @@ +# Copyright (c) 2016: Joey Huchette and contributors +# +# 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. + +struct K1 <: Method + axis_method::Method +end +K1() = K1(Logarithmic()) + +axis_method(method::K1) = method.axis_method + +function formulate_triangle_selection!( + model::JuMP.Model, + λ::Matrix{JuMP.VariableRef}, + triangle_direction::Matrix{Bool}, + method::K1, + structure::K1Triangulation, +) + n_1, n_2 = size(λ) + @assert size(triangle_direction) == (n_1 - 1, n_2 - 1) + counter = model.ext[:PWL].counter + + # TODO: Certify triangulation is "uniform" + + w = JuMP.@variable(model, [1:2], Bin, base_name = "w_$counter") + JuMP.@constraints( + model, + begin + sum( + λ[i, j] for i in 1:n_1, j in 1:n_2 if + mod(i, 2) == mod(j, 2) && (i + j) in 2:4:(n_1+n_2) + ) ≤ w[1] + sum( + λ[i, j] for i in 1:n_1, j in 1:n_2 if + mod(i, 2) == mod(j, 2) && (i + j) in 4:4:(n_1+n_2) + ) ≤ 1 - w[1] + sum( + λ[i, j] for i in 1:n_1, j in 1:n_2 if + mod(i, 2) != mod(j, 2) && (i + j) in 3:4:(n_1+n_2) + ) ≤ w[2] + sum( + λ[i, j] for i in 1:n_1, j in 1:n_2 if + mod(i, 2) != mod(j, 2) && (i + j) in 5:4:(n_1+n_2) + ) ≤ 1 - w[2] + end + ) + + return nothing +end diff --git a/src/methods/bivariate/nine_stencil.jl b/src/methods/bivariate/nine_stencil.jl new file mode 100644 index 0000000..318dc40 --- /dev/null +++ b/src/methods/bivariate/nine_stencil.jl @@ -0,0 +1,69 @@ +# Copyright (c) 2016: Joey Huchette and contributors +# +# 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. + +struct NineStencil <: Method + axis_method::Method +end +NineStencil() = NineStencil(Logarithmic()) + +axis_method(method::NineStencil) = method.axis_method + +# TODO: Unit tests for biclique cover +function formulate_triangle_selection!( + model::JuMP.Model, + λ::Matrix{JuMP.VariableRef}, + triangle_direction::Matrix{Bool}, + method::NineStencil, + _::GridTriangulation, +) + n_1, n_2 = size(λ) + @assert size(triangle_direction) == (n_1 - 1, n_2 - 1) + counter = model.ext[:PWL].counter + + w = JuMP.@variable(model, [1:3, 1:3], Bin, base_name = "w_$counter") + for o_1 in 1:3, o_2 in 1:3 + has_edge_across_stencil = Set{Tuple{Int,Int}}() + for i in o_1:3:n_1, j in o_2:3:n_2 + let k = i + 1, l = j + 1 + if (1 ≤ k ≤ n_1) && (1 ≤ l ≤ n_2) + if triangle_direction[i, j] + push!(has_edge_across_stencil, (k, l)) + end + end + end + let k = i + 1, l = j - 1 + if (1 ≤ k ≤ n_1) && (1 ≤ l ≤ n_2) + if !triangle_direction[i, j-1] + push!(has_edge_across_stencil, (k, l)) + end + end + end + let k = i - 1, l = j - 1 + if (1 ≤ k ≤ n_1) && (1 ≤ l ≤ n_2) + if triangle_direction[i-1, j-1] + push!(has_edge_across_stencil, (k, l)) + end + end + end + let k = i - 1, l = j + 1 + if (1 ≤ k ≤ n_1) && (1 ≤ l ≤ n_2) + if !triangle_direction[i-1, j] + push!(has_edge_across_stencil, (k, l)) + end + end + end + end + JuMP.@constraints( + model, + begin + sum(λ[i, j] for i in o_1:3:n_1, j in o_2:3:n_2) ≤ + 1 - w[o_1, o_2] + sum(λ[i, j] for (i, j) in has_edge_across_stencil) ≤ + w[o_1, o_2] + end + ) + end + return nothing +end diff --git a/src/methods/bivariate/optimal_independent_branching.jl b/src/methods/bivariate/optimal_independent_branching.jl new file mode 100644 index 0000000..4468915 --- /dev/null +++ b/src/methods/bivariate/optimal_independent_branching.jl @@ -0,0 +1,180 @@ +# Copyright (c) 2016: Joey Huchette and contributors +# +# 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. + +# TODO: Generalize to multivariate case. +struct OptimalIndependentBranching <: Method + sub_solver::Any +end + +function formulate_pwl!( + model::JuMP.Model, + input_vars::NTuple{2,VarOrAff}, + output_vars::NTuple{F,VarOrAff}, + pwl::PWLFunctionPointRep{2,F}, + method::OptimalIndependentBranching, + direction::DIRECTION, +) where {F} + initPWL!(model) + counter = model.ext[:PWL].counter + counter += 1 + model.ext[:PWL].counter = counter + + grid = _continuous_gridpoints_or_die(pwl) + + xs, ys = grid.input_vals, grid.output_vals + u_1 = unique(x[1] for x in xs) + u_2 = unique(x[2] for x in xs) + @assert issorted(u_1) + @assert issorted(u_2) + + n_1, n_2 = size(xs) + + x_1_to_i = Dict(u_1[i] => i for i in 1:n_1) + x_2_to_j = Dict(u_2[j] => j for j in 1:n_2) + + λ = JuMP.@variable( + model, + [1:n_1, 1:n_2], + lower_bound = 0, + upper_bound = 1, + base_name = "λ_$counter" + ) + JuMP.@constraint(model, sum(λ) == 1) + JuMP.@constraint( + model, + sum(λ[i, j] * u_1[i] for i in 1:n_1, j in 1:n_2) == input_vars[1] + ) + JuMP.@constraint( + model, + sum(λ[i, j] * u_2[j] for i in 1:n_1, j in 1:n_2) == input_vars[2] + ) + for k in 1:F + rhs = sum(λ[i, j] * ys[i, j][k] for i in 1:n_1, j in 1:n_2) + _constrain_output_var(model, output_vars[k], rhs, direction) + end + + canonical_input_segments = _canonicalize_triangulation(pwl, grid) + # triangle_direction[i,j] = true means that triangle is of the form + # ----- + # | \ | + # ----- + # If false, it is of the form + # ----- + # | / | + # ----- + triangle_direction = fill(false, (n_1 - 1, n_2 - 1)) + for input_seg in canonical_input_segments + # TODO: Remove; assertion redundant with _check_triangulation above. + @assert length(input_seg) == 3 + i_min, i_max = + extrema([input_seg[1][1], input_seg[2][1], input_seg[3][1]]) + j_min, j_max = + extrema([input_seg[1][2], input_seg[2][2], input_seg[3][2]]) + if ((i_min, j_max) in input_seg) && ((i_max, j_min) in input_seg) + triangle_direction[i_min, j_min] = true + else + @assert (i_min, j_min) in input_seg + @assert (i_max, j_max) in input_seg + end + end + + n_1, n_2 = size(λ) + J = Set((i, j) for i in 1:n_1, j in 1:n_2) + x_val = Matrix{Float64}(undef, n_1, n_2) + y_val = Matrix{Float64}(undef, n_1, n_2) + + t = ceil(Int, log2(2 * (n_1 - 1) * (n_2 - 1))) + while true + #@show method.sub_solver + sub_model = JuMP.Model(method.sub_solver) + JuMP.@variable(sub_model, x[1:t, J], Bin) + JuMP.@variable(sub_model, y[1:t, J], Bin) + JuMP.@variable(sub_model, z[1:t, J, J], Bin) + for j in 1:t + for r in J, s in J + # lexicographic ordering on points on grid + if r[1] > s[1] || (r[1] == s[1] && r[2] ≥ s[2]) + continue + end + JuMP.@constraints( + sub_model, + begin + z[j, r, s] ≤ x[j, r] + x[j, s] + z[j, r, s] ≤ x[j, r] + y[j, r] + z[j, r, s] ≤ x[j, s] + y[j, s] + z[j, r, s] ≤ y[j, r] + y[j, s] + z[j, r, s] ≥ x[j, r] + y[j, s] - 1 + z[j, r, s] ≥ x[j, s] + y[j, r] - 1 + end + ) + end + for r in J + JuMP.@constraint(sub_model, x[j, r] + y[j, r] ≤ 1) + end + end + + #@show J + for r in J, s in J + # lexicographic ordering on points on grid + if r[1] > s[1] || (r[1] == s[1] && r[2] ≥ s[2]) + continue + end + edge_is_far_away = LinearAlgebra.norm(r .- s, Inf) > 1 + edge_is_diagonal = + !edge_is_far_away && (abs(r[1] - s[1]) == abs(r[2] - s[2]) == 1) + if edge_is_far_away + JuMP.@constraint(sub_model, sum(z[j, r, s] for j in 1:t) >= 1) + elseif edge_is_diagonal + edge_is_sw_to_ne = + edge_is_diagonal && sign(r[1] - s[1]) == sign(r[1] - s[2]) + subrect_sw_corner = min.(r, s) + triangle_cuts_se_to_nw = + triangle_direction[subrect_sw_corner...] + if (edge_is_sw_to_ne && triangle_cuts_se_to_nw) || ( + edge_is_diagonal && + !edge_is_sw_to_ne && + !triangle_cuts_se_to_nw + ) + JuMP.@constraint( + sub_model, + sum(z[j, r, s] for j in 1:t) >= 1 + ) + else + JuMP.@constraint( + sub_model, + sum(z[j, r, s] for j in 1:t) == 0 + ) + end + else + @assert r[1] == s[1] || r[2] == s[2] + JuMP.@constraint(sub_model, sum(z[j, r, s] for j in 1:t) == 0) + end + end + + JuMP.@objective(sub_model, Min, sum(x) + sum(y)) + JuMP.unset_silent(sub_model) + JuMP.optimize!(sub_model) + if JuMP.primal_status(sub_model) == JuMP.FEASIBLE_POINT + x_val = JuMP.value.(x) + y_val = JuMP.value.(y) + break + else + t += 1 + @show t + end + end + z = JuMP.@variable(model, [1:t], Bin, base_name = "z_$counter") + + for k in 1:t + JuMP.@constraints( + model, + begin + sum(λ[i, j] for (i, j) in J if x_val[k, (i, j)] ≈ 1) ≤ z[k] + sum(λ[i, j] for (i, j) in J if y_val[k, (i, j)] ≈ 1) ≤ 1 - z[k] + end + ) + end + return nothing +end diff --git a/src/methods/bivariate/optimal_triangle_selection.jl b/src/methods/bivariate/optimal_triangle_selection.jl new file mode 100644 index 0000000..e1ee0b6 --- /dev/null +++ b/src/methods/bivariate/optimal_triangle_selection.jl @@ -0,0 +1,117 @@ +# Copyright (c) 2016: Joey Huchette and contributors +# +# 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. + +# TODO: Generalize to multivariate case. +struct OptimalTriangleSelection <: Method + sub_solver::Any + axis_method::Method +end +function OptimalTriangleSelection(sub_solver) + return OptimalTriangleSelection(sub_solver, Logarithmic()) +end + +axis_method(method::OptimalTriangleSelection) = method.axis_method + +function formulate_triangle_selection!( + model::JuMP.Model, + λ::Matrix{JuMP.VariableRef}, + triangle_direction::Matrix{Bool}, + method::OptimalTriangleSelection, + _::GridTriangulation, +) + counter = model.ext[:PWL].counter + n_1, n_2 = size(λ) + J = Set((i, j) for i in 1:n_1, j in 1:n_2) + x_val = Matrix{Float64}(undef, n_1, n_2) + y_val = Matrix{Float64}(undef, n_1, n_2) + + t = 1 + while true + @show method.sub_solver + sub_model = JuMP.Model(method.sub_solver) + JuMP.@variable(sub_model, x[1:t, J], Bin) + JuMP.@variable(sub_model, y[1:t, J], Bin) + JuMP.@variable(sub_model, z[1:t, J, J], Bin) + for j in 1:t + for r in J, s in J + # lexicographic ordering on points on grid + if r[1] > s[1] || (r[1] == s[1] && r[2] ≥ s[2]) + continue + end + JuMP.@constraints( + sub_model, + begin + z[j, r, s] ≤ x[j, r] + x[j, s] + z[j, r, s] ≤ x[j, r] + y[j, r] + z[j, r, s] ≤ x[j, s] + y[j, s] + z[j, r, s] ≤ y[j, r] + y[j, s] + z[j, r, s] ≥ x[j, r] + y[j, s] - 1 + z[j, r, s] ≥ x[j, s] + y[j, r] - 1 + end + ) + end + for r in J + JuMP.@constraint(sub_model, x[j, r] + y[j, r] ≤ 1) + end + end + + for r in J, s in J + # lexicographic ordering on points on grid + if r[1] > s[1] || (r[1] == s[1] && r[2] ≥ s[2]) + continue + end + edge_is_far_away = LinearAlgebra.norm(r .- s, Inf) > 1 + edge_is_diagonal = + !edge_is_far_away && (abs(r[1] - s[1]) == abs(r[2] - s[2]) == 1) + if edge_is_diagonal + edge_is_sw_to_ne = + edge_is_diagonal && sign(r[1] - s[1]) == sign(r[1] - s[2]) + subrect_sw_corner = min.(r, s) + triangle_cuts_se_to_nw = + triangle_direction[subrect_sw_corner...] + if (edge_is_sw_to_ne && triangle_cuts_se_to_nw) || ( + edge_is_diagonal && + !edge_is_sw_to_ne && + !triangle_cuts_se_to_nw + ) + JuMP.@constraint( + sub_model, + sum(z[j, r, s] for j in 1:t) >= 1 + ) + else + JuMP.@constraint( + sub_model, + sum(z[j, r, s] for j in 1:t) == 0 + ) + end + elseif !edge_is_far_away + @assert r[1] == s[1] || r[2] == s[2] + JuMP.@constraint(sub_model, sum(z[j, r, s] for j in 1:t) == 0) + end + end + + JuMP.@objective(sub_model, Min, sum(x) + sum(y)) + JuMP.optimize!(sub_model) + if JuMP.primal_status(sub_model) == JuMP.FEASIBLE_POINT + x_val = JuMP.value.(x) + y_val = JuMP.value.(y) + break + else + t += 1 + end + end + z = JuMP.@variable(model, [1:t], Bin, base_name = "z_$counter") + + for k in 1:t + JuMP.@constraints( + model, + begin + sum(λ[i, j] for (i, j) in J if x_val[k, (i, j)] ≈ 1) ≤ z[k] + sum(λ[i, j] for (i, j) in J if y_val[k, (i, j)] ≈ 1) ≤ 1 - z[k] + end + ) + end + return nothing +end diff --git a/src/methods/bivariate/six_stencil.jl b/src/methods/bivariate/six_stencil.jl new file mode 100644 index 0000000..75c47af --- /dev/null +++ b/src/methods/bivariate/six_stencil.jl @@ -0,0 +1,131 @@ +# Copyright (c) 2016: Joey Huchette and contributors +# +# 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. + +struct SixStencil <: Method + axis_method::Method +end +SixStencil() = SixStencil(Logarithmic()) + +axis_method(method::SixStencil) = method.axis_method + +# TODO: Unit tests for biclique cover +function formulate_triangle_selection!( + model::JuMP.Model, + λ::Matrix{JuMP.VariableRef}, + triangle_direction::Matrix{Bool}, + method::SixStencil, + _::GridTriangulation, +) + n_1, n_2 = size(λ) + @assert size(triangle_direction) == (n_1 - 1, n_2 - 1) + counter = model.ext[:PWL].counter + + # diagonal lines running from SW to NE. Grouped with an offset of 3. + w_sw_ne = JuMP.@variable(model, [0:2], Bin, base_name = "w_sw_ne_$counter") + for o in 0:2 + A_o = Set{Tuple{Int,Int}}() + B_o = Set{Tuple{Int,Int}}() + for off_1 in o:3:(n_1-2) + sw_in_A = true # whether we put the SW corner of the next triangle to cover in set A + for i in (1+off_1):(n_1-1) + j = i - off_1 + if !(1 ≤ i ≤ n_1 - 1) + continue + end + if !(1 ≤ j ≤ n_2 - 1) + continue # should never happen + end + if triangle_direction[i, j] # if we need to cover the edge... + if sw_in_A # figure out which set we need to put it in; this depends on previous triangle in our current line + push!(A_o, (i, j)) + push!(B_o, (i + 1, j + 1)) + else + push!(A_o, (i + 1, j + 1)) + push!(B_o, (i, j)) + end + sw_in_A = !sw_in_A + end + end + end + for off_2 in (3-o):3:(n_2-1) + sw_in_A = true + for j in (off_2+1):(n_2-1) + i = j - off_2 + if !(1 ≤ i ≤ n_1 - 1) + continue + end + if triangle_direction[i, j] + if sw_in_A + push!(A_o, (i, j)) + push!(B_o, (i + 1, j + 1)) + else + push!(A_o, (i + 1, j + 1)) + push!(B_o, (i, j)) + end + sw_in_A = !sw_in_A + end + end + end + JuMP.@constraints( + model, + begin + sum(λ[i, j] for (i, j) in A_o) ≤ w_sw_ne[o] + sum(λ[i, j] for (i, j) in B_o) ≤ 1 - w_sw_ne[o] + end + ) + end + + w_se_nw = JuMP.@variable(model, [0:2], Bin, base_name = "w_se_nw_$counter") + for o in 0:2 + A_o = Set{Tuple{Int,Int}}() + B_o = Set{Tuple{Int,Int}}() + for off_1 in o:3:(n_1-2) + se_in_A = true + for j in 1:(n_2-1) + i = n_1 - j - off_1 + if !(1 ≤ i ≤ n_1 - 1) + continue + end + if !triangle_direction[i, j] + if se_in_A + push!(A_o, (i + 1, j)) + push!(B_o, (i, j + 1)) + else + push!(A_o, (i, j + 1)) + push!(B_o, (i + 1, j)) + end + se_in_A = !se_in_A + end + end + end + for off_2 in (3-o):3:(n_2-1) + se_in_A = true + for j in (off_2+1):(n_2-1) + i = n_1 - j + off_2 + if !(1 ≤ i ≤ n_1 - 1) + continue + end + if !triangle_direction[i, j] + if se_in_A + push!(A_o, (i + 1, j)) + push!(B_o, (i, j + 1)) + else + push!(A_o, (i, j + 1)) + push!(B_o, (i + 1, j)) + end + se_in_A = !se_in_A + end + end + end + JuMP.@constraints( + model, + begin + sum(λ[i, j] for (i, j) in A_o) ≤ w_se_nw[o] + sum(λ[i, j] for (i, j) in B_o) ≤ 1 - w_se_nw[o] + end + ) + end + return nothing +end diff --git a/src/methods/bivariate/union_jack.jl b/src/methods/bivariate/union_jack.jl new file mode 100644 index 0000000..23f3452 --- /dev/null +++ b/src/methods/bivariate/union_jack.jl @@ -0,0 +1,42 @@ +# Copyright (c) 2016: Joey Huchette and contributors +# +# 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. + +struct UnionJack <: Method + axis_method::Method +end +UnionJack() = UnionJack(Logarithmic()) + +axis_method(method::UnionJack) = method.axis_method + +function formulate_triangle_selection!( + model::JuMP.Model, + λ::Matrix{JuMP.VariableRef}, + triangle_direction::Matrix{Bool}, + method::UnionJack, + _::UnionJackTriangulation, +) + n_1, n_2 = size(λ) + @assert size(triangle_direction) == (n_1 - 1, n_2 - 1) + counter = model.ext[:PWL].counter + + # TODO: Certify triangulation is Union Jack + + # j_start = number of triangle segments incident to (1, 1), the lower-left + # most point in the grid. + j_start = triangle_direction[1, 1] ? 1 : 2 + + # diagonal lines running from SW to NE. Grouped with an offset of 3. + z = JuMP.@variable(model, binary = true, base_name = "w_$counter") + + JuMP.@constraints( + model, + begin + sum(λ[i, j] for i in 1:2:n_1, j in j_start:2:n_2) ≤ z + sum(λ[i, j] for i in 2:2:n_1, j in (3-j_start):2:n_2) ≤ 1 - z + end + ) + + return nothing +end diff --git a/src/methods/multivariate/convex_combination.jl b/src/methods/multivariate/convex_combination.jl new file mode 100644 index 0000000..17a588e --- /dev/null +++ b/src/methods/multivariate/convex_combination.jl @@ -0,0 +1,86 @@ +# Copyright (c) 2016: Joey Huchette and contributors +# +# 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. + +struct ConvexCombination <: Method end + +function formulate_pwl!( + model::JuMP.Model, + input_vars::NTuple{D,VarOrAff}, + output_vars::NTuple{F,VarOrAff}, + pwl::PWLFunctionPointRep{D,F}, + method::ConvexCombination, + direction::DIRECTION, +) where {D,F} + # TODO: assert PWL function is continuous + counter = model.ext[:PWL].counter + segments = pwl.segments + z = JuMP.@variable(model, [segments], Bin, base_name = "z_$counter") + JuMP.@constraint(model, sum(z) == 1) + + all_input_vals = Set{NTuple{D,Float64}}() + output_val_map = Dict{NTuple{D,Float64},NTuple{F,Float64}}() + for seg in segments + for (it, input_val) in enumerate(seg.input_vals) + push!(all_input_vals, input_val) + output_val = seg.output_vals[it] + if haskey(output_val_map, input_val) + if output_val_map[input_val] != output_val + error( + "The ConvexCombination method does not currently support discontinuous piecewise linear functions.", + ) + end + else + output_val_map[input_val] = output_val + end + end + end + λ = JuMP.@variable( + model, + [all_input_vals], + lower_bound = 0, + upper_bound = 1, + base_name = "λ_$counter" + ) + JuMP.@constraint(model, sum(λ) == 1) + for j in 1:D + JuMP.@constraint( + model, + input_vars[j] == + sum(λ[input_val] * input_val[j] for input_val in all_input_vals) + ) + end + for j in 1:F + rhs = sum( + λ[input_val] * output_val_map[input_val][j] for + input_val in all_input_vals + ) + _constrain_output_var(model, output_vars[j], rhs, direction) + end + for input_val in all_input_vals + JuMP.@constraint( + model, + λ[input_val] ≤ + sum(z[seg] for seg in segments if input_val in seg.input_vals) + ) + end + return nothing +end + +function formulate_sos2!( + model::JuMP.Model, + λ::Vector{T}, + method::ConvexCombination, +) where {T<:VarOrAff} + counter = model.ext[:PWL].counter + n = length(λ) + z = JuMP.@variable(model, [1:n-1], Bin, base_name = "z_$counter") + JuMP.@constraint(model, sum(z) == 1) + JuMP.@constraint(model, λ[1] ≤ z[1]) + for i in 2:n-1 + JuMP.@constraint(model, λ[i] ≤ z[i-1] + z[i]) + end + JuMP.@constraint(model, λ[n] ≤ z[n-1]) + return nothing +end diff --git a/src/methods/multivariate/disaggregated_logarithmic.jl b/src/methods/multivariate/disaggregated_logarithmic.jl new file mode 100644 index 0000000..cdf5c21 --- /dev/null +++ b/src/methods/multivariate/disaggregated_logarithmic.jl @@ -0,0 +1,63 @@ +# Copyright (c) 2016: Joey Huchette and contributors +# +# 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. + +struct DisaggregatedLogarithmic <: Method end + +function formulate_pwl!( + model::JuMP.Model, + input_vars::NTuple{D,VarOrAff}, + output_vars::NTuple{F,VarOrAff}, + pwl::PWLFunctionPointRep{D,F}, + method::DisaggregatedLogarithmic, + direction::DIRECTION, +) where {D,F} + counter = model.ext[:PWL].counter + + segments = pwl.segments + num_bps = Dict(seg => length(seg.input_vals) for seg in segments) + + γ = JuMP.@variable( + model, + [seg in segments, i in 1:num_bps[seg]], + lower_bound = 0, + upper_bound = 1, + base_name = "γ_$counter" + ) + JuMP.@constraint(model, sum(γ) == 1) + + for j in 1:D + JuMP.@constraint( + model, + input_vars[j] == sum( + sum(γ[seg, i] * seg.input_vals[i][j] for i in 1:num_bps[seg]) + for seg in segments + ) + ) + end + for j in 1:F + rhs = sum( + sum(γ[seg, i] * seg.output_vals[i][j] for i in 1:num_bps[seg]) + for seg in segments + ) + _constrain_output_var(model, output_vars[j], rhs, direction) + end + + r = ceil(Int, log2(length(segments))) + if r == 0 + return nothing + end + _H = _reflected_gray_codes(r) + H = Dict(segments[i] => _H[i] for i in 1:length(segments)) + z = JuMP.@variable(model, [1:r], Bin, base_name = "z_$counter") + for j in 1:r + JuMP.@constraint( + model, + sum( + sum(γ[seg, i] * H[seg][j] for i in 1:num_bps[seg]) for + seg in segments + ) == z[j] + ) + end +end diff --git a/src/methods/multivariate/multiple_choice.jl b/src/methods/multivariate/multiple_choice.jl new file mode 100644 index 0000000..3855dc8 --- /dev/null +++ b/src/methods/multivariate/multiple_choice.jl @@ -0,0 +1,45 @@ +# Copyright (c) 2016: Joey Huchette and contributors +# +# 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. + +struct MultipleChoice <: Method end + +function formulate_pwl!( + model::JuMP.Model, + input_vars::Tuple{VarOrAff}, + output_vars::NTuple{F,VarOrAff}, + pwl::PWLFunctionHyperplaneRep{D,F}, + method::MultipleChoice, + direction::DIRECTION, +) where {D,F} + x_hat = JuMP.@variable(model, [segments, 1:D], base_name = "x_hat_$counter") + y_hat = JuMP.@variable(model, [segments, 1:F], base_name = "y_hat_$counter") + z = JuMP.@variable(model, [segments], Bin, base_name = "z_$counter") + JuMP.@constraint(model, sum(z) == 1) + for i in 1:D + JuMP.@constraint(model, sum(x_hat[:, i]) == x[i]) + end + for i in 1:F + JuMP.@constraint(model, sum(y_hat[:, i]) == y[i]) + end + for seg in segments + for constraint in seg.constraints + coeffs, offset = constraint.coeffs, constraint.offset + JuMP.@constraint( + model, + LinearAlgebra.dot(coeffs, x_hat[seg, :]) + offset * z[seg] ≥ 0 + ) + end + for i in 1:F + output_func = seg.output_funcs[i] + coeffs, offset = output_func.coeffs, output_func.offset + JuMP.@constraint( + model, + y_hat[seg, i] == + LinearAlgebra.dot(coeffs, x_hat[seg, :]) + offset * z[seg] + ) + end + end + return nothing +end diff --git a/src/methods/univariate/incremental.jl b/src/methods/univariate/incremental.jl new file mode 100644 index 0000000..96db5b3 --- /dev/null +++ b/src/methods/univariate/incremental.jl @@ -0,0 +1,46 @@ +# Copyright (c) 2016: Joey Huchette and contributors +# +# 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. + +# TODO: Implement bivariate version of the incremental formulation +struct Incremental <: Method end + +function formulate_pwl!( + model::JuMP.Model, + input_vars::Tuple{VarOrAff}, + output_vars::NTuple{F,VarOrAff}, + pwl::PWLFunctionPointRep{1,F}, + method::Incremental, + direction::DIRECTION, +) where {F} + grid = _continuous_gridpoints_or_die(pwl) + xs, ys = grid.input_vals, grid.output_vals + + counter = model.ext[:PWL].counter + n = length(pwl.segments) + 1 + @assert length(xs) == length(ys) == n + + δ = JuMP.@variable( + model, + [1:n], + lower_bound = 0, + upper_bound = 1, + base_name = "δ_$counter" + ) + z = JuMP.@variable(model, [1:n-1], Bin, base_name = "z_$counter") + JuMP.@constraint( + model, + input_vars[1] == + xs[1][1] + sum(δ[i] * (xs[i+1][1] - xs[i][1]) for i in 1:n-1) + ) + for j in 1:F + rhs = ys[1][j] + sum(δ[i] * (ys[i+1][j] - ys[i][j]) for i in 1:n-1) + _constrain_output_var(model, output_vars[j], rhs, direction) + end + for i in 1:n-1 + JuMP.@constraint(model, δ[i+1] ≤ z[i]) + JuMP.@constraint(model, z[i] ≤ δ[i]) + end + return +end diff --git a/src/methods/univariate/logarithmic_embedding.jl b/src/methods/univariate/logarithmic_embedding.jl new file mode 100644 index 0000000..1c3b357 --- /dev/null +++ b/src/methods/univariate/logarithmic_embedding.jl @@ -0,0 +1,32 @@ +# Copyright (c) 2016: Joey Huchette and contributors +# +# 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. + +struct LogarithmicEmbedding <: Method end + +function formulate_sos2!( + model::JuMP.Model, + λ::Vector{T}, + method::LogarithmicEmbedding, +) where {T<:VarOrAff} + counter = model.ext[:PWL].counter + n = length(λ) + d = n - 1 + if 0 <= d <= 1 + return nothing + end + k = ceil(Int, log2(d)) + if k == 0 + return nothing + end + y = JuMP.@variable(model, [1:k], Bin, base_name = "y_$counter") + _sos2_encoding_constraints!( + model, + λ, + y, + _reflected_gray_codes(k), + _unit_vector_hyperplanes(k), + ) + return nothing +end diff --git a/src/methods/univariate/logarithmic_independent_branching.jl b/src/methods/univariate/logarithmic_independent_branching.jl new file mode 100644 index 0000000..4c549bd --- /dev/null +++ b/src/methods/univariate/logarithmic_independent_branching.jl @@ -0,0 +1,38 @@ +# Copyright (c) 2016: Joey Huchette and contributors +# +# 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. + +struct LogarithmicIndependentBranching <: Method end + +function formulate_sos2!( + model::JuMP.Model, + λ::Vector{T}, + method::LogarithmicIndependentBranching, +) where {T<:VarOrAff} + counter = model.ext[:PWL].counter + n = length(λ) + d = n - 1 + if 0 <= d <= 1 + return nothing + end + k = ceil(Int, log2(d)) + if k == 0 + return nothing + end + z = JuMP.@variable(model, [1:k], Bin, base_name = "y_$counter") + _H = _reflected_gray_codes(k) + H = Dict(i => _H[i] for i in 1:d) + H[0] = H[1] + H[d+1] = H[d] + for j in 1:k + JuMP.@constraints( + model, + begin + sum(λ[i] for i in 1:n if H[i-1][j] == H[i][j] == 1) ≤ z[j] + sum(λ[i] for i in 1:n if H[i-1][j] == H[i][j] == 0) ≤ 1 - z[j] + end + ) + end + return nothing +end diff --git a/src/methods/univariate/native_sos2.jl b/src/methods/univariate/native_sos2.jl new file mode 100644 index 0000000..36d56d3 --- /dev/null +++ b/src/methods/univariate/native_sos2.jl @@ -0,0 +1,15 @@ +# Copyright (c) 2016: Joey Huchette and contributors +# +# 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. + +struct NativeSOS2 <: Method end + +function formulate_sos2!( + model::JuMP.Model, + λ::Vector{T}, + method::NativeSOS2, +) where {T<:VarOrAff} + JuMP.@constraint(model, λ in JuMP.SOS2([k for k in 1:length(λ)])) + return nothing +end diff --git a/src/methods/univariate/sos2_formulation_base.jl b/src/methods/univariate/sos2_formulation_base.jl new file mode 100644 index 0000000..cb63b87 --- /dev/null +++ b/src/methods/univariate/sos2_formulation_base.jl @@ -0,0 +1,168 @@ +# Copyright (c) 2016: Joey Huchette and contributors +# +# 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. + +function _sos2_encoding_constraints!( + m::JuMP.Model, + λ::Vector{T}, + y::Vector{JuMP.VariableRef}, + h::Vector{Vector{Float64}}, + B::Vector{Vector{Float64}}, +) where {T<:VarOrAff} + n = length(λ) - 1 + for b in B + JuMP.@constraints( + m, + begin + 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 +end + +function _reflected_gray_codes(k::Int)::Vector{Vector{Float64}} + if k <= 0 + error("Invalid code length $k") + elseif k == 1 + return [[0.0], [1.0]] + else + codes = _reflected_gray_codes(k - 1) + return vcat( + [vcat(code, 0.0) for code in codes], + [vcat(code, 1.0) for code in reverse(codes)], + ) + end +end + +function _zigzag_codes(k::Int)::Vector{Vector{Float64}} + if k <= 0 + error("Invalid code length $k") + elseif k == 1 + return [[0.0], [1.0]] + else + codes = _zigzag_codes(k - 1) + return vcat( + [vcat(code, 0.0) for code in codes], + [vcat(code, 1.0) for code in codes], + ) + end +end + +function _integer_zigzag_codes(k::Int)::Vector{Vector{Float64}} + if k <= 0 + error("Invalid code length $k") + elseif k == 1 + return [[0.0], [1.0]] + else + codes = _integer_zigzag_codes(k - 1) + offset = [2^(j - 2) for j in k:-1:2] + return vcat( + [vcat(code, 0.0) for code in codes], + [vcat(code .+ offset, 1.0) for code in codes], + ) + end +end + +function _zigzag_hyperplanes(k::Int)::Vector{Vector{Float64}} + hps = Vector{Int}[] + for i in 1:k + hp = zeros(Int, k) + hp[i] = 1 + for j in (i+1):k + hp[j] = 2^(j - i - 1) + end + push!(hps, hp) + end + return hps +end + +function _unit_vector_hyperplanes(k::Int)::Vector{Vector{Float64}} + hps = Vector{Float64}[] + for i in 1:k + hp = zeros(Int, k) + hp[i] = 1 + push!(hps, hp) + end + return hps +end + +# ConvexCombination could also be added to this list, but the multivariate +# implementation works just fine. +const SOS2Method = Union{ + LogarithmicEmbedding, + LogarithmicIndependentBranching, + NativeSOS2, + ZigZagBinary, + ZigZagInteger, +} + +function formulate_pwl!( + model::JuMP.Model, + input_vars::Tuple{VarOrAff}, + output_vars::NTuple{F,VarOrAff}, + pwl::PWLFunctionPointRep{1,F}, + method::SOS2Method, + direction::DIRECTION, +) where {F} + grid = _continuous_gridpoints_or_die(pwl) + λ = _create_convex_multiplier_vars( + model, + grid, + input_vars, + output_vars, + direction, + ) + formulate_sos2!(model, λ, method) + return nothing +end + +function formulate_sos2!( + model::JuMP.Model, + λ::Vector{T}, + method::Method, +) where {T<:VarOrAff} + n = length(λ) + D = 1 + F = n + d = range(0, 1; length = n) + segments = SegmentPointRep{D,F}[] + for i in 1:(n-1) + output_left = ntuple(t -> t == i ? 1.0 : 0.0, n) + output_right = ntuple(t -> t == i + 1 ? 1.0 : 0.0, n) + push!( + segments, + SegmentPointRep{D,F}( + [(d[i],), (d[i+1],)], + [output_left, output_right], + ), + ) + end + dummy_input_var = JuMP.@variable(model, lower_bound = 0, upper_bound = 1) + dummy_pwl = PWLFunctionPointRep{1,F}(segments, Intervals()) + formulate_pwl!( + model, + (dummy_input_var,), + tuple(λ...), + dummy_pwl, + method, + Graph, + ) + return nothing +end diff --git a/src/methods/univariate/zig_zag_binary.jl b/src/methods/univariate/zig_zag_binary.jl new file mode 100644 index 0000000..1405185 --- /dev/null +++ b/src/methods/univariate/zig_zag_binary.jl @@ -0,0 +1,32 @@ +# Copyright (c) 2016: Joey Huchette and contributors +# +# 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. + +struct ZigZagBinary <: Method end + +function formulate_sos2!( + model::JuMP.Model, + λ::Vector{T}, + method::ZigZagBinary, +) where {T<:VarOrAff} + counter = model.ext[:PWL].counter + n = length(λ) + d = n - 1 + if 0 <= d <= 1 + return nothing + end + k = ceil(Int, log2(d)) + if k == 0 + return nothing + end + y = JuMP.@variable(model, [1:k], Bin, base_name = "y_$counter") + _sos2_encoding_constraints!( + model, + λ, + y, + _zigzag_codes(k), + _zigzag_hyperplanes(k), + ) + return nothing +end diff --git a/src/methods/univariate/zig_zag_integer.jl b/src/methods/univariate/zig_zag_integer.jl new file mode 100644 index 0000000..7e9c6f3 --- /dev/null +++ b/src/methods/univariate/zig_zag_integer.jl @@ -0,0 +1,39 @@ +# Copyright (c) 2016: Joey Huchette and contributors +# +# 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. + +struct ZigZagInteger <: Method end + +function formulate_sos2!( + model::JuMP.Model, + λ::Vector{T}, + method::ZigZagInteger, +) where {T<:VarOrAff} + counter = model.ext[:PWL].counter + n = length(λ) + d = n - 1 + if 0 <= d <= 1 + return nothing + end + k = ceil(Int, log2(d)) + if k == 0 + return nothing + end + y = JuMP.@variable( + model, + [i in 1:k], + Int, + lower_bound = 0, + upper_bound = 2^(k - i), + base_name = "y_$counter" + ) + _sos2_encoding_constraints!( + model, + λ, + y, + _integer_zigzag_codes(k), + _unit_vector_hyperplanes(k), + ) + return nothing +end diff --git a/src/methods/util.jl b/src/methods/util.jl new file mode 100644 index 0000000..3326ad8 --- /dev/null +++ b/src/methods/util.jl @@ -0,0 +1,290 @@ +# Copyright (c) 2016: Joey Huchette and contributors +# +# 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. + +function _constrain_output_var( + model::JuMP.Model, + output_var::VarOrAff, + func_val::VarOrAff, + direction::DIRECTION, +) + if direction == Graph + JuMP.@constraint(model, output_var == func_val) + elseif direction == Hypograph + JuMP.@constraint(model, output_var <= func_val) + else + @assert direction == Epigraph + JuMP.@constraint(model, output_var >= func_val) + end + return +end + +struct Grid{D,F} + input_vals::Array{NTuple{D,Float64},D} + output_vals::Array{NTuple{F,Float64},D} + + function Grid( + input_vals::Array{NTuple{D,Float64},D}, + output_vals::Array{NTuple{F,Float64},D}, + ) where {D,F} + if size(input_vals) != size(output_vals) + error( + "Incompatible input/output sizes $(size(input_vals)) and $(size(output_vals)). Must be the same.", + ) + end + return new{D,F}(input_vals, output_vals) + end +end + +function _continuous_gridpoints_or_die( + pwl::PWLFunction{D,F,SegmentPointRep{D,F}}, +) where {D,F} + # Step 1: Collect all input points + xs = Set{NTuple{D,Float64}}() + for segment in pwl.segments + for val in segment.input_vals + push!(xs, val) + end + end + + # Step 2: Verify that function is not multivalued + ys = Dict{NTuple{D,Float64},NTuple{F,Float64}}() + for segment in pwl.segments + for (x, y) in zip(segment.input_vals, segment.output_vals) + if haskey(ys, x) + if y != ys[x] + error( + "Expected piecewise linear function to be continuous.", + ) + end + else + ys[x] = y + end + end + end + + @assert length(xs) == length(ys) + + # Step 3: Build grid, and verify every point is included in PWL function + axes = [sort!(unique(x[i] for x in xs)) for i in 1:D] + + x_grid = collect(Base.Iterators.product(axes...)) + y_grid = map(x -> get(ys, x, nothing), x_grid) + + if length(x_grid) != length(xs) + error("Decomposition is not aligned along a grid.") + end + + # Step 4: Verify that the domain is equal to the convex hull of the gridpoints + # TODO: Implement this. + if D == 1 + if length(pwl.segments) != length(x_grid) - 1 + error( + "Univariate piecewise linear function does not have a connected domain.", + ) + end + elseif D == 2 + if length(pwl.segments) != + 2 * (size(x_grid, 1) - 1) * (size(x_grid, 2) - 1) + error( + "Bivariate piecewise linear function domain is not rectangular.", + ) + end + else + Base.warn( + "Cannot currently verify that multivariate piecewise linear functions have domains equal to entire grid of breakpoints.", + ) + end + + return Grid(x_grid, y_grid) +end + +function _canonicalize_triangulation( + pwl::PWLFunction{D,F,SegmentPointRep{D,F}}, + grid::Grid{D,F}, +) where {D,F} + _check_triangulation(pwl) + + xs = grid.input_vals + U = [unique(x[j] for x in xs) for j in 1:D] + for j in 1:D + @assert issorted(U[j]) + end + x_to_i = [Dict(U[j][i] => i for i in 1:size(xs, j)) for j in 1:D] + + canonical_input_segments = Vector{NTuple{D,Int64}}[] + for segment in pwl.segments + push!( + canonical_input_segments, + [ntuple(j -> x_to_i[j][v[j]], D) for v in segment.input_vals], + ) + end + return canonical_input_segments +end + +function _create_convex_multiplier_vars( + model::JuMP.Model, + grid::Grid{D,F}, + input_vars::NTuple{D,JuMP.VariableRef}, + output_vars::NTuple{F,JuMP.VariableRef}, + direction::DIRECTION, +) where {D,F} + counter = model.ext[:PWL].counter + + # TODO: Name λ variables + λ = similar(Array{JuMP.VariableRef}, axes(grid.input_vals)) + for I in eachindex(λ) + λ[I] = JuMP.@variable(model, lower_bound = 0, upper_bound = 1) + end + + JuMP.@constraint(model, sum(λ) == 1) + for j in 1:D + JuMP.@constraint( + model, + sum(λ[I] * grid.input_vals[I][j] for I in eachindex(λ)) == + input_vars[j] + ) + end + for j in 1:F + _constrain_output_var( + model, + output_vars[j], + sum(λ[I] * grid.output_vals[I][j] for I in eachindex(λ)), + direction, + ) + end + return λ +end + +function _check_triangle(segment::SegmentPointRep{D}) where {D} + if !(length(segment.input_vals) == length(segment.output_vals) == D + 1) + error( + "Encountered something that was not a simplex while expecting a trianglulation.", + ) + end + return +end + +function _check_triangle(segment::SegmentHyperplaneRep{D}) where {D} + if length(segment.breakpoints) != D + 1 + error( + "Encountered something that was not a simplex while expecting a trianglulation.", + ) + end + return +end + +function _check_triangulation(pwl::PWLFunction) where {D,F} + # TODO: Add check that segments lie in a single subrectangle + for segment in pwl.segments + _check_triangle(segment) + end + return +end + +function _generalized_celaya_hyperplanes(k::Int) + return _compute_hyperplanes(_generalized_celaya_codes(k)) +end + +function _symmetric_celaya_hyperplanes(k::Int) + return _compute_hyperplanes(_symmetric_celaya_codes(k)) +end + +# TODO: Figure out assumptions at play here +function _compute_hyperplanes(C::Vector{Vector{T}}) where {T<:Number} + n = length(C) + k = length(C[1]) + + d = zeros(n - 1, k) + for i in 1:n-1 + d[i, :] = C[i+1] - C[i] + end + + if k <= 1 + error("Cannot process codes of length $k") + elseif k == 2 + @assert n == 4 + spanners = Vector{Float64}[] + for i in 1:n-1 + v = canonical!(d[i, :]) + v = [v[2], -v[1]] + push!(spanners, v) + end + indices = [1] + approx12 = isapprox(spanners[1], spanners[2]) + approx13 = isapprox(spanners[1], spanners[3]) + approx23 = isapprox(spanners[2], spanners[3]) + if !approx12 + push!(indices, 2) + end + if !approx13 && !(!approx12 && approx23) + push!(indices, 3) + end + return spanners[indices] + end + + indices = [1, 2] + spanners = Vector{Float64}[] + while !isempty(indices) + if indices == [n - 1] + break + end + if LinearAlgebra.rank(d[indices, :]) == length(indices) && + length(indices) <= k - 1 + if length(indices) == k - 1 + nullsp = LinearAlgebra.nullspace(d[indices, :]) + @assert size(nullsp, 2) == 1 + v = vec(nullsp) + push!(spanners, canonical!(v)) + end + if indices[end] != n - 1 + push!(indices, indices[end] + 1) + else + pop!(indices) + indices[end] += 1 + end + else + if indices[end] != n - 1 + indices[end] += 1 + else + pop!(indices) + indices[end] += 1 + end + end + end + + keepers = [1] + for i in 2:length(spanners) + alreadyin = false + for j in keepers + if isapprox(spanners[i], spanners[j]) + alreadyin = true + break + end + end + if !alreadyin + push!(keepers, i) + end + end + + return spanners[keepers] +end + +function _canonical!(v::Vector{Float64}) + normalize!(v) + for j in 1:length(v) + if abs(v[j]) < 1e-8 + v[j] = 0 + end + end + sgn = sign(v[findfirst([x != 0 for x in v])]) + for j in 1:length(v) + if abs(v[j]) < 1e-8 + v[j] = 0 + else + v[j] *= sgn + end + end + return v +end diff --git a/src/pwlinear.jl b/src/pwlinear.jl new file mode 100644 index 0000000..1e7d8c1 --- /dev/null +++ b/src/pwlinear.jl @@ -0,0 +1,169 @@ +# Copyright (c) 2016: Joey Huchette and contributors +# +# 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. + +function formulate_pwl!( + model::JuMP.Model, + input_vals::Vector{NTuple{D,VarOrAff}}, + output_vals::Vector{NTuple{F,VarOrAff}}, + pwl::PWLFunction, + method::Method, + direction::DIRECTION, +) where {D,F} + return error( + "No support for a R^$D -> R^$F piecewise linear function using the $method method.", + ) +end + +_default_method(::Val{1}) = Logarithmic() +_default_method(::Val{2}) = SixStencil() +# _default_method(::Val) = MultipleChoice() + +""" + piecewiselinear(model, input_vars, pwl::PWLFunction; method, direction, output_vars) + + piecewiselinear(model, input_var, x, f::Function; method, direction, output_var) + piecewiselinear(model, input_var, x, fd; method, direction, output_var) + + piecewiselinear(model, input_var_x, input_var_y, x, y, f::Function; method, direction, output_var, pattern) + + +Formulates a piecewise linear function in the given JuMP model. + +# Arguments +- `model::JuMP.Model`: The JuMP model to which the piecewise linear function will be added. +- `input_vars::NTuple{D,VarOrAff}`: A tuple of input variables or affine expressions. +- `pwl::PWLFunction{D,F,SegmentPointRep{D,F}}`: The piecewise linear function to be added. +- `method::Method`: The method to be used for formulating the piecewise linear function. + Defaults to `Logarithmic()` for 1D and `SixStencil()` for 2D. +- `direction::DIRECTION`: The direction of the piecewise linear function. Defaults to `Graph`. +- `output_vars::Union{Nothing,NTuple{F,VarOrAff}}`: A tuple of output variables or affine expressions. + If `nothing`, new variables will be created. + +# Returns +- `output_vars::NTuple{F,VarOrAff}`: The output variables of the piecewise linear function. + +In addition to the general constructor, there are specialized constructors for univariate +and bivariate problems. + +# Example +```julia +using JuMP +model = Model() +@variable(model, x >= 0) +@variable(model, y >= 0) + +pwl = UnivariatePWLFunction(0:0.1:1, xi -> xi^2) + +# General constructor +output_vars = piecewiselinear(model, (x,), pwl) + +# Specialized constructor for univariate problems +output_var = piecewiselinear(model, x, 0:0.1:1, xi -> xi^2) + +# Specialized constructor for bivariate problems +output_var = piecewiselinear(model, x, y, 0:0.1:1, 0:0.1:1, (xi, yi) -> xi^2 + yi^2; pattern = :BestFit) +``` +""" +function piecewiselinear( + model::JuMP.Model, + input_vars::NTuple{D,VarOrAff}, + pwl::PWLFunction{D,F,SegmentPointRep{D,F}}; + method::Method = _default_method(Val(D)), + direction::DIRECTION = Graph, + output_vars::Union{Nothing,NTuple{F,VarOrAff}} = nothing, +) where {D,F} + initPWL!(model) + counter = model.ext[:PWL].counter + counter += 1 + model.ext[:PWL].counter = counter + + if isempty(pwl.segments) + error( + "I don't know how to handle a piecewise linear function with no breakpoints.", + ) + end + + output_lb = + minimum(minimum(segment.output_vals) for segment in pwl.segments) + output_ub = + maximum(maximum(segment.output_vals) for segment in pwl.segments) + + if output_vars === nothing + output_vars = tuple( + JuMP.@variable( + model, + [i in 1:F], + lower_bound = output_lb[i], + upper_bound = output_ub[i], + base_name = "y_$counter" + )..., + ) + end + + formulate_pwl!(model, input_vars, output_vars, pwl, method, direction) + return output_vars +end + +# Specialization for univariate problems +function piecewiselinear( + model::JuMP.Model, + input_var::VarOrAff, + d, + f::Function; + method::Method = _default_method(Val(1)), + direction::DIRECTION = Graph, + output_var::Union{Nothing,VarOrAff} = nothing, +) + return piecewiselinear( + model, + (input_var,), + UnivariatePWLFunction(d, f); + method = method, + direction = direction, + output_vars = isnothing(output_var) ? nothing : (output_var,), + )[1] +end + +function piecewiselinear( + model::JuMP.Model, + input_var::VarOrAff, + d, + fd; + method::Method = _default_method(Val(1)), + direction::DIRECTION = Graph, + output_var::Union{Nothing,VarOrAff} = nothing, +) + return piecewiselinear( + model, + (input_var,), + UnivariatePWLFunction(d, fd); + method = method, + direction = direction, + output_vars = isnothing(output_var) ? nothing : (output_var,), + )[1] +end + +# Specialization for bivariate problems +function piecewiselinear( + model::JuMP.Model, + input_var_x::VarOrAff, + input_var_y::VarOrAff, + x, + y, + f::Function; + method::Method = _default_method(Val(2)), + direction::DIRECTION = Graph, + output_var::Union{Nothing,VarOrAff} = nothing, + pattern = :K1, +) + return piecewiselinear( + model, + (input_var_x, input_var_y), + BivariatePWLFunction(x, y, f; pattern = pattern); + method = method, + direction = direction, + output_vars = isnothing(output_var) ? nothing : (output_var,), + )[1] +end diff --git a/src/types.jl b/src/types.jl index 5ae0257..fdcced5 100644 --- a/src/types.jl +++ b/src/types.jl @@ -3,131 +3,148 @@ # 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. -struct PWLFunction{D} - x::Vector{NTuple{D,Float64}} - z::Vector{Float64} - T::Vector{Vector{Int}} - meta::Dict -end +abstract type Method end +abstract type UnivariateMethod <: Method end + +@enum DIRECTION Graph Epigraph Hypograph + +# TODO: Make eltypes of input_vals and output_vals a type parameter +abstract type Segment{D,F} end + +struct SegmentPointRep{D,F} <: Segment{D,F} + input_vals::Vector{NTuple{D,Float64}} + output_vals::Vector{NTuple{F,Float64}} -function PWLFunction{D}( - x::Vector{NTuple{D}}, - z::Vector, - T::Vector{Vector}, - meta::Dict, -) where {D} - @assert length(x) == length(z) - for t in T - @assert minimum(t) > 0 && maximum(t) <= length(x) + function SegmentPointRep{D,F}( + input_vals::Vector{NTuple{D,Float64}}, + output_vals::Vector{NTuple{F,Float64}}, + ) where {D,F} + if length(input_vals) != length(output_vals) + error("Must specify the same number of input and output values.") + end + # TODO: Run verifier to ensure this is actually a PWL function + return new{D,F}(input_vals, output_vals) end - return PWLFunction{D}(x, z, T, meta) end -PWLFunction(x, z, T) = PWLFunction(x, z, T, Dict()) +struct AffineFunction{D} + coeffs::NTuple{D,Float64} + offset::Float64 +end + +struct SegmentHyperplaneRep{D,F} <: Segment{D,F} + # Domain given by f_i(x) >= 0 where f_i is i-th constraint in constraints + constraints::Vector{AffineFunction{D}} + funcs::NTuple{F,AffineFunction{D}} +end + +abstract type SegmentStructure{D} end + +struct Intervals <: SegmentStructure{1} end + +abstract type GridTriangulation <: SegmentStructure{2} end +struct UnstructuredTriangulation <: GridTriangulation end +struct K1Triangulation <: GridTriangulation end +struct UnionJackTriangulation <: GridTriangulation end + +struct PWLFunction{D,F,T<:Segment{D,F}} + segments::Vector{T} + structure::SegmentStructure{D} +end + +const PWLFunctionPointRep{D,F} = PWLFunction{D,F,SegmentPointRep{D,F}} +const PWLFunctionHyperplaneRep{D,F} = PWLFunction{D,F,SegmentHyperplaneRep{D,F}} -const UnivariatePWLFunction = PWLFunction{1} +#const UnivariatePWLFunction{F} = PWLFunctionPointRep{1, F} +#const BivariatePWLFunction{F} = PWLFunctionPointRep{2, F} -function UnivariatePWLFunction(x, z) - @assert issorted(x) - return PWLFunction( - Tuple{Float64}[(xx,) for xx in x], - convert(Vector{Float64}, z), - [[i, i + 1] for i in 1:length(x)-1], - ) +const UnivariatePWLFunction = PWLFunctionPointRep{1,1} +const BivariatePWLFunction = PWLFunctionPointRep{2,1} + +function PWLFunctionPointRep{1,1}(x::Vector, z::Vector) + if length(x) != length(z) + error("Mismatch in the number of points and function values") + end + xs = [convert(Float64, xi) for xi in x] + zs = [convert(Float64, zi) for zi in z] + segments = [ + PiecewiseLinearOpt.SegmentPointRep{1,1}( + [(xs[i],), (xs[i+1],)], + [(zs[i],), (zs[i+1],)], + ) for i in 1:length(x)-1 + ] + + return PWLFunction(segments, Intervals()) end -function UnivariatePWLFunction(x, fz::Function) - @assert issorted(x) - return PWLFunction( - Tuple{Float64}[(xx,) for xx in x], - map(t -> convert(Float64, fz(t)), x), - [[i, i + 1] for i in 1:length(x)-1], - ) +function PWLFunctionPointRep{1,1}(x, f::Function) + d = collect(x) + fd = [f(xi) for xi in d] + return PWLFunctionPointRep{1,1}(d, fd) end -const BivariatePWLFunction = PWLFunction{2} +function PWLFunctionPointRep{1,1}(x, z) + return PWLFunctionPointRep{1,1}(collect(x), collect(z)) +end -function BivariatePWLFunction( +function PWLFunctionPointRep{2,1}( x, y, fz::Function; - pattern = :BestFit, + pattern = :K1, seed = hash((length(x), length(y))), ) - @assert issorted(x) - @assert issorted(y) - X = vec(collect(Base.product(x, y))) - # X = vec(Tuple{Float64,Float64}[(_x,_y) for _x in x, _y in y]) - Z = map(t -> convert(Float64, fz(t...)), X) - T = Vector{Vector{Int}}() - m = length(x) - n = length(y) + xs = [convert(Float64, xi) for xi in x] + ys = [convert(Float64, yi) for yi in y] + + segments = SegmentPointRep{2,1}[] + structure = UnstructuredTriangulation() + if pattern == :K1 + structure = K1Triangulation() + elseif pattern == :UnionJack + structure = UnionJackTriangulation() + end + 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] - xL, xU, yL, yU = x[i], x[i+1], y[j], y[j+1] - @assert xL == X[SWt][1] == X[NWt][1] - @assert xU == X[SEt][1] == X[NEt][1] - @assert yL == X[SWt][2] == X[SEt][2] - @assert yU == X[NWt][2] == X[NEt][2] - SW, NW, NE, SE = Z[SWt], Z[NWt], Z[NEt], Z[SEt] - mid1 = 0.5 * (SW + NE) - mid2 = 0.5 * (NW + SE) + for i in 1:length(xs)-1, j in 1:length(ys)-1 + xL, xU, yL, yU = xs[i], xs[i+1], ys[j], ys[j+1] + mid1 = 0.5 * (fz(xL, yL) + fz(xU, yU)) + mid2 = 0.5 * (fz(xL, yU) + fz(xU, yL)) + mid3 = fz(0.5 * (xL + xU), 0.5 * (yL + yU)) + diagonal_nw_se = true if pattern == :Upper - if mid1 > mid2 - t1 = [SWt, NWt, NEt] - t2 = [SWt, NEt, SEt] - else - t1 = [SWt, NWt, SEt] - t2 = [SEt, NWt, NEt] - end + diagonal_nw_se = (mid1 > mid2) elseif pattern == :Lower - if mid1 > mid2 - t1 = [SWt, NWt, SEt] - t2 = [SEt, NWt, NEt] - else - t1 = [SWt, NWt, NEt] - t2 = [SWt, NEt, SEt] - end + diagonal_nw_se = (mid1 < mid2) elseif pattern == :BestFit - mid3 = fz(0.5 * (xL + xU), 0.5 * (yL + yU)) - if abs(mid1 - mid3) < abs(mid2 - mid3) - t1 = [SWt, NWt, NEt] - t2 = [SWt, NEt, SEt] - else - t1 = [SWt, NWt, SEt] - t2 = [SEt, NWt, NEt] - end - elseif pattern == :UnionJack - t1 = [SWt, SEt] - t2 = [NWt, NEt] - if iseven(i + j) - push!(t1, NWt) - push!(t2, SEt) - else - push!(t1, NEt) - push!(t2, SWt) - end + diagonal_nw_se = (abs(mid1 - mid3) < abs(mid2 - mid3)) elseif pattern == :K1 - t1 = [SEt, SWt, NWt] - t2 = [NWt, NEt, SEt] + diagonal_nw_se = false + elseif pattern == :UnionJack + diagonal_nw_se = isodd(i + j) elseif pattern == :Random - if rand(mt, Bool) - t1 = [NWt, NEt, SEt] - t2 = [SEt, SWt, NWt] - else - t1 = [SWt, NWt, NEt] - t2 = [NEt, SEt, SWt] - end + diagonal_nw_se = rand(mt, Bool) + end + + if diagonal_nw_se + corners1 = [(xL, yL), (xL, yU), (xU, yL)] # SW, NW, SE + corners2 = [(xU, yL), (xL, yU), (xU, yU)] # SE, NW, NE else - error("pattern $pattern not currently supported") + corners1 = [(xL, yL), (xU, yU), (xU, yL)] # SW, NE, SE + corners2 = [(xL, yL), (xL, yU), (xU, yU)] # SW, NW, NE end - push!(T, t1) - push!(T, t2) + + push!( + segments, + SegmentPointRep{2,1}(corners1, [(fz(c...),) for c in corners1]), + ) + push!( + segments, + SegmentPointRep{2,1}(corners2, [(fz(c...),) for c in corners2]), + ) end - return PWLFunction{2}(X, Z, T, Dict(:structure => pattern)) + + return PWLFunction(segments, structure) end diff --git a/test/1D-pwl-instances.h5 b/test/1D-pwl-instances.h5 deleted file mode 100644 index 248ab63..0000000 Binary files a/test/1D-pwl-instances.h5 and /dev/null differ diff --git a/test/2D-pwl-instances.h5 b/test/2D-pwl-instances.h5 deleted file mode 100644 index e96e5a1..0000000 Binary files a/test/2D-pwl-instances.h5 and /dev/null differ diff --git a/test/runtests.jl b/test/runtests.jl index 3c89df5..67c98a4 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -4,135 +4,224 @@ # in the LICENSE.md file or at https://opensource.org/licenses/MIT. using PiecewiseLinearOpt +using HiGHS +using JuMP +using LinearAlgebra using Test -import Cbc -import JuMP -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) - -@testset "Univariate tests" begin - @testset "1D: $method" for method in methods_1D - model = JuMP.Model(optimizer) - 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 - @test JuMP.objective_value(model) ≈ 0.90997 rtol = 1e-4 - @test JuMP.objective_value(model) ≈ JuMP.value(z) rtol = 1e-4 - end +const PLO = PiecewiseLinearOpt + +optimizer = optimizer_with_attributes(HiGHS.Optimizer, MOI.Silent() => true) + +const methods_1D = [ + ConvexCombination(), + DisaggregatedLogarithmic(), + Incremental(), + LogarithmicEmbedding(), + LogarithmicIndependentBranching(), + NativeSOS2(), + ZigZagBinary(), + ZigZagInteger(), +] + +@testset "Simple univariate" for method in methods_1D + model = Model(optimizer) + @variable(model, x) + + s1 = PLO.SegmentPointRep{1,1}([(1.0,), (2.0,)], [(2.5,), (3.5,)]) + s2 = PLO.SegmentPointRep{1,1}([(2.0,), (3.0,)], [(3.5,), (1.0,)]) + pwl = PLO.PWLFunction([s1, s2], PLO.Intervals()) + + y = piecewiselinear(model, (x,), pwl; method = method) + @objective(model, Min, y[1]) + + optimize!(model) + @test termination_status(model) == MOI.OPTIMAL + @test value(x) ≈ 3.0 rtol = 1e-4 + @test value(y[1]) ≈ 1.0 rtol = 1e-4 +end + +@testset "Univariate pwlinear" begin + d = 0:0.01:1 + f = (xi -> xi^2) + fd = [f(xi) for xi in d] + pwl = UnivariatePWLFunction(d, f) + + model = Model(optimizer) + @variable(model, x) + y1 = piecewiselinear(model, (x,), pwl) + @constraint(model, x ≤ 0.75) + @objective(model, Max, y1[1]) + optimize!(model) + @test termination_status(model) == MOI.OPTIMAL + @test value(y1[1]) ≈ f(value(x)) rtol = 1e-4 + @test objective_value(model) ≈ 0.5625 rtol = 1e-4 + + model = Model(optimizer) + @variable(model, x) + y2 = piecewiselinear(model, x, d, f) + @constraint(model, x ≤ 0.75) + @objective(model, Max, y2) + optimize!(model) + @test termination_status(model) == MOI.OPTIMAL + @test value(y2) ≈ f(value(x)) rtol = 1e-4 + @test objective_value(model) ≈ 0.5625 rtol = 1e-4 + + model = Model(optimizer) + @variable(model, x) + y3 = piecewiselinear(model, x, d, fd) + @constraint(model, x ≤ 0.75) + @objective(model, Max, y3) + optimize!(model) + @test termination_status(model) == MOI.OPTIMAL + @test value(y3) ≈ f(value(x)) rtol = 1e-4 + @test objective_value(model) ≈ 0.5625 rtol = 1e-4 +end + +@testset "Bivariate pwlinear" begin + d = 0:0.05:1 + f = (xi, yi) -> xi^2 + yi^2 + pwl = BivariatePWLFunction(d, d, f) + + model = Model(optimizer) + @variable(model, x) + @variable(model, y) + z1 = piecewiselinear(model, (x, y), pwl) + @constraint(model, x ≤ 0.75) + @constraint(model, y ≤ 0.75) + @objective(model, Max, z1[1]) + optimize!(model) + @test termination_status(model) == MOI.OPTIMAL + @test value(x) ≈ 0.75 rtol = 1e-4 + @test value(y) ≈ 0.75 rtol = 1e-4 + @test value(z1[1]) ≈ 1.125 rtol = 1e-4 + + model = Model(optimizer) + @variable(model, x) + @variable(model, y) + z2 = piecewiselinear(model, x, y, d, d, f) + @constraint(model, x ≤ 0.75) + @constraint(model, y ≤ 0.75) + @objective(model, Max, z2) + optimize!(model) + @test termination_status(model) == MOI.OPTIMAL + @test value(x) ≈ 0.75 rtol = 1e-4 + @test value(y) ≈ 0.75 rtol = 1e-4 + @test value(z2) ≈ 1.125 rtol = 1e-4 +end + +const sos2_methods = [ + ConvexCombination(), + LogarithmicEmbedding(), + LogarithmicIndependentBranching(), + NativeSOS2(), + ZigZagBinary(), + ZigZagInteger(), +] +const methods_2D_gen = [ + ConvexCombination(), + DisaggregatedLogarithmic(), + #OptimalIndependentBranching(optimizer), + [NineStencil(sos2_method) for sos2_method in methods_1D]..., + [ + OptimalTriangleSelection(optimizer, sos2_method) for + sos2_method in methods_1D + ]..., + [SixStencil(sos2_method) for sos2_method in methods_1D]..., +] +@testset "Simple bivariate" for method in methods_2D_gen + model = Model(optimizer) + @variable(model, x[1:2]) + s1 = PLO.SegmentPointRep{2,1}( + [(0.0, 0.0), (0.0, 1.0), (1.0, 1.0)], + [(0.0,), (1.0,), (2.0,)], + ) + s2 = PLO.SegmentPointRep{2,1}( + [(0.0, 0.0), (1.0, 0.0), (1.0, 1.0)], + [(0.0,), (3.0,), (2.0,)], + ) + pwl = PLO.PWLFunction{2,1,PLO.SegmentPointRep{2,1}}( + [s1, s2], + PLO.UnstructuredTriangulation(), + ) + + y = piecewiselinear(model, (x[1], x[2]), pwl; method = method) + @objective(model, Min, y[1]) + + optimize!(model) + + @test termination_status(model) == MOI.OPTIMAL + @test value(x[1]) ≈ 0.0 rtol = 1e-4 + @test value(x[2]) ≈ 0.0 rtol = 1e-4 + @test value(y[1]) ≈ 0.0 rtol = 1e-4 end -@testset "Bivariate tests " begin - @testset "2D: $method, $pattern" for method in methods_2D, - pattern in patterns_2D - - model = JuMP.Model(optimizer) - JuMP.@variable(model, x[1:2]) - d = range(0; stop = 1, length = 8) - 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 - @test JuMP.value(z) ≈ 0.148753 rtol = 1e-4 - @test JuMP.objective_value(model) ≈ 0.148753 rtol = 1e-3 - @test JuMP.objective_value(model) ≈ JuMP.value(z) rtol = 1e-3 - end +@testset "1D: $method" for method in methods_1D + model = Model(optimizer) + @variable(model, x) + d = 7 + xs = collect(range(1; stop = 2π, length = (d + 1))) + zs = sin.(xs) + y = piecewiselinear(model, x, xs, zs; method = method) + @objective(model, Max, y) + optimize!(model) + @test termination_status(model) == MOI.OPTIMAL + @test value(x) ≈ 1.75474 rtol = 1e-4 + @test value(y) ≈ 0.98313 rtol = 1e-4 + @test objective_value(model) ≈ 0.98313 rtol = 1e-4 + @test objective_value(model) ≈ value(y) rtol = 1e-4 + @constraint(model, x ≤ 1.5y) + optimize!(model) + @test termination_status(model) == MOI.OPTIMAL + @test value(x) ≈ 1.36495 rtol = 1e-4 + @test value(y) ≈ 0.90997 rtol = 1e-4 + @test objective_value(model) ≈ 0.90997 rtol = 1e-4 + @test objective_value(model) ≈ value(y) rtol = 1e-4 end -# println("\nbivariate optimal IB scheme tests") -# @testset "2D: optimal IB, UnionJack" begin -# model = JuMP.Model(JuMP.with_optimizer(Cbc.Optimizer)) -# JuMP.@variable(model, x) -# JuMP.@variable(model, y) -# d = range(0,stop=1,length=3) -# 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 -# @test JuMP.value(z) ≈ 0.222300 rtol=1e-3 -# @test JuMP.objective_value(model) ≈ 0.222300 rtol=1e-3 -# @test JuMP.objective_value(model) ≈ JuMP.value(z) rtol=1e-3 -# end +patterns = [:Upper, :Lower, :BestFit, :K1, :UnionJack, :Random] +method_pattern = vec(collect(Iterators.product(methods_2D_gen, patterns))) +k1_methods = [ + (method, :K1) for method in [K1(sos2_method) for sos2_method in methods_1D] +] +uj_methods = [ + (method, :UnionJack) for + method in [UnionJack(sos2_method) for sos2_method in methods_1D] +] +append!(method_pattern, k1_methods) +append!(method_pattern, uj_methods) + +@testset "2D: $method, $pattern" for (method, pattern) in method_pattern + model = Model(optimizer) + @variable(model, x[1:2]) + d = range(0; stop = 1, length = 8) + f = (x1, x2) -> 2 * (x1 - 1 / 3)^2 + 3 * (x2 - 4 / 7)^4 + z = piecewiselinear( + model, + x[1], + x[2], + d, + d, + f; + method = method, + pattern = pattern, + ) + @objective(model, Min, z) + optimize!(model) + @test termination_status(model) == MOI.OPTIMAL + @test value(x[1]) ≈ 0.285714 rtol = 1e-4 + @test value(x[2]) ≈ 0.571429 rtol = 1e-4 + @test value(z) ≈ 0.004535 rtol = 1e-3 + @test objective_value(model) ≈ 0.004535 rtol = 1e-3 + @test objective_value(model) ≈ value(z) rtol = 1e-3 + + @constraint(model, x[1] ≥ 0.6) + optimize!(model) + + @test termination_status(model) == MOI.OPTIMAL + @test value(x[1]) ≈ 0.6 rtol = 1e-4 + @test value(x[2]) ≈ 0.571428 rtol = 1e-4 + @test value(z) ≈ 0.148753 rtol = 1e-4 + @test objective_value(model) ≈ 0.148753 rtol = 1e-3 + @test objective_value(model) ≈ value(z) rtol = 1e-3 +end