OffsetArrays do not work in let blocks as expected

General Usage
, ,
1

Hi all,

I am using OffsetArrays in my work and recently ran into the issue when using let blocks. MWE without OffsetArrays


let
b = zeros(20,20)
a = ones(20)
global test
function test(p)
  b .= a .* p
end
end
o = ones(20)
@btime for i in 1:1000
  test($o)
end

This runs as expected, and provides 0 allocations.
Now for my actual use case. I have the following code:

const ⊕ = Base.FastMath.add_fast # fast versions for fma
const ⊖ = Base.FastMath.sub_fast
const ⊗ = Base.FastMath.mul_fast
const ⨸ = Base.FastMath.div_fast
const psi = 1.0
const g = 16.0
const k = sqrt(g*psi)
const p = (psi,g,k)
const N = 5

let
dummy = zeros(ComplexF64, 6N+1,2N+1,2N+1)
odummy = OffsetArray(dummy,-3N:3N,-N:N,-N:N)

global C_mat
function C_mat(u::Vector{Complex{Float64}},p::NTuple{3,Float64},t::Float64)

u = OffsetArray(u,-N:N) #if testing, note that dimensionality of u is 2N+1
 
 for q = -3N:1:3N
   for n = -N:1:N
     for m = -N:1:N
       if abs(q-n+m) <= abs(N)
         @inbounds odummy[q,n,m] = (p[2]⨸2)⊗u[n]⊗conj(u[m])⊗u[q-n+m] ::ComplexF64
       end
     end
   end
end

odummy = sum(odummy,dims=2)
final = dropdims(odummy,dims=2)
return parent(final[-N:N,:])
end
end

Now if you call this once, it runs fine. However, if I put this in the loop to evaluate multiple times, I get the stacktrace:

ERROR: BoundsError: attempt to access 31×1×11 OffsetArray(::Array{Complex{Float64},3}, -15:15, -5:-5, -5:5) with eltype Complex{Float64} with indices -15:15×-5:-5×-5:5 at index [-14, -4, 5]
Stacktrace:
 [1] throw_boundserror(::OffsetArray{Complex{Float64},3,Array{Complex{Float64},3}}, ::Tuple{Int64,Int64,Int64}) at ./abstractarray.jl:537
 [2] checkbounds at ./abstractarray.jl:502 [inlined]
 [3] setindex!(::OffsetArray{Complex{Float64},3,Array{Complex{Float64},3}}, ::Complex{Float64}, ::Int64, ::Int64, ::Int64) at /Users/nvv/.julia/packages/OffsetArrays/JNUm0/src/OffsetArrays.jl:157
 [4] macro expansion at /Users/nvv/Desktop/julia/test_mat.jl:26 [inlined]
 [5] macro expansion at ./simdloop.jl:77 [inlined]
 [6] macro expansion at /Users/nvv/Desktop/julia/test_mat.jl:24 [inlined]
 [7] macro expansion at ./simdloop.jl:77 [inlined]
 [8] macro expansion at /Users/nvv/Desktop/julia/test_mat.jl:23 [inlined]
 [9] macro expansion at ./simdloop.jl:77 [inlined]
 [10] C_mat(::Array{Complex{Float64},1}, ::Tuple{Float64,Float64,Float64}, ::Int64) at /Users/nvv/Desktop/julia/test_mat.jl:22
 [11] macro expansion at ./REPL[30]:2 [inlined]
 [12] ##core#396(::Array{Complex{Float64},1}, ::Tuple{Float64,Float64,Float64}, ::Int64) at /Users/nvv/.julia/packages/BenchmarkTools/eCEpo/src/execution.jl:371
 [13] ##sample#397(::BenchmarkTools.Parameters) at /Users/nvv/.julia/packages/BenchmarkTools/eCEpo/src/execution.jl:377
 [14] _run(::BenchmarkTools.Benchmark{Symbol("##benchmark#395")}, ::BenchmarkTools.Parameters; verbose::Bool, pad::String, kwargs::Base.Iterators.Pairs{Symbol,Integer,NTuple{4,Symbol},NamedTuple{(:samples, :evals, :gctrial, :gcsample),Tuple{Int64,Int64,Bool,Bool}}}) at /Users/nvv/.julia/packages/BenchmarkTools/eCEpo/src/execution.jl:405
 [15] (::Base.var"#inner#2"{Base.Iterators.Pairs{Symbol,Integer,NTuple{5,Symbol},NamedTuple{(:verbose, :samples, :evals, :gctrial, :gcsample),Tuple{Bool,Int64,Int64,Bool,Bool}}},typeof(BenchmarkTools._run),Tuple{BenchmarkTools.Benchmark{Symbol("##benchmark#395")},BenchmarkTools.Parameters}})() at ./essentials.jl:715
 [16] #invokelatest#1 at ./essentials.jl:716 [inlined]
 [17] #run_result#37 at /Users/nvv/.julia/packages/BenchmarkTools/eCEpo/src/execution.jl:32 [inlined]
 [18] run(::BenchmarkTools.Benchmark{Symbol("##benchmark#395")}, ::BenchmarkTools.Parameters; progressid::Nothing, nleaves::Float64, ndone::Float64, kwargs::Base.Iterators.Pairs{Symbol,Integer,NTuple{5,Symbol},NamedTuple{(:verbose, :samples, :evals, :gctrial, :gcsample),Tuple{Bool,Int64,Int64,Bool,Bool}}}) at /Users/nvv/.julia/packages/BenchmarkTools/eCEpo/src/execution.jl:94
 [19] #warmup#45 at /Users/nvv/.julia/packages/BenchmarkTools/eCEpo/src/execution.jl:141 [inlined]
 [20] warmup(::BenchmarkTools.Benchmark{Symbol("##benchmark#395")}) at /Users/nvv/.julia/packages/BenchmarkTools/eCEpo/src/execution.jl:141
 [21] top-level scope at /Users/nvv/.julia/packages/BenchmarkTools/eCEpo/src/execution.jl:481

which is not what I expect. I use let block to allocate array statically and avoid array initialization on every call - this is crucial to my implementation. Is there a way of bypassing the error whilst achieving the same goal, or should I start thinking about re-indexing?

Thank you!

2

Here’s a simpler demonstration:

julia> f = let
         z = rand(3)
         function g(x)
           x .= z
           z = [1]
           sum(x)
         end
       end;

julia> f(ones(3))
1.5738779450861395

julia> f(ones(3))
3.0

The variable z is not local to the function, and thus binding it to a new value changes the result of the second run.

If you don’t want this, then you should give the sum a different name, such as odummy2 = sum(odummy,dims=2).

For completeness, the minimal way to get the error from the above is to run this:

C_mat(rand(2N+1).+im, (1.1, 2.2, 3.3), 4.4) # axes(odummy) = (-15:15, -5:5, -5:5)
C_mat(rand(2N+1).+im, (1.1, 2.2, 3.3), 4.4) # axes(odummy) = (-15:15, -5:-5, -5:5)
1 Like
3

Thank you, this worked just fine!