Improve API of `DEDataArray` for fields that are `AbstractArray` but immutable

[ronisbr]

Hi!

I am trying to create a workspace using DEDataArray as follow:

mutable struct Workspace{T} <: DEDataVector{T}
    x::Vector{T} = zeros(12)
    q::Quaternion = Quaternion(1,0,0,0)

where Quaternion is defined in ReferenceFrameRotations.jl. Here, I get the following error:

ERROR: LoadError: setindex! not defined for Quaternion{Float64}
Stacktrace:
 [1] error(::String, ::Type{T} where T) at ./error.jl:42
 [2] error_if_canonical_setindex(::IndexLinear, ::Quaternion{Float64}, ::Int64) at ./abstractarray.jl:1081
 [3] setindex! at ./abstractarray.jl:1072 [inlined]
 [4] copyto!(::Quaternion{Float64}, ::Quaternion{Float64}) at ./multidimensional.jl:962
 [5] recursivecopy!(::Quaternion{Float64}, ::Quaternion{Float64}) at /Users/ronan.arraes/.julia/packages/RecursiveArrayTools/YIS5B/src/utils.jl:25
 [6] macro expansion at /Users/ronan.arraes/.julia/packages/DiffEqBase/8A5uo/src/data_array.jl:60 [inlined]
 [7] recursivecopy!(::Workspace{Float64}, ::Workspace{Float64}) at /Users/ronan.arraes/.julia/packages/DiffEqBase/8A5uo/src/data_array.jl:60
 [8] apply_step!(::OrdinaryDiffEq.ODEIntegrator{Tsit5,true,Workspace{Float64},Nothing,Float64,Tuple{Configuration{Float64},Auxiliary{OrbitPropagatorJ4{Float64},Float64},StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}},Float64,Float64,Float64,Array{Workspace{Float64},1},ODESolution{Float64,2,Array{Workspace{Float64},1},Nothing,Nothing,Array{Float64,1},Array{Array{Workspace{Float64},1},1},ODEProblem{Workspace{Float64},Tuple{Float64,Float64},true,Tuple{Configuration{Float64},Auxiliary{OrbitPropagatorJ4{Float64},Float64},StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}},ODEFunction{true,typeof(dyn),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem},Tsit5,OrdinaryDiffEq.InterpolationData{ODEFunction{true,typeof(dyn),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Array{Workspace{Float64},1},Array{Float64,1},Array{Array{Workspace{Float64},1},1},OrdinaryDiffEq.Tsit5Cache{Workspace{Float64},Workspace{Float64},Workspace{Float64},OrdinaryDiffEq.Tsit5ConstantCache{Float64,Float64}}},DiffEqBase.DEStats},ODEFunction{true,typeof(dyn),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},OrdinaryDiffEq.Tsit5Cache{Workspace{Float64},Workspace{Float64},Workspace{Float64},OrdinaryDiffEq.Tsit5ConstantCache{Float64,Float64}},OrdinaryDiffEq.DEOptions{Float64,Float64,Float64,Float64,typeof(DiffEqBase.ODE_DEFAULT_NORM),typeof(opnorm),CallbackSet{Tuple{},Tuple{DiscreteCallback{typeof(condition_control_loop),typeof(affect_control_loop!),typeof(DiffEqBase.INITIALIZE_DEFAULT)}}},typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN),typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE),typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK),DataStructures.BinaryHeap{Float64,DataStructures.LessThan},DataStructures.BinaryHeap{Float64,DataStructures.LessThan},Nothing,Nothing,Int64,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}},Array{Float64,1},Array{Float64,1}},Workspace{Float64},Float64,Nothing,OrdinaryDiffEq.DefaultInit}) at /Users/ronan.arraes/.julia/packages/OrdinaryDiffEq/V4HnF/src/integrators/integrator_utils.jl:288

