Version 1.2.6 breaks callbacks

[ufechner7]

The following code works with StaticArrays 1.2.5, but fails with 1.2.6 and newer:

using Sundials, StaticArrays
# Tutorial example showing how to use an implicit solver 
# It simulates a falling mass.

const G_EARTH  = [0.0, 0.0, -9.81]

# Type definitions
const SimFloat = Float64
const Vec3     = MVector{3, SimFloat}

# Falling mass.
# State vector y   = mass.pos, mass.vel
# Derivative   yd  = mass.vel, mass.acc
# Residual     res = (y.vel - yd.vel), (yd.acc - G_EARTH)     
function res1(res, yd, y::MVector{S, SimFloat}, p, t) where S
    res[1:3] .= y[4:6] - yd[1:3]
    res[4:6] .= yd[4:6] - G_EARTH
end

vel_0 = Vec3(0.0, 0.0, 50.0)    # Initial velocity
pos_0 = Vec3(0.0, 0.0,  0.0)    # Initial position
acc_0 = Vec3(0.0, 0.0, -9.81)   # Initial acceleration
y0  = vcat(pos_0, vel_0)    # Initial pos, vel
yd0 = vcat(vel_0, acc_0)    # Initial vel, acc

tspan = (0.0, 10.2)         # time span

differential_vars = ones(Bool, 6)
prob = DAEProblem(res1, yd0, y0, tspan, differential_vars=differential_vars)

sol = solve(prob, IDA(), saveat=0.1)

Error message:

julia> include("src/DAE_Example.jl")
ERROR: LoadError: MethodError: no method matching res1(::Vector{Float64}, ::SizedVector{6, Float64, Vector{Float64}}, ::SizedVector{6, Float64, Vector{Float64}}, ::SciMLBase.NullParameters, ::Float64)
Closest candidates are:
  res1(::Any, ::Any, ::MVector{S, Float64}, ::Any, ::Any) where S at /home/ufechner/repos/KiteViewer/src/DAE_Example.jl:15
