From bfa8db554285d0d550e8a88f48b53ebde2385ccf Mon Sep 17 00:00:00 2001 From: Oscar Dowson Date: Wed, 11 Sep 2024 14:26:39 +1200 Subject: [PATCH] Add JuliaFormatter (#53) --- .JuliaFormatter.toml | 9 + .github/workflows/format-check.yml | 30 + src/jump.jl | 1002 ++++++++++++++++++---------- src/types.jl | 92 ++- test/exhaustive-tests.jl | 150 +++-- test/runtests.jl | 58 +- 6 files changed, 900 insertions(+), 441 deletions(-) create mode 100644 .JuliaFormatter.toml create mode 100644 .github/workflows/format-check.yml diff --git a/.JuliaFormatter.toml b/.JuliaFormatter.toml new file mode 100644 index 0000000..f196f10 --- /dev/null +++ b/.JuliaFormatter.toml @@ -0,0 +1,9 @@ +# Configuration file for JuliaFormatter.jl +# For more information, see: https://domluna.github.io/JuliaFormatter.jl/stable/config/ + +always_for_in = true +always_use_return = true +margin = 80 +remove_extra_newlines = true +separate_kwargs_with_semicolon = true +short_to_long_function_def = true diff --git a/.github/workflows/format-check.yml b/.github/workflows/format-check.yml new file mode 100644 index 0000000..4f8e797 --- /dev/null +++ b/.github/workflows/format-check.yml @@ -0,0 +1,30 @@ +name: format-check +on: + push: + branches: + - master + pull_request: + types: [opened, synchronize, reopened] +jobs: + build: + runs-on: ubuntu-latest + steps: + - uses: julia-actions/setup-julia@latest + with: + version: '1' + - uses: actions/checkout@v4 + - name: Format check + shell: julia --color=yes {0} + run: | + using Pkg + # If you update the version, also update the style guide docs. + Pkg.add(PackageSpec(name="JuliaFormatter", version="1")) + using JuliaFormatter + format("."; verbose = true) + out = String(read(Cmd(`git diff`))) + if isempty(out) + exit(0) + end + @error "Some files have not been formatted !!!" + write(stdout, out) + exit(1) diff --git a/src/jump.jl b/src/jump.jl index 431c978..70f1393 100644 --- a/src/jump.jl +++ b/src/jump.jl @@ -15,21 +15,39 @@ function initPWL!(m::JuMP.Model) if !haskey(m.ext, :PWL) m.ext[:PWL] = PWLData() end - return nothing + return end const VarOrAff = Union{JuMP.VariableRef,JuMP.AffExpr} -function piecewiselinear(m::JuMP.Model, x::VarOrAff, d, f::Function; method=defaultmethod()) +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) + return piecewiselinear(m, x, d, fd; method = method) end -piecewiselinear(m::JuMP.Model, x::VarOrAff, d, fd; method=defaultmethod()) = - piecewiselinear(m, x, UnivariatePWLFunction(d, fd); method=method) +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()) +function piecewiselinear( + m::JuMP.Model, + x::VarOrAff, + pwl::UnivariatePWLFunction; + method = defaultmethod(), +) initPWL!(m) counter = m.ext[:PWL].counter counter += 1 @@ -39,71 +57,115 @@ function piecewiselinear(m::JuMP.Model, x::VarOrAff, pwl::UnivariatePWLFunction; n = length(d) if n != length(fd) - error("You provided a different number of breakpoints ($n) and function values at the breakpoints ($(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") + 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") - + 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, 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)) + δ = 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") + 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] + Δ = [(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) + 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) + 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)) + 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]) + 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.@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) + 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) + elseif method in (:LogarithmicIB, :LogIB) sos2_logarithmic_IB_formulation!(m, λ) elseif method == :CC sos2_cc_formulation!(m, λ) - elseif method in (:ZigZag,:ZZB) + elseif method in (:ZigZag, :ZZB) sos2_zigzag_formulation!(m, λ) - elseif method in (:ZigZagInteger,:ZZI) + elseif method in (:ZigZagInteger, :ZZI) sos2_zigzag_general_integer_formulation!(m, λ) elseif method == :GeneralizedCelaya sos2_generalized_celaya_formulation!(m, λ) @@ -115,21 +177,20 @@ function piecewiselinear(m::JuMP.Model, x::VarOrAff, pwl::UnivariatePWLFunction; 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") + 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 nothing + return end # function sos2_mc_formulation!(m::JuMP.Model, λ) # not currently used @@ -142,100 +203,162 @@ end # for i in 1:n-1 # JuMP.@constraint(m, γ[i,i] + γ[i,i+1] ≥ y[i]) # end -# return nothing +# 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 nothing + 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") + 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) + 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 nothing + 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") + 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 nothing + return end function sos2_zigzag_general_integer_formulation!(m::JuMP.Model, λ) counter = m.ext[:PWL].counter - n = length(λ)-1 - k = ceil(Int,log2(n)) + 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 nothing + 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)) + 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 nothing + 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)) + 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") + 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 nothing + return end function sos2_encoding_constraints!(m, λ, y, h, B) - n = length(λ)-1 + 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) + 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 + return end function reflected_gray_codes(k::Int) if k == 0 return Vector{Int}[] elseif k == 1 - return [[0],[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′)]) + 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 @@ -243,11 +366,13 @@ function zigzag_codes(k::Int) if k <= 0 error("Invalid code length $k") elseif k == 1 - return [[0],[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′]) + codes′ = zigzag_codes(k - 1) + return vcat( + [vcat(code, 0) for code in codes′], + [vcat(code, 1) for code in codes′], + ) end end @@ -255,12 +380,14 @@ function integer_zigzag_codes(k::Int) if k <= 0 error("Invalid code length $k") elseif k == 1 - return [[0],[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′]) + 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 @@ -268,18 +395,20 @@ function generalized_celaya_codes(k::Int) if k <= 0 error("Invalid code length $k") elseif k == 1 - return [[0],[1]] + return [[0], [1]] elseif k == 2 - return [[0,0],[1,1],[1,0],[0,1]] + return [[0, 0], [1, 1], [1, 0], [0, 1]] else - codes′ = generalized_celaya_codes(k-1) + codes′ = generalized_celaya_codes(k - 1) n = length(codes′) - hp = Int(n/2) - firstcodes = [codes′[i] for i in 1:hp] + 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]) + 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 @@ -287,16 +416,18 @@ function symmetric_celaya_codes(k::Int) if k <= 0 error("Invalid code length $k") elseif k == 1 - return [[0],[1]] + return [[0], [1]] else - codes′ = generalized_celaya_codes(k-1) + codes′ = generalized_celaya_codes(k - 1) n = length(codes′) - hp = Int(n/2) - firstcodes = [codes′[i] for i in 1:hp] + 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]) + 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 @@ -306,7 +437,7 @@ function zigzag_hyperplanes(k::Int) hp = zeros(Int, k) hp[i] = 1 for j in (i+1):k - hp[j] = 2^(j-i-1) + hp[j] = 2^(j - i - 1) end push!(hps, hp) end @@ -316,70 +447,72 @@ end function unit_vector_hyperplanes(k::Int) hps = Vector{Int}[] for i in 1:k - hp = zeros(Int,k) + hp = zeros(Int, k) hp[i] = 1 push!(hps, hp) end return hps end -generalized_celaya_hyperplanes(k::Int) = compute_hyperplanes(generalized_celaya_codes(k)) +function generalized_celaya_hyperplanes(k::Int) + return compute_hyperplanes(generalized_celaya_codes(k)) +end -symmetric_celaya_hyperplanes(k::Int) = compute_hyperplanes(symmetric_celaya_codes(k)) +function symmetric_celaya_hyperplanes(k::Int) + return compute_hyperplanes(symmetric_celaya_codes(k)) +end -function compute_hyperplanes(C::Vector{Vector{T}}) where T <: Number +function compute_hyperplanes(C::Vector{Vector{T}}) where {T<:Number} n = length(C) k = length(C[1]) - - d = zeros(n-1,k) + d = zeros(n - 1, k) for i in 1:n-1 - d[i,:] = C[i+1] - C[i] + 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 = 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]) + approx12 = isapprox(spanners[1], spanners[2]) + approx13 = isapprox(spanners[1], spanners[3]) + approx23 = isapprox(spanners[2], spanners[3]) if !approx12 - push!(indices,2) + push!(indices, 2) end if !approx13 && !(!approx12 && approx23) - push!(indices,3) + push!(indices, 3) end return spanners[indices] end - - indices = [1,2] + indices = [1, 2] spanners = Vector{Float64}[] while !isempty(indices) - if indices == [n-1] + 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 + 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) + if indices[end] != n - 1 + push!(indices, indices[end] + 1) else pop!(indices) indices[end] += 1 end else - if indices[end] != n-1 + if indices[end] != n - 1 indices[end] += 1 else pop!(indices) @@ -387,7 +520,6 @@ function compute_hyperplanes(C::Vector{Vector{T}}) where T <: Number end end end - keepers = [1] for i in 2:length(spanners) alreadyin = false @@ -401,7 +533,6 @@ function compute_hyperplanes(C::Vector{Vector{T}}) where T <: Number push!(keepers, i) end end - return spanners[keepers] end @@ -412,7 +543,7 @@ function canonical!(v::Vector{Float64}) v[j] = 0 end end - sgn = sign(v[findfirst([x!=0 for x in v])]) + 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 @@ -425,65 +556,61 @@ 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 + E[i, j] = true end end - if haskey(m.ext[:OptimalIBCache], pwl.T) xx, yy = m.ext[:OptimalIBCache][pwl.T] - t = size(xx,1) + 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) - + 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) - + 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) + 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) + 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) + 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) + 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) @@ -494,11 +621,9 @@ function optimal_IB_scheme!(m::JuMP.Model, λ, pwl, subsolver) # IB: independent break end end - m.ext[:OptimalIBCache][pwl.T] = (xx,yy) + 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ʸ) @@ -507,26 +632,56 @@ function optimal_IB_scheme!(m::JuMP.Model, λ, pwl, subsolver) # IB: independent 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) + 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 nothing + return end -piecewiselinear(m::JuMP.Model, x::VarOrAff, y::VarOrAff, dˣ, dʸ, f::Function; method=defaultmethod()) = - piecewiselinear(m, x, y, BivariatePWLFunction(dˣ, dʸ, f); method=method) +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) +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())") + 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 @@ -536,132 +691,206 @@ function piecewiselinear(m::JuMP.Model, x₁::VarOrAff, x₂::VarOrAff, pwl::Biv 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") + 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ʸ])) + 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ˣ])) + 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) + 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 + 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") - + 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) + 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) + 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]] + 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]) + 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] + 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) + @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] + 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]) + @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) + 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") + γ = 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) + 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") + 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]) + 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.@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) - + 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") + 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ʸ) + 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) + 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)])) + 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.") + 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.") + 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ˣ] + 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) + if method in (:Logarithmic, :Log) sos2_logarithmic_formulation!(m, Tx) sos2_logarithmic_formulation!(m, Ty) - elseif method in (:LogarithmicIB,:LogIB) + elseif method in (:LogarithmicIB, :LogIB) sos2_logarithmic_IB_formulation!(m, Tx) sos2_logarithmic_IB_formulation!(m, Ty) - elseif method in (:ZigZag,:ZZB) + elseif method in (:ZigZag, :ZZB) sos2_zigzag_formulation!(m, Tx) sos2_zigzag_formulation!(m, Ty) - elseif method in (:ZigZagInteger,:ZZI) + elseif method in (:ZigZagInteger, :ZZI) sos2_zigzag_general_integer_formulation!(m, Tx) sos2_zigzag_general_integer_formulation!(m, Ty) elseif method == :GeneralizedCelaya @@ -671,10 +900,30 @@ function piecewiselinear(m::JuMP.Model, x₁::VarOrAff, x₂::VarOrAff, pwl::Biv 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.@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 @@ -686,7 +935,9 @@ function piecewiselinear(m::JuMP.Model, x₁::VarOrAff, x₂::VarOrAff, pwl::Biv numT = 0 minˣ, minʸ = uˣ[1], uʸ[1] # find the index of the bottom-left point - idx = findfirst(pwl.x) do w; w == (minˣ, minʸ); end + idx = findfirst(pwl.x) do w + return w == (minˣ, minʸ) + end for t in pwl.T if idx in t numT += 1 @@ -694,44 +945,65 @@ function piecewiselinear(m::JuMP.Model, x₁::VarOrAff, x₂::VarOrAff, pwl::Biv end @assert 1 <= numT <= 2 - w = JuMP.@variable(m, binary=true, base_name="w_$counter") + 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) + 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) + 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.") + 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ʸ] + 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] + 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))) + 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 @@ -739,18 +1011,20 @@ function piecewiselinear(m::JuMP.Model, x₁::VarOrAff, x₂::VarOrAff, pwl::Biv if haskey(m.ext[:OptimalTriSelectCache], E) xx, yy = m.ext[:OptimalTriSelectCache][E] - t = JuMP.size(xx,1) + 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())") + 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) + 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) + 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 @@ -758,27 +1032,46 @@ function piecewiselinear(m::JuMP.Model, x₁::VarOrAff, x₂::VarOrAff, pwl::Biv 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) + 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) + 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) + (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 @@ -789,64 +1082,82 @@ function piecewiselinear(m::JuMP.Model, x₁::VarOrAff, x₂::VarOrAff, pwl::Biv else xx = JuMP.getvalue(xˢᵘᵇ) yy = JuMP.getvalue(yˢᵘᵇ) - m.ext[:OptimalTriSelectCache][E] = (xx,yy) + m.ext[:OptimalTriSelectCache][E] = (xx, yy) break end end end - y = JuMP.@variable(m, [1:t], Bin, base_name="Δselect_$counter") + 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) + 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") + 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 + 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 + 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) + 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 + 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ʸ)] + 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 + if ((im, jM) in IJ) && ((iM, jm) in IJ) + Eⁿᵉ[im, jm] = true else - @assert (im,jm) in IJ && (iM,jM) in IJ + @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") + 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}}() @@ -854,19 +1165,19 @@ function piecewiselinear(m::JuMP.Model, x₁::VarOrAff, x₂::VarOrAff, pwl::Biv 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) + if !(1 ≤ i ≤ nˣ - 1) continue end - if !(1 ≤ j ≤ nʸ-1) + if !(1 ≤ j ≤ nʸ - 1) continue # should never happen end - if Eⁿᵉ[i,j] # if we need to cover the edge... + 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)) + push!(Aᵒ, (i, j)) + push!(Bᵒ, (i + 1, j + 1)) else - push!(Aᵒ, (i+1,j+1)) - push!(Bᵒ, (i ,j )) + push!(Aᵒ, (i + 1, j + 1)) + push!(Bᵒ, (i, j)) end SWinA = !SWinA end @@ -876,47 +1187,50 @@ function piecewiselinear(m::JuMP.Model, x₁::VarOrAff, x₂::VarOrAff, pwl::Biv SWinA = true for j in (offʸ+1):(nʸ-1) i = j - offʸ - if !(1 ≤ i ≤ nˣ-1) + if !(1 ≤ i ≤ nˣ - 1) continue end - if Eⁿᵉ[i,j] + if Eⁿᵉ[i, j] if SWinA - push!(Aᵒ, (i ,j )) - push!(Bᵒ, (i+1,j+1)) + push!(Aᵒ, (i, j)) + push!(Bᵒ, (i + 1, j + 1)) else - push!(Aᵒ, (i+1,j+1)) - push!(Bᵒ, (i ,j )) + 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) + 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") + 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 + # j = offˣ - i + 2 for j in 1:(nʸ-1) i = nˣ - j - offˣ - if !(1 ≤ i ≤ nˣ-1) + if !(1 ≤ i ≤ nˣ - 1) continue end - if !Eⁿᵉ[i,j] + if !Eⁿᵉ[i, j] if SEinA - push!(Aᵒ, (i+1,j )) - push!(Bᵒ, (i ,j+1)) + push!(Aᵒ, (i + 1, j)) + push!(Bᵒ, (i, j + 1)) else - push!(Aᵒ, (i ,j+1)) - push!(Bᵒ, (i+1,j )) + push!(Aᵒ, (i, j + 1)) + push!(Bᵒ, (i + 1, j)) end SEinA = !SEinA end @@ -926,29 +1240,31 @@ function piecewiselinear(m::JuMP.Model, x₁::VarOrAff, x₂::VarOrAff, pwl::Biv SEinA = true for j in (offʸ+1):(nʸ-1) i = nˣ - j + offʸ - if !(1 ≤ i ≤ nˣ-1) + if !(1 ≤ i ≤ nˣ - 1) continue end - if !Eⁿᵉ[i,j] + if !Eⁿᵉ[i, j] if SEinA - push!(Aᵒ, (i+1,j )) - push!(Bᵒ, (i ,j+1)) + push!(Aᵒ, (i + 1, j)) + push!(Bᵒ, (i, j + 1)) else - push!(Aᵒ, (i ,j+1)) - push!(Bᵒ, (i+1,j )) + 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) + 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/types.jl b/src/types.jl index f516220..5ae0257 100644 --- a/src/types.jl +++ b/src/types.jl @@ -10,7 +10,12 @@ struct PWLFunction{D} meta::Dict end -function PWLFunction{D}(x::Vector{NTuple{D}}, z::Vector, T::Vector{Vector}, meta::Dict) where {D} +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) @@ -24,68 +29,83 @@ const UnivariatePWLFunction = PWLFunction{1} 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]) + return PWLFunction( + Tuple{Float64}[(xx,) for xx in x], + convert(Vector{Float64}, z), + [[i, i + 1] for i in 1:length(x)-1], + ) 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]) + 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], + ) end const BivariatePWLFunction = PWLFunction{2} -function BivariatePWLFunction(x, y, fz::Function; pattern=:BestFit, seed=hash((length(x),length(y)))) +function BivariatePWLFunction( + x, + y, + fz::Function; + pattern = :BestFit, + seed = hash((length(x), length(y))), +) @assert issorted(x) @assert issorted(y) - X = vec(collect(Base.product(x,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) + Z = map(t -> convert(Float64, fz(t...)), X) T = Vector{Vector{Int}}() m = length(x) n = length(y) - 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] + 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) - + mid1 = 0.5 * (SW + NE) + mid2 = 0.5 * (NW + SE) if pattern == :Upper if mid1 > mid2 - t1 = [SWt,NWt,NEt] - t2 = [SWt,NEt,SEt] + t1 = [SWt, NWt, NEt] + t2 = [SWt, NEt, SEt] else - t1 = [SWt,NWt,SEt] - t2 = [SEt,NWt,NEt] + t1 = [SWt, NWt, SEt] + t2 = [SEt, NWt, NEt] end elseif pattern == :Lower if mid1 > mid2 - t1 = [SWt,NWt,SEt] - t2 = [SEt,NWt,NEt] + t1 = [SWt, NWt, SEt] + t2 = [SEt, NWt, NEt] else - t1 = [SWt,NWt,NEt] - t2 = [SWt,NEt,SEt] + t1 = [SWt, NWt, NEt] + t2 = [SWt, NEt, SEt] end 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] + 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] + t1 = [SWt, NWt, SEt] + t2 = [SEt, NWt, NEt] end elseif pattern == :UnionJack - t1 = [SWt,SEt] - t2 = [NWt,NEt] - if iseven(i+j) + t1 = [SWt, SEt] + t2 = [NWt, NEt] + if iseven(i + j) push!(t1, NWt) push!(t2, SEt) else @@ -93,23 +113,21 @@ function BivariatePWLFunction(x, y, fz::Function; pattern=:BestFit, seed=hash((l push!(t2, SWt) end elseif pattern == :K1 - t1 = [SEt,SWt,NWt] - t2 = [NWt,NEt,SEt] + t1 = [SEt, SWt, NWt] + t2 = [NWt, NEt, SEt] elseif pattern == :Random if rand(mt, Bool) - t1 = [NWt,NEt,SEt] - t2 = [SEt,SWt,NWt] + t1 = [NWt, NEt, SEt] + t2 = [SEt, SWt, NWt] else - t1 = [SWt,NWt,NEt] - t2 = [NEt,SEt,SWt] + t1 = [SWt, NWt, NEt] + t2 = [NEt, SEt, SWt] end else error("pattern $pattern not currently supported") end - push!(T, t1) push!(T, t2) end - - return PWLFunction{2}(X, Z, T, Dict(:structure=>pattern)) + return PWLFunction{2}(X, Z, T, Dict(:structure => pattern)) end diff --git a/test/exhaustive-tests.jl b/test/exhaustive-tests.jl index 4ccac9b..cda6c79 100644 --- a/test/exhaustive-tests.jl +++ b/test/exhaustive-tests.jl @@ -7,90 +7,162 @@ # solver = CbcSolver(logLevel=0, integerTolerance=1e-9, primalTolerance=1e-9, ratioGap=1e-8) using Gurobi -solver = GurobiSolver(OutputFlag=0) +solver = GurobiSolver(; OutputFlag = 0) # using CPLEX # solver = CplexSolver(CPX_PARAM_SCRIND=0) - using JuMP using PiecewiseLinearOpt using Base.Test using HDF5 -methods_1D = (:CC,:MC,:Logarithmic,:LogarithmicIB,:ZigZag,:ZigZagInteger,:SOS2,:GeneralizedCelaya,:SymmetricCelaya,:Incremental,:DisaggLogarithmic) -methods_2D = (:CC,:MC,:Logarithmic,:LogarithmicIB,:ZigZag,:ZigZagInteger,:SOS2,:GeneralizedCelaya,:SymmetricCelaya,:DisaggLogarithmic) +methods_1D = ( + :CC, + :MC, + :Logarithmic, + :LogarithmicIB, + :ZigZag, + :ZigZagInteger, + :SOS2, + :GeneralizedCelaya, + :SymmetricCelaya, + :Incremental, + :DisaggLogarithmic, +) +methods_2D = ( + :CC, + :MC, + :Logarithmic, + :LogarithmicIB, + :ZigZag, + :ZigZagInteger, + :SOS2, + :GeneralizedCelaya, + :SymmetricCelaya, + :DisaggLogarithmic, +) # methods_2D = (:OptimalIB,) # very slow on all solvers -patterns_2D = (:Upper,:Lower,:UnionJack,:K1,:Random) # not :BestFit because requires function values at midpoints, :OptimalTriangleSelection and :Stencil not supported currently +patterns_2D = (:Upper, :Lower, :UnionJack, :K1, :Random) # not :BestFit because requires function values at midpoints, :OptimalTriangleSelection and :Stencil not supported currently # tests on network flow model with piecewise-linear objective # instance data loaded from .h5 files -instance_data_1D = joinpath(dirname(@__FILE__),"1D-pwl-instances.h5") -instance_data_2D = joinpath(dirname(@__FILE__),"2D-pwl-instances.h5") +instance_data_1D = joinpath(dirname(@__FILE__), "1D-pwl-instances.h5") +instance_data_2D = joinpath(dirname(@__FILE__), "2D-pwl-instances.h5") println("\nunivariate tests") @testset "instance $instance" for instance in ["10104_1_concave_1"] - demand = h5read(instance_data_1D, string(instance,"/demand")) - supply = h5read(instance_data_1D, string(instance,"/supply")) + demand = h5read(instance_data_1D, string(instance, "/demand")) + supply = h5read(instance_data_1D, string(instance, "/supply")) ndem = size(demand, 1) nsup = size(supply, 1) - rawd = h5read(instance_data_1D, string(instance,"/d")) - d = rawd[1,:] - rawfd = h5read(instance_data_1D, string(instance,"/fd")) - fd = rawfd[1,:] + rawd = h5read(instance_data_1D, string(instance, "/d")) + d = rawd[1, :] + rawfd = h5read(instance_data_1D, string(instance, "/fd")) + fd = rawfd[1, :] objval1 = NaN @testset "1D: $method" for method in methods_1D - model = Model(solver=solver) - @variable(model, x[1:nsup,1:ndem] ≥ 0) - @constraint(model, [j in 1:ndem], sum(x[i,j] for i in 1:nsup) == demand[j]) - @constraint(model, [i in 1:nsup], sum(x[i,j] for j in 1:ndem) == supply[i]) - @objective(model, Min, sum(piecewiselinear(model, x[i,j], d, fd, method=method) for i in 1:nsup, j in 1:ndem)) + model = Model(; solver = solver) + @variable(model, x[1:nsup, 1:ndem] ≥ 0) + @constraint( + model, + [j in 1:ndem], + sum(x[i, j] for i in 1:nsup) == demand[j] + ) + @constraint( + model, + [i in 1:nsup], + sum(x[i, j] for j in 1:ndem) == supply[i] + ) + @objective( + model, + Min, + sum( + piecewiselinear(model, x[i, j], d, fd; method = method) for + i in 1:nsup, j in 1:ndem + ) + ) @test solve(model) == :Optimal if isnan(objval1) objval1 = getobjectivevalue(model) else - @test getobjectivevalue(model) ≈ objval1 rtol=1e-4 + @test getobjectivevalue(model) ≈ objval1 rtol = 1e-4 end end end println("\nbivariate tests") # @testset "numpieces $numpieces, variety $variety, objective $objective" for numpieces in [4,8], variety in 1:5, objective in 1:20 -@testset "numpieces $numpieces, variety $variety, objective $objective" for numpieces in [4], variety in 1:1, objective in 1:1 - instance = string(numpieces,"_",variety,"_",objective) +@testset "numpieces $numpieces, variety $variety, objective $objective" for numpieces in + [4], + variety in 1:1, + objective in 1:1 + + instance = string(numpieces, "_", variety, "_", objective) - demand = h5read(instance_data_2D, string(instance,"/demand")) - supply = h5read(instance_data_2D, string(instance,"/supply")) + demand = h5read(instance_data_2D, string(instance, "/demand")) + supply = h5read(instance_data_2D, string(instance, "/supply")) ndem = size(demand, 1) nsup = size(supply, 1) - rawd = h5read(instance_data_2D, string(instance,"/d")) - d = rawd[1,:] - rawfd = h5read(instance_data_2D, string(instance,"/fd")) - fˣ = reshape(rawfd[1,:], length(d), length(d)) + rawd = h5read(instance_data_2D, string(instance, "/d")) + d = rawd[1, :] + rawfd = h5read(instance_data_2D, string(instance, "/fd")) + fˣ = reshape(rawfd[1, :], length(d), length(d)) ˣtoⁱ = Dict(d[p] => p for p in 1:length(d)) - f = (pˣ,pʸ) -> fˣ[ˣtoⁱ[pˣ],ˣtoⁱ[pʸ]] + f = (pˣ, pʸ) -> fˣ[ˣtoⁱ[pˣ], ˣtoⁱ[pʸ]] objval1 = NaN - @testset "2D: $method, $pattern" for method in methods_2D, pattern in patterns_2D - model = Model(solver=solver) - @variable(model, x[1:nsup,1:ndem] ≥ 0) - @constraint(model, [j in 1:ndem], sum(x[i,j] for i in 1:nsup) == demand[j]) - @constraint(model, [i in 1:nsup], sum(x[i,j] for j in 1:ndem) == supply[i]) - @variable(model, y[1:nsup,1:ndem] ≥ 0) - @constraint(model, [j in 1:ndem], sum(y[i,j] for i in 1:nsup) == demand[j]) - @constraint(model, [i in 1:nsup], sum(y[i,j] for j in 1:ndem) == supply[i]) - - @objective(model, Min, sum(piecewiselinear(model, x[i,j], y[i,j], BivariatePWLFunction(d, d, f, pattern=pattern), method=method, subsolver=solver) for i in 1:nsup, j in 1:ndem)) + @testset "2D: $method, $pattern" for method in methods_2D, + pattern in patterns_2D + + model = Model(; solver = solver) + @variable(model, x[1:nsup, 1:ndem] ≥ 0) + @constraint( + model, + [j in 1:ndem], + sum(x[i, j] for i in 1:nsup) == demand[j] + ) + @constraint( + model, + [i in 1:nsup], + sum(x[i, j] for j in 1:ndem) == supply[i] + ) + @variable(model, y[1:nsup, 1:ndem] ≥ 0) + @constraint( + model, + [j in 1:ndem], + sum(y[i, j] for i in 1:nsup) == demand[j] + ) + @constraint( + model, + [i in 1:nsup], + sum(y[i, j] for j in 1:ndem) == supply[i] + ) + + @objective( + model, + Min, + sum( + piecewiselinear( + model, + x[i, j], + y[i, j], + BivariatePWLFunction(d, d, f; pattern = pattern); + method = method, + subsolver = solver, + ) for i in 1:nsup, j in 1:ndem + ) + ) @test solve(model) == :Optimal if isnan(objval1) objval1 = getobjectivevalue(model) else - @test getobjectivevalue(model) ≈ objval1 rtol=1e-4 + @test getobjectivevalue(model) ≈ objval1 rtol = 1e-4 end end end diff --git a/test/runtests.jl b/test/runtests.jl index cbfe74d..3c89df5 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -49,51 +49,65 @@ patterns_2D = ( # :Stencil ) -optimizer=JuMP.optimizer_with_attributes(Cbc.Optimizer, MOI.Silent() => true) +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) + 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 + @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 + @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 end @testset "Bivariate tests " begin - @testset "2D: $method, $pattern" for method in methods_2D, pattern in patterns_2D + @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) + 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 + @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 + @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 end