function MSAi = affine2MS_LTV(Atheta,bound_theta,bound_dtheta)
% function MSAi = affine2MS_LTV(Atheta,bound_theta,bound_dtheta)
% converts a time-varying polynomial matrix given in the affine form, i.e.
% A(theta)= A_0 + theta_1A_1 + ... theta_mA_m,  tmin_i < theta_i < tmax_i
% to a polynomial matrix given in the multi-simplex representation. The following change
% of variables is performed: t theta_i = (tmax_i-tmin_i)alpha_{i1}+tmin_i,  alpha_{i2}=1-alpha_{i1},
% i=1,..,m and subsequently, the resulting matrix is homogenized accordingly. The variables are
% assumed to be time-varying and their variation in time is bounded in the form: 
% dot_theta_i_min <= dot(theta_i) <= dot_theta_i_max. See (Oliveira, Bliman &
% Peres,CDC08, to appear) for details.
%
% inputs: Atheta      -> Matrix coefficients of the affine model [A0 A1 ... Am]
%         bound_theta -> bounds of the parameters: [theta_1_min theta_1_max;
%                                                   theta_2_min theta_2_max;
%                                                       ...             ....
%                                                   theta_m_min theta_m_max];
%         bound_dtheta -> bounds for the time-derivatives of the parameters:  
%                                                   [dot_theta_1_min dot_theta_1_max;
%                                                    dot_theta_2_min dot_theta_2_max;
%                                                       ...             ....
%                                                    dot_theta_m_min dot_theta_m_max];
%                         if the parameter is frozen, [0 0] must be given; if the parameter is
%                         arbitrarily time-varying, [-inf inf] must be given.
%
% outputs: MSAi  -> Homogeneous polynomial matrix in the multi-simplex representation with
%                   time-varying parameters.
%
%
% Example: affine model A(theta)=A0 + A_1theta_1(t) + A_2theta_2(t) + A_3theta_3, with
% A0 = [-2     1; -1    -3];
% A1 = [2     1; 0    -1];
% A2 = [1     1;-1     1 ];
% A3 = [0     0; 0    -1];
% Atheta = [A0 A1 A2 A3];
% bounds_theta = [-1 1; -1 1;-1 1];
% bounds_dtheta = [-1.5 1.5;-inf inf;0 0];
% MSAi = affine2MS_LTV(Atheta,bounds_theta,bounds_dtheta)
% 
% Date: 18/02/2008
% Author: ricfow@dt.fee.unicamp.br
%

N=size(bound_theta,1);
lin = size(Atheta,1);
col = size(Atheta,2)/(N+1);
if rem(col,1) ~= 0
    fprintf('Dimensions of the system vertices do not fit !\n');
    MSAi = [];
    return;
end

Ai_ini = Atheta;
Atheta = Ai_ini(:,1:col);
for i=1:N
    Atheta = Atheta +bound_theta(i,1)*Ai_ini(:,(i+1-1)*col+1:(i+1)*col);
end
for i=1:N
    Atheta = [Atheta (bound_theta(i,2)-bound_theta(i,1))*Ai_ini(:,(i+1-1)*col+1:(i+1)*col)];
end

for i=1:N
    MS{i} = create_MSpolynomial([lin lin],1,i,1,2,'I');    
    j= i+1;    
    coeff.mon{1} = [1; 0];
    coeff.value{1} = Atheta(:,(j-1)*col+1:j*col);
    coeff.mon{2} = [0; 1];
    coeff.value{2} = zeros([lin col]);
    MSs{i} = create_MSpolynomial([lin col],1,i,1,2,'I',coeff);    
end

coeff2.to_all =1;
coeff2.val_all = Atheta(:,1:col);
MSf{1} = create_MSpolynomial([lin col],1,1,1,2,'A0',coeff2);
for i=2:N
    MSf{1} = multiply_MSpoly(MS{i},MSf{1});
end

for i=1:N
    empty = 1;
    for j=1:N
        if j == i
            if empty == 1
                MSf{i+1} = MSs{j};
                empty = 0;
            else
                MSf{i+1}  = multiply_MSpoly(MSf{i+1},MSs{j});
            end
        else
            if empty == 1
                MSf{i+1} = MS{j};
                empty = 0;
            else
                MSf{i+1}  = multiply_MSpoly(MS{j},MSf{i+1});
            end
        end
    end
