Advanced usage

We now cover the options available for the integrate function.

Reduce memory allocations

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):  20.090 μs …  3.253 ms  ┊ GC (min … max): 0.00% … 98.07%
 Time  (median):     21.362 μs              ┊ GC (median):    0.00%
 Time  (mean ± σ):   23.793 μs ± 50.392 μs  ┊ GC (mean ± σ):  5.74% ±  4.34%

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

 Memory estimate: 44.85 KiB, allocs estimate: 9.

While the overhead associated with 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):  20.951 μs … 103.706 μs  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     21.031 μs               ┊ GC (median):    0.00%
 Time  (mean ± σ):   21.629 μs ±   2.356 μs  ┊ GC (mean ± σ):  0.00% ± 0.00%

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

Track convergence progress

The callback keyword argument allows you to monitor the progress of the adaptive integration. The callback function is called after each subdivision and receives:

Here is a practical example to print the history convergence:

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 | %16s | %12s\n", "step", "I", "E")
@printf("  %4s-+-%16s-+-%12s\n", "----", "----------------", "------------")
foreach(history) do h
    @printf("  %4d | %16.12f | %12.4e\n", h.nb_subdiv, h.I, h.E)
end

  step |                I |            E
  -----+------------------+-------------
     0 |   0.639513909998 |   9.1952e-05
     1 |   0.639510345210 |   1.3732e-05
     2 |   0.639510367798 |   6.8358e-06
     3 |   0.639510351940 |   1.5737e-06
     4 |   0.639510351905 |   8.8150e-07
     5 |   0.639510351869 |   1.8933e-07
     6 |   0.639510351891 |   1.1727e-07
     7 |   0.639510351872 |   5.6892e-08
     8 |   0.639510351875 |   4.6596e-08
     9 |   0.639510351871 |   3.6418e-08
    10 |   0.639510351870 |   2.9304e-08
    11 |   0.639510351868 |   2.2191e-08
    12 |   0.639510351870 |   1.6106e-08
    13 |   0.639510351871 |   1.0022e-08
    14 |   0.639510351871 |   8.2339e-09

Choose custom cubature rules

By default, when calling integrate(f, domain), the package uses a default embedded cubature formula for the given domain by calling default_rule. 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 high order 8 and low order 5 (see RadonLaurie). If you want e.g. to use an embedded cubature based on the GrundmannMoeller rule of high order 13 and low order 11, 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; rule = ec)

(3.2580694607121967, 4.863251044928685e-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; rule = 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; rule = 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
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. See the file at Rule/triangle.jl for examples. PRs with new schemes are more than welcome!

Define custom 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}[[0.5, 0.0], [1.0, 0.0], [0.5, 0.5]]), Triangle{Float64}(StaticArraysCore.SVector{2, Float64}[[0.0, 0.5], [0.5, 0.5], [0.0, 1.0]]), Triangle{Float64}(StaticArraysCore.SVector{2, Float64}[[0.5, 0.5], [0.0, 0.5], [0.5, 0.0]]))

But it 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

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

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!