HCubature does not work with ODESolution produced by some ODE solvers #221
Comments
|
@homocomputeris was this resolved on discourse? |
Chris's solution work for me but it doesn't fix the cubature itself. The fact that the cubature need some specific array type is still present. Also
|
|
Looking into the MWE and discourse post, I don't think the issue is with Integrals.jl but with OrdinaryDiffEQ.jl. It is true that HCubature.jl tries to use StaticArrays.jl as the evaluation points, but the caller is responsible for making sure their integrand is compatible with that input. This won't be changed because it is currently the best way to get fast code for small integrands or systems of equations. As discussed on discourse, the patch to fix this is to simply make sure the ODE solver always sees a @ChrisRackauckas this seems like an OrdinaryDiffEq.jl bug, correct? Also, I saw that |
With this code module Mwe
using DifferentialEquations
using Integrals
using Distributions
using LinearAlgebra
using ForwardDiff
# just to speed things up a bit
const TOL = 1e-1
# some parameters; exact values do not matter
tspan = (0, pi)
p = [1, 1]
t = 22 / 7
# limit the domain of interest
rectangle = (-5, 5)
# ODE somewhat similar to the harmonic oscillator
function f(du, u, p, t)
du[1] = u[2]
du[2] = -u[1]
du[3] = 0.0
return du # (or return nothing)
end
# scalar function with vector argument
function g0(x, p, t)
return pdf(truncated(Normal(1, 1), rectangle...), x[1]) *
pdf(truncated(Normal(1, 1), rectangle...), x[2]) *
pdf(truncated(Normal(1, 1), rectangle...), x[3])
end
function g1(x, p, t)
# find the initial condition for a given point
function tmp(y)
prob = ODEProblem(f, y, reverse(tspan), p)
sol = solve(prob, AutoTsit5(Vern9());
abstol = TOL, reltol = TOL)
return sol(tspan[end] - t)
end
# use the found IC in g0
return g0(tmp(x), p, 0.0) * abs(det(ForwardDiff.jacobian(tmp, x)))
end
# Integrate everywhere
function x1(jpdf, x::Number, p, t)
prob = IntegralProblem((tail, p) -> jpdf([x; tail], p, t),
[-5, -Inf],
[5, Inf], p)
sol = solve(prob, HCubatureJL();
reltol = TOL, abstol = TOL)
return sol.u
end
# Truncate in a rectangle
function x1_broken(jpdf, x::Number, p, t)
prob = IntegralProblem((tail, p) -> jpdf([x; tail], p, t),
[-5, rectangle[1]],
[5, rectangle[2]], p)
sol = solve(prob, HCubatureJL();
reltol = TOL, abstol = TOL)
return sol.u
end
@show g0(rand(3), p, t)
@show g1(rand(3), p, t)
# this works
@show x1(g1, 3.0, p, t)
# this doesn't
@show x1_broken(g1, 3.0, p, t)
endg0(rand(3), p, t) = 0.04281316658626982
g1(rand(3), p, t) = 0.0015088150936285213
x1(g1, 3.0, p, t) = 7.758772211499161e-5
ERROR: Initial condition incompatible with functional form.
Detected an in-place function with an initial condition of type Number or SArray.
This is incompatible because Numbers cannot be mutated, i.e.
`x = 2.0; y = 2.0; x .= y` will error.
If using a immutable initial condition type, please use the out-of-place form.
I.e. define the function `du=f(u,p,t)` instead of attempting to "mutate" the immutable `du`.
If your differential equation function was defined with multiple dispatches and one is
in-place, then the automatic detection will choose in-place. In this case, override the
choice in the problem constructor, i.e. `ODEProblem{false}(f,u0,tspan,p,kwargs...)`.
For a longer discussion on mutability vs immutability and in-place vs out-of-place, see:
https://diffeq.sciml.ai/stable/tutorials/faster_ode_example/#Example-Accelerating-a-Non-Stiff-Equation:-The-Lorenz-Equation
Some of the types have been truncated in the stacktrace for improved reading. To emit complete information
in the stack trace, evaluate `TruncatedStacktraces.VERBOSE[] = true` and re-run the code.
Stacktrace:
[1] get_concrete_u0(prob::SciMLBase.ODEProblem{…}, isadapt::Bool, t0::Float64, kwargs::@Kwargs{…})
@ DiffEqBase ~/.julia/packages/DiffEqBase/eLhx9/src/solve.jl:1292
[2] get_concrete_problem(prob::SciMLBase.ODEProblem{…}, isadapt::Bool; kwargs::@Kwargs{…})
@ DiffEqBase ~/.julia/packages/DiffEqBase/eLhx9/src/solve.jl:1148
[3] solve_up(prob::SciMLBase.ODEProblem{…}, sensealg::Nothing, u0::StaticArraysCore.SVector{…}, p::Vector{…}, args::OrdinaryDiffEq.CompositeAlgorithm{…}; kwargs::@Kwargs{…})
@ DiffEqBase ~/.julia/packages/DiffEqBase/eLhx9/src/solve.jl:1052
[4] solve(prob::SciMLBase.ODEProblem{…}, args::OrdinaryDiffEq.CompositeAlgorithm{…}; sensealg::Nothing, u0::Nothing, p::Nothing, wrap::Val{…}, kwargs::@Kwargs{…})
@ DiffEqBase ~/.julia/packages/DiffEqBase/eLhx9/src/solve.jl:981
[5] (::Main.Mwe.var"#tmp#1"{Vector{Int64}, Float64})(y::StaticArraysCore.SVector{3, Float64})
@ Main.Mwe ./Untitled-1:38
[6] g1(x::StaticArraysCore.SVector{3, Float64}, p::Vector{Int64}, t::Float64)
@ Main.Mwe ./Untitled-1:43
[7] #4
@ ./Untitled-1:58 [inlined]
[8] IntegralFunction
@ ~/.julia/packages/SciMLBase/aft1j/src/scimlfunctions.jl:2183 [inlined]
[9] #46
@ ~/.julia/packages/Integrals/d5Wr6/src/Integrals.jl:113 [inlined]
[10] (::HCubature.GenzMalik{…})(f::Integrals.var"#46#48"{…}, a::StaticArraysCore.SVector{…}, b::StaticArraysCore.SVector{…}, norm::typeof(LinearAlgebra.norm))
@ HCubature ~/.julia/packages/HCubature/QvyJW/src/genz-malik.jl:121
[11] hcubature_(f::Integrals.var"#46#48"{…}, a::StaticArraysCore.SVector{…}, b::StaticArraysCore.SVector{…}, norm::typeof(LinearAlgebra.norm), rtol_::Float64, atol::Float64, maxevals::Int64, initdiv::Int64)
@ HCubature ~/.julia/packages/HCubature/QvyJW/src/HCubature.jl:61
[12] hcubature_(f::Integrals.var"#46#48"{…}, a::Vector{…}, b::Vector{…}, norm::Function, rtol::Float64, atol::Float64, maxevals::Int64, initdiv::Int64)
@ HCubature ~/.julia/packages/HCubature/QvyJW/src/HCubature.jl:129
[13] hcubature
@ ~/.julia/packages/HCubature/QvyJW/src/HCubature.jl:178 [inlined]
[14] __solvebp_call(prob::SciMLBase.IntegralProblem{…}, alg::Integrals.HCubatureJL{…}, sensealg::Integrals.ReCallVJP{…}, domain::Tuple{…}, p::Vector{…}; reltol::Float64, abstol::Float64, maxiters::Int64)
@ Integrals ~/.julia/packages/Integrals/d5Wr6/src/Integrals.jl:121
[15] __solvebp_call
@ Integrals ~/.julia/packages/Integrals/d5Wr6/src/Integrals.jl:102 [inlined]
[16] #__solvebp_call#4
@ Integrals ~/.julia/packages/Integrals/d5Wr6/src/common.jl:113 [inlined]
[17] __solvebp_call
@ Integrals ~/.julia/packages/Integrals/d5Wr6/src/common.jl:112 [inlined]
[18] __solvebp
@ Integrals ~/.julia/packages/Integrals/d5Wr6/src/Integrals.jl:65 [inlined]
[19] solve!(cache::Integrals.IntegralCache{…})
@ Integrals ~/.julia/packages/Integrals/d5Wr6/src/common.jl:103
[20] solve(prob::SciMLBase.IntegralProblem{…}, alg::Integrals.HCubatureJL{…}; kwargs::@Kwargs{…})
@ Integrals ~/.julia/packages/Integrals/d5Wr6/src/common.jl:99
[21] x1_broken(jpdf::typeof(Main.Mwe.g1), x::Float64, p::Vector{Int64}, t::Float64)
@ Main.Mwe ./Untitled-1:61
[22] macro expansion
@ show.jl:1181 [inlined]
[23] top-level scope
@ Untitled-1:73
Some type information was truncated. Use `show(err)` to see complete types. |
OK, it would probably be better to call this Now, your MWE makes it clear to me that you are expecting the type of the quadrature points to match the type of the limits of integration, however, you gave mutable limits of type |
No, the ODE solver is doing as expected. If This originates though because HCubature is giving a user StaticArray values even though they never used static arrays, in which case they have to be very careful in their usage of this specific algorithm.
Yes I think that's the key here. It's a bit odd for Integrals.jl to give the user StaticArray values without any warning on it. HCubature.jl does seem to document this (though I had no idea it did this before this issue was opened 😅), but Integrals.jl doesn't document this and it's the only solver that has this behavior. So it's a bit of an interface break and quite odd. In the wrapper we could create our own cache arrays, write into them, and send it to the user, and then have a keyword argument in the algorithm type to simply send the static array on to them, but I can see why it's confusing to a user to have this as the default behavior. |
Describe the bug 🐞
HCubature breaks down when passed an ODESolution obtained wih
Trapezoid()Expected behavior
Calculate the cubature.
Minimal Reproducible Example 👇
Without MRE, we would only be able to help you to a limited extent, and attention to the issue would be limited. to know more about MRE refer to wikipedia and stackoverflow.
Error & Stacktrace⚠️
Environment (please complete the following information):
using Pkg; Pkg.status()using Pkg; Pkg.status(; mode = PKGMODE_MANIFEST)versioninfo()Additional context
See
https://discourse.julialang.org/t/why-does-my-cubature-of-odesolution-break-down-with-finite-limits/108736
Also, the problem in general seems to be non-existent when using a different cubature with
using Cubature.The text was updated successfully, but these errors were encountered: