Docstrings

BernsteinPolynomial(c::AbstractArray, lc, hc)

Create a multidimensional Bernstein polynomial with coefficients c defined on the hyperrectangle [lc[1], hc[1]] × … × [lc[N], hc[N]].

Calling p(x) evaluates the polynomial at the point x = (x[1], …, x[N]) by the formula

\[p(x_1,\dots,x_D)=\sum_{i_j=0}^{d_j}c_{i_1\dots i_D}\prod_{j=1}^D\binom{d_j}{i_j}(x_j-l_j)^{i_j}(r_j-x_j)^{d_j-i_j}\]

where $l_j = lc[j]$ and $r_j = hc[j]$ are the lower and upper bounds of the hyperrectangle, respectively, and $d_j = size(c)[j] - 1$ is the degree of the polynomial in dimension j.

See also berninterp.

BernsteinPolynomial(p::Polynomial, lb, ub)

Convert a Polynomial p with coefficients in the monomial basis to a Bernstein polynomial on the domain defined by the lower corner lb and upper corner ub.

@enum CellType

Enumeration for different types of cells. Options are full_cell, empty_cell, and partial_cell.

struct Config

The Config struct represents the configuration for implicit integration, passed to the integrate function to customize the behavior of the algorithm. It contains the following fields:

• find_zero a function with signature (f, a, b, tol) --> x such that f(x) ≈ 0, a ≤ x ≤ b. The tolerance tol is used to specify the absolute tolerance of the zero approximation (e.g. xatol in Roots).
• quad: a function with signature quad(f, a, b, tol) --> (I,E) such that I approximates the integral of f over [a,b] and E is the estimated error. a and b Tuple(s)/SVector(s) specifying the lower and upper bounds of the integration domain, and tol is the desired absolute tolerance.
• min_vol: a function with signature (tol) --> Float64 used to specify the volume of boxes below which the spatial subdivision stops. low-order method is used to approximate the integral.
• min_qual: a number between 0 and 1 used to specify the minimum quality factor for a height direction to be considered valid for recursion. The quality factor for a direction k is given by |∂ₖϕ| / |∇ϕ|.

GaussLegendre(; order, T = Float64)

Construct a Gauss-Legendre quadrature rule of given order with nodes and weights of type T, callable as quad1d(f, a, b) where f is the function to integrate and a and b are the bounds of the integration interval.

Examples

quad1d = ImplicitIntegration.GaussLegendre(; order = 20)
quad1d(cos, 0, 1) ≈ sin(1)

struct HyperRectangle{N,T<:AbstractFloat}

A struct representing a hyperrectangle in N-dimensional space.

Fields

• lc::SVector{N,T}: The lower corner of the hyperrectangle.
• hc::SVector{N,T}: The upper corner of the hyperrectangle.

struct LogInfo

A structure to store logging information for integration processes.

Fields

• subdivisions::Vector{Int}: A vector containing the subdivisions per dimension used during the integration process.
• loworder::Int: The number of times the low-order method was used.

LogInfo(U::HyperRectangle)

Initialize a LogInfo object for integrating over U.

struct Polynomial{D,T}

D-dimensional polynomial with coefficients of type T. The coefficients are stored as a dense array, and are implicitly associated with a monomial basis; i.e. the c[I] coefficient multiplies the monomial term prod(x.^(I .- 1)), where I is a D-dimensional multi-index.

Passing a Dict{NTuple{N,Int},T} to the constructor will create a polynomial with coefficients given by c[k] for the monomial x₁^k[1] * x₂^k[2] * ... * x_N^k[N], where k is a multi-index of length N corresponding to the key, and c::T is the value associated with that key.

struct Quadrature{N,T}

A collection of coords and weights for integration in N dimensions.

Quadratures are the result of quadgen and can be used to integrate functions through integrate(f,::Quadrature).

const Segment{T} = HyperRectangle{1,T}

A one-dimensional hyperrectangle.

TensorQuadrature(quad1d)

Given a 1D quadrature rule quad1d with signature quad1d(f, a::Number, b::Number), return a tensor product quadrature rule callable through quadnd(f,a::SVector{N}, b::SVector{N}).

