Advanced usage

We now cover the options available for the integrate function.

Buffering

When calling integrate(f, domain), the package allocates memory for storing the various subregions that are generated during the adaptive integration process. Here is what it looks like in practice:

using HAdaptiveIntegration
using BenchmarkTools

t = Triangle((0, 0), (1, 0), (0, 1))
f = x -> 1 / (x[1]^2 + x[2]^2 + 1e-2)
@benchmark integrate($f, $t)

BenchmarkTools.Trial: 10000 samples with 1 evaluation per sample.
 Range (min … max):  27.661 μs …  3.344 ms  ┊ GC (min … max): 0.00% … 97.60%
 Time  (median):     29.014 μs              ┊ GC (median):    0.00%
 Time  (mean ± σ):   33.196 μs ± 50.972 μs  ┊ GC (mean ± σ):  3.13% ±  2.31%

  ▅██▅▂▄▅▂▂▅▄▃▁▂▃▂           ▁▂▁▂▃▄▃▃▂▂▂▁▂▂▁▁                 ▂
  ███████████████████▇▆▆▇▇▇█████████████████████▆▆▆▅▄▄▂▅▅▄▄▅▆ █
  27.7 μs      Histogram: log(frequency) by time      49.8 μs <

 Memory estimate: 44.85 KiB, allocs estimate: 9.

While the overhead associated to these (small) allocations is usually negligible, there are circumstances where one may want to avoid allocations altogether. This can be achieved by passing a buffer to the integrate using allocate_buffer:

using HAdaptiveIntegration: allocate_buffer

buffer = allocate_buffer(f, t)
integrate(f,t; buffer)
b = @benchmark integrate($f, $t; buffer = $buffer)

BenchmarkTools.Trial: 10000 samples with 1 evaluation per sample.
 Range (min … max):  26.369 μs … 69.078 μs  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     26.590 μs              ┊ GC (median):    0.00%
 Time  (mean ± σ):   26.912 μs ±  1.769 μs  ┊ GC (mean ± σ):  0.00% ± 0.00%

  ▆█▆▁  ▁                                          ▁          ▁
  ████▄▅██▅▃▄▃▁▄▄▆▆▅▃▄▁▃▄▁▄▄▄▃▁▁▃▁▁▃▁▄▁▄▃▁▃▃▁▁▁▃▁▅▇██▇▆▅▅▅▅▅▆ █
  26.4 μs      Histogram: log(frequency) by time      34.9 μs <

 Memory estimate: 0 bytes, allocs estimate: 0.

Provided evaluating f does not allocate, and the buffer has a sufficiently large capacity, integrate will not allocate memory during the integration process, as shown in the benchmark above.

Callback

The callback keyword argument allows you to monitor the progress of the adaptive integration. The callback function is called for each estimated value of I and E, including the initial estimate (nb_subdiv=0) and after each subdivision. It receives the current integral I, error estimate E, number of subdivisions nb_subdiv, and the internal buffer.

Here is a simple example that collects the convergence history:

using HAdaptiveIntegration

t = Triangle((0, 0), (1, 0), (0, 1))
f = x -> 1 / (x[1]^2 + x[2]^2 + 0.5)

history = @NamedTuple{I::Float64, E::Float64, nb_subdiv::Int}[]

I, E = integrate(f, t; callback = (I, E, nb_subdiv, _) -> push!(history, (; I, E, nb_subdiv)))

using Printf
@printf("  %4s | %14s | %12s\n", "step", "I", "E")
@printf("  %4s-+-%14s-+-%12s\n", "----", "--------------", "------------")
foreach(history) do h
    @printf("  %4d | %14.10f | %12.4e\n", h.nb_subdiv, h.I, h.E)
end

  step |              I |            E
  -----+----------------+-------------
     0 |   0.6395139100 |   9.1952e-05
     1 |   0.6395103452 |   1.3732e-05
     2 |   0.6395103678 |   6.8358e-06
     3 |   0.6395103519 |   1.5737e-06
     4 |   0.6395103519 |   8.8150e-07
     5 |   0.6395103519 |   1.8933e-07
     6 |   0.6395103519 |   1.1727e-07
     7 |   0.6395103519 |   5.6892e-08
     8 |   0.6395103519 |   4.6596e-08
     9 |   0.6395103519 |   3.6418e-08
    10 |   0.6395103519 |   2.9304e-08
    11 |   0.6395103519 |   2.2191e-08
    12 |   0.6395103519 |   1.6106e-08
    13 |   0.6395103519 |   1.0022e-08
    14 |   0.6395103519 |   8.2339e-09

Embedded cubature formulas

