Add SetProg templates

CLAUDE.md

@@ -70,7 +70,7 @@ src/
 ├── systems/             switched systems, runs, the dual system
 ├── graphs/              queries.jl is THE ADAPTER over HybridSystems; words, predicates,
 │                        languages, dual, de_bruijn, observer, cycles, simulation
-├── templates/           abstract.jl + catalogue/ (six templates; sum_of_squares in ext/)
+├── templates/           abstract.jl + catalogue/ (seven templates; sum_of_squares and set_prog in ext/)
 ├── problems/            abstract.jl + catalogue/ (stability, safety, optimal_control)
 ├── lifts/               COMPARING graphs: abstract.jl, dual.jl, compare.jl + catalogue/
 ├── refinements/         IMPROVING a graph: abstract.jl, dual.jl, strategies.jl, refine.jl + catalogue/

Project.toml

@@ -18,14 +18,17 @@ MathematicalSystems = "0.14"
 Random = "1"
 MultivariatePolynomials = "0.5"
 RecipesBase = "1"
+SetProg = "0.4.4"
 SumOfSquares = "0.7, 0.8"
 julia = "1.10"
 
 [weakdeps]
 MultivariatePolynomials = "102ac46a-7ee4-5c85-9060-abc95bfdeaa3"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
+SetProg = "39881422-4b75-5582-a5c7-0feb14562a65"
 SumOfSquares = "4b9e565b-77fc-50a5-a571-1244f986bda1"
 
 [extensions]
 PathCompleteCertificatesRecipesBaseExt = "RecipesBase"
+PathCompleteCertificatesSetProgExt = "SetProg"
 PathCompleteCertificatesSumOfSquaresExt = ["SumOfSquares", "MultivariatePolynomials"]

docs/Project.toml

@@ -1,24 +1,27 @@
 # The examples are executed as part of the docs build, so this environment needs
-# a solver, a plotting backend and the polynomial stack on top of the
-# documentation tooling. None of them is a dependency of the package, nor of the
-# environment CI instantiates for the tests.
+# a solver, a plotting backend and the polynomial and set-programming stacks on
+# top of the documentation tooling. None of them is a dependency of the package,
+# nor of the environment CI instantiates for the tests.
 [deps]
 Clarabel = "61c947e1-3e6d-4ee4-985a-eec8c727bd6e"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244"
 DynamicPolynomials = "7c1d4256-1411-5781-91ec-d7bc3513ac07"
+HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b"
 HybridSystems = "2207ec0c-686c-5054-b4d2-543502888820"
 JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
 PathCompleteCertificates = "b391c707-6ba7-482f-97f4-86d517cfbcc9"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+SetProg = "39881422-4b75-5582-a5c7-0feb14562a65"
 SumOfSquares = "4b9e565b-77fc-50a5-a571-1244f986bda1"
 
 [compat]
 Documenter = "1"
 DocumenterCitations = "1"
 Literate = "2"
+SetProg = "0.4.4"
 SumOfSquares = "0.7, 0.8"
 
 [sources]

docs/make.jl

@@ -21,6 +21,7 @@ const ORDER = [
     "switched_systems",
     "stability",
     "templates_and_graphs",
+    "set_prog",
     "safety",
     "optimal_control",
     "comparing",
@@ -146,8 +147,7 @@ makedocs(;
             style = :authoryear,
         ),
     ],
-    # The tool paper has six pages including references, so it cannot explain
-    # the package -- these docs have to. `checkdocs = :all` ignores the export
+    # `checkdocs = :all` ignores the export
     # list, so it covers every docstring in the module and *fails the build* on
     # one that no `@autodocs` block under `src/reference/` picks up. Adding a
     # file to `src/` therefore means adding it to a `Pages` list as well.

docs/src/examples/set_prog.jl

