# Howto correctly (hijack) fill an array with formula over indices?
**URL:** https://discourse.julialang.org/t/howto-correctly-hijack-fill-an-array-with-formula-over-indices/25673
**Category:** General Usage
**Tags:** array
**Created:** [June 26, 2019, 2:16am UTC](https://discourse.julialang.org/t/howto-correctly-hijack-fill-an-array-with-formula-over-indices/25673 "2019-06-26T02:16:51Z")
**Posts on this page:** 11
**Page:** 1
### Author:  [@o314](https://discourse.julialang.org/u/o314)
#### Post date: [June 26, 2019, 2:16am UTC](https://discourse.julialang.org/t/howto-correctly-hijack-fill-an-array-with-formula-over-indices/25673/1 "2019-06-26T02:16:51Z")
It seems fill is built to handle static init value only.
I am trying to write something like this
```julia
A = reshape(1:6,(2,3))
fill((i,j)->10i+j,A)
# or, if you prefer
fill(A) do i,j; 10i+j end
```
For now, i have got this one that pass the test
```julia
Base.fill(f::Function, A) = map(x->f(Tuple(x)...),CartesianIndices(A)) # todo fill!
using Test
@test fill((i,j)->10i+j,A) == [11 12 13; 21 22 23]
```
Is there a simpler julian way to do it (especially a builtin in base or stdlib i have missed)?
Is there any pro or con that is better to know with this kind of trick ?
_Of course we could use nested for. But why should we not use high-order function ?_
---
### Author:  [@o314](https://discourse.julialang.org/u/o314)
#### Post date: [June 26, 2019, 3:17am UTC](https://discourse.julialang.org/t/howto-correctly-hijack-fill-an-array-with-formula-over-indices/25673/2 "2019-06-26T03:17:02Z")
_Replying to myself_
I have got this one too (vectorization style)
```julia
Base.fill(f::Function, A) = Base.splat(f).(Tuple.(CartesianIndices(A)))
A = reshape(1:6,(2,3))
@test fill((i,j)->10i+j,A) == [11 12 13; 21 22 23]
```
Not clear yet to me which one is idiomatic …
---
### Author:  [@c42f](https://discourse.julialang.org/u/c42f)
#### Post date: [June 26, 2019, 3:56am UTC](https://discourse.julialang.org/t/howto-correctly-hijack-fill-an-array-with-formula-over-indices/25673/3 "2019-06-26T03:56:52Z")
Perhaps you’re looking for an array comprehension:
```julia
julia> [10i+j for i=1:2, j=1:3]
2×3 Array{Int64,2}:
11 12 13
21 22 23
```
---
### Author:  [@o314](https://discourse.julialang.org/u/o314)
#### Post date: [June 26, 2019, 4:11am UTC](https://discourse.julialang.org/t/howto-correctly-hijack-fill-an-array-with-formula-over-indices/25673/4 "2019-06-26T04:11:17Z")
That’s was the third form i was coming back to append 🙂
A kind of pivoted nested for. I admit it’s very clean for oneliner formula.
However
```julia
[(u=compute_a_long_formula(i,j); v=compute_another_weaponized_math(u); so_finally_got(v)) for i=1:2, j=1:3]
```
can become tricky to read, when
```julia
fill(A) do i,j;
u = compute_a_long_formula(i,j)
v = compute_another_weaponized_math(u)
so_finally_got(v)
end
```
naturally preserve readability.
I was just hoping that was without deserving internals law of good taste and performance.
---
### Author:  [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)
#### Post date: [June 26, 2019, 4:48am UTC](https://discourse.julialang.org/t/howto-correctly-hijack-fill-an-array-with-formula-over-indices/25673/5 "2019-06-26T04:48:32Z")
Don’t define that method for `fill`, it is [type piracy](https://docs.julialang.org/en/latest/manual/style-guide/#Avoid-type-piracy-1).
IMO array comprehension, suggested by @c42f, is your current best option.
There is also `Iterators.product`. See
> [@Outer broadcasting?](https://discourse.julialang.org/t/outer-broadcasting/22701):
>
> Hi, I would like to do something like the [outer product](https://en.wikipedia.org/wiki/Outer_product) for general operations. Basically, I f I have a function f(x1, x2, [...], xn) = [...] and I provide it with n arrays a1 … an as input, then I would like to apply f on all combination of input elements, such that the resulting array has the shape: outer(f, a1, a2, [...], an) |\> size == (size(a1)..., size(a2)..., [...], size(an)...) How can I do that in a clean and performant way? Example: julia\> f(a, b, c) = a+b+c f (generic function …
and for some suggested changes
>
>
> So with broadcasting dots, we've managed to excise nearly all sometimes-scalar/s…ometimes-vectorized functions in Julia to great benefit. There are two notable stragglers: \`getindex\` and \`setindex!\`. Now, since they have their own syntaxes and support so many fancy things it hasn't always been (and really still isn't) completely obvious that it should be deprecated in favor of broadcasting. But they really are just "broadcasting" scalar indexing over the cartesian product of the indices.
>
> \## The advantages
>
> There'd be a number of advantages to spelling nonscalar indexing as some form of broadcasting. First and foremost: fusion. Forget the whole view/copy debate — straight-up fusion could be faster than either. For example, we could define the new syntax \`A.\[I...\]\` to mean \`broadcast((idxs...)-\>A\[idxs...\], I...)\`. Then, to extract the first 10 elements from the first column of \`A\`, you write \`A.\[1:10, 1\]\`. This means that you could fuse in additional operations without any temporaries or views whatsoever with, e.g., \`sqrt.(A.\[1:10, 1\])\`.
>
> The other huge advantage is how this would then generalize these sorts of accesses to \_all\_ data structures, including dictionaries and more. Cf. #24019.
>
> \## The challenges and their potential solutions
>
> Now, there are also some fairly major challenges in making it easy to express all that indexing does with a broadcasting semantic.
> \* \*\*APL Indexing:\*\* The existing \`A\[1:10, 1:10\]\` takes the cartesian product of the two indices, returning a 2-dimensional 10x10 array. I'd expect a broadcasted \`A.\[1:10, 1:10\]\` to use broadcasting's shape-matching semantic and return the diagonal. In general, APL indexing can be thought of as a broadcast where each successive index is "lifted" to a higher dimension than the sum of dimensionalities of all arguments that preceded it. We could potentially have a simple wrapper that flags arguments to be "orthogonalized" before participating in broadcast — options for spelling this wrapper could include \`⟂\` or \`^\`. Thus, \`A\[I, J\]\` becomes \`A.\[I, ^J\]\`. This is a generally useful operation... and for my purposes would completely obviate the pain from the recursive transpose/adjoint fallback removal as you could lift arguments to orthogonal dimensions \_anywhere\_: \`f.(1:10, ^array\_of\_strings)\`. That said, the change in default operation here from APL-like to broadcasting-like may prove to be very painful...
> \* \*\*Index conversions:\*\* The existing indexing API supports \[many types of indices\](https://docs.julialang.org/en/v1/manual/arrays/#man-supported-index-types-1) and even an \[extensible interface for adding more\](https://docs.julialang.org/en/v1/base/arrays/#Base.to\_indices).
> \* For example, \`:\` is actually a \_function\_, but when used as an index it expands out to the entire axis. With broadcast, however, it acts like a function — as a scalar.
> \* Logical indexing with boolean arrays is even worse: there we have an array of trues and falses, but when used as an index it expands out to a vector of the locations of the trues. With broadcast, it behaves just like an array would, yielding the same number of trues and falses as the arrays overall shape.
>
> Resolving this one means unfortunately giving up something fairly major: I don't believe it'll ever be possible to support all these features \_and also\_ support broadcast fusion through to an index computation. So while it would be oh-so-cool to have an expression like \`A.\[clamp.(idx .+ (-N:N), 1, end)\].^2\` (to examine a window about some \`idx\` without worrying about bounds) fuse the whole way down, it simply isn't extensible to the very next thing you'd want to have fuse: \`A.\[A .\> 0\]\`. That would be a nightmare to figure out how to fuse.
> \* \*\*Bounds checking:\*\* In comparison to the other issues, this one feels much simpler, but it's still gonna be a bit of a pain. Non-scalar indexing is able to perform bounds checking at the level of the collections of indices and then perform the scalar indexing with bounds checks off. This is a \_huge\_ win for things like \`:\` (no checks), \`1:10\` (just check endpoints), and logical masks (just check the shape). I think this may be possible to still do — again, if we don't fuse through to the index computations. But I've not completely seen to the end of this tunnel.
> \* \*\*Returned array type:\*\* We’ll need to find a solution to #17533; \`similar\`’s first major use was in defining the output type for nonscalar indexing. Lots of folks specialized it for that reason, and broadcast doesn’t use it.
> \* \*\*Backwards compatibility:\*\* For the most part, this is entirely a new syntax with a straight-forward deprecation. There's one place where I was worried about a clash, though: the \`@.\` macro. That currently leaves indexing expressions alone, but once we introduce the \`A.\[\]\` operator, I'd expect it to dot those bracket operations. It seems like we have sufficient leeway in the macro to do some fancy detection, though, allowing us to insert appropriate handling code if any arguments are arrays and inserting deprecations appropriately.
> \* \*\*Indexed Assignment:\*\* So far I've really only talked about getindex, but all the same applies to \`setindex!\`, too. It actually becomes quite the advantage as the dots can tell the whole story about what's scalar and whats not right where you'd expect (adapted from the sister issue #24086):
> \* \`A\[i\] = x\`: Always sets \`x\` at the scalar index \`i\`.
> \* \`A\[i\] .= X\`: Always broadcasts \`X\` across the element at location \`i\` (always mutating the object in \`A\[i\]\`).
> \* \`A.\[I\] = x\`: Assigns \`x\` to every index in \`I\`.
> \* \`A.\[I\] .= X\`: Always broadcasts \`X\` across the indices selected by \`I\` (always mutating \`A\`).
>
> \## What to do:
>
> So, bringing this all together, I propose we aim to:
>
> \* Introduce the syntax \`A.\[I, J, K\]\` to mean:
> \`\`\`julia
> idxs = to\_indices(A, (I, J, K))
> broadcast((i, j, k)-\>A\[i,j,k\], idxs...)...)
> \`\`\`
> The separation into two statements there is meaningful — that's how it'll behave in fusion. Note, too, that this \_defaults\_ to behaving broadcast-like (and not APL-like). I'm still not entirely sold on it, but I do think it'll unify the language and simplify things. I think this is one of those things we'll have to see how it feels in practice.
> \* Introduce the new "orthogonalization" syntax \`f.(a, ^b, c, ^d)\` to mean:
> \`\`\`julia
> f.(a, Lifted(b, ndims(a)), c, Lifted(d, ndims(a) + ndims(Lifted(b, ndims(a))) + ndims(c))
> \`\`\`
> This may pose some constant-propagation and type-stability challenges, but I hope they're overcome-able.
> \* Apply these semantics to indexed assignments over in #24086.
>
> \## References and prior discussions
>
> This issue brings together lots of iterations of thought that have been spread about through many issues and discourse communications.
> \* \[What to do about multivalue setindex!? (#24086)\](https://github.com/JuliaLang/julia/issues/24086)
> \* \[Broadcast Array Indexing (#19169): \](https://github.com/JuliaLang/julia/issues/19169) discusses if the array itself should participate in the broadcasting (for indexing into many arrays inside an array). \[Discourse: Square Bracket Notation for Broadcasting getindex Across Array of Vectors?\](https://discourse.julialang.org/t/square-bracket-notation-for-broadcasting-getindex-across-array-of-vectors/6227) requests the same thing.
> \* \[Make indexing expressions participate in dot syntax fusion (#22858)\](https://github.com/JuliaLang/julia/issues/22858): started sketching out lots of this, but defaulted to APL indexing semantics. My \[last comment\](https://github.com/JuliaLang/julia/issues/22858#issuecomment-338267870) is fairly negative because lifting arguments to orthogonalize them requires array wrappers — and array wrappers will allocate just like views will. But this omits the fact that many indices won't be heap-allocated (like ranges) and defaulting to broadcasting semantics will further limit the damage here. Upon further reflection I think it's a win.
> \* \[Broadcasting pointwise indexing operator (#2591)\](https://github.com/JuliaLang/julia/issues/2591) is one of the first discussions on this topic.
> \* \[Discourse: understanding view\](https://discourse.julialang.org/t/understanding-view/18286/8?u=mbauman): A discussion about the strange-ness of the sometimes-scalar sometimes-not behaviors of \`getindex\` as compared to \`view\` (which is always nonscalar).
---
### Author:  [@c42f](https://discourse.julialang.org/u/c42f)
#### Post date: [June 26, 2019, 5:11am UTC](https://discourse.julialang.org/t/howto-correctly-hijack-fill-an-array-with-formula-over-indices/25673/6 "2019-06-26T05:11:41Z")
Having `fill(f::Function, dims)` for use with `do` syntax does make some sense. It would technically be a breaking change because it would change the meaning of the following:
```julia
julia> fill(sum, (2,3))
2×3 Array{typeof(sum),2}:
sum sum sum
sum sum sum
```
However, I’m not sure what one could possibly use such an array for (it can only hold objects of the singleton type `typeof(sum)`)!
Keep in mind that in Julia every statement is an expression, so you can create a somewhat readable version of the array comprehension using a block:
```julia
A = [begin
u = compute_a_long_formula(i,j)
v = compute_another_weaponized_math(u)
so_finally_got(v)
end
for i=1:2, j=1:3]
```
Personally I find this kind of odd looking though.
---
### Author:  [@o314](https://discourse.julialang.org/u/o314)
#### Post date: [June 27, 2019, 11:39pm UTC](https://discourse.julialang.org/t/howto-correctly-hijack-fill-an-array-with-formula-over-indices/25673/7 "2019-06-27T23:39:58Z")
> [@Tamas\_Papp](#):
>
> There is also `Iterators.product` .
So one more form
```julia
fil2(f::Function,a::AbstractArray) = Base.splat(f).(Iterators.product(axes(a)...))
@test fil2((i,j)->10i+j, [1 2 3; 4 5 6]) == [11 12 13; 21 22 23]
```
Right. I have forgotten to review the iterators funcs. It’s quite close to CartesianIndex. May be too close to have two.
PRO
CartesianIndex pops here and there (IndexCartesian, CartesianIndex, CartesianIndices, view etc.). Iterators.product seems less used.
CON
Iterators.product may be faster by skipping a call to collect.
* * *
> [@Tamas\_Papp](#):
>
> Don’t define that method for `fill` , it is [type piracy](https://docs.julialang.org/en/latest/manual/style-guide/#Avoid-type-piracy-1).
I disagree. I knew type piracy. It’s a bad practice that had to be fight.
Preventing it by dispatching on func above any, any, any … args in base, is not better.
It’s multidispatch piracy.
There are too much functions foo(x::Any, y::Any) that do rebranch to foo(::Bar, ::Any), foo(::Any, ::Baz), foo(::Qux,::Quux) by bypassing multidispatch via hard-coded spaghetti dispatch.
In very few line, one will inject very hard to grasp, test, optimize, stabilize and check code. This may has been useful to deliver v1 on time. I hope change will come afterward.
* * *
Thanks for the good links. Now i need more time to digest them too 🙂
---
### Author:  [@o314](https://discourse.julialang.org/u/o314)
#### Post date: [June 27, 2019, 11:55pm UTC](https://discourse.julialang.org/t/howto-correctly-hijack-fill-an-array-with-formula-over-indices/25673/8 "2019-06-27T23:55:15Z")
> [@c42f](#):
>
> Having `fill(f::Function, dims)` for use with `do` syntax does make some sense. It would technically be a breaking change because it would change the meaning of the following:
There may be some breaking change but i can not see any major usage of this form too. On the content, Base.fill on function rechecked with compatible type may rebranch to a dynamic fill.
```julia
function Base.fill(f::Function, a::AbstractArray{T,N}) where {T,N} =
isempty(methods(f,NTuple{N,Int})) && error("f") # todo handle IndexStyle
do_fill_with_formula(...)
end
```
* * *
…
> [@c42f](#):
>
> ```julia
> A = [begin
> u = compute_a_long_formula(i,j)
> v = compute_another_weaponized_math(u)
> so_finally_got(v)
> end
> for i=1:2, j=1:3]
>
> ```
is correct, i know.
However what make julia my personal preference is the _“talk like python, run like c, think like lisp”_ motto.
To which i append my personal opinion that is to not loose the _[zen of python](https://www.python.org/dev/peps/pep-0020/)_ in translation.
And too much of those forms hurts it quite a bit imho. so-so.
---
### Author:  [@ElOceanografo](https://discourse.julialang.org/u/ElOceanografo)
#### Post date: [June 27, 2019, 11:59pm UTC](https://discourse.julialang.org/t/howto-correctly-hijack-fill-an-array-with-formula-over-indices/25673/9 "2019-06-27T23:59:50Z")
How about this for zen?
```julia
A = let f(i, j) = so_finally_got(compute_another_weaponized_math(compute_a_long_formula(i,j)))
[f(i, j) for i in 1:2, j in 1:3]
end
```
---
### Author:  [@o314](https://discourse.julialang.org/u/o314)
#### Post date: [June 28, 2019, 12:30am UTC](https://discourse.julialang.org/t/howto-correctly-hijack-fill-an-array-with-formula-over-indices/25673/10 "2019-06-28T00:30:10Z")
A bit hard to read for me. (Readability counts damaged).
3 instruction per lines, nearly identical to
```julia
Base.fill(f::Function, A) = map(x->f(Tuple(x)...),CartesianIndices(A)) # todo fill!
```
Let’s say #ops is ok. What can be improve ?
- Naming - we can shorten them. Those were mine and quite ridiculous.o:
Expert domains bring their own and they [tend to prefer short ones in julia](http://kristofferc.github.io/post/tokenize/).
That may works
- Frequence of use of primitive / patterns. Hoping this topic will help…
…
On a structural viewpoint, I will prefer this one
```julia
A = let f(i, j) =
compute_a_long_formula(i,j) |>
compute_another_weaponized_math |>
so_finally_got;
[f(i,j) for i in 1:2, j in 1:3]
end
```
But i consider julia is better equipped to handle vectorization (v13n) than higher-order functional programming (hours in debugging already lost). So the v13n form with do is what works better for me.
---
### Author:  [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)
#### Post date: [June 28, 2019, 12:01pm UTC](https://discourse.julialang.org/t/howto-correctly-hijack-fill-an-array-with-formula-over-indices/25673/11 "2019-06-28T12:01:31Z")
> [@o314](#):
>
> It’s multidispatch piracy.
Sorry, I don’t understand you.