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