Is there any documentation on "IntegralFunctions"

General Usage
, ,
1

I have been trying to solve two integrals at the same time using the Integrals.jl Package. I am following the example in the documentation of Integrals.jl in the link below.

My function that I am integrating is an in-place function that modifies a vector of size 2.

using LinearAlgebra
using Integrals, Cubature

function func!(y::Vector{Float64}, k::Float64, p::Vector{Float64})

    epsilon = k^2 - p[2];
    energy  = sqrt( epsilon^2 + p[1]^2 );
    tangent = tanh( p[3] * energy ) / energy;

    y[1]    = tangent * k^2 - 1;
    y[2]    = k^2 * ( 1 - tangent * epsilon );

end

prob    = IntegralProblem(func!, 0, 100, [1.0, 1.0, 100.0]; nout = 2);
sol     = solve(prob, HCubatureJL());

println(sol.u)

It tells me that I need to specify that the function I am integrating has two inputs and two outputs with the arument “nout = 2” but when I do so I get a warning telling me that that argument is depricated and that I should use “IntegralsFunctions” instead.

┌ Warning: `nout` and `batch` keywords are deprecated in favor of inplace `IntegralFunction`s or `BatchIntegralFunction`s. See the updated Integrals.jl documentation for details.
â”” @ SciMLBase C:\Users\laura\.julia\packages\SciMLBase\2HZ5m\src\problems\basic_problems.jl:475
[0.13095384303975866, 1.2935022147141437]

The problem is, I can not for the life of me find any documentation for how to use “IntegralFunctions”. This is my first forum post so I hope I am doing this correctly, thanks in advance.

1 Like
2

I would just use an SVector from StaticArrays.jl for this sort of thing. Then you can work out-of-place and it will be fast, and it will probably be faster with HCubatureJL (which calls HCubature.jl) too (since under the hood HCubature.jl doesn’t actually exploit in-place integrand functions).

The problem is, I can not for the life of me find any documentation for how to use “IntegralFunctions”.

Yes, a lack of documentation of IntegralFunction (which is a wrapper around your function) seems to a problem here. To learn how to use it you’d probably need to look inside the package. Or you can keep using nout as you are doing now — the deprecation warning is harmless, if annoying, and they can hardly expect you to switch to a different API that is not documented yet(?).

Maybe file an issue with Integrals.jl on the lack of documentation for IntegralFunction.

2 Likes
3

Thank you, I will probably end up working out of place for this particular problem but I may still want to figure out how to do this for some other more computationally intensive problems I will want to solve in the future.

This was one of the first real programs I tried to make in Julia so it was rather confusing that the first tutorial in the documentation uses the arguments that are depricated and won’t tell me how to fix it.

Do you know how I would go about filing such an issue with Integral.js? My coworkers recently convinced me to try Julia so I am still very new at this.

PS: I did find in the scimlfunctions.lj that I need to specify something called a “integrand_prototype” for the Integral function but I do not understand what that would be.

@doc doc"""

IntegralFunction{iip,specialize,F,T} <: AbstractIntegralFunction{iip}

A representation of an integrand `f` defined by:

math

f(u, p)


For an in-place form of `f` see the `iip` section below for details on in-place or

out-of-place handling.

julia

IntegralFunction{iip,specialize}(f, [integrand_prototype])


Note that only `f` is required, and in the case of inplace integrands a mutable container

`integrand_prototype` to store the result of the integrand. If `integrand_prototype` is

present, `f` is interpreted as in-place, and otherwise `f` is assumed to be out-of-place.

## iip: In-Place vs Out-Of-Place

Out-of-place functions must be of the form ``y = f(u, p)`` and in-place functions of the form

``f(y, u, p)``. Since `f` is allowed to return any type (e.g. real or complex numbers or

arrays), in-place functions must provide a container `integrand_prototype` that is of the

right type for the variable ``y``, and the result is written to this container in-place.

When in-place forms are used, in-place array operations, i.e. broadcasting, may be used by

algorithms to reduce allocations. If `integrand_prototype` is not provided, `f` is assumed

to be out-of-place and quadrature is performed assuming immutable return types.

## specialize

This field is currently unused

## Fields

The fields of the IntegralFunction type directly match the names of the inputs.

"""
4

