Compatibility with IntervalArithmetic v0.22

.github/workflows/CI.yml

@@ -13,7 +13,7 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - '1'
+          - '1.10'
           - 'nightly'
         os:
           - ubuntu-latest

NEWS.md

@@ -1,5 +1,17 @@
 # Updates to `IntervalRootFinding.jl`
 
+## v0.6
+
+Compatibility with IntervalArithmetic v0.22, which fundamentally changes the package.
+- Decorated intervals must be used (which is the default `Interval` in the latest releases of IntervalArithmetic)
+- The signature of the `roots` function changed to `roots(f::Function, X::Union{Interval, AbstractVector} ; kwargs...)`, with the following consequences
+    - Manual derivatives and contractors must now always be passed as keyword arguments. This greatly simplify the internal logic of the `roots` function.
+    - No more `IntervalBox`. Multidimensional problem are specified by returning a vector of intervals, and giving a vector of intervals as initial search region.
+    - Normal vectors of intervals are supported. They are significantly (3x) time slower for a simple 2D problem than SVector, but more convenient.
+- Features of the packages outside the `roots` function are unsupported for the time being (e.g. `Slopes` and quadratic equations). If you need them, please open an issue.
+- `roots` works with `@exact`.
+
+
 ## v0.2
 
 The package has been completely rewritten for v0.2.

Project.toml

@@ -1,23 +1,21 @@
 name = "IntervalRootFinding"
 uuid = "d2bf35a9-74e0-55ec-b149-d360ff49b807"
-version = "0.5.11"
+version = "0.6"
 
 [deps]
 BranchAndPrune = "d3bc4f2e-91e6-11e9-365e-cd067da536ce"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
-Polynomials = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
 [compat]
 BranchAndPrune = "0.2.1"
 ForwardDiff = "0.10"
-IntervalArithmetic = "0.15, 0.16, 0.17, 0.18, 0.19, 0.20"
-Polynomials = "0.5, 0.6, 0.7, 0.8, 1, 2, 3"
+IntervalArithmetic = "0.22.13"
 Reexport = "1"
-StaticArrays = "0.11, 0.12, 1.0"
+StaticArrays = "1"
 julia = "1"
 
 [extras]

README.md

@@ -30,6 +30,8 @@ The basic function is `roots`. Given a standard Julia function and an interval,
 ```julia
 julia> using IntervalArithmetic, IntervalRootFinding
 
+julia> using IntervalArithmetic.Symbols  # to use `..`
+
 julia> f(x) = sin(x) - 0.1*x^2 + 1
 f (generic function with 1 method)
 

benchmark/benchmarks.jl

@@ -49,8 +49,8 @@ sizes = (2, 5, 10)
 
 for n in sizes
     s = S["n = $n"] = BenchmarkGroup()
-    A = Interval.(randn(n, n))
-    b = Interval.(randn(n))
+    A = interval.(randn(n, n))
+    b = interval.(randn(n))
 
     s["Gauss seidel"] = @benchmarkable gauss_seidel_interval($A, $b)
     s["Gauss seidel contractor"] = @benchmarkable gauss_seidel_contractor($A, $b)
@@ -59,7 +59,7 @@ end
 
 
 S = SUITE["Dietmar-Ratz"] = BenchmarkGroup()
-X = Interval(0.75, 1.75)
+X = interval(0.75, 1.75)
 
 for (k, dr) in enumerate(dr_functions)
     s = S["Dietmar-Ratz $k"] = BenchmarkGroup()

docs/make.jl

@@ -10,6 +10,7 @@ makedocs(
     pages = Any[
         "Home" => "index.md",
         "`roots` interface" => "roots.md",
+        "Decorations and guarantees" => "decorations.md",
         "Internals" => "internals.md",
         "Bibliography" => "biblio.md",
         "API" => "api.md"

docs/src/decorations.md

@@ -0,0 +1,78 @@
+# Decorations and guarantees
+
+When you define a simple function, you will notice that the roots
+usually come out with the flag `com` and `NG`.
+
+```jl
+julia> roots(x -> x^2 - 0.1, interval(-5, 5))
+2-element Vector{Root{Interval{Float64}}}:
+ Root([-0.316228, -0.316227]_com_NG, :unique)
+ Root([0.316227, 0.316228]_com_NG, :unique)
+```
+
+In this case, `com` is the *decoration* of the interval,
+and `NG` its *guarantee*.
+See the documentation of `IntervalArithmetic.jl` for a detailed description,
+here we will go on what they mean in the context of finding roots
+of a function.
+
+## Decorations
+
+The decoration tells us whether all operations that affected the interval
+were well-behaved.
+The two crucial ones for us are `def` (the operations were not continuous)
+and `trv` (something horribly wrong happened).
+
+In both cases, it means that the root can not be trusted,
+as the hypothesis of the theory we use are not fulfilled.
+
+## Guarantee
+
+The `NG` flags means that non-interval have been mixed with intervals.
+Since we can not track the source of the non-interval numbers,
+we can not guarantee that they are correct.
+For example, `0.1` is famously not parsed as `0.1`,
+as `0.1` can not be represented exactly as a binary number
+(just like `1/3` can not be represented exactly as a decimal number).
+
+```jl
+julia> big(0.1)
+0.1000000000000000055511151231257827021181583404541015625
+```
+
+If you want the number that you are inputting to be trusted **as is**,
+you can use the `@exact` macro from `IntervalArithmetic.jl`,
+and the `NG` flag will be avoided in most cases.
+
+```jl
+julia> @exact f(x) = x^2 - 0.1
+f (generic function with 1 method)
+
+julia> roots(f, interval(-5, 5))
+2-element Vector{Root{Interval{Float64}}}:
+ Root([-0.316228, -0.316227]_com, :unique)
+ Root([0.316227, 0.316228]_com, :unique)
+```
+
+The `NG` flag can still appear if other computations,
+from a library using non-interval "magic" numbers for example,
+thus indicating that some non-trusted numbers have been used in the computation.
+
+Moreover, the macro work in such a way that you can still use the defined
+function with numbers and get floating point results.
+
+```jl
+julia> f(0.2)
+-0.06
+```
+## Trust
+
+We are doing our best to really give validated and guaranteed results.
+However, in the case that you may make your house explode based on a result
+returned by this package,
+we would like to remind you that you should not trust the package beyond
+the promises of the license:
+
+> THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED
+
+All bug reports and pull requests to fix them are however more than welcome.

docs/src/index.md

@@ -10,19 +10,21 @@ To do so, it uses methods from interval analysis, using interval arithmetic from
 
 ## Basic 1D example
 
-To begin, we need a standard Julia function and an interval in which to search roots of that function. Intervals use the `Interval` type provided by the `IntervalArithmetic.jl` package and are generally constructed using the `..` syntax, `a..b` representing the closed interval $[a, b]$.
+To begin, we need a standard Julia function and an interval in which to search roots of that function. Intervals use the `Interval` type provided by the `IntervalArithmetic.jl` package
+Intervals are generally constructed using the `..` syntax (from the `IntervalArithmetic.Symbols` submodule),
+`a..b` representing the closed interval $[a, b]$.
 
 When provided with this information, the `roots` function will return a vector of all roots of the function in the given interval.
 
 Example:
 
 ```jl
