nchoose(S)

function W = nchoose(S)
% NCHOOSE - all combinations of the elements of a set
%   W = nchoose(S) returns all possible combinations of one or more
%   elements of the set S. In total there are 2^N-1 combinations, where N
%   is the number of elements in S. W is a cell array: each cell holds one
%   of these combination (as a row vector). S can be a cell array, and each
%   cell of W will then contain a cell array.  
%
%   Example:
%      nchoose([2 4 6 8])
%        % -> { [2] ; 
%        %      [4] ; 
%        %      [2 4] ; 
%        %      [6] ;
%        %      ...
%        %      [2 6 8] ;
%        %      [4 6 8] ;
%        %      [2 4 6 8]} ; % in total (2^4)-1 = 15 different combinations
% 
%   Notes: 
%   - For sets containing more than 18 elements a warning is given, as this
%     can take some time. On my PC, a set of 20 elements took 20 seconds.
%     Hit Ctrl-C to intterupt calculations.
%   - If S contain non-unique elements (e.g. S = [1 1 2]), NCHOOSE will
%     return non-unique cells. In other words, NCHOOSE treats all elements
%     of S as being unique. One could use NCHOOSE(UNIQUE(S)) to avoid that.
%   - Loosely speaking, NCHOOSE(S) collects all output of multiple calls to
%     NCHOOSEK(S,K) where K is looping from 1 to the number of elements of
%     S. The implementation of NCHOOSE, however, does rely of a different
%     method and is much faster than such a loop.
%   - By adding an empty cell to the output, the power set of a cell array
%     is formed:
%       S = {1 ; 'hello' ; { 1 2 3}}
%       PowerS = nchoose(S) 
%       PowerS(end+1) = {[]}
%     See: http://en.wikipedia.org/wiki/Power_set
%   
%   See also NCHOOSEK, PERMS
%            COMBN, ALLCOMB on the File Exchange

% for Matlab R13 and up
% version 2.1 (may 2010)
% (c) Jos van der Geest
% email: [e-mail]

% History
% 1.0, may 2008 - inspired by a post on CSSM
% 2.0, may 2008 - improvemed algorithm
% 2.1, mar 2010 - added a note on "power set", updated ML version

% Acknowledgements:
% 2.0: Urs Schwarz, for suggesting a significant speed improvement using bitget


N = numel(S) ; 

% How many unique combinations of 1 to N elements are there
M = (2^N)-1 ;
if N > 18,
    warning('Nchoose:LargeOutput', ...
        'There are %d unique combinations. Please be patient ...',M) ;
end

S = S(:).' ;   % make the set a row vector, for uniform output

% The selection of elements is based on the binary representation of all
% numbers between 1 and M. The binary representation of numbers stored in X
% can be calculated using the formula (see DEC2BIN):
%    B1 = rem(floor(X(:) * pow2(1-N:0)),2) 
% (* means matrix multiplication) in which row B(k,:) represents X(k)
% Based on the suggestion by US, this formula can be rewritten as:
%    B2 = bitget(X2 * (2.^(N-1:-1:0)),N) > 0
% for a specific number X2. Although the bits of B2 are in reversed order
% (compared to B1), this is, however, exactly what we want, as we can now
% select the elements of S directly. We do not need to flip the order
% afterwards. 

% Looping over the numbers is faster than creating a whole matrix at once
W = cell(M,1) ;    % Pre-allocation of output
p2=2.^(N-1:-1:0) ; % This part of the formula can be taken out of the loop

for i=1:M,
    % calculate the (reversed) binary representation of i
    % select the elements of the set based on this representation
    W{i} = S(bitget(i*p2,N) > 0) ; 
end


% Two other options, which are less efficient
% % 1) looping over nchoosek (slow)
% for i=N:-1:1,
%     W{i} = nchoosek(S,i) ;
% end
% % and now W still has to be split into distinct cells ...
% 
% % 2) in one big step (memory inefficient)
% Q = rem(floor([1:M].'*pow2(1-N:0)),2).' ;
% [R,C] = find(Q) ;
% W = mat2cell(S(R),1,sum(Q)).' ;