function [bounds, info] = jsr_norm_maxRelaxation(M, varargin)

% JSR_NORM_MAXRELAXATION Approximates the jsr using Max-Relaxation.
%
%    [BOUNDS, INFO] = JSR_NORM_MAXRELAXATION(M, OPTS)
%      returns lower and upper bounds on the jsr of M using Max-Relaxation
%      scheme (heuristically).
%      M is a cell array with the considered matrices. Uses default values
%      (see below) for all parameters.
%
%    [BOUNDS, INFO] = JSR_NORM_MAXRELAXATION(M, N, AVGFUN)
%      The parameter N corresponds to the number of discretization points.
%      Default 1000.
%      The parameter AVGFUN is either a handler to an averaging function with
%      2 inputs and 1 output, or 'a', 'g', 'h', 'q' for arithmetic,
%      geometric, harmonic or quadratic mean. If not specified, default is 'a'.
%
%    [BOUNDS, INFO] = JSR_NORM_MAXRELAXATION(M, OPTS)
%      Does the same but with the values of parameters specified by the 
%      structure OPTS (see JSRSETTINGS and below for available parameters).
%      
%      BOUNDS contains the lower and upper bounds on the JSR
%
%      INFO is a structure containing various data about the iterations :
%         info.status         - = 0 if normal termination, 1 if maxiter
%                               reached, 2 if stopped in iteration
%                               (CTRL-C or opts.maxTime reached)                              
%         info.elapsedtime 
%         info.niter          - number of iterations
%         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          - evolution of the lower bound  [1x(niter) double]
%         info.allUb          - evolution of the upper bound  [1x(niter) double]
%         info.opts           - the structure of options used
%
%  The field opts.maxrel (generated by jsrsettings) can be used to
%  tune the method :
%
%      maxrel.N              - Number of discretization points. (10000)
%
%      maxrel.avgfun         - either a handler to an averaging function with
%                              2 inputs and 1 output, either 'a', 'g', 'h', 'q' 
%                              for arithmetic, geometric, harmonic or
%                              quadratic mean. ('a')
%
%      maxrel.maxiter        - maximum number of iterations, (500)
%
%      maxrel.tol            - tolerance for stopping condition, stops when
%                              upperBound-lowerBound < tol, (1e-8)
%
%      maxrel.plotBounds     - if 1 plots the evolution of the bounds, (0)
%
%      maxrel.plotTime       - if 1 plots evolution of the bounds w.r.t.
%                              time of computation, (0)
%
% See also JSRSETTINGS
%
% REMARK : This program returns heuristic estimates based on [1]
% 
% REFERENCES
% [1] V.Kozyakin, 
%       "Iterative building of Barabanov norms and computation of the
%        joint spectral radius for matrix sets",
%       arXiv:0810.2154v2 [math.RA], 22 Oct. 2008


if (nargin > 1)
    if (length(varargin) > 1)
        opts = jsrsettings('maxrel.N',varargin{1},'maxrel.avgfun',varargin{2});
    elseif isnumeric(varargin{1})
        opts = jsrsettings('maxrel.N',varargin{1});
    else
        opts = varargin{1};
    end
else
    opts = jsrsettings;
end

avgfun = opts.maxrel.avgfun;
if ischar(avgfun),
    switch (avgfun(1)),
        case 'g',  avgfun = @geomMean;
        case 'h',  avgfun = @harmMean;
        case 'q',  avgfun = @quadMean;
        otherwise, avgfun = @arithMean;
    end
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_maxRelaxation','wt');
        if (logFile == -1)
            warning('Could not open logfile')
        end
    else
        logFile = opts.logfile;
        close =0;
    end
else
    logFile = fopen('log_maxRelaxation','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_norm_maxRelaxation ******** \n \n')

starttime = cputime;

% Parameters
N = opts.maxrel.N;
tol = opts.maxrel.tol;
maxiter = opts.maxrel.maxiter;
m = length(M);
n = size(M{1}, 1);

% Initialization
status = 2;
lowerBounds = zeros(1,maxiter);
upperBounds = zeros(1,maxiter);
timeIter = zeros(1,maxiter);

% Discretization
[X, T] = genDiscretization(N, n);
R = ones(1, N);

