function [bound, details, info] = jsr_prod_lowerBruteForce(M, varargin)

% JSR_PROD_LOWERBRUTEFORCE Gives a lower bound on the jsr using brute force.
%
%    [BOUND, DETAILS, INFO] = JSR_PROD_LOWERBRUTEFORCE(M)
%      returns a lower bound on the jsr of M using brute force by considering
%      products of length 1,3 and 5. M must be a cell array of matrices.
%
%    [BOUND, DETAILS, INFO] = JSR_PROD_LOWERBRUTEFORCE(M, DEPTHS, NDISP)
%      returns a lower bound on the jsr of M using brute force by considering
%      products of length DEPTHS. 
%      DEPTHS is a vector containing the lengths of products one wants to
%      consider. 
%      NDISP is optional, if specified displays the best NDISP products for each depth.
%            The 2*NDISP overall best products will be displayed at the end.
%
%    [BOUND, DETAILS, INFO] = JSR_PROD_LOWERBRUTEFORCE(M, OPTS)
%      does the same but with the values of the parameters specified in the
%      structure OPTS. See below and help JSRSETTINGS for available
%      parameters. 
%
%      BOUND is the best attained lower bound
%
%      DETAILS is a cell array containing the products and their averaged
%              spectral radii.
%
%      INFO is a structure containing various data about the iterations :
%         info.status         - 0 normal termination, 1 some depths could
%                               not be computed due to memory limitations,
%                               2 intermediate results 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 (computed) depth
%         info.allLbprod      - cell array with product indices
%                               corresponding to the lower bounds in        
%                               info.allLb
%         info.opts           - the structure of options used     
%
%  The field options.lowbrut (generated by jsrSettings) can be used to
%  tune the method :
%
%      lowbrut.depths         - vector containing the lengths of products 
%                               one wants to consider, ([1:2:5])
%      lowbrut.ndisp          - displays the best ndisp products for each depth.
%                               The 2*ndisp overall best products will be
%                               displayed at the end., (0)
%      lowbrut.plotBounds     - if 1 plots the lower bound found at
%                               each depth, (0)
%      lowbrut.plotTime       - if 1 plots evolution of the bound found 
%                               w.r.t. time of computation, (0)
%                               
% 
% See also JSRSETTINGS
% 
% REFERENCES
%    R.Jungers, 
%      "The Joint Spectral Radius: Theory and Applications" 
%      Vol. 385 in Lecture Notes in Control and Information
%      Sciences, Springer-Verlag. Berlin Heidelberg, June 2009

