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

% JSR_NORM_BALANCEDCOMPLEXPOLYTOPE Approximates the jsr using b.c.p.'s.
%
%    [BOUNDS, V, INFO] = JSR_NORM_BALANCEDCOMPLEXPOLYTOPE(M, C, K)
%      returns lower and upper bounds on the jsr of M using balanced complex
%      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_BALANCEDCOMPLEXPOLYTOPE(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 best 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.bcp (generated by jsrsettings) can be used to
%  tune the method :
%
%      bcp.maxiter        - max number of iterations, (500)
%
%      bcp.plotBounds     - if 1 plots the evolution of the bounds, (0)
%
%      bcp.plotpopHist    - if 1 plots evolution of the size of the
%                           essential set, (0)
%      bcp.plotTime       - if 1 plots evolution of the bounds w.r.t. time 
%                           of computation, (0)
%
% See also JSRSETTINGS
%
% This is version 1.1 of JSR_NORM_BALANCEDCOMPLEXPOLYTOPE.
% The update fixes a bug where the function returns sometimes 
% an incorrect upper-bound when operating on a non-optimal cycle.
%
% REFERENCES
%    N.Guglielmi and M.Zennaro, 
%      "Finding extremal complex polytope norms for
%       families of real matrices",
%      SIAM J. Matrix Anal. Appl. 31(2):602-620, 2009
%    V.Protasov,
%      "The generalized spectral radius. A geometric approach",
%      Izvestiya Mathematics, 61(5):995-1030, 1997


tolJSR = 1e-6;

if (nargin<4)
    myopts = jsrsettings;
    if (nargin<3)
        error('This method requires a candidate optimal product and its length as input.')
    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_balancedComplexPolytope','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_balancedComplexPolytope ******** \n \n')
starttime = cputime;

% Initialization
opts.disp = 0;
m = length(M);
N = size(M{1}, 1);
[v, d] = eigs(C, 1, 'LM', opts);
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))*(1+tolJSR);
F = cellDivide(M,alpha);
maxiter = myopts.bcp.maxiter;
mumin = Inf;
Vlog = cell(maxiter, 1);
allUb = zeros(1,maxiter);
nv = zeros(1,maxiter);
status = 2;
timeIter = zeros(1,maxiter);

% Preprocessing: create an orthonormal basis using propagation of V0
oldRank = 0;
newRank = 1;
V = [v, zeros(N,N-1)];
while(newRank - oldRank > 0 )
    oldRank = newRank;
    for i = 1:m
        for r =1:newRank
        
            V(:,newRank+1) = M{i}*V(:,r);

            if(rank(V) > newRank)
                newRank = newRank+1;
            end

            if newRank >= N
                break;
            end
        end
    end
end

if(rank(V) < N)
    
    msg(logFile,1,'\n ***********************\n Matrices appear to be jointly block-triangularizable ! \n ***********************\n');
    msg(logFile,1,'\nIn this case the Joint Spectral Radius of the set is equal to max JSR of each set of blocks.\n');
    msg(logFile,1,'The function try to find triangularization with jointTriangul, will terminate and return the sets found as output.\n')
    msg(logFile,1,'\nPlease use quickElim and launch this and/or other method(s) on the relevant sets. \n');
    msg(logFile,1,'\nSee also jointTriangul and quickElim \n')
    msg(logFile,1,'\nREFERENCES')
    msg(logFile,1,'   R.Jungers,')
    msg(logFile,1,'   "The Joint Spectral Radius: Theory and Applications"')
    msg(logFile,1,'   Vol. 385 section 1.2.2.5 in Lecture Notes in Control and Information')
    msg(logFile,1,'   Sciences, Springer-Verlag. Berlin Heidelberg, June 2009')
    
    [TRIANG, BLOCKS, B] = jointTriangul(M);
    
    bounds.blocks = BLOCKS;
    bounds.B = B;
    
    info = 'Matrices appear to be jointly block triangularizable, see printed message.';
    
    return;
end

V = V.';
X = V;
w = size(X, 1);

