function kernels = vkrn(F,G,b,c,N,T,m,P,msk)
%VKRN computes triangular Volterra kernels from the Carleman
%Bilinearization of a system
%
% KERNELS = VKRN(F,G,B,C,N,T,P) computes T-spaced samples up to time N*T of
% the first P<=5 Volterra triangular kernels of the bilinear system
% [F,G,b,c] obtained by Pth order Carleman Bilinearization of a
% linear-analytic sistem. The kernels are returned in vector form.
%
% KERNELS = VKRN(F,G,B,C,N,T,P,msk) uses a cell of P pre-calculated masks
% (see util/mask_P*_N*.m)

% See also CARLBILIN and CARLBL
%
%   CBVK Toolbox
%   JHMG, PMSB 2016

warning('cbvk:vkrn','Consider the use of vkrnsvdx, which is probably faster');


% Vector initialization
h1 = zeros(1,N);
if P > 1
    h2 = zeros(1,nchoosek(N+1,2));
    if P > 2
        h3 = zeros(1,nchoosek(N+2,3));
        if P > 3
            h4 = zeros(1,nchoosek(N+3,4));
            if P > 4
                h5 = zeros(1,nchoosek(N+4,5));
            end
        end
    end
end

% Pre-computed matrices
eF = cell(1,N);
eT = expm(F*T);

%---------------------------------------------------
% This is for greater precision as P grows
eF1=expm(F(1:m,1:m)*T);
if P>2
    M2=m+nchoosek(m+2-1,2);
    eF2=expm(F(1:M2,1:M2)*T);
    if P>3
        M3=M2+nchoosek(m+3-1,3);
        eF3=expm(F(1:M3,1:M3)*T);
        if P>4
            M4=M3+nchoosek(m+4-1,4);
            eF4=expm(F(1:M4,1:M4)*T);
        end
    end
end
if P==2
    eT(1:m,1:m)=eF1;
elseif P==3
    eF2(1:m,1:m)=eF1;
    eT(1:M2,1:M2)=eF2;
elseif P==4
    eF2(1:m,1:m)=eF1;
    eF3(1:M2,1:M2)=eF2;
    eT(1:M3,1:M3)=eF3;
elseif P==5
    eF2(1:m,1:m)=eF1;
    eF3(1:M2,1:M2)=eF2;
    eF4(1:M3,1:M3)=eF3;
    eT(1:M4,1:M4)=eF4;
end
%------------------------------------

eF{1} = eye(size(F,1));
eF{2}=eT;
for k=3:N
  eF{k}=eF{k-1}*eT;
end


% Prepare output
kernels = cell(1,P);

% Kernel indices
i2=1;i3=1;i4=1;i5=1;

% Kernel masks
m2=1;m3=1;m4=1;m5=1;

if nargin <9
    % Main loop  with mask calculation
    for k1=1:N
        % Repetition counter for k indices
        rk = 1;
        v1=c*eF{k1};
        h1(k1)=v1*b;
        if P>1
            v1 = v1*G;
            for k2=k1:N
                if(k2==k1)
                    rk = rk + 1;
                else
                    rk = 1;
                    % Mask for h2
                    m2 = 1/rk;
                end
                v2=v1*eF{k2-k1+1};
                h2(i2) = m2 * v2*b;
                i2=i2+1;
                if P > 2
                    v2 = v2*G;
                    for k3=k2:N
                        if(k3==k2)
                            rk = rk + 1;
                        else
                            rk = 1;
                            % Mask for h3
                            m3 = m2/rk;
                        end
                        v3=v2*eF{k3-k2+1};
                        h3(i3) = m3 * v3*b;
                        i3=i3+1;
                        if P > 3
                            v3 = v3*G;
                            for k4=k3:N
                                if(k4==k3)
                                    rk = rk + 1;
                                else
                                    rk = 1;
                                    m4 = m3/rk;
                                end
                                v4=v3*eF{k4-k3+1};
                                h4(i4) = m4 * v4*b;
                                i4=i4+1;
                                if P > 4
                                    v4=v4*G;
                                    for k5=k4:N
                                        if(k5==k4)
                                            rk = rk + 1;
                                        else
                                            rk = 1;
                                            m5 = m4/rk;
                                        end
                                        v5=m5 * v4*eF{k5-k4+1};
                                        h5(i5) = v5*b;
                                        i5=i5+1;
                                    end
                                end
                            end
                        end
                    end
                end
            end
        end
    end
    kernels{1}=h1;
    if P > 1
        kernels{2} = h2;
        if P > 2
            kernels{3} = h3;
            if P > 3
                kernels{4} = h4;
                if P > 4
                    kernels{5} = h5;
                end
            end
        end
    end

else
    % Main loop without mask calculation
    for k1=1:N
        v1=c*eF{k1};
        h1(k1)=v1*b;
        if P>1
            v1 = v1*G;
            for k2=k1:N
                v2=v1*eF{k2-k1+1};
                h2(i2) = v2*b;
                i2=i2+1;
                if P > 2
                    v2 = v2*G;
                    for k3=k2:N
                        v3=v2*eF{k3-k2+1};
                        h3(i3) = v3*b;
                        i3=i3+1;
                        if P > 3
                            v3 = v3*G;
                            for k4=k3:N
                                v4=v3*eF{k4-k3+1};
                                h4(i4) = v4*b;
                                i4=i4+1;
                                if P > 4
                                    v4=v4*G;
                                    for k5=k4:N
                                        v5=v4*eF{k5-k4+1};
                                        h5(i5) = v5*b;
                                        i5=i5+1;
                                    end
                                end
                            end
                        end
                    end
                end
            end
        end
    end
    kernels{1}=h1;
    if P > 1
        kernels{2} = h2.*msk{2};
        if P > 2
            kernels{3} = h3.*msk{3};
            if P > 3
                kernels{4} = h4.*msk{4};
                if P > 4
                    kernels{5} = h5.*msk{5};
                end
            end
        end
    end
end


    
    

