product/mapgrid/outer

[jariji]

One feature of SplitApplyCombine that's missing from here is product/mapgrid/outer. Like this, but without the intermediate allocation.

julia> foo(f, a, b) = map(splat(f), grid(a,b));

julia> foo(*, 1:5, 10:10:100)
2-dimensional KeyedArray(...) with keys:
↓   5-element UnitRange{Int64}
→   10-element StepRange{Int64,...}
And data, 5×10 Matrix{Int64}:
      (10)  (20)  (30)  (40)  (50)  (60)  (70)  (80)  (90)  (100)
 (1)    10    20    30    40    50    60    70    80    90    100
 (2)    20    40    60    80   100   120   140   160   180    200
 (3)    30    60    90   120   150   180   210   240   270    300
 (4)    40    80   120   160   200   240   280   320   360    400
 (5)    50   100   150   200   250   300   350   400   450    500

[aplavin]

What exists now:

• grid(T, xs...) returns a grid of T constructed from (x, x, x, ...)

julia> grid(SVector, 1:2, 10:10:30)
2-dimensional KeyedArray(...) with keys:
↓   2-element UnitRange{Int64}
→   3-element StepRange{Int64,...}
And data, 2×3 RectiGrids.RectiGridArr{Base.OneTo(2), SVector{2, Int64}, 2, Tuple{Nothing, Nothing}, Tuple{UnitRange{Int64}, StepRange{Int64, Int64}}}:
      (10)         (20)         (30)
 (1)      [1, 10]      [1, 20]      [1, 30]
 (2)      [2, 10]      [2, 20]      [2, 30]

• mapview(T, grid(xs...)) in FlexiMaps returns a (lazy) view, similar to what you need :

julia> mapview(splat(*), grid(1:2, 10:10:30))
2×3 FlexiMaps.MappedArray{Int64, 2, Base.Splat{typeof(*)}, KeyedArray{Tuple{Int64, Int64}, 2, RectiGrids.RectiGridArr{Base.OneTo(2), Tuple{Int64, Int64}, 2, Tuple{Nothing, Nothing}, Tuple{UnitRange{Int64}, StepRange{Int64, Int64}}}, Tuple{UnitRange{Int64}, StepRange{Int64, Int64}}}}:
 10  20  30
 20  40  60

Both can be extended:

• grid(f, xs...) can support functions, not just types. A potential confusion is that functions and types would be handled a bit differently: functions would just be called like f((x, x, ...)), but types are constructed with constructorof(T)(x, x, ...).
• mapview on a RectiGrid (and any other KeyedArray btw) can return a KeyedArray back, instead of a generic MappedArray. Here, we have the common Julian problem of nested array wrappers. I'm not totally sure whether it's more intuitive for mapview(f, array) to always return MappedArray, or to rewrap them into other array types sometimes.

For now, added "KeyedArray propagation" to FlexiMaps. So, mapview(splat(*), grid(1:2, 10:10:30)) is still a MappedArray (not KeyedArray), but supports the KeyedArray API:

julia> A = mapview(splat(*), grid(1:5, 10:10:100))
5×10 FlexiMaps.MappedArray{Int64, 2, Base.Splat{typeof(*)}, KeyedArray{Tuple{Int64, Int64}, 2, RectiGrids.RectiGridArr{Base.OneTo(2), Tuple{Int64, Int64}, 2, Tuple{Nothing, Nothing}, Tuple{UnitRange{Int64}, StepRange{Int64, Int64}}}, Tuple{UnitRange{Int64}, StepRange{Int64, Int64}}}}:
 10   20   30   40   50   60   70   80   90  100
 20   40   60   80  100  120  140  160  180  200
 30   60   90  120  150  180  210  240  270  300
 40   80  120  160  200  240  280  320  360  400
 50  100  150  200  250  300  350  400  450  500

julia> axiskeys(A)
(1:5, 10:10:100)

julia> A(5, 50)
250

julia> A(5, 50..100)
6-element view(::FlexiMaps.MappedArray{Int64, 2, Base.Splat{typeof(*)}, KeyedArray{Tuple{Int64, Int64}, 2, RectiGrids.RectiGridArr{Base.OneTo(2), Tuple{Int64, Int64}, 2, Tuple{Nothing, Nothing}, Tuple{UnitRange{Int64}, StepRange{Int64, Int64}}}, Tuple{UnitRange{Int64}, StepRange{Int64, Int64}}}}, 5, 5:10) with eltype Int64:
 250
 300
 350
 400
 450
 500

Do you have a specific usecase in mind? This can help understanding what the ideal approach would be.
I almost never reach for something like this myself.

[jariji]

grid(f, xs...) can support functions, not just types.

I don't like this because f might be iterable.

mapview(T, grid(xs...)) in FlexiMaps returns a (lazy) view

I was thinking eager, to start with.

Do you have a specific usecase in mind? This can help understanding what the ideal approach would be.
I almost never reach for something like this myself.

