Work on documentation

README.md

@@ -30,7 +30,7 @@ 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> using IntervalArithmetic.Symbols  # to use `..`
 
 julia> f(x) = sin(x) - 0.1*x^2 + 1
 f (generic function with 1 method)

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

docs/src/roots.md

@@ -2,109 +2,99 @@
 
 ## Methods
 
-Three root-finding methods currently available through the `roots` interface are the following:
+Three root-finding contractors currently available through the `roots` interface are the following:
   - `Newton` (default);
   - `Krawczyk`;
   - `Bisection`
 
 Both the Newton and Krawczyk methods can determine if a root is unique in an interval, at the cost of requiring that the function is differentiable. The bisection method has no such requirement, but can never guarantee the existence or uniqueness of a root.
 
-The method used is given as the third (optional) argument of the `roots` function:
+The method used is given using the `contractor` keyword argument:
 
 ```jl
-julia> roots(log, -2..2, Newton)
-1-element Array{Root{Interval{Float64}},1}:
- Root([0.999996, 1.00001], :unique)
+julia> roots(log, -2..2 ; contractor = Newton)
+1-element Vector{Root{Interval{Float64}}}:
+ Root([0.999999, 1.00001]_com, :unique)
 
-julia> roots(log, -2..2, Krawczyk)
-1-element Array{Root{Interval{Float64}},1}:
- Root([0.999984, 1.00002], :unique)
+julia> roots(log, -2..2 ; contractor = Krawczyk)
+1-element Vector{Root{Interval{Float64}}}:
+ Root([0.999999, 1.00001]_com, :unique)
 
-julia> roots(log, -2..2, Bisection)
-1-element Array{Root{Interval{Float64}},1}:
- Root([0.999454, 1.00039], :unknown)
+julia> roots(log, -2..2 ; contractor = Bisection)
+1-element Vector{Root{Interval{Float64}}}:
+ Root([0.999999, 1.00001]_com, :unknown)
 ```
 
 Note that as shown in the example, the `log` function does not complain about being given an interval going outside of its domain. While this may be surprising, this is the expected behavior and no root will ever be found outside the domain of a function.
 
 ## Explicit derivatives
 
-Newton and Krawczyk methods require the function to be differentiable, but the derivative is usually computed automatically using forward-mode automatic differentiation, provided by the `ForwardDiff.jl` package. It is however possible to provide the derivative explicitly for these methods as the second argument of the `roots` function:
+Newton and Krawczyk methods require the function to be differentiable, but the derivative is usually computed automatically using forward-mode automatic differentiation, provided by the `ForwardDiff.jl` package. It is however possible to provide the derivative explicitly using the `derivative` keyword argument:
 
 ```jl
-julia> roots(log, x -> 1/x, -2..2, Newton)
-1-element Array{Root{Interval{Float64}},1}:
- Root([0.999996, 1.00001], :unique)
+julia> roots(log, -2..2 ; contractor = Newton, derivative = x -> 1/x)
+1-element Vector{Root{Interval{Float64}}}:
+ Root([0.999999, 1.00001]_com_NG, :unique)
 
-julia> roots(log, x -> 1/x, -2..2, Krawczyk)
-1-element Array{Root{Interval{Float64}},1}:
- Root([0.999984, 1.00002], :unique)
+julia> roots(log, -2..2 ; contractor = Krawczyk, derivative = x -> 1/x)
+1-element Vector{Root{Interval{Float64}}}:
+ Root([0.999999, 1.00001]_com_NG, :unique)
 ```
 
 When providing the derivative explicitly, the computation is expected to be slightly faster, but the precision of the result is unlikely to be affected.
 
 ```jl
 julia> using BenchmarkTools
 
-julia> @btime roots(log, x -> 1/x, -2..2, Newton)
-  38.600 μs (371 allocations: 27.01 KiB)
-1-element Array{Root{Interval{Float64}},1}:
- Root([0.999996, 1.00001], :unique)
+julia> @btime roots(log, -2..2 ; derivative = x -> 1/x)
+  4.814 μs (53 allocations: 3.16 KiB)
+1-element Vector{Root{Interval{Float64}}}:
+ Root([0.999999, 1.00001]_com_NG, :unique)
 
-julia> @btime roots(log, -2..2, Newton)
-  51.799 μs (373 allocations: 27.20 KiB)
-1-element Array{Root{Interval{Float64}},1}:
- Root([0.999996, 1.00001], :unique)
+julia> @btime roots(log, -2..2)
+  5.767 μs (53 allocations: 3.16 KiB)
+1-element Vector{Root{Interval{Float64}}}:
+ Root([0.999999, 1.00001]_com, :unique)
 ```
 
 This may be useful in some special cases where `ForwardDiff.jl` is unable to compute the derivative of a function. Examples are complex functions and functions whose interval extension must be manually defined (e.g. special functions like `zeta`).
 