mutable struct TreeNode

_find_zeros!(roots, f, U::Segment, config, tol)

Return all zeros of the function f in the Segment U. f should be callable as f(x::SVector{1}).

_integrand_eval(f, phi_vec, s_vec, U, k, config, ::Type{RET_TYPE}, tol, logger)

Return a function f̃ : ℝᴺ⁻¹ -> ℝ that approximates the one-dimensional integral over I(x̃) ⊂ ℝ of the function t -> f(insert(x̃, k, t)), where the integration domain I is defined as I(x̃) = {t ∈ [a,b] : sᵢ*ϕᵢ(insert(̃x,k,t) ≥ 0 ∀ (ϕᵢ,sᵢ) ∈ V)}.

_monomial_to_bernstein(c::Array, lb, ub)

Convert a polynomial with coefficients c in the monomial basis to a Bernstein polynomial on the hyperrectangle defined by lb and ub.

berninterp([T=Float64,] vals::Array, lb, ub)

Construct a Bernstein polynomial of that interpolates the values vals at the points given by uniform_points(size(vals), lb, ub), where lb and ub are the lower and upper corners of the hyperrectangle on which the polynomial is defined.

The optional type parameter T specifies the precision used for computing the polynomial coefficients.

Examples

using StaticArrays
f = (x) -> (1 - x[1])^2 + x[1]^4 + x[2]^5 * x[1]^3
lb = SVector(0.1, -0.3)
ub = SVector(1.2, 1.7)
pts = ImplicitIntegration.uniform_points((5, 6), lb, ub)
vals = f.(pts)
p = ImplicitIntegration.berninterp(vals, lb, ub)
x = lb .+ (ub - lb) .* rand(SVector{2})
f(x) ≈ p(x)

# output

true

berninterp([T=Float64,] f, n, lb, ub)

Construct a Bernstein polynomial of degree n that interpolates the function f at the points given by uniform_points(n, lb, ub), where lb and ub are the lower and upper corners of the hyperrectangle on which the polynomial is defined. The optional type parameter T specifies the precision used for computing the polynomial coefficients.

bernstein_interp(vals, pts, lb, ub)

Construct a Bernstein polynomial on (lb[1], ub[1]) × … × (lb[N], ub[N]) that interpolates the values vals at the points pts. Note that vals and pts...

bound(f, lc, hc) --> (lb, ub)

Return a lower and upper bound for the function f : U → ℝ valid for all x ∈ U.

bound(∇f, lc, hc) --> bnds

Compute a lower and upper bound for the gradient of a function f : U → ℝ valid for all x ∈ U in the sense that bnds[i][1] ≤ ∂f/∂xᵢ(x) ≤ bnds[i][2].

bound(p::BernsteinPolynomial)

Return the bounds of the Bernstein polynomial p as a tuple (m, M), where m is a lower bound and M is an upper bound on the polynomial's values over the hyperrectangle with corners low_corner(p) and high_corner(p).

bounds(rect::HyperRectangle)

Get the lower and upper bounds of a HyperRectangle.

Arguments

• rect::HyperRectangle: The HyperRectangle object.

Returns

• A tuple (lc, hc) representing the lower and upper bounds of the HyperRectangle.

cell_type(ϕ, s, U, surface)

Compute the CellType of a cell defined by the level-set function ϕ, a sign s, and the box U. If surface is true, then the cell is classified as per a surface integral.

derivative(p::BernsteinPolynomial, d::Int)

Compute the derivative along dimension d of the Bernstein polynomial p, returning a new BernsteinPolynomial of the same dimension N.

disable_default_interface()

Throw an error if the default interface is used, to ensure that all methods are implemented for a given function type. This is useful for testing purposes.

enable_default_interface()

Re-enable the default interface, so that the generic bound, gradient, project, and split methods can be used with functions that do not implement these methods.

See also disable_default_interface.

gradient(f)

Compute the gradient function f : ℝᵈ → ℝ. The returned function takes a vector x ∈ ℝᵈ and returns the gradient ∇f(x) ∈ ℝᵈ.

integrate(f, ϕ, lc, hc; tol=1e-8, surface=false, config = Config(), loginfo = false) -->
(val, logger)

