function poly_g = alpha2gamma_path(poly,G)
% function poly_g = alpha2gamma_path(poly,G)
%
% Converts a path-dependent polynomial matrix of arbitrary path length in alpha(k) and
% delta-alpha(k+r), r=1,...,L, to a HPPD matrix in gamma-domain. The relalation between alpha(k) and 
% delta-alpha(k+r), r=1,...,L, and gamma is linear.  This program is an auxiliarly routine necessary for the 
% methods developed in (Oliveira & Peres, Automatica, 2010).
%
% inputs:  poly       : path-dependent polynomial matrix to be converted
%          poly.paths : all possible paths ([path_1 path_2 ... path_last])
%          poly.coef  : coefficients of the paths (monomials) (in cells)
%          poly.N     : number of variables of the polynomial
%          poly.L     : path length
%          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;
L = size(G,2);%path size
M = size(G{1},2);%number of gammas

for i=1:N
    ALPH_SUCC{i,1} = G{1}(i,:);
    for j=2:L
        ALPH_SUCC{i,j} = ALPH_SUCC{i,j-1} + G{j}(i,:);
    end
end

gamma_mon = gen_coefs(M,L);
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.paths,1)%collect the contributions of each alpha-basis monomials
        path_a = poly.paths(j,:);
        c = extract_coefficient(mon_g,path_a,ALPH_SUCC);
        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,path_a,ALPH_SUCC)

for i=1:size(path_a,2)
    H{i} = eye(size(mon_g,2));
end
R = find_solutions(H,mon_g);
c=0;
for i=1:size(R,2)
    ind = R(:,i)';
    cnow = 1;
    for j=1:size(ind,2)        
        cnow = cnow * ALPH_SUCC{path_a(1,j),j}(1,ind(1,j));
    end
    c = c + cnow;    
end

%_____________________________________________________________________________________________________________
function Paths = all_paths(N,L)

Paths = [];
for i=1:L    
    Paths = [Paths kron(ones(1,N^(i-1)),kron(1:N,ones(1,N^(L-i))))'];
end

%_____________________________________________________________________________________________________________
function sols = find_solutions(sets,target)

var = [];
number_of_sets = size(sets,2);
sols = call_recur(number_of_sets,sets,1,var,[],target);
asf=1;

%_____________________________________________________________________________________________________________
function Resp = call_recur(total_sets,H,actual_set,var,Resp,target)

if actual_set <= total_sets    
    for i=1:size(H{actual_set},2)
        var{1}(actual_set,1) = i;        
        var{2}(actual_set,:) = H{actual_set}(:,i)';
        if actual_set == 1
            sum_ = var{2};
        else
            sum_ = sum(var{2});
        end
        if sum(sum_ <= target) == size(target,2)
            Resp = call_recur(total_sets,H,actual_set+1,var,Resp,target);
        end
    end
else
    if size(var{2},1) == 1
        sum_ = var{2};
    else
        sum_ = sum(var{2});
    end
    if sum(sum_ == target) == size(target,2)
        Resp = [Resp var{1}];        
    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
