How to mutate the value (not binding) of a scalar?

General Usage
1

I can do the following to mutate the original memory,

w = ones(3)

 3-element Vector{Float64}:
 1.0
 1.0
 1.0

 x = w

 3-element Vector{Float64}:
 1.0
 1.0
 1.0

 x[2]=2.

 2.0

but I can’t find a way to mutate w if it’s a scalar,

1

 x=w

 1

 x=2

 2

 w

 1

I think I need this because I want to broadcast over a function that takes a scratch or work array, and I need to modify its contents. But maybe there is a different idiom that I can use.

2

There is no way to mutate the value of a scalar such as 1 or 1.0 - they are immutable objects of type Int and Float64 respectively.

Can you expand on what you’re trying to do? There’s certainly a way achieve the same effect without trying to mutate immutables.

3

I’ve settled on this idiom, which is more elegant I think,

geman_mcclure(d2,s) = return one(s) / (one(s)+d2/s^2)^2
    
geman_mcclure!(d2,s,w) = @. w = geman_mcclure(d2,s)
4

I see - in julia, we’d typically have the mutated argument be the first argument of the function, so I’d write it like this:

geman_mcclure(d2, s) = one(s) / (one(s)+d2/s^2)^2
geman_mcclure!(w, d2, s) = geman_mcclure!(d2,s,w) = @. w = geman_mcclure(d2,s)

May I ask, what types do you expect to put into this? Is d2 an array?

5

d2 and w are Vector{T <: AbstractFloat}

6

Then I’d recommend implementing it like this:

function geman_mcclure!(w, d2, s)
    @. w = one(s) / (one(s) + d2 / s^2)^2
end
geman_mcclure(d2, s) = geman_cclure!(similar(d2), d2, s)

If d2 for the non-mutating version is also expected to be a Vector, as otherwise the method won’t work:

julia> gm1(d2, s) = one(s) / (one(s)+d2/s^2)^2
gm1 (generic function with 1 method)

julia> gm1(rand(Float64, 10), rand(Float64))
ERROR: MethodError: no method matching +(::Float64, ::Vector{Float64})
For element-wise addition, use broadcasting with dot syntax: scalar .+ array

Closest candidates are:
  +(::Any, ::Any, ::Any, ::Any...)
   @ Base operators.jl:578
  +(::T, ::T) where T<:Union{Float16, Float32, Float64}
   @ Base float.jl:406
  +(::Union{Float16, Float32, Float64}, ::BigFloat)
   @ Base mpfr.jl:414
  ...

Stacktrace:
 [1] gm1(d2::Vector{Float64}, s::Float64)
   @ Main ./REPL[7]:1
 [2] top-level scope
   @ REPL[8]:1

Though usually, we try to avoid writing an explicitly vectorized function like geman_mcclure! and just directly use broadcasting inline, to allow for loop fusion to occur in longer chains (in which case your original scalar version is perfect already, and no other methods would be needed). It’s just a question of what you’d expect users to pass in.

1 Like
7

Thanks for the tips, but I don’t understand your last comment. In this case, geman_mcclure will be called with vectors of data passed in d2,s. Will broadcasting on the function allow the arithmetic operations to be fused? And can you give an example of what you mean by “broadcasting inline?”

And the non_mutating versions can be scalars. I’m usually writing functions for scalars and broadcast over the function at the caller level. I try not to write vectorized code at the function level, but now I see how your suggestion can accommodate both types of calls. Thanks!

8

Generally, yes - writing geman_mcclure.(d2, s) with d2 a Vector and s a scalar (or even a vector of the same length as d2) is equivalent to

out_arr = similar(d2)
for i in eachindex(d2)
    out_arr[i] = geman_mcclure(d2[i], s) # geman_mcclure(d2[i], s[i]) if `s` is broadcastable
end
out_arr

which is the same as your explicitly vectorized geman_mcclure! function, modulo the additional allocation, which can be done ahead of time manually. If done ahead of time, you can (just with geman_mcclure) then write w .= geman_mcclure.(d2, s) , which is already equivalent to the loop above, except that out_arr does not have to be allocated, because it already exists.

This is a common pattern - instead of writing explicitly vectorized versions of scalar functions like geman_mcclure!, we usually only have the scalar version and broadcast/inplace-broadcast/map as necessary, allowing loop fusion to take place.

In comparison, geman_mcclure! does not allow loop fusion with broadcast outside of itself. That is, this example will not fuse loops:

z .= a .+ geman_mcclure!(w, ds2, s)

(there will first be the inner loop inside of geman_mcclure!, then the outer one over z and a) whereas this will:

z .= a .+ geman_mcclure.(ds2, s)

which is in fact equivalent to

for idx in eachindex(z, a, ds2)
    z[idx] = a[idx] + geman_mcclure(ds2[idx], s)
end

This is also what I mean with “broadcasting inline” - simply using the scalar version in a broadcast is more flexible than having geman_mcclure!, which is required to complete before the outer broadcast can be calculated. The general rule of thumb for when I use explicitly modifying functions is either when the operation is not a scalar operation (i.e., more than a single index is required to perform the calculation, forcing you to write a loop either way) or when there’s a truly persistent modification of a datastructure as a sideeffect happens, e.g. with get!:

help?> get!
search: get! mergewith! get getpid getkey getfield getindex getglobal

  get!(collection, key, default)

  Return the value stored for the given key, or if no mapping for the
  key is present, store key => default, and return default.

  Examples
  ≡≡≡≡≡≡≡≡≡≡

  julia> d = Dict("a"=>1, "b"=>2, "c"=>3);
  
  julia> get!(d, "a", 5)
  1
  
  julia> get!(d, "d", 4)
  4
  
  julia> d
  Dict{String, Int64} with 4 entries:
    "c" => 3
    "b" => 2
    "a" => 1
    "d" => 4
1 Like
9

Thanks for the helpful discussion!