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.111155  0.87656   0.0963009  0.257129   0.11311     0.781497
  2.0  0.386306  0.668226  0.278077   0.52735    0.775634    0.440884
  3.0  0.348826  0.691687  0.964227   0.76448    0.0852624   0.52036
  4.0  0.520184  0.484929  0.013144   0.269905   0.204084    0.117737
  5.0  0.451525  0.247586  0.602631   0.793532   0.842029    0.0683416
  6.0  0.669966  0.640974  0.0171909  0.749553   0.433769    0.452694
  7.0  0.969881  0.130252  0.369225   0.518483   0.0141696   0.714084
  8.0  0.517605  0.93358   0.515445   0.0758114  0.193109    0.150559
  9.0  0.072176  0.474535  0.778122   0.0768485  0.80023     0.91305
 10.0  0.854239  0.735973  0.0102776  0.0299112  0.537046    0.938732

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.0406084
  6.0  0.01576
  7.0  0.104261
  8.0  0.12714
  9.0  0.370194
 10.0  0.48022

Array-like indexing

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

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

(a = 0.2571287355813575, b = 0.925267535912165, c = 0.12713970284423626, d = 0.3258728635315493)

Reducing functions

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

using Statistics

julia> maximum(st)

(a = 0.9698812177371097, b = 0.9982064541308482, c = 0.4802202427553709, d = 0.9919365440188083)

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 stack is fast - indexing a single value return a NamedTuple from all layers is usually 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.037 ns (0 allocations: 0 bytes)
(a = 0.2571287355813575, b = 0.925267535912165, c = 0.12713970284423626, d = 0.3258728635315493)

julia> @btime $st[1, 4]

  4.027 ns (0 allocations: 0 bytes)
(a = 0.2571287355813575, b = 0.925267535912165, c = 0.12713970284423626, d = 0.3258728635315493)