function [bounds, P, info] = jsr_prod_pruningAlgorithm(M, varargin)

% JSR_PROD_PRUNINGALGORITHM Approximates the jsr of nonnegative matrices 
%                           using pruning algorithm.
%
%    [BOUNDS,P,INFO] = JSR_PROD_PRUNINGALGORITHM(M)
%      returns lower and upper bounds on the jsr of a set M of nonnegative
%      matrices using a pruning algorithm and considering products of length up
%      to 4. M must be a cell array of nonnegative matrices.
%
%    [BOUNDS,P,INFO] = JSR_PROD_PRUNINGALGORITHM(M,MAXEVAL,NORMFUN,DOMFUN)
%      returns lower and upper bounds on the JSR of a set M of nonnegative
%      matrices using a pruning algorithm. MAXEVAL limits the number of new
%      products generated.
%      The parameter NORMFUN is a function handle (default: @norm)
%      The parameter DOMFUN is a handle to the domination function:
%          domfun(A, B) is expected to be 1 iff A >= B, and 0 otherwise
%          ( default : domfun(A,B) = all(A*e>=B*e) with e = ones(n,1) )
%
%    [BOUNDS,P,INFO] = JSR_PROD_PRUNINGALGORITHM(M,OPTS)
%      Does the same but with the parameters specified in the structure
%      OPTS. See JSRSETTINGS and below for available parameters and options.
%
%      BOUNDS contains the lower and upper bounds on the JSR
%
%      P is a cell containing the indexes of two products attaining the
%        respective bounds (to be used with buildProduct)
%
%      INFO is a structure containing various data about the iterations :
%         info.status         - = 0 if normal termination, 1 if maxEval 
%                               reached, 2 if stopped in iteration                               
%                               (CTRL-C or opts.maxTime reached)
%         info.elapsedtime
%         info.niter          - max depth reached
%         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 bound at each depth 
%         info.allUb          - upper bound at each depth 
%         info.popHist        - number of candidates along depths
%         info.opts           - the structure of options used
%
%  The field opts.pruning (generated by jsrsettings) can be used to
%  tune the method :
%
%      pruning.maxEval        - maximum number of new products computed,
%                               (1000)
%      pruning.delta          - desired tolerance, stops if  
%                               (upper-lower<delta), (1e-8)
%      pruning.normfun        - function handle to a matrix norm, (@norm)
%                                  
%      pruning.domfun         - function handle to the domination function:
%                               domfun(A, B) is expected to be 1
%                                iff A >= B, and 0 otherwise,
%                               (@defaultdomfun : 
%                                defaultdomfun(A,B) = all(A*e>B*e)
%                                with e = ones(n,1))
%      pruning.plotBounds     - if 1 plots the evolution of the bounds w.r.t.
%                               the depth, (0)
%      pruning.plotpopHist    - if 1 plots the evolution of the number of
%                               candidates, (0)
%      pruning.plotTime       - if 1 plots evolution of the bounds found 
%                               w.r.t. time of computation, (0)
%
% REFERENCES
%    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
%    B.Moision, A.Orlitsky, and P.Siegel. 
%      "On Codes that Avoid Specified Differences",
%      IEEE Transactions on Information Theory 47, 2001
%     

if (nargin>1)
    if (length(varargin)==3 && isnumeric(varargin{1}) && isa(varargin{2},'function_handle') && isa(varargin{2},'function_handle') )
        opts = jsrsettings('pruning.maxEval',varargin{1},'pruning.normfun',varargin{2},'pruning.domfun',varargin{3});
    elseif (length(varargin)== 2 && isnumeric(varargin{1}) && isa(varargin{2},'function_handle'))
        opts = jsrsettings('pruning.maxEval',varargin{1},'pruning.normfun',varargin{2});
    elseif (length(varargin)== 1 && isnumeric(varargin{1}) )
        opts = jsrsettings('pruning.maxEval',varargin{1});
    elseif (length(varargin)== 1 && isstruct(varargin{1}) )
        opts = varargin{1};
    else
        warning('Invalid optional arguments, using default parameters !');
        opts = jsrsettings;
    end
else
    opts = jsrsettings;
end

close = 1;

% logfile opening
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_pruningAlg','wt');
        if (logFile == -1)
            warning('Could not open logfile')
        end
    else
        logFile = opts.logfile;
        close =0;
    end
else
    logFile = fopen('log_pruningAlg','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_pruningAlgorithm ******** \n \n')
starttime = cputime;

% Parameters
maxEval = opts.pruning.maxEval;
normfun  = opts.pruning.normfun;
delta = opts.pruning.delta;
if (isa(opts.pruning.domfun,'function_handle'))
    domfun = opts.pruning.domfun;
else
    domfun = @defaultdomfun;
end

% Initialization
scale = max(max(rho(M)),1e-6);
M = cellDivide(M,scale);
maxTime = 0;
status = 2;
m = length(M); % number of matrices
n = length(M{1}); % dimension of the matrices
D = reshape(eye(n),n*n,1)';
d = 1;
K = [];
k = (1:m)';
lower = scale;
upper = Inf;
P = {[], []};
partial = zeros(2, 200);
npophist = zeros(1,200);
timeIter = zeros(1,200);
eval = 0;
depth = 0;

