function  [Ub, P, info] = jsr_pathcomplete(M,varargin)
% JSR_PATHCOMPLETE  Computes an upper bound for the JSR with (improved) 
%                   bisection method and provides also multiple quadratic 
%                   Lyapunov  functions for the set of matrices M with the 
%                   help of a path-complete automaton.
%
% #########################################################################
% #  This is a basic help. Type help jsr_pathcomplete>fullHelp for an     #
% #  extended help                                                        #
% #########################################################################
%
%   [Ub] = JSR_PATHCOMPLETE(M)
%       M must be a cell array of real square matrices of same dimension.
%       Computes an upper bound of the JSR for the set of matrices M with
%       a de Bruijn graph of dimension 1. 
%       This upper bound is inherent to the graph used, and does not
%       always achieve the real JSR.
%
%   [Ub, P] = JSR_PATHCOMPLETE(M)
%       P is a cell array of positive definite matrices.
%       If there is an edge [i,j,k], then the following inequality holds :
%       M{k}'*P{j}*M{k} <= Ub^2 * P{i}.
%
%   [Ub, P, INFO] = JSR_PATHCOMPLETE(M)
%       INFO is a structure containing various data about the iterations.
%
%   [ ... ] = JSR_PATHCOMPLETE(M,graph)
%       Computes an upper bound for the JSR of M using the specified graph.
%       Graph can be a structure or a tensor.
%       See tens2graph for more information.
%       If this input is empty or equal to zero, then the default graph
%       is generated.
%
%     WARNING : This method only returns a valid upper bound if the entered graph 
%     is path-complete! The method DOES NOT verify the path-completeness of the 
%     graph because it is a computationally expensive problem.
%
%
%   [ ... ] = JSR_PATHCOMPLETE(M, graph, options) 
%   or [ ... ] = JSR_PATHCOMPLETE(M, options)
%       The field opts.pathcomplete (generated by jsrsettings) can be used
%       to tune the method (tolerance, maxTime, etc)
%
%   Warning : the algorithm stops when the length of the bisection
%   interval is smaller than the tolerance. It does NOT mean that
%   (Ub - realJSR) < tol !
% _______________________
%
%   The upper bound found is strongly dependent of the graph used.
%   The field graphOptions can be changed to change the graph.
%   If a graph is specified in input, then it overrides all the options below.
% 
%       graphOptions.type       - String or char.
%                                 Specify the type of graph used.
%                                 Default : 'debruijn'.
%                                 See below for the different type of graph
%                                 and their options.
%
%           # 'DeBruijn' #        >> help debruijn for more information
%
%              debruijn.dimension Non-negative integer.
%                                 Sets the dimension of the de Bruijn graph 
%                                 when graphOptions.type = 'debruijn'. 
%                                 By default, equal to one.
%
%   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.
%
% See also : SOLVE_SEMI_DEFINITE_PROGRAM, TENS2GRAPH, DEBRUIJN, JSRSETTINGS,
%            DEMO3_JSR JSR_PATHCOMPLETE>FULLHELP.


timeBeg = cputime;
graphInInput = false;
fileLoaded = false;
info.status = 2; % Stopped during execution.
info.stopFlag = 5; % Stopped for a unknown reason.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                           PRE PROCESS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if(nargin < 2)
    options = jsrsettings;
end

if(nargin > 1)
    if(length(varargin) == 1) % options only
        options = varargin{1};
    end
    
    if(length(varargin) == 1)  
        if(isfield(varargin{1},'nNodes') && isfield(varargin{1},'edges')) % graph only
            checkGraph(varargin{1})
            options = jsrsettings;
            graphInInput  = true;
            graph = varargin{1}; 
        else                            % options only
            options = varargin{1};
        end
    end
    if(length(varargin) == 2) % options + graph
        graphInInput  = true;
        graph = varargin{1}; 
        options = varargin{2};
    end
end


% Checks if we have to convert graph
if graphInInput && not(isstruct(graph))
    d = size(graph);
    if(length(d) == 3)
        graph = tens2graph(graph);
    else
        error('Input graph must be a tensor or a correct structure.')
    end
end 

