function [bounds, better, x, info] = jsr_conic_linear(M,varargin)

% JSR_CONIC_LINEAR Approximates the jsr of a set of nonnegative matrices 
%                  using the joint conic radius in the positive orthant.
%
%    [BOUNDS, BETTER, X, INFO] = JSR_CONIC_LINEAR(M)
%      returns lower and upper bounds on the jsr of of nonnegative matrices 
%      using the joint conic radius in the positive orthant.
%
%    [BOUNDS, BETTER, X, INFO] = JSR_CONIC_LINEAR(M, LOWER, UPPER) 
%      Uses parameters UPPER and LOWER are starting bounds for bisection
%      method.
%
%    [BOUNDS, BETTER, X, INFO] = JSR_CONIC_LINEAR(M, OPTS)
%      Does the same but with values 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 
%
%      BETTER is 1 iif one could better the initial upper bound. 
%
%      X is the last vector in the cone which growth could be bounded
%
%      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.bissec         - [lower upper] last bisection interval  
%         info.niter          - number of iterations (bisection)
%         info.opts           - the structure of options used
%
%
%  The field opts.linear (generated by jsrsettings) can be used to
%  tune the method :
%
%      linear.initBounds       - = [LOWER UPPER], gives starting bounds for
%                                bisection method, if not specified uses a
%                                few iterations of jsr_prod_bruteForce and 
%                                starts with [lowerBrute, n*upperBrute], ([])                               
%      linear.maxiter          - maximum number of iterations (for bisection),
%                                (1000)
%      linear.tol              - asked precision for bisection. Stopping
%                                condition is 
%                                upper-lower < options.linear.tol*upper
%                                Influences the asked SeDuMi precision; 
%                                params.eps = options.linear.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.Protasov, R.Jungers, V.D.Blondel, 
%      "Joint spectral characteristics of matrices: 
%       a conic programming approach",
%      SIAM J. Matrix Anal. Appl. 31(4):2146-2162, 2010
%
%

if (nargin > 1)
    if (length(varargin) == 2)
        options = jsrsettings('linear.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_linear','wt');
        if (logFile == -1)
            warning('Could not open logfile')
        end
    else
        logFile = options.logfile;
        close = 0;
    end
else
    logFile = fopen('log_linear','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_linear ******** \n \n')
starttime = cputime;

% Parameters
maxiter = options.linear.maxiter;

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

% Starting bounds
if(isempty(options.linear.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 = n*bounds(2);
lower = bounds(1);

else
    lower = options.linear.initBounds(1);
    upper = options.linear.initBounds(2);
end


if (upper-lower <= options.linear.tol*upper),
    x = [];
    
    lowerLinearEst = 1/n*upper;
    bounds = [lowerLinearEst 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;

% Number of nonnegative variables
NB_VAR = n+m*n;
K.l = NB_VAR;

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

% Constraints
A = sparse(m*n+1, NB_VAR);
b = sparse(m*n+1, 1);
b(1,1) = 1;
In1 = ones(n,1);
Inn = spdiags(In1,0,n,n);
Inm1 = ones(m*n,1);
Inmnm = spdiags(Inm1,0,n*m,n*m);
Imult = kron(ones(m,1),Inn);
Mmult = sparse(m*n,n);

for i = 1:m,
    Mmult( (i-1)*n+(1:n),1:n) = M{i};
end

A(1,1:n) = sparse(ones(1,n));
A(2:end,1:n) = sparse(candidate*Imult-Mmult);
A(2:end,n+1:end) = -Inmnm;

options.SDPsolver.solverOptions.fid = 0;
options.SDPsolver.solverOptions.maxiter = 100;
tolSolver = options.linear.tol/10;
options.SDPsolver.solverOptions.bigeps = 1e-3;
stop = 0;

% bisection
msg(logFile,options.verbose>0,'\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,tolSolver,options);

    % Check feasibility
    x = result(1:n);
    for i = 1:m,
        if (sum(sum( (candidate*Inn*x-M{i}*x)<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
    if valid,
        upper = candidate;
        better = 1;
    else
        lower = candidate;
    end

    % Check ending condition
    if (upper-lower < upper*options.linear.tol),
        candidate = upper;
        
        % End of bisection loop
        stop = 1;
    else
        candidate = (upper+lower) / 2;
    end
    
    % Actualise constraints
    A(2:end,1:n) = sparse(candidate*Imult-Mmult);

    % 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;
        
        bounds = [1/n*upper 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 = 0;
end

% Post-processing
result = solve_semi_definite_program(A,b,c,K,tolSolver,options);
x = result(1:n);
msg(logFile,options.verbose>0,'\n> Bounds on the jsr: [%.15g, %.15g]', 1/n*upper, upper);
bounds = [1/n*upper, 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;