From f85a8fc267fe2b82bb8c8c6127ddb85c369abe1f Mon Sep 17 00:00:00 2001 From: Valentin Kaisermayer Date: Tue, 31 Mar 2026 10:15:47 +0200 Subject: [PATCH 1/3] adds docs and tests - adds aqua tests - improves coverage - adds docs - updates CI workflow for docs (not sure if it works) --- .github/workflows/ci.yml | 72 ++++-- .gitignore | 2 + README.md | 1 + docs/Project.toml | 5 + docs/make.jl | 34 +++ docs/src/api.md | 58 +++++ docs/src/index.md | 67 ++++++ docs/src/methods.md | 199 ++++++++++++++++ docs/src/tutorial.md | 222 ++++++++++++++++++ src/PiecewiseLinearOpt.jl | 6 + src/methods/bivariate/k1.jl | 12 + src/methods/bivariate/nine_stencil.jl | 12 + src/methods/bivariate/six_stencil.jl | 12 + src/methods/bivariate/union_jack.jl | 13 + .../multivariate/convex_combination.jl | 8 + .../multivariate/disaggregated_logarithmic.jl | 9 + src/methods/multivariate/multiple_choice.jl | 40 +++- src/methods/univariate/incremental.jl | 9 + .../univariate/logarithmic_embedding.jl | 10 + .../logarithmic_independent_branching.jl | 8 + src/methods/univariate/native_sos2.jl | 8 + src/methods/univariate/zig_zag_binary.jl | 8 + src/methods/univariate/zig_zag_integer.jl | 8 + src/methods/util.jl | 8 +- src/pwlinear.jl | 33 +++ src/types.jl | 95 ++++++++ test/aqua.jl | 3 + test/runtests.jl | 187 +++++++++++++++ 28 files changed, 1115 insertions(+), 34 deletions(-) create mode 100644 docs/Project.toml create mode 100644 docs/make.jl create mode 100644 docs/src/api.md create mode 100644 docs/src/index.md create mode 100644 docs/src/methods.md create mode 100644 docs/src/tutorial.md create mode 100644 test/aqua.jl diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 65db331..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.10', '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 9071367..0917646 100644 --- a/.gitignore +++ b/.gitignore @@ -5,3 +5,5 @@ Manifest.toml *.DS_Store .DS_Store .vscode/settings.json +coverage/ +docs/build/ \ No newline at end of file diff --git a/README.md b/README.md index f4f6c19..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 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 71a5f8d..dc6f231 100644 --- a/src/PiecewiseLinearOpt.jl +++ b/src/PiecewiseLinearOpt.jl @@ -48,6 +48,12 @@ include("methods/univariate/zig_zag_integer.jl") 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 diff --git a/src/methods/bivariate/k1.jl b/src/methods/bivariate/k1.jl index 14a62a7..74b5051 100644 --- a/src/methods/bivariate/k1.jl +++ b/src/methods/bivariate/k1.jl @@ -3,6 +3,18 @@ # 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 diff --git a/src/methods/bivariate/nine_stencil.jl b/src/methods/bivariate/nine_stencil.jl index 318dc40..417aac9 100644 --- a/src/methods/bivariate/nine_stencil.jl +++ b/src/methods/bivariate/nine_stencil.jl @@ -3,6 +3,18 @@ # 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 diff --git a/src/methods/bivariate/six_stencil.jl b/src/methods/bivariate/six_stencil.jl index 75c47af..dcd0b98 100644 --- a/src/methods/bivariate/six_stencil.jl +++ b/src/methods/bivariate/six_stencil.jl @@ -3,6 +3,18 @@ # 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 diff --git a/src/methods/bivariate/union_jack.jl b/src/methods/bivariate/union_jack.jl index 23f3452..fdc61ee 100644 --- a/src/methods/bivariate/union_jack.jl +++ b/src/methods/bivariate/union_jack.jl @@ -3,6 +3,19 @@ # 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 diff --git a/src/methods/multivariate/convex_combination.jl b/src/methods/multivariate/convex_combination.jl index 17a588e..a50b023 100644 --- a/src/methods/multivariate/convex_combination.jl +++ b/src/methods/multivariate/convex_combination.jl @@ -3,6 +3,14 @@ # 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!( diff --git a/src/methods/multivariate/disaggregated_logarithmic.jl b/src/methods/multivariate/disaggregated_logarithmic.jl index cdf5c21..e03d6f0 100644 --- a/src/methods/multivariate/disaggregated_logarithmic.jl +++ b/src/methods/multivariate/disaggregated_logarithmic.jl @@ -3,6 +3,15 @@ # 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!( diff --git a/src/methods/multivariate/multiple_choice.jl b/src/methods/multivariate/multiple_choice.jl index 3855dc8..48222d7 100644 --- a/src/methods/multivariate/multiple_choice.jl +++ b/src/methods/multivariate/multiple_choice.jl @@ -3,41 +3,59 @@ # 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::Tuple{VarOrAff}, + input_vars::NTuple{D,VarOrAff}, output_vars::NTuple{F,VarOrAff}, pwl::PWLFunctionHyperplaneRep{D,F}, method::MultipleChoice, direction::DIRECTION, ) where {D,F} - x_hat = JuMP.@variable(model, [segments, 1:D], base_name = "x_hat_$counter") - y_hat = JuMP.@variable(model, [segments, 1:F], base_name = "y_hat_$counter") - z = JuMP.@variable(model, [segments], Bin, base_name = "z_$counter") + counter = model.ext[:PWL].counter + segments = pwl.segments + S = 1:length(segments) + x_hat = JuMP.@variable(model, [S, 1:D], base_name = "x_hat_$counter") + y_hat = JuMP.@variable(model, [S, 1:F], base_name = "y_hat_$counter") + z = JuMP.@variable(model, [S], Bin, base_name = "z_$counter") JuMP.@constraint(model, sum(z) == 1) for i in 1:D - JuMP.@constraint(model, sum(x_hat[:, i]) == x[i]) + JuMP.@constraint(model, sum(x_hat[:, i]) == input_vars[i]) end for i in 1:F - JuMP.@constraint(model, sum(y_hat[:, i]) == y[i]) + _constrain_output_var( + model, + output_vars[i], + sum(y_hat[:, i]), + direction, + ) end - for seg in segments + 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[seg, :]) + offset * z[seg] ≥ 0 + LinearAlgebra.dot(coeffs, x_hat[s, :]) + offset * z[s] ≥ 0 ) end for i in 1:F - output_func = seg.output_funcs[i] + output_func = seg.funcs[i] coeffs, offset = output_func.coeffs, output_func.offset JuMP.@constraint( model, - y_hat[seg, i] == - LinearAlgebra.dot(coeffs, x_hat[seg, :]) + offset * z[seg] + y_hat[s, i] == + LinearAlgebra.dot(coeffs, x_hat[s, :]) + offset * z[s] ) end end diff --git a/src/methods/univariate/incremental.jl b/src/methods/univariate/incremental.jl index 96db5b3..52d5719 100644 --- a/src/methods/univariate/incremental.jl +++ b/src/methods/univariate/incremental.jl @@ -4,6 +4,15 @@ # 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!( diff --git a/src/methods/univariate/logarithmic_embedding.jl b/src/methods/univariate/logarithmic_embedding.jl index 1c3b357..f82fae9 100644 --- a/src/methods/univariate/logarithmic_embedding.jl +++ b/src/methods/univariate/logarithmic_embedding.jl @@ -3,6 +3,16 @@ # 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!( diff --git a/src/methods/univariate/logarithmic_independent_branching.jl b/src/methods/univariate/logarithmic_independent_branching.jl index 4c549bd..95195ac 100644 --- a/src/methods/univariate/logarithmic_independent_branching.jl +++ b/src/methods/univariate/logarithmic_independent_branching.jl @@ -3,6 +3,14 @@ # 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!( diff --git a/src/methods/univariate/native_sos2.jl b/src/methods/univariate/native_sos2.jl index 36d56d3..87387a0 100644 --- a/src/methods/univariate/native_sos2.jl +++ b/src/methods/univariate/native_sos2.jl @@ -3,6 +3,14 @@ # 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!( diff --git a/src/methods/univariate/zig_zag_binary.jl b/src/methods/univariate/zig_zag_binary.jl index 1405185..753bcf3 100644 --- a/src/methods/univariate/zig_zag_binary.jl +++ b/src/methods/univariate/zig_zag_binary.jl @@ -3,6 +3,14 @@ # 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!( diff --git a/src/methods/univariate/zig_zag_integer.jl b/src/methods/univariate/zig_zag_integer.jl index 7e9c6f3..9b37c3f 100644 --- a/src/methods/univariate/zig_zag_integer.jl +++ b/src/methods/univariate/zig_zag_integer.jl @@ -3,6 +3,14 @@ # 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!( diff --git a/src/methods/util.jl b/src/methods/util.jl index 3326ad8..bd5dd56 100644 --- a/src/methods/util.jl +++ b/src/methods/util.jl @@ -175,7 +175,7 @@ function _check_triangle(segment::SegmentHyperplaneRep{D}) where {D} return end -function _check_triangulation(pwl::PWLFunction) where {D,F} +function _check_triangulation(pwl::PWLFunction) # TODO: Add check that segments lie in a single subrectangle for segment in pwl.segments _check_triangle(segment) @@ -207,7 +207,7 @@ function _compute_hyperplanes(C::Vector{Vector{T}}) where {T<:Number} @assert n == 4 spanners = Vector{Float64}[] for i in 1:n-1 - v = canonical!(d[i, :]) + v = _canonical!(d[i, :]) v = [v[2], -v[1]] push!(spanners, v) end @@ -236,7 +236,7 @@ function _compute_hyperplanes(C::Vector{Vector{T}}) where {T<:Number} nullsp = LinearAlgebra.nullspace(d[indices, :]) @assert size(nullsp, 2) == 1 v = vec(nullsp) - push!(spanners, canonical!(v)) + push!(spanners, _canonical!(v)) end if indices[end] != n - 1 push!(indices, indices[end] + 1) @@ -272,7 +272,7 @@ function _compute_hyperplanes(C::Vector{Vector{T}}) where {T<:Number} end function _canonical!(v::Vector{Float64}) - normalize!(v) + LinearAlgebra.normalize!(v) for j in 1:length(v) if abs(v[j]) < 1e-8 v[j] = 0 diff --git a/src/pwlinear.jl b/src/pwlinear.jl index 1e7d8c1..882d390 100644 --- a/src/pwlinear.jl +++ b/src/pwlinear.jl @@ -106,6 +106,39 @@ function piecewiselinear( 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} + initPWL!(model) + counter = model.ext[:PWL].counter + counter += 1 + model.ext[:PWL].counter = counter + + if isempty(pwl.segments) + error( + "I don't know how to handle a piecewise linear function with no breakpoints.", + ) + end + + if output_vars === nothing + output_vars = tuple( + JuMP.@variable( + model, + [i in 1:F], + base_name = "y_$counter" + )..., + ) + end + + formulate_pwl!(model, input_vars, output_vars, pwl, method, direction) + return output_vars +end + # Specialization for univariate problems function piecewiselinear( model::JuMP.Model, diff --git a/src/types.jl b/src/types.jl index fdcced5..eb8828d 100644 --- a/src/types.jl +++ b/src/types.jl @@ -6,11 +6,30 @@ 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. + +# 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}} @@ -27,11 +46,27 @@ struct SegmentPointRep{D,F} <: Segment{D,F} end end +""" + AffineFunction{D} + +An affine function `f(x) = coeffs ⋅ x + offset` in `D` dimensions. +""" struct AffineFunction{D} coeffs::NTuple{D,Float64} offset::Float64 end +""" + SegmentHyperplaneRep{D,F} + +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}} @@ -47,6 +82,24 @@ 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} @@ -58,7 +111,49 @@ 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. + +# 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) 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 67c98a4..98c7c4c 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -225,3 +225,190 @@ append!(method_pattern, uj_methods) @test objective_value(model) ≈ 0.148753 rtol = 1e-3 @test objective_value(model) ≈ value(z) rtol = 1e-3 end + +# Test OptimalTriangleSelection with one method for coverage +@testset "OptimalTriangleSelection" begin + method = OptimalTriangleSelection(optimizer, ConvexCombination()) + 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 + +# 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,)]) + 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 + +# 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 + + # 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,), pwl; method = FakeMethod()) +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") \ No newline at end of file From 0b750fc520dd72ed6c20f3d6be9412becc67035e Mon Sep 17 00:00:00 2001 From: Valentin Kaisermayer Date: Tue, 31 Mar 2026 20:29:02 +0200 Subject: [PATCH 2/3] Add _next_pwl_id! and _pwl_name helpers --- src/PiecewiseLinearOpt.jl | 22 +++++++++++++++++++ src/methods/bivariate/common.jl | 7 +----- src/methods/bivariate/k1.jl | 4 +--- src/methods/bivariate/nine_stencil.jl | 4 +--- .../optimal_independent_branching.jl | 9 ++------ .../bivariate/optimal_triangle_selection.jl | 3 +-- src/methods/bivariate/six_stencil.jl | 6 ++--- src/methods/bivariate/union_jack.jl | 4 +--- .../multivariate/convex_combination.jl | 8 +++---- .../multivariate/disaggregated_logarithmic.jl | 6 ++--- src/methods/multivariate/multiple_choice.jl | 7 +++--- src/methods/univariate/incremental.jl | 5 ++--- .../univariate/logarithmic_embedding.jl | 3 +-- .../logarithmic_independent_branching.jl | 3 +-- src/methods/univariate/zig_zag_binary.jl | 3 +-- src/methods/univariate/zig_zag_integer.jl | 3 +-- src/methods/util.jl | 6 ++--- src/pwlinear.jl | 17 +++++--------- test/Project.toml | 7 ++++++ 19 files changed, 59 insertions(+), 68 deletions(-) create mode 100644 test/Project.toml diff --git a/src/PiecewiseLinearOpt.jl b/src/PiecewiseLinearOpt.jl index dc6f231..db4da52 100644 --- a/src/PiecewiseLinearOpt.jl +++ b/src/PiecewiseLinearOpt.jl @@ -26,6 +26,28 @@ function initPWL!(m::JuMP.Model) 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") diff --git a/src/methods/bivariate/common.jl b/src/methods/bivariate/common.jl index fd01de7..5fecf46 100644 --- a/src/methods/bivariate/common.jl +++ b/src/methods/bivariate/common.jl @@ -14,11 +14,6 @@ function formulate_pwl!( method::BivariateSOS2Method, direction::DIRECTION, ) where {F} - initPWL!(model) - counter = model.ext[:PWL].counter - counter += 1 - model.ext[:PWL].counter = counter - grid = _continuous_gridpoints_or_die(pwl) xs, ys = grid.input_vals, grid.output_vals @@ -37,7 +32,7 @@ function formulate_pwl!( [1:n_1, 1:n_2], lower_bound = 0, upper_bound = 1, - base_name = "λ_$counter" + base_name = _pwl_name(model, "λ") ) JuMP.@constraint(model, sum(λ) == 1) JuMP.@constraint( diff --git a/src/methods/bivariate/k1.jl b/src/methods/bivariate/k1.jl index 74b5051..a55e3be 100644 --- a/src/methods/bivariate/k1.jl +++ b/src/methods/bivariate/k1.jl @@ -31,11 +31,9 @@ function formulate_triangle_selection!( ) n_1, n_2 = size(λ) @assert size(triangle_direction) == (n_1 - 1, n_2 - 1) - counter = model.ext[:PWL].counter - # TODO: Certify triangulation is "uniform" - w = JuMP.@variable(model, [1:2], Bin, base_name = "w_$counter") + w = JuMP.@variable(model, [1:2], Bin, base_name = _pwl_name(model, "w")) JuMP.@constraints( model, begin diff --git a/src/methods/bivariate/nine_stencil.jl b/src/methods/bivariate/nine_stencil.jl index 417aac9..77aad46 100644 --- a/src/methods/bivariate/nine_stencil.jl +++ b/src/methods/bivariate/nine_stencil.jl @@ -32,9 +32,7 @@ function formulate_triangle_selection!( ) n_1, n_2 = size(λ) @assert size(triangle_direction) == (n_1 - 1, n_2 - 1) - counter = model.ext[:PWL].counter - - w = JuMP.@variable(model, [1:3, 1:3], Bin, base_name = "w_$counter") + 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 diff --git a/src/methods/bivariate/optimal_independent_branching.jl b/src/methods/bivariate/optimal_independent_branching.jl index 4468915..a2417fb 100644 --- a/src/methods/bivariate/optimal_independent_branching.jl +++ b/src/methods/bivariate/optimal_independent_branching.jl @@ -16,11 +16,6 @@ function formulate_pwl!( method::OptimalIndependentBranching, direction::DIRECTION, ) where {F} - initPWL!(model) - counter = model.ext[:PWL].counter - counter += 1 - model.ext[:PWL].counter = counter - grid = _continuous_gridpoints_or_die(pwl) xs, ys = grid.input_vals, grid.output_vals @@ -39,7 +34,7 @@ function formulate_pwl!( [1:n_1, 1:n_2], lower_bound = 0, upper_bound = 1, - base_name = "λ_$counter" + base_name = _pwl_name(model, "λ") ) JuMP.@constraint(model, sum(λ) == 1) JuMP.@constraint( @@ -165,7 +160,7 @@ function formulate_pwl!( @show t end end - z = JuMP.@variable(model, [1:t], Bin, base_name = "z_$counter") + z = JuMP.@variable(model, [1:t], Bin, base_name = _pwl_name(model, "z")) for k in 1:t JuMP.@constraints( diff --git a/src/methods/bivariate/optimal_triangle_selection.jl b/src/methods/bivariate/optimal_triangle_selection.jl index e1ee0b6..5e89663 100644 --- a/src/methods/bivariate/optimal_triangle_selection.jl +++ b/src/methods/bivariate/optimal_triangle_selection.jl @@ -21,7 +21,6 @@ function formulate_triangle_selection!( method::OptimalTriangleSelection, _::GridTriangulation, ) - counter = model.ext[:PWL].counter n_1, n_2 = size(λ) J = Set((i, j) for i in 1:n_1, j in 1:n_2) x_val = Matrix{Float64}(undef, n_1, n_2) @@ -102,7 +101,7 @@ function formulate_triangle_selection!( t += 1 end end - z = JuMP.@variable(model, [1:t], Bin, base_name = "z_$counter") + z = JuMP.@variable(model, [1:t], Bin, base_name = _pwl_name(model, "z")) for k in 1:t JuMP.@constraints( diff --git a/src/methods/bivariate/six_stencil.jl b/src/methods/bivariate/six_stencil.jl index dcd0b98..9e88d69 100644 --- a/src/methods/bivariate/six_stencil.jl +++ b/src/methods/bivariate/six_stencil.jl @@ -32,10 +32,8 @@ function formulate_triangle_selection!( ) n_1, n_2 = size(λ) @assert size(triangle_direction) == (n_1 - 1, n_2 - 1) - counter = model.ext[:PWL].counter - # diagonal lines running from SW to NE. Grouped with an offset of 3. - w_sw_ne = JuMP.@variable(model, [0:2], Bin, base_name = "w_sw_ne_$counter") + 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}}() @@ -89,7 +87,7 @@ function formulate_triangle_selection!( ) end - w_se_nw = JuMP.@variable(model, [0:2], Bin, base_name = "w_se_nw_$counter") + 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}}() diff --git a/src/methods/bivariate/union_jack.jl b/src/methods/bivariate/union_jack.jl index fdc61ee..07be517 100644 --- a/src/methods/bivariate/union_jack.jl +++ b/src/methods/bivariate/union_jack.jl @@ -32,8 +32,6 @@ function formulate_triangle_selection!( ) n_1, n_2 = size(λ) @assert size(triangle_direction) == (n_1 - 1, n_2 - 1) - counter = model.ext[:PWL].counter - # TODO: Certify triangulation is Union Jack # j_start = number of triangle segments incident to (1, 1), the lower-left @@ -41,7 +39,7 @@ function formulate_triangle_selection!( j_start = triangle_direction[1, 1] ? 1 : 2 # diagonal lines running from SW to NE. Grouped with an offset of 3. - z = JuMP.@variable(model, binary = true, base_name = "w_$counter") + z = JuMP.@variable(model, binary = true, base_name = _pwl_name(model, "w")) JuMP.@constraints( model, diff --git a/src/methods/multivariate/convex_combination.jl b/src/methods/multivariate/convex_combination.jl index a50b023..741585c 100644 --- a/src/methods/multivariate/convex_combination.jl +++ b/src/methods/multivariate/convex_combination.jl @@ -22,9 +22,8 @@ function formulate_pwl!( direction::DIRECTION, ) where {D,F} # TODO: assert PWL function is continuous - counter = model.ext[:PWL].counter segments = pwl.segments - z = JuMP.@variable(model, [segments], Bin, base_name = "z_$counter") + 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}}() @@ -49,7 +48,7 @@ function formulate_pwl!( [all_input_vals], lower_bound = 0, upper_bound = 1, - base_name = "λ_$counter" + base_name = _pwl_name(model, "λ") ) JuMP.@constraint(model, sum(λ) == 1) for j in 1:D @@ -81,9 +80,8 @@ function formulate_sos2!( λ::Vector{T}, method::ConvexCombination, ) where {T<:VarOrAff} - counter = model.ext[:PWL].counter n = length(λ) - z = JuMP.@variable(model, [1:n-1], Bin, base_name = "z_$counter") + 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 diff --git a/src/methods/multivariate/disaggregated_logarithmic.jl b/src/methods/multivariate/disaggregated_logarithmic.jl index e03d6f0..bdb63f9 100644 --- a/src/methods/multivariate/disaggregated_logarithmic.jl +++ b/src/methods/multivariate/disaggregated_logarithmic.jl @@ -22,8 +22,6 @@ function formulate_pwl!( method::DisaggregatedLogarithmic, direction::DIRECTION, ) where {D,F} - counter = model.ext[:PWL].counter - segments = pwl.segments num_bps = Dict(seg => length(seg.input_vals) for seg in segments) @@ -32,7 +30,7 @@ function formulate_pwl!( [seg in segments, i in 1:num_bps[seg]], lower_bound = 0, upper_bound = 1, - base_name = "γ_$counter" + base_name = _pwl_name(model, "γ") ) JuMP.@constraint(model, sum(γ) == 1) @@ -59,7 +57,7 @@ function formulate_pwl!( end _H = _reflected_gray_codes(r) H = Dict(segments[i] => _H[i] for i in 1:length(segments)) - z = JuMP.@variable(model, [1:r], Bin, base_name = "z_$counter") + z = JuMP.@variable(model, [1:r], Bin, base_name = _pwl_name(model, "z")) for j in 1:r JuMP.@constraint( model, diff --git a/src/methods/multivariate/multiple_choice.jl b/src/methods/multivariate/multiple_choice.jl index 48222d7..b5bfdc5 100644 --- a/src/methods/multivariate/multiple_choice.jl +++ b/src/methods/multivariate/multiple_choice.jl @@ -23,12 +23,11 @@ function formulate_pwl!( method::MultipleChoice, direction::DIRECTION, ) where {D,F} - counter = model.ext[:PWL].counter segments = pwl.segments S = 1:length(segments) - x_hat = JuMP.@variable(model, [S, 1:D], base_name = "x_hat_$counter") - y_hat = JuMP.@variable(model, [S, 1:F], base_name = "y_hat_$counter") - z = JuMP.@variable(model, [S], Bin, base_name = "z_$counter") + 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]) diff --git a/src/methods/univariate/incremental.jl b/src/methods/univariate/incremental.jl index 52d5719..3961a75 100644 --- a/src/methods/univariate/incremental.jl +++ b/src/methods/univariate/incremental.jl @@ -26,7 +26,6 @@ function formulate_pwl!( grid = _continuous_gridpoints_or_die(pwl) xs, ys = grid.input_vals, grid.output_vals - counter = model.ext[:PWL].counter n = length(pwl.segments) + 1 @assert length(xs) == length(ys) == n @@ -35,9 +34,9 @@ function formulate_pwl!( [1:n], lower_bound = 0, upper_bound = 1, - base_name = "δ_$counter" + base_name = _pwl_name(model, "δ") ) - z = JuMP.@variable(model, [1:n-1], Bin, base_name = "z_$counter") + z = JuMP.@variable(model, [1:n-1], Bin, base_name = _pwl_name(model, "z")) JuMP.@constraint( model, input_vars[1] == diff --git a/src/methods/univariate/logarithmic_embedding.jl b/src/methods/univariate/logarithmic_embedding.jl index f82fae9..bf3f5a7 100644 --- a/src/methods/univariate/logarithmic_embedding.jl +++ b/src/methods/univariate/logarithmic_embedding.jl @@ -20,7 +20,6 @@ function formulate_sos2!( λ::Vector{T}, method::LogarithmicEmbedding, ) where {T<:VarOrAff} - counter = model.ext[:PWL].counter n = length(λ) d = n - 1 if 0 <= d <= 1 @@ -30,7 +29,7 @@ function formulate_sos2!( if k == 0 return nothing end - y = JuMP.@variable(model, [1:k], Bin, base_name = "y_$counter") + y = JuMP.@variable(model, [1:k], Bin, base_name = _pwl_name(model, "y")) _sos2_encoding_constraints!( model, λ, diff --git a/src/methods/univariate/logarithmic_independent_branching.jl b/src/methods/univariate/logarithmic_independent_branching.jl index 95195ac..e4ca907 100644 --- a/src/methods/univariate/logarithmic_independent_branching.jl +++ b/src/methods/univariate/logarithmic_independent_branching.jl @@ -18,7 +18,6 @@ function formulate_sos2!( λ::Vector{T}, method::LogarithmicIndependentBranching, ) where {T<:VarOrAff} - counter = model.ext[:PWL].counter n = length(λ) d = n - 1 if 0 <= d <= 1 @@ -28,7 +27,7 @@ function formulate_sos2!( if k == 0 return nothing end - z = JuMP.@variable(model, [1:k], Bin, base_name = "y_$counter") + 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] diff --git a/src/methods/univariate/zig_zag_binary.jl b/src/methods/univariate/zig_zag_binary.jl index 753bcf3..a5b3fe0 100644 --- a/src/methods/univariate/zig_zag_binary.jl +++ b/src/methods/univariate/zig_zag_binary.jl @@ -18,7 +18,6 @@ function formulate_sos2!( λ::Vector{T}, method::ZigZagBinary, ) where {T<:VarOrAff} - counter = model.ext[:PWL].counter n = length(λ) d = n - 1 if 0 <= d <= 1 @@ -28,7 +27,7 @@ function formulate_sos2!( if k == 0 return nothing end - y = JuMP.@variable(model, [1:k], Bin, base_name = "y_$counter") + y = JuMP.@variable(model, [1:k], Bin, base_name = _pwl_name(model, "y")) _sos2_encoding_constraints!( model, λ, diff --git a/src/methods/univariate/zig_zag_integer.jl b/src/methods/univariate/zig_zag_integer.jl index 9b37c3f..aef91a0 100644 --- a/src/methods/univariate/zig_zag_integer.jl +++ b/src/methods/univariate/zig_zag_integer.jl @@ -18,7 +18,6 @@ function formulate_sos2!( λ::Vector{T}, method::ZigZagInteger, ) where {T<:VarOrAff} - counter = model.ext[:PWL].counter n = length(λ) d = n - 1 if 0 <= d <= 1 @@ -34,7 +33,7 @@ function formulate_sos2!( Int, lower_bound = 0, upper_bound = 2^(k - i), - base_name = "y_$counter" + base_name = _pwl_name(model, "y") ) _sos2_encoding_constraints!( model, diff --git a/src/methods/util.jl b/src/methods/util.jl index bd5dd56..b6053ae 100644 --- a/src/methods/util.jl +++ b/src/methods/util.jl @@ -130,12 +130,10 @@ function _create_convex_multiplier_vars( output_vars::NTuple{F,JuMP.VariableRef}, direction::DIRECTION, ) where {D,F} - counter = model.ext[:PWL].counter - - # TODO: Name λ variables λ = similar(Array{JuMP.VariableRef}, axes(grid.input_vals)) + λ_name = _pwl_name(model, "λ") for I in eachindex(λ) - λ[I] = JuMP.@variable(model, lower_bound = 0, upper_bound = 1) + λ[I] = JuMP.@variable(model, lower_bound = 0, upper_bound = 1, base_name = λ_name) end JuMP.@constraint(model, sum(λ) == 1) diff --git a/src/pwlinear.jl b/src/pwlinear.jl index 882d390..63fb4f7 100644 --- a/src/pwlinear.jl +++ b/src/pwlinear.jl @@ -16,16 +16,15 @@ function formulate_pwl!( ) 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) @@ -74,10 +73,7 @@ function piecewiselinear( direction::DIRECTION = Graph, output_vars::Union{Nothing,NTuple{F,VarOrAff}} = nothing, ) where {D,F} - initPWL!(model) - counter = model.ext[:PWL].counter - counter += 1 - model.ext[:PWL].counter = counter + _next_pwl_id!(model) if isempty(pwl.segments) error( @@ -97,7 +93,7 @@ function piecewiselinear( [i in 1:F], lower_bound = output_lb[i], upper_bound = output_ub[i], - base_name = "y_$counter" + base_name = _pwl_name(model, "y") )..., ) end @@ -114,10 +110,7 @@ function piecewiselinear( direction::DIRECTION = Graph, output_vars::Union{Nothing,NTuple{F,VarOrAff}} = nothing, ) where {D,F} - initPWL!(model) - counter = model.ext[:PWL].counter - counter += 1 - model.ext[:PWL].counter = counter + _next_pwl_id!(model) if isempty(pwl.segments) error( @@ -130,7 +123,7 @@ function piecewiselinear( JuMP.@variable( model, [i in 1:F], - base_name = "y_$counter" + base_name = _pwl_name(model, "y") )..., ) end diff --git a/test/Project.toml b/test/Project.toml new file mode 100644 index 0000000..81ca24d --- /dev/null +++ b/test/Project.toml @@ -0,0 +1,7 @@ +[deps] +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" From 4718d87d54969630f55954c4a8749d49c25c4a4a Mon Sep 17 00:00:00 2001 From: Valentin Kaisermayer Date: Wed, 1 Apr 2026 12:56:25 +0200 Subject: [PATCH 3/3] formatting and speed up tests --- Project.toml | 7 +- src/methods/bivariate/nine_stencil.jl | 7 +- .../bivariate/optimal_triangle_selection.jl | 2 +- src/methods/bivariate/six_stencil.jl | 14 +- .../multivariate/convex_combination.jl | 7 +- src/methods/multivariate/multiple_choice.jl | 6 +- src/methods/util.jl | 7 +- test/Project.toml | 1 + test/runtests.jl | 822 ++++++++++-------- 9 files changed, 491 insertions(+), 382 deletions(-) diff --git a/Project.toml b/Project.toml index 07c3da3..4d21aff 100644 --- a/Project.toml +++ b/Project.toml @@ -12,10 +12,15 @@ Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" HiGHS = "1" JuMP = "1" julia = "1.10" +LinearAlgebra = "1" +Random = "1" +Aqua = "0.8" +Test = "1" [extras] +Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["Test", "HiGHS"] +test = ["Aqua", "Test", "HiGHS"] diff --git a/src/methods/bivariate/nine_stencil.jl b/src/methods/bivariate/nine_stencil.jl index 77aad46..ab3566a 100644 --- a/src/methods/bivariate/nine_stencil.jl +++ b/src/methods/bivariate/nine_stencil.jl @@ -32,7 +32,12 @@ function formulate_triangle_selection!( ) 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")) + 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 diff --git a/src/methods/bivariate/optimal_triangle_selection.jl b/src/methods/bivariate/optimal_triangle_selection.jl index 5e89663..e34e55e 100644 --- a/src/methods/bivariate/optimal_triangle_selection.jl +++ b/src/methods/bivariate/optimal_triangle_selection.jl @@ -28,7 +28,7 @@ function formulate_triangle_selection!( t = 1 while true - @show method.sub_solver + # @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) diff --git a/src/methods/bivariate/six_stencil.jl b/src/methods/bivariate/six_stencil.jl index 9e88d69..9c2564a 100644 --- a/src/methods/bivariate/six_stencil.jl +++ b/src/methods/bivariate/six_stencil.jl @@ -33,7 +33,12 @@ function formulate_triangle_selection!( 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")) + 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}}() @@ -87,7 +92,12 @@ function formulate_triangle_selection!( ) end - w_se_nw = JuMP.@variable(model, [0:2], Bin, base_name = _pwl_name(model, "w_se_nw")) + 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}}() diff --git a/src/methods/multivariate/convex_combination.jl b/src/methods/multivariate/convex_combination.jl index 741585c..b759ac5 100644 --- a/src/methods/multivariate/convex_combination.jl +++ b/src/methods/multivariate/convex_combination.jl @@ -23,7 +23,12 @@ function formulate_pwl!( ) where {D,F} # TODO: assert PWL function is continuous segments = pwl.segments - z = JuMP.@variable(model, [segments], Bin, base_name = _pwl_name(model, "z")) + 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}}() diff --git a/src/methods/multivariate/multiple_choice.jl b/src/methods/multivariate/multiple_choice.jl index b5bfdc5..8090da0 100644 --- a/src/methods/multivariate/multiple_choice.jl +++ b/src/methods/multivariate/multiple_choice.jl @@ -25,8 +25,10 @@ function formulate_pwl!( ) 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")) + 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 diff --git a/src/methods/util.jl b/src/methods/util.jl index b6053ae..cbf90ef 100644 --- a/src/methods/util.jl +++ b/src/methods/util.jl @@ -133,7 +133,12 @@ function _create_convex_multiplier_vars( λ = 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) + λ[I] = JuMP.@variable( + model, + lower_bound = 0, + upper_bound = 1, + base_name = λ_name + ) end JuMP.@constraint(model, sum(λ) == 1) diff --git a/test/Project.toml b/test/Project.toml index 81ca24d..8dae788 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -1,4 +1,5 @@ [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" diff --git a/test/runtests.jl b/test/runtests.jl index 98c7c4c..ef3d86f 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -8,7 +8,6 @@ using HiGHS using JuMP using LinearAlgebra using Test - const PLO = PiecewiseLinearOpt optimizer = optimizer_with_attributes(HiGHS.Optimizer, MOI.Silent() => true) @@ -23,93 +22,6 @@ const methods_1D = [ ZigZagBinary(), ZigZagInteger(), ] - -@testset "Simple univariate" for method in methods_1D - model = Model(optimizer) - @variable(model, x) - - s1 = PLO.SegmentPointRep{1,1}([(1.0,), (2.0,)], [(2.5,), (3.5,)]) - s2 = PLO.SegmentPointRep{1,1}([(2.0,), (3.0,)], [(3.5,), (1.0,)]) - pwl = PLO.PWLFunction([s1, s2], PLO.Intervals()) - - y = piecewiselinear(model, (x,), pwl; method = method) - @objective(model, Min, y[1]) - - optimize!(model) - @test termination_status(model) == MOI.OPTIMAL - @test value(x) ≈ 3.0 rtol = 1e-4 - @test value(y[1]) ≈ 1.0 rtol = 1e-4 -end - -@testset "Univariate pwlinear" begin - d = 0:0.01:1 - f = (xi -> xi^2) - fd = [f(xi) for xi in d] - pwl = UnivariatePWLFunction(d, f) - - model = Model(optimizer) - @variable(model, x) - y1 = piecewiselinear(model, (x,), pwl) - @constraint(model, x ≤ 0.75) - @objective(model, Max, y1[1]) - optimize!(model) - @test termination_status(model) == MOI.OPTIMAL - @test value(y1[1]) ≈ f(value(x)) rtol = 1e-4 - @test objective_value(model) ≈ 0.5625 rtol = 1e-4 - - model = Model(optimizer) - @variable(model, x) - y2 = piecewiselinear(model, x, d, f) - @constraint(model, x ≤ 0.75) - @objective(model, Max, y2) - optimize!(model) - @test termination_status(model) == MOI.OPTIMAL - @test value(y2) ≈ f(value(x)) rtol = 1e-4 - @test objective_value(model) ≈ 0.5625 rtol = 1e-4 - - model = Model(optimizer) - @variable(model, x) - y3 = piecewiselinear(model, x, d, fd) - @constraint(model, x ≤ 0.75) - @objective(model, Max, y3) - optimize!(model) - @test termination_status(model) == MOI.OPTIMAL - @test value(y3) ≈ f(value(x)) rtol = 1e-4 - @test objective_value(model) ≈ 0.5625 rtol = 1e-4 -end - -@testset "Bivariate pwlinear" begin - d = 0:0.05:1 - f = (xi, yi) -> xi^2 + yi^2 - pwl = BivariatePWLFunction(d, d, f) - - model = Model(optimizer) - @variable(model, x) - @variable(model, y) - z1 = piecewiselinear(model, (x, y), pwl) - @constraint(model, x ≤ 0.75) - @constraint(model, y ≤ 0.75) - @objective(model, Max, z1[1]) - optimize!(model) - @test termination_status(model) == MOI.OPTIMAL - @test value(x) ≈ 0.75 rtol = 1e-4 - @test value(y) ≈ 0.75 rtol = 1e-4 - @test value(z1[1]) ≈ 1.125 rtol = 1e-4 - - model = Model(optimizer) - @variable(model, x) - @variable(model, y) - z2 = piecewiselinear(model, x, y, d, d, f) - @constraint(model, x ≤ 0.75) - @constraint(model, y ≤ 0.75) - @objective(model, Max, z2) - optimize!(model) - @test termination_status(model) == MOI.OPTIMAL - @test value(x) ≈ 0.75 rtol = 1e-4 - @test value(y) ≈ 0.75 rtol = 1e-4 - @test value(z2) ≈ 1.125 rtol = 1e-4 -end - const sos2_methods = [ ConvexCombination(), LogarithmicEmbedding(), @@ -123,292 +35,456 @@ const methods_2D_gen = [ DisaggregatedLogarithmic(), #OptimalIndependentBranching(optimizer), [NineStencil(sos2_method) for sos2_method in methods_1D]..., - [ - OptimalTriangleSelection(optimizer, 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 "Simple bivariate" for method in methods_2D_gen - model = Model(optimizer) - @variable(model, x[1:2]) - s1 = PLO.SegmentPointRep{2,1}( - [(0.0, 0.0), (0.0, 1.0), (1.0, 1.0)], - [(0.0,), (1.0,), (2.0,)], - ) - s2 = PLO.SegmentPointRep{2,1}( - [(0.0, 0.0), (1.0, 0.0), (1.0, 1.0)], - [(0.0,), (3.0,), (2.0,)], - ) - pwl = PLO.PWLFunction{2,1,PLO.SegmentPointRep{2,1}}( - [s1, s2], - PLO.UnstructuredTriangulation(), - ) - - y = piecewiselinear(model, (x[1], x[2]), pwl; method = method) - @objective(model, Min, y[1]) - - optimize!(model) - - @test termination_status(model) == MOI.OPTIMAL - @test value(x[1]) ≈ 0.0 rtol = 1e-4 - @test value(x[2]) ≈ 0.0 rtol = 1e-4 - @test value(y[1]) ≈ 0.0 rtol = 1e-4 -end -@testset "1D: $method" for method in methods_1D - model = Model(optimizer) - @variable(model, x) - d = 7 - xs = collect(range(1; stop = 2π, length = (d + 1))) - zs = sin.(xs) - y = piecewiselinear(model, x, xs, zs; method = method) - @objective(model, Max, y) - optimize!(model) - @test termination_status(model) == MOI.OPTIMAL - @test value(x) ≈ 1.75474 rtol = 1e-4 - @test value(y) ≈ 0.98313 rtol = 1e-4 - @test objective_value(model) ≈ 0.98313 rtol = 1e-4 - @test objective_value(model) ≈ value(y) rtol = 1e-4 - @constraint(model, x ≤ 1.5y) - optimize!(model) - @test termination_status(model) == MOI.OPTIMAL - @test value(x) ≈ 1.36495 rtol = 1e-4 - @test value(y) ≈ 0.90997 rtol = 1e-4 - @test objective_value(model) ≈ 0.90997 rtol = 1e-4 - @test objective_value(model) ≈ value(y) rtol = 1e-4 +@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,)], + ) + 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 + + # 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 + + # 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,), + pwl; + method = FakeMethod(), + ) + 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 - -patterns = [:Upper, :Lower, :BestFit, :K1, :UnionJack, :Random] -method_pattern = vec(collect(Iterators.product(methods_2D_gen, patterns))) -k1_methods = [ - (method, :K1) for method in [K1(sos2_method) for sos2_method in methods_1D] -] -uj_methods = [ - (method, :UnionJack) for - method in [UnionJack(sos2_method) for sos2_method in methods_1D] -] -append!(method_pattern, k1_methods) -append!(method_pattern, uj_methods) - -@testset "2D: $method, $pattern" for (method, pattern) in method_pattern - model = Model(optimizer) - @variable(model, x[1:2]) - d = range(0; stop = 1, length = 8) - f = (x1, x2) -> 2 * (x1 - 1 / 3)^2 + 3 * (x2 - 4 / 7)^4 - z = piecewiselinear( - model, - x[1], - x[2], - d, - d, - f; - method = method, - pattern = pattern, - ) - @objective(model, Min, z) - optimize!(model) - @test termination_status(model) == MOI.OPTIMAL - @test value(x[1]) ≈ 0.285714 rtol = 1e-4 - @test value(x[2]) ≈ 0.571429 rtol = 1e-4 - @test value(z) ≈ 0.004535 rtol = 1e-3 - @test objective_value(model) ≈ 0.004535 rtol = 1e-3 - @test objective_value(model) ≈ value(z) rtol = 1e-3 - - @constraint(model, x[1] ≥ 0.6) - optimize!(model) - - @test termination_status(model) == MOI.OPTIMAL - @test value(x[1]) ≈ 0.6 rtol = 1e-4 - @test value(x[2]) ≈ 0.571428 rtol = 1e-4 - @test value(z) ≈ 0.148753 rtol = 1e-4 - @test objective_value(model) ≈ 0.148753 rtol = 1e-3 - @test objective_value(model) ≈ value(z) rtol = 1e-3 -end - -# Test OptimalTriangleSelection with one method for coverage -@testset "OptimalTriangleSelection" begin - method = OptimalTriangleSelection(optimizer, ConvexCombination()) - 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 - -# 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,)]) - 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 - -# 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 - - # 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,), pwl; method = FakeMethod()) -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") \ No newline at end of file