-In dimension greater than one, the function of interest must return a `SVector`, a type provided by the `StaticArrays` package, but otherwise works in the same way as in the 1D case.
+In dimension greater than one, the derivative be given as a function returning the Jacobi matrix:
 
 ```jl
 julia> function f( (x, y) )
-           return SVector(sin(x), cos(y))
+           return [sin(x), cos(y)]
        end
 f (generic function with 1 method)
 
-julia> roots(f, IntervalBox(-3..3, 2))
-2-element Array{Root{IntervalBox{2,Float64}},1}:
- Root([-1.13556e-19, 2.3664e-20] × [1.57079, 1.5708], :unique)
- Root([-7.92188e-19, 3.20973e-19] × [-1.5708, -1.57079], :unique)
-```
-
-When providing the derivative for a multi-dimensional function, this must be given as a function returning the Jacobi matrix of the function as an `SMatrix`:
-
-```jl
 julia> function df( (x, y) )
-           return SMatrix{2, 2}(cos(x), 0, 0, -sin(y))
+           return [cos(x) 0 ; 0 -sin(y)]
        end
-df (generic function with 1 method)
 
-julia> roots(f, df, IntervalBox(-3..3, 2), Newton)
-2-element Array{Root{IntervalBox{2,Float64}},1}:
- Root([-2.35877e-07, 8.22858e-07] × [1.57079, 1.5708], :unique)
- Root([-7.19393e-07, 1.55473e-06] × [-1.5708, -1.57079], :unique)
+julia> roots(f, [-3..3, -3..3] ; derivative = df)
+2-element Vector{Root{Vector{Interval{Float64}}}}:
+ Root(Interval{Float64}[[-1.24409e-21, 1.0588e-22]_com_NG, [-1.5708, -1.57079]_com_NG], :unique)
+ Root(Interval{Float64}[[-1.24409e-21, 1.0588e-22]_com_NG, [1.57079, 1.5708]_com_NG], :unique)
 ```
 
-
 ## Tolerance
 
-An absolute tolerance for the search may be specified as the last argument of the `roots` function, the default being `1e-15`. Currently a method must first be provided in order to be able to choose the tolerance.
+An absolute tolerance for the search may be specified as the `abstol` keyword argument.
+Currently a method must first be provided in order to be able to choose the tolerance.
 
 ```jl
 julia> g(x) = sin(exp(x))
 g (generic function with 1 method)
-
-julia> roots(g, 0..2, Newton)
-2-element Array{Root{Interval{Float64}},1}:
- Root([1.83787, 1.83788], :unique)
- Root([1.14472, 1.14473], :unique)
-
-julia> roots(g, 0..2, Newton, 1e-2)
-2-element Array{Root{Interval{Float64}},1}:
- Root([1.83745, 1.83974], :unique)
- Root([1.14471, 1.14475], :unique)
+        
+julia> roots(g, 0..2)
+2-element Vector{Root{Interval{Float64}}}:
+ Root([1.14472, 1.14474]_com, :unique)
+ Root([1.83787, 1.83788]_com, :unique)
+
+julia> roots(g, 0..2 ; abstol = 1e-1)
+2-element Vector{Root{Interval{Float64}}}:
+ Root([1.14173, 1.15244]_com, :unique)
+ Root([1.78757, 1.84273]_com, :unique)
 ```
 
 A lower tolerance may greatly reduce the computation time, at the cost of an increased number of returned roots having `:unknown` status:
