
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%          NetInfCont
%                      A program for calculating the information
%                      content of a network.
%                      Developed by M. Zanin.
%
%                      V1.01
%
%                      Input arguments:
%                          AM: binary adjacency matrix
%                          directed: 1 if the network is directed,
%                                    0 otherwise
%
%                      Outputs:
%                          Info: the IC of the network
%                          NormInfo: the normalized IC
%
%                      Please cite the following paper:
%                          M. Zanin, P. A. Sousa and E. Menasalvas,
%                          Information content: assessing meso-scale structures in complex networks
%                          EPL 106 (3), (2014) 30001
%


function [Info, NormInfo] = NetInfCont( AM, directed )

numNodes = size(AM, 1);
linkDensity = sum( sum( AM ) ) / ( numNodes * numNodes );
Info = 0;

DistanceMatrix = zeros(numNodes, numNodes);
for n1 = 1 : numNodes
    for n2 = 1 : numNodes
        
        if n1 == n2
            DistanceMatrix(n1, n2) = numNodes;
            continue;
        end
        
        DistanceMatrix(n1, n2) = HD( AM(n1, :), AM(n2, :), AM(:, n1), AM(:, n2) );
                                            
    end
end




for k = 1 : (numNodes - 1)
    [tInfo, AM, DistanceMatrix] = OneIteration(AM, DistanceMatrix);
    Info = Info + tInfo;
end


NormInfo = Info / Normalize( linkDensity, numNodes, directed );


end





function [Info, newAM, newDM] = OneIteration(AM, DistanceMatrix)

numNodes = size(AM, 1);

minDist = min( min( DistanceMatrix ) );
maxDist = max( max( DistanceMatrix ) );

if minDist < ( 2 * numNodes - maxDist )
    [aX, bX] = min( DistanceMatrix, [], 1 );
    [bestHD, cX] = min( aX );
else
    [aX, bX] = max( DistanceMatrix, [], 1 );
    [bestHD, cX] = max( aX );    
end

bestNode = bX(cX);
bestNode2 = cX;

if bestNode == bestNode2
    bestNode = 1;
    bestNode2 = 2;
end



p = bestHD / ( numNodes * 2 );
Info = 0;
if p > 0 && p < 1
    Info = - p * log2( p ) - (1.0 - p) * log2( 1.0 - p );
    Info = Info * numNodes * 2;
end




newDM1 = DistanceMatrix;

for n1 = 1 : numNodes
    for n2 = 1 : numNodes
        
        if n1 == n2
            continue;
        end
        
        newDM1(n1, n2) = newDM1(n1, n2) + double( AM(n1, bestNode) == AM(n2, bestNode) ) - 1;
        newDM1(n1, n2) = newDM1(n1, n2) + double( AM(bestNode, n1) == AM(bestNode, n2) ) - 1;
        
    end
end
offset1 = [ 1 : (bestNode - 1) (bestNode + 1) : numNodes ];
newDM1 = newDM1( offset1, offset1 );



newDM2 = DistanceMatrix;

for n1 = 1 : numNodes
    for n2 = 1 : numNodes
        
        if n1 == n2
            continue;
        end
        
        newDM2(n1, n2) = newDM2(n1, n2) + double( AM(n1, bestNode2) == AM(n2, bestNode2) ) - 1;
        newDM2(n1, n2) = newDM2(n1, n2) + double( AM(bestNode2, n1) == AM(bestNode2, n2) ) - 1;
        
    end
end
offset2 = [ 1 : (bestNode2 - 1) (bestNode2 + 1) : numNodes ];
newDM2 = newDM2( offset2, offset2 );


maxV1 = max( max( newDM1 ) );
maxV2 = max( max( newDM2 ) );
minV1 = min( min( newDM1 ) );
minV2 = min( min( newDM2 ) );

if maxV1 < ( 2 * ( numNodes - 1 ) - minV1 )
    maxV1 = 2 * ( numNodes - 1 ) - minV1;
end
if maxV2 < ( 2 * ( numNodes - 1 ) - minV2 )
    maxV2 = 2 * ( numNodes - 1 ) - minV2;
end

if maxV1 > maxV2
    newAM = AM( offset1, offset1 );
    newDM = newDM1;
else  
    newAM = AM( offset2, offset2 );
    newDM = newDM2;
end


end
        
        
        


function answ = HD( vector1, vector2, vector3, vector4 )

    sizeVector = max( size( vector1 ) );
    answ = sum( double( vector1 == vector2 ) ) + sum( double( vector3 == vector4 ) );
    
    answ = 2 * sizeVector - answ;

end



function approxInfo = approxF( p )

    if p == 0 || p == 1
        approxInfo = 0;
        return;
    end


    approxInfo = - p * log2(p) - (1 - p) * log2(1 - p);

end



function Info = Normalize( p, numNodes, directed )

if directed == 1
    KMax = 0.1625 - 1.822 * ( numNodes + 10.08 ) ^ ( -1.098 );
else
    KMax = 0.1622 - 0.6298 * ( numNodes + 5.905 ) ^ ( -1.132 );
end

alpha = KMax * 4;
effP = alpha * p - alpha * p * p;

Info = 0;

for k = numNodes : -1 : 2
    temp = (2 * k) * approxF( effP );
    Info = Info + temp;
end


end