-julia> using IntervalArithmetic, IntervalRootFinding
+julia> using IntervalArithmetic, IntervalArithmetic.Symbols, IntervalRootFinding
 
 julia> rts = roots(x -> x^2 - 2x, 0..10)
-2-element Array{Root{Interval{Float64}},1}:
- Root([1.99999, 2.00001], :unique)
- Root([0, 4.4724e-16], :unknown)
+2-element Vector{Root{Interval{Float64}}}:
+ Root([0.0, 3.73849e-08]_com_NG, :unknown)
+ Root([1.999999, 2.00001]_com_NG, :unique)
 ```
 
 The roots are returned as `Root` objects, containing an interval and the status of that interval, represented as a `Symbol`. There are two possible types of root status, as shown in the example:
@@ -42,68 +44,89 @@ julia> g(x) = (x^2 - 2)^2 * (x^2 - 3)
 g (generic function with 1 method)
 
 julia> roots(g, -10..10)
-4-element Array{IntervalRootFinding.Root{IntervalArithmetic.Interval{Float64}},1}:
- Root([1.73205, 1.73206], :unique)
- Root([1.41418, 1.4148], :unknown)
- Root([-1.4148, -1.41418], :unknown)
- Root([-1.73206, -1.73205], :unique)
+4-element Vector{Root{Interval{Float64}}}:
+ Root([-1.73206, -1.73205]_com_NG, :unique)
+ Root([-1.41422, -1.41421]_com, :unknown)  
+ Root([1.41421, 1.41422]_com, :unknown)    
+ Root([1.73205, 1.73206]_com_NG, :unique)  
 ```
 
 Here we see that the two double roots are reported as being possible roots without guarantee and the simple roots have been proved to be unique.
 
 
 ## Basic multi-dimensional example
 
