I want to find the multi-dimensional integral of a scalar-valued function with vector argument defined on a rectangle (basically, product of truncated PDFs). The argument of the function comes from solving a system of ODEs.
As the function is defined on a rectangle, it should be sufficient to integrate it in the rectangle and not from [-Inf, Inf]. However, the cubature (and the DifferentialEquations.jl solver) breaks down:
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.
I kinda understand for infinite limits, the cubature changes the variables appropriately. The thing is, whatever the limits are, my code (as far as I understand it) shouldn’t change or cast the ODE solution –whatever flavor of vector x Julia chooses– to Number (why?) or SArray (which I’m not even using).
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 = [u[2], u[1], 0.0]
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)
end