Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Problems evaluating Koopman Expectations over large parameter space #98

Open
joegilkes opened this issue Jan 12, 2023 · 16 comments
Open

Problems evaluating Koopman Expectations over large parameter space #98

joegilkes opened this issue Jan 12, 2023 · 16 comments

Comments

@joegilkes
Copy link

Hi, first time using the package so I apologise if I'm missing something obvious.

I'm trying to propagate parametric uncertainty through a large, stiff system of ODEs representing a chemical reaction network (previous discussion here) using Koopman Expectations, such that I have a normal distribution over each of my ODESystem parameters and I am able to observe the mean and standard deviation over my state vector at every timestep specified by saveat. Here is my setup:

# `rates` represents reaction rate constants, the parameters of the model.
rates, rate_uncertainties = calc_rate_constants(...)
p_dist = [truncated(Normal(r, ru), r-(4*u), r+(4*u)) for (r, u) in zip(rates, rate_uncertainties)]

# `ODEProblem` constructed via Catalyst.jl, shouldn't be relevant.
oprob = ODEProblem(...)

# Want mean and variance for all saved timesteps.
g(sol) = [sol[1,:]; sol[1,:].^2]
gd = GenericDistribution(p_dist...)
# I believe this is the correct expression for `h(x, u, p)` when only applying uncertainty to parameters?
h(x, u, p) = u, x
sm = SystemMap(oprob, Rosenbrock23(); abstol=1e-14, reltol=1e-14, dtmin=1e-16, saveat=0.0:1e6:2e7)
exprob = ExpectationProblem(sm, g, h, gd; nout=2)
sol = solve(exprob, Koopman())

As far as I can tell this is set up correctly, but please correct me if not.

When I run this, it gets into the solve() line and then errors with the following:

