Add Julia Set example (#272)

* Add Julia Set example

* Add ImplicitPlots

* Fixes

* Fixes

* Add disk_integral

* Fix

* Add in_set

* Add LA

* Use sandboxing for examples

* Better viz

* Show solution_summary

* Allow near optimality

Author: blegat

docs/Project.toml

@@ -9,6 +9,7 @@ Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 DynamicPolynomials = "7c1d4256-1411-5781-91ec-d7bc3513ac07"
 GroupsCore = "d5909c97-4eac-4ecc-a3dc-fdd0858a4120"
 HomotopyContinuation = "f213a82b-91d6-5c5d-acf7-10f1c761b327"
+ImplicitPlots = "55ecb840-b828-11e9-1645-43f4a9f9ace7"
 Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9"
 Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
 MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
@@ -18,6 +19,7 @@ MutableArithmetics = "d8a4904e-b15c-11e9-3269-09a3773c0cb0"
 PermutationGroups = "8bc5a954-2dfc-11e9-10e6-cd969bffa420"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
+SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 SumOfSquares = "4b9e565b-77fc-50a5-a571-1244f986bda1"
 SymbolicWedderburn = "858aa9a9-4c7c-4c62-b466-2421203962a2"
 TypedPolynomials = "afbbf031-7a57-5f58-a1b9-b774a0fad08d"

docs/src/tutorials/Systems and Control/julia_set.jl

@@ -0,0 +1,150 @@
+# # Outer approximation of Julia set
+
+#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/generated/Systems and Control/julia_set.ipynb)
+#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/generated/Systems and Control/julia_set.ipynb)
+# **Adapted from**: Section 7.1.3 of [KHJ14]
+#
+# [KHJ14] Milan Korda, Didier Henrion, and Colin N. Jones.
+# *Convex computation of the maximum controlled invariant set for polynomial control systems*.
+# SIAM Journal on Control and Optimization 52.5 (2014): 2944-2969.
+
+# The Julia map is defined as follows:
+
+using Test #src
+function julia_map(point, c)
+    a, b = point
+    return [a^2 - b^2 + real(c), 2a * b + imag(c)]
+end
+
+
+# The *escape radius" is the radius `r` such that `r^2 ≥ r + abs(c)`.
+# Ouside of the circle of that radius, all points diverge so we know the Julia set belongs to that circle.
+
+escape_radius(c) = (1 + √(1 + 4 * abs(c))) / 2
+
+# To check whether a point is in the Julia set, we can iterate and once the point leaves the circle of escape radius,
+# we consider that it is not in the Julia set, if it stays in the set, we consider that it is in the Julia set.
+# This gives an outer approximation that converges to the Julia set when `m` increases.
+
+using LinearAlgebra
+function in_set(x, c, m=2000)
+    r = escape_radius(c)
+    for i in 1:m
+        if norm(x) > r
+            return false
+        end
+        x = julia_map(x, c)
+    end
+    return true
+end
+
+# To sort of minimize a level set of a polynomial we minimize integral of that polynomial.
+# We borrow the following from https://doi.org/10.1080/00029890.2001.11919774
+
+using SpecialFunctions
+using DynamicPolynomials
+β(α) = (α + 1) / 2
+function circle_integral(mono::AbstractMonomial)
+    if any(isodd, exponents(mono))
+        return 0.0
+    else
+        return 2 * prod(gamma ∘ β, exponents(mono)) / gamma(sum(β, exponents(mono)))
+    end
+end
+function disk_integral(mono::AbstractMonomial, r)
+    d = degree(mono) + nvariables(mono)
+    return circle_integral(mono) * r^d / d
+end
+function disk_integral(p::AbstractPolynomialLike, r)
+    return sum(MultivariatePolynomials.coefficient(t) * disk_integral(monomial(t), r) for t in terms(p))
+end
+
+# The following function implements [KHJ14, (8)].
+
+using SumOfSquares
+function outer_approximation(solver, d::Int, c; α = 1/2)
+    @polyvar x[1:2]
+    model = SOSModel(solver)
+    r = escape_radius(c)
+    S = @set sum(x.^2) <= r^2
+    @variable(model, v, Poly(monomials(x, 0:2d)))
+    @variable(model, w0, SOSPoly(monomials(x, 0:d)))
+    @variable(model, w1, SOSPoly(monomials(x, 0:(d - 1))))
+    @constraint(model, α * v(x => julia_map(x, c)) <= v, domain = S)
+    w = w0 + w1 * (r^2 - sum(x.^2))
+    @constraint(model, w >= v + 1, domain = S)
+    @objective(model, Min, disk_integral(w, r))
+    optimize!(model)
+    @test termination_status(model) in [MOI.OPTIMAL, MOI.ALMOST_OPTIMAL] #src
+    @test primal_status(model) in [MOI.FEASIBLE_POINT, MOI.NEARLY_FEASIBLE_POINT] #src
+    if primal_status(model) == MOI.NO_SOLUTION
+        return
+    end
+    return model
+end
+
+# The following function plots the Julia set with the outer approximation.
+
+using ImplicitPlots
+using Plots
+function julia_plot(poly, c, n=200, m=1000; tol=1e-6, res = 1000)
+    r = escape_radius(c)
+    p = implicit_plot(poly; xlims=(-r, r) .* 1.1, ylims=(-r, r), resolution = res, label="")
+    θ = range(0, stop=2π, length=100)
+    points = Vector{Float64}[]
+    as = range(-r, r, length=n)
+    bs = range(-r, r, length=n)
+    for a in as, b in bs
+        point = [a, b]
+        if in_set(point, c, m)
+            push!(points, point)
+        end
+    end
+    xs = [point[1] for point in points]
+    ys = [point[2] for point in points]
+    scatter!(p, xs, ys, label="", markerstrokewidth=0, markersize=1.5, m=:pixel)
+    return p
+end
+
+# We need to pick an SDP solver, see [here](https://jump.dev/JuMP.jl/v1.8/installation/#Supported-solvers) for a list of the available choices.
+
+import CSDP
+solver = optimizer_with_attributes(CSDP.Optimizer, MOI.Silent() => true)
+
+# Let's start with the value of `c` corresponding to the left image of [KHJ14, Figure 3] and with degree 2.
+
+c = -0.7 + 0.2im
+model = outer_approximation(solver, 2, c)
+solution_summary(model)
+
+# We visualize below:
+
+julia_plot(value(model[:v]), c)
+
+# Let's now look at degree 4.
+
+model = outer_approximation(solver, 4, c)
+solution_summary(model)
+
+# We visualize below:
+
+julia_plot(value(model[:v]), c)
+
+# Let's now use the value of `c` corresponding to the right image of [KHJ14, Figure 3] and with degree 2.
+
+c = -0.9 + 0.2im
+model = outer_approximation(solver, 2, c)
+solution_summary(model)
+
+# We visualize below:
+
+julia_plot(value(model[:v]), c)
+
+# Let's now look at degree 4.
+
+model = outer_approximation(solver, 4, c)
+solution_summary(model)
+
+# We visualize below:
+
+julia_plot(value(model[:v]), c)

examples/run_examples.jl

@@ -1,5 +1,18 @@
 using Test
 
+"""
+    _include_sandbox(filename)
+
+Include the `filename` in a temporary module that acts as a sandbox. (Ensuring
+no constants or functions leak into other files.)
+
+This function was taken from `JuMP/docs/make.jl`.
+"""
+function _include_sandbox(filename)
+    mod = @eval module $(gensym()) end
+    return Base.include(mod, filename)
+end
+
 const _TUTORIAL_DIR = joinpath(@__DIR__, "..", "docs", "src", "tutorials")
 
 @testset "run_examples.jl" begin
@@ -9,7 +22,7 @@ const _TUTORIAL_DIR = joinpath(@__DIR__, "..", "docs", "src", "tutorials")
                 continue
             end
             path = joinpath(_TUTORIAL_DIR, dir, filename)
-            include(path)
+            _include_sandbox(path)
         end
     end
     @testset "Chordal" begin