Is there a reason a `MersenneTwister` value is not broadcastable?

General Usage
1

I have multiple functions with the same pattern:

foo(rng::AbstractRNG, x::ValueType) = ...

which are commonly applied to arrays as

foo.(rng::AbstractRNG, a::Array{ValueType})

which is supposed to be equivalent to, essentially

reshape([foo(rng, x) for x in a], size(a))

rng here is a MersenneTwister-type value, but I guess the question applies to any other kind of RNG. For just the code above, I get an error

MethodError: no method matching length(::MersenneTwister)

In order to make it work I have to wrap a MersenneTwister in my own type and define length() and iterate() for it.

Is there some caveat Iā€™m missing which is the reason these methods are not defined in Random?

In fact, since I need this kind of ā€œbroadcastabilityā€ for many other types, I have the following abstract type defined:

abstract type BroadcastableSingleton end
Base.length(::BroadcastableSingleton) = 1
Base.iterate(x::BroadcastableSingleton) = (x, nothing)
Base.iterate(x::BroadcastableSingleton, state) = nothing

This reduces boilerplate code, but wonā€™t be too clear if there are already some complicated abstract type hierarchies in place. I wonder if there is a better way to do it? (Canā€™t help but miss Haskellā€™s data classes here)

2

A Ref is already designed to do exactly what your custom BroadcastableSingleton does. Wrapping something in a Ref will scalarize it in a broadcasted computation:

foo.(Ref(rng), a)
1 Like
3

Thanks, it does work that way, but:

  1. This call pattern is quite common in my code, and is present in the public API as well. It wonā€™t look good to ask a user to wrap every value that needs this in Ref.

  2. Neither the name nor the docstring for Ref give any indication that it may be used for this purpose. It seems more like an undocumented feature that can change at any time.

4

You are mistaken ā€” it is the standard API, and happens to be well-documented.

https://docs.julialang.org/en/v1/manual/arrays/#Broadcasting-1
https://docs.julialang.org/en/v1/base/arrays/#Base.Broadcast.broadcast

1 Like
5

I would disagree.

  1. ā€˜treats any argument that is not an array ā€¦ as a ā€œscalarā€ā€™ in the first link is incorrect; the second link correctly states that everything but certain predetermined types are iterated over.

  2. Moreover, in the first link it says ā€˜and treats any argument that is not an array, tuple or Ref (except for Ptr) as a ā€œscalarā€ā€™, but the Ref arguments are treated as scalars.

  3. The second link seems more consistent with the actual behaviour. But it only mentions Ref as one of the types that can be broadcasted over by default. There is no indication that it is supposed to be used to create singletons (even though it does create them). Neither, again, there is any indication of that in the docstring of Ref itself.

  4. I donā€™t need all these low-level effects Ref provides (ā€œThis type is guaranteed to point to valid, Julia-allocated memory of the correct type. The underlying data is protected from freeing by the garbage collector as long as the Ref itself is referenced.ā€ and so on). I only need a broadcastable singleton. So in addition to not being the intended usage, it seems like a serious overkill.

  5. My point about this being an unnecessary boilerplate (exaggerated by the fact that from the userā€™s standpoint it seems like some kind of black magic) is still valid.

  6. Edit: One more thing. typeof(Ref(x)) != typeof(x), so I would have to explicitly ask for Ref objects in the function signature.

6

Why donā€™t you just fix the docs a little? The problem seems to be that the docs for Ref are incorrect and not lightweight enough in the introductory paragraphs, not that Ref is not lightweight enough.

And Ref arguments are not treated as scalars. In broadcast Ref is unwrapped and its contents passed to the method. A scalar is itself passed to the method unaltered. You donā€™t have to ask for Ref in the method signature because it isnā€™t treated as a scalar.

1 Like
7

They are not incorrect. I actually donā€™t have anything against the docs for Ref, they seem to contain exactly what they need to contain. It could perhaps have been stated that it is an iterable, but that seems like an implementation detail to me, not related to its purpose.

On the other hand, Single- and multi-dimensional Arrays Ā· The Julia Language seems to contain some incorrect statements (not related to Ref). Although I wouldnā€™t mind if it was correct though, that is if every non-iterable object was treated as a scalar.

Yes, I agree, that is technically correct. Refā€™s contents is treated as a scalar.

Yes, I agree with that as well.

