Why does array packing in Table behave like this?

Question

Consider the following example, which generate two identical arrays. Why sometimes the array is packed and sometimes it doesn't? And why they perform differently?

Here the arrays are packed:

ht = 0.002; Ti = 0.4;
a1 = Table[Sin[x], {x, 0.2, 50., 0.02}]; // AbsoluteTiming
a2 = Table[Sin[x], {x, 0.5 Ti, 125 Ti, 0.05 Ti}]; // AbsoluteTiming
Developer`PackedArrayQ[a1]
Developer`PackedArrayQ[a2]
a1 == a2

(*
{0.010927, Null}
{0.335707, Null}
True
True
True
*)

Here one array is packed and the other is not:

a1 = Table[Sin[x], {x, 0.2, 50., 0.2}, {t, 0.002, 25., 0.002}]; // AbsoluteTiming
a2 = Table[Sin[x], {x, 0.5 Ti, 125 Ti, 0.5 Ti}, {t, ht, 12500 ht, ht}]; // AbsoluteTiming
Developer`PackedArrayQ[a1]
Developer`PackedArrayQ[a2]
a1 == a2

(*
{0.371184, Null}
{1.610365, Null}
True
False
True
*)

Since they have huge difference in performance, does that mean we should avoid using variables , instead we should use numbers in the code? Or add Evaluate everywhere we have variables?

Answer 1

I think the deeper reason is really contained in the links that Leonid has given, but probably not as to the point for your specific example as possible. The following cures your problem, without the need to avoid variables:

With[{Ti = Ti},
  a2 = Table[Sin[x], {x, 0.5 Ti, 125 Ti, 0.05 Ti}]; // AbsoluteTiming
]

Edit: as some comments from people more careful than me have pointed out the above doesn't cure the problem (or more precise in the one-iterator case there isn't a problem to be cured). For the case with two iterators solution along the lines of the above doesn't help, persumably because there are still products in the iterator. The following does work as intended (at least for me):

With[{xmin = 0.5 Ti, xmax = 125 Ti, xstep = 0.5 Ti, ht = ht, tmax = 12500 ht},
  a2 = Table[Sin[x], {x, xmin, xmax, xstep}, {t, ht, tmax, ht}]; // AbsoluteTiming
]

Of course this now is much more work for the general case than the methods in Mr. Wizards answer. It also has advantages, though: it does preserve the localization and might be somewhat easier to read and understand for not so advanced users.

What I do is to literally insert the numeric values so that Table now sees that it has to handle constant values and doesn't have to take into account potential complications which could arise from non constant iterators. You might find detailed explanations why that's necessary elsewhere, but sometimes a simple example is easier to understand: Try to predict the result if you use this definition of Ti and you may understand why Table needs to be extra careful whenever it sees symbols in its input:

Module[{x = 10}, Ti := (x++)];
Table[x, {x, 0, 2*Ti, Ti}]

Edit: it is interesting to note that my example is actually identical to the first example in the question, where auto-compilation/packaging did work nevertheless. It looks like for just one iterator the optimizations can be done even when a symbol appears in it, probably because then the generated list always can be assumed to be rectangular (1xn, that is)...

Answer 2

Notice the attributes of Table - HoldAll it means that the arguments (the function to be evaluated, AND the iterators) are not evaluated prior to executing the Table[...]. Therefore, the second loop involves computations in the iterators if you change the second Table to

Table[f[t, x], Evaluate[{x, 0.5 Ti, 125 Ti, 0.5 Ti}], Evaluate[{t, ht, 12500 ht, ht}]]; // AbsoluteTiming

it will enforce only a single evaluation of the values and the rest is similar to the iteration of the first loop, so (almost) identical results in timing

a1 = Table[
    f[t, x], {x, 0.2, 50., 0.2}, {t, 0.002, 25., 
     0.002}]; // AbsoluteTiming
a2 = Table[f[t, x], Evaluate[{x, 0.5 Ti, 125 Ti, 0.5 Ti}], 
Evaluate[{t, ht, 12500 ht, ht}]]; // AbsoluteTiming

results with

{0.673810, Null}

{0.671783, Null}

yehuda

Answer 3

As described in the Q&A's Leonid referenced⁽¹⁾⁽²⁾, Table needs to see explicit iterator values to trigger the internal optimizations. You can pre-evaluate the arguments of Table like this:

packedQ = Developer`PackedArrayQ;

ht = 0.002; Ti = 0.4;

Table @@ {Sin[x], {x, 0.5 Ti, 125 Ti, 0.5 Ti}, {t, ht, 12500 ht, ht}} // packedQ

True

Or if Sin[x] would evaluate incorrectly:

Table[Sin[x], ##] &[{x, 0.5 Ti, 125 Ti, 0.5 Ti}, {t, ht, 12500 ht, ht}] // packedQ

True

Be aware that neither method localizes the iterator variables as Table normally does. You could use Block to do it manually. Alternatively you could inject only the range parameters like this:

x = 5; t = 17; (* global values to demonstrate functioning protection *)

{{0.5 Ti, 125 Ti, 0.5 Ti}, {ht, 12500 ht, ht}} /.
 {{s1__}, {s2__}} :> Table[Sin[x], {x, s1}, {t, s2}] // packedQ

True

Answer 4

[Edit note: This was merged from duplicate question. What are offered below are workarounds for the two-iterator Tables, in response to a comment inquiring about alternatives, so that Table will use internal optimizations. Why there is a need for them is explained in other answers. Here, the solutions address the requirement that the iterator limits have to be evaluated before they are passed to Table in order for packed arrays and the internal optimization to be used.]

Case: iterator symbols undefined

Here is one option if x and t are undefined.

tab[fn_, rest__] := Table[fn, rest];
SetAttributes[tab, HoldFirst];
a3 = tab[f[t, x], {x, 0.5 Ti, 125 Ti, 0.5 Ti}, {t, ht, 12500 ht, ht}]; // AbsoluteTiming
(* {0.541578, Null} *)

Another option is to use With:

a4 = With[{iter1 = {x, 0.5 Ti, 125 Ti, 0.5 Ti}, iter2 = {t, ht, 12500 ht, ht}}, 
    Table[f[t, x], iter1, iter2]]; // AbsoluteTiming
(* {0.552246, Null} *)

Case: blocking symbols that have been defined

If x or t has a defined value, you can use Block:

a5 = Block[{x, t},
       With[{iter1 = {x, 0.5 Ti, 125 Ti, 0.5 Ti}, iter2 = {t, ht, 12500 ht, ht}}, 
        Table[f[t, x], iter1, iter2]]
     ]; // AbsoluteTiming
(* {0.543405, Null} *)

Block can be wrapped around tab, too, but there's a better way.

Automatically blocking iterator symbols

Block can be incorporated in a custom table function:

SetAttributes[tab2, HoldAll];
tab2[fn_, rest__] := 
  With[{iterVars = Cases[Hold[rest], {sym_Symbol, __} :> Hold@sym]},
   Block[iterVars, Table[fn, ##] &[rest]] // ReleaseHold];

If we set x equal to something, the table function still generates the desired table.

x = 5;
a6 = tab2[f[t, x], {x, 0.5 Ti, 125 Ti, 0.5 Ti}, {t, ht, 12500 ht, ht}]; // AbsoluteTiming
(* {0.539895, Null} *)

They all give equivalent results:

a1 == a2 == a3 == a4 == a5 == a6
(* True *)

In the fast versions, because the limits are numbers, Table chooses a fast alternative that, as your question shows, uses packed arrays.