Merge pull request #224 from lxvm/docs

[DOCS] Update examples to use new interface

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)
 ```
 

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"],

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.

docs/src/basics/IntegralFunction.md

@@ -0,0 +1,6 @@
+# [Integral Functions](@id func)
+
+```@docs
+SciMLBase.IntegralFunction
+SciMLBase.BatchIntegralFunction
+```

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,18 +79,20 @@ 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
 ```
 
 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