Change function specfication for approx

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)
 ```

docs/make.jl

@@ -36,7 +36,6 @@ pages = [
 ]
 
 # TODO: Enable modules
-# TODO: Enable doctests
 Documenter.makedocs(
     sitename = "PiecewiseAffineApprox",
     format = Documenter.HTML(;

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
+```
+
+
+
+

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

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 ?

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()