Update StaticArray tests for eltype requirements

[ChrisRackauckas]

No description provided.

[ranocha]

Does this mean that plain Arrays of SVectors will not be supported anymore and we have to switch to VectorOfArray?

[ChrisRackauckas]

See SciML/SciMLBase.jl#561. Basically I am trying to get things to a well-defined interface. This little bit sticks out as sore area of DiffEq where we don't throw any errors about arrays of arrays since 10 out of like 300 methods supports arrays of arrays. But we should be able to support such types of targets more broadly without overhead by just using a small wrapper that makes everything conform to the interface. If this is allowed, then we can throw informative errors as to what's going on and we can vastly clean up a lot of the code to not require hacks for specific methods to avoid assumptions about the lack of Number eltype. Since it's an undocumented bit only used downstream by you and @Datseris, it's a small discussion to be had with many good consequences to most users, and would make things like transitioning a code to implicit methods more readily available since the Vector of Vector form will likely never be able to support linear algebra appropriately.

As you can see from the PR, there is a keyword argument being made to bypass the interface erroring, but I hope that is only used during development and is not a crutch that is actually relied on for these kinds of tasks given that there should be zero overhead to using a more structured and interface-compliant form.

[ranocha]

We use arrays of SVectors at least for some solvers in Trixi.jl with explicit Runge-Kutta methods. If OrdinaryDiffEq.jl transitions to another interface requirement, it would be great to get a breaking release describing this change.

We would also need more interfaces such as resize! and pointer for VectorOfArray.

[ranocha]

As you can see from the PR, there is a keyword argument being made to bypass the interface erroring, but I hope that is only used during development and is not a crutch that is actually relied on for these kinds of tasks given that there should be zero overhead to using a more structured and interface-compliant form.

Which PR do you refer to? I could only see changes to the tests here and new docs at SciML/SciMLBase.jl#561

[ChrisRackauckas]

I'm not sure undocumented side pieces really constitute breaking. We are trying to impose such interfaces so that there is a strong sense of what is actually breaking.

We would also need more interfaces such as resize! and pointer for VectorOfArray.

Add suggestions in the interface PR and I'll happily accept that.

[avik-pal]

In BVDiffEq we were using a AbstractArray{<:AbstractArray} for specifying the initial guess on a mesh. We can upperbound SciMLBase there and cut a breaking release.

[jlchan]

Is bounding SciMLBase compat enough to address the impact of this change?

[avik-pal]

Not SciMLBase my bad. DiffEqBase needs a upperbound. See SciML/BoundaryValueDiffEq.jl#152

[ChrisRackauckas]

In BVDiffEq we were using a AbstractArray{<:AbstractArray} for specifying the initial guess on a mess. We can upperbound SciMLBase there and cut a breaking release.

BVDiffEq we need to specialize there because that's a legitimate use case, that's quite different.

[ranocha]

How do you plan to do that? We may use it as inspiration for us in Trixi.jl.

[ChrisRackauckas]

@AayushSabharwal success!!!!!!!!!!!!!!!!!!! 🎉🎉🎉🎉🎉🎉🎉🎉🎉🎉🎉🎉🎉🎉🎉🎉🎉🎉🎉🎉🎉🎉🎉

As a clarification to everyone, there's a few related things going on here. The first is the SymbolicIndexingInterface v0.3 and RecursiveArrayTools v3 releases. While this nominally was done in order to make sol[x] types of indexing simpler in implementation, better documented, more robust, and more performant, there were a few other things that needed to be addressed along the way. One of the things that was addressed along the way was AbstractVectorOfArray is no longer an AbstractVector. This is beacuse sol[i] = sol.u[i] i.e. it does the vector indexing instead of the array indexing, so sol[i,j] is matrix-like but sol[i] is not and so there's some weird fallback behaviors that can do things that are not correct. We have been debating this for awhile, and thus with the RecursiveArrayTools v3 release we deprecated that indexing so now sol[i] throws a depwarn, and so that is a pretty massive downstream change. This also fixes up some ReverseDiff.jl behavior too though: too many things assume linear indexing of a higher order tensor is the linear index in column order, so we just needed to say it does not conform to the interface.

But with that, a casualty is precisely this test set right now for arrays of static arrays. Now in OrdinaryDiffEq.jl, I've wanted since 2016 to try to support arrays of arrays in multiple ways because it's a nice abstraction for many use cases. However, this change of not relying on the AbstractArray fallback from AbstractVectorOfArray has only exposed the nightmare to which I knew it had become. Basically, vectors of vectors never worked in OrdinaryDiffEq.jl. There are specific cases in SimpleDiffEq.jl where @Datseris added a special dispatch to get it to work for I think two methods for DynamicalSystems.jl, and @ranocha I think got some SSPRK methods working. But almost every algorithm would give up and spit an incomprehensible error message. And now with some of the AbstractArray fallbacks gone, it exposed that there's more work to do in order to actually keep it working. Some things rely on an assumption that "static arrays act like numbers" so some of these cases only worked for Vector of StaticArrays. Yes, it's that icky: a small set of solvers worked with arrays of arrays but only if it's an array of static arrays. That's a documentation nightmare and it's taking a ton of work to keep this side channel open. In practice, it's not really supported but only has just barely enough to keep it alive.

