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 behavior of a DimStack is at times like a NamedTuple of DimArray and, at other times, 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.20259   0.0328886  0.93358    0.515445  0.00387281   0.947461
  2.0  0.717768  0.072176   0.474535   0.778122  0.603823     0.80023
  3.0  0.651877  0.854239   0.735973   0.662632  0.203451     0.537046
  4.0  0.111155  0.87656    0.0963009  0.97164   0.11311      0.781497
  5.0  0.386306  0.668226   0.817575   0.506181  0.775634     0.440884
  6.0  0.348826  0.691687   0.214157   0.76448   0.0852624    0.52036
  7.0  0.520184  0.0170682  0.389215   0.269905  0.204084     0.705701
  8.0  0.451525  0.345667   0.602631   0.793532  0.842029     0.351735
  9.0  0.103636  0.780424   0.0171909  0.749553  0.187655     0.669192
 10.0  0.38291   0.130252   0.369225   0.518483  0.488136     0.714084

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.191224
  6.0  0.943694
  7.0  0.177155
  8.0  0.0406084
  9.0  0.01576
 10.0  0.104261

Array-like indexing

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

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

(a = 0.5154448003579335, b = 0.7410056618197303, c = 0.04060837507113502, d = 0.43235163202114113)

Reducing functions

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

using Statistics

julia> maximum(st)

(a = 0.97164033242648, b = 0.9982064541308482, c = 0.9436941638102098, d = 0.8632634208267141)

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 matrix layers can be rotated

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 a stack is fast - indexing a single value and returning a NamedTuple from all layers is usually measured 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.067 ns (0 allocations: 0 bytes)
(a = 0.5154448003579335, b = 0.7410056618197303, c = 0.04060837507113502, d = 0.43235163202114113)

julia> @btime $st[1, 4]

  4.047 ns (0 allocations: 0 bytes)
(a = 0.5154448003579335, b = 0.7410056618197303, c = 0.04060837507113502, d = 0.43235163202114113)