Integrating a vector valued out-of-place integrand returns a scalar

[Sleort]

Describe the bug 🐞

A 1D integral of a vector valued out-of-place function returns a scalar, even when a vector valued integrand prototype is provided. The equivalent in-place problem works as expected.

Expected behavior

When a vector valued integrand prototype is provided, the solution should be the same for out-of-place and in-place integrands.

Minimal Reproducible Example 👇

This works as expected:

using Integrals

integrand_prototype = zeros(2)

#In-place integrand
f!(y,u,p) = @. y = exp(-u^2)
test! = IntegralFunction(f!, integrand_prototype)
solve(IntegralProblem(test!, -Inf, Inf), QuadGKJL()).u #Vector return value, as expected

while this doesn't:

#Out-of-place integrand
f(u,p) = @. exp(-u^2)
test = IntegralFunction(f, integrand_prototype)
solve(IntegralProblem(test, -Inf, Inf), QuadGKJL()).u #Scalar return value!?

Environment

• Output of using Pkg; Pkg.status()

Status `/tmp/jl_pPHahC/Project.toml`
  [de52edbc] Integrals v4.4.1

• Output of using Pkg; Pkg.status(; mode = PKGMODE_MANIFEST)

Status `/tmp/jl_pPHahC/Manifest.toml`
  [47edcb42] ADTypes v1.2.1
  [7d9f7c33] Accessors v0.1.36
  [79e6a3ab] Adapt v4.0.4
  [66dad0bd] AliasTables v1.1.3
  [4fba245c] ArrayInterface v7.10.0
  [49dc2e85] Calculus v0.5.1
  [861a8166] Combinatorics v1.0.2
  [38540f10] CommonSolve v0.2.4
  [34da2185] Compat v4.15.0
  [a33af91c] CompositionsBase v0.1.2
  [2569d6c7] ConcreteStructs v0.2.3
  [187b0558] ConstructionBase v1.5.5
  [9a962f9c] DataAPI v1.16.0
  [864edb3b] DataStructures v0.18.20
  [e2d170a0] DataValueInterfaces v1.0.0
  [8bb1440f] DelimitedFiles v1.9.1
  [31c24e10] Distributions v0.25.108
  [ffbed154] DocStringExtensions v0.9.3
  [fa6b7ba4] DualNumbers v0.6.8
  [4e289a0a] EnumX v1.0.4
  [e2ba6199] ExprTools v0.1.10
  [1a297f60] FillArrays v1.11.0
  [069b7b12] FunctionWrappers v1.1.3
  [77dc65aa] FunctionWrappersWrappers v0.1.3
  [46192b85] GPUArraysCore v0.1.6
  [19dc6840] HCubature v1.6.0
  [34004b35] HypergeometricFunctions v0.3.23
  [18e54dd8] IntegerMathUtils v0.1.2
  [de52edbc] Integrals v4.4.1
  [3587e190] InverseFunctions v0.1.14
  [92d709cd] IrrationalConstants v0.2.2
  [82899510] IteratorInterfaceExtensions v1.0.0
  [692b3bcd] JLLWrappers v1.5.0
  [73f95e8e] LatticeRules v0.0.1
  [2ab3a3ac] LogExpFunctions v0.3.27
  [1914dd2f] MacroTools v0.5.13
  [e1d29d7a] Missings v1.2.0
  [4886b29c] MonteCarloIntegration v0.2.0
  [77ba4419] NaNMath v1.0.2
  [bac558e1] OrderedCollections v1.6.3
  [90014a1f] PDMats v0.11.31
  [aea7be01] PrecompileTools v1.2.1
  [21216c6a] Preferences v1.4.3
  [27ebfcd6] Primes v0.5.6
  [43287f4e] PtrArrays v1.2.0
  [1fd47b50] QuadGK v2.9.4
  [8a4e6c94] QuasiMonteCarlo v0.3.3
  [3cdcf5f2] RecipesBase v1.3.4
  [731186ca] RecursiveArrayTools v3.19.0
  [189a3867] Reexport v1.2.2
  [ae029012] Requires v1.3.0
  [79098fc4] Rmath v0.7.1
  [7e49a35a] RuntimeGeneratedFunctions v0.5.13
  [0bca4576] SciMLBase v2.39.0
  [c0aeaf25] SciMLOperators v0.3.8
  [53ae85a6] SciMLStructures v1.2.0
  [efcf1570] Setfield v1.1.1
  [ed01d8cd] Sobol v1.5.0
  [a2af1166] SortingAlgorithms v1.2.1
  [276daf66] SpecialFunctions v2.4.0
  [90137ffa] StaticArrays v1.9.4
  [1e83bf80] StaticArraysCore v1.4.2
  [82ae8749] StatsAPI v1.7.0
  [2913bbd2] StatsBase v0.34.3
  [4c63d2b9] StatsFuns v1.3.1
  [2efcf032] SymbolicIndexingInterface v0.3.22
  [3783bdb8] TableTraits v1.0.1
  [bd369af6] Tables v1.11.1
  [efe28fd5] OpenSpecFun_jll v0.5.5+0
  [f50d1b31] Rmath_jll v0.4.2+0
  [0dad84c5] ArgTools v1.1.1
  [56f22d72] Artifacts
  [2a0f44e3] Base64
  [ade2ca70] Dates
  [8ba89e20] Distributed
  [f43a241f] Downloads v1.6.0
  [7b1f6079] FileWatching
  [9fa8497b] Future
  [b77e0a4c] InteractiveUtils
  [b27032c2] LibCURL v0.6.4
  [76f85450] LibGit2
  [8f399da3] Libdl
  [37e2e46d] LinearAlgebra
  [56ddb016] Logging
  [d6f4376e] Markdown
  [a63ad114] Mmap
  [ca575930] NetworkOptions v1.2.0
  [44cfe95a] Pkg v1.10.0
  [de0858da] Printf
  [3fa0cd96] REPL
  [9a3f8284] Random
  [ea8e919c] SHA v0.7.0
  [9e88b42a] Serialization
  [6462fe0b] Sockets
  [2f01184e] SparseArrays v1.10.0
  [10745b16] Statistics v1.10.0
  [4607b0f0] SuiteSparse
  [fa267f1f] TOML v1.0.3
  [a4e569a6] Tar v1.10.0
  [8dfed614] Test
  [cf7118a7] UUIDs
  [4ec0a83e] Unicode
  [e66e0078] CompilerSupportLibraries_jll v1.1.1+0
  [deac9b47] LibCURL_jll v8.4.0+0
  [e37daf67] LibGit2_jll v1.6.4+0
  [29816b5a] LibSSH2_jll v1.11.0+1
  [c8ffd9c3] MbedTLS_jll v2.28.2+1
  [14a3606d] MozillaCACerts_jll v2023.1.10
  [4536629a] OpenBLAS_jll v0.3.23+4
  [05823500] OpenLibm_jll v0.8.1+2
  [bea87d4a] SuiteSparse_jll v7.2.1+1
  [83775a58] Zlib_jll v1.2.13+1
  [8e850b90] libblastrampoline_jll v5.8.0+1
  [8e850ede] nghttp2_jll v1.52.0+1
  [3f19e933] p7zip_jll v17.4.0+2

