function A=multiindex(m,p,combine)
% MULTIINDEX Generate a table of multiindices (block-scheme).
%   A=MULTIINDEX(M,P,COMBINE) generate a table of multiindices using 
%   the standard block scheme i.e. generating all multi-indices up to
%   degree P in all M (random) variables. (Limitation to certain
%   limiters/norms will be added later). If combine is not specified or
%   evaluates to true then the homogeneous multiindices will be combined 
%   into one large (sparse) array A. Otherwise A is a cell array
%   where A{i} represents the multiindices with degree i-1.
%
% Example
%   % To generate the polynomial chaos for 2 random variables up to
%   % polynomial order 4 
%   A=multiindex(2,4); 
%   full(A); % convert from sparse
% 
%   % To generate the polynomial chaos for 5 random variables up to
%   % polynomial order 3, using only the homogeneous chaos of order 3 
%   A=multiindex(5,3,false); 
%   Ahom=full(A{4});
%
% See also MULTIINDEX_ORDER, MULTIINDEX_FACTORIAL

%   Elmar Zander
%   Copyright 2006, Institute of Scientific Computing, TU Braunschweig.
%   $Id: multiindex.m,v 1.1 2006/12/08 08:51:29 krosche Exp $ 

if nargin<3 
    combine=true;
end

% TODO: comment the code, idea of the algorithm

if m==1
    A=cell(1,p+1);
    for d=0:p
        A{d+1}=sparse(1,m,d);
    end
else
    [nonzero,rows]=monostats(m,p);
    A=cell(1,p+1);
    for d=0:p
        A{d+1}=spalloc(rows(d+1),m,nonzero(d+1));
    end

    Am=multiindex(m-1,p,false);
    for d=0:p
        s=0;
        A{d+1}( :, 1:end-1 ) = cat(1,Am{(d+1):-1:1});
        for dm=d:-1:0
            l=size(Am{dm+1},1);
            A{d+1}( (s+1):(s+l), end ) = d-dm;
            s=s+l;
        end
        
        %disp(sprintf('%6d %6d', [nonzero(d+1) nnz(A{d+1}) rows(d+1) size(A{d+1},1)]));
    end
end

if combine
    A=cat(1,A{:});
end




function [nonzero,rows]=monostats(m,p)
rows=ones(1,p+1);
nonzero=[0 ones(1,p)];

for k=2:m
    for d=p:-1:0
        rows(d+1)=sum(rows(1:(d+1)));
        nonzero(d+1)=sum(nonzero(1:(d+1))) + sum(rows(1:(d)));
    end
end

