Define basic FillArrays types in Base

[MasonProtter]

As suggested by @mbauman and @dlfivefifty in #35778 (comment) and #35778 (comment), I have copied some basic FillArrays.jl infrastructure into Base with the expectation that FillArrays.jl will pirate this infrastructure.

This way, we can suggest people use Fill(x) rather than Ref(x) in broadcast code, as well as use that internally which will be more transparent to the compiler due to immutability. For example, this should fix #39151.

I don't have the best understanding of FillArrays.jl (though it is pretty simple), so I'm pretty sure the subset I copied here is not an ideal subset. There's probably stuff I copied that I shouldn't have and vice versa.

Closes #18379

[dlfivefifty]

I believe if you export any new types it would count as “Breaking”. So perhaps better to make a user do Base.Fill for now

[MasonProtter]

I believe Base has added new exported types in the past even though they were technically breaking. I think it'd be good to export Fill at least, but let's see what others have to say.

---

I realize now that I didn't properly attribute the code which I copied. Does anyone have guidance on how to do that properly?

[dlfivefifty]

Probably I wrote most of the code and don’t care too much how you attribute. So probably a comment is fine, though double check the terms of the MIT license

[dlfivefifty]

I think it'd be good to export Fill at least

given that Base.OneTo is not exported I think it should be Base.Fill for consistency.

[MasonProtter]

Well, a big part of the motivation for doing this was to replace Ref for broadcast. If it's not exported, it seems weird to recommend it's use.

We need something immutable we can export from Base that people can use instead of Ref. If this isn't acceptable, we should go back to spicing up Some.

[mbauman]

Do we need this? Why not just define this in getindex and use F[] everywhere?

[dlfivefifty]

That seems inconsistent to me unless you meant something else:

julia> F = fill(1,3)
3-element Vector{Int64}:
 1
 1
 1

julia> F[]
ERROR: BoundsError: attempt to access 3-element Vector{Int64} at index []
Stacktrace:
 [1] throw_boundserror(A::Vector{Int64}, I::Tuple{})
   @ Base ./abstractarray.jl:651
 [2] checkbounds
   @ ./abstractarray.jl:616 [inlined]
 [3] _getindex
   @ ./abstractarray.jl:1196 [inlined]
 [4] getindex(::Vector{Int64})
   @ Base ./abstractarray.jl:1170
 [5] top-level scope
   @ REPL[20]:1

[mbauman]

By F[] I just meant indexing in general. Seems like an unnecessary indirection/complexity.

[dlfivefifty]

It's to allow for unified code for Fill, Zeros, and Ones. Using indexing doesn't work here as one would have to add artificial indices.

[mbauman]

whitespace

[mcabbott]

One quirk about adopting Fill for this is that it specialises broadcasting to avoid repeated evaluations, even along dimensions outside of it:

julia> ones(1,10) .+ rand.(Ref(1:99))
1×10 Matrix{Float64}:
 14.0  73.0  68.0  26.0  70.0  27.0  31.0  36.0  43.0  66.0

julia> ones(1,10) .+ rand.(fill(1:99))
1×10 Matrix{Float64}:
 72.0  36.0  74.0  56.0  30.0  50.0  52.0  27.0  30.0  92.0

julia> using FillArrays

julia> ones(1,10) .+ rand.(Fill(1:99))  # all constant
1×10 Matrix{Float64}:
 13.0  13.0  13.0  13.0  13.0  13.0  13.0  13.0  13.0  13.0

It looks like the implementation here does this, too. Is this desirable?

[dlfivefifty]

This was carried over from FillArrays.jl but I think was a temporary restriction. I'm not entirely sure this restriction should be kept.

[MasonProtter]

Should I just delete that line altogether?

[dlfivefifty]

The check is necessary because the function may be incorrect in the "degenerate" case where a is length 1 and so the shape comes from b, not a. Perhaps it would be better to be rewritten as

if length(a) == 1 
   broadcasted(*, a[1], b)
else
   broadcasted(*, a, _broadcast_getindex_value(b))
end

but this would not be type stable.

In any case some thought should be put into this.

[mbauman]

julia> ones(1,10) .+ rand.(Fill(1:99))  # all constant

Oooof, thanks for bringing this up, @mcabbott. That behavior makes sense for FillArrays, but doesn't make sense for a scalarifier :(

[mcabbott]

Even for Fill it strikes me as strange, why deviate from being equivalent to fill?

