Can you broadcast functions over arbitrary rank? (like J/APL)

Inkq · August 25, 2019, 9:38am

Hello! I’ve been eyeing Julia for a while as a J/APL fan and really want to move some programming I would otherwise do in J into Julia, given they are both array-centric programming languages.

But one of the things that keeps me sticking to J is the ability to easily “broadcast” over arbitrary rank. As an example, if I have a 2 x 10 array a:

    a =. 2 10 $ 1 + i. 20

 1  2  3  4  5  6  7  8  9 10
11 12 13 14 15 16 17 18 19 20

And a 10 x 10 array b:

    b =. 10 10 $ 100 + i. 100

100 101 102 103 104 105 106 107 108 109
110 111 112 113 114 115 116 117 118 119
120 121 122 123 124 125 126 127 128 129
130 131 132 133 134 135 136 137 138 139
140 141 142 143 144 145 146 147 148 149
150 151 152 153 154 155 156 157 158 159
160 161 162 163 164 165 166 167 168 169
170 171 172 173 174 175 176 177 178 179
180 181 182 183 184 185 186 187 188 189
190 191 192 193 194 195 196 197 198 199

Then in J, I can quickly add each row of a to each row of b to produce a 2x10x10 array using J’s rank operator (adding vectors ("1), for each vector in a over b ("1 2)):

   a +"1"1 2 b
101 103 105 107 109 111 113 115 117 119
111 113 115 117 119 121 123 125 127 129
121 123 125 127 129 131 133 135 137 139
131 133 135 137 139 141 143 145 147 149
141 143 145 147 149 151 153 155 157 159
151 153 155 157 159 161 163 165 167 169
161 163 165 167 169 171 173 175 177 179
171 173 175 177 179 181 183 185 187 189
181 183 185 187 189 191 193 195 197 199
191 193 195 197 199 201 203 205 207 209

111 113 115 117 119 121 123 125 127 129
121 123 125 127 129 131 133 135 137 139
131 133 135 137 139 141 143 145 147 149
141 143 145 147 149 151 153 155 157 159
151 153 155 157 159 161 163 165 167 169
161 163 165 167 169 171 173 175 177 179
171 173 175 177 179 181 183 185 187 189
181 183 185 187 189 191 193 195 197 199
191 193 195 197 199 201 203 205 207 209
201 203 205 207 209 211 213 215 217 219

This is a contrived example, but this comes in extremely handy in array programming with J.

In Julia, does a similarly expressive way to do that operation exist, i.e. to broadcast a function over arbitrary specified ranks? I can make the arrays easily enough for this example, but I can’t find a way to add them like above in Julia without writing an explicit for loop. Is there something expressive that would accomplish this, without writing explicit for loops?

Dot syntax works for me for the simple case in which only one dimension needs to be broadcasted, but I’m looking for something as flexible as J’s rank operator.

For reference, these are the same arrays a and b in Julia - I just don’t know how to add them up to produce a 3 dimensional array as above in an easily expressible way, without explicitly looping.

a = collect(1:20) |> x -> reshape(x, (10, 2)) |> transpose |> copy;
b = collect(100:199) |> x -> reshape(x, (10,10)) |> transpose |> copy;

Thank you in advance for any help! I really like Julia so far, and I really want to use it. But I’m having trouble letting go of some of the indulgences that J/APL give users in manipulating arrays easily and expressively.

antoine-levitt · August 25, 2019, 10:51am

I think the easiest 1-liner is a comprehension: [a[i,j]+b[j,k] for i =1:2, j=1:10, k=1:10] which I find very nice and readable for this kind of manipulations (often more than complex broadcasts where I instantly get confused about which dimension is which). It does sometimes suffer from the infamous performance of captured variables in closures · Issue #15276 · JuliaLang/julia · GitHub, though, so careful about its usage in performance-critical code.

Inkq · August 25, 2019, 12:07pm

Thank you - I think that’s an elegant solution. I quite like list comprehensions, but I have to admit I couldn’t think of how to write one for this problem so this is very helpful.

I know broadcasting can get very fiddly, but I have to say I really like how J deals with rank - experience with it has made it feel very intuitive. Definitely wish there was an equivalent in Julia, but perhaps I’ll eventually reach a stage where I don’t need it.

baggepinnen · August 25, 2019, 12:12pm

You might be able to write yourself a macro that accepts a syntax you choose and performs the J functions using list comprehensions

bkamins · August 25, 2019, 12:16pm

For this problem I would use broadcasting in Julia in the following way:

julia> a = reshape(1:20, 1, 10, 2)
1×10×2 reshape(::UnitRange{Int64}, 1, 10, 2) with eltype Int64:
[:, :, 1] =
 1  2  3  4  5  6  7  8  9  10

[:, :, 2] =
 11  12  13  14  15  16  17  18  19  20

julia> b = transpose(reshape(100:199, 10, 10))
10×10 LinearAlgebra.Transpose{Int64,Base.ReshapedArray{Int64,2,UnitRange{Int64},Tuple{}}}:
 100  101  102  103  104  105  106  107  108  109
 110  111  112  113  114  115  116  117  118  119
 120  121  122  123  124  125  126  127  128  129
 130  131  132  133  134  135  136  137  138  139
 140  141  142  143  144  145  146  147  148  149
 150  151  152  153  154  155  156  157  158  159
 160  161  162  163  164  165  166  167  168  169
 170  171  172  173  174  175  176  177  178  179
 180  181  182  183  184  185  186  187  188  189
 190  191  192  193  194  195  196  197  198  199

julia> a .+ b
10×10×2 Array{Int64,3}:
[:, :, 1] =
 101  103  105  107  109  111  113  115  117  119
 111  113  115  117  119  121  123  125  127  129
 121  123  125  127  129  131  133  135  137  139
 131  133  135  137  139  141  143  145  147  149
 141  143  145  147  149  151  153  155  157  159
 151  153  155  157  159  161  163  165  167  169
 161  163  165  167  169  171  173  175  177  179
 171  173  175  177  179  181  183  185  187  189
 181  183  185  187  189  191  193  195  197  199
 191  193  195  197  199  201  203  205  207  209

[:, :, 2] =
 111  113  115  117  119  121  123  125  127  129
 121  123  125  127  129  131  133  135  137  139
 131  133  135  137  139  141  143  145  147  149
 141  143  145  147  149  151  153  155  157  159
 151  153  155  157  159  161  163  165  167  169
 161  163  165  167  169  171  173  175  177  179
 171  173  175  177  179  181  183  185  187  189
 181  183  185  187  189  191  193  195  197  199
 191  193  195  197  199  201  203  205  207  209
 201  203  205  207  209  211  213  215  217  219

EDIT

Now I see that you want the first dimension to be 2 (it was natural for me that the last dimension is 2 for some reason I do not understand ), so you can write: permutedims(a .+ b, (3, 1, 2)) (I reuse my variables; you could also define a and b differently to get the desired effect in one shot).

The accepted answer (using your variables) should be [a[i,k]+b[j,k] for i =1:2, j=1:10, k=1:10] I think (note that I use different indexing).

In general with reshape, permutedims functions and standard Julia broadcasting you should be able to achieve most transformations you expect I think.

Inkq · August 26, 2019, 12:23am

Ahh that’s really awesome, I didn’t know those shapes would be compatible! The order of the dimensions wasn’t so important to me as knowing whether this could be expressed in Julia, so I really like this solution. And I think having the 3D structure expressed in one of the argument arrays is also quite intuitive. Thanks!