I often need gridded versions of continuous functions. For a function with one argument, the usual broad cast is the most elegant way:

``````oned(x) = 2*x # 1D func
xs = 0.0:1.0:10.0 #  x grid
oned_gridded = oned.(xs)
``````

This method is often found in Julia examples.

Then, what about multiple-argument functions? You can of course use the list comprehension, but it’s a bit more verbose and a little bit more error-prone: `[twod(x,y) for x in xs, y in ys]`.

and came up with

``````twod(x,y) = "\$(x) & \$(y)" # 2D func
twod_tuple(tpl) = twod(tpl, tpl)
xs = 0.0:1.0:10.0 #  x grid
ys = 20.0:3.0:50.0 # y grid
twod_gridded = twod_tuple.(Iterators.product(xs,ys))
``````

Is this the most succinct solution today? Is there a standard way to convert a multiple-argument function into a single-tuple-argument function?

By the way, I’ve never been able to remember whether `x in xs, y in ys` or `y in ys, x in xs` is the correct order. That’s one of the reasons why I want to avoid the list comprehension to construct a multidimensional array. For the regular `for` loop,

``````for j in axes(arr,2)
for i in axes(arr,1)
arr[i,j] = func_of(i,j)
``````

is the right order because `i` should change faster. The above is (functionally, at least) equivalent to

``````for j in axes(arr,2), i in axes(arr,1)
arr[i,j] = func_of(i,j)
``````

So far so good. But then, you have to flip the order for the list comprehension

``````arr = [func_of(i,j) for i in is, j in js]
``````

if I’m not mistaken.

In contrast, there is no such doubt in `Iterator.product` because it is obviously the standard outer product in linear algebra.

For this reason, I’d probably abolish the regular `for` loop from my code and write it this way:

``````for (i,j) in Iterators.product(is,js)
arr[i,j] = func_of(i,j)
``````

After all, a nested `for` loop is an outer product.

Then I hope the above loop can be written like

``````for (i,j) in is ⊗ js # outer product
``````
1 Like

You can broadcast over multiple arguments and dimensions (I believe the term actually implies multiple inputs), so you can do

``````twod.(x, y')
``````

Notice the `'`. Now the dimensions of length > 1 will be broadcast (expanded) over the corresponding dimensions in the other argument. So, broadcasting over an Nx1 and 1xM array yields an NxM array. This generalizes to higher dimensions and more inputs.

The broadcasting/dimension expanding behavior is crucial, and distinguishes it from `map`.

4 Likes

That’s indeed what I have been missing! Thanks!

But

This generalizes to higher dimensions and more inputs.

I’d like to know how.

I’m relatively familiar with matrices and so, the usual linear-algebraic operation x y^T, which in julia is `x y'`, makes perfect sense as an argument to the broadcast function application.

But, beyond 2 dimensions, I’m clueless. In my limited experience with tensors, I’ve seen only this type of notation a_{i,j,k} = u_i v_j w_k, which is an analogue of the list comprehension.

To apply the `x y'` notation to 3-dimensions, perhaps we would first construct `xy = x y'` and then serialize (1-dimensionalize) `xy` into one dimensional column vector and then `xyz = xy_1d z'` and then somehow recover the first two dimensions. . . .

[Edit:] I’ve just been able to implement the idea:

``````threed(x,y,z) = "\$(x),\$(y),\$(z)" # 3D func
tuplify(f) = tpl -> f(tpl...)
xs = 1:4 # x grid
ys = 10:10:30 # y grid
zs = 100:100:200 # z grid
gridded1 = (tuplify(threed)).(Iterators.product(xs,ys,zs)) # my original idea.
gridded2 = threed.(reshape(xs, (size(xs,1),1,1) ),
reshape(ys, (1, size(ys,1),1)),
reshape(zs, (1, 1, size(zs,1))) )
``````

(I’ve also improved the way to convert a multi-argument function to a single-tuple function.)
As in the above, you can use `reshape` to convert the vectors into L × 1 × 1 array, 1 × M × 1 array, and 1 × 1 × N array.

Now, I wonder how to do this succinctly.

Broadcasting automatically extends with trailing singleton dimensions and you can use the colon feature of reshape:

``````gridded3 = threed.(xs, reshape(ys, 1, :), reshape(zs, 1, 1, :))
``````

Still somewhat too cumbersome for comfort though.

1 Like

Your “tuplify” function is already in Base and called `splat`.

1 Like

I like RectiGrids.jl for this.

``````julia> G = grid(20:30, [:a, :b, :c])
2-dimensional KeyedArray(...) with keys:
↓   11-element UnitRange{Int64}
→   3-element Vector{Symbol}
And data, 11×3 RectiGrids.RectiGridArr{Base.OneTo(2), Tuple{Int64, Symbol}, 2, Tuple{Nothing, Nothing}, Tuple{UnitRange{Int64}, Vector{Symbol}}}:
(:a)        (:b)        (:c)
(20)    (20, :a)    (20, :b)    (20, :c)
(21)    (21, :a)    (21, :b)    (21, :c)
(22)    (22, :a)    (22, :b)    (22, :c)
(23)    (23, :a)    (23, :b)    (23, :c)
(24)    (24, :a)    (24, :b)    (24, :c)
(25)    (25, :a)    (25, :b)    (25, :c)
(26)    (26, :a)    (26, :b)    (26, :c)
(27)    (27, :a)    (27, :b)    (27, :c)
(28)    (28, :a)    (28, :b)    (28, :c)
(29)    (29, :a)    (29, :b)    (29, :c)
(30)    (30, :a)    (30, :b)    (30, :c)
``````

The syntax may be simplified somewhat by defining a function to perform the reshape:

``````julia> ↑(a, ::Val{N}) where {N} = reshape(a, ntuple(_->1, N-1)..., :)
↑ (generic function with 1 method)

julia> xs = 1:2; ys = 10:10:20; zs = 100:100:200;

julia> threed.(xs, ↑(ys, Val(2)), ↑(zs, Val(3)))
2×2×2 Array{String, 3}:
[:, :, 1] =
"1,10,100"  "1,20,100"
"2,10,100"  "2,20,100"

[:, :, 2] =
"1,10,200"  "1,20,200"
"2,10,200"  "2,20,200"
``````

You can also mimic numpy’s `newaxis` in Julia, which generalizes to more dimensions:

``````const newaxis = [CartesianIndex()]
gridded = threed.(xs, ys[newaxis,:], zs[newaxis,newaxis,:])
``````
2 Likes

Wow! That has made it perfectly clear to me what broadcast is. Before reading the quoted sentence above, it had been quite nebulous to me what broadcast is. Thank you for the great explanation!

Aha! Probably that’s the best method using only core features (broadcast and `reshape`). I didn’t know `reshape` has the same convenience as `view`. (Speaking of that, can we replace `reshape` with `view` in the above? . . . Never mind, I’ll investigate.)

I guess we should advocate this. I’ve been kind of surprised that I keep seeing examples on the Internet like

``````gridded = zeros(nx, ny, nz) # Off topic, but I guess the programmer wrote >>
#  >> this because Array{Float64}(undef, nx, ny, nz) is cumbersome to write.
for k in 1:nz
for j in 1:ny
for i in 1:nx
gridded[i,j,k] = threed(xs[i], ys[j], zs[k])
``````

I’m not saying that explicit looping is necessarily bad. A lot of algorithms are more simply written with indices. But the present example doesn’t belong to that category.