From 332139d70dae9e45cbadb750db6fbaa7e599115f Mon Sep 17 00:00:00 2001
From: Lars Hellemo <Lars.Hellemo@sintef.no>
Date: Mon, 21 Oct 2024 21:48:17 +0200
Subject: [PATCH 1/8] Substitute @info with @debug

---
 src/kazda_li.jl     | 6 +++---
 src/magnani_boyd.jl | 4 ++--
 2 files changed, 5 insertions(+), 5 deletions(-)

diff --git a/src/kazda_li.jl b/src/kazda_li.jl
index cd139ac..e7639c1 100644
--- a/src/kazda_li.jl
+++ b/src/kazda_li.jl
@@ -94,13 +94,13 @@ function convexify(
 
     obj_val = objective_value(m)
     if obj_val == 0
-        @info "Data points are $(curv_string)"
+        @debug "Data points are $(curv_string)"
         return f
     end
 
     pts_adj = count(v -> v > 0, value.(s⁺ + s⁻))
 
-    @info "Data are not $(curv_string), dev = $(round(obj_val; digits=3)), $(pts_adj) point(s) adjusted"
+    @debug "Data are not $(curv_string), dev = $(round(obj_val; digits=3)), $(pts_adj) point(s) adjusted"
     values = [f.values[l] + value(s⁺[l]) - value(s⁻[l]) for l ∈ 1:K]
     return FunctionEvaluations(f.points, values)
 end
@@ -280,7 +280,7 @@ function _progressive_pwa(
         p += 1
         pwa_red = _allocation_improvement(f, p, U, optimizer)
     end
-    @info "Fitting finished, error = $(round(_approx_error(f, pwa_red, metric); digits = 3)), p = $p"
+    @debug "Fitting finished, error = $(round(_approx_error(f, pwa_red, metric); digits = 3)), p = $p"
     return pwa_red
 end
 
diff --git a/src/magnani_boyd.jl b/src/magnani_boyd.jl
index 614812c..3bf1331 100644
--- a/src/magnani_boyd.jl
+++ b/src/magnani_boyd.jl
@@ -54,7 +54,7 @@ function _convex_linearization_mb(X::Matrix, z::Vector, options::Cluster)
     strict = options.strict
     optimizer = options.optimizer
 
-    @info "Starting heuristic search "
+    @debug "Starting heuristic search "
     approxes = collect(
         fetch.(
             @spawn _convex_linearization_mb_single(
@@ -70,7 +70,7 @@ function _convex_linearization_mb(X::Matrix, z::Vector, options::Cluster)
     )
     min_error, pwa_best = argmin(first, approxes)
 
-    @info "Terminating search - best approximation error = $(min_error) ($metric)"
+    @debug "Terminating search - best approximation error = $(min_error) ($metric)"
     return pwa_best
 end
 

From 3c7884279f29e95cad4b053a53dbbcf6a32c6f2b Mon Sep 17 00:00:00 2001
From: Lars Hellemo <Lars.Hellemo@sintef.no>
Date: Mon, 21 Oct 2024 22:09:23 +0200
Subject: [PATCH 2/8] Prep v.0.6.2

---
 Project.toml | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/Project.toml b/Project.toml
index 70f1f43..e569615 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "PiecewiseAffineApprox"
 uuid = "029c5f1c-0dc4-4551-8b87-a10b73fe3713"
 authors = ["Truls Flatberg <Truls.Flatberg@sintef.no>, Lars Hellemo <Lars.Hellemo@sintef.no>, Thiago Silva <Thiago.Silva.sintef.no>"]
-version = "0.6.1"
+version = "0.6.2"
 
 [deps]
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
@@ -24,7 +24,7 @@ MakieExt = "Makie"
 Combinatorics = "1"
 Distributions = "0.25"
 JuMP = "1"
-julia = "^1.9"
+julia = "1.10"
 Makie = "0.21"
 StableRNGs = "1"
 Statistics = "1"

From 27bd5b363082263b1d117e22f38b51dfa115dfdd Mon Sep 17 00:00:00 2001
From: Lars Hellemo <Lars.Hellemo@sintef.no>
Date: Mon, 21 Oct 2024 22:13:59 +0200
Subject: [PATCH 3/8] Update CI

---
 .github/workflows/ci.yml | 11 ++++-------
 1 file changed, 4 insertions(+), 7 deletions(-)

diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index d15d7ac..361100a 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -16,18 +16,15 @@ jobs:
         # Since PiecewiseAffineApprox doesn't have binary dependencies, only test on a subset
         # of possible platforms.
         include:
-          - version: '1'  # The latest point-release (Linux)
+          - version: ['lts', '1']  # LTS + latest point-release (Linux)
             os: ubuntu-latest
             arch: x64
-          - version: '1'  # The latest point-release (Windows)
+          - version: ['lts', '1']  # LTS + latest point-release (Windows)
             os: windows-latest
             arch: x64
-          - version: '1.9'  # 1.9, not 1.6 LTS (64-bit Linux)
-            os: ubuntu-latest
-            arch: x64
     steps:
-      - uses: actions/checkout@v2
-      - uses: julia-actions/setup-julia@v1
+      - uses: actions/checkout@v3
+      - uses: julia-actions/setup-julia@v2
         with:
           version: ${{ matrix.version }}
           arch: ${{ matrix.arch }}

From 61a00e0c7a5de5856763d6980c521b5a5f5e86f9 Mon Sep 17 00:00:00 2001
From: Lars Hellemo <Lars.Hellemo@sintef.no>
Date: Mon, 21 Oct 2024 22:23:57 +0200
Subject: [PATCH 4/8] fix test matrix

---
 .github/workflows/ci.yml | 12 +++---------
 1 file changed, 3 insertions(+), 9 deletions(-)

diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 361100a..97b4da3 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -13,15 +13,9 @@ jobs:
     strategy:
       fail-fast: false
       matrix:
-        # Since PiecewiseAffineApprox doesn't have binary dependencies, only test on a subset
-        # of possible platforms.
-        include:
-          - version: ['lts', '1']  # LTS + latest point-release (Linux)
-            os: ubuntu-latest
-            arch: x64
-          - version: ['lts', '1']  # LTS + latest point-release (Windows)
-            os: windows-latest
-            arch: x64
+        version: ['lts', '1']
+        os: [ubuntu-latest, windows-latest]
+        arch: [x64]
     steps:
       - uses: actions/checkout@v3
       - uses: julia-actions/setup-julia@v2

From ba6812865228e4d01610b5811a0f9d5d0e5177c0 Mon Sep 17 00:00:00 2001
From: Lars Hellemo <Lars.Hellemo@sintef.no>
Date: Tue, 22 Oct 2024 00:22:58 +0200
Subject: [PATCH 5/8] Add yield() to _reduced_order_pwa to fix hang

---
 src/kazda_li.jl | 590 ++++++++++++++++++++++++------------------------
 1 file changed, 295 insertions(+), 295 deletions(-)

diff --git a/src/kazda_li.jl b/src/kazda_li.jl
index e7639c1..4c2c075 100644
--- a/src/kazda_li.jl
+++ b/src/kazda_li.jl
@@ -1,295 +1,295 @@
-#=
-This file holds an implementation of a progressive fitting heuristic
-inspired by the algorithm in the paper "A linear programming approach to
-difference-of-convex piecewise linear approximation" by Kazda and Li (2024).
-
-The current implementation requires that the function evaluations are
-samples of a convex function, i.e. it is not robust with regard to noisy
-and non-convex data.
-
-The main algorithm consists of the following steps:
-1. While error < tolerance
-    a. Introduce an extra plane in the pwa approximation
-    b. Find local optimum for allocating points to plane (using linear programming)
-=#
-
-"""
-    enforce_curvature(f::FunctionEvaluations, curvature::Curvature, optimizer, metric = :l1)
-
-Create a slightly perturbed version of the function values
-to ensure that the data points can be interpolated by a convex/concave
-piecewise affine function.
-
-Enforcing the curvature is performed using a linear optimization problem
-that adjust the function values and tries to minimize the total deviation.
-The total deviation can be measured in different metrics specified
-by the `metric` parameter.
-This function can be useful as a pre-step for running the full-order and progressive
-fitting heuristics that require data points that are convex/concave in this sense.
-"""
-function enforce_curvature(
-    f::FunctionEvaluations,
-    c::Convex,
-    optimizer,
-    metric = :l1,
-)
-    return convexify(f, optimizer, metric)
-end
-
-function enforce_curvature(
-    f::FunctionEvaluations,
-    c::Concave,
-    optimizer,
-    metric = :l1,
-)
-    fneg = FunctionEvaluations(f.points, -f.values)
-    fconvex = convexify(fneg, optimizer, metric, "concave")
-    return FunctionEvaluations(f.points, -fconvex.values)
-end
-
-function convexify(
-    f::FunctionEvaluations{D},
-    optimizer,
-    metric = :l1,
-    curv_string = "convex",
-) where {D}
-    K = length(f)
-
-    m = Model(optimizer)
-    @variable(m, a[1:K, 1:D])
-    @variable(m, b[1:K])
-    @variable(m, s⁺[1:K] ≥ 0)
-    @variable(m, s⁻[1:K] ≥ 0)
-
-    for (l, (x, z)) ∈ enumerate(f)
-        for k ∈ 1:K
-            @constraint(m, z + s⁺[l] - s⁻[l] ≥ dot(a[k, :], x) + b[k])
-        end
-        @constraint(m, z + s⁺[l] - s⁻[l] ≤ dot(a[l, :], x) + b[l])
-    end
-    obj = AffExpr()
-    if metric == :l2 || metric == :rms
-        obj = sum(s⁺[l]^2 + s⁻[l]^2 for l ∈ 1:K)
-    elseif metric == :max
-        @variable(m, smax)
-        obj = smax
-        for l ∈ 1:K
-            @constraint(m, smax ≥ s⁺[l] + s⁻[l])
-        end
-    elseif metric == :l1
-        obj = sum(s⁺[l] + s⁻[l] for l ∈ 1:K)
-    elseif metric == :none
-        # No regularization term added
-    else
-        error("Unrecognized/unsupported metric type $(metric)")
-    end
-    @objective(m, Min, obj)
-
-    optimize!(m)
-
-    if !is_solved_and_feasible(m)
-        #JuMP.write_to_file(m, "convexify.lp")
-        error("Unable to find a feasible solution")
-    end
-
-    obj_val = objective_value(m)
-    if obj_val == 0
-        @debug "Data points are $(curv_string)"
-        return f
-    end
-
-    pts_adj = count(v -> v > 0, value.(s⁺ + s⁻))
-
-    @debug "Data are not $(curv_string), dev = $(round(obj_val; digits=3)), $(pts_adj) point(s) adjusted"
-    values = [f.values[l] + value(s⁺[l]) - value(s⁻[l]) for l ∈ 1:K]
-    return FunctionEvaluations(f.points, values)
-end
-
-# Create a full order convex approximation for the given data points,
-# i.e. using a separate plane for each data point.
-function _full_order_pwa(f::FunctionEvaluations{D}, optimizer, metric) where {D}
-    K = length(f)
-
-    m = Model(optimizer)
-    @variable(m, a[1:K, 1:D])
-    @variable(m, b[1:K])
-    @variable(m, s[1:K, 1:K] ≥ 0)
-    for (l, (x, z)) ∈ enumerate(f)
-        for k ∈ 1:K
-            @constraint(m, z - s[l, k] == dot(a[k, :], x) + b[k])
-        end
-        @constraint(m, z ≤ dot(a[l, :], x) + b[l])
-    end
-    obj = AffExpr()
-    if metric == :l2 || metric == :rms
-        obj = sum(s .^ 2)
-    elseif metric == :max
-        @variable(m, smax)
-        obj = smax
-        @constraint(m, s .≤ smax)
-    elseif metric == :l1
-        obj = sum(s)
-    elseif metric == :none
-        # No regularization added to m
-    else
-        error("Unrecognized/unsupported metric type $(metric)")
-    end
-    @objective(m, Min, obj)
-    optimize!(m)
-
-    if !is_solved_and_feasible(m)
-        error(
-            "No full order convex approximation exists. Function evaluations should be " *
-            "from a convex function.",
-        )
-    end
-
-    aᶠ = value.(a)
-    bᶠ = value.(b)
-    pwa = PWAFunc([Plane(aᶠ[k, :], bᶠ[k]) for k ∈ 1:K], Convex())
-    return pwa
-end
-
-# Create a pwa approximation using `p` planes where `U` provides a mapping from
-# data points indices to the indices of planes. The planes are required to lie below
-# all data points. Points are allowed to lie above their associated plane with
-# a metric that is minimized in the optimization.
-function _reduced_order_pwa(
-    f::FunctionEvaluations{D},
-    p::Int,
-    U,
-    optimizer,
-    μ = 1e-4,
-) where {D}
-    K = length(f)
-
-    m = Model(optimizer)
-    @variable(m, ared[1:p, 1:D])
-    @variable(m, bred[1:p])
-    @variable(m, s[1:K] ≥ 0)
-    @variable(m, sʳ[1:K, 1:p] ≥ 0)
-    for (l, (x, z)) ∈ enumerate(f)
-        for k ∈ 1:p
-            @constraint(m, z - sʳ[l, k] == dot(ared[k, :], x) + bred[k])
-        end
-    end
-    for (l, k) ∈ U
-        z = f.values[l]
-        x = f.points[l]
-        @constraint(m, z ≤ dot(ared[k, :], x) + bred[k] + s[l])
-    end
-    @objective(m, Min, sum(s) + μ / p * sum(sʳ))
-    optimize!(m)
-    obj = objective_value(m)
-
-    aʳ = value.(ared)
-    bʳ = value.(bred)
-    pwa_red = PWAFunc([Plane(aʳ[k, :], bʳ[k]) for k ∈ 1:p], Convex())
-
-    return obj, pwa_red
-end
-
-# The hyperplane being active at the point x for a convex pwa function
-function _active_plane(pwl::PWAFunc{Convex,D}, x) where {D}
-    return argmax(collect(evaluate(p, x) for p ∈ pwl.planes))
-end
-
-# Map data points to the active plane for the piecewise approximation
-function _update_allocation(f::FunctionEvaluations, pwa_red)
-    U = [(l, _active_plane(pwa_red, f.points[l])) for l ∈ 1:length(f)]
-    return U
-end
-
-# Do a local improvement of point allocation to planes until error
-# is not further reduced
-function _allocation_improvement(f::FunctionEvaluations, p, Uᴵ, optimizer)
-    U = Uᴵ
-    prev_obj = Inf64
-    obj, pwa_red = _reduced_order_pwa(f, p, U, optimizer)
-    while obj < prev_obj
-        prev_obj = obj
-        U = _update_allocation(f, pwa_red)
-        obj, pwa_red = _reduced_order_pwa(f, p, U, optimizer)
-    end
-    return pwa_red
-end
-
-# Add an extra plane to the pwa by identifying the point lying
-# further away and use the corresponding plane in the full order pwa
-# to update the allocation of points to planes.
-function _increase_order(f::FunctionEvaluations, pwa_red, pwa, used)
-    s = [z - evaluate(pwa_red, x) for (x, z) ∈ f]
-    imax = 0
-    smax = -1
-    for i ∈ 1:length(f)
-        if s[i] > smax && !(i in used)
-            imax = i
-            smax = s[i]
-        end
-    end
-    push!(pwa_red.planes, pwa.planes[imax])
-    U = _update_allocation(f, pwa_red)
-    return U, vcat(used, imax)
-end
-
-# Calculate the approximation error between input data points and the piecewise
-# affine approximation in different metrics
-function _approx_error(f::FunctionEvaluations, pwa::PWAFunc, metric = :l1)
-    err = 0.0
-    for (x, z) ∈ f
-        v = evaluate(pwa, x)
-        if metric == :l1
-            err += abs(v - z)
-        elseif metric == :l2 || metric == :rms
-            err += (v - z)^2
-        elseif metric == :max
-            err = max(err, abs(v - z))
-        end
-    end
-    if metric == :rms
-        err = err / length(f.points)
-    end
-    if metric == :l2 || metric == :rms
-        err = sqrt(err)
-    end
-    return err
-end
-
-# Main algorithm
-function _progressive_pwa(
-    f::FunctionEvaluations,
-    optimizer,
-    δᵗᵒˡ = 1e-3,
-    metric = :max,
-    full_order_metric = :l1,
-)
-
-    # Find the full order convex approximation, will error if the points are
-    # not convex
-    pwa = _full_order_pwa(f, optimizer, full_order_metric)
-
-    # Start with an approximation with one plane and allocate all points to that segment
-    p = 1
-    U = [(l, 1) for l ∈ 1:length(f)]
-
-    # Increase the number of planes until the required tolerance is met
-    pwa_red = _allocation_improvement(f, p, U, optimizer)
-    used = []
-    while _approx_error(f, pwa_red, metric) > δᵗᵒˡ && p < length(f.points)
-        U, used = _increase_order(f, pwa_red, pwa, used)
-        p += 1
-        pwa_red = _allocation_improvement(f, p, U, optimizer)
-    end
-    @debug "Fitting finished, error = $(round(_approx_error(f, pwa_red, metric); digits = 3)), p = $p"
-    return pwa_red
-end
-
-# Map algorithm structure to correct parameters
-function _progressive_pwa(f::FunctionEvaluations, options::Progressive)
-    return _progressive_pwa(
-        f,
-        options.optimizer,
-        options.tolerance,
-        options.metric,
-    )
-end
+#=
+This file holds an implementation of a progressive fitting heuristic
+inspired by the algorithm in the paper "A linear programming approach to
+difference-of-convex piecewise linear approximation" by Kazda and Li (2024).
+
+The current implementation requires that the function evaluations are
+samples of a convex function, i.e. it is not robust with regard to noisy
+and non-convex data.
+
+The main algorithm consists of the following steps:
+1. While error < tolerance
+    a. Introduce an extra plane in the pwa approximation
+    b. Find local optimum for allocating points to plane (using linear programming)
+=#
+
+"""
+    enforce_curvature(f::FunctionEvaluations, curvature::Curvature, optimizer, metric = :l1)
+
+Create a slightly perturbed version of the function values
+to ensure that the data points can be interpolated by a convex/concave
+piecewise affine function.
+
+Enforcing the curvature is performed using a linear optimization problem
+that adjust the function values and tries to minimize the total deviation.
+The total deviation can be measured in different metrics specified
+by the `metric` parameter.
+This function can be useful as a pre-step for running the full-order and progressive
+fitting heuristics that require data points that are convex/concave in this sense.
+"""
+function enforce_curvature(
+    f::FunctionEvaluations,
+    c::Convex,
+    optimizer,
+    metric = :l1,
+)
+    return convexify(f, optimizer, metric)
+end
+
+function enforce_curvature(
+    f::FunctionEvaluations,
+    c::Concave,
+    optimizer,
+    metric = :l1,
+)
+    fneg = FunctionEvaluations(f.points, -f.values)
+    fconvex = convexify(fneg, optimizer, metric, "concave")
+    return FunctionEvaluations(f.points, -fconvex.values)
+end
+
+function convexify(
+    f::FunctionEvaluations{D},
+    optimizer,
+    metric = :l1,
+    curv_string = "convex",
+) where {D}
+    K = length(f)
+
+    m = Model(optimizer)
+    @variable(m, a[1:K, 1:D])
+    @variable(m, b[1:K])
+    @variable(m, s⁺[1:K] ≥ 0)
+    @variable(m, s⁻[1:K] ≥ 0)
+
+    for (l, (x, z)) ∈ enumerate(f)
+        for k ∈ 1:K
+            @constraint(m, z + s⁺[l] - s⁻[l] ≥ dot(a[k, :], x) + b[k])
+        end
+        @constraint(m, z + s⁺[l] - s⁻[l] ≤ dot(a[l, :], x) + b[l])
+    end
+    obj = AffExpr()
+    if metric == :l2 || metric == :rms
+        obj = sum(s⁺[l]^2 + s⁻[l]^2 for l ∈ 1:K)
+    elseif metric == :max
+        @variable(m, smax)
+        obj = smax
+        for l ∈ 1:K
+            @constraint(m, smax ≥ s⁺[l] + s⁻[l])
+        end
+    elseif metric == :l1
+        obj = sum(s⁺[l] + s⁻[l] for l ∈ 1:K)
+    elseif metric == :none
+        # No regularization term added
+    else
+        error("Unrecognized/unsupported metric type $(metric)")
+    end
+    @objective(m, Min, obj)
+
+    optimize!(m)
+
+    if !is_solved_and_feasible(m)
+        #JuMP.write_to_file(m, "convexify.lp")
+        error("Unable to find a feasible solution")
+    end
+
+    obj_val = objective_value(m)
+    if obj_val == 0
+        @debug "Data points are $(curv_string)"
+        return f
+    end
+
+    pts_adj = count(v -> v > 0, value.(s⁺ + s⁻))
+
+    @debug "Data are not $(curv_string), dev = $(round(obj_val; digits=3)), $(pts_adj) point(s) adjusted"
+    values = [f.values[l] + value(s⁺[l]) - value(s⁻[l]) for l ∈ 1:K]
+    return FunctionEvaluations(f.points, values)
+end
+
+# Create a full order convex approximation for the given data points,
+# i.e. using a separate plane for each data point.
+function _full_order_pwa(f::FunctionEvaluations{D}, optimizer, metric) where {D}
+    K = length(f)
+
+    m = Model(optimizer)
+    @variable(m, a[1:K, 1:D])
+    @variable(m, b[1:K])
+    @variable(m, s[1:K, 1:K] ≥ 0)
+    for (l, (x, z)) ∈ enumerate(f)
+        for k ∈ 1:K
+            @constraint(m, z - s[l, k] == dot(a[k, :], x) + b[k])
+        end
+        @constraint(m, z ≤ dot(a[l, :], x) + b[l])
+    end
+    obj = AffExpr()
+    if metric == :l2 || metric == :rms
+        obj = sum(s .^ 2)
+    elseif metric == :max
+        @variable(m, smax)
+        obj = smax
+        @constraint(m, s .≤ smax)
+    elseif metric == :l1
+        obj = sum(s)
+    elseif metric == :none
+        # No regularization added to m
+    else
+        error("Unrecognized/unsupported metric type $(metric)")
+    end
+    @objective(m, Min, obj)
+    optimize!(m)
+
+    if !is_solved_and_feasible(m)
+        error(
+            "No full order convex approximation exists. Function evaluations should be " *
+            "from a convex function.",
+        )
+    end
+
+    aᶠ = value.(a)
+    bᶠ = value.(b)
+    pwa = PWAFunc([Plane(aᶠ[k, :], bᶠ[k]) for k ∈ 1:K], Convex())
+    return pwa
+end
+
+# Create a pwa approximation using `p` planes where `U` provides a mapping from
+# data points indices to the indices of planes. The planes are required to lie below
+# all data points. Points are allowed to lie above their associated plane with
+# a metric that is minimized in the optimization.
+function _reduced_order_pwa(
+    f::FunctionEvaluations{D},
+    p::Int,
+    U,
+    optimizer,
+    μ = 1e-4,
+) where {D}
+    yield()
+    K = length(f)
+    m = Model(optimizer)
+    @variable(m, ared[1:p, 1:D])
+    @variable(m, bred[1:p])
+    @variable(m, s[1:K] ≥ 0)
+    @variable(m, sʳ[1:K, 1:p] ≥ 0)
+    for (l, (x, z)) ∈ enumerate(f)
+        for k ∈ 1:p
+            @constraint(m, z - sʳ[l, k] == dot(ared[k, :], x) + bred[k])
+        end
+    end
+    for (l, k) ∈ U
+        z = f.values[l]
+        x = f.points[l]
+        @constraint(m, z ≤ dot(ared[k, :], x) + bred[k] + s[l])
+    end
+    @objective(m, Min, sum(s) + μ / p * sum(sʳ))
+    optimize!(m)
+    obj = objective_value(m)
+
+    aʳ = value.(ared)
+    bʳ = value.(bred)
+    pwa_red = PWAFunc([Plane(aʳ[k, :], bʳ[k]) for k ∈ 1:p], Convex())
+
+    return obj, pwa_red
+end
+
+# The hyperplane being active at the point x for a convex pwa function
+function _active_plane(pwl::PWAFunc{Convex,D}, x) where {D}
+    return argmax(collect(evaluate(p, x) for p ∈ pwl.planes))
+end
+
+# Map data points to the active plane for the piecewise approximation
+function _update_allocation(f::FunctionEvaluations, pwa_red)
+    U = [(l, _active_plane(pwa_red, f.points[l])) for l ∈ 1:length(f)]
+    return U
+end
+
+# Do a local improvement of point allocation to planes until error
+# is not further reduced
+function _allocation_improvement(f::FunctionEvaluations, p, Uᴵ, optimizer)
+    U = Uᴵ
+    prev_obj = Inf64
+    obj, pwa_red = _reduced_order_pwa(f, p, U, optimizer)
+    while obj < prev_obj
+        prev_obj = obj
+        U = _update_allocation(f, pwa_red)
+        obj, pwa_red = _reduced_order_pwa(f, p, U, optimizer)
+    end
+    return pwa_red
+end
+
+# Add an extra plane to the pwa by identifying the point lying
+# further away and use the corresponding plane in the full order pwa
+# to update the allocation of points to planes.
+function _increase_order(f::FunctionEvaluations, pwa_red, pwa, used)
+    s = [z - evaluate(pwa_red, x) for (x, z) ∈ f]
+    imax = 0
+    smax = -1
+    for i ∈ 1:length(f)
+        if s[i] > smax && !(i in used)
+            imax = i
+            smax = s[i]
+        end
+    end
+    push!(pwa_red.planes, pwa.planes[imax])
+    U = _update_allocation(f, pwa_red)
+    return U, vcat(used, imax)
+end
+
+# Calculate the approximation error between input data points and the piecewise
+# affine approximation in different metrics
+function _approx_error(f::FunctionEvaluations, pwa::PWAFunc, metric = :l1)
+    err = 0.0
+    for (x, z) ∈ f
+        v = evaluate(pwa, x)
+        if metric == :l1
+            err += abs(v - z)
+        elseif metric == :l2 || metric == :rms
+            err += (v - z)^2
+        elseif metric == :max
+            err = max(err, abs(v - z))
+        end
+    end
+    if metric == :rms
+        err = err / length(f.points)
+    end
+    if metric == :l2 || metric == :rms
+        err = sqrt(err)
+    end
+    return err
+end
+
+# Main algorithm
+function _progressive_pwa(
+    f::FunctionEvaluations,
+    optimizer,
+    δᵗᵒˡ = 1e-3,
+    metric = :max,
+    full_order_metric = :l1,
+)
+
+    # Find the full order convex approximation, will error if the points are
+    # not convex
+    pwa = _full_order_pwa(f, optimizer, full_order_metric)
+
+    # Start with an approximation with one plane and allocate all points to that segment
+    p = 1
+    U = [(l, 1) for l ∈ 1:length(f)]
+
+    # Increase the number of planes until the required tolerance is met
+    pwa_red = _allocation_improvement(f, p, U, optimizer)
+    used = []
+    while _approx_error(f, pwa_red, metric) > δᵗᵒˡ && p < length(f.points)
+        U, used = _increase_order(f, pwa_red, pwa, used)
+        p += 1
+        pwa_red = _allocation_improvement(f, p, U, optimizer)
+    end
+    @debug "Fitting finished, error = $(round(_approx_error(f, pwa_red, metric); digits = 3)), p = $p"
+    return pwa_red
+end
+
+# Map algorithm structure to correct parameters
+function _progressive_pwa(f::FunctionEvaluations, options::Progressive)
+    return _progressive_pwa(
+        f,
+        options.optimizer,
+        options.tolerance,
+        options.metric,
+    )
+end

From ba8191bb586a9c06cba4a924f24481ebcde74a79 Mon Sep 17 00:00:00 2001
From: Lars Hellemo <Lars.Hellemo@sintef.no>
Date: Tue, 22 Oct 2024 09:15:42 +0200
Subject: [PATCH 6/8] Add comment explaining the yield()

---
 src/kazda_li.jl | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/kazda_li.jl b/src/kazda_li.jl
index 4c2c075..dfd77b9 100644
--- a/src/kazda_li.jl
+++ b/src/kazda_li.jl
@@ -161,7 +161,8 @@ function _reduced_order_pwa(
     optimizer,
     μ = 1e-4,
 ) where {D}
-    yield()
+    yield() # Avoid blocking of main thread causing hangs in some cases, see
+            # https://github.com/sintefore/PiecewiseAffineApprox.jl/issues/11
     K = length(f)
     m = Model(optimizer)
     @variable(m, ared[1:p, 1:D])

From e8b5f0b7543c6a56c7c5ccc810a7e8e4d61364a3 Mon Sep 17 00:00:00 2001
From: Lars Hellemo <Lars.Hellemo@sintef.no>
Date: Tue, 22 Oct 2024 09:48:46 +0200
Subject: [PATCH 7/8] Attempt to fix warnings on documentation workflow

---
 .github/workflows/documentation.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.github/workflows/documentation.yml b/.github/workflows/documentation.yml
index 76aff3f..bc20d3d 100644
--- a/.github/workflows/documentation.yml
+++ b/.github/workflows/documentation.yml
@@ -9,7 +9,7 @@ jobs:
   build:
     runs-on: ubuntu-latest
     steps:
-      - uses: actions/checkout@v2
+      - uses: actions/checkout@v3
       - uses: julia-actions/setup-julia@latest
         with:
           version: '1'

From 4454299491d3bbcf6b41db171ba69fd090474f64 Mon Sep 17 00:00:00 2001
From: Lars Hellemo <Lars.Hellemo@sintef.no>
Date: Tue, 22 Oct 2024 09:57:16 +0200
Subject: [PATCH 8/8] Remove defunct codecov workflow

---
 .github/workflows/codecov.yml | 2 --
 1 file changed, 2 deletions(-)
 delete mode 100644 .github/workflows/codecov.yml

diff --git a/.github/workflows/codecov.yml b/.github/workflows/codecov.yml
deleted file mode 100644
index 3217852..0000000
--- a/.github/workflows/codecov.yml
+++ /dev/null
@@ -1,2 +0,0 @@
-codecov:
-  token: 3c305759-ecc3-419d-aee0-20ab9a53d1e9
\ No newline at end of file