I just figured out how to do it, while I was copying text to reply I noticed that the integrand prototype is just the variable that is appended by the in-place function.

So the following bit of code:

prob    = IntegralProblem(func!, 0, 100, [1.0, 1.0, 100.0]; nout = 2);

Can be replaced by:

intfunc = IntegralFunction(func!, [0.0, 0.0])
prob    = IntegralProblem(intfunc, 0, 100, [1.0, 1.0, 100.0]);

Yay :partying_face:

1 Like
5

Go to Issues · SciML/Integrals.jl · GitHub and click “New issue” (assuming you have a github.com account … if not sign up for an account and log in first), then create an issue with a title like “missing documentation for IntegralFunction”. Then, in the body of the issue, point out something like:

IntegralFunction is mentioned in the manual for Integrals.jl, and is also suggested by the nout deprecation message, but doesn’t seem to be documented in the manual. Please include it in the manual, ideally a couple of examples to illustrate how it is intended to be used.

You might also link to this discourse discussion.

An example of using IntegralFunction to define an in-place integrand can be found in one of the Integrals.jl tests, which looks like:

IntegralFunction((y, x, p) -> f_helper!(f, y, x, p), zeros(1))

Apparently, the second argument is simply a sample vector, produced by e.g. zeros, giving the shape of the desired output vector. That is, it should be zeros(nout) in your case.

Again, they should really document this — this is not your fault.

2 Likes
6

Thank you very much, I will report the issue tomorrow. :slightly_smiling_face:

Edit: Done Missing documentation for IntegralFunction and BatchIntegralFunction · Issue #218 · SciML/Integrals.jl · GitHub

7

Yes, that’s an issue on our part. When I saw your issue open I thought I’d just point to the doc page / tutorial with it, and realized it’s just completely missing. We should fix that.

8

I’m not sure whether this is worth generating a new thread for but how would you then go about batching that function? From what I understand batching allows you to evaluate the integral in chunks on different cores but the sparing information Ic an find suggests the requirement of matrices and always seems to involve solving multiple different integrals at once. Is it possible to simply parallelise a vector integral?

9

Hi, I should have seen this sooner since I am responsible for all of the updates and missing web documentation. I dedicated more time to writing good tests and neglected the documentation, but there have always been complete docstrings, e.g.

julia> using Integrals

julia>?BatchIntegralFunction

and I will put them in the web documentation.

What the BatchIntegralFunction does is evaluate the same function at multiple quadrature nodes, so if f(x::Number)::Number is your integrand, assuming it is a scalar function of one variable, a simple implementation of a batched integrand is to evaluate it “point-wise” using broadcasting: bf(bx::Vector)::Vector = f.(bx). Here I used the letter b to indicate a batched variable since when you use this interface the input bx = [x1, x2, ..., xN] is a list of evaluation points that could be of any length. In practice, this interface exists so you can parallelize the evaluation of your function, for example with threads as shown in the quadgk documentation.

The interface also depends on if the integral is multidimensional or vector-valued. For the sake of compatibility with C libraries, the rule of thumb is that the inputs and outputs of the batch integrand should be arrays whose last dimension is the batching dimension, i.e. f(bx[:, i]) == bf(bx)[:, ..., :, i], where each index on that axis corresponds to a unique evaluation point/function value. In summary, the cases are:

  • univariate, scalar function: f(x::Number)::Number becomes bf(bx::Vector)::Vector
  • multivariate, scalar function : f(x::Vector)::Number becomes bf(bx::Matrix)::Vector
  • univariate, vector function: f(x::Number)::Vector becomes bf(bx::Vector)::Matrix
  • multi-variate, vector function: f(x::Vector)::Vector becomes bf(bx::Matrix)::Matrix

In general, any array shapes should just work, i.e. f(x)::Array{Number,N} becomes bf(bx)::Array{Number,N+1}, and to use the inplace interface you have to supply a prototype of the output array with the correct size, although the length of the batching dimension can be zero. In fact for the out-of-place interface with algorithms implemented in C, your integrand will be given a length-zero batch array and has to return the intended type and shape of array correctly. I hope this helps!

1 Like