If you have f : A -> B -> C and a : Vector{A} and b : Vector{B} then the two obvious ways to combine them are inner(f, a, b) (aka map/zipWith) and outer(f, a, b). outer is just the array way of writing a double loop. The current grid function is outer(tuple, a, b).

Good demos and analysis are in Outer Product as an Introduction to APL and a Pretty Cool Thing in General.

Examples:

Identity matrix

julia> outer(==, 1:5, 1:5)
2-dimensional KeyedArray(...) with keys:
↓   5-element UnitRange{Int64}
→   5-element UnitRange{Int64}
And data, 5×5 Matrix{Bool}:
      (1)  (2)  (3)  (4)  (5)
 (1)    1    0    0    0    0
 (2)    0    1    0    0    0
 (3)    0    0    1    0    0
 (4)    0    0    0    1    0
 (5)    0    0    0    0    1

Upper triangular

julia> outer(≤, 1:5, 1:5)
2-dimensional KeyedArray(...) with keys:
↓   5-element UnitRange{Int64}
→   5-element UnitRange{Int64}
And data, 5×5 Matrix{Bool}:
      (1)  (2)  (3)  (4)  (5)
 (1)    1    1    1    1    1
 (2)    0    1    1    1    1
 (3)    0    0    1    1    1
 (4)    0    0    0    1    1
 (5)    0    0    0    0    1

Multiplication table

julia> outer(*, 1:5, 1:5)
2-dimensional KeyedArray(...) with keys:
↓   5-element UnitRange{Int64}
→   5-element UnitRange{Int64}
And data, 5×5 Matrix{Int64}:
      (1)  (2)  (3)  (4)  (5)
 (1)    1    2    3    4    5
 (2)    2    4    6    8   10
 (3)    3    6    9   12   15
 (4)    4    8   12   16   20
 (5)    5   10   15   20   25

Menu

julia> meat = ["chicken", "pork", "beef"]
       style = ["fried rice", "lo mein", "broccoli", "with black bean sauce"]
       outer((m,s)->"$m $s", meat, style)
2-dimensional KeyedArray(...) with keys:
↓   3-element Vector{String}
→   4-element Vector{String}
And data, 3×4 Matrix{String}:
               ("fried rice")          ("lo mein")          ("broccoli")          ("with black bean sauce")
  ("chicken")    "chicken fried rice"    "chicken lo mein"    "chicken broccoli"    "chicken with black bean sauce"
  ("pork")       "pork fried rice"       "pork lo mein"       "pork broccoli"       "pork with black bean sauce"
  ("beef")       "beef fried rice"       "beef lo mein"       "beef broccoli"       "beef with black bean sauce"

Pascal's triangle

julia> outer(binomial, 0:12, 0:12)
2-dimensional KeyedArray(...) with keys:
↓   13-element Vector{Int64}
→   13-element Vector{Int64}
And data, 13×13 Matrix{Int64}:
       (0)  (1)  (2)  (3)  (4)  (5)  (6)  (7)  (8)  (9)  (10)  (11)  (12)
  (0)    1    0    0    0    0    0    0    0    0    0     0     0     0
  (1)    1    1    0    0    0    0    0    0    0    0     0     0     0
  (2)    1    2    1    0    0    0    0    0    0    0     0     0     0
  (3)    1    3    3    1    0    0    0    0    0    0     0     0     0
  (4)    1    4    6    4    1    0    0    0    0    0     0     0     0
  (5)    1    5   10   10    5    1    0    0    0    0     0     0     0
  (6)    1    6   15   20   15    6    1    0    0    0     0     0     0
  (7)    1    7   21   35   35   21    7    1    0    0     0     0     0
  (8)    1    8   28   56   70   56   28    8    1    0     0     0     0
  (9)    1    9   36   84  126  126   84   36    9    1     0     0     0
 (10)    1   10   45  120  210  252  210  120   45   10     1     0     0
 (11)    1   11   55  165  330  462  462  330  165   55    11     1     0
 (12)    1   12   66  220  495  792  924  792  495  220    66    12     1

Preserves associativity of f

julia> outer(*, outer(*, 'a':'d', 'A':'D'), 'α':'δ') == outer(*, 'a':'d', outer(*, 'A':'D'), 'α':'δ') # true 

[aplavin]

I was thinking eager, to start with.

Then map(f, grid(...)) as you say, it's the straightforward solution.

Like this, but without the intermediate allocation.

It should only allocate the final result of map.

[jariji]

I think that isn't associative for associative f like *.

julia> outer(f, a, b) = map(splat(f), grid(a,b));

julia> outer(*, 1:3, 10:10:30)
2-dimensional KeyedArray(...) with keys:
↓   3-element UnitRange{Int64}
→   3-element StepRange{Int64,...}
And data, 3×3 Matrix{Int64}:
      (10)  (20)  (30)
 (1)    10    20    30
 (2)    20    40    60
 (3)    30    60    90

julia> outer(*, outer(*, 1:3, 10:10:30), 100:100:300)
ERROR: MethodError: no method matching grid(::KeyedArray{Int64, 2, Matrix{Int64}, Tuple{UnitRange{Int64}, StepRange{Int64, Int64}}}, ::StepRange{Int64, Int64})

