diff --git a/README.md b/README.md index de8d989..fcb6ce2 100644 --- a/README.md +++ b/README.md @@ -25,7 +25,32 @@ using JuMP, PiecewiseAffineApprox, HiGHS m = Model(HiGHS.Optimizer) @variable(m, x) # Create a piecewise linear approximation to x^2 on the interval [-1, 1] -pwa = approx(x -> x[1]^2, [(-1, 1)], Convex(), MILP(optimizer = HiGHS.Optimizer, planes=5)) +pwa = approx(x -> x^2, [(-1, 1)], Convex(), MILP(optimizer = HiGHS.Optimizer, planes=5)) + +cluster = Cluster(optimizer = HiGHS.Optimizer, planes=5) + +# Alternative formulation with defined evaluation points +pwa_alt = approx(x -> x^2, -1:0.1:1, Convex(), cluster) +# 2 dimensional variant +pwa_2d = approx((x, y) -> x^2 + y^2, -1:0.1:1, -1:0.1:1, Convex(), cluster) + +# Another variant with explicit points +pwa_pts = approx(LinRange(0, 1, 10), [1 + i^2 for i in 1:10], Convex(), cluster) + +# Input as a matrix +fv = [1 2 3 4 5 6 + 1 4 9 16 25 36] +pwa_mat = approx(fv, Convex(), cluster) +fv_2d = [1 1 1 1 2 2 2 2 3 3 3 3 4 4 4 4 + 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4 + 2 5 10 17 5 8 13 20 10 13 18 25 17 20 25 32] +pwa_mat = approx(fv_2d, Convex(), cluster) + +# Input from a csv-file with row-based observed values +using DataFrames +data = DataFrame(CSV.File("test/observations.csv")) +pwa_data = approx(Matrix(data), Convex(), cluster) + # Add the pwa function to the model z = pwaffine(m, x, pwa) # Minimize @@ -45,9 +70,9 @@ The following demonstrates how this can be achieved: using PiecewiseAffineApprox, CairoMakie, HiGHS x = LinRange(0, 1, 20) -f(x) = first(x)^2 +f(x) = x^2 pwa = approx(f, [(0, 1)], Convex(), MILP(optimizer = HiGHS.Optimizer, planes = 3)) -p = plot(x, f.(x), pwa) +p = Makie.plot(x, f.(x), pwa) save("approx.svg", p) ``` diff --git a/docs/approx.svg b/docs/approx.svg index b0b2415..4f0c2d9 100644 --- a/docs/approx.svg +++ b/docs/approx.svg @@ -1,168 +1,155 @@ - + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - + + - - - - - - - - - - - + + + + + + + + - - - - - - - - - + + + + + + + + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - - + + + + + + + + + + + + + + + + + + + + + + + - - - - - - - - + + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + diff --git a/docs/make.jl b/docs/make.jl index 95ba769..53f612c 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -36,7 +36,6 @@ pages = [ ] # TODO: Enable modules -# TODO: Enable doctests Documenter.makedocs( sitename = "PiecewiseAffineApprox", format = Documenter.HTML(; diff --git a/docs/src/howto/approximate_points.md b/docs/src/howto/approximate_points.md index 47b8766..7dd4f12 100644 --- a/docs/src/howto/approximate_points.md +++ b/docs/src/howto/approximate_points.md @@ -24,6 +24,27 @@ length(pwa.planes) 5 ``` +You may find it more convenient to use Matrix input for the data points rather than wrapping the data in the FunctionEvaluations, consider the alternative version of the example above: + +```jldoctest +using PiecewiseAffineApprox, JuMP, HiGHS +optimizer = optimizer_with_attributes(HiGHS.Optimizer, MOI.Silent()=>true) + +x = collect(range(-1, 1; length = 10)) +z = x .^ 2 +m = hcat(x, z)' # Collect data in a single Matrix + +pwa = approx( + m, + Convex(), + MILP(;optimizer, metric = :l1, planes = 5), +) + +length(pwa.planes) +# output +5 +``` + # Bi-variate functions (3D) A dataset can also be used as input for 3D piecewise affine approximations. The following code creates a uniformly sampled domain X within (-1,1) and calculates the corresponding function values for a concave function z_concave. @@ -56,3 +77,40 @@ length(pwa.planes) # output 17 ``` + +Similarly, the data can be passed as a Matrix if that is more convenient: + +```jldoctest +using PiecewiseAffineApprox, JuMP, HiGHS +optimizer = optimizer_with_attributes(HiGHS.Optimizer, MOI.Silent()=>true) + +xg = [i for i ∈ -1:0.5:1] +yg = [j for j ∈ -1:0.5:1] + +X = [repeat(xg, inner = [size(yg, 1)]) repeat(yg, outer = [size(xg, 1)])] + +z = X[:, 1] .^ 2 + X[:, 2] .^ 2 + +m = hcat(X, z)' + +np = 17 + +pwa = approx( + m, + Convex(), + MILP( + ;optimizer, + planes = np, + strict = :outer, + metric = :l1, + ), +) +length(pwa.planes) + +# output +17 +``` + + + + diff --git a/docs/src/reference/api.md b/docs/src/reference/api.md index e409f24..3a1d503 100644 --- a/docs/src/reference/api.md +++ b/docs/src/reference/api.md @@ -5,19 +5,18 @@ CurrentModule = PiecewiseAffineApprox ## Algorithms ```@docs -#FullOrderFitting -#Heuristic -#Interpol -#Optimized -#ProgressiveFitting +Cluster +Interpol +MILP +Progressive +FullOrder ``` ## Other types ```@docs - -#Convex # Missing docstring -#Concave # Missing docstring +# Convex # Not found when building docs +# Concave # Not found when building docs Plane FunctionEvaluations PWAFunc diff --git a/src/convexapprox.jl b/src/convexapprox.jl index f659989..8b3c541 100644 --- a/src/convexapprox.jl +++ b/src/convexapprox.jl @@ -5,13 +5,47 @@ defaulttimelimit() = 60 """ approx(input::FunctionEvaluations{D}, c::Curvature, a::Algorithm) -Return PWAFunc{Convex,D} or PWAFunc{Concave,D} depending on `c`, approximating the `input` points in `D` dimensions +Return PWAFunc{Convex,D} or PWAFunc{Concave,D} depending on curvature `c`, approximating the `input` points in `D` dimensions + +## Arguments +- input +- c::Curvature `Convex` or `Concave` +- a::Algorithm `MILP`, `Cluster`, `Interpol`, `Progressive` or `FullOrder` """ function approx(input, c::Concave, a::Algorithm) cv = approx(FunctionEvaluations(input.points, -input.values), Convex(), a;) return PWAFunc{Concave,dims(cv)}(cv.planes) end + +""" + approx(fv::AbstractMatrix, c::Curvature, a::Algorithm;) + +Return approximation calculated by point estimates `fv` in Matrix. + +## Arguments +- fv::AbstractMatrix dimensions 1:m-1 are coordinats, dimension m is function values +- c::Curvature `Convex` or `Concave` +- a::Algorithm `MILP`, `Cluster`, `Interpol`, `Progressive` or `FullOrder` +""" +function approx(fv::AbstractMatrix, c::Curvature, a::Algorithm;) + m, n = size(fv) + + # Transpose if necessary, assuming that there are more points than values + if m > n + @warn("Transposing input matrix assuming row-based format") + fv = fv' + m, n = size(fv) + end + + xvals = fv[1:m-1, :] + fvals = fv[m, :] + pts = [Tuple(xi for xi in x) for x in eachcol(xvals)] + return approx(FunctionEvaluations(pts, fvals), c, a) +end + + + dims(pwa::PWAFunc{C,D}) where {C,D} = D # Using dispatch for specializing on dimensions. If performance were a concern, @@ -74,7 +108,14 @@ function approx(input::FunctionEvaluations, c::Convex, a::FullOrder;) return _full_order_pwa(input, a.optimizer, a.metric) end -# Optimal convex approximation using mixed integer optimization +""" + approx()input::FunctionEvaluations{D}, + c::Convex, + options::MILP, +) where {D} + +Approximate using MILP and data points sampled from a Convex function. +""" function approx( input::FunctionEvaluations{D}, c::Convex, @@ -180,7 +221,7 @@ function _sample_uniform(f::Function, bbox::Vector{<:Tuple}, nsamples) ) x = vec(collect(it)) end - y = [f(xx) for xx ∈ x] + y = [f(xx...) for xx ∈ x] return FunctionEvaluations(x, y) end @@ -198,6 +239,24 @@ function approx(f::Function, bbox::Vector{<:Tuple}, c::Curvature, a::Algorithm;) return approx(_sample_uniform(f, bbox, samples), c, a) end +function approx(f::Function, xvals, c::Curvature, a::Algorithm;) + pts = [(x,) for x in xvals] + yvals = [f(x) for x ∈ xvals] + return approx(FunctionEvaluations(pts, yvals), c, a) +end + +function approx(f::Function, xvals, yvals, c::Curvature, a::Algorithm;) + pts = [(x, y) for (x, y) in zip(xvals, yvals)] + yvals = [f(x, y) for (x, y) ∈ pts] + return approx(FunctionEvaluations(pts, yvals), c, a) +end + +function approx(xvals, fvals, c::Curvature, a::Algorithm;) + pts = [(x,) for x in xvals] + return approx(FunctionEvaluations(pts, collect(fvals)), c, a) +end + + # Utility function to find the value of a hyperplane at the point x function _plane_f(x, normal, d) return abs(normal[3]) > 1e-4 ? diff --git a/src/types.jl b/src/types.jl index 6b72333..6e89c0f 100644 --- a/src/types.jl +++ b/src/types.jl @@ -1,6 +1,20 @@ -# Input types +""" + Curvature + +Type parameter for PWAFunc to indicate curvature +""" abstract type Curvature end +""" + Concave + +Type parameter for PWAFunc when curvature is concave. +""" struct Concave <: Curvature end +""" + Convex + +Type parameter for PWAFunc when curvature is convex. +""" struct Convex <: Curvature end abstract type Algorithm end @@ -12,6 +26,15 @@ Note that this algoritm computes multiple approximations and selects the best. If julia is started with multiple threads, these are computed in parallel. Consider how many threads will be beneficial, particularly when using a commercial solver where the license may restrict the number of simultanous solver instances. + + ## Arguments +- planes::Int Number of planes to use in approximation +- metric::Symbol Allowed values are :l1, :l2, and :max +- trials::Int Number of randomized starts +- itlim::Int Number of iterations to try improving approximation +- strict::Symbol Allowed values are :none, :outer, and :inner +- maxtime::Int Time limit in seconds +- optimizer::T Optimizer to use to calculate approximation """ @kwdef struct Cluster{T} <: Algorithm planes::Int = defaultplanes() @@ -26,6 +49,11 @@ end Interpol Compute affine approximation by method proposed by Flatberg. Only available for 1D. + +## Arguments +- planes::Int Number of planes to use in approximation +- metric::Symbol Allowed values are :l1, :l2, and :max +- optimizer::T Optimizer to use to calculate approximation """ @kwdef struct Interpol{T} <: Algorithm planes::Int = defaultplanes() @@ -39,6 +67,13 @@ Compute affine approximation using a variation of the method proposed by Toriell Note that the resulting approximation is sensitive to the selection of the Big-M value used when solving the optimization problem. + +## Arguments +- planes::Int Number of planes to use in approximation +- metric::Symbol Allowed values are :l1, :l2, and :max +- strict::Symbol Allowed values are :none, :outer, and :inner +- maxtime::Int Time limit in seconds +- optimizer::T Optimizer to use to calculate approximation """ @kwdef struct MILP{T} <: Algorithm planes::Int = defaultplanes() @@ -53,6 +88,11 @@ end Compute affine approximation based on a variation of the method of Kazda and Li (2024) specialized to convex approximations. + +## Arguments +- tolerance::Float64 TODO: Check if relative gap +- metric::Symbol Allowed values are :l1, :l2, and :max +- optimizer::T Optimizer to use to calculate approximation """ @kwdef struct Progressive{T} <: Algorithm tolerance::Float64 = 0.01 @@ -65,6 +105,11 @@ end Compute affine approximation based on a variation of the method of Kazda and Li (2024) specialized to convex approximations with full order approximation (no reduction of number of planes). + +## Arguments +- metric::Symbol Allowed values are :l1, :l2, and :max +- optimizer::T Optimizer to use to calculate approximation + """ @kwdef struct FullOrder{T} <: Algorithm metric::Symbol = defaultmetric() diff --git a/test/observations.csv b/test/observations.csv new file mode 100644 index 0000000..28db574 --- /dev/null +++ b/test/observations.csv @@ -0,0 +1,122 @@ +x,y,value +0.0, 0.0, 0.054638237027137215 +0.0, 0.1, 0.067630837595478 +0.0, 0.2, 0.06397225788628423 +0.0, 0.3, 0.17949183411348948 +0.0, 0.4, 0.21735888944008888 +0.0, 0.5, 0.28586798898185645 +0.0, 0.6, 0.44321605368239025 +0.0, 0.7, 0.5672744690405076 +0.0, 0.8, 0.6574109301389564 +0.0, 0.9, 0.8143421958896492 +0.0, 1.0, 1.0522341457992452 +0.1, 0.0, 0.09052210158556362 +0.1, 0.1, 0.1009520216061986 +0.1, 0.2, 0.10171532988599352 +0.1, 0.3, 0.13722616611324548 +0.1, 0.4, 0.19718294890244215 +0.1, 0.5, 0.26755230776393085 +0.1, 0.6, 0.42994214465458563 +0.1, 0.7, 0.5255749028209615 +0.1, 0.8, 0.7054404944382558 +0.1, 0.9, 0.8195068120824216 +0.1, 1.0, 1.0986800254920859 +0.2, 0.0, 0.07391924013411105 +0.2, 0.1, 0.08143005606383755 +0.2, 0.2, 0.1331007397675232 +0.2, 0.3, 0.13154620011587023 +0.2, 0.4, 0.21591357533769745 +0.2, 0.5, 0.3396374379061799 +0.2, 0.6, 0.4256001526659258 +0.2, 0.7, 0.5260797691897607 +0.2, 0.8, 0.6554675903724397 +0.2, 0.9, 0.8199487228774796 +0.2, 1.0, 1.0305368537989132 +0.3, 0.0, 0.09065658438028218 +0.3, 0.1, 0.010718818140496137 +0.3, 0.2, 0.13725083785741932 +0.3, 0.3, 0.14769062988342788 +0.3, 0.4, 0.16867615523135612 +0.3, 0.5, 0.34100841301260043 +0.3, 0.6, 0.45602397556766017 +0.3, 0.7, 0.5272890111755907 +0.3, 0.8, 0.6697053819387533 +0.3, 0.9, 0.8346036090905168 +0.3, 1.0, 1.0261608803455926 +0.4, 0.0, 0.07562067903495283 +0.4, 0.1, 0.03778027951177275 +0.4, 0.2, 0.04641740213439717 +0.4, 0.3, 0.16419673637538534 +0.4, 0.4, 0.22312679707699493 +0.4, 0.5, 0.2630522620543643 +0.4, 0.6, 0.4371176818339449 +0.4, 0.7, 0.5723323115053733 +0.4, 0.8, 0.7080551108616417 +0.4, 0.9, 0.8279367483021247 +0.4, 1.0, 1.0566861425416065 +0.5, 0.0, 0.0313041526807114 +0.5, 0.1, 0.09261693968504361 +0.5, 0.2, 0.1366180358573751 +0.5, 0.3, 0.17997862719957242 +0.5, 0.4, 0.16191580585589166 +0.5, 0.5, 0.31836821823205896 +0.5, 0.6, 0.42233901538281543 +0.5, 0.7, 0.5600346360197604 +0.5, 0.8, 0.6992482325761759 +0.5, 0.9, 0.862258365787062 +0.5, 1.0, 1.0378091868392207 +0.6, 0.0, 0.0988068010181403 +0.6, 0.1, 0.03773387710194175 +0.6, 0.2, 0.10176653979089584 +0.6, 0.3, 0.14887220116919433 +0.6, 0.4, 0.25252497778669736 +0.6, 0.5, 0.3055300354904886 +0.6, 0.6, 0.3917340248489845 +0.6, 0.7, 0.5525997648069417 +0.6, 0.8, 0.6696021355843519 +0.6, 0.9, 0.8988929370024403 +0.6, 1.0, 1.0981723152875262 +0.7, 0.0, 0.013481301649133704 +0.7, 0.1, 0.10330920184602638 +0.7, 0.2, 0.08754342206405685 +0.7, 0.3, 0.17203402265475137 +0.7, 0.4, 0.21355416005595512 +0.7, 0.5, 0.32423032390355605 +0.7, 0.6, 0.4004330635685615 +0.7, 0.7, 0.5766931867856282 +0.7, 0.8, 0.6528643584003283 +0.7, 0.9, 0.8429171716399838 +0.7, 1.0, 1.0616458248980662 +0.8, 0.0, 0.033106729829805406 +0.8, 0.1, 0.023291980971513665 +0.8, 0.2, 0.07185018965121617 +0.8, 0.3, 0.12165196581499371 +0.8, 0.4, 0.2448963321019051 +0.8, 0.5, 0.2707808035796927 +0.8, 0.6, 0.43494661906860166 +0.8, 0.7, 0.5014935299134817 +0.8, 0.8, 0.7352837904543503 +0.8, 0.9, 0.848317112179193 +0.8, 1.0, 1.0429381820645305 +0.9, 0.0, 0.09130352242808197 +0.9, 0.1, 0.06657183579297654 +0.9, 0.2, 0.11223658384915326 +0.9, 0.3, 0.1364353366977988 +0.9, 0.4, 0.19712838734472587 +0.9, 0.5, 0.3051728121945353 +0.9, 0.6, 0.3749660409901386 +0.9, 0.7, 0.5190051514823596 +0.9, 0.8, 0.6774512640919589 +0.9, 0.9, 0.8310034058242509 +0.9, 1.0, 1.0689673662128705 +1.0, 0.0, 0.009062810309745007 +1.0, 0.1, 0.05327474257878392 +1.0, 0.2, 0.12975922343712176 +1.0, 0.3, 0.11834110855807276 +1.0, 0.4, 0.1980660815725522 +1.0, 0.5, 0.346034727499414 +1.0, 0.6, 0.4300237539110094 +1.0, 0.7, 0.5368342913405174 +1.0, 0.8, 0.695356382862654 +1.0, 0.9, 0.8624470450842356 +1.0, 1.0, 1.053732774959407 \ No newline at end of file diff --git a/test/test_kazda_li.jl b/test/test_kazda_li.jl index 59b498c..e209cec 100644 --- a/test/test_kazda_li.jl +++ b/test/test_kazda_li.jl @@ -18,19 +18,18 @@ end # 2D - g(x) = x[1]^2 + x[2]^2 + g(x, y) = x^2 + y^2 vals = PiecewiseAffineApprox._sample_uniform(g, [(-1, 1), (-1, 1)], 10) pwa_red = approx( vals, Convex(), Progressive(optimizer = optimizer, tolerance = 0.2, metric = :max), ) - # Disabled as it gives inconsistent results locally and in CI - # @test length(pwa_red.planes) == 10 + @test length(pwa_red.planes) == 10 @test evaluate(pwa_red, (0, 0)) ≈ 0.024 atol = 0.001 @testset "Nonconvex" begin - h(x) = sin(5 * x[1]) * (x[1]^2 + x[2]^2) + h(x, y) = sin(5 * x) * (x^2 + y^2) vals = PiecewiseAffineApprox._sample_uniform(h, [(-1, 1), (-1, 1)], 10) @test_throws ErrorException approx( @@ -89,7 +88,7 @@ end end # 2D - g(x) = x[1]^2 + x[2]^2 + g(x, y) = x^2 + y^2 vals = PiecewiseAffineApprox._sample_uniform(g, [(-1, 1), (-1, 1)], 10) pwa_red = approx(vals, Convex(), FullOrder(optimizer = optimizer, metric = :max)) @@ -97,7 +96,7 @@ end @test evaluate(pwa_red, (0, 0)) ≈ 0.024 atol = 0.001 @testset "Nonconvex" begin - h(x) = sin(5 * x[1]) * (x[1]^2 + x[2]^2) + h(x, y) = sin(5 * x) * (x^2 + y^2) vals = PiecewiseAffineApprox._sample_uniform(h, [(-1, 1), (-1, 1)], 10) @test_throws ErrorException approx(