diff --git a/README.md b/README.md index bf8d4c73..e8155a7a 100644 --- a/README.md +++ b/README.md @@ -57,10 +57,10 @@ multithreading with Cubature.jl: using Integrals, Cubature, Base.Threads function f(dx, x, p) Threads.@threads for i in 1:size(x, 2) - dx[i] = sum(sin.(@view(x[:, i]))) + dx[i] = sum(sin, @view(x[:, i])) end end -prob = IntegralProblem(f, ones(2), 3ones(2), batch = 2) +prob = IntegralProblem(BatchIntegralFunction(f, zeros(0)), ones(2), 3ones(2)) sol = solve(prob, CubatureJLh(), reltol = 1e-3, abstol = 1e-3) ``` diff --git a/docs/pages.jl b/docs/pages.jl index 99ef1719..d8c8a7d6 100644 --- a/docs/pages.jl +++ b/docs/pages.jl @@ -2,6 +2,7 @@ pages = ["index.md", "Tutorials" => Any["tutorials/numerical_integrals.md", "tutorials/differentiating_integrals.md"], "Basics" => Any["basics/IntegralProblem.md", + "basics/IntegralFunction.md", "basics/SampledIntegralProblem.md", "basics/solve.md", "basics/FAQ.md"], diff --git a/docs/src/basics/FAQ.md b/docs/src/basics/FAQ.md index c996e4b9..baaff810 100644 --- a/docs/src/basics/FAQ.md +++ b/docs/src/basics/FAQ.md @@ -1,3 +1,63 @@ # Frequently Asked Questions -Ask more questions. +## How should I use the in-place interface? + +The in-place interface allows evaluating vector-valued integrands without +allocating an output array. This can be beneficial for reducing allocations when +integrating many functions simultaneously or to make use of existing in-place +code. However, note that not all algorithms use in-place operations under the +hood, i.e. `HCubatureJL()`, and may still allocate. + +You can construct an `IntegralFunction(f, prototype)`, where `f` is of the form +`f(y, u, p)` where `prototype` is of the desired type and shape of `y`. + +For small array outputs of a known size, consider using StaticArrays.jl for the +return value of your integrand. + +## How should I use the batch interface? + +The batch interface allows evaluating one (or more) integrals simultaneously at +different points, which maximizes the parallelism for a given algorithm. + +You can construct an out-of-place `BatchIntegralFunction(bf)` where `bf` is of +the form `bf(u, p) = stack(x -> f(x, p), eachslice(u; dims=ndims(u)))`, where +`f` is the (unbatched) integrand. + +You can construct an in-place `BatchIntegralFunction(bf, prototype)`, where `bf` +is of the form `bf(y, u, p) = foreach((y,x) -> f(y,x,p), eachslice(y, dims=ndims(y)), eachslice(x, dims=ndims(x)))`. + +Note that not all algorithms use in-place batched operations under the hood, +i.e. `QuadGKJL()`. + +## What should I do if my solution is not converged? + +Certain algorithms, such as `QuadratureRule` used a fixed number of points to +calculate an integral and cannot provide an error estimate. In this case, you +have to increase the number of points and check the convergence yourself, which +will depend on the accuracy of the rule you choose. + +For badly-behaved integrands, such as (nearly) singular and highly oscillatory +functions, most algorithms will fail to converge and either throw an error or +silently return the incorrect result. In some cases Integrals.jl can provide an +error code when things go wrong, but otherwise you can always check if the error +estimate for the integral is less than the requested tolerance, e.g. `sol.resid +< max(abstol, reltol*norm(sol.u))`. Sometimes using a larger tolerance or higher +precision arithmetic may help. + +## How can I integrate arbitrarily-spaced data? + +See `SampledIntegralProblem`. + +## How can I integrate on arbitrary geometries? + +You can't, since Integrals.jl currently supports integration on hypercubes +because that is what lower-level packages implement. + +## I don't see algorithm X or quadrature scheme Y ? + +Fixed quadrature rules from other packages can be used with `QuadratureRule`. +Otherwise, feel free to open an issue or pull request. + +## Can I take derivatives with respect to the limits of integration? + +Currently this is not implemented. \ No newline at end of file diff --git a/docs/src/basics/IntegralFunction.md b/docs/src/basics/IntegralFunction.md new file mode 100644 index 00000000..eee17094 --- /dev/null +++ b/docs/src/basics/IntegralFunction.md @@ -0,0 +1,6 @@ +# [Integral Functions](@id func) + +```@docs +SciMLBase.IntegralFunction +SciMLBase.BatchIntegralFunction +``` \ No newline at end of file diff --git a/docs/src/tutorials/numerical_integrals.md b/docs/src/tutorials/numerical_integrals.md index 740e79cd..fe9cf90c 100644 --- a/docs/src/tutorials/numerical_integrals.md +++ b/docs/src/tutorials/numerical_integrals.md @@ -39,15 +39,15 @@ For example, we also want to evaluate: ```@example integrate2 using Integrals f(u, p) = [sum(sin.(u)), sum(cos.(u))] -prob = IntegralProblem(f, ones(3), 3ones(3); nout = 2) +prob = IntegralProblem(f, ones(3), 3ones(3)) sol = solve(prob, HCubatureJL(); reltol = 1e-3, abstol = 1e-3) sol.u ``` -The keyword `nout` now has to be specified equal to the number of integrals we are are calculating, 2. Another way to think about this is that the integrand is now a vector valued function. -The default value for the keyword `nout` is 1, -and thus it does not need to be specified for scalar valued functions. +In general, we should be able to integrate any type that is in a vector space +and supports addition and scalar multiplication, although Integrals.jl allows +scalars and arrays. In the above example, the integrand was defined out-of-position. This means that a new output vector is created every time the function `f` is called. If we do not want these allocations, we can also define `f` in-position. @@ -58,14 +58,16 @@ function f(y, u, p) y[1] = sum(sin.(u)) y[2] = sum(cos.(u)) end -prob = IntegralProblem(f, ones(3), 3ones(3); nout = 2) +prototype = zeros(2) +prob = IntegralProblem(IntegralFunction(f, prototype), ones(3), 3ones(3)) sol = solve(prob, CubatureJLh(); reltol = 1e-3, abstol = 1e-3) sol.u ``` -where `y` is a cache to store the evaluation of the integrand. +where `y` is a cache to store the evaluation of the integrand and `prototype` is +an instance of `y` with the desired type and shape. We needed to change the algorithm to `CubatureJLh()` -because `HCubatureJL()` does not support in-position. +because `HCubatureJL()` does not support in-position under the hood. `f` evaluates the integrand at a certain point, but most adaptive quadrature algorithms need to evaluate the integrand at multiple points in each step of the algorithm. @@ -77,11 +79,12 @@ For example, here we do allocation-free multithreading with Cubature.jl: using Integrals, Cubature, Base.Threads function f(y, u, p) Threads.@threads for i in 1:size(u, 2) - y[1, i] = sum(sin.(@view(u[:, i]))) - y[2, i] = sum(cos.(@view(u[:, i]))) + y[1, i] = sum(sin, @view(u[:, i])) + y[2, i] = sum(cos, @view(u[:, i])) end end -prob = IntegralProblem(f, ones(3), 3ones(3); nout = 2, batch = 2) +prototype = zeros(2, 0) +prob = IntegralProblem(BatchIntegralFunction(f, prototype), ones(3), 3ones(3)) sol = solve(prob, CubatureJLh(); reltol = 1e-3, abstol = 1e-3) sol.u ``` @@ -89,6 +92,7 @@ sol.u Both `u` and `y` changed from vectors to matrices, where each column is respectively a point the integrand is evaluated at or the evaluation of the integrand at the corresponding point. +The `prototype` now has an extra dimension for batching that can be of size zero. Try to create yourself an out-of-position version of the above problem. For the full details of the batching interface, see the [problem page](@ref prob). @@ -109,7 +113,7 @@ The [solvers page](@ref solvers) gives an overview of which arguments each algor ## One-dimensional integrals -Integrals.jl also has specific solvers for integrals in a single dimension, such as `QuadGKLJ`. +Integrals.jl also has specific solvers for integrals in a single dimension, such as `QuadGKJL`. For example, we can create our own sine function by integrating the cosine function from 0 to x. ```@example integrate6