Stacktrace:
  [1] (::DAEFunction{true, typeof(res1), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing})(::Vector{Float64}, ::Vararg{Any, N} where N)
    @ SciMLBase ~/.julia/packages/SciMLBase/mndcy/src/scimlfunctions.jl:369
  [2] (::Sundials.var"#104#108"{DAEProblem{MVector{6, Float64}, MVector{6, Float64}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, DAEFunction{true, typeof(res1), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, Vector{Bool}}, Tuple{Int64}, Tuple{Int64}})(out::Vector{Float64}, du::SizedVector{6, Float64, Vector{Float64}}, u::SizedVector{6, Float64, Vector{Float64}}, p::SciMLBase.NullParameters, t::Float64)
    @ Sundials ~/.julia/packages/Sundials/sg9S0/src/common_interface/solve.jl:1111
  [3] __init(prob::DAEProblem{MVector{6, Float64}, MVector{6, Float64}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, DAEFunction{true, typeof(res1), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, Vector{Bool}}, alg::IDA{:Dense, Nothing, Nothing}, timeseries::Vector{Any}, ts::Vector{Any}, ks::Vector{Any}; verbose::Bool, dt::Nothing, dtmax::Float64, save_on::Bool, save_start::Bool, callback::Nothing, abstol::Float64, reltol::Float64, saveat::Float64, tstops::Vector{Float64}, maxiters::Int64, timeseries_errors::Bool, dense_errors::Bool, save_everystep::Bool, dense::Bool, save_timeseries::Nothing, save_end::Bool, progress::Bool, progress_name::String, progress_message::typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), advance_to_tstop::Bool, stop_at_next_tstop::Bool, userdata::Nothing, kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ Sundials ~/.julia/packages/Sundials/sg9S0/src/common_interface/solve.jl:1312
  [4] __solve(prob::DAEProblem{MVector{6, Float64}, MVector{6, Float64}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, DAEFunction{true, typeof(res1), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, Vector{Bool}}, alg::IDA{:Dense, Nothing, Nothing}, timeseries::Vector{Any}, ts::Vector{Any}, ks::Vector{Any}, recompile::Type{Val{true}}; kwargs::Base.Iterators.Pairs{Symbol, Float64, Tuple{Symbol}, NamedTuple{(:saveat,), Tuple{Float64}}})
    @ Sundials ~/.julia/packages/Sundials/sg9S0/src/common_interface/solve.jl:13
  [5] solve_call(_prob::DAEProblem{MVector{6, Float64}, MVector{6, Float64}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, DAEFunction{true, typeof(res1), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, Vector{Bool}}, args::IDA{:Dense, Nothing, Nothing}; merge_callbacks::Bool, kwargs::Base.Iterators.Pairs{Symbol, Float64, Tuple{Symbol}, NamedTuple{(:saveat,), Tuple{Float64}}})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/8q10H/src/solve.jl:61
  [6] solve_up(prob::DAEProblem{MVector{6, Float64}, MVector{6, Float64}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, DAEFunction{true, typeof(res1), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, Vector{Bool}}, sensealg::Nothing, u0::MVector{6, Float64}, p::SciMLBase.NullParameters, args::IDA{:Dense, Nothing, Nothing}; kwargs::Base.Iterators.Pairs{Symbol, Float64, Tuple{Symbol}, NamedTuple{(:saveat,), Tuple{Float64}}})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/8q10H/src/solve.jl:87
  [7] solve(prob::DAEProblem{MVector{6, Float64}, MVector{6, Float64}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, DAEFunction{true, typeof(res1), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, Vector{Bool}}, args::IDA{:Dense, Nothing, Nothing}; sensealg::Nothing, u0::Nothing, p::Nothing, kwargs::Base.Iterators.Pairs{Symbol, Float64, Tuple{Symbol}, NamedTuple{(:saveat,), Tuple{Float64}}})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/8q10H/src/solve.jl:73
  [8] top-level scope
    @ ~/repos/KiteViewer/src/DAE_Example.jl:31
  [9] include(fname::String)
    @ Base.MainInclude ./client.jl:444
 [10] top-level scope
    @ REPL[1]:1
in expression starting at /home/ufechner/repos/KiteViewer/src/DAE_Example.jl:31

I tested with Julia 1.6.4 on Ubuntu Linux, but the same problem happens with Julia 1.7

[ufechner7]

I guess there is an automatic conversion from MVector to SizedVector happening, which should NOT happen.

[mateuszbaran]

I think this should be fine if you change res1(res, yd, y::MVector{S, SimFloat}, p, t) to res1(res, yd, y::StaticVector{S, SimFloat}, p, t)? There is no automatic conversion from MVector to SizedVector but there was a change in reshaping: #933 .

[mateuszbaran]

Anyway, there is indeed a bug with reshaping, vec shouldn't change MVector to SizedVector.
BTW, this: https://github.com/SciML/Sundials.jl/blob/595bf1d75e068c03099c60e2d91607c0693fc219/src/common_interface/solve.jl#L1091 , even with correct reshaping in StaticArrays, can turn MArray into SizedVector. Reshaping (including vec) should return a view and thus SizedArray sometimes has to appear.

[ufechner7]

Well, with a SizedArray my code is 100 times slower than with an MVector. I mean, my real code, not this toy example. So I think there is a bug that should be fixed. I need to be able to make callbacks using MVectors.

[ufechner7]

Just tested it. It works with the newest version if I change MVector to SizedVector, but 100 times slower ...

[mateuszbaran]

I don't recall a single operation that is so much slower for SizedVector than for MVector -- IIRC the performance difference is mostly about creation of new MVector being about 2-3 times faster than SizedVector.

Anyway, your performance problems are not caused by that StaticArrays.jl bug but rather Sundials.jl isn't doing the right thing and it was actually a bug in StaticArrays.jl that caused it to work before.

[ufechner7]

Well, with an MVector I have zero memory allocations, and a lot with SizedVector:

# With  KVec3     = SizedVector{3, SimFloat}
# Time  (mean ± σ):   52.215 μs ± 37.302 μs  Memory estimate: 16.92 KiB, allocs estimate: 614.
#
# KVec3     = MVector{3, SimFloat}
#  Time  (mean ± σ):   509.464 ns ±  52.723 ns  Memory estimate: 0 bytes, allocs estimate: 0.

[ufechner7]

Can you explain a little bit more what was the bug in StaticArrays and what should be fixed in Sundials?

I afraid if I tell the maintainer from Sundials.jl I upgraded StaticArrays, things stopped to work and therefore he should fix Sundials he might think that I am crazy...

[mateuszbaran]

Well, with an MVector I have zero memory allocations, and a lot with SizedVector:

# With  KVec3     = SizedVector{3, SimFloat}
# Time  (mean ± σ):   52.215 μs ± 37.302 μs  Memory estimate: 16.92 KiB, allocs estimate: 614.
#
# KVec3     = MVector{3, SimFloat}
#  Time  (mean ± σ):   509.464 ns ±  52.723 ns  Memory estimate: 0 bytes, allocs estimate: 0.

It's hard to tell what can cause that. Maybe there is a method with too narrow type bounds, or maybe with MVector Julia's escape analysis can elide some allocations. In any case it's perfectly fine to keep working with MVectors.

Can you explain a little bit more what was the bug in StaticArrays and what should be fixed in Sundials?

I afraid if I tell the maintainer from Sundials.jl I upgraded StaticArrays, things stopped to work and therefore he should fix Sundials he might think that I am crazy...

You can just report that things like vec(copy(prob.du0)) cannot be expected to return an MVector when prob.du0 is an MArray with more than one dimension, and reference this issue and #933 . Briefly, vec is supposed to return an array sharing the same memory as the original, which AFAIK can't be done on MArray other than MVector, so we need to make a view of the original MArray which can then be wrapped in SizedVector.

[mateuszbaran]

#982 might be enough to fix your example.

[ufechner7]

Thanks for the quick commit! Lets see if it helps...

[ufechner7]

In the latest master version the issue is fixed. Shall I leave this open until there is a new release?

[mateuszbaran]

Yes, this is fixed on the master branch so this issue can be closed.