@@ -0,0 +1,154 @@
+# # Piecewise semi-ellipsoids: a template from SetProg
+#
+# [SetProg](https://github.com/blegat/SetProg.jl) writes optimisation problems
+# over sets: ellipsoids, sublevel sets of polynomials, polytopes, and sets that
+# are one of these on each cone of a partition. [`SetProgTemplate`](@ref)
+# makes such a set the node function of a node, the function whose unit
+# sublevel set it is, and the stability inequality
+# ``V_d(A_i x) \le \gamma^d V_s(x)`` becomes the inclusion
+# ``(A_i / \gamma)\, S_s \subseteq S_d``, which SetProg imposes. Three of its
+# sets are templates this package already has, and the example checks that;
+# the fourth, the piecewise semi-ellipsoid, is new.
+
+import PathCompleteCertificates as PCC
+import Clarabel, HiGHS, JuMP
+using SetProg
+using DynamicPolynomials
+using LinearAlgebra
+using Plots
+
+const OPTIMIZER = Clarabel.Optimizer
+
+# The rotation and shear of the other examples, on the one-node graph, so that
+# the template alone decides the rate.
+
+A = [0.69 * [0.0 1.0; -1.0 0.0], 0.69 * [1.0 1.0; 0.0 1.0]]
+problem = PCC.StabilityProblem(PCC.switched_system(A))
+graph = PCC.common_graph(problem.system)
+
+bound(template) =
+    PCC.jsr_bound(template, graph, problem; optimizer = OPTIMIZER, rtol = 1e-3)
+
+# A partition is a polytope with the origin inside: its pieces are the cones
+# over its facets. The polytope is built with Polyhedra on a linear programming
+# solver, which SetProg works on the pieces with; a simplex solver such as
+# HiGHS keeps that cheap.
+
+LIBRARY = Polyhedra.DefaultLibrary{Float64}(
+    JuMP.optimizer_with_attributes(HiGHS.Optimizer, JuMP.MOI.Silent() => true),
+)
+
+polygon(normals) = polyhedron(hrep([HalfSpace(a, norm(a)) for a in normals]), LIBRARY)
+
+octagon = polygon([[2, 1], [1, 2], [-1, 2], [-2, 1], [-2, -1], [-1, -2], [1, -2], [2, -1]])
+
+# ## Three sets that are templates already
+#
+# An ellipsoid is a quadratic form, the sublevel set of a polynomial is a sum
+# of squares, and a polytope that is linear on each cone of a partition is the
+# free-facet polyhedral template on the same cones; the cones over the facets
+# of this octagon are those of [`planar_conic_partition`](@ref)`(1)`. Each
+# pair certifies the same rate.
+
+@polyvar x[1:2]
+
+pairs = [
+    (
+        "quadratic",
+        PCC.QuadraticTemplate(),
+        PCC.SetProgTemplate(Ellipsoid(; symmetric = true)),
+    ),
+    (
+        "sum of squares, degree 4",
+        PCC.SumOfSquaresTemplate(2, x),
+        PCC.SetProgTemplate(PolySet(; symmetric = true, degree = 4, variables = x)),
+    ),
+    (
+        "free facets, 4 cones",
+        PCC.ConicPolyhedralTemplate([PCC.planar_conic_partition(1)]),
+        PCC.SetProgTemplate(Polytope(; symmetric = true, piecewise = octagon)),
+    ),
+]
+
+[
+    (name, round(bound(ours).rate; digits = 4), round(bound(theirs).rate; digits = 4)) for
+    (name, ours, theirs) in pairs
+]
+
+# ## The piecewise semi-ellipsoid
+#
+# `Ellipsoid(symmetric = true, piecewise = P)` is a function quadratic on each
+# cone of `P` and continuous across them, a template none of the others gives.
+# Refining the partition can only help, since a function on the coarser one is
+# a function on the finer one too.
+
+square = polygon([[1, 0], [0, 1], [-1, 0], [0, -1]])
+sixteen = polygon([
+    [1, 0],
+    [2, 1],
+    [1, 1],
+    [1, 2],
+    [0, 1],
+    [-1, 2],
+    [-1, 1],
+    [-2, 1],
+    [-1, 0],
+    [-2, -1],
+    [-1, -1],
+    [-1, -2],
+    [0, -1],
+    [1, -2],
+    [1, -1],
+    [2, -1],
+])
+
+partitions = [("square", square), ("octagon", octagon), ("16-gon", sixteen)]
+
+piecewise = [
+    (name, bound(PCC.SetProgTemplate(Ellipsoid(; symmetric = true, piecewise = P)))) for
+    (name, P) in partitions
+]
+
+[
+    (name, nhalfspaces(P), round(certificate.rate; digits = 4)) for
+    ((name, P), (_, certificate)) in zip(partitions, piecewise)
+]
+
+# On the four cones of the square the piecewise function certifies the rate of
+# the plain ellipsoid; the octagon's eight cones improve on it, and sixteen
+# cones of degree two beat the polynomial of degree four above.
+
+# ## Drawing it
+#
+# The sublevel set ``\{x : V(x) \le 1\}`` of each function, normalised as in
+# [Templates and graphs: the two axes](@ref). The piecewise one is a union of
+# arcs of ellipses, meeting along the rays of the partition.
+
+function sublevel_boundary(certificate; samples = 721)
+    degree = PCC.homogeneity_degree(PCC.template(certificate))
+    angles = range(0, 2pi; length = samples)
+    radii = [certificate([cos(a), sin(a)])^(-1 / degree) for a in angles]
+    radii ./= maximum(radii)
+
+    return radii .* cos.(angles), radii .* sin.(angles)
+end
+
+figure = plot(;
+    title = "one node, three partitions",
+    aspect_ratio = :equal,
+    legend = :outerbottom,
+    framestyle = :origin,
+)
+
+for (name, certificate) in piecewise
+    xs, ys = sublevel_boundary(certificate)
+    plot!(
+        figure,
+        xs,
+        ys;
+        label = "$name, rate $(round(certificate.rate; digits = 3))",
+        lw = 2,
+    )
+end
+
+figure

