New subtype for state space set that contains normal Vectors

[Datseris]

For performance reasons, and because often one works in low dimensional cases, StateSpaceSet always has the points as SVector. This becomes limiting in high dimensional systems with 100s or more dimensions.

We should allow for another type that contains the points in normal Julia vectors, but is still type parameterized for its dimension.

See also JuliaDynamics/Attractors.jl#76 and we have the same problem in ComplexityMeasures.jl when trying to compute histograms in 100-dimensional spaces (cc @kahaaga ).

[Datseris]

It should be a direct duplicate of StateSpaceSet with a SizedArray.

[ikottlarz]

Hi! I'd like to contribute here, as this would also help me a great deal to implement the simulation of spatiotemporal systems. Before I do anything, I have a question though:

Currently, many methods for AbstractStateSpaceSet are hardcoded so that they return SVectors (such as minima, maxima, Base.hcat, ...). I was thinking I would make a new subtype HighDimStateSpaceSet (or some other name) of AbstractStateSpaceSet, but now this would be dispatched to the methods that return SVectors just like StateSpaceSet. So I see two possibilities to deal with this:

1. Change the expected type from AbstractStateSpaceSet to StateSpaceSet in the existing methods, make HighDimStateSpaceSet a subtype of AbstractStateSpaceSet and overload all methods with versions that take HighDimStateSpaceSet as an argument and return a SizedVector instead of an SVector. This would end up in a lot of duplicate code, since essentially both statespaceset.jl and statespaceset_concrete.jl would have to be duplicated with just some types changed
2. Leave the expected types of methods in statespaceset.jl as AbstractStateSpaceSet and add an if condition in each method that checks with which subtype we're dealing.
3. Something entirely else

I feel like 1. would be more conform with Julia as a multiple dispatch language, but as I said above, I see this leading to a large amount of essentially duplicated code.

Also: Currently trajectory automatically returns a StateSpaceSet. Do we want the user to be able to decide whether they want a StateSpaceSet or HighDimStateSpaceSet or should this be a purely automated decision based on the dimensionality of the dynamical system (in any case this would have to be implemented in a separate PR in DynamicalSystemsBase)?

[Datseris]

The very first thing is to define a new function containertype or similar name that returns the element type of the container. For StateSpaceSet this will be SVector{D,T} while for the new Set it will be Vector{T}. Then the abstract interface will be changed to utilize this type and return the correct appropriate type instead of hardcoding SVector everywhere.

Then, as you suggested, a new type HighDimStateSpaceSet or VectorBasedSSSet can be made that subtypes the abstract type.

if conditions are not necessary to distinguish types. You can utilize multiple dispatch and this new function to write code that applies to both cases, such as

function f(ssset)
V = containertype(ssset)
return V(x) # or whatever
end

Currently trajectory automatically returns a StateSpaceSet. Do we want the user to be able to decide whether they want a StateSpaceSet or HighDimStateSpaceSet or should this be a purely automated decision based on the dimensionality of the dynamical system (in any case this would have to be implemented in a separate PR in DynamicalSystemsBase)?

A decision based on the dimensionality would not be type stable. Instead, we can add one more keyword to trajectory that decides whether the output should be returned as StateSpaceSet (default) or HighDimStateSpaceSet. Again in the source code of trajectory the containertype function can be used to make code that works for both cases. And yes, this needs to be a separate PR.

[kahaaga]

This issue ties in with another issue I've had: it should be possible to use MVectors as the underlying data point type for StateSpaceSet too. I opened an issue on this way back, but can't locate it atm after the ecosystem restructuring.

One possibility to control the return type is to parameterize StateSpaceSet such that the type of the internal container is a type parameter. We can then use a (keyword?) argument in methods defined for AbstractStateSpaceSet to dictate the return type, i.e., hcat(x, y; return_type::T = SVector). Then return_type = SVector would reproduce the behaviour of the current StateSpaceSet, T = MVector would return a StateSpaceSet with mutable points, T = Vector would use regular vectors, T = Array for spatiotemporal systems, etc.

[Datseris]

Using keywords to change the container type in hcat is not something I would advise. We can satisfy both @kahaaga and @ikottlarz by altering the internals of StateSpaceSet to include a third type parameter V that is the type of the container. Actually, this is by far the simplest and cleanest solution. We would need to slightly expand some code, but all other functions such as hcat or whatever would "inherit" the third type parameter V and ensure they use vectors of type V instead of forcing SVectors.

[Datseris]

All constructors of StateSPaceSet can otain V as a first optional argument, which will default to SVector, breaking nothing. This is how Julia handles optional return types, like round(Int, x) versus round(x).

[kahaaga]

Haha, you beat me to the finishing line by seconds, @Datseris!

We can satisfy both @kahaaga and @ikottlarz by altering the internals of StateSpaceSet to include a third type parameter V that is the type of the container. Actually, this is by far the simplest and cleanest solution. We would need to slightly expand some code, but all other functions such as hcat or whatever would "inheritthe third type parameterVand ensure they use vectors of typeV` instead of forcing SVectors.

I also think this is the way to go, but:

Using keywords to change the container type in hcat is not something I would advise.

What would the return type be when doing e.g. hcat(x::StateSpaceSet{SVector}, y::StateSpaceSet{MVector}? I guess we'd need to use promotion rules, or allow some other mechanism to control it?

[kahaaga]

All constructors of StateSPaceSet can otain V as a first optional argument, which will default to SVector, breaking nothing. This is how Julia handles optional return types, like round(Int, x) versus round(x).

Agreed. Could we potentially also use the same optional argument to control the return type in the case of potential ambiguities, e.g. hcat(::T, ::StateSpaceSet{SVector}, ::StateSpaceSet{MVector})?, where T dictates the return type?

[Datseris]

What would the return type be when doing e.g. hcat(x::StateSpaceSet{SVector}, y::StateSpaceSet{MVector}? I guess we'd need to use promotion rules, or allow some other mechanism to control it?

We use the same promotion rules as in base Julia if they exist, other wise error. There is a promotion rule for when doing SVector + MVector (the function promote_something) and we should use exactly the same.

[Datseris]

V = promote_type(containertype(A), containertype(B))

is the function.

[kahaaga]

We use the same promotion rules as in base Julia if they exist, other wise error. There is a promotion rule for when doing SVector + MVector (the function promote_something) and we should use exactly the same.

Excellent. I guess StaticArrays.jl has also defined relevant promotion rules for the case of SArrays and regular AbstractArray too, so it should be possible to fall back on that.

[kahaaga]

There actually not that many instances of promote_rule defined in StaticArrays.jl, so we'll have to define them ourselves. We'll have to make a few choices when it comes to type promotion.

I'm working on a PR atm, so will present a proposal shortly.

[Datseris]

What do you mean? There should only be one instance of promote rule defined, maybe two, one for vector + any_vector_from_staticarrays and one for svector + mvector

[Datseris]

@kahaaga or @ikottlarz did either of you have had any progress here?

[kahaaga]

@kahaaga or @ikottlarz did either of you have had any progress here?

Yes, but it's not finished. I'm focusing on getting the ComplexityMeasures.jl stuff ready before I finish my work on this issue.

[Datseris]

if you have some code maybe open a PR and someone else can finish it? (likely, me)

[Datseris]

I am doing this now