BoundsError when trying to index certain dimensions of an array with a vector containing indices

A = rand(3,3)
b = findall(vec(A).>0.5) # contains indices
c = rand(2,length(b)) # some values that we want to have at the corresponding indexes stored in b, in array B.
B = rand(2,3,3)
B[:,b] = c # this line returns an error

ERROR: BoundsError: attempt to access 2Γ—3Γ—3 Array{Float64, 3} at index [1:2, [1, 4, 5, 7, 8, 9]]
Stacktrace:
 [1] throw_boundserror(A::Array{Float64, 3}, I::Tuple{Base.Slice{Base.OneTo{Int64}}, Vector{Int64}})
   @ Base ./abstractarray.jl:651
 [2] checkbounds
   @ ./abstractarray.jl:616 [inlined]
 [3] _setindex!
   @ ./multidimensional.jl:886 [inlined]
 [4] setindex!(::Array{Float64, 3}, ::Matrix{Float64}, ::Function, ::Vector{Int64})
   @ Base ./abstractarray.jl:1267
 [5] top-level scope
   @ REPL[15]:1

In the last line, I would like that the value in c fill B in its 3x3 part (second and third dimension of B) at the locations contained in b for a 3x3 matrix (locations obtained from A). How should I write the last line ?

FYI, it would be much better to make your B 3x3x2.


Your A is a 3x3 mask (based on > 0.5 or not), your B has 2 β€œlayers” of 3x3 matrix.

My understanding is that you want to assign random values to each of those two layers at location according to mask A.

julia> for i = CartesianIndices(A)
           if A[i] > 0.5
               @. B[:, i] = rand()
           end
       end

this should do?

1 Like

The solution I was looking for is

A = rand(3,3)
b = findall(vec(A).>0.5) # contains indices
c = rand(2,length(b)) # some values that we want to have at the corresponding indexes stored in b, in array B.
B = rand(2,3,3)
b = CartesianIndices(A)[b] # ADDING THIS LINE SOLVES THE PROBLEM
B[:,b] = c # b as CartesianIndices is legal here

Vector b as CartesianIndices is a legal input in this scenario (instead of β€˜LinearIndices’).

@jling Thank you for your suggestion with the loop, but I will avoid using a loop if there is an alternative.

why, Julia is not MATLAB/Python/R, loops are fast in Julia and allocate less than your β€œvectorzied” code.

julia> A = rand(3,3);

julia> function f(A,B)
           b = findall(vec(A).>0.5) 
           c = rand(2,length(b))
           b = CartesianIndices(A)[b]
           B[:,b] = c
           B
       end

julia> @benchmark f($A, B) setup=(B=rand(2,3,3))
BenchmarkTools.Trial: 10000 samples with 376 evaluations.
 Range (min … max):  256.104 ns …   9.701 ΞΌs  β”Š GC (min … max): 0.00% … 93.84%
 Time  (median):     263.614 ns               β”Š GC (median):    0.00%
 Time  (mean Β± Οƒ):   302.790 ns Β± 381.676 ns  β”Š GC (mean Β± Οƒ):  8.49% Β±  6.50%

  β–…β–ˆβ–…β–‚β–                                                         ▁
  β–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–†β–†β–‡β–ˆβ–ˆβ–‡β–‡β–‡β–‡β–†β–…β–†β–‡β–‡β–‡β–‡β–†β–‡β–†β–‡β–ˆβ–‡β–‡β–†β–…β–…β–„β–„β–„β–„β–†β–„β–…β–…β–ƒβ–…β–„β–„β–„β–„β–…β–…β–…β–…β–„β–…β–†β–†β–†β–†β–…β–…β–„β–ƒβ–„ β–ˆ
  256 ns        Histogram: log(frequency) by time        554 ns <

 Memory estimate: 608 bytes, allocs estimate: 7.

julia> function g(A,B)
           for i = CartesianIndices(A)
                  if A[i] > 0.5
                      @. B[:, i] = rand()
                  end
              end
        B
       end
g (generic function with 1 method)

julia> @benchmark g($A, B) setup=(B=rand(2,3,3))
BenchmarkTools.Trial: 10000 samples with 988 evaluations.
 Range (min … max):  51.816 ns … 381.994 ns  β”Š GC (min … max): 0.00% … 0.00%
 Time  (median):     52.579 ns               β”Š GC (median):    0.00%
 Time  (mean Β± Οƒ):   53.107 ns Β±   4.857 ns  β”Š GC (mean Β± Οƒ):  0.00% Β± 0.00%

  β–…β–ƒβ–ˆβ–   ▁                                                     ▁
  β–ˆβ–ˆβ–ˆβ–ˆβ–‡β–‡β–†β–ˆβ–‡β–†β–…β–ƒβ–…β–…β–‡β–†β–…β–ƒβ–ƒβ–β–ƒβ–„β–β–„β–„β–ƒβ–β–β–β–ƒβ–β–„β–„β–β–β–ƒβ–ƒβ–…β–…β–„β–ƒβ–„β–…β–…β–…β–„β–†β–†β–„β–…β–†β–…β–„β–„β–†β–…β–…β–ƒβ–†β–„ β–ˆ
  51.8 ns       Histogram: log(frequency) by time      72.2 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

So my version is 4x faster from micro benchmark and has 0 allocation?

1 Like