Integrate the function f over an implicit domain defined by:

• Ω = {lc ≤ 𝐱 ≤ hc: ϕ(𝐱) < 0} if surface = false
• Γ = {lc ≤ 𝐱 ≤ hc: ϕ(𝐱) = 0} if surface = true

where lc and hc denote the lower and upper corners of the bounding box.

tol specifies the desired (absolute) tolerance for the approximation.

The function returns a named tuple (val, logger) where val is the approximated value, and logger is a LogInfo object containing information about the integration process if loginfo = true; otherwise, logger is nothing.

For a finer control over the integration process, pass a config object (see Config).

Note that both f and ϕ must be callable with a single argument 𝐱 of type SVector. Furthemore, ϕ is expected to return a real value.

See also quadgen if you want to generate a quadrature instead of directly computing the value of the integral.

By default, ImplicitIntegration uses ForwardDiff to compute gradients and IntervalArithmetic to compute bounds, both of which are needed for the algorithm to work. While these work reasonably well in most cases, you may want to overload the following methods for the type of your input function ϕ:

• ϕ(x::SVector{N,<:Real}) -> Real to evaluate the level-set function at x.
• ϕ(xI::SVector{N,<:Interval{<:Real}}) -> Interval{<:Real} to evaluate a bound on ϕ on the interval xI.
• ϕ(xD::SVector{N,Dual{N,<:Real}}) -> Dual{N,<:Real} to evaluate the level-set function and its gradient at x.
• ϕ(xDI::SVector{N,Dual{N,<:Interval{<:Real}}}) -> Dual{N,<:Interval{<:Real}} to evaluate a bound on ϕ and its gradient on the interval xDI.

You may need to overload the methods above if typeof(ϕ) is not supported by ForwardDiff and/or IntervalArithmetic, or if you have a better/faster implementation; see the main docs.

Examples

To compute the area of a quarter of a disk of radius 1.0:

a, b = (0.0, 0.0), (1.5, 1.5)
ϕ = (x) -> x[1]^2 + x[2]^2 - 1
f = (x) -> 1.0
res = integrate(f, ϕ, a, b) # area of quarter of a disk
res.val ≈ π / 4

To compute the perimeter of a quarter of a circle of radius 1.0:

a, b = (0.0, 0.0), (1.5, 1.5)
ϕ = (x) -> x[1]^2 + x[2]^2 - 1
f = (x) -> 1.0
res = integrate(x -> 1.0, ϕ, a, b; surface = true) # perimeter of quarter of a circle
res.val ≈ 2π / 4

integrate(f, Q::Quadrature)

Shorthand for ∑ᵢf(xᵢ)wᵢ, where xᵢ and wᵢ are the nodes and weights of the Quadrature.

lift_all_dims(root::TreeNode)

Lift all nodes in root to the same dimension as root.

lift_dim(node, dim)

Lift all boxes of dimension dim in the tree rooted at node to the same dimension as their parent.

lift_node_and_children(node::TreeNode, dim, lb, ub)

Lift each box in the tree rooted at node by appending lb and ub to the dim-th dimension of the box.

lower_restrict(p::BernsteinPolynomial{D}, d::Integer) where {D}

Restricts the given BernsteinPolynomial p to the lower face (i.e., where the d-th coordinate is at its lower bound) along the specified dimension d. Returns a new BernsteinPolynomial of dimension D-1 representing the restriction.

 project(f, k, v)

Given a function f : ℝᵈ → ℝ, a value v ∈ ℝ and an integer 1 ≤ k ≤ d, return the function f̃ : ℝᵈ⁻¹ → ℝ defined by projecting f onto the hyperplane xₖ = v; i.e. f̃(x) = f(x₁, ..., x_{k-1}, v, x_{k}, ..., x_d).

quadgen(ϕ, lc, hc; order, surface, kwargs...)

Similar to integrate, but generate a Quadrature over an implict domain defined by:

• Ω = {lc ≤ 𝐱 ≤ hc: ϕ(𝐱) < 0} if surface = false
• Γ = {lc ≤ 𝐱 ≤ hc: ϕ(𝐱) = 0} if surface = true

