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,'w');
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','w');
if (logFile == -1)
warning('Could not open logfile')
end
else
logFile = opts.logfile;
close = 0;
end
else
logFile = fopen('log_lowerBruteForce','w');
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
