function [bounds, prodOpt, allProd, info] = jsr_prod_Gripenberg(M,varargin)
%
% JSR_PROD_GRIPENBERG(M) Computes bounds by a branch and bound method on
%                        products of growing length.
%
% [BOUNDS, PRODOPT, ALLPROD, INFO] = JSR_PROD_GRIPENBERG(M)
%     For M a cell array of matrices, computes products of growing length
%     and keeps only the products that could produce an upper bound higher
%     than current_lower + delta. See default values below.
%
% [BOUNDS, PRODOPT, ALLPROD, INFO] = JSR_PROD_GRIPENBERG(M,OPTS) or
% [BOUNDS, PRODOPT, ALLPROD, INFO] = JSR_PROD_GRIPENBERG(M,OPTSNAME,OPTSVALUE)
%     Uses the values of the parameters in the structure OPTS generated by
%     jsrsettings. Or the values given by pair of OPTSNAME, OPTSVALUE
%     arguments. OPTSNAME can be 'delta', 'maxEval' or 'normfun' for
%     instance.
%
%   BOUNDS  is a vector with the bounds : [lower, upper]
%
%   PRODOPT contains the indices of a product whose average spectral radius
%           attains the lower bound (use with buildProduct)
%   ALLPROD contains on each row the indices of the products kept
%
%   INFO is a structure containing various data about the iterations :
%         info.status         - = 0 if delta was achieved, 1 if maxEval was
%                               reached,2 if stopped in iteration
%                               (CTRL-C or OPTS.maxTime reached)
%         info.elapsedtime
%         info.timeIter       - cputtime - start time in s. at the end of
%                               each iteration. Note : diff(info.timeIter)
%                               gives the time taken by each iteration.
%         info.allLb          - lower bounds at each depth
%         info.allUb          - upper bounds at each depth
%         info.popHist        - number of candidates along depths
%         info.opts           - the structure of options used
%
%  The field OPTS.grip (generated by jsrsettings) can be used to
%  set parameters and options :
%
%      grip.delta         - aimed difference between bounds, influences the
%                           the number of kept products and hence the
%                           speed, (1e-2)
%      grip.maxEval       - maximum number of eigenvalues evaluation and
%                           norms calculated, (1000)
%      grip.normfun       - function handle to a matrix norm,
%                           (@norm)
%      grip.plotBounds    - if 1 plots the evolution of the bounds w.r.t.
%                           the length, (0)
%      grip.plotpopHist   - if 1 plots the evolution of the number of
%                           candidates, (0)
%      grip.plotTime      - if 1 plots evolution of the bounds found
%                           w.r.t. time of computation, (0)
%
% See also jsrsettings
%
% REFERENCES
%    G.Gripenberg,
%      "Computing the joint spectral radius",
%      Linear Algebra and its Applications, 234:43-60, 1996
%    R.Jungers,
%      "The Joint Spectral Radius: Theory and Applications"
%      Vol. 385 section 2.3.3 in Lecture Notes in Control and Information
%      Sciences, Springer-Verlag. Berlin Heidelberg, June 2009
%

% When scale is too low, the problem is hill-conditionned. We need to bound
% it below.


opts = jsrsettings;

if (nargin>1)
    iarg = 1;
    while (iarg<=length(varargin))
        if (ischar(varargin{iarg}))
            opts = jsrsettings(opts,['grip.' varargin{iarg}],varargin{iarg+1});
            iarg = iarg+2;
            
        elseif (isstruct(varargin{iarg}))
            opts = varargin{iarg};
            iarg = iarg+1;
        end
        
    end
end

% logfile opening
close = 1;
if (ischar(opts.logfile) )
    logFile = fopen(opts.logfile,'wt');
    if (logFile == -1)
        warning(sprintf('Could not open file %s',opts.logfile));
    end
elseif isnumeric(opts.logfile)
    if (opts.logfile==0)
        logFile = -1;
    elseif (opts.logfile == 1)
        logFile = fopen('log_grip','wt');
        if (logFile == -1)
            warning('Could not open logfile')
        end
    else
        logFile = opts.logfile;
        close =0;
    end
else
    logFile = fopen('log_grip','wt');
    if (logFile == -1)
        warning('Could not open logfile')
    end
end

if (logFile~=-1)
    fprintf(logFile,[datestr(now) '\n\n']);
end

msg(logFile,opts.verbose>1,'\n \n******** Starting jsr_prod_Gripenberg ******** \n \n')
tic;

% Parameters
m = length(M);
n = length(M{1});
normfun = opts.grip.normfun;
delta = opts.grip.delta;
maxeval = opts.grip.maxEval;

% Initialisation
status = 2;
neval = m;
l = 1;
nt = m;
prods(1:m,1) = (1:m)';
rhoM = rho(M);
[alpha, prodOpt] = max(rhoM);

% Sometimes, scale can be equal to zero
if(alpha < eps)
    scale = eps;
else
    scale = alpha;
end
beta = 0;
P = zeros(1,m);
popHist = zeros(1,50);
popHist(1) = m;
timeIter = zeros(1,50);
allLb = timeIter;
allUb = timeIter;
% Note: comparison between for and cellfun was performed, for is quicker in
% majority of tests, see testCellfun_for.m
for im=1:m
    normi = feval(normfun,full(M{im}));
    if (normi>beta)
        beta = normi;
    end
    P(im) = normi;
end

M = cellDivide(M,scale);
T = M;