ERROR: LoadError: OutOfMemoryError()
Stacktrace:
  [1] Array
    @ ./boot.jl:459 [inlined]
  [2] signcombos(k::Int64, λ::Float64, #unused#::Val{36})
    @ HCubature ~/.julia/packages/HCubature/QvyJW/src/genz-malik.jl:34
  [3] HCubature.GenzMalik(v::Val{36}, ::Type{Float64})
    @ HCubature ~/.julia/packages/HCubature/QvyJW/src/genz-malik.jl:99
  [4] cubrule
    @ ~/.julia/packages/HCubature/QvyJW/src/HCubature.jl:37 [inlined]
  [5] hcubature_(f::Integrals.var"#19#21"{IntegralProblem{false, RecursiveArrayTools.ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}, SciMLExpectations.var"#14#15"{SciMLExpectations.GenericDistribution{SciMLExpectations.var"#pdf_func#4"{NTuple{36, Distributions.Truncated{Distributions.Normal{Float64}, Distributions.Continuous, Float64}}}, SciMLExpectations.var"#rand_func#6"{NTuple{36, Distributions.Truncated{Distributions.Normal{Float64}, Distributions.Continuous, Float64}}}, StaticArraysCore.SVector{36, Float64}, StaticArraysCore.SVector{36, Float64}}, AutoBreakdown.var"#h#45", AutoBreakdown.var"#g#44", SciMLExpectations.SystemMap{ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, SciMLBase.AutoSpecialize, ModelingToolkit.var"#f#487"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x11954d43, 0xfeb7cb16, 0x36117848, 0xb62dc5d6, 0x86beea97)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xf7d03d71, 0x4a56f8e0, 0x6b4c3b18, 0xa7e9f6b2, 0x184e3749)}}, UniformScaling{Bool}, Nothing, Nothing, ModelingToolkit.var"#_jac#491"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xeb9be0d9, 0x805797f6, 0x35efa81e, 0xf7ba3bad, 0x26a49a9b)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x6e53ff16, 0x22f787c8, 0xc15faf6f, 0x2bf5f454, 0x2964e826)}}, Nothing, Nothing, SparseArrays.SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#503#generated_observed#495"{Bool, ModelingToolkit.ODESystem, Dict{Any, Any}}, Nothing, ModelingToolkit.ODESystem}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tuple{Rosenbrock23{0, true, Nothing, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}}, Base.Pairs{Symbol, Any, NTuple{7, Symbol}, NamedTuple{(:saveat, :reltol, :progress_steps, :dtmin, :abstol, :progress, :kwargshandle), Tuple{Vector{Float64}, Float64, Int64, Float64, Float64, Bool, DataType}}}}}, StaticArraysCore.SVector{36, Float64}, StaticArraysCore.SVector{36, Float64}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}}, a::StaticArraysCore.SVector{36, Float64}, b::StaticArraysCore.SVector{36, Float64}, norm::typeof(norm), rtol_::Float64, atol::Float64, maxevals::Int64, initdiv::Int64)
    @ HCubature ~/.julia/packages/HCubature/QvyJW/src/HCubature.jl:56
  [6] #hcubature#3
    @ ~/.julia/packages/HCubature/QvyJW/src/HCubature.jl:179 [inlined]
  [7] __solvebp_call(::IntegralProblem{false, RecursiveArrayTools.ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}, SciMLExpectations.var"#14#15"{SciMLExpectations.GenericDistribution{SciMLExpectations.var"#pdf_func#4"{NTuple{36, Distributions.Truncated{Distributions.Normal{Float64}, Distributions.Continuous, Float64}}}, SciMLExpectations.var"#rand_func#6"{NTuple{36, Distributions.Truncated{Distributions.Normal{Float64}, Distributions.Continuous, Float64}}}, StaticArraysCore.SVector{36, Float64}, StaticArraysCore.SVector{36, Float64}}, AutoBreakdown.var"#h#45", AutoBreakdown.var"#g#44", SciMLExpectations.SystemMap{ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, SciMLBase.AutoSpecialize, ModelingToolkit.var"#f#487"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x11954d43, 0xfeb7cb16, 0x36117848, 0xb62dc5d6, 0x86beea97)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xf7d03d71, 0x4a56f8e0, 0x6b4c3b18, 0xa7e9f6b2, 0x184e3749)}}, UniformScaling{Bool}, Nothing, Nothing, ModelingToolkit.var"#_jac#491"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xeb9be0d9, 0x805797f6, 0x35efa81e, 0xf7ba3bad, 0x26a49a9b)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x6e53ff16, 0x22f787c8, 0xc15faf6f, 0x2bf5f454, 0x2964e826)}}, Nothing, Nothing, SparseArrays.SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#503#generated_observed#495"{Bool, ModelingToolkit.ODESystem, Dict{Any, Any}}, Nothing, ModelingToolkit.ODESystem}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tuple{Rosenbrock23{0, true, Nothing, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}}, Base.Pairs{Symbol, Any, NTuple{7, Symbol}, NamedTuple{(:saveat, :reltol, :progress_steps, :dtmin, :abstol, :progress, :kwargshandle), Tuple{Vector{Float64}, Float64, Int64, Float64, Float64, Bool, DataType}}}}}, StaticArraysCore.SVector{36, Float64}, StaticArraysCore.SVector{36, Float64}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, ::Integrals.HCubatureJL, ::Integrals.ReCallVJP{Integrals.ZygoteVJP}, ::StaticArraysCore.SVector{36, Float64}, ::StaticArraysCore.SVector{36, Float64}, ::RecursiveArrayTools.ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}; reltol::Float64, abstol::Float64, maxiters::Int64, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ Integrals ~/.julia/packages/Integrals/KYdjn/src/Integrals.jl:260
  [8] #__solvebp#13
    @ ~/.julia/packages/Integrals/KYdjn/src/Integrals.jl:221 [inlined]
  [9] #solve#12
    @ ~/.julia/packages/Integrals/KYdjn/src/Integrals.jl:217 [inlined]
 [10] integrate(quadalg::Integrals.HCubatureJL, adalg::SciMLExpectations.NonfusedAD, f::Function, lb::StaticArraysCore.SVector{36, Float64}, ub::StaticArraysCore.SVector{36, Float64}, p::RecursiveArrayTools.ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}; nout::Int64, batch::Int64, kwargs::Base.Pairs{Symbol, Real, Tuple{Symbol, Symbol, Symbol}, NamedTuple{(:reltol, :abstol, :maxiters), Tuple{Float64, Float64, Int64}}})
    @ SciMLExpectations ~/.julia/packages/SciMLExpectations/DWytJ/src/expectation.jl:148
 [11] solve(::SciMLExpectations.ExpectationProblem{SciMLExpectations.SystemMap{ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, SciMLBase.AutoSpecialize, ModelingToolkit.var"#f#487"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x11954d43, 0xfeb7cb16, 0x36117848, 0xb62dc5d6, 0x86beea97)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xf7d03d71, 0x4a56f8e0, 0x6b4c3b18, 0xa7e9f6b2, 0x184e3749)}}, UniformScaling{Bool}, Nothing, Nothing, ModelingToolkit.var"#_jac#491"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xeb9be0d9, 0x805797f6, 0x35efa81e, 0xf7ba3bad, 0x26a49a9b)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x6e53ff16, 0x22f787c8, 0xc15faf6f, 0x2bf5f454, 0x2964e826)}}, Nothing, Nothing, SparseArrays.SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#503#generated_observed#495"{Bool, ModelingToolkit.ODESystem, Dict{Any, Any}}, Nothing, ModelingToolkit.ODESystem}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tuple{Rosenbrock23{0, true, Nothing, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}}, Base.Pairs{Symbol, Any, NTuple{7, Symbol}, NamedTuple{(:saveat, :reltol, :progress_steps, :dtmin, :abstol, :progress, :kwargshandle), Tuple{Vector{Float64}, Float64, Int64, Float64, Float64, Bool, DataType}}}}, AutoBreakdown.var"#g#44", AutoBreakdown.var"#h#45", SciMLExpectations.GenericDistribution{SciMLExpectations.var"#pdf_func#4"{NTuple{36, Distributions.Truncated{Distributions.Normal{Float64}, Distributions.Continuous, Float64}}}, SciMLExpectations.var"#rand_func#6"{NTuple{36, Distributions.Truncated{Distributions.Normal{Float64}, Distributions.Continuous, Float64}}}, StaticArraysCore.SVector{36, Float64}, StaticArraysCore.SVector{36, Float64}}, RecursiveArrayTools.ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}}, ::SciMLExpectations.Koopman{SciMLExpectations.NonfusedAD}; maxiters::Int64, batch::Int64, quadalg::Integrals.HCubatureJL, ireltol::Float64, iabstol::Float64, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ SciMLExpectations ~/.julia/packages/SciMLExpectations/DWytJ/src/expectation.jl:134
 [12] solve(::SciMLExpectations.ExpectationProblem{SciMLExpectations.SystemMap{ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, SciMLBase.AutoSpecialize, ModelingToolkit.var"#f#487"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x11954d43, 0xfeb7cb16, 0x36117848, 0xb62dc5d6, 0x86beea97)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xf7d03d71, 0x4a56f8e0, 0x6b4c3b18, 0xa7e9f6b2, 0x184e3749)}}, UniformScaling{Bool}, Nothing, Nothing, ModelingToolkit.var"#_jac#491"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xeb9be0d9, 0x805797f6, 0x35efa81e, 0xf7ba3bad, 0x26a49a9b)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x6e53ff16, 0x22f787c8, 0xc15faf6f, 0x2bf5f454, 0x2964e826)}}, Nothing, Nothing, SparseArrays.SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#503#generated_observed#495"{Bool, ModelingToolkit.ODESystem, Dict{Any, Any}}, Nothing, ModelingToolkit.ODESystem}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tuple{Rosenbrock23{0, true, Nothing, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}}, Base.Pairs{Symbol, Any, NTuple{7, Symbol}, NamedTuple{(:saveat, :reltol, :progress_steps, :dtmin, :abstol, :progress, :kwargshandle), Tuple{Vector{Float64}, Float64, Int64, Float64, Float64, Bool, DataType}}}}, AutoBreakdown.var"#g#44", AutoBreakdown.var"#h#45", SciMLExpectations.GenericDistribution{SciMLExpectations.var"#pdf_func#4"{NTuple{36, Distributions.Truncated{Distributions.Normal{Float64}, Distributions.Continuous, Float64}}}, SciMLExpectations.var"#rand_func#6"{NTuple{36, Distributions.Truncated{Distributions.Normal{Float64}, Distributions.Continuous, Float64}}}, StaticArraysCore.SVector{36, Float64}, StaticArraysCore.SVector{36, Float64}}, RecursiveArrayTools.ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}}, ::SciMLExpectations.Koopman{SciMLExpectations.NonfusedAD})
    @ SciMLExpectations ~/.julia/packages/SciMLExpectations/DWytJ/src/expectation.jl:125

