i’m writing a program where i’m using my own array type, and i’d like to use DifferentialEquations to solve some ODEs that show up. i was under the impression that it should work with any `AbstractArray`, but this is either not the case or (most likely) i’m doing something wrong.

below i illustrate the issue with a minimal working example.

``````using DifferentialEquations

struct my_type{A,N} <: AbstractArray{A,N}
data :: A
name :: String
end
my_type(data::AbstractArray{T,N}, name::String) where {T,N} =
my_type{typeof(data),N}(data, name)

Base.size(var::my_type) = size(var.data)

Base.getindex(var::my_type, i::Int) = var.data[i]
Base.getindex(var::my_type, I::Vararg{Int,N}) where {N} = var.data[I...]
Base.getindex(var::my_type, ::Colon) = var.data[:]
Base.getindex(var::my_type, kr::AbstractRange) = var.data[kr]

Base.setindex!(var::my_type, v, i::Int) = (var.data[i] = v)
Base.setindex!(var::my_type, v, I::Vararg{Int,N}) where {N} = (var.data[I...] = v)
Base.setindex!(var::my_type, v, ::Colon) = (var.data[:] .= v)
Base.setindex!(var::my_type, v, kr::AbstractRange) = (var.data[kr] .= v)

function rhs_test(f, p, t)
f
end

xmin   = -2.0*pi
xmax   =  2.0*pi
xnodes =  600
hx     = (xmax - xmin) / xnodes

xx = range(xmin, stop=xmax, length=600)

x0   = 0
w    = 0.4
A    = 1

f0    = A * exp.( -((xx .- x0) ./ w).^2 )
foo   = my_type(f0, "foo")

tspan = (0.0, 1.0)

# this works fine
prob  = ODEProblem(rhs_test, f0, tspan)
sol   = solve(prob, RK4())

# this does not
prob2 = ODEProblem(rhs_test, foo, tspan)
sol2  = solve(prob2, RK4())
``````

basically i’m defining my own subtype of `AbstractArray` and define the `getindex` and `setindex!` methods as specified in the documentation, but this doesn’t seem to be enough as i get the error

``````ERROR: LoadError: MethodError: Cannot `convert` an object of type Float64 to an
object of type Array{Float64,1}
``````

am i missing some method, is this not possible at all, or is it the wrong way of approaching the problem?

Custom number types is super easy. You just need the standard arithmetic (`+`, `-`, `/`, `*`, and then `sqrt` if you don’t use a custom `internalnorm`). Custom array types take a little bit more. They need to use a compatible number type, be validly defined, define similar (see the AbstractArray interface, it mentions the `similar` methods you need), handle scalar multiples, and compatible additions. You need compatibility with your chosen linear solver or Jacobian if using an implicit method. You need a valid in-place broadcast for in-place methods (usually this can happen automatically if you define indexing).

``````struct my_type{A,N} <: AbstractArray{A,N}
data :: A
name :: String
end
``````

what you’re saying here is that the element type of your array is `A`, which in your example is `eltype(foo) == Vector{Float64}`. Thus when it’s doing its setup, it thinks that the “number type” you’re using is `Vector{Float64}` and errors. What you meant to say was that your element type was `Float64`, as in:

``````struct my_type{A,N} <: AbstractArray{A,N}
data :: Vector{A}
name :: String
end
my_type(data::AbstractArray{T,N}, name::String) where {T,N} =
my_type{eltype(data),N}(data, name)
``````

Then your array type now makes sense. However, it’s not closed under arithmetic. For example:

``````typeof(2foo) # Vector{Float64}
typeof(foo + foo) # Vector{Float64}
``````

How is the integrator supposed to know how to propagate your `my_type`? We need some kind of definition choice here. For fun, I’m going to say that we only take the left strings. So let’s define:

``````Base.:+(x::my_type,y::my_type) = my_type(x.data+y.data,x.name)
Base.:*(x::Number,y::my_type) = my_type(x*y.data,y.name)
Base.:/(x::my_type,y::Number) = my_type(x.data/y,x.name)
Base.similar(foo::my_type) = my_type(similar(foo.data),foo.name)
``````

