

% The following code is property of Aimdyn, Inc.
% Version of July 27, 2011
% All rights reserved

% Inputs: 
% matr - weighted adjacency matrix

% Outputs:
% sum_1 - the first term in the complexity function
% sum_2 - the second term in the complexity 1 function
% sum_2_second - the second term in the complexity 2 function
% metrics - complexity 1
% metrics_second - complexity 2
% sinks - sinks
% kmin - number of groups not including sinks as a group
% kmin1 - number of groups including sinks as a group
% answer_all - groups
% ncomp_all size of groups
% cycles_ind - number of generating eigenvalues
% nodes - original nodes


function [sum_1, sum_2, sum_2_second, metrics, metrics_second, sinks, kmin, kmin1, answer_all,...
          ncomp_all, cycles_ind, nodes] = ComputeComplexityAndClustering(matr)

    showplots = false;
    I=sqrt(-1);

    % removing sources
    nodes = 1:size(matr,1);
    csum = sum(matr,1);
    rsum = sum(matr,2);
    rmv = (csum == 0) & (rsum' ~= 0);
    while sum(rmv)>0;
        matr = matr(~rmv,~rmv);
        nodes(rmv) = [];
        csum = sum(matr,1);
        rsum = sum(matr,2);
        rmv = (csum == 0) & (rsum' ~= 0);
    end

    % finding sinks
    matr0 = matr;
    notsinks = 1:size(matr0,1);
    csum = sum(matr0,1);
    rsum = sum(matr0,2);
    rmv = (rsum == 0) & (csum' ~= 0);
    while sum(rmv)>0;
        matr0 = matr0(~rmv,~rmv);
        notsinks(rmv) = [];
        csum = sum(matr0,1);
        rsum = sum(matr0,2);
        rmv = (rsum == 0) & (csum' ~= 0);
    end
    sinks = setdiff(1:size(matr,1),notsinks);

    % normalizing adjacency matrix
    size_new=size(matr,1);
    matr1=zeros(size_new,size_new);
    for j=1:size_new
        sum1=sum(matr(j,:));
        if (sum1>0)
            matr1(j,:)=matr(j,:)/sum1;
        else
            matr1(j,j)=1;
        end
    end

    % plotting graph
    if (showplots)
        nnames = cell(size_new,1);
        for i2 = 1:size_new
            nnames{i2} = ['Node ' int2str(nodes(i2))];
        end
        bg2 = biograph(matr1,nnames);
        view(bg2)
    end

    % computing eigenvalues and eigenvectors
    [V,Evd] = eig(matr1);
    eigval = zeros(size(Evd,1),1);
    for i=1:size(Evd,1)
        eigval(i) = Evd(i,i);
    end
    notzeros=1:size_new;
    notzeros(abs(eigval)<1e-8)=[];

    % plotting eigenvalues
    if (showplots)
        figure;
        hold on
        circle = rsmak('circle');
        fnplt(circle)
        plot(real(eigval(notzeros)),imag(eigval(notzeros)),'.r','MarkerSize',10)
        daspect([1 1 1]);
        hold off
        title('Nonzero Eigenvalues')
    end

    % computing first term in complexity formula
    size1=length(notzeros);
    maxfun=zeros(size1,1);
    cou=1;
    for j=notzeros
        maxfun(cou)=(1-abs(eigval(j)));
        cou=cou+1;
    end
    if (size1>=1)
        sum_1=sum(maxfun)/size1;
    else
        sum_1=0;
    end

    % computing second term in complexity formula (with 2 methods)
    fe=zeros(size_new,1);
    for z=1:size_new
        fe(z)=atan2(imag(eigval(z)),real(eigval(z)));
    end
    fe(fe<0) = fe(fe<0)+2*pi;

    filim=zeros(size1,1);
    filim(1)=0.5*2*pi/size1;
    for z=2:size1
       filim(z)=filim(z-1)+2*pi/size1;
    end

    i1=zeros(size1,1);
    for z=1:size1
       if (z~=size1)
           for t=notzeros
               if ((fe(t)<=filim(z+1))&&(fe(t)>filim(z)))
                   i1(z)=i1(z)+1;
               end
           end
       else
           for t=notzeros
               if ((fe(t)>filim(z))&&(fe(t)<=2*pi)) || ((fe(t)>=0)&&(fe(t)<=filim(1)))
                   i1(z)=i1(z)+1;
               end
           end
       end
    end

    maxlog=zeros(size1,1);
    for j=1:size1
        if (size1>1)
            if (i1(j)==0)
                maxlog(j)=0;
            else
                maxlog(j)=(i1(j)/size1)*log(i1(j)/size1);
            end
        else
            maxlog(j)=0;
        end
    end

    if (size1>1)
        sum_2=-sum(maxlog)/log(size1);
    else
        sum_2=0;
    end

    metrics=sum_1+sum_2;

    if (size1>0)
        sum_2_second=1-i1(size1)/size1;
    else
        sum_2_second=0;
    end

    metrics_second=sum_1+sum_2_second;

    % compute clustering
    disj_comp=[];
    for i1=1:size_new
        if ((sum(matr1(:,i1))==1)&&(matr1(i1,i1)==1))
            disj_comp=[disj_comp i1];
        end
    end

    if (length(disj_comp)<size_new)

        s = zeros(1,size1-1);
        sm = zeros(1,size1-1);
        evs1 = eigval(notzeros);

        angl = atan2(imag(evs1),real(evs1));
        angl(angl<0) = angl(angl<0)+2*pi;

        for k=2:size1
            for t=2:k
                angidx = find(angl >= (2*pi/k)*(t-1.5) & angl <= (2*pi/k)*(t-0.5));
                if isempty(angidx)
                    sm(t-1)=1;
                else
                    sm(t-1)=min(abs(exp(I*(2*pi/k)*(t-1))-evs1(angidx)));
                end
            end
            s(k-1)=(1/(k-1))*sum(sm(1:(k-1)));
        end
        [smin kt] = min(s);
        % number of clusters
        if (size1<2)
            kmin = 1;
        else
            kmin=kt+1;
        end

        % find generating eigenvalue(s)
        angidx1 = angl >= (2*pi/kmin)*0.5 & angl <= (2*pi/kmin)*1.5;
        evs2 = evs1;
        evs2(~angidx1) = Inf;
        arr2=abs(evs2-exp(I*(2*pi/kmin)));
        arr2_min=min(arr2);
        [a_m sj1_m]=sort(arr2);
        sj1_min = sj1_m(a_m<arr2_min+1e-4);

        size_sj1_min=length(sj1_min);
        cycles_ind=size_sj1_min;

        eigvec_sj1=zeros(size_new,size_sj1_min);
        for q1=1:size_new
            eigvec_sj1(q1,:)=V(q1,notzeros(sj1_min));
        end

        notsinksdisj = setdiff(notsinks,disj_comp);
        for q_sj1=1:size_sj1_min

            answer=cell(1,0);
            ncomp=zeros(1,0);

            % plot eigenvector
            if (showplots)
                figure;
                hold on
                circle = rsmak('circle');
                fnplt(circle)
                plot(real(eigvec_sj1(notsinks,q_sj1)),imag(eigvec_sj1(notsinks,q_sj1)),'.r','MarkerSize',10)
                daspect([1 1 1]);
                hold off
                title('Notsink components of Eigenvector')
            end

            fi=zeros(size_new,1);
            for z=1:size_new
               fi(z)=atan2(imag(eigvec_sj1(z,q_sj1)),real(eigvec_sj1(z,q_sj1)));
            end
            fi(fi<0) = fi(fi<0)+2*pi;

            filim=zeros(1,kmin);
            filim(1)=0.5*2*pi/kmin;
            for z=2:kmin
               filim(z)=filim(z-1)+2*pi/kmin;
            end

            for z=1:kmin
                i1=1;
                comp=[];
                if (z~=kmin)
                   for t=notsinksdisj
                       if ((fi(t)<=filim(z+1))&&(fi(t)>filim(z)))
                           comp(i1)=t;
                           i1=i1+1;
                       end
                   end
                else
                   for t=notsinksdisj
                       if ((fi(t)>filim(z))&&(fi(t)<=2*pi)) || ((fi(t)>=0)&&(fi(t)<=filim(1)))
                           comp(i1)=t;
                           i1=i1+1;
                       end
                   end
                end
                if (~isempty(comp))
                    answer = [answer {comp}];
                    ncomp = [ncomp length(comp)];
                end
            end
            answer_all{q_sj1}=answer;
            ncomp_all{q_sj1}=ncomp;
        end

        if (~isempty(disj_comp))
            kmin=kmin+length(disj_comp);
            for q_sj1=1:size_sj1_min
                sizesize=zeros(1,size_sj1_min);
                sizesize(q_sj1)=size(answer_all{q_sj1},2);
                for i2=1:length(disj_comp)
                    answer_all{q_sj1}{sizesize(q_sj1)+i2}=disj_comp(i2);
                    ncomp_all{q_sj1}=[ncomp_all{q_sj1} 1];
                end
            end
        end

        kmin1=kmin;
        if (~isempty(sinks))
            kmin1=kmin1+1;
        end

    else
        answer={1};
        ncomp=1;
        for i1=2:size_new
            comp=i1;
            answer = [answer {comp}];
            ncomp = [ncomp length(comp)];
        end
        answer_all{1}=answer;
        ncomp_all{1}=ncomp;
        kmin=size_new;
        kmin1=kmin;
        cycles_ind=1;
    end
end
