Suppose in Julia you have a function which outputs a vector, like this:
foo(x,y) = rand(3)
And you broadcast it over some arrays
xs = [i for i=1:4, j=1:4]
ys = [j for i=1:4, j=1:4]
A = foo.(xs,ys)
Now A is an array of arrays. How can I turn this into a 3 dimensional array? More generally, is there a solution for higher dimensional situations? E.g. for a similar function foo that outputs an n dimensional array, broadcast over m dimensional arrays xs and ys, how can one convert the broadcasted output to an n m dimensional array?
I've tried this solution, but vcat and hcat only work on one dimension. I've tried the Destruct.jl package, but again it only works on one dimension. I've tried all the obvious indexing tricks I can think of, but e.g. (foo.(xs,ys))[:,:][1] just returns (foo.(xs,ys))[1], which is not what I want.
CodePudding user response：
You can use SplitApplyCombine.jl like this:
julia> using SplitApplyCombine
julia> A
4×4 Matrix{Vector{Float64}}:
[0.718738, 0.911636, 0.0113524] [0.356858, 0.454262, 0.599563] … [0.455822, 0.475587, 0.567987]
[0.126141, 0.435478, 0.854416] [0.621569, 0.210589, 0.404811] [0.682072, 0.245013, 0.962997]
[0.361148, 0.753464, 0.752801] [0.0643847, 0.193872, 0.820067] [0.676852, 0.942302, 0.093561]
[0.0223175, 0.0972373, 0.215464] [0.737345, 0.737404, 0.996896] [0.183609, 0.335617, 0.720999]
julia> combinedims(A)
3×4×4 Array{Float64, 3}:
[:, :, 1] =
0.718738 0.126141 0.361148 0.0223175
0.911636 0.435478 0.753464 0.0972373
0.0113524 0.854416 0.752801 0.215464
[:, :, 2] =
0.356858 0.621569 0.0643847 0.737345
0.454262 0.210589 0.193872 0.737404
0.599563 0.404811 0.820067 0.996896
[:, :, 3] =
0.135787 0.851451 0.930372 0.498012
0.74397 0.552808 0.960158 0.506592
0.502612 0.0145137 0.915655 0.538791
[:, :, 4] =
0.455822 0.682072 0.676852 0.183609
0.475587 0.245013 0.942302 0.335617
0.567987 0.962997 0.093561 0.720999
Is this what you wanted to get?
CodePudding user response：
I assume that you want to concatenate such that for the result B, we have
B[i,j,:] == A[i,j]
Here's how I'd write this by hand, with a couple of test cases:
julia> function stack_inner(A)
m = ndims(A)
n = ndims(eltype(A))
s_outer = size(A)
s_inner = size(A[begin])
T = eltype(eltype(A))
B = similar(A, T, s_outer..., s_inner...)
.. = fill(:, n)
for ix in CartesianIndices(A)
view(B, ix, (..)...) .= A[ix]
end
return B
end
stack_inner (generic function with 1 method)
julia> stack_inner(foo.(xs, ys))
4×4×3 Array{Float64, 3}:
[:, :, 1] =
0.146902 0.454526 0.191392 0.0322454
0.824091 0.482875 0.700646 0.431301
0.701653 0.0762824 0.194861 0.62421
0.663212 0.0853607 0.313588 0.368867
[:, :, 2] =
0.828077 0.74424 0.255279 0.666977
0.854976 0.302373 0.649691 0.41975
0.759374 0.200208 0.502984 0.886694
0.315596 0.683564 0.956973 0.170769
[:, :, 3] =
0.765305 0.327369 0.824123 0.0537041
0.642428 0.595402 0.235029 0.53
0.603463 0.87867 0.913007 0.548221
0.058201 0.0320288 0.636 0.39045
julia> A = foo.(xs,ys)
4×4 Matrix{Vector{Float64}}:
[0.369792, 0.692929, 0.330885] [0.986054, 0.628871, 0.604634] [0.367734, 0.974091, 0.621425] [0.848115, 0.76681, 0.070687]
[0.00585058, 0.0253985, 0.0470831] [0.10664, 0.373489, 0.111656] [0.719105, 0.126512, 0.660547] [0.999209, 0.0836153, 0.56231]
[0.88527, 0.745378, 0.380452] [0.861579, 0.252228, 0.303043] [0.506468, 0.645717, 0.443472] [0.322553, 0.80937, 0.90342]
[0.783752, 0.553846, 0.830212] [0.868647, 0.0431845, 0.868717] [0.533789, 0.247143, 0.968839] [0.813371, 0.78052, 0.0166259]
julia> stack_inner(A)
4×4×3 Array{Float64, 3}:
[:, :, 1] =
0.369792 0.986054 0.367734 0.848115
0.00585058 0.10664 0.719105 0.999209
0.88527 0.861579 0.506468 0.322553
0.783752 0.868647 0.533789 0.813371
[:, :, 2] =
0.692929 0.628871 0.974091 0.76681
0.0253985 0.373489 0.126512 0.0836153
0.745378 0.252228 0.645717 0.80937
0.553846 0.0431845 0.247143 0.78052
[:, :, 3] =
0.330885 0.604634 0.621425 0.070687
0.0470831 0.111656 0.660547 0.56231
0.380452 0.303043 0.443472 0.90342
0.830212 0.868717 0.968839 0.0166259
julia> xs = rand(2,2,3); ys = rand(2,2,3);
julia> A = foo.(xs, ys)
2×2×3 Array{Vector{Float64}, 3}:
[:, :, 1] =
[0.992239, 0.0603648, 0.687885] [0.896318, 0.825111, 0.31629]
[0.348628, 0.580808, 0.365436] [0.31947, 0.126727, 0.364692]
[:, :, 2] =
[0.714268, 0.0538692, 0.404262] [0.1527, 0.556172, 0.922746]
[0.521115, 0.383689, 0.731707] [0.663383, 0.764024, 0.61838]
[:, :, 3] =
[0.23384, 0.388472, 0.413886] [0.119036, 0.117612, 0.365978]
[0.208551, 0.875924, 0.783887] [0.444344, 0.466899, 0.523953]
julia> stack_inner(A)
2×2×3×3 Array{Float64, 4}:
[:, :, 1, 1] =
0.992239 0.896318
0.348628 0.31947
[:, :, 2, 1] =
0.714268 0.1527
0.521115 0.663383
[:, :, 3, 1] =
0.23384 0.119036
0.208551 0.444344
[:, :, 1, 2] =
0.0603648 0.825111
0.580808 0.126727
[:, :, 2, 2] =
0.0538692 0.556172
0.383689 0.764024
[:, :, 3, 2] =
0.388472 0.117612
0.875924 0.466899
[:, :, 1, 3] =
0.687885 0.31629
0.365436 0.364692
[:, :, 2, 3] =
0.404262 0.922746
0.731707 0.61838
[:, :, 3, 3] =
0.413886 0.365978
0.783887 0.523953
julia> foo(x,y) = rand(3, 3)
foo (generic function with 1 method)
julia> xs = rand(1,2); ys = rand(1,2);
julia> A = foo.(xs, ys)
1×2 Matrix{Matrix{Float64}}:
[0.886281 0.898241 0.0377659; 0.720996 0.401784 0.567878; 0.670619 0.681678 0.457421] [0.0157844 0.23003 0.155695; 0.698028 0.535197 0.458348; 0.500317 0.880123 0.970028]
julia> stack_inner(A)
1×2×3×3 Array{Float64, 4}:
[:, :, 1, 1] =
0.886281 0.0157844
[:, :, 2, 1] =
0.720996 0.698028
[:, :, 3, 1] =
0.670619 0.500317
[:, :, 1, 2] =
0.898241 0.23003
[:, :, 2, 2] =
0.401784 0.535197
[:, :, 3, 2] =
0.681678 0.880123
[:, :, 1, 3] =
0.0377659 0.155695
[:, :, 2, 3] =
0.567878 0.458348
[:, :, 3, 3] =
0.457421 0.970028
(The funny .. is just a hack to do what one could use EllipsisNotation.jl for.)