% Iteration
for iter = 1:maxiter
    Xsize = 4;
    Xused = 0;
    Xnew = zeros(Xsize, N);
    mu = 1;
    
    sizeV = size(V, 1);
    A = zeros(2*N,3*w);
    for i = 1:sizeV
        A(: ,3*i-1) = [ real(V(i,:))' ; imag(V(i,:))' ];
        A(: ,3*i  ) = [-imag(V(i,:))' ; real(V(i,:))' ];
    end
    
    c = repmat([1 0 0],1,sizeV);
    K.q = ones(sizeV,1)*3;
    for i = 1:m,
        Y = F{i}*X.';
        for j = 1:w,
            
            b = [real(Y(:,j)) ; imag(Y(:,j))];
            [trash1,trash2,feas,stopFlagSDP,mutemp] = solve_semi_definite_program(A,b,c,K,1e-8);
                
            if stopFlagSDP == 2 || feas == false,
                mutemp = Inf;
            end
            
            mu = max(mu, mutemp);
            
            if (mutemp > 1)
                % not feasible
                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);
    
    if (iter>1000)
        if (mod(iter,100)==0)
            msg(logFile,myopts.verbose>0,'> Iteration %d - current bounds: [%.15g, %.15g]', iter, alpha, alpha*mumin*(1+tolJSR));
        end
    elseif (iter>150)
        if (mod(iter,50)==0)
            msg(logFile,myopts.verbose>0,'> Iteration %d - current bounds: [%.15g, %.15g]', iter, alpha, alpha*mumin*(1+tolJSR));
        end
    elseif (iter>20)
        if (mod(iter,10)==0)
            msg(logFile,myopts.verbose>0,'> Iteration %d - current bounds: [%.15g, %.15g]', iter, alpha, alpha*mumin*(1+tolJSR));
        end
    else
        msg(logFile,myopts.verbose>0,'> Iteration %d - current bounds: [%.15g, %.15g]', iter, alpha, alpha*mumin*(1+tolJSR));
    end
    
    allUb(iter) = alpha*mumin*(1+tolJSR);
    Xnew = Xnew(1:Xused, :);
    if isempty(Xnew)
        msg(logFile,myopts.verbose>0,'> Polyedra is proven to be invariant, up to precision %f.\n', tolJSR);
        timeIter(iter) = cputime-starttime;
        break;

    end
    X = unique(Xnew, 'rows');
    V = unique([V ; X], 'rows');
    if (size(V, 1) > size(V, 2))
        % Purge points in the interior
        i = 1;
        w = size(V, 1)-1;
        while (i <= w+1),
             A = zeros(2*size(V,2)+1,3*size(V,1));
            for k = 1:size(V,1)
                A(1     ,3*k-2) = 1;
                A(2:end ,3*k-1) = [ real(V(k,:))' ; imag(V(k,:))' ];
                A(2:end ,3*k  ) = [-imag(V(k,:))' ; real(V(k,:))' ];
            end
            A(:,(3*i-2):(3*i)) = [];
            b = [1 ; real(V(i,:))' ; imag(V(i,:))'];
            c = zeros(3*size(V,1)-3,1);
            K.q = ones(size(V,1)-1,1)*3;
            
            [trash1,trash2,feas,stopFlagSDP,mutemp] = solve_semi_definite_program(A,b,c,K,1e-8); %#ok<*ASGLU>
            
            if stopFlagSDP == 2 || feas == false
                mutemp = Inf;
            end
            if (mutemp <= 1),
                V(i, :) = [];
                w = w-1;
            else
                i = i+1;
            end
        end
    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*(1+tolJSR)];
        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 = allUb(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 = iter;
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*(1+tolJSR));
bounds = [alpha, alpha*mumin*(1+tolJSR)];
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(1:indLastV);
info.allLb = ones(1,iter)*alpha;
info.allUb = allUb(1:iter);
info.opts = myopts;

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

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

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

if (myopts.bcp.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('BCP : Evolution of the bounds w.r.t. to time')
    legend('Lower bound','Upper bound')   
end

end
