Resizing ODE integrator with matrix or VectorOfArray initial condition

General Usage
,
1

Consider the following:

using OrdinaryDiffEq
using ElasticArrays
function sysde!(du, u, p, t)
    du[:, 1] .= rand(2)
    du[:, 2] .= rand(2)
    du[:, 3] .= rand(2)
    nothing 
end
u0 = ElasticArray{Float64}([0.0 0.0 0.0; 0.0 0.0 0.0])
tspan=(0.0,1.0)
prob=ODEProblem(sysde!, u0, tspan)
integrator = init(prob,Tsit5())
resize!(integrator, 4) # fails


ERROR: MethodError: no method matching resize!(::ElasticMatrix{Float64, Vector{Float64}}, ::Int64)
Closest candidates are:
  resize!(::OrdinaryDiffEq.ODEIntegrator, ::Int64) at C:\Users\licer\.julia\packages\OrdinaryDiffEq\c2Dwp\src\integrators\integrator_interface.jl:117
  resize!(::SciMLBase.DEIntegrator, ::Int64) at C:\Users\licer\.julia\packages\SciMLBase\cJ8FF\src\integrator_interface.jl:20
  resize!(::OrdinaryDiffEq.AbstractNLSolver, ::Any, ::Int64) at C:\Users\licer\.julia\packages\OrdinaryDiffEq\c2Dwp\src\nlsolve\utils.jl:407
  ...