The problem is that Quaternion (like SVector) is an immutable AbstractArray. In the current code, this is handle as follows:

        if Tf <: Union{Number,SArray}
            expressions[i] = :( dest.$f = getfield( src, $qf ) )
        elseif Tf <: AbstractArray
            expressions[i] = :( recursivecopy!(dest.$f, getfield( src, $qf ) ) )
        else
            expressions[i] = :( dest.$f = deepcopy( getfield( src, $qf ) ) )
        end

Thus, an element of type SVector is copied as b = a. However, Quaternion will be copied using recursivecopy!. This, in turn, will call copyto! that is not defined.

Thus, I have two suggestions:

1. Create a simple API like DiffEqBase.ismutablearray(::Some_type) = true. If this is true, then proceed with recursivecopy!. Otherwise, copy the element using b = a.
2. (I think this is better) Check if the method setindex! is implemented for the type. If it is, then proceed with recursivecopy!. Otherwise, copy the element using b = a.

P.S.: If this is accepted, I can submit a PR.

[devmotion]

IMO we should probably just use ArrayInterface.ismutable (see https://github.com/SciML/ArrayInterface.jl).

[ChrisRackauckas]

yup that's what it was made for. This code just originates from before ArrayInterface.jl

[ronisbr]

Something like this?

         if !ArrayInterface.ismutable(f)
            expressions[i] = :( dest.$f = getfield( src, $qf ) )
        elseif Tf <: AbstractArray
            expressions[i] = :( recursivecopy!(dest.$f, getfield( src, $qf ) ) )
        else
            expressions[i] = :( dest.$f = deepcopy( getfield( src, $qf ) ) )
        end

[devmotion]

It should be !ArrayInterface.ismutable(Tf) AFAICT but that's basically what I had in mind.

[ChrisRackauckas]

yup

[ronisbr]

I think I am missing something (I always have problem with @generated functions). With the modification, the following code is still failing:

struct Quaternion{T} <: AbstractVector{T}
    q0::T
    q1::T
    q2::T
    q3::T
end

Base.size(::Quaternion) = 4
ArrayInterface.ismutable(::Type{<:Quaternion}) = false

mutable struct SimWorkspace3{T} <: DEDataVector{T}
    x::Vector{T}
    q::Quaternion{T}
end

a = SimWorkspace3([1.0,2.0,3.0], Quaternion(cosd(15), 0.0, 0.0, sind(15)))
b = SimWorkspace3([0.0,0.0,0.0], Quaternion(1.0, 0.0, 0.0, 0.0))
recursivecopy!(b,a)

ArrayInterface.ismutable(Tf) is returning true inside the generated function but false in REPL. Any ideas what I am doing wrong?

EDIT: I added the following print just before the if:

        println("Type: $Tf")
        println("ismutable: $(ArrayInterface.ismutable(Tf))")

Then I copied and pasted the previous code in REPL. I am getting this:

julia> using DiffEqBase, RecursiveArrayTools, ArrayInterface

julia> struct Quaternion{T} <: AbstractVector{T}
           q0::T
           q1::T
           q2::T
           q3::T
       end

julia> Base.size(::Quaternion) = 4

julia> ArrayInterface.ismutable(::Type{<:Quaternion}) = false

julia> mutable struct SimWorkspace3{T} <: DEDataVector{T}
           x::Vector{T}
           q::Quaternion{T}
       end

julia> a = SimWorkspace3([1.0,2.0,3.0], Quaternion(cosd(15), 0.0, 0.0, sind(15)))
3-element SimWorkspace3{Float64}:
 1.0
 2.0
 3.0

julia> b = SimWorkspace3([0.0,0.0,0.0], Quaternion(1.0, 0.0, 0.0, 0.0))
3-element SimWorkspace3{Float64}:
 0.0
 0.0
 0.0

julia> recursivecopy!(b,a)
Type: Array{Float64,1}
ismutable: true
Type: Quaternion{Float64}
ismutable: true
ERROR: MethodError: no method matching checkbounds_indices(::Type{Bool}, ::Base.OneTo{Int64}, ::NTuple{4,Int64})
Closest candidates are:
  checkbounds_indices(::Type{Bool}, ::Tuple{}, ::Tuple) at abstractarray.jl:531
  checkbounds_indices(::Type{Bool}, ::Tuple, ::Tuple) at abstractarray.jl:527
  checkbounds_indices(::Type{Bool}, ::Tuple{}, ::Tuple{}) at abstractarray.jl:535
  ...
Stacktrace:
 [1] checkbounds(::Type{Bool}, ::Quaternion{Float64}, ::Int64, ::Int64, ::Vararg{Int64,N} where N) at ./abstractarray.jl:481
 [2] checkbounds(::Quaternion{Float64}, ::Int64, ::Int64, ::Int64, ::Vararg{Int64,N} where N) at ./abstractarray.jl:502
 [3] copyto!(::Quaternion{Float64}, ::Quaternion{Float64}) at ./multidimensional.jl:959
 [4] recursivecopy!(::Quaternion{Float64}, ::Quaternion{Float64}) at /Users/ronan.arraes/.julia/packages/RecursiveArrayTools/YIS5B/src/utils.jl:25
 [5] macro expansion at /Users/ronan.arraes/tmp/DiffEqBase.jl/src/data_array.jl:60 [inlined]
 [6] recursivecopy!(::SimWorkspace3{Float64}, ::SimWorkspace3{Float64}) at /Users/ronan.arraes/tmp/DiffEqBase.jl/src/data_array.jl:60
 [7] top-level scope at REPL[16]:1

julia> ArrayInterface.ismutable(Quaternion{Float64})
false

No idea why...

[ChrisRackauckas]

is it supposed to be ArrayInterface.ismutable(::Quaternion) = false?

[ronisbr]

Same problem. I just used the same structure in the declaration of StaticArrays in ArrayInterfaces.jl:

https://github.com/SciML/ArrayInterface.jl/blob/65fff23126664e7bad07b15fc7bf77667de7c393/src/ArrayInterface.jl#L491

However, if I do ArrayInterface.ismutable(::Quaternion) = false, then ArrayInterface.ismutable(Quaternion{Float64}) will be true but ArrayInterface.ismutable(q), where q is a Quaternion will be false.

[ChrisRackauckas]

try a new REPL session?

[ronisbr]

Yes, I tried in a new REPL session.

[ronisbr]

Aha! It only works if I define ArrayInterfaces(::Type{<:Quaternion}) before using DiffEqBase. Is it expected?

julia> using ArrayInterface

julia> struct Quaternion{T} <: AbstractVector{T}
           q0::T
           q1::T
           q2::T
           q3::T
       end

julia> Base.size(::Quaternion) = 4

julia> ArrayInterface.ismutable(::Type{<:Quaternion}) = false

julia> using DiffEqBase, RecursiveArrayTools


julia> mutable struct SimWorkspace3{T} <: DEDataVector{T}
           x::Vector{T}
           q::Quaternion{T}
       end

julia> a = SimWorkspace3([1.0,2.0,3.0], Quaternion(cosd(15), 0.0, 0.0, sind(15)))
3-element SimWorkspace3{Float64}:
 1.0
 2.0
 3.0

julia> b = SimWorkspace3([0.0,0.0,0.0], Quaternion(1.0, 0.0, 0.0, 0.0))
3-element SimWorkspace3{Float64}:
 0.0
 0.0
 0.0

julia> recursivecopy!(b,a)
Type: Array{Float64,1}
ismutable: true
Type: Quaternion{Float64}
ismutable: false
3-element SimWorkspace3{Float64}:
 1.0
 2.0
 3.0

julia> a
3-element SimWorkspace3{Float64}:
 1.0
 2.0
 3.0

julia> b
3-element SimWorkspace3{Float64}:
 1.0
 2.0
 3.0

[ChrisRackauckas]

Yes, it's a generated function that's made during using time.

[ronisbr]

hum, I see. I thought that even being a generated function, at compile time (first time it is called), then it would be able to call the correct method.

I will prepare a PR.