function [F,G,b,c] = carlbl(f,g,h,x,names,values,P,vb,wtylr,D,Dg,K)
%CARLBL computes the non-redundant Carleman Bilinearization of a
%linear-analytic system
%
% [F,G,B,C] = CARLBL(f,g,H,X,NAMES,VALUES,P,VB) computes the non-redundant
% Pth order Carleman Bilinearization v' = Fv + Gvu + bu, y = cv of the
% linear-analytic system x'= f(x) + g(x)u, y = hx. The procedure takes
% advantage of the sparsity of the matrices involved.
% 
% NAMES is a structure with the names of the parameters of f, g and h and
% VALUES is a structure (or a cell of structures with different sets) of
% parameter values. If VALUES is a cell then F,G and B are cells of
% matrices.
% 
% VB ~= 0 sets verbose mode on (default off)
%
% [F,G,B,C] = CARLBL(f,g,H,X,NAMES,VALUES,WTYLR) uses pre-calculated Taylor
% expansion coefficients
%
% [F,G,B,C] = CARLBL(f,g,H,X,NAMES,VALUES,WTYLR,D,Dg,K) uses also
% pre-calculated matrices in the bilinearization
%
% Modify this function accordingly if only a certain operation on the
% bilinear model for each set of parameters is required, and not returning
% the models themselves

% CBVK Toolbox
% JHMG, PMSB 2016

% ----------------------------------------------------

m=length(x);

if nargin<8
    vb=0;
end


% Taylor expansion weights
if nargin < 9
    if vb
        disp('Taylor expansion weights')
    end
    wtylr = mtaylorweights(m,P);
end

% CB matrices
if nargin < 10
    if vb
        disp('CB matrices')
    end
    [D,Dg,K] = precbsp(m,P);
end

% Non-redundant literal Taylor expansion using sparsity
if vb
    disp('Literal taylor expansion')
end
[F1i,F1c,F11] = multaylorsp(f,x,P);
[G1i,G1c,G11] = multaylorsp(g,x,P-1);


% Number of parameter sets
if iscell(values)
    NS = length(values);
else
    NS = 1;
end

MP=nchoosek(m+P,P)-1;
c=zeros(1,MP);
c(1:m)=h;


% Loop of parameter sets
for ns=1:NS
    
    if vb
        disp(['Parameter set ' num2str(ns)]);
        disp('Numerical evaluation of Taylor expansion')
    end
    % Numerical evaluation of Taylor expansion
    if ~iscell(values)
        [F1cn,F11n,G1cn,G11n]=evtaylorsp(m,F1c,F1i,F11,G1c,G1i,G11,names,values,wtylr);
    else
        [F1cn,F11n,G1cn,G11n]=evtaylorsp(m,F1c,F1i,F11,G1c,G1i,G11,names,values{ns},wtylr);
    end
    
    % Carleman Bilinearization
    if vb
        disp('Carleman Bilinearization')
    end
    [Fi,Fc,Gi,Gc] = carlbilinsp(F1i,F1cn,G1i,G1cn,G11n,P,K,D,Dg);

    % Return to matrix form
    % (C routine would probably be faster)
    if vb
        disp('Return to matrix form')
    end
    
    Fi=Fi';
    Fc=Fc';
    Gi=Gi';
    Gc=Gc';
    Fm=zeros(MP);
    Gm=zeros(MP);
    
    
    for i=1:MP
        ind=(Fi(i,:)~=0);
        Fm(i,Fi(i,ind))=Fc(i,ind);
        ind=(Gi(i,:)~=0);
        Gm(i,Gi(i,ind))=Gc(i,ind);
    end
    bm=Gm(:,1);
    Gm=[Gm(:,2:MP) zeros(MP,1)];
    
    if ~iscell(values)
        F=Fm; G=Gm; b=bm;
    else
        F{ns}=Fm; G{ns}=Gm; b{ns}=bm;
    end
    
    
end



