diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 3743a4b..4ddbeeb 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -1,39 +1,73 @@ name: CI on: push: - branches: [master] + branches: + - master + tags: ["*"] pull_request: - types: [opened, synchronize, reopened] -# needed to allow julia-actions/cache to delete old caches that it has created -permissions: - actions: write - contents: read + workflow_dispatch: +concurrency: + # Skip intermediate builds: always. + # Cancel intermediate builds: only if it is a pull request build. + group: ${{ github.workflow }}-${{ github.ref }} + cancel-in-progress: ${{ startsWith(github.ref, 'refs/pull/') }} jobs: test: - name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }} + name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} runs-on: ${{ matrix.os }} + timeout-minutes: 60 + permissions: # needed to allow julia-actions/cache to proactively delete old caches that it has created + actions: write + contents: read strategy: fail-fast: false matrix: - version: ['1.6', '1'] - os: [ubuntu-latest, macOS-latest, windows-latest] - arch: [x64] - include: - # Also test against 32-bit Linux. - - version: '1' - os: ubuntu-latest - arch: x86 + version: + - "1.10" + - "1.12" + - "pre" + os: + - ubuntu-latest + arch: + - x64 steps: - - uses: actions/checkout@v4 + - uses: actions/checkout@v6 - uses: julia-actions/setup-julia@v2 with: version: ${{ matrix.version }} arch: ${{ matrix.arch }} - - uses: julia-actions/cache@v1 + - uses: julia-actions/cache@v2 + - run: sudo apt-get update && sudo apt-get install -y xorg-dev mesa-utils xvfb libgl1 freeglut3-dev libxrandr-dev libxinerama-dev libxcursor-dev libxi-dev libxext-dev xsettingsd x11-xserver-utils - uses: julia-actions/julia-buildpkg@v1 - uses: julia-actions/julia-runtest@v1 - uses: julia-actions/julia-processcoverage@v1 - - uses: codecov/codecov-action@v4 + - uses: codecov/codecov-action@v5 with: - file: lcov.info + files: lcov.info token: ${{ secrets.CODECOV_TOKEN }} + fail_ci_if_error: false + docs: + name: Documentation + runs-on: ubuntu-latest + permissions: + actions: write # needed to allow julia-actions/cache to proactively delete old caches that it has created + contents: write + statuses: write + steps: + - uses: actions/checkout@v6 + - uses: julia-actions/setup-julia@v2 + with: + version: "1" + - uses: julia-actions/cache@v2 + - uses: julia-actions/julia-buildpkg@v1 + - uses: julia-actions/julia-docdeploy@v1 + env: + GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} # For authentication with GitHub Actions token + DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }} # For authentication with SSH deploy key + - name: Run doctests + shell: julia --project=docs --color=yes {0} + run: | + using Documenter: DocMeta, doctest + using PiecewiseLinearOpt + DocMeta.setdocmeta!(PiecewiseLinearOpt, :DocTestSetup, :(using PiecewiseLinearOpt); recursive=true) + doctest(PiecewiseLinearOpt) diff --git a/.gitignore b/.gitignore index b1fb063..0917646 100644 --- a/.gitignore +++ b/.gitignore @@ -3,3 +3,7 @@ *.jl.mem Manifest.toml *.DS_Store +.DS_Store +.vscode/settings.json +coverage/ +docs/build/ \ No newline at end of file diff --git a/Project.toml b/Project.toml index 50d9b8c..4d21aff 100644 --- a/Project.toml +++ b/Project.toml @@ -6,18 +6,21 @@ 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" +LinearAlgebra = "1" +Random = "1" +Aqua = "0.8" +Test = "1" [extras] -Cbc = "9961bab8-2fa3-5c5a-9d89-47fab24efd76" +Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" +HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["Cbc", "Test"] +test = ["Aqua", "Test", "HiGHS"] diff --git a/README.md b/README.md index abfab63..298fc46 100644 --- a/README.md +++ b/README.md @@ -2,6 +2,7 @@ [![Build Status](https://github.com/jump-dev/PiecewiseLinearOpt.jl/workflows/CI/badge.svg)](https://github.com/jump-dev/PiecewiseLinearOpt.jl/actions?query=workflow%3ACI) [![codecov](https://codecov.io/gh/jump-dev/PiecewiseLinearOpt.jl/branch/master/graph/badge.svg)](https://codecov.io/gh/jump-dev/PiecewiseLinearOpt.jl) +[![Aqua QA](https://juliatesting.github.io/Aqua.jl/dev/assets/badge.svg)](https://github.com/JuliaTesting/Aqua.jl) [PiecewiseLinearOpt.jl](https://github.com/jump-dev/PiecewiseLinearOpt.jl) is a JuMP extension for modeling optimization problems containing piecewise linear @@ -33,37 +34,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 +57,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 +78,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/docs/Project.toml b/docs/Project.toml new file mode 100644 index 0000000..4e6fbc0 --- /dev/null +++ b/docs/Project.toml @@ -0,0 +1,5 @@ +[deps] +Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" +HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b" +JuMP = "4076af6c-e467-56ae-b986-b466b2749572" +PiecewiseLinearOpt = "0f51c51e-adfa-5141-8a04-d40246b8977c" diff --git a/docs/make.jl b/docs/make.jl new file mode 100644 index 0000000..675bc9a --- /dev/null +++ b/docs/make.jl @@ -0,0 +1,34 @@ + +# julia --project=docs docs/make.jl +using PiecewiseLinearOpt +using Documenter + +DocMeta.setdocmeta!( + PiecewiseLinearOpt, + :DocTestSetup, + :(using PiecewiseLinearOpt, JuMP); + recursive = true, +) + +makedocs(; + modules = [PiecewiseLinearOpt], + sitename = "PiecewiseLinearOpt.jl", + authors = "Joey Huchette and contributors", + format = Documenter.HTML(; + canonical = "https://jump-dev.github.io/PiecewiseLinearOpt.jl", + edit_link = "master", + assets = String[], + ), + pages = [ + "Home" => "index.md", + "Tutorial" => "tutorial.md", + "Formulation Methods" => "methods.md", + "API Reference" => "api.md", + ], + warnonly = true, +) + +deploydocs(; + repo = "github.com/jump-dev/PiecewiseLinearOpt.jl", + devbranch = "master", +) diff --git a/docs/src/api.md b/docs/src/api.md new file mode 100644 index 0000000..b0bf6cb --- /dev/null +++ b/docs/src/api.md @@ -0,0 +1,58 @@ +# API Reference + +## Main function + +```@docs +piecewiselinear +``` + +## PWL function types + +```@docs +PWLFunction +UnivariatePWLFunction +BivariatePWLFunction +``` + +## Direction enum + +```@docs +PiecewiseLinearOpt.DIRECTION +``` + +## Univariate methods + +```@docs +Logarithmic +LogarithmicEmbedding +LogarithmicIndependentBranching +Incremental +NativeSOS2 +ZigZagBinary +ZigZagInteger +``` + +## Bivariate methods + +```@docs +K1 +SixStencil +NineStencil +UnionJack +``` + +## Multivariate methods + +```@docs +ConvexCombination +DisaggregatedLogarithmic +MultipleChoice +``` + +## Internal types + +```@docs +PiecewiseLinearOpt.SegmentPointRep +PiecewiseLinearOpt.SegmentHyperplaneRep +PiecewiseLinearOpt.AffineFunction +``` diff --git a/docs/src/index.md b/docs/src/index.md new file mode 100644 index 0000000..4b68d9b --- /dev/null +++ b/docs/src/index.md @@ -0,0 +1,67 @@ +```@meta +CurrentModule = PiecewiseLinearOpt +``` + +# PiecewiseLinearOpt.jl + +[PiecewiseLinearOpt.jl](https://github.com/jump-dev/PiecewiseLinearOpt.jl) is a +[JuMP](https://github.com/jump-dev/JuMP.jl) extension for modeling optimization +problems containing piecewise linear functions. It provides MIP (mixed-integer +programming) formulations for embedding piecewise linear functions into +optimization models. + +This package is an accompaniment to the paper: + +> J. Huchette and J.P. Vielma, [_Nonconvex piecewise linear functions: Advanced +> formulations and simple modeling tools_](https://arxiv.org/abs/1708.00050), +> Operations Research, 2019. + +## Features + +- **Univariate piecewise linear functions**: Approximate any univariate function + with a piecewise linear function defined on a set of breakpoints. +- **Bivariate piecewise linear functions**: Approximate bivariate functions on + triangulated grids with multiple triangulation patterns. +- **Multiple MIP formulations**: Choose from a variety of formulations with + different trade-offs in terms of the number of binary variables, continuous + variables, and constraints. +- **Graph, epigraph, and hypograph modes**: Model the graph of the function + exactly, or relax to epigraph/hypograph constraints. + +## Installation + +Install PiecewiseLinearOpt using the Julia package manager: + +```julia +import Pkg +Pkg.add("PiecewiseLinearOpt") +``` + +## Getting help + +If you need help, please ask a question on the +[JuMP community forum](https://jump.dev/forum). + +If you have a reproducible example of a bug, please open a +[GitHub issue](https://github.com/jump-dev/PiecewiseLinearOpt.jl/issues/new). + +## License + +`PiecewiseLinearOpt.jl` is licensed under the +[MIT license](https://github.com/jump-dev/PiecewiseLinearOpt.jl/blob/master/LICENSE.md). + +## Quick start + +```julia +using JuMP, PiecewiseLinearOpt, HiGHS + +model = Model(HiGHS.Optimizer) +set_silent(model) +@variable(model, x) +# Approximate sin(x) on [0, 2π] with breakpoints every 0.5 +z = piecewiselinear(model, x, 0:0.5:2pi, sin) +@objective(model, Max, z) +optimize!(model) +value(x) # ≈ π/2 +value(z) # ≈ 1.0 +``` diff --git a/docs/src/methods.md b/docs/src/methods.md new file mode 100644 index 0000000..12fa2b2 --- /dev/null +++ b/docs/src/methods.md @@ -0,0 +1,199 @@ +# Formulation Methods + +PiecewiseLinearOpt.jl provides a variety of MIP (mixed-integer programming) +formulations for embedding piecewise linear functions into optimization models. +Each formulation makes different trade-offs in terms of the number of binary +variables, continuous variables, and constraints. The choice of formulation can +significantly affect solver performance. + +For a detailed theoretical treatment, see: + +> J. Huchette and J.P. Vielma, [_Nonconvex piecewise linear functions: Advanced +> formulations and simple modeling tools_](https://arxiv.org/abs/1708.00050), +> Operations Research, 2019. + +## Univariate methods + +These methods are designed for one-dimensional piecewise linear functions. Given +a function with ``n`` breakpoints (i.e., ``n - 1`` linear segments), the methods +differ in how they enforce the SOS2 (Special Ordered Set of Type 2) condition on +the convex combination weights. + +| Method | Binary variables | Type | +|:-------|:-----------------|:-----| +| [`Logarithmic`](@ref) / [`LogarithmicEmbedding`](@ref) | ``\lceil \log_2(n-1) \rceil`` | Embedding | +| [`LogarithmicIndependentBranching`](@ref) | ``\lceil \log_2(n-1) \rceil`` | Independent branching | +| [`ZigZagBinary`](@ref) | ``\lceil \log_2(n-1) \rceil`` | Encoding | +| [`ZigZagInteger`](@ref) | ``\lceil \log_2(n-1) \rceil`` | Encoding (integer) | +| [`Incremental`](@ref) | ``n - 1`` | Incremental | +| [`NativeSOS2`](@ref) | 0 (uses solver SOS2) | Native | +| [`ConvexCombination`](@ref) | ``n - 1`` | Big-M | + +### Logarithmic / LogarithmicEmbedding + +```julia +Logarithmic() # or equivalently, LogarithmicEmbedding() +``` + +The default method for univariate functions. Uses a logarithmic number of binary +variables based on reflected Gray codes. This is often the best choice for +general-purpose use. + +### LogarithmicIndependentBranching + +```julia +LogarithmicIndependentBranching() +``` + +Similar to `Logarithmic`, but uses an independent branching scheme. Also uses +a logarithmic number of binary variables. + +### ZigZagBinary + +```julia +ZigZagBinary() +``` + +Uses a zig-zag encoding with binary variables. The number of binary variables is +logarithmic in the number of breakpoints. + +### ZigZagInteger + +```julia +ZigZagInteger() +``` + +A variant of the zig-zag formulation that uses general integer variables instead +of purely binary variables. The number of integer variables is logarithmic. + +### Incremental + +```julia +Incremental() +``` + +Uses one binary variable per segment. The formulation is based on an incremental +representation of the piecewise linear function. Simple but uses more binary +variables than the logarithmic methods. + +### NativeSOS2 + +```julia +NativeSOS2() +``` + +Delegates the SOS2 constraint directly to the solver using JuMP's built-in SOS2 +support. This method does not add any binary variables itself, but relies on the +solver to handle the SOS2 constraint natively. Not all solvers support SOS2 +constraints. + +## Bivariate methods + +These methods are designed for two-dimensional piecewise linear functions defined +on triangulated rectangular grids. They combine SOS2 constraints along each axis +with a triangle selection mechanism. + +Each bivariate method accepts an optional `axis_method` argument that specifies +which univariate formulation to use for the axis-aligned SOS2 constraints. The +default is `Logarithmic()`. + +| Method | Triangle selection variables | Grid requirement | +|:-------|:---------------------------|:-----------------| +| [`SixStencil`](@ref) | 6 binary | Any grid triangulation | +| [`NineStencil`](@ref) | 9 binary | Any grid triangulation | +| [`K1`](@ref) | 2 binary | K1 triangulation | +| [`UnionJack`](@ref) | 1 binary | Union Jack triangulation | + +### SixStencil + +```julia +SixStencil() # uses Logarithmic() for axis SOS2 +SixStencil(Incremental()) # uses Incremental() for axis SOS2 +``` + +The default method for bivariate functions. Works with any grid triangulation +and uses 6 binary variables for triangle selection. A good general-purpose +choice. + +### NineStencil + +```julia +NineStencil() +NineStencil(Incremental()) +``` + +Works with any grid triangulation. Uses 9 binary variables for triangle +selection. + +### K1 + +```julia +K1() +K1(Incremental()) +``` + +An efficient method that uses only 2 binary variables for triangle selection, but +requires the grid to use a K1 triangulation (set `pattern = :K1` when +constructing the bivariate function, which is the default). + +### UnionJack + +```julia +UnionJack() +UnionJack(Incremental()) +``` + +Uses only 1 binary variable for triangle selection, but requires a Union Jack +triangulation (set `pattern = :UnionJack` when constructing the bivariate +function). + +## Multivariate methods + +These methods work for piecewise linear functions of any dimension. + +### ConvexCombination + +```julia +ConvexCombination() +``` + +A general-purpose method based on the convex combination (or "lambda") formulation. +Uses one binary variable per segment. Works for any dimensionality and can also +be used as a univariate method. + +### DisaggregatedLogarithmic + +```julia +DisaggregatedLogarithmic() +``` + +A disaggregated version of the logarithmic formulation that works for general +piecewise linear functions of any dimension. Uses a logarithmic number of binary +variables in the number of segments. + +### MultipleChoice + +```julia +MultipleChoice() +``` + +A multiple choice formulation for piecewise linear functions specified in +hyperplane representation (i.e., each segment is described by affine functions +and linear constraints defining its domain). This is the only method that +supports the [`PWLFunction`](@ref) with +[`SegmentHyperplaneRep`](@ref PiecewiseLinearOpt.SegmentHyperplaneRep) segments. + +## Choosing a method + +For **univariate** problems, the default [`Logarithmic`](@ref) method is +generally a good choice. It provides a compact formulation with a logarithmic +number of binary variables. If the solver supports native SOS2 constraints, +[`NativeSOS2`](@ref) may be faster. + +For **bivariate** problems, the default [`SixStencil`](@ref) works with any +triangulation. If you can use a specific triangulation pattern, [`K1`](@ref) +or [`UnionJack`](@ref) provide more compact formulations. + +For **general multivariate** problems, [`DisaggregatedLogarithmic`](@ref) or +[`ConvexCombination`](@ref) are the main options, with +[`DisaggregatedLogarithmic`](@ref) being more compact. diff --git a/docs/src/tutorial.md b/docs/src/tutorial.md new file mode 100644 index 0000000..c8a8057 --- /dev/null +++ b/docs/src/tutorial.md @@ -0,0 +1,222 @@ +# Tutorial + +This tutorial demonstrates how to use PiecewiseLinearOpt.jl to model piecewise +linear functions in JuMP optimization models. + +## Univariate piecewise linear functions + +The simplest use case is approximating a univariate function with a piecewise +linear function. You specify a set of breakpoints and either a function or +explicit function values at those breakpoints. + +### Using a function + +```@example univariate +using JuMP, PiecewiseLinearOpt, HiGHS + +model = Model(HiGHS.Optimizer) +set_silent(model) +@variable(model, x) + +# Approximate x^2 on [0, 2] with breakpoints at 0, 0.5, 1.0, 1.5, 2.0 +z = piecewiselinear(model, x, 0:0.5:2, xi -> xi^2) + +@objective(model, Min, z) +optimize!(model) +println("x = $(value(x)), z = $(value(z))") +``` + +### Using explicit values + +You can also pass the function values directly instead of a function: + +```@example univariate +model = Model(HiGHS.Optimizer) +set_silent(model) +@variable(model, x) + +breakpoints = [0.0, 1.0, 2.0, 3.0] +values = [0.0, 1.0, 4.0, 9.0] # x^2 at each breakpoint +z = piecewiselinear(model, x, breakpoints, values) + +@objective(model, Max, z) +@constraint(model, x <= 2.5) +optimize!(model) +println("x = $(value(x)), z = $(value(z))") +``` + +### Choosing a formulation method + +PiecewiseLinearOpt.jl supports several MIP formulation methods. The default for +univariate functions is [`Logarithmic`](@ref) (an alias for +[`LogarithmicEmbedding`](@ref)), which uses ``\log_2(n)`` binary variables for +``n`` breakpoints. You can choose a different method using the `method` keyword: + +```@example univariate +model = Model(HiGHS.Optimizer) +set_silent(model) +@variable(model, x) + +# Use the Incremental formulation (n-1 binary variables) +z = piecewiselinear(model, x, 0:0.5:2, xi -> xi^2; method = Incremental()) + +@objective(model, Min, z) +optimize!(model) +println("x = $(value(x)), z = $(value(z))") +``` + +See [Formulation Methods](@ref) for a full list of available methods. + +### Epigraph and hypograph constraints + +By default, `piecewiselinear` models the graph of the function, i.e., it +enforces `z == f(x)`. You can relax this to an epigraph (`z >= f(x)`) or +hypograph (`z <= f(x)`) constraint using the `direction` keyword: + +```@example univariate +model = Model(HiGHS.Optimizer) +set_silent(model) +@variable(model, x) + +# Epigraph: z >= f(x) +z = piecewiselinear( + model, x, 0:0.5:2, xi -> xi^2; + direction = PiecewiseLinearOpt.Epigraph, +) + +@objective(model, Min, z) # minimize z subject to z >= f(x) +optimize!(model) +println("x = $(value(x)), z = $(value(z))") +``` + +The three options are: +- `PiecewiseLinearOpt.Graph` (default): `z == f(x)` +- `PiecewiseLinearOpt.Epigraph`: `z >= f(x)` +- `PiecewiseLinearOpt.Hypograph`: `z <= f(x)` + +## Bivariate piecewise linear functions + +PiecewiseLinearOpt.jl also supports bivariate piecewise linear functions defined +on triangulated rectangular grids. + +```@example bivariate +using JuMP, PiecewiseLinearOpt, HiGHS + +model = Model(HiGHS.Optimizer) +set_silent(model) +@variable(model, x) +@variable(model, y) + +# Approximate exp(x + y) on [0,1] × [0,1] +z = piecewiselinear( + model, x, y, + 0:0.25:1, 0:0.25:1, + (xi, yi) -> exp(xi + yi), +) + +@objective(model, Min, z) +optimize!(model) +println("x = $(value(x)), y = $(value(y)), z = $(value(z))") +``` + +### Triangulation patterns + +The `pattern` keyword controls how each rectangular cell in the grid is split +into two triangles. The available patterns are: + +- `:K1` (default): A standard triangulation compatible with the [`K1`](@ref) method. +- `:UnionJack`: A triangulation compatible with the [`UnionJack`](@ref) method + where the diagonal direction alternates. +- `:BestFit`: Chooses the diagonal that best approximates the function at the + midpoint. +- `:Upper`: Chooses the diagonal that overestimates the function. +- `:Lower`: Chooses the diagonal that underestimates the function. +- `:Random`: Randomly chooses the diagonal direction. + +```@example bivariate +model = Model(HiGHS.Optimizer) +set_silent(model) +@variable(model, x) +@variable(model, y) + +z = piecewiselinear( + model, x, y, + 0:0.25:1, 0:0.25:1, + (xi, yi) -> xi^2 + yi^2; + pattern = :BestFit, + method = SixStencil(), +) + +@objective(model, Min, z) +optimize!(model) +println("x = $(value(x)), y = $(value(y)), z = $(value(z))") +``` + +## Using the general API + +For more control, you can construct a [`PWLFunction`](@ref) object directly and +pass it to [`piecewiselinear`](@ref) with a tuple of input variables. + +### Univariate + +```@example general +using JuMP, PiecewiseLinearOpt, HiGHS + +# Construct a UnivariatePWLFunction explicitly +pwl = UnivariatePWLFunction([0.0, 1.0, 2.0, 3.0], [0.0, 1.0, 0.5, 1.5]) + +model = Model(HiGHS.Optimizer) +set_silent(model) +@variable(model, x) + +# Pass input as a tuple +output = piecewiselinear(model, (x,), pwl) +# output is a 1-tuple; extract the variable +z = output[1] + +@objective(model, Max, z) +optimize!(model) +println("x = $(value(x)), z = $(value(z))") +``` + +### Bivariate + +```@example general +pwl = BivariatePWLFunction( + 0:0.5:1, 0:0.5:1, + (xi, yi) -> xi * yi; + pattern = :K1, +) + +model = Model(HiGHS.Optimizer) +set_silent(model) +@variable(model, x) +@variable(model, y) + +output = piecewiselinear(model, (x, y), pwl) +z = output[1] + +@objective(model, Max, z) +optimize!(model) +println("x = $(value(x)), y = $(value(y)), z = $(value(z))") +``` + +## Providing an existing output variable + +If you already have a variable you want to use as the output, you can pass it +using the `output_var` (univariate/bivariate) or `output_vars` (general) +keyword: + +```@example univariate +model = Model(HiGHS.Optimizer) +set_silent(model) +@variable(model, x) +@variable(model, z) + +piecewiselinear(model, x, 0:0.5:2, xi -> xi^2; output_var = z) + +@objective(model, Max, z) +@constraint(model, x <= 1.5) +optimize!(model) +println("x = $(value(x)), z = $(value(z))") +``` diff --git a/src/PiecewiseLinearOpt.jl b/src/PiecewiseLinearOpt.jl index 4327bec..db4da52 100644 --- a/src/PiecewiseLinearOpt.jl +++ b/src/PiecewiseLinearOpt.jl @@ -8,12 +8,96 @@ 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 + +""" + _next_pwl_id!(model::JuMP.Model) -> Int + +Atomically initialize the PWL extension (if needed), increment the counter, +and return the new unique ID. Call this exactly once per `piecewiselinear()` +invocation.""" +function _next_pwl_id!(model::JuMP.Model) + initPWL!(model) + model.ext[:PWL].counter += 1 + return model.ext[:PWL].counter +end + +""" + _pwl_name(model::JuMP.Model, varname::String) -> String + +Return a prefixed variable name like `"pwl3_λ"` using the current PWL counter. +Use this in all formulation methods to generate descriptive variable names. +""" +function _pwl_name(model::JuMP.Model, varname::String) + return "pwl$(model.ext[:PWL].counter)_$(varname)" +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 +""" + Logarithmic + +Alias for [`LogarithmicEmbedding`](@ref). The default formulation method for +univariate piecewise linear functions. +""" +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..5fecf46 --- /dev/null +++ b/src/methods/bivariate/common.jl @@ -0,0 +1,102 @@ +# 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} + 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 = _pwl_name(model, "λ") + ) + 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..a55e3be --- /dev/null +++ b/src/methods/bivariate/k1.jl @@ -0,0 +1,60 @@ +# 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. + +""" + K1() + K1(axis_method::Method) + +Bivariate formulation for K1-triangulated grids. + +Uses 2 binary variables for triangle selection plus the binary variables from +`axis_method` for the axis-aligned SOS2 constraints. Requires the piecewise +linear function to use a K1 triangulation (`pattern = :K1`). + +The default `axis_method` is [`Logarithmic`](@ref). +""" +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) + # TODO: Certify triangulation is "uniform" + + w = JuMP.@variable(model, [1:2], Bin, base_name = _pwl_name(model, "w")) + 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..ab3566a --- /dev/null +++ b/src/methods/bivariate/nine_stencil.jl @@ -0,0 +1,84 @@ +# 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. + +""" + NineStencil() + NineStencil(axis_method::Method) + +Bivariate formulation using a nine-stencil biclique cover for triangle selection. + +Uses 9 binary variables for triangle selection plus the binary variables from +`axis_method` for the axis-aligned SOS2 constraints. Works with any grid +triangulation pattern. + +The default `axis_method` is [`Logarithmic`](@ref). +""" +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) + w = JuMP.@variable( + model, + [1:3, 1:3], + Bin, + base_name = _pwl_name(model, "w") + ) + 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..a2417fb --- /dev/null +++ b/src/methods/bivariate/optimal_independent_branching.jl @@ -0,0 +1,175 @@ +# 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} + 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 = _pwl_name(model, "λ") + ) + 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 = _pwl_name(model, "z")) + + 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..e34e55e --- /dev/null +++ b/src/methods/bivariate/optimal_triangle_selection.jl @@ -0,0 +1,116 @@ +# 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, +) + 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 = _pwl_name(model, "z")) + + 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..9c2564a --- /dev/null +++ b/src/methods/bivariate/six_stencil.jl @@ -0,0 +1,151 @@ +# 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. + +""" + SixStencil() + SixStencil(axis_method::Method) + +Bivariate formulation using a six-stencil biclique cover for triangle selection. + +Uses 6 binary variables for triangle selection plus the binary variables from +`axis_method` for the axis-aligned SOS2 constraints. Works with any grid +triangulation pattern. This is the default method for bivariate functions. + +The default `axis_method` is [`Logarithmic`](@ref). +""" +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) + # diagonal lines running from SW to NE. Grouped with an offset of 3. + w_sw_ne = JuMP.@variable( + model, + [0:2], + Bin, + base_name = _pwl_name(model, "w_sw_ne") + ) + 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 = _pwl_name(model, "w_se_nw") + ) + 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..07be517 --- /dev/null +++ b/src/methods/bivariate/union_jack.jl @@ -0,0 +1,53 @@ +# 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. + +""" + UnionJack() + UnionJack(axis_method::Method) + +Bivariate formulation for Union Jack-triangulated grids. + +Uses only 1 binary variable for triangle selection plus the binary variables +from `axis_method` for the axis-aligned SOS2 constraints. Requires the +piecewise linear function to use a Union Jack triangulation +(`pattern = :UnionJack`). + +The default `axis_method` is [`Logarithmic`](@ref). +""" +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) + # 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 = _pwl_name(model, "w")) + + 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..b759ac5 --- /dev/null +++ b/src/methods/multivariate/convex_combination.jl @@ -0,0 +1,97 @@ +# 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. + +""" + ConvexCombination() + +The convex combination ("lambda") formulation for piecewise linear functions. + +Uses one binary variable per segment. Works for any input dimension. Can also +be used as a univariate SOS2 formulation method. +""" +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 + segments = pwl.segments + z = JuMP.@variable( + model, + [segments], + Bin, + base_name = _pwl_name(model, "z") + ) + 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 = _pwl_name(model, "λ") + ) + 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} + n = length(λ) + z = JuMP.@variable(model, [1:n-1], Bin, base_name = _pwl_name(model, "z")) + 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..bdb63f9 --- /dev/null +++ b/src/methods/multivariate/disaggregated_logarithmic.jl @@ -0,0 +1,70 @@ +# 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. + +""" + DisaggregatedLogarithmic() + +The disaggregated logarithmic formulation for piecewise linear functions. + +Uses `⌈log₂(S)⌉` binary variables where `S` is the number of segments. +Works for any input dimension. More compact than [`ConvexCombination`](@ref) +for functions with many segments. +""" +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} + 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 = _pwl_name(model, "γ") + ) + 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 = _pwl_name(model, "z")) + 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..8090da0 --- /dev/null +++ b/src/methods/multivariate/multiple_choice.jl @@ -0,0 +1,64 @@ +# 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. + +""" + MultipleChoice() + +The multiple choice formulation for piecewise linear functions in hyperplane +representation. + +This is the only method that supports `PWLFunction` objects with +[`SegmentHyperplaneRep`](@ref PiecewiseLinearOpt.SegmentHyperplaneRep) segments, +where each segment is defined by affine constraints and affine output functions. +""" +struct MultipleChoice <: Method end + +function formulate_pwl!( + model::JuMP.Model, + input_vars::NTuple{D,VarOrAff}, + output_vars::NTuple{F,VarOrAff}, + pwl::PWLFunctionHyperplaneRep{D,F}, + method::MultipleChoice, + direction::DIRECTION, +) where {D,F} + segments = pwl.segments + S = 1:length(segments) + x_hat = + JuMP.@variable(model, [S, 1:D], base_name = _pwl_name(model, "x_hat")) + y_hat = + JuMP.@variable(model, [S, 1:F], base_name = _pwl_name(model, "y_hat")) + z = JuMP.@variable(model, [S], Bin, base_name = _pwl_name(model, "z")) + JuMP.@constraint(model, sum(z) == 1) + for i in 1:D + JuMP.@constraint(model, sum(x_hat[:, i]) == input_vars[i]) + end + for i in 1:F + _constrain_output_var( + model, + output_vars[i], + sum(y_hat[:, i]), + direction, + ) + end + for (s, seg) in enumerate(segments) + for constraint in seg.constraints + coeffs, offset = constraint.coeffs, constraint.offset + JuMP.@constraint( + model, + LinearAlgebra.dot(coeffs, x_hat[s, :]) + offset * z[s] ≥ 0 + ) + end + for i in 1:F + output_func = seg.funcs[i] + coeffs, offset = output_func.coeffs, output_func.offset + JuMP.@constraint( + model, + y_hat[s, i] == + LinearAlgebra.dot(coeffs, x_hat[s, :]) + offset * z[s] + ) + 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..3961a75 --- /dev/null +++ b/src/methods/univariate/incremental.jl @@ -0,0 +1,54 @@ +# 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 +""" + Incremental() + +The incremental MIP formulation for univariate piecewise linear functions. + +Uses `n - 1` binary variables for a function with `n` breakpoints. The formulation +is based on an incremental representation where each binary variable indicates +whether the function has "passed" a given breakpoint. +""" +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 + + n = length(pwl.segments) + 1 + @assert length(xs) == length(ys) == n + + δ = JuMP.@variable( + model, + [1:n], + lower_bound = 0, + upper_bound = 1, + base_name = _pwl_name(model, "δ") + ) + z = JuMP.@variable(model, [1:n-1], Bin, base_name = _pwl_name(model, "z")) + 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..bf3f5a7 --- /dev/null +++ b/src/methods/univariate/logarithmic_embedding.jl @@ -0,0 +1,41 @@ +# 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. + +""" + LogarithmicEmbedding() + +The logarithmic embedding formulation for SOS2 constraints. + +Uses `⌈log₂(n-1)⌉` binary variables for `n` breakpoints, based on reflected +Gray codes. This is the default method for univariate piecewise linear functions. + +`Logarithmic` is an alias for `LogarithmicEmbedding`. +""" +struct LogarithmicEmbedding <: Method end + +function formulate_sos2!( + model::JuMP.Model, + λ::Vector{T}, + method::LogarithmicEmbedding, +) where {T<:VarOrAff} + 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 = _pwl_name(model, "y")) + _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..e4ca907 --- /dev/null +++ b/src/methods/univariate/logarithmic_independent_branching.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. + +""" + LogarithmicIndependentBranching() + +The logarithmic independent branching formulation for SOS2 constraints. + +Uses `⌈log₂(n-1)⌉` binary variables for `n` breakpoints. Similar to +[`LogarithmicEmbedding`](@ref) but uses an independent branching scheme. +""" +struct LogarithmicIndependentBranching <: Method end + +function formulate_sos2!( + model::JuMP.Model, + λ::Vector{T}, + method::LogarithmicIndependentBranching, +) where {T<:VarOrAff} + 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 = _pwl_name(model, "y")) + _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..87387a0 --- /dev/null +++ b/src/methods/univariate/native_sos2.jl @@ -0,0 +1,23 @@ +# 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. + +""" + NativeSOS2() + +Delegates the SOS2 constraint directly to the solver using JuMP's built-in +SOS2 support. No additional binary variables are introduced. + +Requires a solver that supports SOS2 constraints natively. +""" +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..a5b3fe0 --- /dev/null +++ b/src/methods/univariate/zig_zag_binary.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. + +""" + ZigZagBinary() + +The zig-zag binary encoding formulation for SOS2 constraints. + +Uses `⌈log₂(n-1)⌉` binary variables for `n` breakpoints with a zig-zag +encoding scheme. +""" +struct ZigZagBinary <: Method end + +function formulate_sos2!( + model::JuMP.Model, + λ::Vector{T}, + method::ZigZagBinary, +) where {T<:VarOrAff} + 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 = _pwl_name(model, "y")) + _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..aef91a0 --- /dev/null +++ b/src/methods/univariate/zig_zag_integer.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. + +""" + ZigZagInteger() + +The zig-zag integer encoding formulation for SOS2 constraints. + +Uses `⌈log₂(n-1)⌉` general integer variables for `n` breakpoints. A variant +of [`ZigZagBinary`](@ref) that uses integer variables instead of binary. +""" +struct ZigZagInteger <: Method end + +function formulate_sos2!( + model::JuMP.Model, + λ::Vector{T}, + method::ZigZagInteger, +) where {T<:VarOrAff} + 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 = _pwl_name(model, "y") + ) + _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..cbf90ef --- /dev/null +++ b/src/methods/util.jl @@ -0,0 +1,293 @@ +# 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} + λ = similar(Array{JuMP.VariableRef}, axes(grid.input_vals)) + λ_name = _pwl_name(model, "λ") + for I in eachindex(λ) + λ[I] = JuMP.@variable( + model, + lower_bound = 0, + upper_bound = 1, + base_name = λ_name + ) + 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) + # 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}) + 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 diff --git a/src/pwlinear.jl b/src/pwlinear.jl new file mode 100644 index 0000000..63fb4f7 --- /dev/null +++ b/src/pwlinear.jl @@ -0,0 +1,195 @@ +# 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 methods for dimensions 1 and 2 +_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} + _next_pwl_id!(model) + + 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 = _pwl_name(model, "y") + )..., + ) + end + + formulate_pwl!(model, input_vars, output_vars, pwl, method, direction) + return output_vars +end + +function piecewiselinear( + model::JuMP.Model, + input_vars::NTuple{D,VarOrAff}, + pwl::PWLFunction{D,F,SegmentHyperplaneRep{D,F}}; + method::Method = MultipleChoice(), + direction::DIRECTION = Graph, + output_vars::Union{Nothing,NTuple{F,VarOrAff}} = nothing, +) where {D,F} + _next_pwl_id!(model) + + if isempty(pwl.segments) + error( + "I don't know how to handle a piecewise linear function with no breakpoints.", + ) + end + + if output_vars === nothing + output_vars = tuple( + JuMP.@variable( + model, + [i in 1:F], + base_name = _pwl_name(model, "y") + )..., + ) + 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..eb8828d 100644 --- a/src/types.jl +++ b/src/types.jl @@ -3,131 +3,243 @@ # 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 + +""" + DIRECTION + +Enum controlling how the output variable relates to the piecewise linear function value. + +- `Graph`: The output equals the function value, i.e., `z == f(x)`. +- `Epigraph`: The output is an upper bound, i.e., `z >= f(x)`. +- `Hypograph`: The output is a lower bound, i.e., `z <= f(x)`. +""" +@enum DIRECTION Graph Epigraph Hypograph + +# TODO: Make eltypes of input_vals and output_vals a type parameter +abstract type Segment{D,F} end + +""" + SegmentPointRep{D,F} + +A segment of a piecewise linear function represented by its vertices (point representation). +Each segment is a simplex defined by `D + 1` vertices in `D`-dimensional input space. -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) +# Fields +- `input_vals::Vector{NTuple{D,Float64}}`: The input coordinates of the segment vertices. +- `output_vals::Vector{NTuple{F,Float64}}`: The output values at each vertex. +""" +struct SegmentPointRep{D,F} <: Segment{D,F} + input_vals::Vector{NTuple{D,Float64}} + output_vals::Vector{NTuple{F,Float64}} + + 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()) +""" + AffineFunction{D} + +An affine function `f(x) = coeffs ⋅ x + offset` in `D` dimensions. +""" +struct AffineFunction{D} + coeffs::NTuple{D,Float64} + offset::Float64 +end -const UnivariatePWLFunction = PWLFunction{1} +""" + SegmentHyperplaneRep{D,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], - ) +A segment of a piecewise linear function represented by affine constraints (hyperplane +representation). The domain of the segment is defined by `f_i(x) >= 0` for each +constraint, and the function values are given by affine functions. + +# Fields +- `constraints::Vector{AffineFunction{D}}`: Linear constraints defining the segment domain. +- `funcs::NTuple{F,AffineFunction{D}}`: Affine functions giving the output on this segment. +""" +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 -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], - ) +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 + +""" + PWLFunction{D,F,T} + +A piecewise linear function from ℝ^D to ℝ^F, composed of a collection of +segments of type `T`. + +`D` is the input dimension, `F` is the output dimension, and `T` is the +segment representation type (either [`SegmentPointRep`](@ref PiecewiseLinearOpt.SegmentPointRep) +or [`SegmentHyperplaneRep`](@ref PiecewiseLinearOpt.SegmentHyperplaneRep)). + +# Fields +- `segments::Vector{T}`: The segments composing the piecewise linear function. +- `structure::SegmentStructure{D}`: Metadata about the structure of the segments (e.g., grid triangulation type). + +# Type aliases +- `UnivariatePWLFunction = PWLFunction{1,1,SegmentPointRep{1,1}}`: A univariate PWL function. +- `BivariatePWLFunction = PWLFunction{2,1,SegmentPointRep{2,1}}`: A bivariate PWL function. +""" +struct PWLFunction{D,F,T<:Segment{D,F}} + segments::Vector{T} + structure::SegmentStructure{D} end -const BivariatePWLFunction = PWLFunction{2} +const PWLFunctionPointRep{D,F} = PWLFunction{D,F,SegmentPointRep{D,F}} +const PWLFunctionHyperplaneRep{D,F} = PWLFunction{D,F,SegmentHyperplaneRep{D,F}} + +#const UnivariatePWLFunction{F} = PWLFunctionPointRep{1, F} +#const BivariatePWLFunction{F} = PWLFunctionPointRep{2, F} + +""" + UnivariatePWLFunction + +A univariate piecewise linear function (ℝ → ℝ) in point representation. + +# Constructors + + UnivariatePWLFunction(x::Vector, z::Vector) + +Construct from breakpoints `x` and corresponding function values `z`. + + UnivariatePWLFunction(x, f::Function) + +Construct from breakpoints `x` and a function `f` evaluated at those points. -function BivariatePWLFunction( +# Examples +```jldoctest +julia> pwl = UnivariatePWLFunction([0.0, 1.0, 2.0], [0.0, 1.0, 0.0]); + +julia> pwl = UnivariatePWLFunction(0:0.5:2, x -> x^2); +``` +""" +const UnivariatePWLFunction = PWLFunctionPointRep{1,1} + +""" + BivariatePWLFunction + +A bivariate piecewise linear function (ℝ² → ℝ) in point representation, +defined on a triangulated rectangular grid. + +# Constructor + + BivariatePWLFunction(x, y, fz::Function; pattern=:K1, seed=...) + +Construct from grid breakpoints `x` and `y`, and a function `fz(xi, yi)`. +The `pattern` keyword controls how grid rectangles are triangulated: +- `:K1` (default), `:UnionJack`, `:BestFit`, `:Upper`, `:Lower`, `:Random` + +# Examples +```jldoctest +julia> pwl = BivariatePWLFunction(0:0.5:1, 0:0.5:1, (x, y) -> x + y); +``` +""" +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 PWLFunctionPointRep{1,1}(x, f::Function) + d = collect(x) + fd = [f(xi) for xi in d] + return PWLFunctionPointRep{1,1}(d, fd) +end + +function PWLFunctionPointRep{1,1}(x, z) + return PWLFunctionPointRep{1,1}(collect(x), collect(z)) +end + +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/Project.toml b/test/Project.toml new file mode 100644 index 0000000..8dae788 --- /dev/null +++ b/test/Project.toml @@ -0,0 +1,8 @@ +[deps] +Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" +HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b" +JuMP = "4076af6c-e467-56ae-b986-b466b2749572" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" +PiecewiseLinearOpt = "0f51c51e-adfa-5141-8a04-d40246b8977c" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" diff --git a/test/aqua.jl b/test/aqua.jl new file mode 100644 index 0000000..cf87a25 --- /dev/null +++ b/test/aqua.jl @@ -0,0 +1,3 @@ +using PiecewiseLinearOpt +using Aqua +Aqua.test_all(PiecewiseLinearOpt) diff --git a/test/runtests.jl b/test/runtests.jl index 3c89df5..ef3d86f 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -4,135 +4,487 @@ # in the LICENSE.md file or at https://opensource.org/licenses/MIT. using PiecewiseLinearOpt +using HiGHS +using JuMP +using LinearAlgebra using Test +const PLO = PiecewiseLinearOpt -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, +optimizer = optimizer_with_attributes(HiGHS.Optimizer, MOI.Silent() => true) + +const methods_1D = [ + ConvexCombination(), + DisaggregatedLogarithmic(), + Incremental(), + LogarithmicEmbedding(), + LogarithmicIndependentBranching(), + NativeSOS2(), + ZigZagBinary(), + ZigZagInteger(), +] +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 is tested in a dedicated testset below with a + # coarser grid to keep CI fast (the selection MIP is expensive). + [SixStencil(sos2_method) for sos2_method in methods_1D]..., +] + +@testset "PiecewiseLinearOpt.jl" begin + @testset "Simple univariate" begin + @testset "$method" 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 + 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 + + @testset "Simple bivariate" begin + @testset "$method" 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 + end + + @testset "1D" begin + @testset "$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 + end + + @testset "2D" begin + 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 "$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 + end + + # Test OptimalTriangleSelection with multiple SOS2 methods on a small problem + # (the small 2-triangle problem is fast; the expensive part is the grid-based test) + @testset "OptimalTriangleSelection" begin + @testset "$sos2" for sos2 in [ + ConvexCombination(), + LogarithmicEmbedding(), + NativeSOS2(), + ] + method = OptimalTriangleSelection(optimizer, sos2) + 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 + end + + # Test OptimalTriangleSelection on a grid-based 2D problem with a coarser grid + # to keep CI fast (the selection MIP scales poorly with grid size). + @testset "OptimalTriangleSelection 2D grid" begin + @testset "$pattern" for pattern in [:Upper, :K1] + method = OptimalTriangleSelection(optimizer, ConvexCombination()) + model = Model(optimizer) + @variable(model, x[1:2]) + # Use a coarser grid (length=4) vs the main 2D tests (length=8) + d = range(0; stop = 1, length = 4) + 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(z) ≈ objective_value(model) rtol = 1e-3 + end + end + + # Test OptimalIndependentBranching (small problem only, the MIP can be very slow) + @testset "OptimalIndependentBranching" begin + method = OptimalIndependentBranching(optimizer) + 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,)], ) - 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 + 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 end -end -@testset "Bivariate tests " begin - @testset "2D: $method, $pattern" for method in methods_2D, - pattern in patterns_2D + # Test MultipleChoice (requires SegmentHyperplaneRep) + @testset "MultipleChoice" begin + method = MultipleChoice() + + # 1D: Two segments: f(x) = x + 1.5 on [1,2], f(x) = -2.5x + 8.5 on [2,3] + # Segment 1: domain x >= 1 and x <= 2, output y = x + 1.5 + s1 = PLO.SegmentHyperplaneRep{1,1}( + [ + PLO.AffineFunction{1}((1.0,), -1.0), + PLO.AffineFunction{1}((-1.0,), 2.0), + ], + (PLO.AffineFunction{1}((1.0,), 1.5),), + ) + # Segment 2: domain x >= 2 and x <= 3, output y = -2.5x + 8.5 + s2 = PLO.SegmentHyperplaneRep{1,1}( + [ + PLO.AffineFunction{1}((1.0,), -2.0), + PLO.AffineFunction{1}((-1.0,), 3.0), + ], + (PLO.AffineFunction{1}((-2.5,), 8.5),), + ) + pwl = PLO.PWLFunction{1,1,PLO.SegmentHyperplaneRep{1,1}}( + [s1, s2], + PLO.Intervals(), + ) + + model = Model(optimizer) + @variable(model, 1 <= x <= 3) + y = piecewiselinear(model, (x,), pwl; method = method) + @objective(model, Min, y[1]) + optimize!(model) + @test termination_status(model) == MOI.OPTIMAL + # Min of y = x + 1.5 on [1,2] is 2.5 at x=1 + # Min of y = -2.5x + 8.5 on [2,3] is 1.0 at x=3 + @test value(y[1]) ≈ 1.0 rtol = 1e-4 + @test value(x) ≈ 3.0 rtol = 1e-4 + end + + # Test direction parameter (Epigraph and Hypograph) + @testset "Direction tests" begin + 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()) + + model = Model(optimizer) + @variable(model, x) + y_epi = piecewiselinear(model, (x,), pwl; direction = PLO.Epigraph) + @objective(model, Min, y_epi[1]) + optimize!(model) + @test termination_status(model) == MOI.OPTIMAL + + model = Model(optimizer) + @variable(model, x) + y_hypo = piecewiselinear(model, (x,), pwl; direction = PLO.Hypograph) + @objective(model, Max, y_hypo[1]) + optimize!(model) + @test termination_status(model) == MOI.OPTIMAL + end + + # Test error cases + @testset "Error cases" begin + # Mismatched input/output in SegmentPointRep + @test_throws ErrorException PLO.SegmentPointRep{1,1}( + [(1.0,), (2.0,)], + [(2.5,)], + ) + + # Empty segments error + empty_pwl = PLO.PWLFunction(PLO.SegmentPointRep{1,1}[], PLO.Intervals()) + model = Model(optimizer) + @variable(model, x) + @test_throws ErrorException piecewiselinear(model, (x,), empty_pwl) + + # Test output_vars kwarg + model = Model(optimizer) + @variable(model, x) + @variable(model, y_out) + 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; output_vars = (y_out,)) + @test y[1] === y_out + + # Test univariate output_var kwarg with function + model = Model(optimizer) + @variable(model, x) + @variable(model, y_out2) + d = 0:0.1:1 + f = xi -> xi^2 + y = piecewiselinear(model, x, d, f; output_var = y_out2) + @test y === y_out2 + + # Test univariate output_var with fd (vector) version + model = Model(optimizer) + @variable(model, x) + @variable(model, y_out3) + fd = [xi^2 for xi in d] + y = piecewiselinear(model, x, d, fd; output_var = y_out3) + @test y === y_out3 - 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( + # Test bivariate output_var kwarg + model = Model(optimizer) + @variable(model, x) + @variable(model, y) + @variable(model, z_out) + d2 = 0:0.5:1 + f2 = (xi, yi) -> xi + yi + z = piecewiselinear(model, x, y, d2, d2, f2; output_var = z_out) + @test z === z_out + + # Test formulate_pwl! fallback error + struct FakeMethod <: PLO.Method end + 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()) + @test_throws MethodError piecewiselinear( model, - x[1], - x[2], - BivariatePWLFunction(d, d, f; pattern = pattern); - method = method, + (x,), + pwl; + method = FakeMethod(), ) - 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 -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 + # Test internal utility functions + @testset "Utility functions" begin + # Test _canonical! + v = [3.0, 4.0] + result = PLO._canonical!(v) + @test result[1] ≈ 0.6 atol = 1e-6 + @test result[2] ≈ 0.8 atol = 1e-6 + + # Test with negative leading element + v2 = [-3.0, 4.0] + result2 = PLO._canonical!(v2) + @test result2[1] > 0 # should be positive after sign flip + + # Test with near-zero elements + v3 = [1e-10, 1.0] + result3 = PLO._canonical!(v3) + @test result3[1] == 0.0 + @test result3[2] ≈ 1.0 atol = 1e-6 + + # Test _compute_hyperplanes with k=2 + C = [[0, 0], [1, 0], [1, 1], [0, 1]] + hyperplanes = PLO._compute_hyperplanes(C) + @test length(hyperplanes) >= 1 + + # Test _compute_hyperplanes with k > 2 + C3 = [[0, 0, 0], [1, 0, 0], [1, 1, 0], [1, 1, 1], [0, 1, 1]] + hyperplanes3 = PLO._compute_hyperplanes(C3) + @test length(hyperplanes3) >= 1 + + # Test _compute_hyperplanes error with k <= 1 + @test_throws ErrorException PLO._compute_hyperplanes([[1], [2], [3]]) + + # Test _continuous_gridpoints_or_die error: discontinuous function + s1 = PLO.SegmentPointRep{1,1}([(1.0,), (2.0,)], [(2.5,), (3.5,)]) + s2 = PLO.SegmentPointRep{1,1}([(2.0,), (3.0,)], [(999.0,), (1.0,)]) + pwl_disc = PLO.PWLFunction([s1, s2], PLO.Intervals()) + @test_throws ErrorException PLO._continuous_gridpoints_or_die(pwl_disc) + + # Test Grid size mismatch error + input_vals = Array{NTuple{1,Float64},1}([(1.0,), (2.0,)]) + output_vals = Array{NTuple{1,Float64},1}([(1.0,)]) + @test_throws ErrorException PLO.Grid(input_vals, output_vals) + end + + include("aqua.jl") +end