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):  37.891 μs …  2.913 ms  ┊ GC (min … max): 0.00% … 97.16%
 Time  (median):     39.073 μs              ┊ GC (median):    0.00%
 Time  (mean ± σ):   42.738 μs ± 48.309 μs  ┊ GC (mean ± σ):  2.96% ±  4.61%

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

 Memory estimate: 58.30 KiB, allocs estimate: 8.

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):  37.530 μs … 64.731 μs  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     38.141 μs              ┊ GC (median):    0.00%
 Time  (mean ± σ):   38.517 μs ±  1.868 μs  ┊ GC (mean ± σ):  0.00% ± 0.00%

   ▂▇██▅▁                                           ▁▁▁       ▂
  ▄█████████▄▃▃▁▄▄▇▇▅▅▄▁▁▁▁▁▃▄▁▄▁▃▃▃▃▁▄▃▃▃▄▃▃▃▃▁▁▄▇█████▇▆▅▆▅ █
  37.5 μs      Histogram: log(frequency) by time      46.6 μ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.

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.258069460712197, 4.863253362519249e-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!

CUBE_BE115
CUBE_BE65
CUBE_GM33
GenzMalik
GrundmannMoeller
RadonLaurie
SQUARE_CH21
SQUARE_CH25
SQUARE_GM17
TETRAHEDRON_GM35
TRIANGLE_GM19
TRIANGLE_RL19

To add a custom embedded quadrature 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.2580694603139104, 4.787767539718435e-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!