function [bounds, P, info] = jsr_conic_ellipsoid(M,varargin)

% JSR_CONIC_ELLIPSOID Approximates the jsr using ellipsoidal norms.
%
%
%    [BOUNDS, P, INFO] = JSR_CONIC_ELLIPSOID(M)
%      returns lower and upper bounds on the jsr of M using ellipsoid
%      norms. The set of matrices M must be a cell array of matrices.
%
%    [BOUNDS, P, INFO] = JSR_CONIC_ELLIPSOID(M, LOWER, UPPER) 
%      Uses the parameters UPPER and LOWER as starting bounds for the
%      bisection.
%
%    [BOUNDS, P, INFO] = JSR_CONIC_ELLIPSOID(M, OPTS)
%      Does the same but with value of the parameters defined in the
%      structure OPTS. See below and help JSRSETTINGS for available     
%      parameters. 
%
%      BOUNDS contains the lower and upper bounds on the JSR
%
%      P is the matrix defining the optimal norm
%
%      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 (ctr-c
%                               or maxTime reached)
%         info.elapsedtime 
%         info.bissec         - [lower upper] last bisection interval  
%         info.niter          - number of iterations (bisection)
%         info.opts           - the structure of options used
%
%
%  The field opts.ellips (generated by jsrsettings) can be used to
%  tune the method :
%
%      ellips.initBounds       - = [LOWER UPPER], gives starting bounds for
%                                bisection method, if not specified uses a
%                                few iterations of jsr_prod_bruteForce
%      ellips.maxiter          - maximum number of iterations (for bisection),
%                                (1000)
%      ellips.tol              - asked precision for bisection. Stopping
%                                condition is 
%                                upper-lower < options.ellips.tol*upper
%                                Influences the asked SeDuMi precision; 
%                                params.eps = options.ellips.tol/10, (1e-8)
%
% See also JSRSETTINGS, SOLVE_SEMI_DEFINITE_PROGRAM
%
%   NOTE :  This function uses the interface solve_semi_definite_program 
%           in order to call a LMI solver (default : SeDuMi). 
%           See help solve_semi_definite_program for more information.
%
% REFERENCES
%    V.D.Blondel, Y.Nesterov, J.Theys, 
%      "On the accuracy of the ellipsoid norm approximation
%       of the joint spectral radius",
%      Linear Algebra and its Applications, 394(1):91-107, 2005
%    T.Ando, M.-h. Shih, 
%      "Simultaneous contractibility", 
%      SIAM J. Matrix Anal. Appl. 19(2):487-498, 1998

if (nargin > 1)
    if (length(varargin) == 2)
        options = jsrsettings('ellips.initBounds',[varargin{1} varargin{2}]);
    elseif isstruct(varargin{1})
        options = varargin{1};
    else
        warning('Invalid optional arguments, using default parameters !');
        options = jsrsettings;
    end
else
    options = jsrsettings;
end

% logfile opening
close =1;
if (ischar(options.logfile) )    
    logFile = fopen(options.logfile,'wt');
    if (logFile == -1)
        warning(sprintf('Could not open file %s',options.logfile));
    end
elseif isnumeric(options.logfile)
    if (options.logfile==0)
        logFile = -1;
    elseif options.logfile==1
        logFile = fopen('log_ellipsoid','wt');
        if (logFile == -1)
            warning('Could not open logfile')
        end
    else
        logFile = options.logfile;
        close =0;
    end
else
    logFile = fopen('log_ellipsoid','wt');
    if (logFile == -1)
        warning('Could not open logfile')
    end
end

if (logFile~=-1)
    fprintf(logFile,[datestr(now) '\n\n']);
end

msg(logFile,options.verbose>1,'\n \n******** Starting jsr_conic_ellipsoid ******** \n \n')
starttime = cputime;

% Parameters
maxiter = options.ellips.maxiter;

% Initialization
m = length(M); % number of matrices
n = length(M{1}); % dimension of the matrices
opts.disp = 0;
iter = 0;
status = 2;

% Starting bounds
if(isempty(options.ellips.initBounds))
    maxdepth = ceil(log(42)/log(m));
    optsbrute = jsrsettings('verbose',0,'logfile',0,'bruteforce.maxdepth',maxdepth);
    bounds = jsr_prod_bruteForce(M, optsbrute);
    msg(logFile,options.verbose>1,'Bounds from bruteforce with depth %2.0f : [%.15g, %.15g]',maxdepth,bounds(1),bounds(2));
    upper = bounds(2);
    lower = bounds(1);
else
    lower = options.ellips.initBounds(1);
    upper = options.ellips.initBounds(2);
end


if (upper-lower <= options.ellips.tol*upper),
    P = [];
 
    bounds = [lower upper];
    msg(logFile,options.verbose>0,'No bisection needed, initial bounds are tight enough');
    
    elapsedtime = cputime - starttime;
    info.status = 0;
    info.elapsedtime = elapsedtime;
    info.niter = iter;
    info.opts = options;
    return;
