Unique combinations of matrices

Given a cell array that consists of several matrices of various sizes, I'd like to find all unique combinations of these matrices. I am aware that this sort of question has been asked many times over, but all queries seem to seek combinations of (row) vectors. I am also aware that this problem can explode to be computationally expensive very fast, but I am hoping to extend this at manageable levels.

Here is a minimum working example:

T = cell(1,4);
T{1} = [0 0];
T{2} = [1 0; 0 1];
T{3} = [1];
T{4} = [1 1 0; 1 0 1; 0 1 1];

So, T is a 1x4 cell array consisting of four matrices of various sizes. I would like to find the matrix, M, that is the concatenation of all combinations of these four matrices. Each of these four matrices is considered one item when these combinations are considered. Thus, M would be of size

[prod(a), sum(b)]

where

a = cellfun('size', T, 1);
b = cellfun('size', T, 2);

In this example, M would consist of 6-rows and 8-columns (because the total number of permuted rows is 1*2*1*3 = 6 = 3!, and because the total number of combined columns of the matrices of T sums to 8). Via concatenation, M would be equal to:

M = [0 0 1 0 1 1 1 0; 0 0 1 0 1 1 0 1; 0 0 1 0 1 0 1 1; 0 0 0 1 1 1 1 0; 0 0 0 1 1 1 0 1; 0 0 0 1 1 0 1 1];
% or
M =

     0     0     1     0     1     1     1     0
     0     0     1     0     1     1     0     1
     0     0     1     0     1     0     1     1
     0     0     0     1     1     1     1     0
     0     0     0     1     1     1     0     1
     0     0     0     1     1     0     1     1

Using the results of previous vector-combinations queries (see: allcomb and Generate a matrix containing all combinations of elements taken from n vectors), I have tried:

C = allcomb(S{:}); % returns not enough columns and non-unique rows (size is 72x4)
C = unique(allcomb(S{:}), 'rows'); % returns not enough columns or rows (size is 4x4)

and

nCells = numel(S);
combs = cell(1,nCells);
[combs{end:-1:1}] = ndgrid(S{end:-1:1});
combs = cat(nCells+1, combs{:});
combs = reshape(combs, [], nCells); % size here is 74x4
combs = unique(combs, 'rows'); % size here is 4x4

Neither of these attempts returns the correct number of columns or the correct number of ("unique") rows. I believe that I have made a good start and that I am very close. What am I missing? What am I overlooking? Thanks!

Answer 1

[Guillaume]

T = {[0 0], [1 0; 0 1], [1], [1 1 0; 1 0 1; 0 1 1]}

rowset = cellfun(@(t) 1:size(t, 1), T, 'UniformOutput', false);  %row indices of each matrix
[rowset{:}] = ndgrid(rowset{:});                                 %all combinations of row indices
M = cell2mat(cellfun(@(row, t) t(row, :), rowset, T, 'UniformOutput', false))

Note: if you want to visualise what is in rowset (the row indices combination for each matrix):

reshape(cat(numel(T)+1, rowset{:}), [], numel(T))

However, keeping rowset as a cell array makes the next step easier.

[Matthew]

@Guillaume: I like this answer, but it doesn't seem to give me the same result that I am anticipating. Using these commands, I get:

M =

     0     0     1     0     1     1     1     0
     0     0     0     1     1     1     1     0
     0     0     1     0     1     1     0     1
     0     0     0     1     1     1     0     1
     0     0     1     0     1     0     1     1
     0     0     0     1     1     0     1     1

which is not the same as (see above):

M =
   0     0     1     0     1     1     1     0
   0     0     1     0     1     1     0     1
   0     0     1     0     1     0     1     1
   0     0     0     1     1     1     1     0
   0     0     0     1     1     1     0     1
   0     0     0     1     1     0     1     1

I realize that this may just be a sorting issue, but I wonder if there is a way to update your sequence of commands so that the combinations/permutations of the matrix of T are executed from left to right?

In other words, so that the matrix M is constructed by finding and concatenating the permutations of T{1}, and then T{2}, and then T{3}, and then T{4} (etc)?

[Guillaume]

That's due to the order in which ndgrid compute the cartesian product. You just have to swap the order of the inputs to ndgrid and swap back the outputs to get it the other way round. Possibly, the most compact way of doing this:

rowset = fliplr(cellfun(@(t) 1:size(t, 1), T, 'UniformOutput', false));
[rowset{:}] = ndgrid(rowset{:});
M = cell2mat(cellfun(@(row, t) t(row, :), fliplr(rowset), T, 'UniformOutput', false))