DimStacks

An AbstractDimStack represents a collection of AbstractDimArray layers that share some or all dimensions. For any two layers, a dimension of the same name must have the identical lookup - in fact only one is stored for all layers to enforce this consistency.

julia> using DimensionalData

julia> x, y = X(1.0:10.0), Y(5.0:10.0)

↓ X 1.0:1.0:10.0,
→ Y 5.0:1.0:10.0

julia> st = DimStack((a=rand(x, y), b=rand(x, y), c=rand(y), d=rand(x)))

╭───────────────╮
│ 10×6 DimStack │
├───────────────┴──────────────────────────────────────────── dims ┐
  ↓ X Sampled{Float64} 1.0:1.0:10.0 ForwardOrdered Regular Points,
  → Y Sampled{Float64} 5.0:1.0:10.0 ForwardOrdered Regular Points
├────────────────────────────────────────────────────────── layers ┤
  :a eltype: Float64 dims: X, Y size: 10×6
  :b eltype: Float64 dims: X, Y size: 10×6
  :c eltype: Float64 dims: Y size: 6
  :d eltype: Float64 dims: X size: 10
└──────────────────────────────────────────────────────────────────┘

The behaviour of a DimStack is at times like a NamedTuple of DimArray and, others an AbstractArray of NamedTuple.

NamedTuple-like indexing

Layers can be accessed with .name or [:name]

julia> st.a

╭────────────────────────────╮
│ 10×6 DimArray{Float64,2} a │
├────────────────────────────┴─────────────────────────────── dims ┐
  ↓ X Sampled{Float64} 1.0:1.0:10.0 ForwardOrdered Regular Points,
  → Y Sampled{Float64} 5.0:1.0:10.0 ForwardOrdered Regular Points
└──────────────────────────────────────────────────────────────────┘
  ↓ →  5.0        6.0        7.0        8.0       9.0        10.0
  1.0  0.42647    0.109647   0.352971   0.413182  0.727907    0.428786
  2.0  0.932088   0.638534   0.871102   0.844364  0.978476    0.799307
  3.0  0.168171   0.0218458  0.556944   0.716964  0.940252    0.419048
  4.0  0.274651   0.311835   0.803492   0.689708  0.0278592   0.446461
  5.0  0.607209   0.0118438  0.800428   0.4462    0.816343    0.923608
  6.0  0.0495992  0.607637   0.442958   0.466858  0.721096    0.172443
  7.0  0.481653   0.991939   0.848509   0.918917  0.873178    0.23117
  8.0  0.508449   0.845379   0.0408718  0.683317  0.0191855   0.881938
  9.0  0.557033   0.213854   0.32464    0.244024  0.162406    0.00802652
 10.0  0.259676   0.401827   0.468547   0.980044  0.822858    0.359875

julia> st[:c]

╭─────────────────────────────────╮
│ 6-element DimArray{Float64,1} c │
├─────────────────────────────────┴───────────────────────── dims ┐
  ↓ Y Sampled{Float64} 5.0:1.0:10.0 ForwardOrdered Regular Points
└─────────────────────────────────────────────────────────────────┘
  5.0  0.830848
  6.0  0.386346
  7.0  0.180918
  8.0  0.762112
  9.0  0.617365
 10.0  0.396534

Array-like indexing

Indexing with a scalar returns a NamedTuple of values, one for each layer:

julia> st[X=1, Y=4]

(a = 0.413182007748902, b = 0.6500838767926648, c = 0.7621122643241707, d = 0.7762823999009395)

Reducing functions

Base functions like mean, maximum, reverse are applied to all layers of the stack.

julia> maximum(st)

(a = 0.9919389861816956, b = 0.9703359276120026, c = 0.8308481472259899, d = 0.893325868995048)

julia> maximum(st; dims=Y)

╭───────────────╮
│ 10×1 DimStack │
├───────────────┴──────────────────────────────────────────── dims ┐
  ↓ X Sampled{Float64} 1.0:1.0:10.0 ForwardOrdered Regular Points,
  → Y Sampled{Float64} 7.5:6.0:7.5 ForwardOrdered Regular Points
├────────────────────────────────────────────────────────── layers ┤
  :a eltype: Float64 dims: X, Y size: 10×1
  :b eltype: Float64 dims: X, Y size: 10×1
  :c eltype: Float64 dims: Y size: 1
  :d eltype: Float64 dims: X size: 10
└──────────────────────────────────────────────────────────────────┘

broadcast_dims broadcasts functions over any mix of AbstractDimStack and AbstractDimArray returning a new AbstractDimStack with layers the size of the largest layer in the broadcast. This will work even if dimension permutation does not match in the objects.

Only atrix layers can be rotaed

julia> rotl90(st[(:a, :b)])

╭───────────────╮
│ 6×10 DimStack │
├───────────────┴───────────────────────────────────────────── dims ┐
  ↓ Y Sampled{Float64} 10.0:-1.0:5.0 ReverseOrdered Regular Points,
  → X Sampled{Float64} 1.0:1.0:10.0 ForwardOrdered Regular Points
├─────────────────────────────────────────────────────────── layers ┤
  :a eltype: Float64 dims: Y, X size: 6×10
  :b eltype: Float64 dims: Y, X size: 6×10
└───────────────────────────────────────────────────────────────────┘

julia> rotl90(st[(:a, :b)], 2)

╭───────────────╮
│ 10×6 DimStack │
├───────────────┴───────────────────────────────────────────── dims ┐
  ↓ X Sampled{Float64} 10.0:-1.0:1.0 ReverseOrdered Regular Points,
  → Y Sampled{Float64} 10.0:-1.0:5.0 ReverseOrdered Regular Points
├─────────────────────────────────────────────────────────── layers ┤
  :a eltype: Float64 dims: X, Y size: 10×6
  :b eltype: Float64 dims: X, Y size: 10×6
└───────────────────────────────────────────────────────────────────┘

Performance

Indexing stack is fast - indexing a single value return a NamedTuple from all layers is usally measures in nanoseconds, and no slower than manually indexing into each parent array directly.

There are some compilation overheads to this though, and stacks with very many layers can take a long time to compile.

julia> using BenchmarkTools

julia> @btime $st[X=1, Y=4]

  4.027 ns (0 allocations: 0 bytes)
(a = 0.413182007748902, b = 0.6500838767926648, c = 0.7621122643241707, d = 0.7762823999009395)

julia> @btime $st[1, 4]

  4.028 ns (0 allocations: 0 bytes)
(a = 0.413182007748902, b = 0.6500838767926648, c = 0.7621122643241707, d = 0.7762823999009395)