[aplavin]

Ah, I see what you mean! For now only 1d arrays are supported as grid components.
Potentially, nd are also in scope, and I'll support a PR adding those. Semantics aren't totally clear for now: for example, how the result should look like for named grid components?

[aplavin]

This works

julia> outer(*, outer(*, 1:3, 10:10:30) |> vec, 100:100:300)
2-dimensional KeyedArray(...) with keys:
↓   9-element Vector{Int64}
→   3-element StepRange{Int64,...}
And data, 9×3 Matrix{Int64}:
       (100)  (200)  (300)
 (10)   1000   2000   3000
 (20)   2000   4000   6000
 (30)   3000   6000   9000
 (20)   2000   4000   6000
 (40)   4000   8000  12000
 (60)   6000  12000  18000
 (30)   3000   6000   9000
 (60)   6000  12000  18000
 (90)   9000  18000  27000

but probably you are looking for 3d array, right?

[jariji]

but probably you are looking for 3d array, right?

Right.

[jariji]

julia> using AxisKeys

julia> outer(f, a, b) =
           KeyedArray(collect(Iterators.map(splat(f), Iterators.product(a, b)));
               NamedTuple(dimnames(a) .=> axiskeys(a))...,
               NamedTuple(dimnames(b) .=> axiskeys(b))...,
           )
outer (generic function with 1 method)

julia> ka1 = KeyedArray(1:2, i='a':'b')
1-dimensional KeyedArray(NamedDimsArray(...)) with keys:
↓   i ∈ 2-element StepRange{Char,...}
And data, 2-element UnitRange{Int64}:
 ('a')  1
 ('b')  2

julia> ka2 = KeyedArray(3:5, j='A':'C')
1-dimensional KeyedArray(NamedDimsArray(...)) with keys:
↓   j ∈ 3-element StepRange{Char,...}
And data, 3-element UnitRange{Int64}:
 ('A')  3
 ('B')  4
 ('C')  5

julia> ka3 = KeyedArray(6:9, k='α':'δ')
1-dimensional KeyedArray(NamedDimsArray(...)) with keys:
↓   k ∈ 4-element StepRange{Char,...}
And data, 4-element UnitRange{Int64}:
 ('α')  6
 ('β')  7
 ('γ')  8
 ('δ')  9

julia> @assert outer(*, outer(*, ka1, ka2), ka3) == outer(*, ka1, outer(*, ka2, ka3))

julia> outer(*, outer(*, ka1, ka2), ka3)
3-dimensional KeyedArray(NamedDimsArray(...)) with keys:
↓   i ∈ 2-element StepRange{Char,...}
→   j ∈ 3-element StepRange{Char,...}
◪   k ∈ 4-element StepRange{Char,...}
And data, 2×3×4 Array{Int64, 3}:
[:, :, 1] ~ (:, :, 'α'):
         ('A')  ('B')  ('C')
  ('a')  18     24     30
  ('b')  36     48     60

[:, :, 2] ~ (:, :, 'β'):
         ('A')  ('B')  ('C')
  ('a')  21     28     35
  ('b')  42     56     70

[:, :, 3] ~ (:, :, 'γ'):
         ('A')  ('B')  ('C')
  ('a')  24     32     40
  ('b')  48     64     80

[:, :, 4] ~ (:, :, 'δ'):
         ('A')  ('B')  ('C')
  ('a')  27     36     45
  ('b')  54     72     90

julia> outer(vcat, outer(vcat, ka1, ka2), ka3)
3-dimensional KeyedArray(NamedDimsArray(...)) with keys:
↓   i ∈ 2-element StepRange{Char,...}
→   j ∈ 3-element StepRange{Char,...}
◪   k ∈ 4-element StepRange{Char,...}
And data, 2×3×4 Array{Vector{Int64}, 3}:
[:, :, 1] ~ (:, :, 'α'):
         ('A')        ('B')        ('C')
  ('a')    [1, 3, 6]    [1, 4, 6]    [1, 5, 6]
  ('b')    [2, 3, 6]    [2, 4, 6]    [2, 5, 6]

[:, :, 2] ~ (:, :, 'β'):
         ('A')        ('B')        ('C')
  ('a')    [1, 3, 7]    [1, 4, 7]    [1, 5, 7]
  ('b')    [2, 3, 7]    [2, 4, 7]    [2, 5, 7]

[:, :, 3] ~ (:, :, 'γ'):
         ('A')        ('B')        ('C')
  ('a')    [1, 3, 8]    [1, 4, 8]    [1, 5, 8]
  ('b')    [2, 3, 8]    [2, 4, 8]    [2, 5, 8]

[:, :, 4] ~ (:, :, 'δ'):
         ('A')        ('B')        ('C')
  ('a')    [1, 3, 9]    [1, 4, 9]    [1, 5, 9]
  ('b')    [2, 3, 9]    [2, 4, 9]    [2, 5, 9]