At dot of function of five arguments against five vectors to five dimensional array

I’ve seen many examples that generate a two dimensional array (matrix) from a function of two arguments against two vectors:

x = 1:100
y = 1:100
f(x,y) = x*y
z = @. f(x', y)

I want to do something analagous, but five-dimensional:

v = 1:100
w = 1:100
x = 1:100
y = 1:100
z = 1:100
g(v,w,x,y,z) = v*w*x*y*z
u = @. g(reshape(v,(1,1,1,1,100)),reshape(w,(1,1,1,100,1)),reshape(x,(1,1,100,1,1)),reshape(y,(1,100,1,1,1,1)),reshape(z,(100,1,1,1,1)))

But when I try to do the latter I get

ERROR: DimensionMismatch: arrays could not be broadcast to a common size; got a dimension with lengths 100 and 5

(I also tried nonsense like v'''' but that just results in v again because it limits itself to 2 dimensions)

Is there a good way to do what I’m trying to do?

@. is too powerful for what you want to do, it fails to broadcast the reshape calls. Try

u = g.(reshape(v,(1,1,1,1,100)),reshape(w,(1,1,1,100,1)),reshape(x,(1,1,100,1,1)),reshape(y,(1,100,1,1,1,1)),reshape(z,(100,1,1,1,1)))

instead.

That gives me an out of memory error, which is promising! Thanks! (Same as I’d get for collect(1:100^5) which makes sense)

1 Like

Perhaps adding a utility function can make the expression simpler. Here is an example:

alongdim(v,d,D) = reshape(v,ntuple(t->t==d ? D : 1))
alongdim(v,d) = alongdim(v,d,d)

Using alongdim in this context:

julia> u = (sin∘prod∘vcat).(alongdim.(Ref(1:2),1:4)...)
2×2×2×2 Array{Float64, 4}:
[:, :, 1, 1] =
 0.841471   0.909297
 0.909297  -0.756802

[:, :, 2, 1] =
  0.909297  -0.756802
 -0.756802   0.989358

[:, :, 1, 2] =
  0.909297  -0.756802
 -0.756802   0.989358

[:, :, 2, 2] =
 -0.756802   0.989358
  0.989358  -0.287903

alongdim(v,d,D) would return an array of dimension D with the vector v laid out a long dimension d.

u = g.(reshape(v,(1,1,1,1,100)),reshape(w,(1,1,1,100,1)),reshape(x,(1,1,100,1,1)),reshape(y,(1,100,1,1,1,1)),reshape(z,(100,1,1,1,1)))

can be rewritten as:

u = g.(
  alongdim(v,5,5),
  alongdim(w,4,5),
  alongdim(x,3,5), 
  alongdim(y,2,5), 
  alongdim(z,1,5)
)

or even shorter:

u = g.(alongdim.((v,w,x,y,z),5:-1:1,5)...)
1 Like
Base.splat(g).(Iterators.product(v,w,x,y,z))
g.([v;],[w...;;],[x...;;;],[y...;;;;],[z...;;;;;])

a=[v,w,x,y,z]

args=ntuple(i->cat(a[i]...,dims=i),5)

g.(args...)
nn=ntuple(i->ntuple(j->j==i ? (:) : 1,5),5)
g.(reshape.(a,nn)...)
1 Like
u = [g(V, W, X, Y, Z) for Z in z, Y in y, X in x, W in w, V in v]
2 Likes

Writing this using TensorCast.jl is compact and has the merit of clarity in tracking indices:

using TensorCast
@cast u[i,j,k,l,m] := g(v[i], w[j], x[k], y[l], z[m])
3 Likes