Generating a list with elements based conditionals on other lists in Mathematica

Question

I want to calculate an 2D array "tocalc" in which the elements are calculated based on tests on three other lists (z,b1,b2).

(*example data*)
z = Range[0, 100, 10];
x = Range[10];
b1 = ConstantArray[0., Length[x]];
tocalc = ConstantArray[0, {Length[x], Length[z]}];
b2 = {0, 20, 30, 40, 50, 40, 30, 20, 10, 0};

one solution to this would be

(*simple but slow solution*)
Do[
 Do[
   If[z[[k]] <= b2[[i]] && z[[k]] >= b1[[i]], 
    tocalc[[i, k]] = (b2[[i + 1]] - b2[[i - 1]])],
   {k, 1, Length[z]}];,
 {i, 2, Length[x] - 1}]

with the result

{{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {30, 30, 30, 0, 0, 0, 0, 0, 0, 0, 
  0}, {20, 20, 20, 20, 0, 0, 0, 0, 0, 0, 0}, {20, 20, 20, 20, 20, 0, 
  0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
  0}, {-20, -20, -20, -20, -20, 0, 0, 0, 0, 0, 
  0}, {-20, -20, -20, -20, 0, 0, 0, 0, 0, 0, 0}, {-20, -20, -20, 0, 0,
   0, 0, 0, 0, 0, 0}, {-20, -20, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0,
   0, 0, 0, 0, 0, 0, 0, 0}}

The Question: How can this be done efficiently in Mathematica?

If this is evaluated 10000 times then it would take then 3.66 seconds. While in Matlab this takes 0.04sec so Matlab is almost 100 times faster. I know that the solution with the two Do loops is not perfect for Mathematica, so i tried several other solutions such as with MapIndexed, Table, functions, Conditionals, and so on. but all are not really faster maybe even slower the the two Do loops. Here is one example with MapIndexed:

tocalc = ConstantArray[0, {Length[x], Length[z]}];
MapIndexed[
  If[z[[Part[#2, 2]]] <= b2[[Part[#2, 1]]] && 
     z[[Part[#2, 2]]] >= b1[[Part[#2, 1]]] && Part[#2, 1] >= 2 && 
     Part[#2, 1] <= Length[x] - 1, 
    tocalc[[Part[#2, 1], Part[#2, 2]]] = (b2[[Part[#2, 1] + 1]] - 
       b2[[Part[#2, 1] - 1]]), 0.] &, tocalc, {2}];

The ideal solution should work for larger matrices and real numbers as well and also for more complicated conditionals.

---edit:

Since it looks some solutions to this are even slower in my real problem here is one example of it:

the Real-World Problem

b2 = {0.`, 0.`, 0.`, 990.3440201085594`, 1525.7589030785484`, 
   1897.6531659202747`, 2191.6073263357594`, 2433.0441988616717`, 
   2630.6658409463894`, 2799.347578394955`, 2944.656306810331`, 
   3070.718467691769`, 3179.485627984329`, 3272.3788096129415`, 
   3346.199103579602`, 3405.384848015466`, 3346.199103579602`, 
   3272.3788096129415`, 3179.485627984329`, 3070.718467691769`, 
   2944.656306810331`, 2799.347578394955`, 2630.6658409463894`, 
   2433.0441988616717`, 2191.6073263357594`, 1897.6531659202747`, 
   1525.7589030785484`, 990.3440201085594`, 0.`, 0.`, 0.`};
z = {0.`, 250.`, 500.`, 750.`, 1000.`, 1250.`, 1500.`, 1750.`, 2000.`,
   2250.`, 2500.`, 2750.`, 3000.`, 3250.`, 
  3500.`}; (*z(k)*)
imax = 31; (*number of x(i)*)
b1 = ConstantArray[0., imax]; (*lower boundary, can be different form 0*)
deltax = 50000.`;
mmax = 10000.; (*number of calculations*)
A00 = 1.127190283243198`*^-12; (*somefactor*)
n = 3;

one solution:

f2C = Compile[{{b2, _Real, 1}, {z, _Real, 1}, {b1, _Real, 1}},
   With[{zeros = {ConstantArray[0., Length[z]]}},
    Join[zeros, 
     Table[If[
       b1[[i]] <= z[[k]] <= 
        b2[[i]], -(A00*(Abs[(b2[[i + 1]] - b2[[i - 1]])/(2.0*
                 deltax)])^(n - 
              1.0)*(b2[[i]]^(n + 1.) - (b2[[i]] - z[[k]])^(n + 
                1.)))*((b2[[i + 1]] - b2[[i - 1]])/(2.0*deltax))
       , 0.],
      {i, 2, Length[b2] - 1}, {k, Length[z]}
      ], zeros]]
   , CompilationTarget -> "C"];

The Result is

Timing[Do[f2C[b2, z, b1];, {mmax}]]
Out[85]= {81.3544, Null}

Thanks!

When you say "If this is evaluated 10000 times then it would take more then 30 seconds.", do you mean wrapping your first two code blocks in Do[ ... , {10000}] // Timing? This runs in 7 seconds on a 6.5 year old mobile CPU here. Or maybe you mean 100000 times? Just trying to make sure I was benchmarking the correct thing. — Szabolcs
ok. thanks you are right I actually used 100000 times. I edited the text. — gogoolplex
How long does it take on MATLAB? Just so I know when to stop trying to optimize it :) "Vectorizing" it pushed the timing from 7 seconds to between 0.5 and 1 (i.e. around 10x speedup), without the need to compile. But the code was not pretty and unintuitive. For larger matrices the speedup would be greater. — Szabolcs
I now made the same example in Matlab and 10^4 evaluations. It took 0.04 sec and in Mathematica 3.66. So Matlab is about 91.5 times faster. (On a different machine but the Mathematica one is even newer). — gogoolplex

Daniel Lichtblau Daniel Lichtblau · Accepted Answer · 2012-05-30T13:45:31

You can do something like below. You will need to figure out how you want to handle the boundaries though (where b2[[i+1]] or b2[[i-1]] is not defined).

f[x_, y_] := If[x[[1]] <= y <= x[[2]], x[[4]] - x[[3]], 0]

Here I restrain the level to which Outer goes, so that I do not need to change the head (as I was doing in the original response).

In[1309]:= Outer[f, 
 Transpose[{b1, b2, RotateRight[b2], RotateLeft[b2]}], z, 1]

Out[1309]= {{20, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {30, 30, 30, 0, 0, 0, 
  0, 0, 0, 0, 0}, {20, 20, 20, 20, 0, 0, 0, 0, 0, 0, 0}, {20, 20, 20, 
  20, 20, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
  0}, {-20, -20, -20, -20, -20, 0, 0, 0, 0, 0, 
  0}, {-20, -20, -20, -20, 0, 0, 0, 0, 0, 0, 0}, {-20, -20, -20, 0, 0,
   0, 0, 0, 0, 0, 0}, {-20, -20, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {-10, 0, 
  0, 0, 0, 0, 0, 0, 0, 0, 0}}

Speed check:

In[1298]:= Timing[
 Do[Outer[f, 
   Apply[list, 
    Transpose[{b1, b2, RotateRight[b2], RotateLeft[b2]}], {1}], 
   z], {10^4}]]

Out[1298]= {2.68, Null}

We can compile the function to get better speed.

fC = Compile[{{x, _Integer, 1}, {y, _Integer}}, 
   If[x[[1]] <= y <= x[[2]], x[[4]] - x[[3]], 0]];

In[1306]:= Timing[
 Do[Outer[fC, Transpose[{b1, b2, RotateRight[b2], RotateLeft[b2]}], z,
    1], {10^4}]]

Out[1306]= {0.8, Null}

--- edit ---

Variants include compiling the entire routine. Here is one such.

ff = Compile[{{b1, _Integer, 1}, {b2, _Integer, 1}, {z, _Integer, 
     1}},
   With[{lc = 
      RotateRight[ListConvolve[{1, 0, -1}, b2, {-1, -1}, 0]]},
    Table[
     If[b1[[i]] <= z[[k]] <= b2[[i]], lc[[i]], 0], {i, 
      Length[b2]}, {k, Length[z]}
     ]]];
In[385]:= Timing[Do[ff[b1, b2, z], {10^4}]]

Out[385]= {0.24, Null}

If I add CompilationTarget -> "C" then it gets around twice as fast.

Another variant, in C code, gets under 0.1 seconds.

In[441]:= 
ff2C = Compile[{{b1, _Integer, 1}, {b2, _Integer, 1}, {z, _Integer, 
     1}},
   With[{zeros = {ConstantArray[0, Length[z]]}},
    Join[zeros, Table[
      If[b1[[i]] <= z[[k]] <= b2[[i]], b2[[i + 1]] - b2[[i - 1]], 
       0], {i, 2, Length[b2] - 1}, {k, Length[z]}
      ], zeros]], CompilationTarget -> "C"];

In[442]:= Timing[Do[ff2C[b1, b2, z], {10^4}]]

Out[442]= {0.04, Null}

In[443]:= ff2C[b1, b2, z]

Out[443]= {{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {30, 30, 30, 0, 0, 0, 0,
   0, 0, 0, 0}, {20, 20, 20, 20, 0, 0, 0, 0, 0, 0, 0}, {20, 20, 20, 
  20, 20, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
  0}, {-20, -20, -20, -20, -20, 0, 0, 0, 0, 0, 
  0}, {-20, -20, -20, -20, 0, 0, 0, 0, 0, 0, 0}, {-20, -20, -20, 0, 0,
   0, 0, 0, 0, 0, 0}, {-20, -20, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0,
   0, 0, 0, 0, 0, 0, 0, 0}}

I would guess there are variants that may be faster still.

--- end edit ---

--- edit 2 ---

Of course, if you have global variables (that is, defined outside of your Compile), then there is a bit more work to do. I am aware of two possibilities. Prior to version 8 one would suck in constants using a With[] around the Compile, as below.

f2C = With[{n = n, deltax = deltax, A00 = A00}, Compile[{{b2, _Real, 1}, {z, _Real, 1}, {b1, _Real, 1}}, With[{zeros = {ConstantArray[0., Length[z]]}}, Join[zeros, Table[If[ b1[[i]] <= z[[k]] <= b2[[i]], -(A00*(Abs[(b2[[i + 1]] - b2[[i - 1]])/(2.0* deltax)])^(n - 1.0)(b2[[i]]^(n + 1.) - (b2[[i]] - z[[k]])^(n + 1.)))((b2[[i + 1]] - b2[[i - 1]])/(2.0*deltax)), 0.], {i, 2, Length[b2] - 1}, {k, Length[z]}], zeros]], CompilationTarget -> "C"]];

In version 8 the following achieves the same effect.

f2Cb = Compile[{{b2, _Real, 1}, {z, _Real, 1}, {b1, _Real, 1}}, With[{zeros = {ConstantArray[0., Length[z]]}}, Join[zeros, Table[If[ b1[[i]] <= z[[k]] <= b2[[i]], -(A00*(Abs[(b2[[i + 1]] - b2[[i - 1]])/(2.0* deltax)])^(n - 1.0)(b2[[i]]^(n + 1.) - (b2[[i]] - z[[k]])^(n + 1.)))((b2[[i + 1]] - b2[[i - 1]])/(2.0*deltax)), 0.], {i, 2, Length[b2] - 1}, {k, Length[z]}], zeros]], CompilationTarget -> "C", CompilationOptions -> {"InlineExternalDefinitions" -> True}];

With either I get a result on the more realistic example in around 0.7 seconds, whereas my machine would take over 100 seconds without those globals being defined inside the Compile.

A more general approach might be to pass them as parameters (if they were likely to change rather than be constants). That would lead to a slightly slower run time though.

Regarding that option approach, you might have a look at ref/CompilationOptions in the Cocumentation Center

--- end edit 2 ---

Generating a list with elements based conditionals on other lists in Mathematica

1 Answers