similar(v) broken for v::VectorOfArray

[Gregstrq]

similar(v) for v::VectorOfArray{T, N, A} gives you an instance of AbstractArray{T,N} instead of another VectorOfArray.
MWE:

a = VectorOfArray([[1,2,3,4,5],[6,7,8,9,10]])
similar(a)

gives you a 5×2 Array{Int64,2}.

I guess you can write something like

@inline Base.similar(VA::AbstractVectorOfArray) = VectorOfArray([similar(VA[i]) for i in eachindex(VA)])

to deal with this problem.

[ChrisRackauckas]

I actually thought we already had this implemented. Yes, that is how I'd do it (along with one for DiffEq array that just passes along the .t). Would you mind putting in the PR?

[Gregstrq]

No problem.
I am just a bit busy nowadays, so it will take several days.

[Gregstrq]

It seems that VectorOfArray breaks ODE_DEFAULT_NORM function from DiffEqBase.jl/src/common_defaults (it dispatches on AbstractArray and tries to apply UNITLESS_ABS2 to subarrays of VectorOfArray).
And we can also have ArrayPartition and VectorOfArray as partition of ArrayPartition.

Seems logical to treat this cases in the same manner as for AbstractArray{T,N} where {T<:AbstractArray,N}. Both ArrayPartition and VectorOfArray are composite arrays, essentially.

An easy workaround would be to write something like

@inline DiffEqBase.ODE_DEFAULT_NORM(u::RecursiveArrayTools.ArrayPartition) = sqrt(sum(DiffEqBase.ODE_DEFAULT_NORM, u.x)/length(u.x))
@inline DiffEqBase.ODE_DEFAULT_NORM(u::RecursiveArrayTools.AbstractVectorOfArray) = sqrt(sum(DiffEqBase.ODE_DEFAULT_NORM, u)/length(u))

But there is still a problem when u is, for example, ArrayPartition(Array, VectorOfArray).
In this case we have an array from the first partition and arrays from VectorOfArray. So if we use the piece of code which I wrote, we will first average over each of the arrays, then we will average over averages for arrays of VectorOfArray; finally, we average over partitions of ArrayPartition.
Long story short, we don't treat arrays coming VectorOfArray and arrays coming from other partitions of ArrayPartition equally.

Are the algorithms which use the norm are actually very sensitive to the choice of this norm?

[ChrisRackauckas]

Seems logical to treat this cases in the same manner as for AbstractArray{T,N} where {T<:AbstractArray,N}. Both ArrayPartition and VectorOfArray are composite arrays, essentially.

That can't be done because then eltype would have to be the internal array type, which would then make the indexing not work or throw a bunch of errors if it comes out as numbers.

An easy workaround would be to write something like

Yeah, those look good.

Are the algorithms which use the norm are actually very sensitive to the choice of this norm?

Not at all. This is choice the norm for the error in the adaptive timestepping. The Hairer 2-seminorm is just standard. But even then, Hairer's dopri5 and dop853 had an error in the norm calculation (in the initial dt part) for years and nobody noticed, so that shows it really doesn't matter all that much. And some packages choose other defaults, like a maximum norm.

But there is still a problem when u is, for example, ArrayPartition(Array, VectorOfArray).

I think it's a problem with the iterator interface.

A = (rand(5),rand(5))
p = ArrayPartition(A)
for i in p
  @show i
end

iterates over the values, while VectorOfArray iterates over the internal arrays. If ArrayPartition was overloaded to iterate over p.x, then it should recurse the iteration and work.

[Gregstrq]

I think it's a problem with the iterator interface.

Are you proposing to completely redo the iteration for ArrayPartition?
Won't it break something?

[ChrisRackauckas]

Won't it break something?

It might, but that's what the tests are for. This should be run against OrdinaryDiffEq.jl as well to make sure it would still work over there. But it looks like we can fix a lot of stuff just by overloading any, all and maybe reduce.

[ChrisRackauckas]

fixed