function [bounds, V, info] = jsr_norm_balancedRealPolytope(M, C, k, myopts)

% JSR_NORM_BALANCEDREALPOLYTOPE Approximates the jsr using b.r.p.'s.
%
%    [BOUNDS,V, INFO] = JSR_NORM_BALANCEDREALPOLYTOPE(M, C, K)
%      returns lower and upper bounds on the jsr of M using balanced real
%      polytope norms.
%      The set M must be a cell array containing the considered matrices.
%      The matrix C is the candidate optimal product.
%      The parameter K corresponds to the length of the candidate product.
%
%    [BOUNDS,V, INFO] = JSR_NORM_BALANCEDREALPOLYTOPE(M, C, K, OPTS)
%      Uses the values of the parameters defined in the structure OPTS.
%      See below and JSRSETTINGS for the available parameters and options.
%
%      BOUNDS are the last lower and upper bounds found.
%
%      V is a matrix from wich the rows correspond to an essential 
%      system of vertices for the polytope.
%
%      INFO is a structure containing various data about the iterations :
%         info.status         - = 0 if normal termination, 1 if maxiter has
%                               been reached, 2 if stopped in iteration
%                               (CTRL-C or opts.maxTime reached)
%         info.elapsedtime 
%         info.niter          - number of iterations performed
%         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.Vhist          - Cell containing the sets of vectors at each
%                               iteration
%         info.popHist        - Number of vectors in essential set 
%                               at each iteration
%         info.allLb          - lower bounds, (constant = rho(C)^(1/K) )  [1xniter double]
%         info.allUb          - evolution of the upper bound              [1xniter double]
%         info.opts           - the structure of options used
%
%  The field opts.brp (generated by jsrsettings) can be used to
%  tune the method :
%
%      brp.maxiter        - max number of iterations, (1000)
%
%      brp.plotBounds     - if 1 plots the evolution of the bounds, (0)
%
%      brp.plotpopHist    - if 1 plots evolution of the size of the
%                           essential set, (0)
%      brp.plotTime       - if 1 plots evolution of the bounds w.r.t. time of
%                           computation, (0)
%
% See also JSRSETTINGS
%
% REFERENCES
%    A.Cicone, N.Guglielmi, S.Serra-Capizzano and M.Zennaro, 
%     "Finiteness property of pairs of 2x2 sign-matrices via real
%      extremal polytope norms",
%     Linear Algebra and its Applications, 432(2-3):796-816, 2010
%    V.Protasov,
%      "The generalized spectral radius. A geometric approach" ,
%      Izvestiya Mathematics, 61(5):995-1030, 1997


if (nargin < 4)
    myopts = jsrsettings;
    if (nargin <3)
        error('This method requires a candidate optimal product and its length')
    end
end

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

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

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

% Initialization
status = 2;
optseigs.disp = 0;
options = optimset('Display', 'off', 'TolFun', 1e-12);
m = length(M);
N = size(M{1}, 1);
F = M;
[v, d] = eigs(C, 1, 'LM', optseigs);
d = abs(d);
if(d < 1e-4)
    if(~jsr_zeroJsr(M))
        error('The spectral radius of the products is equal to zero, but the output of jsr_zeroJsr says that the JSR of the set is larger than zero. Please give another product with a positive spectral radius.');
    else
        msg(logFile,myopts.verbose>0,'The spectral radius has been proven to be equal to zero (see jsr_zeroJsr).\n The method will stop and returns [0,0] as lower and upper bounds.')
        bounds = [0,0];
        V = [];
        info.status = 0;
        info.elapsedtime = 0;
        info.niter = 0;
        info.timeIter = 0;
        info.Vhist = [];
        info.popHist = 0;
        info.allLb = 0;
        info.allUb = 0;
        info.opts = myopts;
        return;
    end
end
alpha = abs(d)^(1/k);
F = cellDivide(F,alpha);
V = v';
X = v';
w = size(X, 1);
maxiter = myopts.brp.maxiter;
mumin = Inf;
Vlog = cell(maxiter, 1);
nv = zeros(1,maxiter);
upperBounds = zeros(1,maxiter);
timeIter = zeros(1,maxiter);