Watching my system memory usage at the time indicates that I really shouldn't be running out of memory (~7 GB / 32 GB used when it errors), so I'm really not sure what to make of the issue. I'd like this to work on systems of >2000 ODES with > 7000 parameters, but right now this crash is happening on a much smaller system of 20 ODEs with only 36 parameters (hence the Rosenbrock23 solver, I usually use QNDF).
@agerlach
Copy link
Contributor

I'm not sure if this is driving the memory issue, but the default quadrature alg is probably only good up to 10 dimensions. You can set quadalg kwarg during solve to anything supported in Integrals.jl. For the size system you are interested in you will have to rely on the Cuba quasi-mc integration methods, in particular I would try CubaDivonne and CubaSUAVE.

@agerlach
Copy link
Contributor

Your h function looks correct.

@joegilkes
Copy link
Author

Ah thanks, that's good to know. Changing from the default quadalg seems to get rid of the memory issue, but both of the Cuba quadratures seem to have independent issues with the dimensionality of the system.

Changing to CubaDivonne yields

ERROR: LoadError: ArgumentError: In Divonne ndim must be ≥ 2
Stacktrace:
 [1] divonne(integrand::IntegralsCuba.var"#4#20"{IntegralProblem{false, RecursiveArrayTools.ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}, SciMLExpectations.var"#14#15"{SciMLExpectations.GenericDistribution{SciMLExpectations.var"#pdf_func#4"{Tuple{Distributions.Truncated{Distributions.Normal{Float64}, Distributions.Continuous, Float64}}}, SciMLExpectations.var"#rand_func#6"{Tuple{Distributions.Truncated{Distributions.Normal{Float64}, Distributions.Continuous, Float64}}}, StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SVector{1, Float64}}, AutoBreakdown.var"#h#49", AutoBreakdown.var"#g#48", SciMLExpectations.SystemMap{ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, SciMLBase.AutoSpecialize, ModelingToolkit.var"#f#487"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x11954d43, 0xfeb7cb16, 0x36117848, 0xb62dc5d6, 0x86beea97)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xf7d03d71, 0x4a56f8e0, 0x6b4c3b18, 0xa7e9f6b2, 0x184e3749)}}, UniformScaling{Bool}, Nothing, Nothing, ModelingToolkit.var"#_jac#491"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xeb9be0d9, 0x805797f6, 0x35efa81e, 0xf7ba3bad, 0x26a49a9b)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x6e53ff16, 0x22f787c8, 0xc15faf6f, 0x2bf5f454, 0x2964e826)}}, Nothing, Nothing, SparseArrays.SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#503#generated_observed#495"{Bool, ModelingToolkit.ODESystem, Dict{Any, Any}}, Nothing, ModelingToolkit.ODESystem}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tuple{Rosenbrock23{0, true, Nothing, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}}, Base.Pairs{Symbol, Any, NTuple{7, Symbol}, NamedTuple{(:saveat, :reltol, :progress_steps, :dtmin, :abstol, :progress, :kwargshandle), Tuple{Vector{Float64}, Float64, Int64, Float64, Float64, Bool, DataType}}}}}, StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SVector{1, Float64}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SVector{1, Float64}, RecursiveArrayTools.ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}}, ndim::Int64, ncomp::Int64; nvec::Int64, rtol::Float64, atol::Float64, flags::Int64, seed::Int64, minevals::Int64, maxevals::Int64, key1::Int64, key2::Int64, key3::Int64, maxpass::Int64, border::Float64, maxchisq::Float64, mindeviation::Float64, ngiven::Int64, ldxgiven::Int64, xgiven::Matrix{Float64}, nextra::Int64, peakfinder::Ptr{Nothing}, statefile::String, spin::Ptr{Nothing}, reltol::Missing, abstol::Missing, userdata::Missing)
   @ Cuba ~/.julia/packages/Cuba/aaNMS/src/divonne.jl:129
 [2] __solvebp_call(::IntegralProblem{false, RecursiveArrayTools.ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}, SciMLExpectations.var"#14#15"{SciMLExpectations.GenericDistribution{SciMLExpectations.var"#pdf_func#4"{Tuple{Distributions.Truncated{Distributions.Normal{Float64}, Distributions.Continuous, Float64}}}, SciMLExpectations.var"#rand_func#6"{Tuple{Distributions.Truncated{Distributions.Normal{Float64}, Distributions.Continuous, Float64}}}, StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SVector{1, Float64}}, AutoBreakdown.var"#h#49", AutoBreakdown.var"#g#48", SciMLExpectations.SystemMap{ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, SciMLBase.AutoSpecialize, ModelingToolkit.var"#f#487"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x11954d43, 0xfeb7cb16, 0x36117848, 0xb62dc5d6, 0x86beea97)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xf7d03d71, 0x4a56f8e0, 0x6b4c3b18, 0xa7e9f6b2, 0x184e3749)}}, UniformScaling{Bool}, Nothing, Nothing, ModelingToolkit.var"#_jac#491"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xeb9be0d9, 0x805797f6, 0x35efa81e, 0xf7ba3bad, 0x26a49a9b)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x6e53ff16, 0x22f787c8, 0xc15faf6f, 0x2bf5f454, 0x2964e826)}}, Nothing, Nothing, SparseArrays.SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#503#generated_observed#495"{Bool, ModelingToolkit.ODESystem, Dict{Any, Any}}, Nothing, ModelingToolkit.ODESystem}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tuple{Rosenbrock23{0, true, Nothing, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}}, Base.Pairs{Symbol, Any, NTuple{7, Symbol}, NamedTuple{(:saveat, :reltol, :progress_steps, :dtmin, :abstol, :progress, :kwargshandle), Tuple{Vector{Float64}, Float64, Int64, Float64, Float64, Bool, DataType}}}}}, StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SVector{1, Float64}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, ::IntegralsCuba.CubaDivonne, ::Integrals.ReCallVJP{Integrals.ZygoteVJP}, ::StaticArraysCore.SVector{1, Float64}, ::StaticArraysCore.SVector{1, Float64}, ::RecursiveArrayTools.ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}; reltol::Float64, abstol::Float64, maxiters::Int64, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
   @ IntegralsCuba ~/.julia/packages/IntegralsCuba/CPNxI/src/IntegralsCuba.jl:100
 [3] #__solvebp#13
   @ ~/.julia/packages/Integrals/KYdjn/src/Integrals.jl:221 [inlined]
 [4] #solve#12
   @ ~/.julia/packages/Integrals/KYdjn/src/Integrals.jl:217 [inlined]
 [5] integrate(quadalg::IntegralsCuba.CubaDivonne, adalg::SciMLExpectations.NonfusedAD, f::Function, lb::StaticArraysCore.SVector{1, Float64}, ub::StaticArraysCore.SVector{1, Float64}, p::RecursiveArrayTools.ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}; nout::Int64, batch::Int64, kwargs::Base.Pairs{Symbol, Real, Tuple{Symbol, Symbol, Symbol}, NamedTuple{(:reltol, :abstol, :maxiters), Tuple{Float64, Float64, Int64}}})
   @ SciMLExpectations ~/.julia/packages/SciMLExpectations/DWytJ/src/expectation.jl:148
 [6] solve(::SciMLExpectations.ExpectationProblem{SciMLExpectations.SystemMap{ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, SciMLBase.AutoSpecialize, ModelingToolkit.var"#f#487"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x11954d43, 0xfeb7cb16, 0x36117848, 0xb62dc5d6, 0x86beea97)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xf7d03d71, 0x4a56f8e0, 0x6b4c3b18, 0xa7e9f6b2, 0x184e3749)}}, UniformScaling{Bool}, Nothing, Nothing, ModelingToolkit.var"#_jac#491"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xeb9be0d9, 0x805797f6, 0x35efa81e, 0xf7ba3bad, 0x26a49a9b)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x6e53ff16, 0x22f787c8, 0xc15faf6f, 0x2bf5f454, 0x2964e826)}}, Nothing, Nothing, SparseArrays.SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#503#generated_observed#495"{Bool, ModelingToolkit.ODESystem, Dict{Any, Any}}, Nothing, ModelingToolkit.ODESystem}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tuple{Rosenbrock23{0, true, Nothing, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}}, Base.Pairs{Symbol, Any, NTuple{7, Symbol}, NamedTuple{(:saveat, :reltol, :progress_steps, :dtmin, :abstol, :progress, :kwargshandle), Tuple{Vector{Float64}, Float64, Int64, Float64, Float64, Bool, DataType}}}}, AutoBreakdown.var"#g#48", AutoBreakdown.var"#h#49", SciMLExpectations.GenericDistribution{SciMLExpectations.var"#pdf_func#4"{Tuple{Distributions.Truncated{Distributions.Normal{Float64}, Distributions.Continuous, Float64}}}, SciMLExpectations.var"#rand_func#6"{Tuple{Distributions.Truncated{Distributions.Normal{Float64}, Distributions.Continuous, Float64}}}, StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SVector{1, Float64}}, RecursiveArrayTools.ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}}, ::SciMLExpectations.Koopman{SciMLExpectations.NonfusedAD}; maxiters::Int64, batch::Int64, quadalg::IntegralsCuba.CubaDivonne, ireltol::Float64, iabstol::Float64, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
   @ SciMLExpectations ~/.julia/packages/SciMLExpectations/DWytJ/src/expectation.jl:134