The one argument I saw was that sparse arrays do something similar. The justification there is that you'd like abs2.(mat) to preserve sparsity, and it does. Broadcasting a function which does not preserve zero gives you a full sparse array, which I think you never ever want, so you're expected to watch out. This means they are further from being a carefree replacement for dense arrays than (I think) Fill wants to be.

Also, what sparse arrays do is not so strong. It seems to only apply along some dimensions which they create in the broadcast, not along others:

using SparseArrays
r99(ignore) = rand(1:99);

# almost like above?
ones(1,10) .+ r99.(spzeros(1))  # all random

# alone
r99.(spzeros(10))     # all constant
r99.(spzeros(10, 10)) # all constant

# broadcasting with special vectors
ones(1,10) .+ r99.(fill("zero", 10))  # all random
ones(1,10) .+ r99.(spzeros(10))       # columns are constant
ones(1,10) .+ r99.(Fill("zero", 10))  # all constant

ones(10,1) .+ r99.(spzeros(10)')      # all random

# ... and when those are extended
ones(10,10) .+ r99.(fill("zero", 10))  # all random
ones(10,10) .+ r99.(spzeros(10))       # all random
ones(10,10) .+ r99.(Fill("zero", 10))  # all constant

# broadcasting with a sparse matrix
ones(1,1) .+ r99.(spzeros(10, 10))     # columns are constant
ones(1,10) .+ r99.(spzeros(10, 10))    # columns are constant
ones(10,1) .+ r99.(spzeros(10, 10))    # all constant

[dlfivefifty]

It was a decision made out of practicality, since this was the only way to allow broadcasting, and yes SparseArrays provided precedence for this decision.

So there are a few options:

1. Live with the surprising behaviour
2. Forgo this optimised broadcasting behaviour (we can always keep around FillArrays.jl for packages that need this, which include InfiniteArrays.jl).
3. Make it possible to recognise if a function is pure at compile time and only have this special broadcasting behaviour in that case. Though I think this would be a substantial change in the compiler so is not realistic.
4. Special case rand and variants to not use the specialised broadcasting

[mcabbott]

It would be nice more generally if broadcasting had some knowledge of purity, as this could be exploited any time a dimension is extended:

julia> @btime (1:1000)' ./ sqrt.(1:20);
  19.291 μs (2 allocations: 156.33 KiB)

julia> @btime (1:1000)' ./ map(sqrt, 1:20);  # 0.1% as many sqrt calls
  3.203 μs (3 allocations: 156.56 KiB)

No idea how hard that would be, though -- for some loop orders there will be a tradeoff with storage, which the Fill case doesn't have. But perhaps whatever the level of exploitation, it could be something functions opt in to, rather than being something the compiler infers? Broadcast.broadcastpure(::typeof(sqrt)) = identity or something.

[mbauman]

So Ranges do some O(1) things like this, but they're explicitly enabled for a select group of operators that are known to be linear and range-preserving, like:

julia/base/broadcast.jl

Line 1124 in 29d5158

broadcasted(::DefaultArrayStyle{1}, ::typeof(+), r::OrdinalRange) = r

The trouble here is that we have a much larger group of functions (pure instead of just linear). Trying to do the same will be a never-ending task (or require that bigger purity system). That said, in many cases, it's as simple as manually hoisting the Fill (at a user-level) to get this behavior. E.g., instead of rand.(Fill(1:99, 10)) it's Fill(rand(1:99), 10). Yes, that doesn't work if you're passing around Fills and such, but I can live with that.

[dlfivefifty]

I don't think the Range comparison is fair since the number of functions where structure is preserved is much narrower than Fill.

[dlfivefifty]

@MasonProtter Can you please provide the reason for closure for future reference?

I'm assuming its because it was not being actively maintained which is fine.

[MasonProtter]

@dlfivefifty, I stopped working on this years ago because it seemed there was a fundamental mismatch between how this Base would want to use this, and how FillArrays.jl would use it.

For this to be a useful immutable replacement for Ref in Base, the different broadcasting behaviour is a non-starter, but that was apparently important for FillArrays.jl to have. So that seemed to imply that Base and FillArrays wouldn't end up 'sharing' this type, but if they end up not sharing it, then adding this to Base seems kinda antagonistic towards FillArrays because the name does need to be exported from Base to be reasonably useful.

All the options laid out in #39184 (comment) seemed worse than just not doing this. And now it's been a few years and the PR is quite stale anyways.