allLb(1) = alpha;
allUb(1) = beta;

if (m==1)
    bounds = [alpha, alpha];
    prodOpt = 1;
    allProd = [];
    info.status = 0;
    info.elapsedtime = toc;
    info.timeIter = 0;
    info.allLb = allLb(1:l);
    info.allUb = allUb(1:l);
    info.popHist = popHist(1:l);
    info.opts = opts;
    return;
end

% Main loop
while (neval < maxeval && beta > alpha + delta)
    l = l+1;
    if (l>length(popHist))
        popHist = [popHist zeros(1,50)];
        timeIter = [timeIter zeros(1,50)];
        allLb = [allLb zeros(1,50)];
        allUb = [allUb zeros(1,50)];
    end
    
    newT = cell(1,m*nt);
    newProds = sparse(m*nt,l);
    newP = sparse(1,m*nt);
    newnt = 0;
    maxNewP = alpha;
    
    msg(logFile,opts.verbose>1,'Depth %d, starting computation of %d new products',l,m*nt);
    alphaprev=alpha;
    for it = 1:nt
        
        for im = 1:m
            MT = M{im}*T{it};
            normi = feval(normfun,full(MT))^(1/l)*scale;
            rhoMT = rho(MT)^(1/l)*scale;
            
            ilast = find(prods(it,:),1,'last');
            if isempty(ilast)
                ilast = 0;
            end
            prodi = prods(it,1:ilast);
            
            newPi = min(P(it),normi);
            
            % Actualise alpha
            if (rhoMT>alpha)
                alpha = rhoMT;
                prodOpt = [prodi im];
                
            end
            
            % Keep or not
            if (newPi>alpha+delta)
                newnt = newnt+1;
                newT{newnt} = MT;
                
                % Add to newProds
                newProds(newnt,1:ilast+1) = [prodi im];
                % Actualise newP
                newP(newnt) = newPi;
                
                % Keep max newP to actualise beta
                if (newP(newnt)>maxNewP)
                    maxNewP = full(newP(newnt));
                end
                
            end
        end
    end
    %  % Rescale the different sets
    if (alpha>alphaprev)
        newScale = alpha;
        if(newScale < eps)
        % Sometimes scale can be equal to zero
            newScale = beta;
        end
        
        M = cellDivide(M,newScale/scale);
        T = cellDivide(T,(newScale/scale)^(l-1));
        newT = cellDivide(newT,(newScale/scale)^(l));
        scale = newScale;
    end
    
    msg(logFile,opts.verbose>1,'Number of products kept : %d',newnt);
    neval = neval + m*nt;
    nt = newnt;
    T = newT(1:nt);
    prods = newProds(1:nt,:);
    P = newP(1:nt);
    beta = min(beta,max(alpha+delta,maxNewP));
    
    popHist(l) = nt;
    timeIter(l) = toc;
    allLb(l) = alpha;
    allUb(l) = beta;
    
    if (l>999)
        if (mod(l,100)==0)
            msg(logFile,opts.verbose>0,'\n>  Depth %3d  -  [%15.12g, %15.12g]  -  %d candidates\n', l, alpha, beta, nt);
        end
    elseif (l>150)
        if (mod(l,50)==0)
            msg(logFile,opts.verbose>0,'\n>  Depth %3d  -  [%15.12g, %15.12g]  -  %d candidates\n', l, alpha, beta, nt);
        end
    elseif (l>20)
        if (mod(l,10)==0)
            msg(logFile,opts.verbose>0,'\n>  Depth %3d  -  [%15.12g, %15.12g]  -  %d candidates\n', l, alpha, beta, nt);
        end
    else
        msg(logFile,opts.verbose>0,'\n>  Depth %3d  -  [%15.12g, %15.12g]  -  %d candidates\n', l, alpha, beta, nt);
    end
end

if (beta-alpha)<=delta
    status = 0;
else
    status=1;
end

elapsedtime = toc;
msg(logFile,opts.verbose>0,'\n> Bounds on the jsr : [%.15g, %.15g]', alpha, beta);

msg(logFile,opts.verbose>1,'\nEnd of algorithm after %5.2f s',elapsedtime)

if (logFile~=-1 && close)
    fclose(logFile);
end

% Post-processing
bounds = [alpha, beta];
prodOpt = full(prodOpt);
allProd = full(prods);
info.status = status;
info.elapsedtime = elapsedtime;
info.timeIter = timeIter(1:l);
info.allLb = allLb(1:l);
info.allUb = allUb(1:l);
info.popHist = popHist(1:l);
info.opts = opts;

if(opts.grip.plotBounds)
    figure
    plot(1:l,info.allLb,'-+g',1:l,info.allUb,'-*r')
    title('grip : Bounds found at each depth')
    xlabel('Length of products')
end

if opts.grip.plotpopHist
    figure
    plot(1:l,info.popHist,'-o')
    title('grip : Evolution of number of candidates')
    xlabel('Length of products')
end


if (opts.grip.plotTime)
    figure
    if (info.timeIter(end) > 600)
        time = info.timeIter/60;
        plot(time,info.allLb,'-+g',time,info.allUb,'-*r','MarkerSize',10);
        xlabel('time from start in min')
    else
        plot(info.timeIter,info.allLb,'-+g',info.timeIter,info.allUb,'-*r','MarkerSize',10)
        xlabel('time from start in s')
    end
    title('grip : Bounds at each depth w.r.t. to time')
    legend('Upper bound','Lower bound')
end



end