Meanwhile, changing to CubaSUAVE yields

ERROR: LoadError: DimensionMismatch: array could not be broadcast to match destination
Stacktrace:
  [1] check_broadcast_shape
    @ ./broadcast.jl:540 [inlined]
  [2] check_broadcast_axes
    @ ./broadcast.jl:543 [inlined]
  [3] check_broadcast_axes
    @ ./broadcast.jl:546 [inlined]
  [4] instantiate
    @ ./broadcast.jl:284 [inlined]
  [5] materialize!
    @ ./broadcast.jl:871 [inlined]
  [6] materialize!(dest::Vector{Float64}, bc::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1}, Nothing, typeof(*), Tuple{Vector{Float64}, Float64}})
    @ Base.Broadcast ./broadcast.jl:868
  [7] (::IntegralsCuba.var"#4#20"{IntegralProblem{false, RecursiveArrayTools.ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}, SciMLExpectations.var"#14#15"{SciMLExpectations.GenericDistribution{SciMLExpectations.var"#pdf_func#4"{Tuple{Distributions.Truncated{Distributions.Normal{Float64}, Distributions.Continuous, Float64}}}, SciMLExpectations.var"#rand_func#6"{Tuple{Distributions.Truncated{Distributions.Normal{Float64}, Distributions.Continuous, Float64}}}, StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SVector{1, Float64}}, AutoBreakdown.var"#h#49", AutoBreakdown.var"#g#48", SciMLExpectations.SystemMap{ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, SciMLBase.AutoSpecialize, ModelingToolkit.var"#f#487"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x11954d43, 0xfeb7cb16, 0x36117848, 0xb62dc5d6, 0x86beea97)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xf7d03d71, 0x4a56f8e0, 0x6b4c3b18, 0xa7e9f6b2, 0x184e3749)}}, UniformScaling{Bool}, Nothing, Nothing, ModelingToolkit.var"#_jac#491"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xeb9be0d9, 0x805797f6, 0x35efa81e, 0xf7ba3bad, 0x26a49a9b)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x6e53ff16, 0x22f787c8, 0xc15faf6f, 0x2bf5f454, 0x2964e826)}}, Nothing, Nothing, SparseArrays.SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#503#generated_observed#495"{Bool, ModelingToolkit.ODESystem, Dict{Any, Any}}, Nothing, ModelingToolkit.ODESystem}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tuple{Rosenbrock23{0, true, Nothing, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}}, Base.Pairs{Symbol, Any, NTuple{7, Symbol}, NamedTuple{(:saveat, :reltol, :progress_steps, :dtmin, :abstol, :progress, :kwargshandle), Tuple{Vector{Float64}, Float64, Int64, Float64, Float64, Bool, DataType}}}}}, StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SVector{1, Float64}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SVector{1, Float64}, RecursiveArrayTools.ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}})(x::Vector{Float64}, dx::Vector{Float64})
    @ IntegralsCuba ~/.julia/packages/IntegralsCuba/CPNxI/src/IntegralsCuba.jl:40
  [8] generic_integrand!(ndim::Int32, x_::Ptr{Float64}, ncomp::Int32, f_::Ptr{Float64}, func!::IntegralsCuba.var"#4#20"{IntegralProblem{false, RecursiveArrayTools.ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}, SciMLExpectations.var"#14#15"{SciMLExpectations.GenericDistribution{SciMLExpectations.var"#pdf_func#4"{Tuple{Distributions.Truncated{Distributions.Normal{Float64}, Distributions.Continuous, Float64}}}, SciMLExpectations.var"#rand_func#6"{Tuple{Distributions.Truncated{Distributions.Normal{Float64}, Distributions.Continuous, Float64}}}, StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SVector{1, Float64}}, AutoBreakdown.var"#h#49", AutoBreakdown.var"#g#48", SciMLExpectations.SystemMap{ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, SciMLBase.AutoSpecialize, ModelingToolkit.var"#f#487"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x11954d43, 0xfeb7cb16, 0x36117848, 0xb62dc5d6, 0x86beea97)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xf7d03d71, 0x4a56f8e0, 0x6b4c3b18, 0xa7e9f6b2, 0x184e3749)}}, UniformScaling{Bool}, Nothing, Nothing, ModelingToolkit.var"#_jac#491"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xeb9be0d9, 0x805797f6, 0x35efa81e, 0xf7ba3bad, 0x26a49a9b)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x6e53ff16, 0x22f787c8, 0xc15faf6f, 0x2bf5f454, 0x2964e826)}}, Nothing, Nothing, SparseArrays.SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#503#generated_observed#495"{Bool, ModelingToolkit.ODESystem, Dict{Any, Any}}, Nothing, ModelingToolkit.ODESystem}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tuple{Rosenbrock23{0, true, Nothing, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}}, Base.Pairs{Symbol, Any, NTuple{7, Symbol}, NamedTuple{(:saveat, :reltol, :progress_steps, :dtmin, :abstol, :progress, :kwargshandle), Tuple{Vector{Float64}, Float64, Int64, Float64, Float64, Bool, DataType}}}}}, StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SVector{1, Float64}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SVector{1, Float64}, RecursiveArrayTools.ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}})
    @ Cuba ~/.julia/packages/Cuba/aaNMS/src/Cuba.jl:72
  [9] dointegrate!
    @ ~/.julia/packages/Cuba/aaNMS/src/suave.jl:27 [inlined]
 [10] dointegrate
    @ ~/.julia/packages/Cuba/aaNMS/src/Cuba.jl:209 [inlined]
 [11] suave(integrand::IntegralsCuba.var"#4#20"{IntegralProblem{false, RecursiveArrayTools.ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}, SciMLExpectations.var"#14#15"{SciMLExpectations.GenericDistribution{SciMLExpectations.var"#pdf_func#4"{Tuple{Distributions.Truncated{Distributions.Normal{Float64}, Distributions.Continuous, Float64}}}, SciMLExpectations.var"#rand_func#6"{Tuple{Distributions.Truncated{Distributions.Normal{Float64}, Distributions.Continuous, Float64}}}, StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SVector{1, Float64}}, AutoBreakdown.var"#h#49", AutoBreakdown.var"#g#48", SciMLExpectations.SystemMap{ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, SciMLBase.AutoSpecialize, ModelingToolkit.var"#f#487"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x11954d43, 0xfeb7cb16, 0x36117848, 0xb62dc5d6, 0x86beea97)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xf7d03d71, 0x4a56f8e0, 0x6b4c3b18, 0xa7e9f6b2, 0x184e3749)}}, UniformScaling{Bool}, Nothing, Nothing, ModelingToolkit.var"#_jac#491"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xeb9be0d9, 0x805797f6, 0x35efa81e, 0xf7ba3bad, 0x26a49a9b)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x6e53ff16, 0x22f787c8, 0xc15faf6f, 0x2bf5f454, 0x2964e826)}}, Nothing, Nothing, SparseArrays.SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#503#generated_observed#495"{Bool, ModelingToolkit.ODESystem, Dict{Any, Any}}, Nothing, ModelingToolkit.ODESystem}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tuple{Rosenbrock23{0, true, Nothing, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}}, Base.Pairs{Symbol, Any, NTuple{7, Symbol}, NamedTuple{(:saveat, :reltol, :progress_steps, :dtmin, :abstol, :progress, :kwargshandle), Tuple{Vector{Float64}, Float64, Int64, Float64, Float64, Bool, DataType}}}}}, StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SVector{1, Float64}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SVector{1, Float64}, RecursiveArrayTools.ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}}, ndim::Int64, ncomp::Int64; nvec::Int64, rtol::Float64, atol::Float64, flags::Int64, seed::Int64, minevals::Int64, maxevals::Int64, nnew::Int64, nmin::Int64, flatness::Float64, statefile::String, spin::Ptr{Nothing}, reltol::Missing, abstol::Missing, userdata::Missing)
    @ Cuba ~/.julia/packages/Cuba/aaNMS/src/suave.jl:89
 [12] __solvebp_call(::IntegralProblem{false, RecursiveArrayTools.ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}, SciMLExpectations.var"#14#15"{SciMLExpectations.GenericDistribution{SciMLExpectations.var"#pdf_func#4"{Tuple{Distributions.Truncated{Distributions.Normal{Float64}, Distributions.Continuous, Float64}}}, SciMLExpectations.var"#rand_func#6"{Tuple{Distributions.Truncated{Distributions.Normal{Float64}, Distributions.Continuous, Float64}}}, StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SVector{1, Float64}}, AutoBreakdown.var"#h#49", AutoBreakdown.var"#g#48", SciMLExpectations.SystemMap{ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, SciMLBase.AutoSpecialize, ModelingToolkit.var"#f#487"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x11954d43, 0xfeb7cb16, 0x36117848, 0xb62dc5d6, 0x86beea97)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xf7d03d71, 0x4a56f8e0, 0x6b4c3b18, 0xa7e9f6b2, 0x184e3749)}}, UniformScaling{Bool}, Nothing, Nothing, ModelingToolkit.var"#_jac#491"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xeb9be0d9, 0x805797f6, 0x35efa81e, 0xf7ba3bad, 0x26a49a9b)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x6e53ff16, 0x22f787c8, 0xc15faf6f, 0x2bf5f454, 0x2964e826)}}, Nothing, Nothing, SparseArrays.SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#503#generated_observed#495"{Bool, ModelingToolkit.ODESystem, Dict{Any, Any}}, Nothing, ModelingToolkit.ODESystem}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tuple{Rosenbrock23{0, true, Nothing, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}}, Base.Pairs{Symbol, Any, NTuple{7, Symbol}, NamedTuple{(:saveat, :reltol, :progress_steps, :dtmin, :abstol, :progress, :kwargshandle), Tuple{Vector{Float64}, Float64, Int64, Float64, Float64, Bool, DataType}}}}}, StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SVector{1, Float64}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, ::IntegralsCuba.CubaSUAVE, ::Integrals.ReCallVJP{Integrals.ZygoteVJP}, ::StaticArraysCore.SVector{1, Float64}, ::StaticArraysCore.SVector{1, Float64}, ::RecursiveArrayTools.ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}; reltol::Float64, abstol::Float64, maxiters::Int64, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ IntegralsCuba ~/.julia/packages/IntegralsCuba/CPNxI/src/IntegralsCuba.jl:96
 [13] #__solvebp#13
    @ ~/.julia/packages/Integrals/KYdjn/src/Integrals.jl:221 [inlined]
 [14] #solve#12
    @ ~/.julia/packages/Integrals/KYdjn/src/Integrals.jl:217 [inlined]
 [15] integrate(quadalg::IntegralsCuba.CubaSUAVE, adalg::SciMLExpectations.NonfusedAD, f::Function, lb::StaticArraysCore.SVector{1, Float64}, ub::StaticArraysCore.SVector{1, Float64}, p::RecursiveArrayTools.ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}; nout::Int64, batch::Int64, kwargs::Base.Pairs{Symbol, Real, Tuple{Symbol, Symbol, Symbol}, NamedTuple{(:reltol, :abstol, :maxiters), Tuple{Float64, Float64, Int64}}})
    @ SciMLExpectations ~/.julia/packages/SciMLExpectations/DWytJ/src/expectation.jl:148
 [16] solve(::SciMLExpectations.ExpectationProblem{SciMLExpectations.SystemMap{ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, SciMLBase.AutoSpecialize, ModelingToolkit.var"#f#487"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x11954d43, 0xfeb7cb16, 0x36117848, 0xb62dc5d6, 0x86beea97)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xf7d03d71, 0x4a56f8e0, 0x6b4c3b18, 0xa7e9f6b2, 0x184e3749)}}, UniformScaling{Bool}, Nothing, Nothing, ModelingToolkit.var"#_jac#491"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xeb9be0d9, 0x805797f6, 0x35efa81e, 0xf7ba3bad, 0x26a49a9b)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x6e53ff16, 0x22f787c8, 0xc15faf6f, 0x2bf5f454, 0x2964e826)}}, Nothing, Nothing, SparseArrays.SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#503#generated_observed#495"{Bool, ModelingToolkit.ODESystem, Dict{Any, Any}}, Nothing, ModelingToolkit.ODESystem}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tuple{Rosenbrock23{0, true, Nothing, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}}, Base.Pairs{Symbol, Any, NTuple{7, Symbol}, NamedTuple{(:saveat, :reltol, :progress_steps, :dtmin, :abstol, :progress, :kwargshandle), Tuple{Vector{Float64}, Float64, Int64, Float64, Float64, Bool, DataType}}}}, AutoBreakdown.var"#g#48", AutoBreakdown.var"#h#49", SciMLExpectations.GenericDistribution{SciMLExpectations.var"#pdf_func#4"{Tuple{Distributions.Truncated{Distributions.Normal{Float64}, Distributions.Continuous, Float64}}}, SciMLExpectations.var"#rand_func#6"{Tuple{Distributions.Truncated{Distributions.Normal{Float64}, Distributions.Continuous, Float64}}}, StaticArraysCore.SVector{1, Float64}, StaticArraysCore.SVector{1, Float64}}, RecursiveArrayTools.ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}}, ::SciMLExpectations.Koopman{SciMLExpectations.NonfusedAD}; maxiters::Int64, batch::Int64, quadalg::IntegralsCuba.CubaSUAVE, ireltol::Float64, iabstol::Float64, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ SciMLExpectations ~/.julia/packages/SciMLExpectations/DWytJ/src/expectation.jl:134

