function k=loop_complexity(M)
% first, convert automaton M to graph G
nQ = length(M.Q);
G = zeros(nQ);
for q=M.Q
    [qq1,ss1] = get_ddelta(M,q);
    for q1=qq1
        G(q,q1) = 1;
    end
end
k = lc_rec(G,1:nQ);
return

function k = lc_rec(G,nodes)
% definition taken from
% http://www.liafa.jussieu.fr/~lombardy/publi/ICWLC.pdf
BB = get_balls(G,nodes);
if isempty(BB)
    k = 0;
    return
end
if (length(BB)==1)
    % G is a ball
    K = zeros(size(nodes));
    for j=1:length(nodes)
        K(j) = lc_rec(G,setdiff(nodes,nodes(j)));
    end
    k = 1+min(K);
else
    K = zeros(1,length(BB));
    for j=1:length(BB)
        K(j) = lc_rec(G,BB{j});
    end
    k = max(K);
end
return


function BB = get_balls(G,nodes)
% a ball is a strongly connected component of G
% that contains an arc
[n,n] = size(G);
Id = eye(n);
C = ((Id+G)^(n-1)>0);
% C(i,j) = 1 iff there is a directed path from i to j in G
%nodes = 1:n;
BB = {};
while ~isempty(nodes)
    v = nodes(1);
    B = sc_comp(C,v);
    nodes = setdiff(nodes,B);
    if (G(v,v) | (length(B)>1))
        BB{end+1}=B;
    end
end
return

function B = sc_comp(C,v)
% B is the sc-component containing v
C = C.*C';
B = find(C(v,:));
return