end
candidate = (upper+lower) / 2;

% Semidefinite variables
K.s = ones(1, m+1) * n;
NB_VAR = (m+1) * n^2;

% Objective function
c = sparse(NB_VAR, 1);

% Constraints
A = sparse(m*n*(n+1)/2, NB_VAR);
b = sparse(m*n*(n+1)/2, 1);

idx = 0;
for i = 1:m,
    for j = 1:n,
        for k = j:n,
            idx = idx + 1;
            A(idx, i*n^2 + (j-1)*n + k) = -1;
            A(idx, 1:n^2) = - vec(M{i}(:, k) * M{i}(:, j)')';
            A(idx, j*n-n+k) = A(idx, j*n-n+k) + candidate^2;
        end
    end
end

stop = 0;

% bisection
msg(logFile,options.verbose>1,'\nStarting bisection method');
while (stop == 0 && iter < maxiter)
    iter = iter+1;
    
    valid = 0;
    
    if (iter>1000)
        if (mod(iter,100)==0)
            msg(logFile,options.verbose>0,'\n> Current upper-bound on JSR: %.15g (testing %.15g)',upper, lower, upper, candidate);
        end
    elseif (iter>150)
        if (mod(iter,50)==0)
            msg(logFile,options.verbose>0,'\n> Current upper-bound on JSR: %.15g (testing %.15g)',upper, lower, upper, candidate);
        end
    elseif (iter>20)
        if (mod(iter,10)==0)
            msg(logFile,options.verbose>0,'\n> Current upper-bound on JSR: %.15g (testing %.15g)',upper, lower, upper, candidate);
        end
    else
        msg(logFile,options.verbose>0,'\n> Current upper-bound on JSR: %.15g (testing %.15g)',upper, lower, upper, candidate);
    end
    
    [result,dual,feas,stopFlag] = solve_semi_definite_program(A,b,c,K,options.ellips.tol/10,options);
    
    % Check feasibility
    P = mat(result(1:n^2), n);
    for i = 1:m,
        if (min(eig(candidate^2 * P - M{i}' * P * M{i})) < 0),
            valid = -1;
            break;
        end
    end
    
    if (valid == -1),
        msg(logFile,options.verbose>1,'  Failure - Infeasible solution');
        valid = 0;
    elseif (feas == 0),
        msg(logFile,options.verbose>1,'  Failure - Infeasible problem');
    elseif (stopFlag ~= 0),
        msg(logFile,options.verbose>1,'  Failure - Numerical problems');
    else
        msg(logFile,options.verbose>1,'  Success - Found feasible solution');
        valid = 1;
    end

    % Update the bounds
    old = candidate;
    if valid,
        upper = candidate;
    else
        lower = candidate;
    end

    % Update the constraints
    if (upper-lower < upper*options.ellips.tol),
        candidate = upper;
        
        % End of bisection loop
        stop = 1;
    else
        candidate = (upper+lower) / 2;
    end
    
    idx = 0;
    for i = 1:m,
        for j = 1:n,
            for k = j:n,
                idx = idx + 1;
                A(idx, j*n-n+k) = A(idx, j*n-n+k) - old^2 + candidate^2;
            end
        end
    end

    % Save to file option
    if (ischar(options.saveinIt))
        if(iter==maxiter); status = 1; end;
        elapsedtime = cputime-starttime;
        
        info.status = status;
        info.elapsedtime = elapsedtime;
        info.bissec = [lower upper];
        info.niter = iter;
        info.opts = options;
        
        tempLowerEllipsEst = max(1/sqrt(n)*upper,1/sqrt(m)*upper);
        bounds = [tempLowerEllipsEst upper];
        
        save([options.saveinIt '.mat'],'bounds','info'); 
    end
    
    if ( (cputime-starttime)>=options.maxTime ) 
        stop = 1;
        msg(logFile,options.verbose>0,'\nopts.maxTime reached')
    end
    
end

if (iter==maxiter)
    status = 1;
elseif (cputime-starttime)>=options.maxTime
    status = 2;
else
    status = 0;
end

% Post-processing
result = solve_semi_definite_program(A,b,c,K,options.ellips.tol/10,options);
P = mat(result(1:n^2), n);
lowerEllipsEst = max(1/sqrt(n)*upper,1/sqrt(m)*upper);
msg(logFile,options.verbose>0,'\n> Bounds on the jsr: [%.15g, %.15g]', lowerEllipsEst, upper);
bounds = [lowerEllipsEst, upper];
elapsedtime = cputime - starttime;

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

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

info.status = status;
info.elapsedtime = elapsedtime;
info.bissec = [lower, upper];
info.niter = iter;
info.opts = options;