Technically, according to the AbstractArray interface you should add a few more similar methods, but this works. Together, the code is:

``````using DifferentialEquations

struct my_type{A,N} <: AbstractArray{A,N}
data :: Vector{A}
name :: String
end
my_type(data::AbstractArray{T,N}, name::String) where {T,N} =
my_type{eltype(data),N}(data, name)

Base.size(var::my_type) = size(var.data)

Base.getindex(var::my_type, i::Int) = var.data[i]
Base.getindex(var::my_type, I::Vararg{Int,N}) where {N} = var.data[I...]
Base.getindex(var::my_type, ::Colon) = var.data[:]
Base.getindex(var::my_type, kr::AbstractRange) = var.data[kr]

Base.setindex!(var::my_type, v, i::Int) = (var.data[i] = v)
Base.setindex!(var::my_type, v, I::Vararg{Int,N}) where {N} = (var.data[I...] = v)
Base.setindex!(var::my_type, v, ::Colon) = (var.data[:] .= v)
Base.setindex!(var::my_type, v, kr::AbstractRange) = (var.data[kr] .= v)

function rhs_test(f, p, t)
f
end

xmin   = -2.0*pi
xmax   =  2.0*pi
xnodes =  600
hx     = (xmax - xmin) / xnodes

xx = range(xmin, stop=xmax, length=600)

x0   = 0
w    = 0.4
A    = 1

f0    = A * exp.( -((xx .- x0) ./ w).^2 )
foo   = my_type(f0, "foo")

tspan = (0.0, 1.0)

Base.:+(x::my_type,y::my_type) = my_type(x.data+y.data,x.name)
Base.:*(x::Number,y::my_type) = my_type(x*y.data,y.name)
Base.similar(foo::my_type) = my_type(similar(foo.data),foo.name)

prob  = ODEProblem(rhs_test, foo, tspan)
sol   = solve(prob, RK4())
``````

(Note that to work on master you also need to define `/`, but I just put in a quick PR to remove that requirement). Now if you want to make your type work with in-place, then you don’t need the arithmetic since that’s defined by broadcast. Here you just need one of your other missing `similar` methods:

``````Base.similar(foo::my_type,::Type{T}) where T = my_type(similar(foo.data,T),foo.name)
``````

So a full working code for that is:

``````using DifferentialEquations

struct my_type{A,N} <: AbstractArray{A,N}
data :: Vector{A}
name :: String
end
my_type(data::AbstractArray{T,N}, name::String) where {T,N} =
my_type{eltype(data),N}(data, name)

Base.size(var::my_type) = size(var.data)

Base.getindex(var::my_type, i::Int) = var.data[i]
Base.getindex(var::my_type, I::Vararg{Int,N}) where {N} = var.data[I...]
Base.getindex(var::my_type, ::Colon) = var.data[:]
Base.getindex(var::my_type, kr::AbstractRange) = var.data[kr]

Base.setindex!(var::my_type, v, i::Int) = (var.data[i] = v)
Base.setindex!(var::my_type, v, I::Vararg{Int,N}) where {N} = (var.data[I...] = v)
Base.setindex!(var::my_type, v, ::Colon) = (var.data[:] .= v)
Base.setindex!(var::my_type, v, kr::AbstractRange) = (var.data[kr] .= v)

function rhs_test(f, p, t)
f
end

xmin   = -2.0*pi
xmax   =  2.0*pi
xnodes =  600
hx     = (xmax - xmin) / xnodes

xx = range(xmin, stop=xmax, length=600)

x0   = 0
w    = 0.4
A    = 1

f0    = A * exp.( -((xx .- x0) ./ w).^2 )
foo   = my_type(f0, "foo")

tspan = (0.0, 1.0)

