function [bounds, allBounds, details] = itMeth(M,opts)
%
% ITMETH(M) Computes all products, of matrices in M, of length 1 to a
% certain limit, possibly takes kroenecker power and iterates
% asked methods.
%
% [BOUNDS, ALLBOUNDS, DETAILS] = ITMETH(M)
% Computes products of length 1 to 4 and launches jsr_conic_ellipsoid
% on the sets of products of each length.
%
% [BOUNDS, ALLBOUNDS, DETAILS] = ITMETH(M,OPTS)
% Uses the parameters values and options specified in OPTS generated
% by jsrsettings.
% Computes all products of length 1 to max(OPTS.itMeth.length). For
% the length's contained in the vector OPTS.itMeth.length, launches the
% methods specified in OPTS.itMeth.meth and keeps the best bounds. Can
% also take the kroenecker powers of the matrices in the set.
% Stops when all the length's have been computed or when
% (upper-lower)<upper*OPTS.itMeth.tol.
% See below for default values of the fields in OPTS.
%
% BOUNDS is a vector with the best attained bounds
%
% AllBOUNDS contains the best bounds found at each length
%
% DETAILS is a structure with different fields containing information,
% (generated only if nargout>2);
%
% DETAILS.elapsedtime
%
% DETAILS.timeIter - time elapsed at the end of each iteration
%
% DETAILS.opts - the option structure used
%
% DETAILS.<method> - contains the bounds found at each length for
% the method <method>
%
% The field OPTS.itMeth (generated by jsrsettings) can be used to
% tune the method :
%
% itMeth.lift - = d, a positive integer, will take the d-th
% kroenecker power of each matrix at each
% length, (1)
% itMeth.length - if is a positive integer, itMeth will call
% the specified methods for products of
% length 1 to itMeth.length(end). If is a
% vector of increasing positive integers,
% launches only the methods on products of
% asked length, (4)
% itMeth.meth - cell array containing strings referring
% to the methods to use.
% Values can be :
% 'conitope' | 'ellips' | 'grip' |
% 'linear' | 'linrel' | 'maxrel'|'maxrel2D'
% 'pruning' | 'semidef'| 'sos'
% ( {'ellips'} )
% itMeth.tol - stopping condition, stops the method if
% upper-lower < OPTS.itMeth.tol*upper
% (1e-8)
%
% REMARKS : 1) There are m^k different products of length k
% and the d-th kroenecker power of an nxn matrix is
% (n^(d))x(n^(d)). Too high values of these parameters might
% cause Matlab to stop for memory issues.
% 2) The adequation between a method and the properties of the
% matrices is not checked by the algorithm, for instance no
% warning will be thrown if jsr_prod_pruningAlgorithm is
% called on matrices with negative entries.
%
% --------------------------------------------
% Example :
%
% opts = jsrsettings('itMeth.lift',2,'itMeth.length',[1, 5],...
% 'itMeth.meth',[{'ellips'} {'linear'}]);
%
% bounds = itMeth(M,opts)
%
% Will compute all products of length 1 (M), take their kroenecker
% power ( M2{i} = kron(M{i},M{i}) ), launch methods jsr_conic_ellipsoid and
% jsr_conic_linear on the resulting set and deduce bounds on the jsr of M.
% Then will repeat this operation on all products of length 5.
% --------------------------------------------
%
%
% See also jsrsettings
%
if nargin<2
opts = jsrsettings;
end
if (length(opts.itMeth.length)>1)
lengths = opts.itMeth.length;
maxLength = lengths(end);
nLengths = length(lengths);
else
maxLength = opts.itMeth.length;
lengths = 1:maxLength;
nLengths = maxLength;
end
% logfile opening
close =1;
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_itMeth','w');
if (logFile == -1)
warning('Could not open logfile')
end
else
logFile = opts.logfile;
close = 0;
end
else
logFile = fopen('log_itMeth','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 itMeth ******** \n \n')
starttime = cputime;
% Initialisation
d = opts.itMeth.lift;
meth = opts.itMeth.meth;
m = length(M);
n = length(M{1});
scale = ones(1,maxLength);
bounds = [0 Inf];
allBounds = [zeros(nLengths,1) Inf(nLengths,1)];
if (nargout>2)
for imeth =1:length(meth); eval(sprintf('bounds_%s = allBounds;',meth{imeth})); end
end
k = 0;
cont=1;
timeIt = zeros(1,nLengths);
Mk = M;
while (k < maxLength && cont)
k=k+1;
if k>1
Mksup = appMat(M,Mk);
scalek = max(rho(Mksup));
scale(k) = scale(k-1)*scalek;
else
Mksup = M;
scalek = 1;
end
% Prevents numeric explosion
Mksup = cellDivide(Mksup,scalek);
if (sum(lengths==k))
% Take d-th kroenecker power
if (d>1)
MksupLift = kronLift(Mksup,d);
else
MksupLift = Mksup;
end
msg(logFile,opts.verbose>0,'\n>>Starting methods for %d products of length %d, matrices are %dx%d\n\n',...
length(Mksup),k,size(MksupLift{1},1),size(MksupLift{1},2));
elapsedtime = cputime-starttime;
optsMeth = jsrsettings('verbose',opts.verbose,'logfile',logFile);
for imeth = 1:length(meth) % Executes each method
timeMeth = (opts.maxTime - elapsedtime)/((length(meth)-imeth+1)*nLengths);
optsMeth = jsrsettings(optsMeth,'maxTime',timeMeth);
switch lower(meth{imeth})
case 'conitope'
boundsk = jsr_norm_conitope(MksupLift,optsMeth);
case 'ellips'
boundsk = jsr_conic_ellipsoid(MksupLift,optsMeth);
case 'grip'
boundsk = jsr_prod_Gripenberg(MksupLift,optsMeth);
case 'linear'
boundsk = jsr_conic_linear(MksupLift,optsMeth);
case 'linrel'
boundsk = jsr_norm_linearRelaxation2D(MksupLift,optsMeth);
case 'maxrel'
boundsk = jsr_norm_maxRelaxation(MksupLift,optsMeth);
case 'maxrel2D'
boundsk = jsr_norm_maxRelaxation2D(MksupLift,optsMeth);
case 'pruning'
boundsk = jsr_prod_pruningAlgorithm(MksupLift,optsMeth);
case 'semidef'
boundsk = jsr_lift_semidefinite(MksupLift,optsMeth);
case 'sos'
boundsk = jsr_opti_sos(MksupLift,optsMeth);
otherwise
error(['Method %s in opts.itMeth.meth unknown'...
'\n Values can be :'...
'\n ''conitope'' | ''ellips''| ''Grip'' | ''linear'' '...
'\n ''linrel'' | ''maxrel'' |''maxrel2D'' '...
'\n ''pruning'' |''semidef''| ''sos'''],meth{imeth});
end
% Factor from kroenecker power
boundsk = boundsk.^(1/d);
% Factor from length and scaling
boundsk = (boundsk*scale(k)).^(1/k);
allBounds(lengths==k,1) = max(allBounds(lengths==k,1),boundsk(1));
allBounds(lengths==k,2) = min(allBounds(lengths==k,2),boundsk(2));
if (nargout>2)
eval(sprintf('bounds_%s(lengths==k,:) = boundsk;',meth{imeth}));
end
end
bounds(1) = max(allBounds(:,1));
bounds(2) = min(allBounds(:,2));
msg(logFile,opts.verbose>0,'\n>> ITMETH Length %d with kroenecker power %d terminated. Bounds : [%.15g, %.15g]',k,d,bounds(1),bounds(2))
if ((bounds(2)-bounds(1))<bounds(2)*opts.itMeth.tol)
cont = 0;
end
timeIt(lengths==k) = cputime-starttime;
if(timeIt>= opts.maxTime)
msg(-1,opts.verbose>0,'>>itMeth opts.maxTime reached')
cont = 0;
end
% Save to file option
if (ischar(opts.saveinIt))
if (nargout<3)
save([options.saveinIt '.mat'],'bounds','allBounds');
else
elapsedtime = cputime-starttime;
details.elapsedTime = elapsedtime;
details.timeIter = timeIt;
details.opts = opts;
for imeth = 1:length(meth)
eval(sprintf('details.%s = bounds_%s;',meth{imeth},meth{imeth}));
end
save([options.saveinIt '.mat'],'bounds','allBounds','details');
clear elapsedtime details
end
end
end
Mk = Mksup;
end
elapsedtime = cputime-starttime;
msg(logFile,opts.verbose>0,'\n>> ITMETH Bounds on the jsr : [%.15g, %.15g]',bounds(1),bounds(2));
msg(logFile,opts.verbose>1,'\n>> ITMETH End of algorithm after %5.2f s',elapsedtime);
if(logFile~=-1 && close)
fclose(logFile);
end
% Post-Processing
allBounds = allBounds(lengths<=k,:);
timeIt = timeIt(lengths<=k);
if nargout > 2
details.elapsedTime = elapsedtime;
details.timeIter = timeIt;
details.opts = opts;
for imeth = 1:length(meth)
eval(sprintf('details.%s = bounds_%s(lengths<=k,:);',meth{imeth},meth{imeth}));
end
end
% Save to file option
if (ischar(opts.saveEnd))
if (nargout<3)
save([options.saveEnd '.mat'],'bounds','allBounds');
else
save([options.saveEnd '.mat'],'bounds','allBounds','details');
end
end
% PlotBounds
if (opts.itMeth.plotBounds)
figure
plot(lengths(lengths<=k),allBounds(:,2),'-*r',lengths(lengths<=k),allBounds(:,1),'-g+')
title('itMeth : Evolution of the bounds on the JSR')
legend('Upper bound','Lower bound')
xlabel('Length of products')
end
% PlotTime
if (opts.itMeth.plotTime)
figure
if (timeIt(end) > 600)
time = timeIt/60;
plot(time,allBounds(:,1),'-+g',time,allBounds(:,2),'-*r','MarkerSize',10);
xlabel('time from start in min')
else
plot(timeIt,allBounds(:,1),'-+g',timeIt,allBounds(:,2),'-*r','MarkerSize',10);
xlabel('time from start in s')
end
title('itMeth : Evolution of the bounds w.r.t. time')
legend('Lower bound','Upper bound')
end
end
function Mksup = appMat(M,Mk)
m = length(M);
mk = length(Mk);
mksup = m*mk;
Mksup = cell(1,mksup);
for im=1:m
for imk = 1:mk
Mksup{mk*(im-1)+imk} = M{im}*Mk{imk};
end
end
end
function Mlift = kronLift(M,d)
m = length(M);
pow=1;
Mlift = M;
while(pow<d)
pow = pow+1;
for im = 1:m
Mlift{im} = kron(M{im},Mlift{im});
end
end
end