-For dimensions $n > 1$, the function passed to `roots` must currently return an `SVector` from the `StaticArrays.jl` package.
+For dimensions $n > 1$, the function passed to `roots` must take an array as
+argument and return an array.
+The initial search region is an array of interval.
 
 Here we give a 3D example:
 
 ```jl
 julia> function g( (x1, x2, x3) )
-           return SVector(x1^2 + x2^2 + x3^2 - 1,
-                          x1^2 + x3^2 - 0.25,
-                          x1^2 + x2^2 - 4x3
-                         )
+          return [
+              x1^2 + x2^2 + x3^2 - 1,
+              x1^2 + x3^2 - 0.25,
+              x1^2 + x2^2 - 4x3
+          ]
        end
 g (generic function with 1 method)
 
 julia> X = -5..5
-[-5, 5]
-
-julia> rts = roots(g, X × X × X)
-4-element Array{Root{IntervalBox{3,Float64}},1}:
- Root([0.440762, 0.440763] × [0.866025, 0.866026] × [0.236067, 0.236068], :unique)
- Root([0.440762, 0.440763] × [-0.866026, -0.866025] × [0.236067, 0.236068], :unique)
- Root([-0.440763, -0.440762] × [0.866025, 0.866026] × [0.236067, 0.236068], :unique)
- Root([-0.440763, -0.440762] × [-0.866026, -0.866025] × [0.236067, 0.236068], :unique)
+[-5.0, 5.0]_com
+
+julia> rts = roots(g, [X, X, X])
+4-element Vector{Root{Vector{Interval{Float64}}}}:
+ Root(Interval{Float64}[[-0.440763, -0.440762]_com_NG, [-0.866026, -0.866025]_com_NG, [0.236067, 0.236069]_com_NG], :unique)
+ Root(Interval{Float64}[[-0.440763, -0.440762]_com_NG, [0.866025, 0.866026]_com_NG, [0.236067, 0.236069]_com_NG], :unique)
+ Root(Interval{Float64}[[0.440762, 0.440763]_com_NG, [-0.866026, -0.866025]_com_NG, [0.236067, 0.236069]_com_NG], :unique)
+ Root(Interval{Float64}[[0.440762, 0.440763]_com_NG, [0.866025, 0.866026]_com_NG, [0.236067, 0.236069]_com_NG], :unique)
 ```
 
-Thus the system admits four unique roots in the box $[-5, 5]^3$. We have used the unicode character `×` (typed as `\times<tab>`) to compose several intervals into a multidimensional box.
+Thus, the system admits four unique roots in the box $[-5, 5]^3$.
+
+Moreover, the package is compatible with `StaticArrays.jl`.
+Usage of static arrays is recommended to increase performance.
+```jl
+julia> using StaticArrays
+
+julia> h((x, y)) = SVector(x^2 - 4, y^2 - 16)
+h (generic function with 1 method)
+
+julia> roots(h, SVector(X, X))
+4-element Vector{Root{SVector{2, Interval{Float64}}}}:
+ Root(Interval{Float64}[[-2.00001, -1.999999]_com_NG, [-4.00001, -3.999999]_com_NG], :unique)
+ Root(Interval{Float64}[[-2.00001, -1.999999]_com_NG, [3.999999, 4.00001]_com_NG], :unique)
+ Root(Interval{Float64}[[1.999999, 2.00001]_com_NG, [-4.00001, -3.999999]_com_NG], :unique)
+ Root(Interval{Float64}[[1.999999, 2.00001]_com_NG, [3.999999, 4.00001]_com_NG], :unique)
+```
 
 ## Stationary points
 
 Stationary points of a function $f:\mathbb{R}^n \to \mathbb{R}$ may be found as zeros of the gradient of $f$.
-The package exports the `∇` operator to calculate gradients using `ForwardDiff.jl`:
+The gradient can be computed using `ForwardDiff.jl`:
 
 ```jl
+julia> using ForwardDiff: gradient
+
 julia> f( (x, y) ) = sin(x) * sin(y)
 f (generic function with 1 method)
 
-julia> ∇f = ∇(f)  # gradient operator from the package
+julia> ∇f(x) = gradient(f, x)  # gradient operator from the package
 (::#53) (generic function with 1 method)
 
 julia> rts = roots(∇f, IntervalBox(-5..6, 2), Newton, 1e-5)
-25-element Array{IntervalRootFinding.Root{IntervalArithmetic.IntervalBox{2,Float64}},1}:
- Root([4.71238, 4.71239] × [4.71238, 4.71239], :unique)
- Root([4.71238, 4.71239] × [1.57079, 1.5708], :unique)
+25-element Vector{Root{Vector{Interval{Float64}}}}:
+ Root(Interval{Float64}[[-4.7124, -4.71238]_com, [-4.7124, -4.71238]_com], :unique)
+ Root(Interval{Float64}[[-3.1416, -3.14159]_com, [-3.1416, -3.14159]_com], :unique)
  ⋮
- [output snipped for brevity]
+ [output skipped for brevity]
 ```
 
 Now let's find the midpoints and plot them:
 
 ```jl
-midpoints = mid.(interval.(rts))
+midpoints = [mid.(root_region(rt)) for rt in rts]
 
 xs = first.(midpoints)
 ys = last.(midpoints)
@@ -114,6 +137,4 @@ surface(-5:0.1:6, -6:0.1:6, (x,y) -> f([x, y]))
 scatter!(xs, ys, f.(midpoints))
 ```
 
-The result is the following:
-
 ![stationary points](stationary_points.png)