docs/src/index.md

@@ -64,7 +64,9 @@ In reading order:
 
 The two other problems have an example each:
 [Safety: a barrier certificate](@ref) and
-[Optimal control: a feedback per node](@ref).
+[Optimal control: a feedback per node](@ref). A template from another
+package has one too:
+[Piecewise semi-ellipsoids: a template from SetProg](@ref).
 
 ## Scope
 

docs/src/manual/support.md

@@ -5,7 +5,7 @@ CurrentModule = PathCompleteCertificates
 # What works with what
 
 Every template composes with every problem in principle. In practice,
-stability runs on all six templates, and the two other problems on the
+stability runs on all seven templates, and the two other problems on the
 quadratic one.
 
 |                                | [`StabilityProblem`](@ref) | [`SafetyProblem`](@ref) | [`OptimalControlProblem`](@ref) |
@@ -16,11 +16,12 @@ quadratic one.
 | [`PolyhedralTemplate`](@ref)        | ✓ | ✗ | ✗ |
 | [`ConicPolyhedralTemplate`](@ref)   | ✓ | ✗ | ✗ |
 | [`SumOfSquaresTemplate`](@ref)      | ✓ ² | ✗ | ✗ |
+| [`SetProgTemplate`](@ref)           | ✓ ² | ✗ | ✗ |
 
 ¹ Positive systems only: entrywise nonnegative mode matrices.
 
-² Behind a package extension: `using SumOfSquares` first, or its primitives are
-not defined.
+² Behind a package extension: `using SumOfSquares` first, or `using SetProg`
+for [`SetProgTemplate`](@ref), or its primitives are not defined.
 
 A ✗ throws an `ArgumentError` before any model is built. Optimal control needs
 a **complete** graph, not merely a path-complete one. Stability accepts a
@@ -42,8 +43,9 @@ how the loop reads where a certificate is tight:
 | [`PolyhedralTemplate`](@ref) | ✓ | identity; max when the node matrices coincide |
 | [`ConicPolyhedralTemplate`](@ref) | ✓ | identity |
 | [`SumOfSquaresTemplate`](@ref) | ✗ ³ | identity, sum; composition if invertible |
+| [`SetProgTemplate`](@ref) | ✗ ³ | identity; sum when the nodes share the set variable |
 
-³ No `domination_slack` yet, so the loop cannot read where that template is
+³ No `domination_slack` yet, so the loop cannot read where these templates are
 tight.
 
 ## The gaps

docs/src/manual/templates.md

@@ -17,7 +17,7 @@ callable.
     data, [`PolyhedralTemplate`](@ref) holds one matrix per node, and a type
     has nowhere to put it.
 