Could this be something to do with only having uncertainty over the parameters, but none over the initial conditions?

@agerlach
Copy link
Contributor

I think your nout is wrong. It should actually be 2*length(0.0:1e6:2e7) since you are computing the mean and 2nd raw moment at each time point. Before trying that though, try something simpler, like nout=2 and

g(sol) = [sol[1,end]; sol[1,end]^2]

@joegilkes
Copy link
Author

joegilkes commented Jan 12, 2023 •
edited
Loading

I seem to have it running now, I actually had an issue in my calc_rate_constants() function that was causing p_dist to be formed incorrectly, hence the issues with dimensionality from CubaDivonne. Thanks for the simplification idea, it helped narrow things down.

I'm now trying to get the full observable out, so the mean and 2nd moment for every state at every saved timestep. Is it valid to do something like g(sol, p) = [sol[:, :], sol[:, :].^2] with nout = 2, making the two outputs the individual arrays, or do the solution arrays need to be fully unrolled with nout = 2 * length(0.0:1e6:2e7) * n_states so that every observable is individually represented, and then split/reshaped to obtain the desired arrays afterwards?

@agerlach
Copy link
Contributor

I think you will need to unroll. I doubt Integrals.jl can support what you are asking. You could possible try some fancy vector type like ComponentArrays.jl to keep everything partitioned.

