Preventing collect on broadcasting expressions

General Usage
, ,
1

In writing efficient code for image processing, views are incredibly useful. They can nicely be chained for example to implement a phasor in optics propagation:

using IndexFunArrays
k0 = 5; phasor = exp.(1im .* 2pi .*sqrt.(k0^2 .- rr2((7,7))))
7×7 Matrix{ComplexF64}:

But the the phasor is again just a nother builden block in further expressions and it might be nice to intentionally keep it from being evaluated at definition. In other words I would like the broadcast / loop fusion operation to not automatically evaluate but rather encapsulate it in a view. I tried view(myexpression)
Yet this yields ERROR: BoundsError: attempt to access 7×7 Matrix{ComplexF64} at index []. Is there any known convenient way to prevent loop evaluation?
One may think that view() might do the job. But there are two problems: 1. view cannot be called with only a single array element, but always seems to require ranges which means at least the number of dimension needs be known. 2. The argument of view is also fused and evaluated if it contains such a fusion operation:

@time v = view(rr2((3000,3000)) .+ 0,:,:);
  0.046590 seconds (9 allocations: 68.665 MiB, 43.28% gc time)
 @time v = view(rr2((3000,3000)),:,:);
  0.000023 seconds (7 allocations: 336 bytes)

Any ideas how to prevent loop fusion upon assignment?

2

See eg BroadcastArray in

3

Thanks for the quick answer. I might be missing something here, but in my hands this never worked:

@time q = BroadcastArray(+, ones((3000,3000)), 0.0);
  0.039350 seconds (3 allocations: 68.665 MiB, 46.17% gc time)

It works if there are allocate operations, but the broadcasting seems incabable of dealing with on-the-fly generators like ones or rr2.

4

BroadcastArray still evaluates its arguments, so it can’t remove the allocation of ones. For that you could in theory use FillArrays.jl, but there appears to be some awkward interaction between how the two packages alter broadcasting:

julia> using LazyArrays, FillArrays

julia> o3k = ones((3000,3000));

julia> @time q = BroadcastArray(@~ o3k .+ 1.0);
  0.000028 seconds (6 allocations: 208 bytes)

julia> Fill(1.0, 3000, 3000) .+ 1.0
3000×3000 Fill{Float64, 2, Tuple{OneTo{Int64}, OneTo{Int64}}} = 2.0

julia> @time q = BroadcastArray(+, Fill(1.0, 3000, 3000), 1.0);
ERROR: type Fill has no field f
5

Thanks. But the whole point of this thread was to be able to chain operations that are capable to generate data on the fly (like rr2() or also ones()) together across multiple lines. Not as a default behaviour but on purpose. FillArrays seem to be a special case, as they essentially host only a single value and broadcasting needs to be applied to this one value only. So this is probably the way they achieve it there, but this breaks down quickly:

a = Fill(1.0,3000,3000) .* ones(3000,3000)
3000×3000 Matrix{Float64}:

Theoretically it is still one value but the mechanism is naturally incapable of seeing this.

6

@RainerHeitzmann, this discussion is well above my level but luck is a factor helping beginners, sometimes.

Does this function composition make sense?

using IndexFunArrays
k0 = 5; 
phasor(r) = exp.(1im .* 2pi .*sqrt.(k0^2 .- r))

n = 7; m = 7;
ph = IndexFunArray(phasor ∘ rr2,(n,m))

7×7 IndexFunArray{Matrix{ComplexF64}, 2, ComposedFunction{typeof(phasor), typeof(rr2)}}:
7

I am not sure what “mechanism” you are expecting here. In any case,

julia> Fill(1.0,3000,3000) .* Fill(1.0, 3000,3000)
3000×3000 Fill{Float64}: entries equal to 1.0

works. Once the bug mentioned by @mcabbott is fixed, it should work fine for your purposes.