Base.similar(foo::my_type) = my_type(similar(foo.data),foo.name)
Base.similar(foo::my_type,::Type{T}) where T = my_type(similar(foo.data,T),foo.name)
function rhs_test2(df, f, p, t)
df.data .= f.data
end
prob  = ODEProblem(rhs_test2, foo, tspan)
sol   = solve(prob, RK4())
``````

So even for arrays it’s not bad: only two methods not on the `AbstractArray` interface were needed for out-of-place to define arithmetic, and in-place only needed the AbstractArray interface (but needed to make sure you define the `similar` part!).

Note that for these “extra” features we have a little bit more to do in order to get back to where we were in Julia v0.6. Since the definition of an AbstractArray and Broadcast has changed 3 times over the last 2 years, we haven’t been able to sit down and properly define our interface here. However, with the stability of v1.0 I hope we can finally do so soon .

5 Likes

When defining a wrapper type like your custom array type, I often find the macro `@forward` from Lazy.jl very convenient. It allows you to “forward” a set of function on your special type to the wrapped field inside, this way you do not have to manually define all those functions to simply operate on the field.

Example

``````julia> using Lazy

julia> struct MM
x
end

julia> @forward MM.x (Base.length, Base.iterate, Base.getindex)

julia> m = MM([1,2,3])
MM([1, 2, 3])

julia> length(m)
3

julia> m[3]
3
``````
6 Likes

this is it, thanks!

Is it possible to do this with an immutable type too? I.e. with a type:

``````using StaticArrays
struct my_type{N,A} <: AbstractVector{A}
data :: SVector{N,A}
name :: String
end
my_type(data::AbstractVector{T}, name::String) where {T} =
my_type{length(data), eltype(data)}(SVector{length(data)}(data), name)
``````

In this case `setindex!` cannot be defined and also `similar` will have to return a different type. The fact that StaticArrays can be used in DiffEq suggests that it is possible. But a naive modification of above (not defining `setindex!` and `similar`) does not work.

Then you also need to change to an out-of-place function, i.e. omit the 1st argument and return the value:

``````function rhs_test2_oop(f, p, t)
f # return value is the derivative
end
prob  = ODEProblem(rhs_test2_oop, foo, tspan)
sol   = solve(prob, RK4())
``````

Yes, I messed something up before, The example then runs but has `typeof(sol[1])==Vector` and not `Mtype`.

And thus when trying to carry a parameter around, it does not work:

``````using DifferentialEquations, StaticArrays

struct Mtype{N,A} <: AbstractVector{A}
data :: SVector{N,A}
name :: String
p::Float64
end
Mtype(data::AbstractVector{T}, name::String, p) where {T} =
Mtype{length(data), eltype(data)}(SVector{length(data)}(data), name, p)

Base.size(var::Mtype{N}) where N = (N,)

Base.getindex(var::Mtype, i::Int) = var.data[i]
Base.IndexStyle(::Mtype) = IndexLinear()

rhs_test(f, p, t) = f*f.p

xmin   = -2.0*pi
xmax   =  2.0*pi
xnodes =  10
hx     = (xmax - xmin) / xnodes

xx = range(xmin, stop=xmax, length=10)

x0   = 0
w    = 0.4
A    = 1

f0    = A * exp.( -((xx .- x0) ./ w).^2 )
foo   = Mtype(f0, "foo", 2.0)

tspan = (0.0, 1.0)

Base.:+(x::Mtype,y::Mtype) = Mtype(x.data+y.data,x.name)
Base.:*(x::Number,y::Mtype) = Mtype(x*y.data,y.name)

prob  = ODEProblem(rhs_test, foo, tspan)
sol   = solve(prob, RK4())
``````

errors with `ERROR: LoadError: type Array has no field p`. What seems to happen is that `similar(::Mtype)` is used, which creates an `Array` and thus has no field `foo.p`.

``````struct Mtype{N,A} <: AbstractVector{A}
data :: SVector{N,A}
name :: String
p::Float64
end
Mtype(data::AbstractVector{T}, name::String, p) where {T} =
Mtype{length(data), eltype(data)}(SVector{length(data)}(data), name, p)