I would test that using Integrals.jl directly with a MWE

@joegilkes
Copy link
Author

joegilkes commented Jan 13, 2023 •
edited
Loading

Thanks, I've got it working now by doing

function g(sol, p)
    solvec = reduce(vcat, sol.u)
    return [solvec; solvec.^2]
end

to concatenate all the solution values together, and then

nout = 2 * length(saveat) * n_states

after which I just split and reshape the resulting sol.u. However, I'm not getting the expected results. Using 'fake' uncertainties such that in each normal distribution $\sigma = \frac{\mu}{100}$ (just to ensure that my provided uncertainties aren't the issue), my resulting expectations are all ~ 1.6e-315, when they should be on the order of 0.1. Even reducing to $\sigma = \frac{\mu}{1e8}$ so that the distributions are essentially just sampling the mean values should return the original solution, but does pretty much the same, with the $\sigma$ resulting from the 2nd central moment sitting around 4e-158 and occasionally bouncing up to 1e153.

I've attached the full trace of the results to illustrate (top), as well as the trace that I should be getting (bottom, as a result of doing a regular solve on just the ODEProblem).

expectation_test_1000K
rads_small_certain_1000K

The above was with CubaSUAVE, although CubaDivonne does pretty similarly. Could this be a problem with e.g. integral solver tolerances? Happy to provide a MWE if it helps, I'll just need to package up the ODEProblem since it's a pain to recreate without providing all of the files behind it as well.

