Allow StaticArray or FieldVector with extra fields

[mauro3]

It would be nice if custom defined StaticArrays or FieldVectors could have extra fields which would not form part of the array. Example:

struct Point3D <: FieldVector{3, Float64}
    x::Float64
    y::Float64
    z::Float64
    proj::String # Proj4 or similar string
end

This feature would allow to use arbitrary fields of a type to be exposed as a vector, and then be used in e.g. diff-eq solver.

I don't think this can work with the current similar_type interface. If similar_type would take an instance instead of a type, then that could be used to make an appropriate constructor.

[andyferris]

Hmm, yes we've definitely thought about this before.

My current thoughts are: you want similar_type to return a constructor. It doesn't necessarily need to be a type, or "data-less". For what you want you might want the signature to be more like similar_type(old_instance, new_eltype, new_size). Not sure how painful that would be in implementation. (Of course, we should also rename the function so it doesn't refer to "type").

There are some more, deeper considerations. For example, the AbstractArray interface defines precisely how isequal, hash, isless and so-on are defined - they depend on the elements but never the containers. So these "extra fields" couldn't participate in comparisons or for distinguishing dictionary keys and so-on. At some point this kind of caveat becomes painful, and you might actually be better off not being an AbstractArray (remember in Julia, a subtype should obey the entire interface, and if you want only a subset of it you should use a different part of the type tree. Hopefully interfaces/traits/multiple inheretence will land one day to make the different parts of the interface much more modular).

[andyferris]

By the way, for geographic coordinates I have been thinking a non-AbstractVector could be quite nice in the end. Really, you want to be able to do:

• point - point --> vector
• point + vector --> point
• point - vector --> point

But not

• point + point
• scalar * point
• etc

So these really aren't very AbstractVector-like...

Does that make sense?

[mauro3]

For example, the AbstractArray interface defines precisely how isequal, hash, isless and so-on are defined - they depend on the elements but never the containers.

That's a very good point. Thanks for the explanation. So, probably this is a no-starter.

My current motivation that I made this nested, isbits type encoding a machine. Now I want to simulate it as an ODE. For that a vector-like thing of the variables to be simulated (which are a bunch of the fields of the full type) needs to be passed off to DiffEq. So, as far as DiffEq is concerned that type should be a vector but otherwise, say for hashing/dict, it would have to be breaking the AbstractArray interface.

[andyferris]

Yes.

The only solution I've been satisfied with so far is: carefully transform all your inputs to the same, (approximately) Euclidean CRS first, do the DiffEq / optimization style work with (static) vectors, and carefully transform your answers to the correct CRS when saving.

Ideally we wouldn't need to do this but the autodiff tools seem to well with arrays rather than opaque structs.

[c42f]

For the original example, it might be ok to have the CRS in the type (as a symbol?) But generally I agree with Andy that this provides limited value in practice, and it might be best to carefully bring all data into a shared CRS once at load time and once again before saving the results. It's not perfect but the alternative of having a huge proliferation of incompatible CRS types isn't great.

It's a bit harder to comment on your example of the machine without knowing specifics of how you've modeled the problem. Generally I suppose you want to be able to split the full state into two parts, one vector to be optimized and another vector of fixed hyperparameters. Do you have the case where you want to flexibly select which parameters are fixed and which are to be optimized? Thinking out loud — maybe you could solve that in general by defining a tagged scalar type ("tagged for optimization" vs "tagged as fixed") with the tag in the type?

[mauro3]

What I'm after here is a immutable version of the DEDataArray of DiffEqBase (docs, code). Its docs describe the motivation for it quite well, in particular in this subsection. Note that the DEDataArray also suffers the equality/hashing problem:

julia> using DiffEqBase

julia> mutable struct SimType{T} <: DEDataVector{T}
           x::Array{T,1}
           f1::T
       end

julia> SimType([1,2], 3) == SimType([1,2], 4)
true

I think for the use-case with DiffEq, that odd-equality is ok (although it should probably be mentioned in the docs). Whereas for an example like the Point3D it would not be ok as two point with the same coordinates in different projections are definitely not the same point.

Anyway, the use-case with DiffEq is quite compelling, maybe @ChrisRackauckas has some insights?

Last, this came from a attempt of mine to implement this by hand: https://discourse.julialang.org/t/juliadiffeq-with-custom-types/16294/17

I'll need to ponder @c42f's suggestion. Edit: Yes, this would be neat to automatically select which fields to solve for. But to use with DiffEq, it would still need to be turned into a vector.

[mauro3]

Here an example, now I use a wrapper around the "machine" which is supposed to be the vector-like thing:

# The machine
struct Plane
    liftratio::Float64
    fuelweight::Float64
    dryweight::Float64
    area::Float64
    maxthrust::Float64
    thrust::Float64
    skintemp::Float64
    coords::SVector{3,Float64}
    vel::SVector{3,Float64}
end

# Wrapper which should act as a Vector to be used in DiffEq
struct PlaneDynamicVars <: StaticVector{8,Float64}
    p::plane
end

# Indexing
function Base.getindex(pd::PlaneDynamicVars, i)
    if 1<=i<=3
        return pd.p.coords[i]
    elseif 4<=i<=6
        return pd.p.vel[i-3]
    elseif i==7
        # keep track of the fuel
        pd.p.fuelweight
    elseif i==8
        # keep track of the skin temperature, only relevant for
        # very fast planes
        pd.p.skintemp
    else
        BoundsError(pd, i)
    end
end
Base.size(::PlaneDynamicVars) = (8,)
Base.IndexStyle(::PlaneDynamicVars) = IndexLinear()

# Broadcasting is hard thus it would be nice to just hook into StaticArray

Thus, ideally, it should be possible to map arbitrary components of the type to the indices.

[c42f]

I just looked at this again trying to figure out what is actionable. Possible ways forward:

• We could have some tools (maybe a macro) for making it easier to define your PlaneDynamicVars wrapper.
• We could test / document that similar_type can carry variables along via a constructor-like object, as @andyferris suggested.

@mauro3 has PlaneDynamicVars been a good workaround for your application? An alternative approach which might work just as well (and be more efficient in general) would be if you had helper functions to split Plane into a flat vector of dynamic vars (SVector{8}) and static vars (...) to be passed as u0 and p parameters to ODEProblem. This seems like a common enough problem that @ChrisRackauckas might already have a great solution to it somewhere in the DiffEq ecosystem.

[juliohm]

Bumping this issue as it is related to JuliaGeometry/CoordinateTransformations.jl#75

I am trying to define a new custom static vector called Cartesian to represent coordinates in a geographic system. What is the minimum amount of code I need to write to inherit all the niceties of linear algebra with SVector instances? Basically I am trying to define a new SVector with the exact same behavior to avoid recompilations downstream when users load 3rd-party modules with functions defined on SVector.

Is there a base abstract type that I can use to inherit behavior of a new "SVector" type?