% Iteration
for iter = 1:maxiter,
    Xsize = 4;
    Xused = 0;
    Xnew = zeros(Xsize, N);
    mu = 1;
    for i = 1:m,
        Y = F{i}*X';
        for j = 1:w,
            objective = ones(2*w, 1);
            constraintAeq = zeros(N, 2*w);
            constraintAeq(:, 1:2:2*w-1) = X';
            constraintAeq(:, 2:2:2*w) = -X';
            constraintBeq = Y(:, j);
            [dummy, mutemp, exitflag] = linprog(objective, [], [], constraintAeq, constraintBeq, zeros(2*w, 1), Inf*ones(2*w, 1), [], options);
            if (exitflag == -2) || (exitflag == 0) || isempty(mutemp),
                mutemp = Inf;
            end
            mu = max(mu, mutemp);
            if (mutemp > 1),
                if (Xused == Xsize),
                    Xsize = Xsize*2;
                    Xnew(Xsize, N) = 0;
                end
                Xused = Xused + 1;
                Xnew(Xused, :) = Y(:, j)';
            end
        end
    end
    mumin = min(mumin, mu);
    upperBounds(iter) = alpha*mumin;
    
    if (iter>1000)
        if (mod(iter,100)==0)
            msg(logFile,myopts.verbose>0,'> Iteration %d - current bounds: [%.15g, %.15g]', iter, alpha, alpha*mumin);
        end
    elseif (iter>150)
        if (mod(iter,50)==0)
            msg(logFile,myopts.verbose>0,'> Iteration %d - current bounds: [%.15g, %.15g]', iter, alpha, alpha*mumin);
        end
    elseif (iter>20)
        if (mod(iter,10)==0)
            msg(logFile,myopts.verbose>0,'> Iteration %d - current bounds: [%.15g, %.15g]', iter, alpha, alpha*mumin);
        end
    else
        msg(logFile,myopts.verbose>0,'> Iteration %d - current bounds: [%.15g, %.15g]', iter, alpha, alpha*mumin);
    end
            
    Xnew = Xnew(1:Xused, :);
    if isempty(Xnew),
        timeIter(iter) = cputime-starttime;
        break;
    end
    X = unique(Xnew, 'rows');
    V = unique([V ; -V ; X], 'rows');
    if (size(V, 1) > size(V, 2))
        W = unique(convhulln(V, {'QJ', 'Pp'}));
        V = V(W, :);
    end
    Vlog{iter} = V;
    nv(iter) = size(V,1);
    X = intersect(X, V, 'rows');
    w = size(X, 1);
    
    timeIter(iter) = cputime-starttime;
    
    % Save to file option
    if (ischar(myopts.saveinIt))
        if (iter==maxiter);status=1;end
        bounds = [alpha, alpha*mumin];
        elapsedtime = cputime - starttime;

        info.status = status;
        info.elapsedtime = elapsedtime;
        info.niter = iter;
        info.timeIter = timeIter(1:iter);
        info.Vhist = Vlog(1:iter);
        info.popHist = nv(1:iter);
        info.allLb = ones(1,iter)*alpha;
        info.allUb = upperBounds(1:iter);
        info.opts = myopts;

        save([myopts.saveinIt '.mat'],'bounds','V','info')
    end
    
    if (timeIter(iter)>=myopts.maxTime)
        msg(logFile,myopts.verbose>0,'\nopts.maxTime reached\n');
        break;
    end
end

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

% Post-processing
msg(logFile,myopts.verbose>0,'\n> Bounds on the jsr: [%.15g, %.15g]', alpha, alpha*mumin);
Vlog = Vlog(1:indLastV);
nv = nv(1:indLastV);
bounds = [alpha, alpha*mumin];
elapsedtime = cputime - starttime;

msg(logFile,myopts.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.niter = iter;
info.timeIter = timeIter(1:iter);
info.Vhist = Vlog(1:indLastV);
info.popHist = nv;
info.allLb = ones(1,iter)*alpha;
info.allUb = upperBounds(1:iter);
info.opts = myopts;

% Save output option
if (ischar(myopts.saveEnd))
   save([myopts.saveEnd '.mat'],'bounds','V','info')
end

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

if (myopts.brp.plotpopHist)
   figure
   plot(1:indLastV,nv,'-*');
   title('brp : Evolution of the size of the essential system');
   xlabel('Iteration')
   xlim([1 indLastV])
end

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

end