Base.size(var::Mtype{N}) where N = (N,)

Base.getindex(var::Mtype, i::Int) = var.data[i]
Base.IndexStyle(::Mtype) = IndexLinear()

rhs_test(f, p, t) = f*f.p

xmin   = -2.0*pi
xmax   =  2.0*pi
xnodes =  10
hx     = (xmax - xmin) / xnodes

xx = range(xmin, stop=xmax, length=10)

x0   = 0
w    = 0.4
A    = 1

f0    = A * exp.( -((xx .- x0) ./ w).^2 )
foo   = Mtype(f0, "foo", 2.0)

foo .+ foo
``````

This works:

``````using DiffEqBase, OrdinaryDiffEq, StaticArrays
using Base: @propagate_inbounds

struct Mtype{N,T} <: AbstractVector{T}
data :: SVector{N,T}
name :: String
p::Float64
end
Mtype(data::AbstractVector{T}, name::String, p) where {T} =
Mtype{length(data), T}(SVector{length(data)}(data), name, p)
Base.copy(m::Mtype) = deepcopy(m)
Base.zero(m::Mtype{N,T}) where {N,T} = Mtype{N,T}(zero(SVector{N,T}), m.name, m.p)

# AbstractVector interface
Base.size(var::Mtype{N}) where N = (N,)
@propagate_inbounds Base.getindex(var::Mtype, i::Int) = var.data[i]
Base.IndexStyle(::Mtype) = IndexLinear()

# Whenever you subtype AbstractArrayStyle, you also need to define rules for combining
# dimensionalities, by creating a constructor for your style that takes a Val(N) argument.
MtypeStyle{N,T}(::Val{0}) where {N,T} = MtypeStyle{N,T}()
MtypeStyle{N,T}(::Val{1}) where {N,T} = MtypeStyle{N,T}()

# traverse Broadcasted tree to find Mtype instances
name, p = get_name_p(bc)

# TODO: this is what causes the slowness:
Mtype{N,T}(SVector(tuple((v for v in bc)...)), name, p)
# e.g compared with this:
# Mtype{N,T}(zeros(SVector{10,Float64}), name, p)
end

"Retrieve the fields name and p from the Broadcasted instance and make sure they all match"
name::Union{String, Nothing}=nothing
p::Union{Float64, Nothing}=nothing
for a in bc.args
if a isa Mtype
if name===nothing
name = a.name
p = a.p
else
name==a.name || error("Broadcasting of Mtypes with different name, p not supported.")
p==a.p || error("Broadcasting of Mtypes with different name, p not supported.")
end
nn, pp = get_name_p(a)
if name===nothing
name = nn
p = pp
elseif nn!==nothing
name==nn || error("Broadcasting of Mtypes with different name, p not supported.")
p==pp || error("Broadcasting of Mtypes with different name, p not supported.")
end
end
end
return name, p
end

# Test it
rhs_test(f::Mtype, p, t) = f*f.p

f0    = rand(10)
sf0   = SVector{length(f0)}(f0)
foo   = Mtype(sf0, "foo", 2.0)

tspan = (0.0, 100.0)

prob  = ODEProblem(rhs_test, foo, tspan)
sol   = solve(prob, RK4())
@time solve(prob, RK4()) # 0.109459 seconds (1.95 M allocations: 71.409 MiB, 13.93% gc time)

# test against straight SVector version
rhs_test2(f, p, t) = f*2
prob2  = ODEProblem(rhs_test2, sf0, tspan)
sol2   = solve(prob2, RK4())
@time solve(prob2, RK4()) # 0.000486 seconds (3.62 k allocations: 427.281 KiB)

# against Vector version
prob2  = ODEProblem(rhs_test2, f0, tspan)
sol2   = solve(prob2, RK4())
@time solve(prob2, RK4()) # 0.005367 seconds (52.43 k allocations: 7.910 MiB)

using Test
for i = 1:length(sol)
@test sol[i]==sol2[i]
end
``````