Nevertheless, I still feel like I failed to put the issue across. Instead of Ref(x) I could have used (x,), [x], Set(x) or whatever else built-in iterable container to create a singleton (in fact, Iā€™m surprised that Ref was proposed instead of a tuple - what makes it a better choice? At least with a tuple you can guess that the author of the code wanted to make an iterable). But all of these are clutches. They create noise in the code, they are unintuitive. I think that behaving like a singleton in broadcasting is a natural behaviour for anything that is not an iterable.

8

Sure, I meant either of those docs, some improvements connecting Ref and broadcast. But I also agree that it seems arbitrary which wrapper to use. (,) does seem more intuitive.

Treating everything as an iterable by default could probably be achieved with:

length(x) = 1

but there is probably a reason thatā€™s a bad idea, broadcast has a few more moving parts than thisā€¦

9 
julia> (1,) .+ 2
(3,)

julia> Ref(1) .+ 2
3
1 Like
10

Youā€™d have to define iterate as well. Alternatively, you could broadcast as singleton anything that doesnā€™t define length(), but I suspect it would affect performance (reflection like that is probably too slow).

11

I see. I feel like weā€™re getting offtopic, but this behaviour of Ref seems weird to me. Itā€™s a container like a tuple, why isnā€™t it preserved?

Edit: Actually, sorry, I see that thereā€™s a difference, but I donā€™t see why it makes Ref a better choice.

12 
julia> axes((1,))
(Base.OneTo(1),)

julia> axes(Ref(1))
()

Basically, a tuple is one dimensional while a Ref is zero dimensional.

1 Like
13

As mentioned by others, broadcast defaults to treating its arguments as collections. Types can opt in to being treated as ā€œscalarsā€ by defining Base.broadcastable(x::MyType) = Ref(x), or you can opt-in at the call site with an explicit Ref.

That being said, I agree that RNG types should opt-in to being treated as scalars. Fortunately, it is a one-line non-breaking change to add this to any existing type:

https://github.com/JuliaLang/julia/issues/31071

3 Likes
14

Thanks for that, it seems like a better way to customise broadcasting behaviour than overriding length and iterate.

Iā€™m still a bit weirded out by using Ref in this context though. Namely, I donā€™t understand why it having zero dimensions makes it a better choice; and, more importantly, why Ref incidentally having zero dimensions (which is kind of implied by Ref objects are dereferenced (loaded or stored) with []. in its docstring, but seems like an implementation detail, because it can easily get a custom getindex() in future) means that it should be used here. Would it really be clear for anyone with some Julia programming experience why Ref is used in that line?

15

I think that historically Ref precedes the current broadcasting implementation, and was kind of hijacked for the purpose. It is innocuous, but I agree that it can be confusing. See

16

Itā€™s still an off topic, but let me mention that a syntax &x for making x a scalar in broadcast is proposed:

https://github.com/JuliaLang/julia/issues/27563

If this happens I think Ref just becomes a lower-level detail for most of users.

2 Likes
17

Honestly, I donā€™t really see a qualitative difference between it and Ref. Itā€™s a contraction, which is nice, but & is just as low-level. It is still associated with ā€œreferenceā€, not with ā€œscalarā€. My point (which I think people keep missing) was that I, or a user of my library, shouldnā€™t have to do anything to mark that something that is not an iterable (and the compiler knows that) should be treated as a scalar. @stevengj above agreed that RNG objects should be treated as scalars; why doesnā€™t it apply to any other non-iterable? Is there some conceptual problem with that, or is it just impossible to implement at this point without breaking existing code?

18

Itā€™s not currently possible to identify ā€œnon-iterableā€ reliably at compile-time, because ā€œiterableā€ is not a type.

There was a basic design choice that needed to be made for broadcast: should it default to treating arguments as scalars or as collections, with the opposite behavior being opt-in? Initially, it defaulted to treating things as scalars, but it was changed to default to collections after a long discussion (https://github.com/JuliaLang/julia/issues/18618 and https://github.com/JuliaLang/julia/pull/26435). At this point, the default is set in stone ā€” it cannot be changed in 1.x without breaking backwards compatibility.

There were reasonable arguments on both sides. One advantage of defaulting to a collection is that when a type like AbstractRNG arises where this is the wrong choice, we can update it in 1.x without breaking anything, which is not true of the opposite default.

2 Likes