function MSpoly = create_MSpolynomial(dim,m,var,deg,N,name_var,coeff)
% MSpoly = create_MSpolynomial(dim,m,var,deg,N,name_var,coeff)
% Creates a homogeneous polynomial matrix with parameters belonging to a multi-simplex.
%
% inputs:  dim           -> dimensions of the coefficient matrices.
%          m             -> number of simplexes inside the multi-simplex.
%          var           -> letters to be used for the variabels of the polynomial: 1-> a,
%                           2-> b, 3-> c, etc.
%          deg           -> partial degrees of the polynomial [g1 g2 ... gm].
%          N             -> number of variables in each simplex [N1 N2 ... Nm]. 
%          name_var      -> letter representing the polynomial matrix, ex: 'P' for P(alpha).
%          coeff         -> value to be assigned to the coefficients. The monomials that
%                           are not modified are considered as the identity matrix.
%          coeff.mon{i}  -> a monomial, given as a vector containing the powers of the
%                           variables; the alignment is as follows: [a b c ... ]=[a_1 .. 
%                           a_N1 b_1 ... b_N2 c_1 .... c_N3 ....]
%                           for instance: a_1^3a_2b_1^2b_2^3c_1c_2^4 -> [3 1 2 3 1 4].
%          coeff.value{i}-> coefficient value of the i-th monomial.
%          coeff.to_all  -> if the field 'to_all' exists, one value will be assigned 
%                           to all the coefficients of all the monomials.
%          coeff.val_all -> value to be assigned to all coefficients.
%
% outputs: MSpoly        -> resulting polynomial matrix
%
% Example: P(alpha) = 2*a_1b_1 + 3*a_1b_2 + 4*a_2b_1 + 5*a_2b_2
%                     
% with coefficients of order 1.
% dim = [1 1];
% m = 2;
% var = [1 2]; % a and b
% deg = [1 1];% affine in both variables
% N = [2 2];
% name_var = 'P';
% coeff.mon{1} = [1 0 1 0]'; % a_1b_1
% coeff.value{1} = 2;
% coeff.mon{2} = [1 0 0 1]'; % a_1b_1
% coeff.value{2} = 3;
% coeff.mon{3} = [0 1 1 0]'; % a_2b_1
% coeff.value{3} = 4;
% coeff.mon{4} = [0 1 0 1]'; % a_2b_2
% coeff.value{4} = 5;
% MSpoly = create_MSpolynomial(dim,m,var,deg,N,name_var,coeff)
%
% Date: 18/02/2008
% Author: ricfow@dt.fee.unicamp.br
%
 
MSpoly.m = m;
MSpoly.var = var;
MSpoly.deg = deg;
MSpoly.nxm = dim;
MSpoly.name = name_var;
MSpoly.N = N;
for i=1:MSpoly.m
    mon{i} = gen_coefs(MSpoly.N(1,i),MSpoly.deg(1,i))';
end
MSpoly.mon = all_possible_comb(MSpoly.m,mon);
MSpoly.nmon = size(MSpoly.mon,2);

for i=1:MSpoly.nmon
    
    if nargin == 7        
        if isfield(coeff,'to_all')
            if ~isfield(coeff,'val_all')
                fprintf('the general value for all the coefficients (coeff.val_all) should have been given !');
                MSpoly = [];
                return;
            end
            MSpoly.coef{i}.coef_val = coeff.val_all;
        else
            MSpoly.coef{i}.coef_val = get_coefficient(coeff,MSpoly.mon(:,i));
        end
    else
        if MSpoly.nxm(1) == MSpoly.nxm(2)
            c = eye(MSpoly.nxm(1));
        else
            c = ones(MSpoly.nxm);
        end
        MSpoly.coef{i}.coef_val = c*get_multinomial_coefficient(MSpoly.mon(:,i),MSpoly.N,MSpoly.deg);
    end
end

%________________________________________________________________________________________________________
function [comb] = gen_coefs(vertices,degree)
%generates all the solutions of: a_1 + a_2 + ... + a_vertices = degree

name = strcat('comb_N',int2str(vertices),'d',int2str(degree),'.mat');
if exist(name,'file') == 2
    load(name);
else
    comb = [];
    if vertices == 1
        comb = degree;
    else
        for i = 0:degree
            temp = gen_coefs(vertices-1,degree-i);
            comb   = [ comb
                i*ones(size(temp,1),1) temp ];
        end
    end
end

%________________________________________________________________________________________________________
function Resp = all_possible_comb(n,H)

var = [];
Resp = call_recur(n,H,1,var,[],1);

function Resp = call_recur(n,H,na,var,Resp,pos_ini)

if na <= n
    tam = size(H{na},1);
    for i=1:size(H{na},2)
        var(pos_ini:pos_ini+tam-1,1) = H{na}(:,i);
        Resp = call_recur(n,H,na+1,var,Resp,pos_ini+tam);
    end
else
    Resp = [Resp var];
end

%________________________________________________________________________________________________________
function value = get_coefficient(all_coefs,mon)

for i=1:size(all_coefs.mon,2)    
    if sum(all_coefs.mon{i} == mon) == size(all_coefs.mon{i},1)
        value = all_coefs.value{i};
        return;
    end
end

%________________________________________________________________________________________________________
function value = get_multinomial_coefficient(monomial,N,degs)

value = 1;
pos = 1;
for i=1:size(N,2)
    tuple_now = monomial(pos:pos+N(i)-1,:);
    v = factorial(degs(1,i))/prod(factorial(tuple_now));
    value = value*v;
    pos = pos + N(i);
end