@@ -113,32 +103,32 @@ A lower tolerance may greatly reduce the computation time, at the cost of an inc
 julia> h(x) = cos(x) * sin(1 / x)
 h (generic function with 1 method)
 
-julia> @btime roots(h, 0.05..1, Newton)
-  1.316 ms (9676 allocations: 202.21 KiB)
-6-element Array{Root{Interval{Float64}},1}:
- Root([0.106103, 0.106104], :unique)
- Root([0.318309, 0.31831], :unique)
- Root([0.0795774, 0.0795775], :unique)
- Root([0.0636619, 0.063662], :unique)
- Root([0.0530516, 0.0530517], :unique)
- Root([0.159154, 0.159155], :unique)             
-
-julia> @btime roots(h, 0.05..1, Newton, 1e-2)
-   475.500 μs (4171 allocations: 94.00 KiB)
- 6-element Array{Root{Interval{Float64}},1}:
-  Root([0.317179, 0.319299], :unique)
-  Root([0.157209, 0.165989], :unknown)
-  Root([0.104739, 0.107542], :unknown)
-  Root([0.0515382, 0.0531049], :unique)
-  Root([0.0785458, 0.0797755], :unknown)
-  Root([0.0570253, 0.0641614], :unknown)
-
-julia> @btime roots(h, 0.05..1, Newton, 1e-1)
-  151.101 μs (1382 allocations: 32.86 KiB)
-3-element Array{Root{Interval{Float64}},1}:
- Root([0.107541, 0.165989], :unknown)
- Root([0.283803, 0.3555], :unknown)
- Root([0.0499999, 0.107542], :unknown)     
+julia> @btime roots(h, 0.05..1)
+  79.500 μs (301 allocations: 13.97 KiB)
+6-element Vector{Root{Interval{Float64}}}:
+ Root([0.0530516, 0.0530517]_com_NG, :unique)
+ Root([0.0636619, 0.063662]_com_NG, :unique)
+ Root([0.0795774, 0.0795775]_com_NG, :unique)
+ Root([0.106103, 0.106104]_com_NG, :unique)
+ Root([0.159154, 0.159156]_com_NG, :unique)
+ Root([0.318309, 0.31831]_com_NG, :unique)
+
+julia> @btime roots(h, 0.05..1 ; abstol = 1e-2)
+  48.500 μs (265 allocations: 11.61 KiB)
+6-element Vector{Root{Interval{Float64}}}:
+ Root([0.0514445, 0.0531087]_com_NG, :unique)
+ Root([0.0570253, 0.0641615]_com, :unknown)
+ Root([0.0785458, 0.0797797]_com_NG, :unknown)
+ Root([0.104754, 0.107542]_com_NG, :unknown)
+ Root([0.157236, 0.165989]_com_NG, :unknown)
+ Root([0.31716, 0.319318]_com_NG, :unique)
+
+julia> @btime roots(h, 0.05..1 ; abstol = 1e-1)
+  17.300 μs (114 allocations: 5.16 KiB)
+3-element Vector{Root{Interval{Float64}}}:
+ Root([0.04999999, 0.107542]_com, :unknown)
+ Root([0.107541, 0.165989]_com, :unknown)
+ Root([0.283803, 0.356099]_com_NG, :unknown)
 ```
 
 The last example shows a case where the tolerance was too large to be able to isolate the roots in distinct regions.