if (nargin>1)
    if (length(varargin)==2 && isvector(varargin{1}) && isscalar(varargin{2}) )
        opts = jsrsettings('lowbrut.depths',varargin{1},'lowbrut.ndisp',varargin{2});
    elseif (length(varargin)==1 && isnumeric(varargin{1}))
        opts = jsrsettings('lowbrut.depths',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_lowerBruteForce','wt');
        if (logFile == -1)
            warning('Could not open logfile')
        end
    else
        logFile = opts.logfile;
        close = 0;
    end
else
    logFile = fopen('log_lowerBruteForce','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_LowerBruteForce ******** \n \n')
starttime = cputime;

% Initialization
maxtime = 0;
status = 2;
depths = opts.lowbrut.depths;
ndisp = opts.lowbrut.ndisp;
bound = 0;
s = length(M);
n = length(depths);
details = cell(n, 2);
bestsval = zeros(n*ndisp, 1);
bestsarg = cell(n*ndisp, 1);
impDepths = zeros(1,n);
bestDepth = zeros(1,n);
bestDepthprod = cell(1,n);
timeIter = zeros(1,n);

% Scan through all products ignoring cyclic permutations
for i = 1:n,
    p = depths(i);
    try
        msg(logFile,opts.verbose>1,'\nGenerate all products of length %d',p);
        details{i, 1} = 1+genNecklaces(p, s); % Generates all products to test
        D = details{i, 1};
        m = size(D, 1);
        sr = zeros(m, 1);
        for j = 1:m,
            X = eye(size(M{1}, 1));
            for k = 1:p,
                X = X * M{ D(j, k) };
            end
            sr(j) = rho(X);
        end
        details{i, 2} = sr.^(1/p);
        [besti indbesti] = max(details{i,2});
        bestDepth(i) = besti;
        bestDepthprod{i} = details{i,1}(indbesti,:);
        
        msg(logFile,opts.verbose>1 && ~ndisp,'Best lower bound at length %d : %.15g \n', p,besti);
        
        if (besti>bound)
            bound = besti;
            bestProd = details{i,1}(indbesti,:);
        end

        % Display top <ndisp> products
        if ndisp,
            [sorted, idx] = sort(details{i, 2}, 1, 'descend');
            msg(logFile,1,'Best products of length %d:', p);
            for j = 1:min(ndisp, length(idx)),
                prodval = details{i, 2}(idx(j));
                prodarg = details{i, 1}(idx(j), :);
                msg(logFile,1,'  #%2d : bound of %15.9g with product %s', j, prodval, num2str(prodarg));
                bestsval((i-1)*ndisp + j) = prodval;
                bestsarg{(i-1)*ndisp + j} = prodarg;
            end
        end
        
    catch ME
        if(strcmp(ME.identifier,'MATLAB:pmaxsize') )
            msg(logFile,opts.verbose>0,'\n**********************\n Impossible for depth %d : \n %s \n**********************\n ',p,ME.message)
            impDepths(i) = 1;
            status = 1;
        else
            rethrow(ME);
        end
    end
    
    timeIter(i) = cputime-starttime;
    % Save to file option
    if (ischar(opts.saveinIt))
        elapsedtime = cputime - starttime;

        info.status = status
        info.elapsedtime = elapsedtime;
        info.timeIter = timeIter(1:i);
        info.allLb = bestDepth;
        info.allLbprod = bestDepthprod;
        info.opts = opts;

        save([opts.saveinIt '.mat'],'bound','details','info')
    end
    
    if (timeIter(i)>= opts.maxTime)
        maxtime = 1;
        msg(logFile,opts.verbose>0,'\nopts.maxTime reached\n');
        break;
    end
        
end

% Display overall top 2*<ndisp> products
if ndisp,
    msg(logFile,1,'\nOverall best products:');
    [sorted, idx] = sort(bestsval, 1, 'descend');
    bestsval = bestsval(idx);
    bestsarg = bestsarg(idx);
    
    % Remove equivalent products
    idx = 0;
    new = 1;
    for i = 1:n*ndisp,
        if isempty(bestsarg{i}),
            break;
        end
        if new || abs(bestsval(i) - bestsval(idx)) >= 1e-10,
            new = 0;
            idx = i;
        else
            r = length(bestsarg{i});
            update = 0;
            for j = idx:i-1,
                if isempty(bestsarg{j}),
                    continue;
                end
                if ~update,
                    update = 1;
                    idx = j;
                end
                q = length(bestsarg{j});
                s = lcm(q, r);
                oldtest = repmat(bestsarg{j}, 1, s/q);
                newtest = repmat(bestsarg{i}, 1, s/r);
                if all(oldtest == newtest),
                    if (q < r),
                        bestsval(i) = 0;
                        bestsarg{i} = [];
                        break;
                    else
                        bestsval(j) = 0;
                        bestsarg{j} = [];
                        idx = idx + 1;
                    end
                end
            end
        end
    end
    
    % Display best products
    [sorted, idx] = sort(bestsval, 1, 'descend');
    bestsval = bestsval(idx);
    bestsarg = bestsarg(idx);
    for j = 1:min(2*ndisp, length(idx)),
        if isempty(bestsarg{j}),
            break;
        end
        msg(logFile,1,'  #%2d : bound of %15.9g with product %s', j, bestsval(j), num2str(bestsarg{j}));
    end
end

if(sum(impDepths))
    status = 1;
elseif (maxtime)
    status = 2;
else
    status = 0;
end

% Post-processing
elapsedtime = cputime - starttime;

info.status = status;
info.elapsedtime = elapsedtime;
info.timeIter = timeIter(1:n);
info.allLb = bestDepth;
info.allLbprod = bestDepthprod;
info.opts = opts;

msg(logFile,opts.verbose>0,'\n> Best lower bound on the jsr: %.15g', bound);

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

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

if(ischar(opts.saveEnd))
    save([opts.saveEnd '.mat'],'bound','details','info')
end

% Figure
if(opts.lowbrut.plotBounds)
    figure
    plot(depths(~impDepths),bestDepth(~impDepths),'-+g','MarkerSize',5)
    title('lowerBruteForce : Best lower bounds found at each depth')
    xlabel('Depth')
end

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

end