Getting directly a packed array of repeated packed arrays

Question

Consider this code:

Needs["Developer`"];
N@Range[5];
PackedArrayQ@%
Table[%%, 5];
PackedArrayQ@%

which prints True, False. One sees that the trivial way to create a list out of repetitions of some list does not give you automatically packed array. Is there any way to accomplish this on the fly?

Answer 1

Personally, I would prefer march's solution of using ConstantArray[] instead of the tweak I am about to show. As I noted, one should not be (re)setting system options willy-nilly, and especially if you don't know what you're doing. Thus, here is how one can localize the setting change I mentioned in the comments:

With[{copt = SystemOptions["CompileOptions"]},
     Internal`WithLocalSettings[SetSystemOptions["CompileOptions" ->
                                                 {"TableCompileLength" -> 1}],
     arr = Table[N @ Range[5], {5}],
     SetSystemOptions[copt]]];

Check:

Developer`PackedArrayQ[arr]
   True

Answer 2

Mr.Wizard has shown how memory performance is affected. I wish to illustrate various aspects of timing performance.

First, Developer`ToPackedArray@Table[] takes about the same time as ConstantArray[]:

v = Range@1*^6;  (* all examples with Mr.Wizard's (packed)  v  *)

With[{copt = SystemOptions["CompileOptions"]}, 
 Internal`WithLocalSettings[
  SetSystemOptions["CompileOptions" -> {"TableCompileLength" -> Infinity}], 
  Benchmark[{ConstantArray[v, {#}] &, Developer`ToPackedArray@Table[v, {#}] &}, 2^# &],
  SetSystemOptions[copt]]]
(*
  {{{1, 0.0058}, {2, 0.012}, {4, 0.045},  {6, 0.324}, {8, 1.2742}},
   {{1, 0.0058}, {2, 0.012}, {4, 0.0442}, {6, 0.31},  {8, 1.250}}}
*)

Second, this is because Table[v, {n}] is really fast, since it is only copying n pointers:

With[{copt = SystemOptions["CompileOptions"]}, 
 Internal`WithLocalSettings[
  SetSystemOptions["CompileOptions" -> {"TableCompileLength" -> Infinity}], 
  Benchmark[{ConstantArray[v, {#}] &, Table[v, {#}] &}, 2^# &],
  SetSystemOptions[copt]]]
(*
{{{1, 0.0058},    {2, 0.011},      {4, 0.043},     {6, 0.317},     {8, 1.4}},
 {{1, 9.4*10^-7}, {2, 9.89*10^-7}, {3, 1.2*10^-6}, {5, 2.3*10^-6}, {8, 0.000016}, ...}}
*)

Third, Table[v, {n}] is about twice as slow as ConstantArray[v, {n}], if n is greater than or equal to "TableCompileLength":

With[{copt = SystemOptions["CompileOptions"]}, 
 Internal`WithLocalSettings[
  SetSystemOptions["CompileOptions" -> {"TableCompileLength" -> 1}], 
  Benchmark[{ConstantArray[v, {#}] &, Table[v, {#}] &}, 
   2^# &, {1, 2, 4, 6, 8},
   TimeConstraint -> 20.],
  SetSystemOptions[copt]]]
(*
  {{{1, 0.0059}, {2, 0.012}, {4, 0.043}, {6, 0.312}, {8, 1.258}},
   {{1, 0.013}, {2, 0.023}, {4, 0.086}, {6, 0.61}, {8, 3.}}}
*)

Fourth, for I/O, there's no difference in read speed, but the write speed is a bit strange (to me). I put in some Developer`PackedArrayQ calls and turn on unpacking warnings to make sure there wasn't something I was missing. I timed both random and systematic access, and the results were opposite for the completely packed p and the top-level unpacked up.

read[arrayfn_] := Module[{d, a},
   a = arrayfn[];
   d = Dimensions[a];
   BlockRandom[
    SeedRandom[0, Method -> "MersenneTwister"];
    {AbsoluteTiming[
       Do[a[[RandomInteger[{1, d[[1]]}], RandomInteger[{1, d[[2]]}]]], {10^4}]; 
       Developer`PackedArrayQ@a],
     AbsoluteTiming[       (* different offsets j*k in each row *)
       Do[a[[j, j*k]], {k, 10^3}, {j, 10}]; Developer`PackedArrayQ@a]}
    ]
   ];

write[arrayfn_] := Module[{d, a},
  a = arrayfn[];
  d = Dimensions[a];
  BlockRandom[
   SeedRandom[0, Method -> "MersenneTwister"];
   {AbsoluteTiming[       (* different offsets j*k in each row *)
      Do[a[[RandomInteger[{1, d[[1]]}], RandomInteger[{1, d[[2]]}]]] = -1, {10^4}]; 
      Developer`PackedArrayQ@a],

    a = arrayfn[];    (* reset array *)
    AbsoluteTiming[
      Do[a[[j, j*k]] = -2, {k, 10^3}, {j, 10}]; 
      Developer`PackedArrayQ@a]}
   ]
  ]

Tests:

On["Packing"]  (* just to double check *)

read[ConstantArray[v, {50}] &]  (* packed === p *)
read[Table[v, {50}] &]          (* unpacked of packed === up *)
(*
  {{0.043588, True},  {0.009283, True}}
  {{0.043868, False}, {0.009309, False}}
*)

write[ConstantArray[v, {50}] &]  (* packed *)
write[Table[v, {50}] &]          (* unpacked of packed *)
(*
  {{0.056558, True},  {0.261569, True}}
  {{0.243272, False}, {0.037485, False}}
*)

Off["Packing"]

Why is the systematic overwriting of p so slow? (It makes no difference if you switch the orders of the iterators j, k.) If you comment out the reset-array line, then the timings are as follows:

write[ConstantArray[v, {50}] &]  (* packed *)
write[Table[v, {50}] &]          (* unpacked of packed *)
(*
  {{0.05759, True},   {0.014992, True}}
  {{0.244438, False}, {0.014401, False}}
*)

They're both the same, and twice as fast as the fastest with the array reset. If I switch the order of the systematic and random writes, then the packed array timing for the random writes depends on whether the array is reset. Maybe I'm doing something wrong, or maybe it's worth a separate question. Note that there is no difference in read if the array is reset; but there is a difference in speed between the random and systematic reads, which I guess is due to caching by the computer and not due to Mathematica.

Answer 3

Functions

march in a comment proposed ConstantArray, and indeed this works:

pQ = Developer`PackedArrayQ;
toP = Developer`ToPackedArray;
v = Range[5];

ConstantArray[v, 7] // pQ                    (* True *)

Padding functions will return packed arrays with some syntax but not others:

PadRight[{v}, {7, Automatic}, v] // pQ       (* True *)

PadRight[{{}}, {7, 5}, v]        // pQ       (* False *)

PadRight[{v}, 7, {v}]            // pQ       (* False *)

ArrayPad[{v}, {{0, 6}, 0}, v]    // pQ       (* True *)

Memory

Something to consider is that expressions are shared in memory and because of this performance may be quite different from what you expect!

Starting with a fresh kernel here is the memory performance of ConstantArray:

v = Range @ 1*^6;
ByteCount[v]
MaxMemoryUsed[]

p = ConstantArray[v, 50];
ByteCount[p]
MaxMemoryUsed[]

Developer`PackedArrayQ[p]

8000144
42129680

400000152
442130368

True

Our packed vector v takes ~8MB, and a packed array containing fifty copies of it takes ~400MB. The maximum memory use was ~442MB.

So how will Table used in a fashion that produces an unpacked outer list compare?

(* start with a fresh kernel *)

v = Range @ 1*^6;
ByteCount[v]
MaxMemoryUsed[]

u = Table[v, {50}];
ByteCount[u]
MaxMemoryUsed[]

Developer`PackedArrayQ[u]

8000144
42130176

400007648
42134816

False

Notice anything odd?

ByteCount again reports that our array is ~400MB, but the maximum memory used in the session is only ~42MB!

The actual memory use of the unpacked list (of packed vectors) can be a fraction of that of the equivalent fully packed array.

Of course if one starts modifying the individual packed elements of that list memory sharing may be compromised, but there can be very real advantages to not making a fully packed array in this context.

See also Addendum of:

• Fast method to select matrix elements based on a vector of positions