I think the straw that really broke the camel's back was that I was running a workshop and people got hung up on this due to a bad error message, and then I saw the same thing on Discourse (https://discourse.julialang.org/t/solving-a-system-of-matrix-differential-equations-using-differentialequations-jl/107663/4). So this piece that doesn't really work was holding up any possibility of allowing most people to have the good experience of a clear error message. And, I'm having trouble fixing this side channel for the RecursiveArrayTools v3 release because it's essentially relying on some assumptions that break what's common in the interfaces (for example, that the element type of the array is supposed to give a number).

So I was left with a few choices:

1. Fix the loophole and continue like nothing has happened. We've been doing this for years and it's a ton of work, so I don't think this is a great option.
2. Plug the whole with a nice error message (Throw a nicer error if eltype is not a number DiffEqBase.jl#976) and tell everyone to use VectorOfArrays if they want to do this. Make sure VectorOfArray is then the optimal way of handling it, battle tested with good performance.
3. Do (2) but add a trait system so that certain algorithms can opt into not throwing the error message, and then do (1) so that the status quo keeps alive.
4. Automatically wrap the user's u0 in a VectorOfArray

I don't think I could feasibly do (3) because that would take a lot more time than I have, and we as an organization have been seeking to plug up all of our interface holes in order to make things more robust and documented anyways. So those two factors together said, it's time to rip off the band-aid and do (2). And to be real, it's been a long time coming, and it's best for the vast majority of users. (4) has many potential issues since it wouldn't dispatch as the user thinks, so I didn't want to open that can of worms.

But yes, there are these two codes that do use this side channel, I don't think it would be anyone else because it would be fairly difficult to even find out that these algorithms were exceptions for these specific cases 😅. It's not documented but you could argue it's public API, but I'd rather do a much more targeted check here. Again, we could do a trait system to opt specific algorithms out here https://github.com/SciML/DiffEqBase.jl/blob/master/src/solve.jl#L591, but I'm not convinced that's even a good idea in the long term anyways. In the long term, now that VectorOfArray does everything that is needed here, I think the solution is just, if you want to do Vector of Vectors just slap a VectorOfArray on it, and if you don't we throw an informative error message telling you the easy solution. And with that, now everything is optimized, everyone gets a clear error message, and everything conforms to a vastly simpler interface.

I'm definitely open for opinions though

[ChrisRackauckas]

BVDiffEq we need to specialize there because that's a legitimate use case, that's quite different.

In https://github.com/SciML/DiffEqBase.jl/blob/master/src/solve.jl#L591 we should add a check that it's not a BoundaryValueProblem. @avik-pal can you follow up with that?

How do you plan to do that? We may use it as inspiration for us in Trixi.jl.

I'm not so sure how that's applicable @ranocha . The issue is that for a BVP u0 being a Vector of Vectors is a valid choice since that's used to give a time-dependent initial condition. It's just not a valid object in the same way for an ODE.

[ranocha]

Well, our ODE is a spatial semidiscretization of a system of PDEs. From the semidiscretization's point of view, an array of vectors (for each point in space) is a natural approach.
Anyway, if SciML drops support of this quality-of-life improvement, we have to find a way around it with some additional wrapping. It may not be trivial, though - we also have to support StructArrays.jl wrapping StaticArrays.jl with a different data layout than classics arrays of static arrays.

Minor clarification: It's not only a handful of SSPRK methods - all explicit RK methods not doing anything special (like ROCK methods) have been supporting arrays of static arrays for years. Everything required for explicit RK methods is a basic implementation of a vector space interface.

[ChrisRackauckas]

Minor clarification: It's not only a handful of SSPRK methods - all explicit RK methods not doing anything special (like ROCK methods) have been supporting arrays of static arrays for years. Everything required for explicit RK methods is a basic implementation of a vector space interface.

Are you sure about that? Because broadcast doesn't work well on arrays of arrays. zero is an operation that's not defined IIRC. similar. At least if most did, I was unaware of it and it wasn't tested.

Anyway, if SciML drops support of this quality-of-life improvement, we have to find a way around it with some additional wrapping. It may not be trivial, though - we also have to support StructArrays.jl wrapping StaticArrays.jl with a different data layout than classics arrays of static arrays.

Does VectorOfArray wrapping not work for you case? It's covered everything that we've tried so far at least.

[ranocha]

Minor clarification: It's not only a handful of SSPRK methods - all explicit RK methods not doing anything special (like ROCK methods) have been supporting arrays of static arrays for years. Everything required for explicit RK methods is a basic implementation of a vector space interface.

Are you sure about that? Because broadcast doesn't work well on arrays of arrays. zero is an operation that's not defined IIRC. similar. At least if most did, I was unaware of it and it wasn't tested.

It did indeed work and was tested with Euler(), Tsit5() and Vern9() - you have just removed the tests in this PR.

[ChrisRackauckas]

But were there AD tests on it? I am going back through the AD testing and don't see any of those cases checked, nor do I see how they would've worked in general with reverse mode. I get the very basic "it can solve" was checked, but I don't see any other interface compatibility checks on this case beyond that, and from what I could tell those mostly failed.

What exactly do you need? I still don't quite understand why VectorOfArray would add any complexity or performance issues.

[ranocha]

We just needed the basic "it can solve" capability, so we never tested AD etc.

We need for example the option to use StructArrays of StaticArrays. Is there an official and stable way to get the array wrapped by a VectorOfArrays? I would have expected something like

(jl_YFxZqo) pkg> add StaticArrays StructArrays RecursiveArrayTools

(jl_YFxZqo) pkg> add StaticArrays StructArrays RecursiveArrayTools
...
  [731186ca] + RecursiveArrayTools v3.4.2
  [90137ffa] + StaticArrays v1.9.0
  [09ab397b] + StructArrays v0.6.16
...


julia> using StaticArrays, StructArrays, RecursiveArrayTools

julia> u = VectorOfArray(StructArray(randn(SVector{2, Float64}, 10)))
VectorOfArray{Float64,2}:
10-element StructArray(::Vector{Float64}, ::Vector{Float64}) with eltype SVector{2, Float64}:
 [-0.7640292999038704, 0.6553189479381939]
 [0.8378337172716237, 0.07143533596999121]
 [-0.5413059411379175, 1.545169798144512]
 [-0.1140800570149823, 2.0271290114354232]
 [-0.09736760972644991, 0.9762634582911948]
 [-0.9241989885704037, -0.7643537491391698]
 [-0.3234905544876145, 0.035062794764888564]
 [0.6533985389766392, 1.0330990755657259]
 [-1.7927415052914886, -0.8398035991158214]
 [-0.5444953510216021, 0.28585092256556105]

julia> parent(u)
ERROR: MethodError: no method matching parent(::VectorOfArray{Float64, 2, StructVector{SVector{2, Float64}, Tuple{Vector{Float64}, Vector{Float64}}, Int64}})

Closest candidates are:
  parent(::SizedArray)
   @ StaticArrays ~/.julia/packages/StaticArrays/PLKkM/src/SizedArray.jl:95
  parent(::StructArrays.GroupPerm)
   @ StructArrays ~/.julia/packages/StructArrays/0h2SD/src/sort.jl:14
  parent(::PermutedDimsArray)
   @ Base permuteddimsarray.jl:48
  ...

Stacktrace:
 [1] top-level scope
   @ REPL[5]:1

It looks like we could use

julia> u.u
10-element StructArray(::Vector{Float64}, ::Vector{Float64}) with eltype SVector{2, Float64}:
 [-0.7640292999038704, 0.6553189479381939]
 [0.8378337172716237, 0.07143533596999121]
 [-0.5413059411379175, 1.545169798144512]
 [-0.1140800570149823, 2.0271290114354232]
 [-0.09736760972644991, 0.9762634582911948]
 [-0.9241989885704037, -0.7643537491391698]
 [-0.3234905544876145, 0.035062794764888564]
 [0.6533985389766392, 1.0330990755657259]
 [-1.7927415052914886, -0.8398035991158214]
 [-0.5444953510216021, 0.28585092256556105]

but I couldn't find it in the docs.

[ranocha]

What exactly do you need? I still don't quite understand why VectorOfArray would add any complexity or performance issues.

I just took a look at implementing parent as discussed above. While doing this, I noticed that we really need/use multi-dimensional arrays of SVectors. This is not supported by VectorOfArray. Which alternative would you recommend?

[jlchan]

I'll second this; my unknowns are Matrix{SVector}. VectorOfArray doesn't handle this.

[ChrisRackauckas]

Does https://github.com/JuliaArrays/ArraysOfArrays.jl work for this?

[jlchan]

Some broadcast operations (like adding constants to an array) work with VectorOfArray but don't work with ArraysOfArrays, so maybe? I also remember that it allocated unnecessarily in some operations, but this was several years ago now.

I tried extending VectorOfArray to use a Matrix parent array. The broadcasting machinery seems to work since VectorOfArray just relies on linearly indexing of the parent. Do you see any downsides to this approach? If not, I can make a PR.

[ChrisRackauckas]

I tried extending VectorOfArray to use a Matrix parent array. The broadcasting machinery seems to work since VectorOfArray just relies on linearly indexing of the parent. Do you see any downsides to this approach? If not, I can make a PR.

The RecursiveArrayTools.jl one? Yeah it sounds fine to generalize that. I'd take a PR if you have code for it.

[jlchan]

Great, will add when I'm back