%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% construct extraction matrix for multiple tensor product               %
%% input: number of variables                                            %
%%        degree of tensor product                                       %
%% output: Matrix A such that A*kron(y,y...,y) has no repeated elements  %
%% example:                                                              %
%% >> syms a b c                                                         %
%% >> v=[a;b;c];                                                         %
%% >> A=construct_matrix(3,4);                                           %
%% >> expand(ones(1,15)*A*kron(v,kron(v,kron(v,v)))-(a+b+c)^4)           %
%% ans= 0;                                                               %
%% giovanni.lombardo@ecb.int @2010                                       %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function A=construct_matrix(nvar,order)

y=[];
for jj=1:nvar;
    y=[y;sym(char(65+jj))];
end
y3=kron(y,y);
for jj=1:order-2
y3=kron(y3,y);
end
% order=length(y3)^(1/nvar);
sz=nvar+(order-1)*nvar+1;
A=zeros(sz,length(y3));

[ii,jj]=unique(y3);
ysint=y3(sort(jj));
obj1=regexprep(sym2cell(y3),'\*','_');
obj1=regexprep(obj1,'\^','#');
obj2=regexprep((sym2cell(ysint)),'\*','_');
obj2=regexprep(obj2,'\^','#');
for kk=1:length(obj2);
[x,yy]=find_incell(obj1,obj2{kk});
A(kk,yy)=1;
end
return
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% EXAMPLE
% syms a b c d
% y=[a;b;c];
% y3=kron(kron(y,y),y);
% A*y3
%  
% ans =
%  
%      a^3
%  3*a^2*b
%  3*a*b^2
%      b^3
%  3*a^2*c
%  6*a*b*c
%  3*b^2*c
%  3*a*c^2
%  3*b*c^2
%      c^3