order specifies the degree of exactness of the quadrature rule; that is, the quadrature rule will integrate exactly polynomials of degree up to order, but not order+1. A GaussLegendre quadrature rule is used.

The function returns a named tuple (quad, logger) where quad contains the generated Quadrature, and logger is a LogInfo object containing information about the integration process.

See integrate for more information on the available kwargs.

Examples

To compute the area of a quarter of a disk of radius 1.0:

a, b = (0.0, 0.0), (1.5, 1.5)
ϕ = (x) -> x[1]^2 + x[2]^2 - 1
f = (x) -> 1.0
out = quadgen(ϕ, a, b; order = 20)
integrate(f, out.quad) ≈ π / 4 # area of quarter of a disk

To compute the perimeter of a quarter of a circle of radius 1.0:

a, b = (0.0, 0.0), (1.5, 1.5)
ϕ = (x) -> x[1]^2 + x[2]^2 - 1
f = (x) -> 1.0
out = quadgen(ϕ, a, b; order = 20, surface = true)
Q = out.quad
integrate(f, Q) ≈ 2π / 4 # perimeter of quarter of a circle

reference_vandermonde_matrix([T=Float64,] n)

Vandermond matrix for Bernstein interpolation on uniform_points(n, lb, ub), with lb = (0, … , 0) and ub = (1, … , 1). The ordering of the basis elements is the same as in the evaluation of the Bernstein polynomial. The optional type parameter T specifies the element type of the matrix.

remove_dimension(rect::HyperRectangle, k)

Remove a dimension from a HyperRectangle by deleting the k-th element from the lower and upper corners.

Arguments

• rect::HyperRectangle: The input hyperrectangle.
• k: The index of the dimension to be removed.

Returns

A new HyperRectangle with the k-th dimension removed.

sgn(m, s, S::Bool, σ)

Helper function to compute the sign of lower and upper restrictions of a level-set function ϕᵢ in a box along a given height direction k. Here m is sign of ∂ₖϕᵢ, which is assume not to change throughout the box since k is assumed to be a height direction, s is the sign of ϕᵢ in the box, S is a flag to indicate whether we are in the special case of a surface integral.

split(f, lb, ub)

Given a function f : ℝᵈ → ℝ and lower and upper bounds lb, ub ∈ ℝᵈ, return the restriction of f to the hyperrectangle defined by lb and ub.

By default this function simply returns f, but it computing sharper bounds on the restricted function requires a more sophisticated implementation.

split(U::HyperRectangle, dir)

Split a hyperrectangle U along the specified direction dir.

Arguments

• U::HyperRectangle: The hyperrectangle to be split.
• dir: The direction along which to split the hyperrectangle.

Returns

• Uₗ: The left half of the split hyperrectangle.
• Uᵣ: The right half of the split hyperrectangle.

split(p::BernsteinPolynomial, d::Integer, α = 0.5) --> pₗ, pᵣ

Split the Bernstein polynomial p along dimension d at lc[d] + (hc[d] - lc[d]) * α, where lc and hc are the lower and upper corners of the hyperrectangle on which p is defined. Returns two new Bernstein polynomials pₗ and pᵣ representing the left and right polynomials.

uniform_points(n[, lc, hc])

Generate n[1] × … × n[D] uniformly spaced points in the hyperrectangle defined by (lc[1], hc[1]) × … × (lc[D], hc[D]), where D = length(n) = length(lc) = length(hc). By default lc = (0, … , 0) and hc = (1, … , 1), generating n[1] × … × n[D] points in the unit hypercube [0, 1]^D.

The points are returned as an array, of size n, containing SVector of dimension D.

upper_restrict(p::BernsteinPolynomial{D}, d::Integer) where {D}

Same as lower_restrict, but restricts to the upper face (i.e., where the d-th coordinate is at its upper bound)

use_heuristic_bounds(F::Type, n = 10)

Overload interface methods for type F to use a sample-based heuristic when bounding functions of type F and its gradient.

The bounds are obtained by sampling the function (or its gradient) on an n × ... × n grid and taking the maximum/minimum of the attained values. This is obviously a heuristic, and may fail in practice.