@joegilkes
Copy link
Author

It's also probably worth noting that I'm pinning the truncation of my parameters' normal distributions to be positive by doing

p_dist = [truncated(Normal(r, ru), max(r-(4*u), 0.0), r+(4*u)) for (r, u) in zip(rates, rate_uncertainties)]

as none of these rate constants should be negative, but this shouldn't be affecting things with the $\sigma = \frac{\mu}{1e8}$ example at all since the width of the distribution is so many orders of magnitude below its mean.

@joegilkes joegilkes changed the title OutOfMemoryError when constructing cubature Problems evaluating Koopman Expectations over large parameter space Jan 13, 2023
@agerlach
Copy link
Contributor

agerlach commented Jan 13, 2023 •
edited
Loading

I will admit that I've never tried this technique on a system with both a large number of uncertain dimensions AND a large observable vector.

It is most likely a matter of tuning the integration parameters. Note that when you consider a vector valued observable like you are the tolerencing is with respect to some norm of the observation vector. Some of the integration methods allow you to change the norm and others don't. Intergrals.jl just wraps other integration packages into a common interface. You may want to check the upstream docs.

Just to try to isolate the issue I would try something like form an integrand of nout = 2 * length(saveat) * n_states normal pdfs and try integrating it directly using Integrals.jl. It should give 1 everywhere, but it may help you get a feel for the play between the norm, tolerances, and other integration parameters

@joegilkes
Copy link
Author

Okay thanks, will give that a go. Just to check, do the abstol and reltol that go into SystemMap (and therefore go into the ODE solver call) also affect the reltol and abstol of the integration method, or can these be set separately?

@agerlach
Copy link
Contributor

This is exactly why SystemMap exists; to prevent kwarg collision. So, ODE kwargs to to SystemMap and Integral.jl kwargs go to solve

@agerlach
Copy link
Contributor

agerlach commented Jan 13, 2023 •
edited
Loading 

solve(exprob, Koopman(); quadalg=CubatureJLh(),
    ireltol=1e-3, iabstol=1e-3)

For setting spatial integration tols. Note, it looks like they are 1e-2 by default

SciMLExpectations.jl/src/expectation.jl

Line 137 in 4e8cbaf

