Changing size of ODE system: resize! works but deleteat! doesn't

I am playing around with the capacity of changing the size of an ODE system by tweaking the example presented here. For some reason, shrinking the model by resize!()ing the integrator works as expected, but the (equivalent?) deleteat!() fails with a “DimensionMismatch” or a “BoundsError”.

using DifferentialEquations
using Plots

const α = -0.3
function f(du,u,p,t)
  for i in 1:length(u)
    du[i] = α*u[i]
  end
end

function condition(u,t,integrator) # Event when event_f(u,t) == 0
  minimum(u)-0.01
end

function affect!(integrator)
  u = integrator.u
  resize!(integrator,length(u)-1) # This works
  # deleteat!(integrator, length(u)) # This errors
  nothing
end

callback = ContinuousCallback(condition,affect!)
u0 = [5, 3, 1]
tspan = (0.0,20.0)
prob = ODEProblem(f,u0,tspan)
sol = solve(prob,callback=callback)

plot(sol.t,map((x)->length(x),sol[:]),lw=3,
     ylabel="Number of Cells",xlabel="Time")

Can anybody help me understand what’s going on here?

I think what is going on is the following:

integrator is a specific structure, which has scalar and a vector. You want to resize the vector (u). Apparently, what means to resize the integrator structure correctly is implemented for that type of structure somewhere. This is not evident. For example:

julia> struct A
          a
          x
       end

julia> a = A(1,[1,2])
A(1, [1, 2])

julia> resize!(a,1)
ERROR: MethodError: no method matching resize!(::A, ::Int64)
Closest candidates are:
  resize!(::Array{T,1} where T, ::Integer) at array.jl:1082
  resize!(::BitArray{1}, ::Integer) at bitarray.jl:799
Stacktrace:
 [1] top-level scope at REPL[3]:1

There is no “resize!” function defined for that structure. Very likely, somewhere within the DifferentialEquations package there is a proper definition of resize! to deal with the the type of integrator, which handles things correctly for what it will be used, and there is no similar definition of deleteat!.

Simply resizing the integrator.u vector does not work either, it seems, therefore changing the size of that vector is not enough to keep things working under the hood.

Well, from the docs, it seems that both should work and handle correctly what is to be handled. Maybe that is a bug

From: Integrator Interface · DifferentialEquations.jl

  • resize!(integrator,k) : Resizes the DE to a size k . This chops off the end of the array, or adds blank values at the end, depending on whether k>length(integrator.u) .
  • resize_non_user_cache!(integrator,k) : Resizes the non-user facing caches to be compatible with a DE of size k . This includes resizing Jacobian caches. Note that in many cases, resize! simple resizes user_cache variables and then calls this function. This finer control is required for some AbstractArray operations.
  • deleteat_non_user_cache!(integrator,idxs) : deleteat! s the non-user facing caches at indices idxs . This includes resizing Jacobian caches. Note that in many cases, deleteat! simple deleteat! s user_cache variables and then calls this function. This finer control is required for some AbstractArray operations.
  • addat_non_user_cache!(integrator,idxs) : addat! s the non-user facing caches at indices idxs . This includes resizing Jacobian caches. Note that in many cases, addat! simple addat! s user_cache variables and then calls this function. This finer control is required for some AbstractArray operations.
  • deleteat!(integrator,idxs) : Shrinks the ODE by deleting the idxs components.
  • addat!(integrator,idxs) : Grows the ODE by adding the idxs components. Must be contiguous indices.

Thank you for the reply, I will look into it. And yes from reading the docs I also had the impression it should work, maybe a bug indeed.

It works when you specify the algorithm:

using DifferentialEquations
using Plots

const α = -0.3
function f(du,u,p,t)
  for i in 1:length(u)
    du[i] = α*u[i]
  end
end

function condition(u,t,integrator) # Event when event_f(u,t) == 0
  minimum(u)-0.01
end

function affect!(integrator)
  u = integrator.u
  #resize!(integrator,length(u)-1) # This works
  deleteat!(integrator, length(u)) # This errors
  nothing
end

callback = ContinuousCallback(condition,affect!)
u0 = [5, 3, 1]
tspan = (0.0,20.0)
prob = ODEProblem(f,u0,tspan)
sol = solve(prob,Tsit5(),callback=callback)

plot(sol.t,map((x)->length(x),sol[:]),lw=3,
     ylabel="Number of Cells",xlabel="Time")

The issue seems to be that if you use an automatic stiffness detection algorithm then, since they have some shared cache variables, it’s deleting from those shared cache variables twice. Could you file a bug report so we can track and get this fixed? It shouldn’t effect your usage if you select an algorithm, but it would be nice to fix it sooner or later.

1 Like

Excellent thank you, specifying an algorithm indeed does fix the problem. I opened an issue on github.