% load after each iterations
if(options.pathcomplete.loadIt)
    fileLoaded = true;
    load(options.pathcomplete.loadItFile,'M','graph','A','b','c','K','Agamma');
end


if not(iscell(M))
    error('M must be a cell array!');
end

% If M is not a cell vector, reshape
if(size(M,1) ~= 1 && size(M,2) ~= 1)
    warning('Reshape M to a vector of cells.')
    M = reshape(M,1,numel(M)); % matrix to vec.
end

if iscellcomplex(M)
    error('Works only with real matrix.')
end

sizeSetM = length(M);
sizeMatrix = size(M{1},1);

% Checks dimensions
for i=1:sizeSetM
    if(size(M{i},1) ~=sizeMatrix || size(M{i},2) ~=sizeMatrix)
        error('M must be a cell array of square matrices of same dimension.')
    end
end

% logfile opening
close =1;
oldLogFile = options.logfile;
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_pathcomplete','wt');
        if (logFile == -1)
            warning('Could not open logfile')
        end
    else
        logFile = options.logfile;
        close =0;
    end
else
    logFile = fopen('log_pathcomplete','wt');
    if (logFile == -1)
        warning('Could not open logfile')
    end
end

if (logFile~=-1)
    fprintf(logFile,[datestr(now) '\n\n']);
end
options.logfile = logFile; % To give logfile to subroutine.

msg(logFile,options.verbose>1,'\n \n******** Starting jsr_pathcomplete ******** \n \n')


% Checks if we have to generate a graph
nLabels = length(M);
if graphInInput || fileLoaded
    checkGraph(graph,nLabels)
else
    msg(logFile,options.verbose>1,'pathcomplete : create graph with generate_graph');
    graph = generate_graph(nLabels,options);
    msg(logFile,options.verbose>1,'generate_graph : Done.\n');
end

% Generates SDP if there is no file to load.
if not(fileLoaded)
    msg(logFile,options.verbose>1,'pathcomplete : create SDP problem with generate_pathcomplete_sdp');
    [A,Agamma,b,c,K] = generate_pathcomplete_sdp(M,graph,options); % See help for more info.
    msg(logFile,options.verbose>1,'generate_pathcomplete_sdp : Done.\n');
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                           MAIN LOOP
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
msg(logFile,options.verbose>1,'\njsr_pathcomplete : try to find optimal UB.\n')
try
    [Ub, P, trash, info] = find_pathcomplete_lyapunov(M,A,Agamma,b,c,K,graph,options);
catch errorStruct
    % If there is an error, then we throw it and restore options.logfile
    options.logfile = oldLogFile;
    msg(logFile,options.verbose>0,'jsr_pathcomplete : Error occured while executing find_pathcomplete_lyapunov.\n');
    msg(logFile,options.verbose>0,['Error msg : ',errorStruct.message]);
    errorfile = errorStruct.stack.file;
    errorlineNum = errorStruct.stack.line;
    errorline = num2str(errorlineNum);
    msg(logFile,options.verbose>0,['In file : ',errorfile]);
    msg(logFile,options.verbose>0,['At line : ',errorline]);
    errorStruct.rethrow
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                           Post process
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%�

msg(logFile,options.verbose>1,'find_pathcomplete_lyapunov : Done.\n')

info.status = 0;

if(info.stopFlag==3) % max ite reached
    msg(logFile,options.verbose>0,'\njsr_findpathcomplete : Maximum iterations reached. The solution is correct but may be inacurate\n')
    info.status = 1;
end
if(info.stopFlag == 4) % maxtime reached
    msg(logFile,options.verbose>0,'\njsr_findpathcomplete : Maximum time reached. The solution is correct but may be inacurate\n')
    info.status = 2;
end

info.opts = options;
options.logfile = oldLogFile;

if ischar(options.saveEnd)
    save([options.saveEnd,'.mat'],'Bounds','P','info');
end

info.elapsedtime = cputime - timeBeg;
msg(logFile,options.verbose>1,'\njsr_pathcomplete : End of the algorithm after %5.2f seconds.\n',info.elapsedtime)

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



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                           INTERNAL FUNCTIONS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%
% CHECK GRAPH %
%%%%%%%%%%%%%%%