function DiffEqBase.solve(prob::ExpectationProblem, expalg::Koopman, args...;

@joegilkes
Copy link
Author

That's exactly what I was looking for, thanks. I'll have a play around with those and see if I can get it working properly.

@joegilkes
Copy link
Author

I'm still testing things out, but just a point of note. My initial condition u0 = zeros(Float64, n_states); u0[1] = 1.0 is non-zero for the first state in my observable state vector, and yet with the problematic solves above, the expectation of this value at t = 0.0 always equals ~ 1.6e-315. Further, in the final expectation matrix it seems like there's a pattern:

200×28 adjoint(::Matrix{Float64}) with eltype Float64:
 1.56768e-315  5.0e-324  1.56768e-315  5.0e-324  1.56766e-315  5.0e-324  1.56766e-315  5.0e-324  1.56766e-315  5.0e-324  1.56766e-315  …  1.56765e-315  5.0e-324  1.56764e-315  5.0e-324  1.56764e-315  5.0e-324  1.56764e-315  5.0e-324  1.56763e-315  5.0e-324
 1.56763e-315  5.0e-324  1.56763e-315  5.0e-324  1.56763e-315  5.0e-324  1.56744e-315  5.0e-324  1.56744e-315  5.0e-324  1.56743e-315     1.56742e-315  5.0e-324  1.56742e-315  5.0e-324  1.56742e-315  5.0e-324  1.56742e-315  5.0e-324  1.56741e-315  5.0e-324
 1.56741e-315  5.0e-324  1.56741e-315  5.0e-324  1.5674e-315   5.0e-324  1.5674e-315   5.0e-324  1.5674e-315   5.0e-324  1.5674e-315      1.56739e-315  5.0e-324  1.56738e-315  5.0e-324  1.56738e-315  5.0e-324  1.56738e-315  5.0e-324  1.56737e-315  5.0e-324
 1.56737e-315  5.0e-324  1.56737e-315  5.0e-324  1.56737e-315  5.0e-324  1.56736e-315  5.0e-324  1.56736e-315  5.0e-324  1.56736e-315     1.56735e-315  5.0e-324  1.56735e-315  5.0e-324  1.56734e-315  5.0e-324  1.56734e-315  5.0e-324  1.56734e-315  5.0e-324
 1.56733e-315  5.0e-324  1.56733e-315  5.0e-324  1.56733e-315  5.0e-324  1.56733e-315  5.0e-324  1.56732e-315  5.0e-324  1.56732e-315     1.56731e-315  5.0e-324  1.56731e-315  5.0e-324  1.56731e-315  5.0e-324  1.5673e-315   5.0e-324  1.5673e-315   5.0e-324
 1.5673e-315   5.0e-324  1.56729e-315  5.0e-324  1.56729e-315  5.0e-324  1.56729e-315  5.0e-324  1.56729e-315  5.0e-324  1.56728e-315  …  1.56727e-315  5.0e-324  1.56727e-315  5.0e-324  1.56727e-315  5.0e-324  1.56726e-315  5.0e-324  1.56726e-315  5.0e-324
 1.56726e-315  5.0e-324  1.56726e-315  5.0e-324  1.56725e-315  5.0e-324  1.56725e-315  5.0e-324  1.56725e-315  5.0e-324  1.56725e-315     1.56724e-315  5.0e-324  1.56723e-315  5.0e-324  1.56723e-315  5.0e-324  1.56723e-315  5.0e-324  1.56722e-315  5.0e-324
 1.56722e-315  5.0e-324  1.56722e-315  5.0e-324  1.56722e-315  5.0e-324  1.56721e-315  5.0e-324  1.56721e-315  5.0e-324  1.56721e-315     1.5672e-315   5.0e-324  1.56719e-315  5.0e-324  1.56719e-315  5.0e-324  1.56719e-315  5.0e-324  1.56719e-315  5.0e-324
 ⋮                                                             ⋮                                                         ⋮             ⋱                          ⋮                                                             ⋮                       
 1.56714e-315  5.0e-324  1.56714e-315  5.0e-324  1.56713e-315  5.0e-324  1.56713e-315  5.0e-324  1.56713e-315  5.0e-324  1.56712e-315     1.56711e-315  5.0e-324  1.56711e-315  5.0e-324  1.56711e-315  5.0e-324  1.56711e-315  5.0e-324  1.5671e-315   5.0e-324
 1.5671e-315   5.0e-324  1.5671e-315   5.0e-324  1.5671e-315   5.0e-324  1.56709e-315  5.0e-324  1.56709e-315  5.0e-324  1.56709e-315     1.56708e-315  5.0e-324  1.56707e-315  5.0e-324  1.56707e-315  5.0e-324  1.56707e-315  5.0e-324  1.56707e-315  5.0e-324
 1.56706e-315  5.0e-324  1.56706e-315  5.0e-324  1.56706e-315  5.0e-324  1.56706e-315  5.0e-324  1.56705e-315  5.0e-324  1.56705e-315     1.56704e-315  5.0e-324  1.56704e-315  5.0e-324  1.56703e-315  5.0e-324  1.56703e-315  5.0e-324  1.56703e-315  5.0e-324
 1.56703e-315  5.0e-324  1.56702e-315  5.0e-324  1.56702e-315  5.0e-324  1.56702e-315  5.0e-324  1.56701e-315  5.0e-324  1.56701e-315  …  1.567e-315    5.0e-324  1.567e-315    5.0e-324  1.567e-315    5.0e-324  1.56699e-315  5.0e-324  1.56699e-315  5.0e-324
 1.56699e-315  5.0e-324  1.56699e-315  5.0e-324  1.56698e-315  5.0e-324  1.56698e-315  5.0e-324  1.56698e-315  5.0e-324  1.56697e-315     1.56696e-315  5.0e-324  1.56696e-315  5.0e-324  1.56696e-315  5.0e-324  1.56696e-315  5.0e-324  1.56695e-315  5.0e-324
 1.56695e-315  5.0e-324  1.56695e-315  5.0e-324  1.56694e-315  5.0e-324  1.56694e-315  5.0e-324  1.56694e-315  5.0e-324  1.56694e-315     1.56693e-315  5.0e-324  1.56692e-315  5.0e-324  1.56692e-315  5.0e-324  1.56692e-315  5.0e-324  1.56692e-315  5.0e-324
 1.56691e-315  5.0e-324  1.56691e-315  5.0e-324  1.56691e-315  5.0e-324  1.5669e-315   5.0e-324  1.5669e-315   5.0e-324  1.5669e-315      1.56689e-315  5.0e-324  1.56689e-315  5.0e-324  1.56688e-315  5.0e-324  1.56688e-315  5.0e-324  1.56687e-315  5.0e-324
 1.56687e-315  5.0e-324  1.56671e-315  5.0e-324  1.5667e-315   5.0e-324  1.5667e-315   5.0e-324  1.5667e-315   5.0e-324  1.5667e-315      1.56669e-315  5.0e-324  1.56668e-315  5.0e-324  1.56668e-315  5.0e-324  1.56668e-315  5.0e-324  1.56667e-315  5.0e-324

Each column is a state, and each row is a saved timestep. The values for every other column just repeat, within floating point error. Can poor integrator tolerances alone explain this kind of behaviour, or could something else be causing a problem?

Also I switched to using non-truncated distributions, as in #74, and can confirm that it seems to be working (at least I'm getting the same results that I was getting before anyway).

@agerlach
Copy link
Contributor

agerlach commented Jan 17, 2023 •
edited
Loading

Also I switched to using non-truncated distributions, as in #74, and can confirm that it seems to be working (at least I'm getting the same results that I was getting before anyway).

"Working" as in gives the correct answer? In the end I think you need a simpler MWE to determine if this is a bug and where it might reside.

@joegilkes
Copy link
Author

"Working" as in functioning without a crash or error and giving the same results as above, i.e. unchanged from what I was getting with the truncated distributions. Just thought I'd include it as a point of note on the back of the other issue.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees
No one assigned
Labels
None yet
Projects
None yet
Milestone
No milestone
Development

No branches or pull requests
2 participants
@agerlach @joegilkes