-## The six templates
+## The seven templates
 
 | Template | ``V_s(x)`` | Homogeneity degree | Solved for | Source |
 | :-- | :-- | :-: | :-- | :-- |
@@ -27,12 +27,20 @@ callable.
 | [`PolyhedralTemplate`](@ref) | ``\max_k \lvert (G_s x)_k \rvert / w_{s,k}`` | 1 | weights ``w_s``, facets ``G_s`` fixed | [athanasopoulos2019polyhedral](@cite) |
 | [`ConicPolyhedralTemplate`](@ref) | ``\max_i \lvert p_{s,i}^\top x \rvert`` | 1 | facets ``p_{s,i}`` | unpublished ¹ |
 | [`SumOfSquaresTemplate`](@ref) | ``z(x)^\top Q_s z(x)``, ``z`` the monomials of degree ``d`` | ``2d`` | ``Q_s \succeq 0`` | [parrilo2008approximation](@cite) ² |
+| [`SetProgTemplate`](@ref) | the gauge of a [SetProg](https://github.com/blegat/SetProg.jl) set | the set's | the set | ³ |
 
 ¹ The free-facet construction is under submission.
 
 ² Behind a package extension: load SumOfSquares first. Raising ``d`` tightens
 the bound monotonically, and ``d = 1`` is exactly [`QuadraticTemplate`](@ref).
 
+³ Behind a package extension: load SetProg first. An ellipsoid is the
+quadratic template, a polynomial set the sum-of-squares one and a piecewise
+polytope the free-facet one on the same cones; what it adds are the piecewise
+semi-ellipsoids, quadratic on each cone of a partition, which improve on the
+quadratic template as the partition is refined;
+[Piecewise semi-ellipsoids: a template from SetProg](@ref) shows all four.
+
 The **homogeneity degree**, [`homogeneity_degree`](@ref), is the ``d`` in
 ``V(cx) = c^d V(x)``. Stability scales its edge inequality by ``\gamma^d``, so
 a rate means the same thing whatever the template, and bounds from different
@@ -69,6 +77,7 @@ hold for it. A template declares its closures with [`is_closed_under`](@ref):
 | [`DualCopositiveTemplate`](@ref) | | ✓ | | |
 | [`PolyhedralTemplate`](@ref) | | ✓ if the node matrices coincide | | |
 | [`ConicPolyhedralTemplate`](@ref) | | | | |
+| [`SetProgTemplate`](@ref) | ✓ if the nodes share the set variable | | | |
 
 ¹ Not closed as a set of functions, yet the min lift is valid for it on
 positive linear systems ([debauche2024thesis](@cite), Thm. 7.43), the one
@@ -85,15 +94,15 @@ was given.
 
 [`has_dual`](@ref), [`dual`](@ref)`(template)` and `dual(template, V)`, the
 dual norm of a fitted function. Quadratic forms are self-dual and the
-copositive templates are each other's dual; the polyhedral and sum-of-squares
-templates have none here. A closure declared on either of a dual pair serves
-both; see [Duality](@ref).
+copositive templates are each other's dual; the polyhedral, sum-of-squares and
+SetProg templates have none here. A closure declared on either of a dual pair
+serves both; see [Duality](@ref).
 
 ## Per-node data
 
 A template whose data is attached to the nodes, such as the facet matrices of
-[`PolyhedralTemplate`](@ref), must say what a node copied by a refinement
-inherits. [`reindex`](@ref)`(template, origins)` gives every node of the new
+[`PolyhedralTemplate`](@ref) or the set variables of [`SetProgTemplate`](@ref),
+must say what a node copied by a refinement inherits. [`reindex`](@ref)`(template, origins)` gives every node of the new
 graph the data of the node it came from; a data-free template returns itself.
 
 ## Writing your own

docs/src/reference/templates.md

@@ -58,3 +58,13 @@ required before this template can be solved with.
 Modules = [PathCompleteCertificates]
 Pages   = ["templates/catalogue/sum_of_squares.jl"]
 ```
+
+## SetProg sets
+
+The primitives are defined in a package extension, so `using SetProg` is
+required before this template can be solved with.
+
+```@autodocs
+Modules = [PathCompleteCertificates]
+Pages   = ["templates/catalogue/set_prog.jl"]
+```