Arad = zeros(m, N);
Asph = cell(m, 1);
msg(logFile,opts.verbose>1,'Starting to compute the images of the points');
for i = 1:m,
    Y = X*M{i}';
    Arad(i, :) = sqrt(sum(Y.^2, 2))';
    Asph{i} = Y ./ (Arad(i, :)' * ones(1, n));
end
idx = dsearchn(X, T, [cell2mat(Asph) ; 1, zeros(1, n-1)]);
Nearest = reshape(idx(1:end-1), N, m)';
NearestZero = idx(end);

% Iteration
msg(logFile,opts.verbose>1,'Starting iteration \n');
for iter = 1:maxiter,
    Rall = Arad .* R(Nearest);
    Rstar = max(Rall, [], 1);
    RHOup = max(Rstar./R);
    RHOdown = min(Rstar./R);
    invGamma = 1/avgfun(RHOup, RHOdown);
    R = max(R, Rstar * invGamma);
    R = R / R(NearestZero);
    lowerBounds(iter) = RHOdown;
    upperBounds(iter) = RHOup;
    timeIter(iter) = cputime-starttime;
    
    if (iter>1000)
        if (mod(iter,100)==0)
    msg(logFile,opts.verbose>0,'> Iteration %3.0f - current bounds: [%.15g, %.15g] \n', iter, RHOdown, RHOup);
        end
    elseif (iter>150)
        if (mod(iter,50)==0)
    msg(logFile,opts.verbose>0,'> Iteration %3.0f - current bounds: [%.15g, %.15g] \n', iter, RHOdown, RHOup);
        end
    elseif (iter>20)
        if (mod(iter,10)==0)
    msg(logFile,opts.verbose>0,'> Iteration %3.0f - current bounds: [%.15g, %.15g] \n', iter, RHOdown, RHOup);
        end
    else
        msg(logFile,opts.verbose>0,'> Iteration %3.0f - current bounds: [%.15g, %.15g] \n', iter, RHOdown, RHOup);
    end
    
    % Save to file option
    if (ischar(opts.saveinIt))
        if (iter==maxiter);status=1;end
        bounds = [RHOdown, RHOup];
        elapsedtime = cputime - starttime;

        info.status = status;
        info.elapsedtime = elapsedtime;
        info.niter = iter;
        info.timeIter = timeIter(1:iter);
        info.R = R;
        info.X = X;
        info.allLb = lowerBounds(1:iter);
        info.allUb = upperBounds(1:iter);
        info.opts  = opts;
        
        save([opts.saveinIt '.mat'],'bounds','info')
    end
    
    if (RHOup - RHOdown < tol),
        break;
    end
    
    if (timeIter(iter)>=opts.maxTime)
        msg(logFile,opts.verbose>0,'\nopts.maxTime reached\n');
        break;
    end
end

if(iter==maxiter)
    status = 1;
elseif timeIter(iter)>=opts.maxTime
    status = 2;
else
    status = 0;
end

% Post-processing
msg(logFile,opts.verbose>0,'\n> Bounds on the jsr: [%.15g, %.15g]', RHOdown, RHOup);
bounds = [RHOdown, RHOup];
elapsedtime = cputime - starttime;

msg(logFile,opts.verbose>1,'\nEnd of algorithm after %5.2f s.',elapsedtime);
if (logFile~=-1 && close)
    fclose(logFile);
end

info.status = status;
info.elapsedtime = elapsedtime;
info.niter = iter;
info.timeIter = timeIter(1:iter);
info.R = R;
info.X = X;
info.allLb = lowerBounds(1:iter);
info.allUb = upperBounds(1:iter);
info.opts  = opts;

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

% Figures
if(opts.maxrel.plotBounds)
    figure
    it = 1:iter;
    plot(it,info.allUb,'-*r',it,info.allLb,'-g+')
    title('maxRel : Evolution of the bounds on the JSR')
    legend('Upper bound','Lower bound')
    xlabel('Iterations')
end

if (opts.maxrel.plotTime)
    figure
    if (info.timeIter(end) > 600)
        time = info.timeIter/60;
        plot(time,info.allLb,'-+g',time,info.allUb,'-*r','MarkerSize',5);
        xlabel('time from start in min')
    else
        plot(info.timeIter,info.allLb,'-+g',info.timeIter,info.allUb,'-*r','MarkerSize',5)
        xlabel('time from start in s')
    end
    title('maxRel : Evolution of the bounds w.r.t. to time')
    legend('Upper bound','Lower bound')   
end
end

%%%%%%%%%%%%%%%%%%%%%%%%%
% AVERAGING FUNCTIONS   %
%%%%%%%%%%%%%%%%%%%%%%%%%

function z = arithMean(x, y)
z = (x+y)/2;
end

function z = geomMean(x, y)
z = sqrt(x*y);
end

function z = harmMean(x, y)
z = 2*x*y/(x+y);
end

function z = quadMean(x, y)
z = sqrt((x*x + y*y)/2);
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DISCRETIZATION FUNCTION   %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [X, T] = genDiscretization(N, n)

% http://sci.tech-archive.net/Archive/sci.math.num-analysis/2005-04/msg00203.html
% Knuth, _The Art of Computer Programming_, vol. 2, 3rd. ed., section 3.4.1.E.6. 
X = randn(N, n);
X = X./(sqrt(sum(X.^2, 2))*ones(1, n));
T = delaunayn(X);
end
