Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I'd like to select (offset) indexes of a number of array points, that match a selection criteria. I can do this by creating a Table of all the coordinates, and then selecting all the values:

data = Array[ RandomInteger[1] &, {4, 4}]
{wx, wy} = data // Dimensions ;
t = Flatten[
  Table[ {x - wx/2, y - wy/2, data[[x]][[y]]}, {x, 1, wx}, {y, 1,
    wy}], 1]
Select[ t, #[[3]] == 1 &]

This resulted in the respective values:

{{1, 1, 1, 1}, {1, 0, 0, 1}, {1, 0, 1, 0}, {1, 0, 0, 1}}
{{-1, -1, 1}, {-1, 0, 1}, {-1, 1, 1}, {-1, 2, 1}, {0, -1, 1}, {0, 0,
  0}, {0, 1, 0}, {0, 2, 1}, {1, -1, 1}, {1, 0, 0}, {1, 1, 1}, {1, 2,
  0}, {2, -1, 1}, {2, 0, 0}, {2, 1, 0}, {2, 2, 1}}
{{-1, -1, 1}, {-1, 0, 1}, {-1, 1, 1}, {-1, 2, 1}, {0, -1, 1}, {0, 2,
  1}, {1, -1, 1}, {1, 1, 1}, {2, -1, 1}, {2, 2, 1}}

However, my real data is big (800x800) and sparse. This method creates a big data set and then throws most of it away. What is a more efficient way to do this?

share|improve this question
    
Is your data a SparseArray? Can you make it a SparseArray? –  Szabolcs May 28 at 14:27
    
BTW you can write RandomInteger[1, {4,4}] and data[[x,y]]. –  Szabolcs May 28 at 14:37

4 Answers 4

up vote 2 down vote accepted

Please consider:

{wx, wy} = {4, 4};

SeedRandom[1];
data = RandomInteger[1, {wx, wy}]

{#1 - wx/2, #2 - wy/2, data[[##]]} & @@@ Position[data, 1]
{{-1, -1, 1}, {-1, 0, 1}, {-1, 2, 1}, {0, 2, 1}, {1, 0, 1}}

For better performance you can use SparseArray to directly compute the list of positions, then use the "NonzeroPositions" Property to extract them. For example:

SeedRandom[1]
m = RandomInteger[5, {10, 10}];

Optimized code to find the positions of all elements equal to four:

sa = SparseArray[Unitize @ Subtract[m, 4], Automatic, 1];

pos = sa["NonzeroPositions"]
{{1, 1}, {1, 3}, {2, 6}, {2, 7}, {3, 1}, {3, 3}, {3, 7},
 {3, 9}, {5, 5}, {6, 8}, {7, 4}, {8, 5}, {9, 6}, {10, 3}}

(The long form of Subtract is used because peculiarly it is faster.)

You could then operate on these positions as I did above. Better still, set all non-selection elements in the array to the same value, and use that as the background of the SparseArray, then also use the "NonzeroValues" Property as well. For example let's select all elements in the interval $[2, 4]$:

{wx, wy} = {14, 20};

SeedRandom[1]
m = RandomInteger[99, {wx, wy}];

sa = SparseArray[Clip[m, {2, 4}, {0, 0}]];

pos = sa["NonzeroPositions"]
val = sa["NonzeroValues"]
{{1, 5}, {1, 13}, {3, 5}, {3, 9}, {4, 14}, {9, 5}, {9, 20}, {11, 10}, {14, 16}}

{3, 4, 4, 2, 3, 2, 2, 2, 2}

We can then transform positions using fast vector operations:

pos2 = (pos\[Transpose] - {wx/2, wy/2})\[Transpose]  (* this looks much nicer in a Notebook *)
{{-6, -5}, {-6, 3}, {-4, -5}, {-4, -1}, {-3, 4}, {2, -5}, {2, 10}, {4, 0}, {7, 6}}

And combine with the values using Join, which does not unpack:

Join[pos2, val ~Partition~ 1, 2]

% // Developer`PackedArrayQ
{{-6, -5, 3}, {-6, 3, 4}, {-4, -5, 4}, {-4, -1, 2}, {-3, 4, 3},
 {2, -5, 2}, {2, 10, 2}, {4, 0, 2}, {7, 6, 2}}

True

Note that using SparseArray properties with packed/packable data is much faster than using ArrayRules:

{wx, wy} = {14000, 20000};

SeedRandom[1]
m = RandomInteger[99, {wx, wy}];

sa = SparseArray[Clip[m, {2, 4}, {0, 0}]];

Timing[
 pos = sa["NonzeroPositions"];
 val = sa["NonzeroValues"];
]

Timing[
 rules = Most @ ArrayRules[sa];
]
{0.0088, Null}

{4.649, Null}
share|improve this answer

Here's some sparse sample data of size 800 by 800:

am = AdjacencyMatrix@RandomGraph[{800, 1000}];
data = am.am;
data -= DiagonalMatrix@Diagonal[data];

It contains 1s and 2s (other than zeros).

Let's find the positions of 2s:

Cases[Most@ArrayRules[data], HoldPattern[position_ -> 2] :> position]

These positions can be easily offset after they're found.

You're probably looking something more like this variation:

Cases[Most@ArrayRules[data], HoldPattern[position_ -> value_ /; value == 2] :> Append[position, value]]

This approach is best when you have sparse arrays to work with.

share|improve this answer

It seems that, depending on how dense the input matrix is, Position may perform better than the Cases/ ArrayRules combination:

am = AdjacencyMatrix@RandomGraph[{800, 1000}];
data = am.am;
data -= DiagonalMatrix@Diagonal[data];

ClearSystemCache[];
out1 = Cases[Most@ArrayRules[data], HoldPattern[position_ -> 2] :> position]; // Timing
(* {0.015625,Null} *)
ClearSystemCache[];
out2 = Position[Normal@data, 2, {2}]; // Timing
(* {0.03125,Null} *)
Length@out1/(800  800) // N
(* 0.0000375 *)
Sort@out1 == Sort@out2
(* True *)

For a denser matrix

am = AdjacencyMatrix@RandomGraph[{800, 10000}];
data = am.am;
data -= DiagonalMatrix@Diagonal[data];
ClearSystemCache[];
out1 = Cases[Most@ArrayRules[data], HoldPattern[position_ -> 2] :> position]; // Timing
(* {0.3125,Null} *)
ClearSystemCache[];
out2 = Position[Normal@data, 2, {2}]; // Timing
(* {0.046875,Null} *)
Length@out1/(800  800) // N
(* 0.140584375` *)
Sort@out1 == Sort@out2 
(* True *)
share|improve this answer

Use functions Sow and Reap

t = Reap[
  Do[
   If[data[[x]][[y]] == 1, Sow[{x - wx/2, y - wy/2, data[[x]][[y]]}]],
   {x, 1, wx}, {y, 1, wy}
   ]
  ]
t = t[[2,1]]
share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.