% Just create error if there is something wrong.

    function checkGraph(graph,nLabels) % check if graph is a correct structure
        
        if(  not(isfield(graph,'nNodes')) )
            error('The field nNodes does not exist in the structure "graph"');
        end
        
        if( graph.nNodes < 1 )
            error('The variable graph does not have a positive value for nNodes field.');
        end
        
        if( not(isfield(graph,'edges')))
            error('The field edges does not exist in the structure "graph"');
        end
        
        if( size(graph.edges,2) ~= 3 )
            error('The field edges is not a nEdgesx3 matrix');
        end
        
        if(max(max(graph.edges(:,1:2))) > graph.nNodes)
            error('There exist a label of node which is greater than nNodes');
        end
        
        if(sum(graph.edges(:,:)<=0) > 0)
            error('Field graph.edges must contain only positive entries.')
        end
        
        if sum(not(round(graph.edges(:,:))==graph.edges(:,:)))>0
            error('Field graph.edges must contain only integers.')
        end
        
        if(sum(graph.edges(:,:)<=0) > 0)
            error('Field graph.edges must contain only positive entries.')
        end
        
        if(max(graph.edges(:,3)) > nLabels )
            error('There exist a edge label larger than nLabels');
        end
    end


end

function fullHelp()
% JSR_PATHCOMPLETE  Computes an upper bound for the JSR with bisection 
%                   method and provides also multiple quadratic Lyapunov 
%                   functions for the set of matrices M with the help
%                   of a path-complete automaton. 
%                   The labels of this automaton are matrices of M, 
%                   and the nodes are symetric positive definite 
%                   matrices used to build the Lyapunov functions.
%                   See demo3_jsr for a quick start guide.
%
% #########################################################################
% #  Full text for jsr_pathcomplete                                       #
% #########################################################################
%
%   [Ub] = JSR_PATHCOMPLETE(M)
%       M must be a cell array of real square matrices of same dimension.
%       Computes an upper bound of the JSR for the set of matrices M with
%       a de Bruijn graph of dimension 1. 
%       This upper bound is inherent to the graph used, and does not
%       always achieve the real JSR.
%
%   [Ub, P] = JSR_PATHCOMPLETE(M)
%       Returns also matrices present in multiple quadratic Lyapunov
%       functions. P is a cell array of positive definite matrices.
%       If there is an edge [i,j,k], then the following inequality holds :
%       M{k}'*P{j}*M{k} <= Ub^2 * P{i}.
%
%   [Ub, P, info] = JSR_PATHCOMPLETE(M)
%       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    Total elapsed time
%
%           info.bisec          [lower upper] last bisection interval
%
%           info.niter          Number of iterations (bisection)
%
%           info.iterTime       Vector containing the elapsed
%                               time of each iteration
%
%           info.opts           The structure of options used
%
%           info.iterGamma      Value of gamma tested at each
%                               iterations (see [1]), with gamma = 1/Ub.
%
%           info.iterFeas       Binary vector. The kth entry is true
%                               when info.iterGamma(k) is feasible
%
%           info.primalVar      Primal variables of the last feasible SDP
%
%           info.dualVar        Dual variables of the last feasible SDP
%
%           info.iterObj        Vector with all values of the objective 
%                               function at each iteration
%
%           info.stopFlag       More detailed than info.status.
%                           = 0 if the algorithm converges without any
%                               numerical problem
%                           = 1 if Ub is found with the required tolerance
%                               but not P (you can change the tolerance in
%                               options.solveroptions, see help
%                               solve_semi_definite_program)
%                           = 2 if there is a complete numerical failure
%                               (the solver is not able to solve the SDP
%                               problem properly)
%                           = 3 if the number of iterations is greater than
%                               maxiter value (see options below)
%                           = 4 if maxTime is reached (see options below)
%                           = 5 If the method stops for an unknown reason
%                           =-1 if the given or computed UB is not
%                               feasible. In this case, this UB will be
%                               returned but P will be empty
%
%
%   [ ... ] = JSR_PATHCOMPLETE(M,graph)
%       Computes an upper bound for the JSR of M using the specified graph.
%       Graph can be a structure or a tensor.
%       See tens2graph for more information.
%       If this input is empty or equal to zero, then the default graph
%       is generated.
%
%     WARNING : This method only returns a valid upper bound if the entered graph 
%     is path-complete! The method DOES NOT verify the path-completeness of the 
%     graph because it is a computationally expensive problem.
%
%
%   [ ... ] = JSR_PATHCOMPLETE(M, graph, options) 
%   or [ ... ] = JSR_PATHCOMPLETE(M, options)
%       The field opts.pathcomplete (generated by jsrsettings) can be used
%       to tune the method.
%
%       pathcomplete.LbBisec    - Real positive number.
%                                 Gives starting LB for bisection method.
%                                 If not specified, takes the biggest
%                                 spectral radius of each matrices.
%       pathcomplete.UbBisec    - Real positive number.
%                                 Gives starting UB for bisection method.
%                                 If not specified, takes the biggest
%                                 norm 2 of each matrices.
%       pathcomplete.maxiter    - Positive integer.
%                                 Maximum number of iterations (for
%                                 bisection), (inf).
%       pathcomplete.reltol     - Real positive number.
%                                 Maximum relative length for bisection interval. 
%                                 Stopping criterion is
%                                 UbBisec-LbBisec < options.pathcomplete.reltol*UbBisec
%                                 Default : 1e-6.
%       pathcomplete.abstol     - Real positive number.
%                                 Maximum absolute length for bisection interval. 
%                                 Stopping criterion is
%                                 UbBisec-LbBisec < options.pathcomplete.abstol
%                                 Default : +inf.
%   Warning : the algorithm stops when the length of the bisection
%   interval is smaller than the tolerance. It does NOT mean that
%   (Ub - realJSR) < tol !
%
%       pathcomplete.testUb     - Binary value.
%                                 If testUb =1, then the method checks if
%                                 the upper bound for the JSR is feasible.
%                                 By default, testUb =1. 
%       pathcomplete.testLb     - Binary value.
%                                 If testLb =1, then the method checks if
%                                 the lower bound for the JSR is infeasible.
%                                 By default, testLb =1. 
%       pathcomplete.loadIt     - Binary value.
%                                 If loadIt =1, then the  method loads the
%                                 file with name [pathcomplete.loadItFile]
%                                 and resumes from the last finished 
%                                 iteration.
%                                 See options 'saveinIt' in jsrsettings to
%                                 know how to save data after each iterations.
%                                 In case where the problem is big, it is
%                                 strongly recommended to enable the
%                                 options saveinIt.
%                                 By default, loadIt = 0. 
%       pathcomplete.loadItFile - String.
%                                 By default, loadItFile = 'saveinIt.mat'
%                                 See pathcomplete.loadIt.
%
%   The upper bound found is strongly dependent of the graph used.
%   The field graphOptions can be changed to change the graph.
%   If a graph is specified in input, then it overrides all the options below.
% 
%       graphOptions.type       - String or char.
%                                 Specify the type of graph used.
%                                 Default : 'debruijn'.
%                                 See below for the different type of graph
%                                 and their options.
%
%           # 'DeBruijn' #        >> help debruijn for more information
%
%              debruijn.dimension Non-negative integer.
%                                 Sets the dimension of the DeBruijn graph 
%                                 when graphOptions.type = 'debruijn'. 
%                                 By default, equal to one.
%
%   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.
% _________________________________________________________________________
%
%   We can deduce a lower bound with the output of this algorithm (see [1]),
%   but this lower bound is quite inacurate and is not helpful.
% _________________________________________________________________________
%
% REFERENCES
%   [1] Ahmadi, Jungers, Parrilo and Roozbehani,
%   "Joint spectral radius and path-complete graph Lyapunov functions"
%   Vol. 52, No1, pp. 687-717 in SIAM J. CONTROL OPTIM, 2014.
%
% See also : SOLVE_SEMI_DEFINITE_PROGRAM, TENS2GRAPH, DEBRUIJN, JSRSETTINGS,
%            DEMO3_JSR.
end
