The Cartesian product of a sequence of m lists
L1, ..., Ln,
with list i having ni elements,
is given by the sequence of Πi ni lists
[L11,...,Lm1], ...,
[L1n1,...,Lmnm],
where Ljk is the k-th element of the j-th list.
If the right-most element in the lists varies the fastest,
the sequence is in row-major order,
if the left-most elements varies the fastest,
the sequence is in column-major order.
With one exception (noted) all the routines given below use
row-major order.
Because a Cartesian product can be large,
Maple provides a procedure, combinat[cartprod],
that generates an iterator for the Cartesian product.
Each call to this iterator returns the next value in the sequence.
Here is a procedure that uses cartprod
to iterate through each list in the Cartesian product of its arguments,
passing each generated list to the function f.
Note:
This procedure, as well as most of the following,
uses the seq parameter modifier,
a mechanism introduced in Maple 11.
The parameter L is assigned
an expression sequence consisting of the following lists.
Here it is equivalent to using the older
args[2..-1], provided the remaining arguments
are lists.
use_combinat := proc(f, L::seq(list))
local T;
uses combinat;
T := cartprod([L]);
while not T['finished'] do
f(T['nextvalue']());
end do;
NULL;
end proc:
By assigning f to print its argument,
we can demonstrate that this works as intended:
f := curry(printf," %a"):
lsts := [$1..3],[a,b]:
use_combinat(f,lsts): printf("\n"):
[1, a] [1, b] [2, a] [2, b] [3, a] [3, b]
Using combinat[cartprod] is a straight-forward and
effective way to step through the entire set without having to
allocate memory for the set. The iterator itself, however, requires
some time to execute and if the loop itself is otherwise
fast—say it calls a compiled procedure—then the time spent
evaluating the iterator can be significant.
The Maple built-in procedure rtable_scanblock
applies a function to each index of an rtable.
If each dimension of the rtable is a range that corresponds to
the number of elements in each of the lists,
then there is a simple mapping from the indices of the rtable
to the sublists in the Cartesian product.
We can use this, then, to effectively iterate over the
Cartesian product. Because the rtable is
created with the storage=sparse option,
the memory allocation is minimal.
Here is a test procedure that implements this technique.
use_scanblock := proc(f)
local a,dims,n,L;
L := args[2..-1]:
n := nargs-1;
dims := seq(1..nops(a), a in L);
rtable_scanblock(rtable(dims, 'storage=sparse')
, [dims]
, proc(val,indx,data)
local i;
f([seq(L[i][indx[i]], i=1..n)]);
end proc
);
NULL;
end proc:
use_scanblock(f,lsts): printf("\n"):
[1, a] [1, b] [2, a] [2, b] [3, a] [3, b]
The following tabulates, for each procedure,
the time and memory required to iterate through each element
of a Cartesian product of lists.
Four different inputs are tested.
In each case, the rtable_scanblock
technique is faster (about 2×) and uses less memory
than combinat[cartprod].
Iterate through the Cartesian product of N lists of E elements,
generating E^N iterations. At each step, call the (unassigned)
global function f with the generated list.
gcfreq = 10^7
+-------- timing -------+--------- memory ---------+
256^2 | maple (s) | total (s) | used (MiB) | alloc (MiB) | gctimes
---------------+-----------+-----------+------------+-------------+--------
use_combinat | 1.77 | 1.90 | 29.98 | 23.75 | 0
use_scanblock | 0.93 | 0.95 | 19.23 | 13.87 | 0
---------------+-----------+-----------+------------+-------------+--------
10^5 | | | | |
---------------+-----------+-----------+------------+-------------+--------
use_combinat | 3.38 | 3.47 | 47.84 | 31.74 | 1
use_scanblock | 1.84 | 1.93 | 28.49 | 22.68 | 0
---------------+-----------+-----------+------------+-------------+--------
3^10 | | | | |
---------------+-----------+-----------+------------+-------------+--------
use_combinat | 2.73 | 2.81 | 37.70 | 31.49 | 0
use_scanblock | 1.40 | 1.45 | 21.51 | 16.06 | 0
---------------+-----------+-----------+------------+-------------+--------
2^16 | | | | |
---------------+-----------+-----------+------------+-------------+--------
use_combinat | 3.99 | 4.13 | 49.50 | 32.62 | 1
use_scanblock | 2.06 | 2.08 | 25.63 | 21.00 | 0
---------------------------------------------------------------------------
While iterators are useful for reducing the memory footprint, at
times it is necessary to instantiate the Cartesian product as a single
entity. Maple, alas, does not provide a function for doing so
(this might make be a nice addition to the ListTools package).
It is not difficult to create a procedure, the interesting part is making
it fast.
The following
site
has a few examples of implementations; I have used them below with minor
modifications, along with some of my own contributions.
The first example uses combinat[cartprod],
called Πi ni times,
to generate the Cartesian product.
CartProdDeiterated := proc(L::seq(list))
local C,i,S;
C := combinat['cartprod']([L]);
[seq(C['nextvalue'](), i=1..mul(nops(S), S=[L]))]
end proc:
The next example uses a modification of the rtable_scanblock
technique. Instead of calling rtable_scanblock,
it generates an rtable with an initializer that
maps each index to the corresponding entry in the Cartesian product.
The elements of the rtable are then converted to a list by using
seq.
CartRtable := proc(L::seq(list))
local i,a,n;
n := nargs;
[seq(a, a in rtable(seq(1..nops(a), a in args)
,() -> [seq(L[i][args[i]],i=1..n)]
,'order = C_order'
))];
end proc:
Here is an offering from Carl DeVore (aka Carl Love). I've
modified it slightly to use the seq parameter
modifier previously mentioned.
This routine produces the Cartesian product in column-major order;
that is, the left-most element in the lists varies the fastest.
All other routines produce the elements in row-major order.
CartProdAsList := proc(L::seq(list))
option `Copyright (C) 2002, Carl DeVore. All rights reserved.`;
local q,i,N;
N := map(nops, [L]);
[seq([seq(L[i][1+irem(q,N[i],'q')], i= 1..nargs)], q=0..`*`(N[])-1)]
end proc:
Here is a suggestion from Helmut Kahovec, also slightly modified.
CartProdRecursive := proc(L::seq(list))
local i,j;
option `Copyright (C) 2001, Helmut Kahovec. All rights reserved.`;
if nargs=2 then
[ seq(seq([`if`(i::list,i[],i),
`if`(j::list,j[],j)
] , j in args[2] )
, i in L[1] ) ];
elif nargs>2 then
procname( procname(args[1..2]), args[3..-1] );
else
args;
end if;
end proc:
The previous code has a few minor flaws and permits a few
optimizations. Here is a non-recursive version that corrects the
flaws and improves its performance.
CartProdNonrecursive := proc(L::seq(list))
local i,j,X,lst;
option `Copyright (C) 2007, Joseph Riel. All rights reserved.`;
if nargs=0 then return [];
elif nargs=1 then return map(`[]`,L);
end if;
X := [seq(seq([i,j], j in L[2] ), i in L[1] )] ;
for lst in L[3..-1] do
X := [seq(seq([i[],j], j in lst ), i in X )];
end do;
X;
end proc:
Here is an interesting version that I originally devised
as a way to write an add-like procedure
that summed over multiple dimensions, converted here to
use seq to generate the Cartesian product.
Because foldl is not a built-in,
it is possible to improve it further, however, the resulting
speed gains were trivial (less than five percent).
CartProdSeq := proc(L::seq(list))
local Seq,i,j;
option `Copyright (C) 2007, Joseph Riel. All rights reserved.`;
eval([subs(Seq=seq, foldl(Seq
, [cat(i,1..nargs)]
, seq(cat(i,j)=L[j],j=nargs..1,-1)
))]);
end proc:
To better understand what it does, run it without the eval
statement, or do
test := subs(eval=''eval'',op(CartProdSeq)):
test([1,2,3],[a,b]); %;
eval([seq(seq([i1, i2], i2 = [a, b]), i1 = [1, 2, 3])])
[[1, a], [1, b], [2, a], [2, b], [3, a], [3, b]]
The following tabulates, for each procedure, the time and memory
required to generate the Cartesian product of a sequence of lists.
Four different inputs are tested.
For all cases the CartProdSeq is the fastest
and uses the least memory.
Compute the cartesian product of N lists of E elements
to get a list of E^N sublists of N elements
gcfreq = 10^7
+-------- timing -------+--------- memory ---------+
256^2 | maple (s) | total (s) | used (MiB) | alloc (MiB) | gctimes
---------------------+-----------+-----------+------------+-------------+--------
CartProdDeiterated | 1.30 | 1.33 | 23.06 | 19.50 | 0
CartProdRtable | 0.54 | 0.55 | 10.00 | 7.50 | 0
CartProdAsList | 0.77 | 0.86 | 14.81 | 10.81 | 0
CartProdRecursive | 0.25 | 0.26 | 7.79 | 5.75 | 0
CartProdNonrecursive | 0.13 | 0.13 | 5.77 | 3.75 | 0
CartProdSeq | 0.12 | 0.13 | 5.76 | 3.75 | 0
---------------------+-----------+-----------+------------+-------------+--------
10^5 | | | | |
---------------------+-----------+-----------+------------+-------------+--------
CartProdDeiterated | 2.82 | 2.94 | 43.31 | 32.74 | 1
CartProdRtable | 1.30 | 1.36 | 19.85 | 14.93 | 0
CartProdAsList | 1.94 | 1.98 | 27.67 | 21.93 | 0
CartProdRecursive | 0.74 | 0.74 | 18.41 | 13.56 | 0
CartProdNonrecursive | 0.47 | 0.49 | 14.78 | 9.94 | 0
CartProdSeq | 0.32 | 0.35 | 13.53 | 9.37 | 0
---------------------+-----------+-----------+------------+-------------+--------
3^10 | | | | |
---------------------+-----------+-----------+------------+-------------+--------
CartProdDeiterated | 2.28 | 2.34 | 31.04 | 27.99 | 0
CartProdRtable | 1.08 | 1.10 | 12.54 | 10.56 | 0
CartProdAsList | 1.85 | 1.87 | 24.33 | 20.68 | 0
CartProdRecursive | 0.66 | 0.73 | 18.49 | 13.37 | 0
CartProdNonrecursive | 0.43 | 0.46 | 14.54 | 10.19 | 0
CartProdSeq | 0.24 | 0.27 | 10.56 | 8.25 | 0
---------------------+-----------+-----------+------------+-------------+--------
2^16 | | | | |
---------------------+-----------+-----------+------------+-------------+--------
CartProdDeiterated | 3.67 | 3.80 | 44.03 | 35.24 | 1
CartProdRtable | 1.64 | 1.67 | 17.02 | 14.56 | 0
CartProdAsList | 3.06 | 3.18 | 35.65 | 30.74 | 0
CartProdRecursive | 1.15 | 1.31 | 28.41 | 23.31 | 0
CartProdNonrecursive | 0.78 | 0.83 | 23.88 | 19.25 | 0
CartProdSeq | 0.44 | 0.48 | 18.21 | 14.00 | 0
---------------------------------------------------------------------------------