Stacktrace:
 [1] resize!(integrator::OrdinaryDiffEq.ODEIntegrator{Tsit5, true, ElasticMatrix{Float64, Vector{Float64}}, Nothing, Float64, SciMLBase.NullParameters, Float64, Float64, Float64, Float64, Vector{ElasticMatrix{Float64, Vector{Float64}}}, ODESolution{Float64, 3, Vector{ElasticMatrix{Float64, Vector{Float64}}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{ElasticMatrix{Float64, Vector{Float64}}}}, ODEProblem{ElasticMatrix{Float64, Vector{Float64}}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, ODEFunction{true, SciMLBase.AutoSpecialize, typeof(sysde!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5, OrdinaryDiffEq.InterpolationData{ODEFunction{true, SciMLBase.AutoSpecialize, typeof(sysde!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Vector{ElasticMatrix{Float64, Vector{Float64}}}, Vector{Float64}, Vector{Vector{ElasticMatrix{Float64, Vector{Float64}}}}, OrdinaryDiffEq.Tsit5Cache{ElasticMatrix{Float64, Vector{Float64}}, ElasticMatrix{Float64, Vector{Float64}}, ElasticMatrix{Float64, Vector{Float64}}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}}}, DiffEqBase.DEStats}, ODEFunction{true, SciMLBase.AutoSpecialize, typeof(sysde!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, OrdinaryDiffEq.Tsit5Cache{ElasticMatrix{Float64, Vector{Float64}}, ElasticMatrix{Float64, Vector{Float64}}, ElasticMatrix{Float64, Vector{Float64}}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}}, OrdinaryDiffEq.DEOptions{Float64, Float64, Float64, Float64, PIController{Rational{Int64}}, typeof(DiffEqBase.ODE_DEFAULT_NORM), typeof(opnorm), Nothing, CallbackSet{Tuple{}, Tuple{}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, Nothing, Nothing, Int64, Tuple{}, Tuple{}, Tuple{}}, ElasticMatrix{Float64, Vector{Float64}}, Float64, Nothing, OrdinaryDiffEq.DefaultInit}, i::Int64)
   @ OrdinaryDiffEq C:\Users\licer\.julia\packages\OrdinaryDiffEq\c2Dwp\src\integrators\integrator_interface.jl:123
 [2] top-level scope
   @

I want to be able to resize the integrator so that u0 has an extra column. Is there any way to do that with resize!, or how can I define a new method to do it?

I initially wanted to use a vector of vectors, but that has the same issue:

using OrdinaryDiffEq
using RecursiveArrayTools
function sysde!(du, u, p, t)
    du[1] .= rand(2)
    du[2] .= rand(2)
    du[3] .= rand(2)
    nothing
end
u0 = VectorOfArray([[0.0,0.0],[0.0,0.0],[0.0,0.0]])
tspan = (0.0, 1.0)
prob = ODEProblem(sysde!, u0, tspan)
integrator = init(prob, Tsit5())
resize!(integrator, 4) # fails

ERROR: MethodError: no method matching resize!(::VectorOfArray{Float64, 2, Vector{Vector{Float64}}}, ::Int64)
Closest candidates are:
  resize!(::OrdinaryDiffEq.ODEIntegrator, ::Int64) at C:\Users\licer\.julia\packages\OrdinaryDiffEq\c2Dwp\src\integrators\integrator_interface.jl:117
  resize!(::SciMLBase.DEIntegrator, ::Int64) at C:\Users\licer\.julia\packages\SciMLBase\cJ8FF\src\integrator_interface.jl:20
  resize!(::OrdinaryDiffEq.AbstractNLSolver, ::Any, ::Int64) at C:\Users\licer\.julia\packages\OrdinaryDiffEq\c2Dwp\src\nlsolve\utils.jl:407
  ...
Stacktrace:
 [1] resize!(integrator::OrdinaryDiffEq.ODEIntegrator{Tsit5, true, VectorOfArray{Float64, 2, Vector{Vector{Float64}}}, Nothing, Float64, SciMLBase.NullParameters, Float64, Float64, Float64, Float64, Vector{VectorOfArray{Float64, 2, Vector{Vector{Float64}}}}, ODESolution{Float64, 3, Vector{VectorOfArray{Float64, 2, Vector{Vector{Float64}}}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{VectorOfArray{Float64, 2, Vector{Vector{Float64}}}}}, ODEProblem{VectorOfArray{Float64, 2, Vector{Vector{Float64}}}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, ODEFunction{true, SciMLBase.AutoSpecialize, typeof(sysde!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, 
Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5, OrdinaryDiffEq.InterpolationData{ODEFunction{true, SciMLBase.AutoSpecialize, typeof(sysde!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Vector{VectorOfArray{Float64, 2, Vector{Vector{Float64}}}}, Vector{Float64}, Vector{Vector{VectorOfArray{Float64, 2, Vector{Vector{Float64}}}}}, OrdinaryDiffEq.Tsit5Cache{VectorOfArray{Float64, 2, Vector{Vector{Float64}}}, VectorOfArray{Float64, 2, Vector{Vector{Float64}}}, VectorOfArray{Float64, 2, Vector{Vector{Float64}}}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}}}, DiffEqBase.DEStats}, ODEFunction{true, SciMLBase.AutoSpecialize, typeof(sysde!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, OrdinaryDiffEq.Tsit5Cache{VectorOfArray{Float64, 2, Vector{Vector{Float64}}}, VectorOfArray{Float64, 2, Vector{Vector{Float64}}}, VectorOfArray{Float64, 2, Vector{Vector{Float64}}}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}}, OrdinaryDiffEq.DEOptions{Float64, Float64, Float64, Float64, PIController{Rational{Int64}}, typeof(DiffEqBase.ODE_DEFAULT_NORM), typeof(opnorm), Nothing, CallbackSet{Tuple{}, Tuple{}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, Nothing, Nothing, 
Int64, Tuple{}, Tuple{}, Tuple{}}, VectorOfArray{Float64, 2, Vector{Vector{Float64}}}, Float64, Nothing, OrdinaryDiffEq.DefaultInit}, i::Int64)
   @ OrdinaryDiffEq C:\Users\licer\.julia\packages\OrdinaryDiffEq\c2Dwp\src\integrators\integrator_interface.jl:123
 [2] top-level scope
   @

(If a fix is possible, I’d much prefer the ElasticArray approach since that (for whatever reason) works 10x faster than the VectorOfArray approach for my application.)

2

I didn’t run any code, but could it be that you need to call resize! with two size arguments (e.g. for two-dimensional arrays instead of one-dimensional vectors)? For example resize!(integrator, 2, 4), just guessing.

1 Like
3

I considered that initially, e.g.

resize!(integrator, (4, 5))
resize!(integrator, 4, 5)

Doesn’t work, unfortunately (stacktrace is essentially the same).

1 Like
4

Actually, I think this is the solution - I should’ve been using

resize!(integrator, (2, i))

since ElasticArrays needs the number of rows in the first part of that tuple:

using OrdinaryDiffEq
using ElasticArrays
function sysde!(du, u, p, t)
    du[:, 1] .= rand(2)
    du[:, 2] .= rand(2)
    du[:, 3] .= rand(2)
    nothing 
end
u0 = ElasticArray{Float64}([0.0 0.0 0.0; 0.0 0.0 0.0])
tspan=(0.0,1.0)
prob=ODEProblem(sysde!, u0, tspan)
integrator = init(prob,Tsit5())
resize!(integrator, (2, 4)) # works

Not sure about the VectorOfArray case, but this works.

3 Likes