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?
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.
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.
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.