function poly_g = alpha2gamma_HPPD(poly,G)
% function poly_g = alpha2gamma_HPPD(poly,G)
%
% Converts a homogeneous polynomial parameter-dependent (HPPD) matrix of arbitrary degree in alpha and/or 
% delta-alpha domain to a HPPD matrix in gamma-domain. The relalation between alpha, delta-alpha and gamma is 
% linear.  This program is an auxiliarly routine necessary for the methods developed in (Oliveira & Peres, 
% Automatica, 2010).
%
% inputs:  poly     : homogeneous polynomial to be converted
%          poly.mon : monomials of the homogeneous polynomial ([mon_1;mon_2;...mon_last])
%          poly.coef: coefficients of the monomials (in cells)
%          poly.N   : number of variables of the polynomial
%          poly.deg : degree of the polynomial
%          G        : matrix defining the linear relation between alpha and gamma        
%
% outputs: poly_g   : converted homogeneous polynomial

% Date: 26/06/2008
% Author: ricfow@dt.fee.unicamp.br


N=poly.N;
M = size(G,2);

for i=1:N
    for j=1:poly.deg
        C{i,j} = alpha_i_power_g(M,N,G(i,:),j);
    end
end

gamma_mon = gen_coefs(M,poly.deg);
for i=1:size(gamma_mon,1)%for each monomial of the gamma-basis polynomial
    mon_g = gamma_mon(i,:);
    
    poly_g.coef{i} = 0;
    for j=1:size(poly.mon,1)%collect the contributions of each alpha-basis monomials
        mon_a = poly.mon(j,:);
        c = extract_coefficient(mon_g,mon_a,C);
        poly_g.coef{i} = poly_g.coef{i} + c*poly.coef{j};
    end    
end
poly_g.mon = gamma_mon;

%_____________________________________________________________________________________________________________
function c = extract_coefficient(mon_g,mon_a,C)

cont=1;
cas = [];
for i=1:size(mon_a,2)
    if mon_a(1,i) > 0
        other_mons{cont} = C{i,mon_a(1,i)};
        cont=cont + 1;
        cas = [cas i];
    end
end
R = find_solutions(mon_g,other_mons);
c=0;
for i=1:size(R,2)
    ind = R{i}.ind';
    cnow = 1;
    for j=1:size(ind,2)
        cnow = cnow * C{cas(j),mon_a(1,cas(j))}.coef(ind(1,j),1);
    end
    c = c + cnow;    
end


%_____________________________________________________________________________________________________________
function R = find_solutions(mon,other_mons)

nsol = 1;
n_mon = size(other_mons,2);
for i=1:n_mon
    H{i} = 1:size(other_mons{i}.mon,1);
end
C = all_possible_comb(n_mon,H);
for i=1:size(C,2)
    sum_now = zeros(1,size(mon,2));
    sol = [];
    for j=1:size(C,1)
        sum_now = sum_now + other_mons{j}.mon(C(j,i),:);
        sol = [sol;other_mons{j}.mon(C(j,i),:)];
    end
    if sum(sum_now==mon) == size(mon,2)
        R{nsol}.mon = sol;
        R{nsol}.ind = C(:,i);
        nsol = nsol + 1;
    end
end

%_____________________________________________________________________________________________________________
function C = alpha_i_power_g(M,N,G_i,deg)

C.mon = gen_coefs(M,deg);
C.coef = [];
for i=1:size(C.mon,1)
    mon = C.mon(i,:);
    c = factorial(deg)/prod(factorial(mon));
    for j=1:M
        c = c*my_pow(G_i(1,j),mon(j));
    end    
    C.coef = [C.coef;c];
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 r = my_pow(a,e)

if e == 0
    r=1;
else
    r=a;
    for i=2:e
        r=r*a;
    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