end

coeff2.val_all = zeros([lin col]);
MSAi = create_MSpolynomial([lin col],N,[1:N],ones(1,N),2*ones(1,N),'A',coeff2);

for i=1:MSAi.nmon
    coef_val = zeros([lin col]);
    mon_atual = MSAi.mon(:,i);
    for j=1:N+1
        for k=1:MSf{j}.nmon
            mon_comp = MSf{j}.mon(:,k);
            if sum(mon_atual == mon_comp) == size(mon_atual,1)
                coef_val = coef_val + MSf{j}.coef{k}.coef_val;
                break;
            end
        end
    end
    MSAi.coef{i}.coef_val = coef_val;
end

%arranging the time-derivatives bounds
for i=1:N
    MSAi.bounds{i}.H = [];
    MSAi.bounds{i}.columns = 0;
    if bound_dtheta(i,1) == 0 &  bound_dtheta(i,2) == 0
        MSAi.bounds{i}.H = zeros(MSAi.N(1,i),1);% convention for time-invariant parameter
        MSAi.bounds{i}.columns = 1;
    elseif bound_dtheta(i,1) == -inf &  bound_dtheta(i,2) == inf 
        MSAi.bounds{i}.H = inf*ones(MSAi.N(1,i),1);% convention for arbitrarily fast parameter
        MSAi.bounds{i}.columns = 1;        
    else        
        dot_am = [bound_dtheta(i,1)/(bound_theta(i,2)-bound_theta(i,1)) bound_dtheta(i,2)/(bound_theta(i,2)-bound_theta(i,1))];    
        H = gen_base(2,[dot_am;dot_am]);
        MSAi.bounds{i}.columns = 2;
        MSAi.bounds{i}.H = H;
    end
end

%________________________________________________________________________________________________________
function [H] = gen_base(N,bounds)
%function [H] = gen_base(N,bounds)
% This function generates the representation base H  used by the diffential
% inclusion method
% inputs: N       -> number of vertices of the polytope
%         bounds  -> time-derivative bounds. ex.: [[-1 1];[-1 1]...[-1 1]];
%
% output: H -> representation base

H = solve_plane(N,bounds,0);

%________________________________________________________________________________________________________
function H = solve_plane(N,limites,plano)

H = [];
for i=1:N
    for k=1:N
        if i == k
            continue;
        end
        sobra = [];
        pos = [];
        if N == 2
            if -limites(i,1) >= limites(k,1) & -limites(i,1) <= limites(k,2)
                vec_final(1,i) = limites(i,1);
                vec_final(1,k) = -limites(i,1);
                if ja_contem(H,vec_final') == 0
                    H = [H vec_final'];
                end
            end
            if -limites(i,2) >= limites(k,1) & -limites(i,2) <= limites(k,2)
                vec_final(1,i) = limites(i,2);
                vec_final(1,k) = -limites(i,2);
                if ja_contem(H,vec_final') == 0
                    H = [H vec_final'];
                end
            end
        else
            for o=1:N
                if o ~= i & o ~= k
                    pos = o;
                    sobra = limites(o,:);
                    break;
                end
            end
            for m=1:size(sobra',1);
                soma1 = -sum(sobra(:,m)) - limites(i,1) - plano;
                soma2 = -sum(sobra(:,m)) - limites(i,2) - plano;
                if soma1 >= limites(k,1) & soma1 <= limites(k,2)
                    vec_final(1,i) = limites(i,1);
                    vec_final(1,k) = soma1;
                    vec_final(1,pos) = sobra(:,m);

                    if ja_contem(H,vec_final') == 0
                        H = [H vec_final'];
                    end
                end
                if soma2 >= limites(k,1) & soma2 <= limites(k,2)
                    vec_final(1,i) = limites(i,2);
                    vec_final(1,k) = soma2;
                    vec_final(1,pos) = sobra(:,m);
                    if ja_contem(H,vec_final') == 0
                        H = [H vec_final'];
                    end
                end
            end
        end
    end
end

%________________________________________________________________________________________________________
function R = ja_contem(H,h)

R = 0;
for i=1:size(H,2)
    if sum(H(:,i) == h) == size(h,1)
        R = 1;
        return;
    end
end