8

I have a few large arrays (100 times larger than in the example) and I need to re-arrange them. Basically using tests on b and c, I need to pick elements of a.

a = RandomReal[1, 10^7];
b = RandomReal[1, 10^7];
c = RandomReal[1, 10^7];

sel = MapThread[(Abs[#1] <= 0.1) ∧ (Abs[#2] > 0.9)&, {b, c}]; // AbsoluteTiming

{10.9262,Null}

d = Pick[a, sel]; // AbsoluteTiming

{0.584831,Null}

MapThread and Pick are the only ways to do this, as far as I found, I can't re-arrange the code any more. The selection generation takes awhile, many times longer than on other platforms, is there a faster way to speed this up?

Alexey Popkov
  • 61,809
  • 7
  • 149
  • 368
Anatoly
  • 409
  • 2
  • 8

2 Answers2

14

You can use UnitStep to construct your selector array:

SeedRandom[1]
a = RandomReal[1, 10^7];
b = RandomReal[1, 10^7];
c = RandomReal[1, 10^7];

sel = MapThread[(Abs[#1] <= 0.1) ∧ (Abs[#2] > 0.9) &, {b, c}]; //
   AbsoluteTiming // First

15.332532

d = Pick[a, sel]; // AbsoluteTiming // First

0.605384

sel2 = UnitStep[.1 - Abs[b]](1- UnitStep[.9 - Abs[c]]); // AbsoluteTiming // First

0.630973

d2 = Pick[a, sel2, 1]; // AbsoluteTiming // First

0.015618

d == d2

True

If sel is already created, you can convert it a packed array. Additional processing time to do this is more than compensated for by speed gains in Pick:

sel3 = Developer`ToPackedArray[With[{True = 1, False = 0}, Evaluate[sel]]]; // 
  AbsoluteTiming // First

0.37475

d3 = Pick[a, sel3, 1]; // AbsoluteTiming // First

0.031245

d == d3

True

See this answer by Mr.Wizard for the origin of the With[{True = 1, False = 0}, ...] trick.

Update: Slightly faster way to construct the selector array (thanks: Fred Simons):

With[{sel3 = UnitStep[.1 - Abs[b]] - UnitStep[0.9 - Abs[c]]}, 
 Pick[a, sel3, 1]]
kglr
  • 394,356
  • 18
  • 477
  • 896
  • That's exactly what I was going to post. Always use vectorized functions if you can and always stay with machine integers/reals if you can. Furthermore, sel is not a packed array because MapThread doesn't pack automatically for some reason. That makes MapThread even less appealing. – Sjoerd Smit Oct 26 '18 at 20:38
  • One issue: you are effectively using Abs[c] >= 0.9 and not Abs[c] > 0.9 as OP requested. – Szabolcs Oct 26 '18 at 21:22
  • @Sjoerd If you use this technique, check out BoolEval (see my answer). – Szabolcs Oct 26 '18 at 21:22
  • @Szabolcs, thank you. I will update with a fix in a moment. – kglr Oct 26 '18 at 21:23
  • Replacing True/False with 1/0 is more easily done with Boole. It works on lists too. – Szabolcs Oct 26 '18 at 21:27
  • @Szabolcs, yes but Mr.Wizard's method is much faster than Boole. – kglr Oct 26 '18 at 21:29
  • Thank you @kglr, this works ! Still not sure which functions are vectorized,and exactly how packing works, but this makes sense – Anatoly Oct 26 '18 at 21:31
  • A slight simplification of the UnitStep construction, speeding it up: With[{sel=UnitStep[.1-Abs[b]]+UnitStep[Abs[c]-0.9]}, Pick[a, sel,2]]. – Fred Simons Nov 06 '18 at 15:55
  • Thank you @FredSimons. I updated with variant (to have Abs[c] > 0.9, not Abs[c] >= 0.9 ) of your suggestion. – kglr Nov 06 '18 at 21:46
10

@kglr already showed how to do this with vectorization, which offers much better performance.

However, the syntax with all those UnitSteps is not very friendly. My BoolEval package offers much better syntax for the very same operation.

Demo:

a = RandomReal[1, 10^7];
b = RandomReal[1, 10^7];
c = RandomReal[1, 10^7];

sel = MapThread[(Abs[#1] <= 0.1) ∧ (Abs[#2] > 0.9) &, {b, c}]; // AbsoluteTiming
(* {10.2691, Null} *)

d = Pick[a, sel]; // AbsoluteTiming
(* {0.460785, Null} *)

<< BoolEval`

e = BoolPick[a, Abs[b] <= 0.1 && Abs[c] > 0.9]; // AbsoluteTiming
(* {0.186719, Null} *)

e == d
(* True *)

We went from ~11 seconds to only 0.18 s, while shortening the code to a very readable one-liner. I hope this already convinced you that you really need BoolEval ;-)

I should also point out that @kglr's solution actually uses Abs[c] >= 0.9 and not Abs[c] > 0.9 as you requested. To achieve a true >, we could use 1 - UnitStep[0.9 - Abs[c]], which would make the code even less readable. BoolEval takes care of this for you.

One more thing to note:

Pick is much faster if the second argument is a packed array and the third argument is a simple number (not the default True).

Szabolcs
  • 234,956
  • 30
  • 623
  • 1,263
  • Yes, matlab notation is what I am used to, this will help keep my mind straight when I go between the two, thank you!! – Anatoly Oct 26 '18 at 21:32