% Iteration
while (eval < maxEval && (upper-lower)>delta)
    eval = eval+d*m;
    depth = depth +1;
    
    if(depth>length(npophist))
        partial = [partial zeros(2,200)];
        npophist = [npophist zeros(1,200)];
        timeIter = [timeIter zeros(1,200)];
    end
    
    % Build the new products
    newD = zeros(d*m, n^2);
    for i = 1:d,
        current = reshape(D(i, :), n, n);
        for j = 1:m,
            newD(m*i+j-m, :) = reshape(M{j}*current,n*n,1)';
        end
    end
    newK = zeros(d*m, depth);
    for i = 1:m,
        newK(i:m:end, 1:depth-1) = K;
    end
    newK(:, depth) = repmat(k, d, 1);
    
    % Remove dominated products
    [D, K] = prune(domfun, newD, newK, n);
    d = size(D, 1);
    npophist(depth) = d;
    holdupper = 0;
    holdupperarg = [];
    
    % Analyse dominating products
    for i = 1:d,
        current = reshape(D(i, :), n,n);
        
        sublower = (abs(max(eig(current))))^(1/depth);
      
        partial(1, depth) = max(partial(1, depth), scale*sublower);
        if (sublower > 1),
            lower = sublower*scale;
            scale=lower;
            rr=abs(max(eig(current)));
            D=D/rr;
            
            M = cellDivide(M,sublower);
            P{1} = K(i, :);
        end
        subupper = (feval(normfun,current))^(1/depth)*scale;
        partial(2, depth) = max(partial(2, depth), subupper);
        if (subupper > holdupper),
            holdupper = subupper;
            holdupperarg = K(i, :);
        end
    end
    if (holdupper < upper),
        upper = holdupper;
        P{2} = holdupperarg;
    end
    
    if (depth>1000)
        if (mod(depth,100)==0)
            msg(logFile,opts.verbose>0,'  Depth %3d  -  [%15.12g, %15.12g]  -  %d candidates', depth, lower, upper, d);
        end
    elseif (depth>150)
        if (mod(depth,50)==0)
            msg(logFile,opts.verbose>0,'  Depth %3d  -  [%15.12g, %15.12g]  -  %d candidates', depth, lower, upper, d);
        end
    elseif (depth>20)
        if (mod(depth,10)==0)
            msg(logFile,opts.verbose>0,'  Depth %3d  -  [%15.12g, %15.12g]  -  %d candidates', depth, lower, upper, d);
        end
    else
        msg(logFile,opts.verbose>0,'  Depth %3d  -  [%15.12g, %15.12g]  -  %d candidates', depth, lower, upper, d);
    end
    
    

    
    timeIter(depth) = cputime-starttime;
    
    % Save to file option
    if (ischar(opts.saveinIt))
        elapsedtime = cputime - starttime;
        bounds = [lower, upper];

        info.status = status;
        info.elapsedtime = elapsedtime;
        info.niter = depth;
        info.timeIter = timeIter(1:depth);
        info.allLb = (partial(1,1:depth));
        info.allUb = (partial(2,1:depth));
        info.popHist = npophist;
        info.opts  = opts;
        
        save([opts.saveinIt '.mat'],'bounds','P','info')
    end
    
    % Check projected time for next iteration
    if(depth>1)
        limProj = (cputime + timeIter(depth)-timeIter(depth) - starttime)>opts.maxTime;
    else
        limProj = 0;
    end
    
    if (timeIter(depth)>=opts.maxTime || limProj)
        maxTime = 1;
        msg(logFile,opts.verbose>0,'\nopts.maxTime reached\n');
        break;
    end
        
end

if (upper-lower<=delta)
    status = 0;
elseif (eval>= maxEval)
    status = 1;
elseif (maxTime)
    status = 2;
end

% Postprocessing
msg(logFile,opts.verbose>0,'\n> Bounds on the j.s.r.: [%.15g, %.15g]', lower, upper);
bounds = [lower, upper];
elapsedtime = cputime - starttime;

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

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

info.status = status;
info.elapsedtime = elapsedtime;
info.niter = depth;
info.timeIter = timeIter(1:depth);
info.allLb = (partial(1,1:depth));
info.allUb = (partial(2,1:depth));
info.popHist = npophist(1:depth);
info.opts  = opts;

if (ischar(opts.saveEnd))
    save([opts.saveEnd '.mat'],'bounds','P','info');
end

% Figures
if(opts.pruning.plotBounds)
    figure
    plot(1:info.niter,info.allLb,'-+g',1:info.niter,info.allUb,'-*r')
    title('pruningAlg : Bounds found at each depth')
    legend('Lower bound','Upper bound')
    xlabel('Depth')
end

if(opts.pruning.plotpopHist)
    figure
    plot(1:info.niter,info.popHist,'-o')
    title('pruningAlg : Evolution of number of candidates')
    xlabel('Depth')
end

if (opts.pruning.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('pruning : Bounds at each depth w.r.t. to time')
    legend('Lower Bound','Upper bound')   
end

end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [D, K] = prune(domfun, D, K, n)

[D, key] = unique(D, 'rows');
K = K(key, :);
m = size(D, 1);
i = 1;
while i < m,
    A = reshape(D(i, :),n,n);
    j = i+1;
    while j <= m,
        B = reshape(D(j, :), n,n);
        if feval(domfun,A, B),
            D = D([1:j-1, j+1:m], :);
            K = K([1:j-1, j+1:m], :);
            m = m-1;
            j = j-1;
        elseif feval(domfun,B, A),
            D = D([1:i-1, i+1:m], :);
            K = K([1:i-1, i+1:m], :);
            m = m-1;
            i = i-1;
            break;
        end
        j = j+1;
    end
    i = i+1;
end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function b = defaultdomfun(A, B)
n = length(A);
e = ones(n,1);

b = all(A*e>B*e);

% Chia-Tche's default dominating function :
% b = all(all(A'*A >= B'*B));
end
