API Reference

bisect(f, grid; threaded=false, monotonic=false, iterations=5, verbose=false)

Perform the multidimensional bisection algorithm of function f on an initial grid of evaluation points.

Arguments

Positional Arguments

• f: function to be evaluated
• grid: a tuple of ranges describing the initial evaluation grid

Keyword Arguments

• threaded: set to true for multithreaded evaluation
• monotonic: set to true if f is known to be monotonic in each dimension
• iterations: number of times to iterate the algorithm
• verbose: set to true to print the number of new function evaluations at each iteration

Examples

julia> grid = (0.0:1.0, 0.0:1.0);

julia> bisect(x -> 1 - sum(abs2, x), grid)
BisectionGrid{Float64, 2}
       Domain: (0.0:0.0625:1.0, 0.0:0.0625:1.0)
  Grid points: 289
  Evaluations: 112
        Edges: 32


julia> bisect(x -> x[2] - x[1], grid)
BisectionGrid{Float64, 2}
       Domain: (0.0:0.0625:1.0, 0.0:0.0625:1.0)
  Grid points: 289
  Evaluations: 112
        Edges: 32


julia> bisect(x -> x[2] - x[1], grid; monotonic=true)
BisectionGrid{Float64, 2}
       Domain: (0.0:0.0625:1.0, 0.0:0.0625:1.0)
  Grid points: 289
  Evaluations: 67
        Edges: 32

Extended help

Algorithm

The algorithm begins with an initial evaluation grid of $N$-dimensional hypercubes (lines, squares, cubes, etc.) and proceeds by splitting it in "half" at each stage, but only in regions where the function is known to change sign. After the function is evaluated, each hypercube is split in half (a line into two lines, a square into four squares, a cube into eight cubes, etc.) and the function sign at each vertex is checked. If not all vertices share the same sign, the function must change sign somewhere within the hypercube, so the vertices of the subcubes (e.g. the midpoint of a line) will be marked for evaluation in the next stage; if the vertices all share the same sign, the function is presumed to have the same sign everywhere throughout the hypercube, so the sign is inferred at all interior points and the function is not evaluated in the hypercube again.

By passing the keyword monotonic=true, the algorithm further assumes that the function cannot change sign in any one direction between two points with the same sign. This means that the sign of the function should be checked at the vertices of each edge of the hypercube; if the signs are the same, the function is assumed to have the same sign at the midpoint along the edge and will not be evaluated.

Multithreading

Set the keyword argument threaded=true to make evaluation of the function at each iteration multithreaded. For functions that are very cheap to evaluate, the overhead from setting up threads will actually make each iteration longer. However, this effect disappears quickly.

Example

julia> ffast(x) = 1 - sum(abs2, x);

julia> grid = (0.0:1.0, 0.0:1.0);

julia> @btime bisect(ffast, $grid; threaded=false);
  49.720 μs (40 allocations: 3.59 KiB)

julia> @btime bisect(ffast, $grid; threaded=true);
  146.173 μs (163 allocations: 18.23 KiB)

julia> fslow(x) = (sleep(0.1); ffast(x))
fslow (generic function with 1 method)

julia> @btime bisect(fslow, $grid; threaded=false);
  11.397 s (605 allocations: 20.83 KiB)

julia> @btime bisect(fslow, $grid; threaded=true);
  3.170 s (736 allocations: 35.45 KiB)

Monotonicity

If the function f is known to be monotonic, then the number of evaluations can be reduced by inferring the behavior of f along edges where both vertices have the same sign.

Example

julia> bisect(x -> 1 - sum(abs2, x), (0.0:1.0, 0.0:1.0); monotonic=false)
BisectionGrid{Float64, 2}
       Domain: (0.0:0.0625:1.0, 0.0:0.0625:1.0)
  Grid points: 289
  Evaluations: 112
        Edges: 32


julia> bisect(x -> 1 - sum(abs2, x), (0.0:1.0, 0.0:1.0); monotonic=true)
BisectionGrid{Float64, 2}
       Domain: (0.0:0.0625:1.0, 0.0:0.0625:1.0)
  Grid points: 289
  Evaluations: 80
        Edges: 32