but is slow, about 200x slower than StaticArray (times are in the comments). The reason is the
line `Mtype{N,T}(SVector(tuple((v for v in bc)...)), name, p)` in the `copy` function.

You might want to make this `<: StaticArray`. That might make it have some better dispatches. There are likely some other ways to speed this up too.

The problem with subtyping from StaticVector is that then I need to define `similar_type` https://juliaarrays.github.io/StaticArrays.jl/stable/pages/api.html#Implementing-your-own-types-1. But that assumes, as far as I can tell, that all the fields are part of the array, i.e. there cannot be any metadata `name`, `p`.

Otherwise that would be perfect (or similarly working with FieldVector).

This recovers the fastness of Static arrays, using Setfield.jl

``````function Base.copy(bc::Broadcasted{<:MtypeStyle{N,T}}) where {N,T}
# traverse Broadcasted tree to find Mtype instances
name, p = get_name_p(bc)
# This trick works:
# https://discourse.julialang.org/t/constructing-svector-with-a-loop/15372/7
s = SVector{N,T}(1:N)
for (i,v) in enumerate(bc)
@set s[i] = v
end
Mtype{N,T}(s, name, p)
end
``````
2 Likes

Would something like this be a good addition to the docs? Where would it go?

I think you need to use `s = @set s[i] = v` or its shortcut `@set! s[i] = v`.

It would be a good addition to the docs, probably an FAQ on how define state types, or a discussion on the problem interface page

1 Like

Oops, you’re right. This makes timing worse, but only about a factor 20x slower instead of 200x (and still ~2x slower than vector).

Aright, improving… The `Mtype` was not `isbits`, which obviously hurt performance a lot. Now it’s only 1.5x slower than the SVector version (when skipping `get_name_p`). There is an ugly hack with the `_copy` methods. If anyone can tell me how to improve on that, that be awesome.

I’ll try to put a PR for the docs together.

``````using DiffEqBase, OrdinaryDiffEq, StaticArrays
using Base: @propagate_inbounds
using Setfield
using Test, BenchmarkTools

struct Mtype{N,T} <: AbstractVector{T}
data::SVector{N,T}
name::Int
p::Float64
end
Mtype(data::AbstractVector{T}, name::Int, p) where {T} =
Mtype{length(data), T}(SVector{length(data)}(data), name, p)
Base.copy(m::Mtype) = deepcopy(m)
Base.zero(m::Mtype{N,T}) where {N,T} = Mtype{N,T}(zero(SVector{N,T}), m.name, m.p)
# NOTE: that only the array-parts are hashed for <:AbstractVector
tmp1   = Mtype(1:2, -1, 2.0)
tmp2   = Mtype(1:2, -1, 3.0)
@test_broken hash(tmp1)!=hash(tmp2)
@test isbits(tmp1)

# AbstractVector interface
Base.size(var::Mtype{N}) where N = (N,)
@propagate_inbounds Base.getindex(var::Mtype, i::Int) = var.data[i]
Base.IndexStyle(::Mtype) = IndexLinear()

# Whenever you subtype AbstractArrayStyle, you also need to define rules for combining
# dimensionalities, by creating a constructor for your style that takes a Val(N) argument.
MtypeStyle{N,T}(::Val{0}) where {N,T} = MtypeStyle{N,T}()
MtypeStyle{N,T}(::Val{1}) where {N,T} = MtypeStyle{N,T}()

return _copy(bc, bc.args...)
end
# these separate _copy functions speed things up by about 2x.
# But they sure are ugly.
@inline function _copy(bc::Broadcasted{<:MtypeStyle{N,T}}, a1::Mtype{N,T}, a2::Mtype{N,T}) where {N,T}
f = bc.f
a1.name==a2.name || error("Broadcasting of Mtypes with different name, p not supported.")
a1.p==a2.p || error("Broadcasting of Mtypes with different name, p not supported.")
end
@inline _copy(bc::Broadcasted{<:MtypeStyle{N,T}}, a1::Mtype{N,T}, a2) where {N,T} =
@inline _copy(bc::Broadcasted{<:MtypeStyle{N,T}}, a1, a2::Mtype{N,T}) where {N,T} =
@inline function _copy(bc::Broadcasted{<:MtypeStyle{N,T}}, args...) where {N,T}
# traverse Broadcasted tree to find Mtype instances
name, p = -1, 2.0 # get_name_p(bc) # this uses about 120μs and 0.4MiB

