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

% JSR_LIFT_SEMIDEFINITE Approximates the jsr using semidefinite liftings.
%
%    [BOUNDS, INFO] = JSR_LIFT_SEMIDEFINITE(M)
%      returns lower and upper bounds on the jsr of M using semidefinite
%      liftings (recursive procedure described in [1], theorem 5).
%      M must be a cell array containing the considered matrices.
%      Uses default values for the parameters (see below).
%
%    [BOUNDS, INFO] = JSR_LIFT_SEMIDEFINITE(M, MAXDEPTH)
%      where MAXDEPTH is an integer, goes only to MAXDEPTH
%      lifts.
%
%    [BOUNDS, INFO] = JSR_LIFT_SEMIDEFINITE(M, OPTS)
%      With OPTS an option structure created by jsrsettings uses the
%      parameters in OPTS. 
%
%  The function may abort the procedure before reaching MAXDEPTH
%  (opts.semidef.maxDepth) if there is not enough memory
%  (software of hardware limitation).
%
%  NOTE : if activated, the general option 'saveinIt' (see help jsrsettings) 
%         for saving at each iteration uses memory and
%         causes the procedure to abort earlier.         
%      
%      BOUNDS are the best attained bounds
%
%      INFO is a structure containing various data about the iterations :
%         info.status         - = 0 if normal termination, 1 if aborted
%                               due to hardware software limitations, 2
%                               if stopped in iteration (CTRL-C or
%                               opts.maxTime reached)
%         info.elapsedtime 
%         info.niter          - number of iterations (Depth reached)
%         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  [1xniter double]
%         info.allUb          - evolution of the upper bound  [1xniter double]
%         info.opts           - the structure of options used
%
%
%  The field options.semidef (generated by jsrSettings) can be used to
%  tune the method :
%
%      semidef.maxDepth       - max number of iterations (lifts), (6)
%
%      semidef.plotBounds     - if 1 plots the evolution of the bounds, (0)
%
%      semidef.plotTime       - if 1 plots evolution of the bounds w.r.t. time of
%                               computation, (0)
%
% See also JSRSETTINGS
%
% REFERENCES
%  [1]  V.D.Blondel and Y.Nesterov,
%         "Computationally efficient approximations of the joint spectral radius",
%         SIAM J. of Matrix Analysis, 27(1):256-272, 2005
%  [2]  V.D.Blondel, R. Jungers and V.Protasov,
%         "On the complexity of computing the capacity of codes that avoid forbidden 
%          difference patterns" 
%         IEEE Transactions on Information Theory, 52(11):5122-5127, 2006
%  [3]  R.Jungers, 
%         "The Joint Spectral Radius: Theory and Applications" 
%         Vol. 385 section 2.3.6 in Lecture Notes in Control and Information
%         Sciences, Springer-Verlag. Berlin Heidelberg, June 2009

if (nargin == 1)
    opts = jsrsettings;
end

if (nargin==2)
   if(isstruct(varargin{1}))
    opts = varargin{1};
   elseif(isnumeric(varargin{1}))
       opts = jsrsettings('semidef.maxDepth',varargin{1});
   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_semidefinite','wt');
        if (logFile == -1)
            warning('Could not open logfile')
        end
    else
        logFile = opts.logfile;
        close =0;
    end
else
    logFile = fopen('log_semidefinite','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_lift_semidefinite ******** \n \n')
starttime = cputime;

% Initialization
maxdepth = opts.semidef.maxDepth;
m = length(M);
n = size(M{1}, 1);
allBounds = zeros(maxdepth, 2);
allBounds(:, 2) = Inf;
timeIter = zeros(1,maxdepth);

% Iteration
try
    for i = 1:maxdepth,
       
        numEntries = length(M{1})^4; % number of entries in Mj below  % Correction

        numBytesNeeded = numEntries * 8; % supposes double (8 Bytes) 

        if numBytesNeeded >  0.85*available_memory % limit to 85% of physical memory
            error('JSR_LIFT_SEMIDEF:OUTOFMEMORY','too much memory would be needed')
        end
        
        msg(logFile,opts.verbose>1,...
            '\n Starting lift number %2.0f \n', i);
        MM = cell(m, 1);
        nn = n*(n+1)/2;
        S = zeros(nn, nn);
        keep = true(n^2, 1);
        
        % Build representative matrices
        for j = 1:m,
            Mj = kron(M{j}, M{j});
            for k = 1:n,
                keep(k*n-n+k+1:k*n) = 0;
                Mj(:, k*n+k:n:n*n-n+k) = Mj(:, k*n+k:n:n*n-n+k) + Mj(:, k*n-n+k+1:k*n);
            end
            MM{j} = Mj(keep, keep);
            S = S + MM{j};
        end
        allBounds(i, 2) = max(abs(eig(S)))^(1/2^i);
        allBounds(i, 1) = 1/(m^(1/2^i)) * allBounds(i, 2);
        
        msg(logFile,opts.verbose>0,...
            '\n > Iteration %d - current bounds: [%.15g, %.15g] \n', i, allBounds(i,1), allBounds(i,2));
        
        M = MM;
        n = nn;
        timeIter(i) = cputime-starttime;
        
        % Save to file option
        if (ischar(opts.saveinIt))
            status = 2;
            bounds = [max(allBounds(:, 1)), min(allBounds(:, 2))];
            elapsedtime = cputime - starttime;

            info.status = status;
            info.elapsedtime = elapsedtime;
            info.niter = i;
            info.timeIter = timeIter(1:i); 
            info.allLb = allBounds(1:i,1)';
            info.allUb = allBounds(1:i,2)';
            info.opts  = opts;
            
            save([opts.saveinIt '.mat'],'bounds','info') 
            clear status elapsedtime bounds info
        end
        
        if (timeIter(i))>= opts.maxTime
            msg(logFile,opts.verbose>0,'\nopts.maxTime reached');
            break;
        end
    end
catch ME
        msg(logFile,opts.verbose>0,'Aborting at depth %d...\n', i);
        i = i-1;
end

% Post-processing
if (i==maxdepth)
    status = 0;
elseif (timeIter(i) >= opts.maxTime)
    status = 2;
else
    status = 1;
end

bounds = [max(allBounds(:, 1)), min(allBounds(:, 2))];
elapsedtime = cputime - starttime;

msg(logFile,opts.verbose>0,'\n> Bounds on the jsr: [%.15g, %.15g]', bounds(1), bounds(2));

msg(logFile,opts.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 = i;
info.timeIter = timeIter(1:i);
info.allLb = allBounds(1:i,1)';
info.allUb = allBounds(1:i,2)';
info.opts  = opts; 

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

% Figures
if (opts.semidef.plotBounds)
    figure
    plot(1:i,info.allLb,'-+g',1:i,info.allUb,'-*r')
    title('semidefinite : Evolution of the bounds on the JSR')
    xlabel('Number of lifts')
    legend('Lower bound','Upper bound')
end

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

end