Iterations

For finer grids, repeat the algorithm for more iterations.

Example

julia> bisect(x -> 1 - sum(abs2, x), (0.0:1.0, 0.0:1.0))
BisectionGrid{Float64, 2}
       Domain: (0.0:0.0625:1.0, 0.0:0.0625:1.0)
  Grid points: 289
  Evaluations: 112
        Edges: 32


julia> bisect(x -> 1 - sum(abs2, x), (0.0:1.0, 0.0:1.0); iterations=3)
BisectionGrid{Float64, 2}
       Domain: (0.0:0.25:1.0, 0.0:0.25:1.0)
  Grid points: 25
  Evaluations: 22
        Edges: 8


julia> bisect(x -> 1 - sum(abs2, x), (0.0:1.0, 0.0:1.0); iterations=7)
BisectionGrid{Float64, 2}
       Domain: (0.0:0.015625:1.0, 0.0:0.015625:1.0)
  Grid points: 4225
  Evaluations: 490
        Edges: 128

Logging

For longer-running bisections, it may be useful to report the number of new evaluations as a kind of thread-safe progress meter.

Example

julia> bisect(x -> 1 - sum(abs2, x), (0.0:1.0, 0.0:1.0); verbose=true)
[ Info: Iteration 1: 4 gridpoints (4 evaluations)
[ Info: Iteration 2: 9 gridpoints (5 evaluations)
[ Info: Iteration 3: 25 gridpoints (13 evaluations)
[ Info: Iteration 4: 81 gridpoints (29 evaluations)
[ Info: Iteration 5: 289 gridpoints (61 evaluations)
BisectionGrid{Float64, 2}
       Domain: (0.0:0.0625:1.0, 0.0:0.0625:1.0)
  Grid points: 289
  Evaluations: 112
        Edges: 32


julia> bisect(x -> 1 - sum(abs2, x), (0.0:1.0, 0.0:1.0); verbose=true, iterations=7)
[ Info: Iteration 1: 4 gridpoints (4 evaluations)
[ Info: Iteration 2: 9 gridpoints (5 evaluations)
[ Info: Iteration 3: 25 gridpoints (13 evaluations)
[ Info: Iteration 4: 81 gridpoints (29 evaluations)
[ Info: Iteration 5: 289 gridpoints (61 evaluations)
[ Info: Iteration 6: 1089 gridpoints (125 evaluations)
[ Info: Iteration 7: 4225 gridpoints (253 evaluations)
BisectionGrid{Float64, 2}
       Domain: (0.0:0.015625:1.0, 0.0:0.015625:1.0)
  Grid points: 4225
  Evaluations: 490
        Edges: 128

corners(CI::CartesianIndices)

Construct an iterator over the corners of the hypercube defined by CI.

Examples

julia> CI = CartesianIndices((1:3, 1:4))
CartesianIndices((1:3, 1:4))

julia> collect(CI)
3×4 Matrix{CartesianIndex{2}}:
 CartesianIndex(1, 1)  CartesianIndex(1, 2)  …  CartesianIndex(1, 4)
 CartesianIndex(2, 1)  CartesianIndex(2, 2)     CartesianIndex(2, 4)
 CartesianIndex(3, 1)  CartesianIndex(3, 2)     CartesianIndex(3, 4)

julia> corners(CI)
CartesianIndices((1:2:3, 1:3:4))

julia> collect(corners(CI))
2×2 Matrix{CartesianIndex{2}}:
 CartesianIndex(1, 1)  CartesianIndex(1, 4)
 CartesianIndex(3, 1)  CartesianIndex(3, 4)

domainindex(X, I::CartesianIndex)

Like getindex, except X is a tuple of ranges instead of an array.

Examples

julia> X = (1.0:5.0, 0.5:2.5);

julia> domainindex(X, CartesianIndex(1, 1))
(1.0, 0.5)

julia> domainindex(X, CartesianIndex(2, 1))
(2.0, 0.5)

julia> domainindex(X, CartesianIndex(1, 2))
(1.0, 1.5)

julia> domainindex(X, last(CartesianIndices(eachindex.(X))))
(5.0, 2.5)

edgebounds(edge)

Return the bounds of the edge in its varying dimension.

Examples

julia> edgebounds(((0.5, 1.0), (1.0, 1.0)))
(0.5, 1.0)

julia> edgebounds(((0.5, 1.0), (0.5, 1.5)))
(1.0, 1.5)

julia> edgebounds(((1.0, 1.0, 1.0), (1.0, 1.0, 2.0)))
(1.0, 2.0)

edgedim(edge)

Return the dimension of the edge (i.e. the only direction in which the vertices vary).

Examples

julia> edgedim(((0.5, 1.0), (1.0, 1.0)))
1

julia> edgedim(((0.5, 1.0), (0.5, 1.5)))
2

julia> edgedim(((1.0, 1.0, 1.0), (1.0, 1.0, 2.0)))
3

edgeindices(BG::BisectionGrid{T, N}) where {T, N}

Find the CartesianIndices of all edges in the current grid domain.

edgeroot(f, edge, M; kwargs...)

Find the root along edge with a call to Roots.find_zero. Both M and kwargs are passed directly to find_zero.

Examples

julia> f(x) = 1 - sum(abs2, x); # unit circle

julia> BG = bisect(f, (0.0:1.0, 0.0:1.0); iterations=3);

julia> E = edges(BG)
8-element Vector{Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}:
 ((0.75, 0.0), (1.0, 0.0))
 ((0.75, 0.25), (1.0, 0.25))
 ((0.75, 0.5), (1.0, 0.5))
 ((0.75, 0.5), (0.75, 0.75))
 ((0.0, 0.75), (0.0, 1.0))
 ((0.25, 0.75), (0.25, 1.0))
 ((0.5, 0.75), (0.75, 0.75))
 ((0.5, 0.75), (0.5, 1.0))

julia> tracks = Roots.Tracks();

julia> root = edgeroot(f, E[3], Roots.ITP(); tracks)
(0.8660254037844387, 0.5)

julia> f(root)
0.0

julia> tracks
Results of univariate zero finding:

* Converged to: 0.8660254037844387
* Algorithm: Roots.ITP{Float64, Int64}(0.2, 2, 1)
* iterations: 7
* function evaluations ≈ 9
* stopped as f(x_n) = 0

Trace:
(a₀, b₀) = ( 0.75, 1 )
(a₁, b₁) = ( 0.75, 0.86964285714285705 )
(a₂, b₂) = ( 0.86290337975046705, 0.86964285714285705 )
(a₃, b₃) = ( 0.86290337975046705, 0.8660279692952968 )
(a₄, b₄) = ( 0.86602344653979335, 0.8660279692952968 )
(a₅, b₅) = ( 0.86602344653979335, 0.86602540378563075 )
(a₆, b₆) = ( 0.86602540378367243, 0.86602540378563075 )
(a₇, b₇) = ( 0.86602540378443871, 0.86602540378563075 )

edges(BG::BisectionGrid)

Find all hypercube edges where the sign of the function changes. Returns a Vector of 2-tuples of points in the domain.

edgetuple(edge)

Return a function converting a value along the edge dimension to a point.

Examples

julia> tup = edgetuple(((0.5, 1.0), (1.0, 1.0)));

julia> tup(0.75)
(0.75, 1.0)

julia> tup(0.9)
(0.9, 1.0)

julia> edgetuple(((0.5, 1.0), (0.5, 1.5)))(1.2)
(0.5, 1.2)

edgevalues!(f, BG::BisectionGrid{T, N}; threaded=false, verbose=false) where {T, N}

Evaluate f at all edge vertices where it has not been evaluated.

expandzeros(A)

Expand the array A, inserting a zero between each element.

Examples

julia> A = reshape(1:4, 2, 2)
2×2 reshape(::UnitRange{Int64}, 2, 2) with eltype Int64:
 1  3
 2  4

julia> expandzeros(A)
3×3 Matrix{Int64}:
 1  0  3
 0  0  0
 2  0  4

julia> B = reshape(1:6, 3, 2)
3×2 reshape(::UnitRange{Int64}, 3, 2) with eltype Int64:
 1  4
 2  5
 3  6

julia> expandzeros(B)
5×3 Matrix{Int64}:
 1  0  4
 0  0  0
 2  0  5
 0  0  0
 3  0  6

findhypercubes(A; dist=2)

Find the starting vertex of all hypercubes in the array A.

Examples

julia> A = ones(Bool, 5, 5)
5×5 Matrix{Bool}:
 1  1  1  1  1
 1  1  1  1  1
 1  1  1  1  1
 1  1  1  1  1
 1  1  1  1  1

julia> cubes = findhypercubes(A)
CartesianIndices((1:2:3, 1:2:3))

julia> collect(cubes)
2×2 Matrix{CartesianIndex{2}}:
 CartesianIndex(1, 1)  CartesianIndex(1, 3)
 CartesianIndex(3, 1)  CartesianIndex(3, 3)

julia> cubes3 = findhypercubes(ones(Bool, 5, 5, 5))
CartesianIndices((1:2:3, 1:2:3, 1:2:3))

julia> collect(cubes3)
2×2×2 Array{CartesianIndex{3}, 3}:
[:, :, 1] =
 CartesianIndex(1, 1, 1)  CartesianIndex(1, 3, 1)
 CartesianIndex(3, 1, 1)  CartesianIndex(3, 3, 1)

[:, :, 2] =
 CartesianIndex(1, 1, 3)  CartesianIndex(1, 3, 3)
 CartesianIndex(3, 1, 3)  CartesianIndex(3, 3, 3)

forwardinds(d)

Get CartesianIndex steps forward in each of d dimensions.

Examples

julia> forwardinds(1)
(CartesianIndex(1,),)

julia> forwardinds(2)
(CartesianIndex(1, 0), CartesianIndex(0, 1))

julia> forwardinds(3)
(CartesianIndex(1, 0, 0), CartesianIndex(0, 1, 0), CartesianIndex(0, 0, 1))

interpolate(rootfinder, BG::BisectionGrid{T, N}; threaded=false) where {T, N}

Compute the roots along all bracketing edges of BG.

The function rootfinder can be any method that receives an edge and returns a root, such as edgeroot, linearroot, or marchingsquares.

Examples

julia> f(x) = 1 - sum(abs2, x); # unit circle

julia> BG = bisect(f, (0.0:1.0, 0.0:1.0); iterations=3);

julia> interpolate(marchingsquares, BG)
8-element Vector{Tuple{Float64, Float64}}:
 (0.875, 0.0)
 (0.875, 0.25)
 (0.875, 0.5)
 (0.75, 0.625)
 (0.0, 0.875)
 (0.25, 0.875)
 (0.625, 0.75)
 (0.5, 0.875)

julia> interpolate(edge -> linearroot(f, edge), BG)
8-element Vector{Tuple{Float64, Float64}}:
 (1.0, 0.0)
 (0.9642857142857143, 0.25)
 (0.8571428571428571, 0.5)
 (0.75, 0.65)
 (0.0, 1.0)
 (0.25, 0.9642857142857143)
 (0.65, 0.75)
 (0.5, 0.8571428571428571)

julia> interpolate(edge -> edgeroot(f, edge, Roots.ITP()), BG)
8-element Vector{Tuple{Float64, Float64}}:
 (1.0, 0.0)
 (0.9682458365518543, 0.25)
 (0.8660254037844387, 0.5)
 (0.75, 0.6614378277661477)
 (0.0, 1.0)
 (0.25, 0.9682458365518543)
 (0.6614378277661477, 0.75)
 (0.5, 0.8660254037844387)

interpolate(BG::BisectionGrid{T, N}) where {T, N}

Compute the roots along all bracketing edges of BG with a linear approximation using precomputed function evaluations. Equivalent to interpolate(linearroot(f), BG), but does not require re-evaluating the function at each edge.

Examples

```jldoctest julia> f(x) = 1 - sum(abs2, x); # unit circle

julia> BG = bisect(f, (0.0:1.0, 0.0:1.0); iterations=3);

julia> interpolate(linearroot(f), BG) 8-element Vector{Tuple{Float64, Float64}}: (1.0, 0.0) (0.9642857142857143, 0.25) (0.8571428571428571, 0.5) (0.75, 0.65) (0.0, 1.0) (0.25, 0.9642857142857143) (0.65, 0.75) (0.5, 0.8571428571428571)

julia> interpolate(BG) 8-element Vector{Tuple{Float64, Float64}}: (1.0, 0.0) (0.9642857142857143, 0.25) (0.8571428571428571, 0.5) (0.75, 0.65) (0.0, 1.0) (0.25, 0.9642857142857143) (0.65, 0.75) (0.5, 0.8571428571428571)

linearroot(f, edge)

Construct a linear interpolation of f through edge and return its root. Requires two function evaluations.

Examples

julia> f(x) = 1 - sum(abs2, x); # unit circle

julia> BG = bisect(f, (0.0:1.0, 0.0:1.0); iterations=3);

julia> E = edges(BG)
8-element Vector{Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}:
 ((0.75, 0.0), (1.0, 0.0))
 ((0.75, 0.25), (1.0, 0.25))
 ((0.75, 0.5), (1.0, 0.5))
 ((0.75, 0.5), (0.75, 0.75))
 ((0.0, 0.75), (0.0, 1.0))
 ((0.25, 0.75), (0.25, 1.0))
 ((0.5, 0.75), (0.75, 0.75))
 ((0.5, 0.75), (0.5, 1.0))

julia> linearroot.(f, E)
8-element Vector{Tuple{Float64, Float64}}:
 (1.0, 0.0)
 (0.9642857142857143, 0.25)
 (0.8571428571428571, 0.5)
 (0.75, 0.65)
 (0.0, 1.0)
 (0.25, 0.9642857142857143)
 (0.65, 0.75)
 (0.5, 0.8571428571428571)

marchingsquares(edge)

Return the midpoint of the edge, in a manner similar to the marching squares algorithm. Does not require any function evaluations.

Examples

julia> f(x) = 1 - sum(abs2, x); # unit circle

julia> BG = bisect(f, (0.0:1.0, 0.0:1.0); iterations=3);


julia> E = edges(BG)
8-element Vector{Tuple{Tuple{Float64, Float64}, Tuple{Float64, Float64}}}:
 ((0.75, 0.0), (1.0, 0.0))
 ((0.75, 0.25), (1.0, 0.25))
 ((0.75, 0.5), (1.0, 0.5))
 ((0.75, 0.5), (0.75, 0.75))
 ((0.0, 0.75), (0.0, 1.0))
 ((0.25, 0.75), (0.25, 1.0))
 ((0.5, 0.75), (0.75, 0.75))
 ((0.5, 0.75), (0.5, 1.0))

julia> marchingsquares.(E)
8-element Vector{Tuple{Float64, Float64}}:
 (0.875, 0.0)
 (0.875, 0.25)
 (0.875, 0.5)
 (0.75, 0.625)
 (0.0, 0.875)
 (0.25, 0.875)
 (0.625, 0.75)
 (0.5, 0.875)

splitsign(BG::BisectionGrid)

Split the evaluated points of the bisection grid into a vector of positive points and a vector of negative points.

Examples

julia> BG = bisect(x -> 1 - sum(abs2, x), (0.0:1.0, 0.0:1.0); iterations=2)
BisectionGrid{Float64, 2}
       Domain: (0.0:0.5:1.0, 0.0:0.5:1.0)
  Grid points: 9
  Evaluations: 9
        Edges: 4

julia> posx, negx = splitsign(BG);


julia> posx
4-element Vector{Tuple{Float64, Float64}}:
 (0.0, 0.0)
 (0.5, 0.0)
 (0.0, 0.5)
 (0.5, 0.5)

julia> negx
5-element Vector{Tuple{Float64, Float64}}:
 (1.0, 0.0)
 (1.0, 0.5)
 (0.0, 1.0)
 (0.5, 1.0)
 (1.0, 1.0)