lxvm
commented
9 months ago @homocomputeris was this resolved on discourse?
Open homocomputeris opened 9 months ago
homocomputeris
commented
9 months ago @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
Open an issue on Integrals.jl so we can discuss this more.
lxvm
commented
9 months ago 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 Vector, e.g. with prob = IntegralProblem((tail, p) -> jpdf(collect([x; tail]), p, t), ..., but the underlying problem is probably some mishandling of static array types in the function OrdinaryDiffEq.calc_W whose result disagrees with the pre-computed cache type used in OrdinaryDiffEq.update_W! that is causing the conversion error
ERROR: MethodError: Cannot `convert` an object of type LinearAlgebra.LU{Float64, Matrix{Float64}, Vector{Int64}} to an object of type StaticArrays.LU{LinearAlgebra.LowerTriangular{Float64, StaticArraysCore.SMatrix{3, 3, Float64, 9}}, LinearAlgebra.UpperTriangular{Float64, StaticArraysCore.SMatrix{3, 3, Float64, 9}}, StaticArraysCore.SVector{3, Int64}}@ChrisRackauckas this seems like an OrdinaryDiffEq.jl bug, correct?
Also, I saw that Mwe.x1 is broken as well, and there is no difference between (in)finite limits.
homocomputeris
commented
9 months ago Also, I saw that
Mwe.x1is broken as well, and there is no difference between (in)finite limits.
With this code x1 works but x1_broken doesn't.
It's probably due to rather stupid mistake AutoTsit5(Vern9()) (which was probably supposed to be AutoVern9(...) but the fact still holds: for some integration limits it works, for some it doesn't.
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.With this code
x1works butx1_brokendoesn't.
OK, it would probably be better to call this Mwe2, since there is a silent error, i.e. x1 should also break with HCubatureJL() for the same reason x1_broken does. It looks like the infinity transformation can change the type of the input point, which it shouldn't, so I will open a pr to correct that.
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 Vector and instead are getting immutable points of type SVector. I agree it would be reasonable for Integrals.jl to enforce consistent types, so it would be good to add this to the list in #215.
ChrisRackauckas
commented
9 months ago @ChrisRackauckas this seems like an OrdinaryDiffEq.jl bug, correct?
No, the ODE solver is doing as expected. If u0 is a StaticArray, then your f should be returning StaticArrays. I think we could add an earlier error throw on this that's more explicit, but the issue is that u0 is a StaticArray but f is returning Vectors, which is not a behavior that is supported and something expected to error.
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.
I agree it would be reasonable for Integrals.jl to enforce consistent types, so it would be good to add this to the list in https://github.com/SciML/Integrals.jl/issues/215.
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.