# This is slow:
# Mtype{N,T}(SVector(tuple((v for v in bc)...)), name, p)
# thus use Setfield.jl trick:
# https://discourse.julialang.org/t/constructing-svector-with-a-loop/15372/7
s = SVector{N,T}(1:N)
for (i,v) in enumerate(bc)
@set! s[i] = v
end
Mtype{N,T}(s, name, p)
end

"""
Retrieve the fields name and p from the Broadcasted
instance and make sure they all match

TODO: a wee bit slow, this needs to be done better.
"""
name::Union{Int, Nothing}=nothing
p::Union{Float64, Nothing}=nothing
for a in bc.args
if a isa Mtype
if name===nothing
name = a.name
p = a.p
else
name==a.name || error("Broadcasting of Mtypes with different name, p not supported.")
p==a.p || error("Broadcasting of Mtypes with different name, p not supported.")
end
nn, pp = get_name_p(a)
if name===nothing
name = nn
p = pp
elseif nn!==nothing
name==nn || error("Broadcasting of Mtypes with different name, p not supported.")
p==pp || error("Broadcasting of Mtypes with different name, p not supported.")
end
end
end
return name::Int, p::Float64
end

# Test it
rhs_test(f::Mtype, p, t) = f*f.p

f0    = rand(10)
sf0   = SVector{length(f0)}(f0)
foo   = Mtype(sf0, -1, 2.0)

tspan = (0.0, 100)

prob  = ODEProblem(rhs_test, foo, tspan)
sol   = solve(prob, RK4())
@btime solve(prob, RK4()) #  568.960 μs (6553 allocations: 1.32 MiB)

# test against straight SVector version
rhs_test2(f, p, t) = f*2
prob2  = ODEProblem(rhs_test2, sf0, tspan)
sol2   = solve(prob2, RK4())
@btime solve(prob2, RK4()) # 366.731 μs (3615 allocations: 426.84 KiB)

# against Vector version
prob2  = ODEProblem(rhs_test2, f0, tspan)
sol2   = solve(prob2, RK4())
@btime solve(prob2, RK4()) # 2.246 ms (52423 allocations: 7.91 MiB)

@test length(sol)==length(sol2)
for i = 1:length(sol)
@test sol[i]==sol2[i]
end
``````

On this topic, there’s currently an issue here that `ArrayPartition` is not compatible with esp. stiff solvers since linear algebra is not implemented for it, (it’s mentioned that `PartitionedODEProblem` use them internally so maybe it’s also affected? btw, `PartitionedODEProblem not defined` appears to be undefined, looks like code is absent)

Is there a roadmap about what needs to be done if one wishes to help with it (overloading the functions with `ArrayPartition` or copy it in a linearized way to a cache)? It can be very useful for a set of coupled odes that included matrices, vectors and numbers, e.g.

v’ = fv(m_i,v,s)
m_i’ = fm(m_j,v,s)
s’ = fs(m_i,v,s)
where v is a vector, m_{1…n} are matrices and s is a number.

right now, I’m handling it through reshaping v,m,s to 1d and vcat them together, speed is okay for small systems but for bigger systems it may require preallocating and passing the preallocated arrays via the `p` parameter (may still generate some allocations).

By the way, I suppose `s = SVector{N,T}(ntuple(i -> (@inbounds bc[i]), Val(N)))` is equivalent to the `@set!` solution in terms of type stability?

1 Like

For this case we would really want to do blocked linear algebra, but I am not sure of a good solution without something like Swizzles. We may just have to densify until we have more linear algebra tools.

1 Like