• Output of versioninfo()

Julia Version 1.10.3
Commit 0b4590a5507 (2024-04-30 10:59 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 16 × Intel(R) Core(TM) i7-10875H CPU @ 2.30GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-15.0.7 (ORCJIT, skylake)
Threads: 1 default, 0 interactive, 1 GC (on 16 virtual cores)

[lxvm]

Hi,

To me it looks like your integrand is behaving as expected because QuadGK.jl gives your integrand a scalar, u, and so your broadcasting function f will also return a scalar because there is nothing to give it a shape, unlike f! which gets its shape from y. Does this help?

[Sleort]

Hmm...Yeah... I guess I just misunderstood the documentation saying that

If integrand_prototype is present for either in-place or out-of-place integrands it is used to infer the return type of the integrand.

In other words, I imagined integrand_prototype would provide the shape...

Maybe cases like this one, i.e. when typeof(solution(...).u) != typeof(integrand_prototype), should error instead?

[lxvm]

Yes, we use integrand_prototype to determine the type and shape of the integrand in some cases, in particular for C libraries or to pre-allocate workspaces with correct types for algorithms. Also, we currently assume the output of the user's function is consistent with the prototype, but we don't assert it because the type of the solution may depend on the library used. Any pr to check for correctness would be welcome as long as the type assertions don't incur runtime overhead for type-stable use-cases.