List of indices of matrix satisfying certain criterion

[David P. Sanders]

Hi,

Given a matrix, which will be large (say 10^5 x 10^5), I need to extract the list of indices (i.e. the pairs (x,y) of positions) of those places in the matrix

where the value stored satisfies a certain condition. For a minimal example, the condition can just be that the value is greater than 0.5.

The code below achieves this, but seems inefficient, since it constructs the whole array of indices, even though there may be only a few places where the condition is satisfied.

Is there a more efficient way of doing this? 

Is there an "if" clause in an array comprehension? -- adding this to the definition of 'indices' would seem to do what I want, but I have not found this syntax in the manual for arrays.

Thanks,

David.

julia> L = 5

5

julia> r= rand(L, L)

5x5 Array{Float64,2}:

 0.705585   0.534721  0.158935  0.343876  0.624299

 0.0624089  0.525414  0.131139  0.590439  0.554686

 0.190085   0.557751  0.591916  0.485526  0.6307

 0.365398   0.943102  0.575083  0.858705  0.105142

 0.047924   0.116424  0.756757  0.576293  0.461547

julia> indices = [(x,y) for x in 1:L, y in 1:L]

5x5 Array{(Any,Any),2}:

 (1,1)  (1,2)  (1,3)  (1,4)  (1,5)

 (2,1)  (2,2)  (2,3)  (2,4)  (2,5)

 (3,1)  (3,2)  (3,3)  (3,4)  (3,5)

 (4,1)  (4,2)  (4,3)  (4,4)  (4,5)

 (5,1)  (5,2)  (5,3)  (5,4)  (5,5)

julia> where = indices[r .> 0.5]

14-element Array{(Any,Any),1}:

 (1,1)

 (1,2)

 (2,2)

 (3,2)

 (4,2)

 (3,3)

 (4,3)

 (5,3)

 (2,4)

 (4,4)

 (5,4)

 (1,5)

 (2,5)

 (3,5)

[Mauro]

> Given a matrix, which will be large (say 10^5 x 10^5), I need to extract
> the list of indices (i.e. the pairs (x,y) of positions) of those places in
> the matrix
> where the value stored satisfies a certain condition. For a minimal
> example, the condition can just be that the value is greater than 0.5.

Give you a logical array which you can use for indexing:
r.>.5
and to get the linear indices
find(r.>.5)

> The code below achieves this, but seems inefficient, since it constructs
> the whole array of indices, even though there may be only a few places
> where the condition is satisfied.
> Is there a more efficient way of doing this?
> Is there an "if" clause in an array comprehension? -- adding this to the
> definition of 'indices' would seem to do what I want, but I have not found
> this syntax in the manual for arrays.

There is no if clause for comprehensions but you can always just write
the loops.

--

[Ivar Nesje]

There is an open issue for adding a filter/if to a list comprehension. https://github.com/JuliaLang/julia/issues/550

for now you will need to do

L = 5

r= rand(L, L)

where = (Int,Int)[]

for x = 1:L, y= 1:L

    if r[x,y] < 0.5

        push!(where, (x,y))

    end

end

Preferably you will wrap that in a function, so that you do not use global variables, as they are currently slow.

Ivar

[David P. Sanders]

El viernes, 21 de febrero de 2014 08:25:00 UTC-6, Mauro escribió:

> Given a matrix, which will be large (say 10^5 x 10^5), I need to extract
> the list of indices (i.e. the pairs (x,y) of positions) of those places in
> the matrix
> where the value stored satisfies a certain condition. For a minimal
> example, the condition can just be that the value is greater than 0.5.

Give you a logical array which you can use for indexing:
  r.>.5
and to get the linear indices
  find(r.>.5)

Great, didn't know about find.

Maybe I should just rewrite everything thinking in terms of 1D arrays. (This is something I have taught in several courses, 

but for some reason it didn't occur to me!)

Thanks,

David.

[David P. Sanders]

El viernes, 21 de febrero de 2014 08:31:20 UTC-6, Ivar Nesje escribió:

There is an open issue for adding a filter/if to a list comprehension. https://github.com/JuliaLang/julia/issues/550

for now you will need to do

L = 5

r= rand(L, L)

where = (Int,Int)[]

for x = 1:L, y= 1:L

    if r[x,y] < 0.5

        push!(where, (x,y))

    end

end

Preferably you will wrap that in a function, so that you do not use global variables, as they are currently slow.

OK, thanks.

The version in a function can be profiled with `@time`.

How can I profile the version with globals (just to convince myself!)?

David.

[Mauro]

On Sun, 2014-02-23 at 05:51, dpsa...@gmail.com wrote:
> El viernes, 21 de febrero de 2014 08:25:00 UTC-6, Mauro escribió:
>>
>> > Given a matrix, which will be large (say 10^5 x 10^5), I need to extract
>> > the list of indices (i.e. the pairs (x,y) of positions) of those places
>> in
>> > the matrix
>> > where the value stored satisfies a certain condition. For a minimal
>> > example, the condition can just be that the value is greater than 0.5.
>>
>> Give you a logical array which you can use for indexing:
>> r.>.5
>> and to get the linear indices
>> find(r.>.5)
>>
>
> Great, didn't know about find.
> Maybe I should just rewrite everything thinking in terms of 1D arrays.
> (This is something I have taught in several courses,
> but for some reason it didn't occur to me!)

No need to rewrite in terms of 1D arrays. You can index both with a
linear index or with a 2D index (or higher D, depending on the array):

julia> a = [1 2; 3 4]
2x2 Array{Int64,2}:
1 2
3 4

julia> a[2]
3

julia> a[2,1]
3

If you want to translate between the two use ind2sub and sub2ind:

julia> ind = ind2sub(size(a), 2)
(2,1)

julia> a[ind...]
3

[Ivar Nesje]

@time begin
#code
end

[David P. Sanders]

Excellent, thanks! 

[David P. Sanders]

El domingo, 23 de febrero de 2014 15:56:35 UTC-6, Ivar Nesje escribió:

@time begin
   #code
end

Should have guessed that, thanks!