By default, when calling integrate(f, domain), the package uses a default embedded cubature formula for the given domain by calling default_embedded_cubature. Although these are generally good choices, you can also specify a custom embedded cubature formula by passing it as a keyword argument to integrate. For example, in the case of a triangle, the package defaults to a Radon-Laurie embedded cubature formula of orders 5 and 7 (see the rules_simplex.jl file for more details). If you want e.g. to use an embedded cubature based on the GrundmannMoeller rule of order 13, you can do

using HAdaptiveIntegration
using HAdaptiveIntegration.Rule: GrundmannMoeller, embedded_cubature

t = Triangle((0, 0), (1, 0), (0, 1))
f = x -> 1 / (x[1]^2 + x[2]^2 + 1e-2)
ec = embedded_cubature(GrundmannMoeller{2}(13, 11))
I, E = integrate(f, t; embedded_cubature = ec)

(3.2580694607121954, 4.863250004094599e-9)

Which cubature rule is best depends on the function being integrated, as well as on the desired accuracy; as a rule of thumb, higher-order cubature rules will perform better for globally smooth functions f or higher accuracy requirements. Here is a short study on the number of function evaluations required to achieve a given accuracy for the default Radon-Laurie cubature and the GrundmannMoeller cubature rule above:

const cc = Ref(0) # a counter for the number of function evaluations
g = x -> (cc[] += 1; f(x))
rtol = 1e-2
cc[] = 0; integrate(g, t; rtol); counter_default = cc[]
cc[] = 0; integrate(g, t; embedded_cubature = ec, rtol); counter_custom = cc[]
counter_default, counter_custom

(171, 395)

For rtol = 1e-2, we see that the default cubature rule requires fewer function evaluations. However, decreasing rtol changes the balance:

rtol = 1e-8
cc[] = 0; integrate(g, t; rtol); counter_default = cc[]
cc[] = 0; integrate(g, t; embedded_cubature = ec, rtol); counter_custom = cc[]
counter_default, counter_custom

(8455, 2291)

This example illustrates that testing is necessary to determine which cubature rule is best for your specific application!

Available embedded cubature formulas

The list of available embedded cubature formulas is:

CUBE_BE115
CUBE_BE65
GenzMalik
GrundmannMoeller
RadonLaurie
SEGMENT_GK15
SEGMENT_GK31
SEGMENT_GK7
SQUARE_CH21
SQUARE_CH25

To add a custom embedded cubature for a given domain, you must write a constructor e.g. my_custom_cubature(args...) that returns a valid EmbeddedCubature object (see the function embedded_cubature in the file Rule/triangle.jl for some examples on how this is done). PRs with new schemes are more than welcome!

Subdivision strategies

The package uses a default subdivision strategy for the given domain by calling default_subdivision. For example, by default triangles are subdivided into 4 smaller triangles by connecting the midpoints of the edges:

using HAdaptiveIntegration

t = Triangle((0, 0), (1, 0), (0, 1))
subdiv_algo = HAdaptiveIntegration.default_subdivision(t)

subdivide_triangle (generic function with 1 method)

Here are the subdivided triangles:

subdiv_algo(t)

(Triangle{Float64}(StaticArraysCore.SVector{2, Float64}[[0.0, 0.0], [0.5, 0.0], [0.0, 0.5]]), Triangle{Float64}(StaticArraysCore.SVector{2, Float64}[[1.0, 0.0], [0.5, 0.5], [0.5, 0.0]]), Triangle{Float64}(StaticArraysCore.SVector{2, Float64}[[0.0, 1.0], [0.0, 0.5], [0.5, 0.5]]), Triangle{Float64}(StaticArraysCore.SVector{2, Float64}[[0.5, 0.0], [0.5, 0.5], [0.0, 0.5]]))

But is also possible (and maybe desirable) to split the triangle into 2 smaller triangles instead. The following function accomplishes this:

using StaticArrays

function subdivide_triangle2(t::Triangle{T}) where {T}
        a, b, c = t.vertices
        bc = (b + c) / 2
        return (Triangle{T}(SVector(bc, a, b)), Triangle{T}(SVector(bc, c, a)))
    end

subdivide_triangle2 (generic function with 1 method)

Passing subdivide_triangle2 as the subdiv_algo to integrate will use this instead of the default:

f = x -> 1 / (x[1]^2 + x[2]^2 + 1e-2)
I, E = integrate(f, t; subdiv_algo = subdivide_triangle2)

(3.258069460313913, 4.7877675435348266e-8)

Which subdivision strategy is best depends on the function being integrated; for the example presented above, it turns out the default strategy is more efficient!