Java tutorial
/******************************************************************************* * Copyright (c) 2016 Pablo Pavon-Marino. * All rights reserved. This program and the accompanying materials * are made available under the terms of the GNU Lesser Public License v2.1 * which accompanies this distribution, and is available at * http://www.gnu.org/licenses/lgpl.html * * Contributors: * Pablo Pavon-Marino - Jose-Luis Izquierdo-Zaragoza, up to version 0.3.1 * Pablo Pavon-Marino - from version 0.4.0 onwards ******************************************************************************/ package com.net2plan.libraries; import cern.colt.matrix.tdouble.DoubleFactory1D; import cern.colt.matrix.tdouble.DoubleFactory2D; import cern.colt.matrix.tdouble.DoubleMatrix1D; import cern.colt.matrix.tdouble.DoubleMatrix2D; import cern.colt.matrix.tdouble.algo.DenseDoubleAlgebra; import cern.colt.matrix.tdouble.algo.decomposition.DenseDoubleEigenvalueDecomposition; import cern.jet.math.tdouble.DoubleFunctions; import com.net2plan.interfaces.networkDesign.Link; import com.net2plan.interfaces.networkDesign.Node; import com.net2plan.libraries.GraphUtils.JGraphTUtils; import com.net2plan.libraries.GraphUtils.JUNGUtils; import com.net2plan.utils.CollectionUtils; import com.net2plan.utils.Constants; import com.net2plan.utils.DoubleUtils; import edu.uci.ics.jung.algorithms.scoring.BetweennessCentrality; import edu.uci.ics.jung.graph.Graph; import org.apache.commons.collections15.Transformer; import org.jgrapht.DirectedGraph; import org.jgrapht.alg.EdmondsKarpMaximumFlow; import org.jgrapht.alg.StrongConnectivityInspector; import java.util.*; /** * <p>Class to deal with graph-theory metrics computation.</p> * * <p><b>Important</b>: Internal computations (like shortest-paths) are cached in order to improve efficiency.</p> * * @author Pablo Pavon-Marino, Jose-Luis Izquierdo-Zaragoza */ @SuppressWarnings("unchecked") public class GraphTheoryMetrics { private final List<Node> nodes; private final List<Link> linkMap; private Map<Link, Double> costMap; private final int N; private final int E; private DoubleMatrix2D adjacencyMatrix; private double[] adjacencyMatrixEigenvalues; private double averageSPLength, diameter, heterogeneity; private DirectedGraph<Node, Link> graph_jgrapht; private Graph<Node, Link> graph_jung; private DoubleMatrix2D incidenceMatrix; private DoubleMatrix2D laplacianMatrix; private double[] laplacianMatrixEigenvalues; private Transformer<Link, Double> nev; private DoubleMatrix1D linkBetweenessCentrality; private DoubleMatrix1D nodeBetweenessCentrality; private DoubleMatrix1D outNodeDegree; /** * Default constructor * @param nodes List of odes * @param links List of links * @param linkCostMap Cost per link, where the key is the link identifier and the value is the cost of traversing the link. No special iteration-order (i.e. ascending) is required */ public GraphTheoryMetrics(List<Node> nodes, List<Link> links, Map<Link, Double> linkCostMap) { this.nodes = nodes; this.linkMap = links; this.N = nodes.size(); this.E = links.size(); adjacencyMatrix = null; adjacencyMatrixEigenvalues = null; graph_jgrapht = null; graph_jung = null; incidenceMatrix = null; laplacianMatrix = null; laplacianMatrixEigenvalues = null; linkBetweenessCentrality = null; nev = null; nodeBetweenessCentrality = null; outNodeDegree = null; configureLinkCostMap(linkCostMap); } private void computeBetweenessCentrality() { Graph<Node, Link> aux_graph = getGraph_JUNG(); BetweennessCentrality bc = new BetweennessCentrality(aux_graph, getCostTransformer()); nodeBetweenessCentrality = DoubleFactory1D.dense.make(N); for (Node node : this.nodes) nodeBetweenessCentrality.set(node.getIndex(), bc.getVertexScore(node)); linkBetweenessCentrality = DoubleFactory1D.dense.make(E); for (Link link : linkMap) linkBetweenessCentrality.set(link.getIndex(), bc.getEdgeScore(link)); } private void computeSPDistanceMetrics() { diameter = 0; averageSPLength = 0; heterogeneity = 0; Graph<Node, Link> aux_graph = getGraph_JUNG(); Transformer<Link, Double> aux_nev = getCostTransformer(); /* Compute network diameter using nave Floyd-Warshall algorithm */ double[][] costMatrix = new double[N][N]; for (int n = 0; n < N; n++) { Arrays.fill(costMatrix[n], Double.MAX_VALUE); costMatrix[n][n] = 0; } for (Link edge : aux_graph.getEdges()) { int a_e = edge.getOriginNode().getIndex(); int b_e = edge.getDestinationNode().getIndex(); double newCost = aux_nev.transform(edge); if (newCost < costMatrix[a_e][b_e]) costMatrix[a_e][b_e] = newCost; } for (int k = 0; k < N; k++) { for (int i = 0; i < N; i++) { if (i == k) continue; for (int j = 0; j < N; j++) { if (j == k || j == i) continue; double newValue = costMatrix[i][k] + costMatrix[k][j]; if (newValue < costMatrix[i][j]) costMatrix[i][j] = newValue; } } } int numPaths = 0; double sum = 0; double M = 0.0; double S = 0.0; for (int i = 0; i < N; i++) { for (int j = i + 1; j < N; j++) { double dist_ij = costMatrix[i][j]; if (dist_ij < Double.MAX_VALUE) { sum += dist_ij; numPaths++; double tmpM = M; M += (dist_ij - tmpM) / numPaths; S += (dist_ij - tmpM) * (dist_ij - M); if (dist_ij > diameter) diameter = dist_ij; } double dist_ji = costMatrix[j][i]; if (dist_ji < Double.MAX_VALUE) { sum += dist_ji; numPaths++; double tmpM = M; M += (dist_ji - tmpM) / numPaths; S += (dist_ji - tmpM) * (dist_ji - M); if (dist_ji > diameter) diameter = dist_ji; } } } if (numPaths == 0) return; averageSPLength = numPaths == 0 ? 0 : sum / numPaths; heterogeneity = averageSPLength == 0 ? 0 : Math.sqrt(S / numPaths) / averageSPLength; } /** * Re-configures link cost setting. Related information, such as shortest paths, is cleared. * * @param linkCostMap Cost per link, where the key is the link identifier and the value is the cost of traversing the link. No special iteration-order (i.e. ascending) is required */ public void configureLinkCostMap(Map<Link, Double> linkCostMap) { if (linkCostMap == null) this.costMap = null; else this.costMap = new LinkedHashMap<Link, Double>(linkCostMap); averageSPLength = -1; diameter = -1; heterogeneity = -1; nev = null; } /** * Returns the adjacency matrix of the network. The adjacency matrix is a * <i>NxN</i> matrix (where <i>N</i> is the number of nodes in the network), * where each position (<i>i</i>,<i>j</i>) represents the number of directed * links from <i>i</i> to <i>j</i>. * * @return Adjacency matrix */ private DoubleMatrix2D getAdjacencyMatrix() { if (adjacencyMatrix == null) { adjacencyMatrix = DoubleFactory2D.sparse.make(N, N); for (Link link : linkMap) { int a_e = link.getOriginNode().getIndex(); int b_e = link.getDestinationNode().getIndex(); adjacencyMatrix.setQuick(a_e, b_e, adjacencyMatrix.get(a_e, b_e) + 1); } } return adjacencyMatrix; } /** * Returns the eigenvalues of the adjacency matrix. * * @return Eigenvalues of the adjacency matrix */ private double[] getAdjacencyMatrixEigenvalues() { if (adjacencyMatrixEigenvalues == null) { DoubleMatrix2D A_nn = getAdjacencyMatrix().copy(); A_nn.assign(A_nn.viewDice(), DoubleFunctions.max); DenseDoubleAlgebra alg = new DenseDoubleAlgebra(); DenseDoubleEigenvalueDecomposition eig = alg.eig(A_nn); adjacencyMatrixEigenvalues = eig.getRealEigenvalues().toArray(); DoubleUtils.sort(adjacencyMatrixEigenvalues, Constants.OrderingType.ASCENDING); } return adjacencyMatrixEigenvalues; } /** * <p>Returns the algebraic connectivity of the network. The algebraic connectivity * is equal to the second smallest eigenvalue of the laplacian matrix.</p> * * <p>For symmetric (or undirected) networks, if the algebraic connectivity * is different from zero, it is ensured that the network is connected, that is, * it is possible to find a path between each node pair. * * @return Algebraic connectivity */ public double getAlgebraicConnectivity() { double[] eig = getLaplacianMatrixEigenvalues(); return eig[1]; } /** * Returns the assortativity of the network. * * @return Assortativity */ public double getAssortativity() { if (E == 0) return 0; DoubleMatrix1D aux_outNodeDegree = getOutNodeDegree(); double a = 0; double b = 0; double y = 0; for (Link link : linkMap) { Node originNode = link.getOriginNode(); Node destinationNode = link.getDestinationNode(); int j_e = (int) aux_outNodeDegree.get(originNode.getIndex()); int k_e = (int) aux_outNodeDegree.get(destinationNode.getIndex()); y += j_e + k_e; a += j_e * k_e; b += j_e * j_e + k_e * k_e; } y /= 2.0D * E; y *= y; a /= E; b /= 2.0D * E; return (a - y) / (b - y); } /** * Returns the average neighbor connectivity. * * @return Average neighbor connectivity */ public double getAverageNeighborConnectivity() { if (E == 0) return 0; DoubleMatrix1D aux_outNodeDegree = getOutNodeDegree(); int maxNodeDegree = aux_outNodeDegree.size() == 0 ? 0 : (int) aux_outNodeDegree.getMaxLocation()[0]; double[] knn = new double[maxNodeDegree + 1]; DoubleMatrix2D m = DoubleFactory2D.sparse.make(maxNodeDegree + 1, maxNodeDegree + 1); for (Link link : linkMap) { Node originNode = link.getOriginNode(); Node destinationNode = link.getDestinationNode(); int degree_k1 = (int) aux_outNodeDegree.get(originNode.getIndex()); int degree_k2 = (int) aux_outNodeDegree.get(destinationNode.getIndex()); m.set(degree_k1, degree_k2, m.get(degree_k1, degree_k2) + 1); } for (int k_1 = 1; k_1 <= maxNodeDegree; k_1++) { knn[k_1] = k_1 * m.viewRow(k_1).zSum(); } return DoubleUtils.averageNonZeros(knn) / (E - 1); } /** * Returns the average number of outgoing links per node. * * @return Average number of outgoing links per node */ public double getAverageOutNodeDegree() { return getOutNodeDegree().size() == 0 ? 0 : getOutNodeDegree().zSum() / getOutNodeDegree().size(); } /** * Returns the average shortest path distance among all node-pair shortest paths. * * @return Average shortest path distance */ public double getAverageShortestPathDistance() { if (averageSPLength == -1) computeSPDistanceMetrics(); return averageSPLength; } /** * Returns the average two-term reliability (A2TR) of the network. A2TR is computed * as the ratio between the number of node-pair for which a path can be found * and the same number when the network is connected (<i>Nx(N-1)</i>, where * <i>N</i> is the number of nodes in the network). The value is in range [0, 1]. * * @return Average two-term reliability */ public double getAverageTwoTermReliability() { if (E == 0) return 0; DirectedGraph<Node, Link> graph = getGraph_JGraphT(); StrongConnectivityInspector<Node, Link> ci = new StrongConnectivityInspector<Node, Link>(graph); List<Set<Node>> connectedComponents = ci.stronglyConnectedSets(); double sum = 0; Iterator<Set<Node>> it = connectedComponents.iterator(); while (it.hasNext()) { int componentSize = it.next().size(); sum += componentSize * (componentSize - 1); } return sum / (N * (N - 1)); } /** * Returns the clustering coefficient of the network. * * @return Clustering coefficient */ public double getClusteringCoefficient() { if (E == 0) return 0; DoubleMatrix1D aux_outNodeDegree = getOutNodeDegree(); Map<Node, Double> clusteringCoefficient = new LinkedHashMap<Node, Double>(); for (Node node : nodes) { switch ((int) aux_outNodeDegree.get(node.getIndex())) { case 0: break; case 1: clusteringCoefficient.put(node, 1.0); break; default: Collection<Node> neighbors = getNeighbors(node); int aux = 0; for (Node i : neighbors) { Collection<Node> aux_neighbors = getNeighbors(i); for (Node j : neighbors) { if (i.equals(j)) continue; if (CollectionUtils.contains(aux_neighbors, j)) aux++; } } clusteringCoefficient.put(node, (double) aux / neighbors.size()); break; } } return DoubleUtils.average(clusteringCoefficient); } private Transformer<Link, Double> getCostTransformer() { if (nev == null) nev = JUNGUtils.getEdgeWeightTransformer(costMap); return nev; } /** * Returns the density of the network. The density represents the ratio * between the number of links in the network and the number of links needed * to build a full-mesh network (<i>Nx(N-1)</i>, where <i>N</i> is the number of * nodes in the network). * * @return Density */ public double getDensity() { if (N == 0) return 0; return (double) E / (N * (N - 1)); } /** * Returns the diameter of the network. The diameter is the longest path distance * among all node-pair shortest paths. * * @return Network diameter */ public double getDiameter() { if (diameter == -1) computeSPDistanceMetrics(); return diameter; } private DirectedGraph<Node, Link> getGraph_JGraphT() { if (graph_jgrapht == null) graph_jgrapht = (DirectedGraph<Node, Link>) JGraphTUtils.getGraphFromLinkMap(nodes, linkMap); return graph_jgrapht; } private Graph<Node, Link> getGraph_JUNG() { if (graph_jung == null) graph_jung = JUNGUtils.getGraphFromLinkMap(nodes, linkMap); return graph_jung; } /** * Returns the heterogeneity of the network. The heterogeneity is equal to the * standard deviation of all node-pair shortest paths divided by the average * shortest path distance. * * @return Heterogeneity */ public double getHeterogeneity() { if (heterogeneity == -1) computeSPDistanceMetrics(); return heterogeneity; } /** * Returns the incidence matrix of the network. The incidence matrix is a * <i>NxE</i> matrix (where <i>N</i> and <i>E</i> are the number of nodes * and links in the network, respectively), where each position (<i>i</i>, * <i>j</i>) is equal to: '1', if node <i>i</i> is the origin node of the link * <i>j</i>; '-1', if node <i>i</i> is the destination node of the link * <i>j</i>; and '0', otherwise. * * @return Incidence matrix */ private DoubleMatrix2D getIncidenceMatrix() { if (incidenceMatrix == null) { incidenceMatrix = DoubleFactory2D.sparse.make(N, E); for (Link link : linkMap) { int e = link.getIndex(); int a_e = link.getOriginNode().getIndex(); int b_e = link.getDestinationNode().getIndex(); incidenceMatrix.setQuick(a_e, e, 1); incidenceMatrix.setQuick(b_e, e, -1); } } return incidenceMatrix; } /** * Returns the laplacian matrix of the network. The laplacian matrix is equal * to the product of the incidence matrix by its transpose matrix. * * @return Laplacian matrix */ private DoubleMatrix2D getLaplacianMatrix() { if (laplacianMatrix == null) { DoubleMatrix2D A_ne = getIncidenceMatrix().copy(); laplacianMatrix = A_ne.zMult(A_ne.viewDice(), null); } return laplacianMatrix; } /** * Returns the eigenvalues of the laplacian matrix of the network. * * @return Eigenvalues of the laplacian matrix */ private double[] getLaplacianMatrixEigenvalues() { if (laplacianMatrixEigenvalues == null) { DenseDoubleAlgebra alg = new DenseDoubleAlgebra(); DenseDoubleEigenvalueDecomposition eig = alg.eig(getLaplacianMatrix()); laplacianMatrixEigenvalues = eig.getRealEigenvalues().toArray(); DoubleUtils.sort(laplacianMatrixEigenvalues, Constants.OrderingType.ASCENDING); } return laplacianMatrixEigenvalues; } /** * <p>Returns the betweeness centrality of each link. The betweeness * centrality of a link is equal to the number of node-pair shortest paths which traverses * the link.</p> * * <p>Internally it makes use of the Brandes' algorithm.</p> * * @return Betweeness centrality of each link */ public DoubleMatrix1D getLinkBetweenessCentrality() { if (linkBetweenessCentrality == null) computeBetweenessCentrality(); return linkBetweenessCentrality; } /** * <p>Returns the link connectivity. The link connectivity is equal to the smallest * number of link-disjoint paths between each node pair.</p> * * <p>Internally it makes use of the Edmonds-Karp algorithm to compute the maximum * flow between each node pair, assuming a link capacity equal to one for every link.</p> * * @return Link connectivity */ public int getLinkConnectivity() { if (E == 0) return 0; DirectedGraph<Node, Link> graph = getGraph_JGraphT(); EdmondsKarpMaximumFlow<Node, Node> ek = new EdmondsKarpMaximumFlow(graph); int k = Integer.MAX_VALUE; for (Node originNode : nodes) { for (Node destinationNode : nodes) { if (originNode.equals(destinationNode)) continue; ek.calculateMaximumFlow(originNode, destinationNode); k = Math.min(k, ek.getMaximumFlowValue().intValue()); if (k == 0) break; } } return k == Integer.MAX_VALUE ? 0 : k; } /** * Returns the set of nodes reachable from a given node. * * @param node Node * @return Collection of reachable nodes */ public Collection<Node> getNeighbors(Node node) { return Collections.unmodifiableCollection(getGraph_JUNG().getSuccessors(node)); } /** * <p>Returns the betweeness centrality of each node. The betweeness * centrality of a node is equal to the number of node-pair shortest paths which traverses * the node.</p> * * <p>Internally it makes use of the Brandes' algorithm.</p> * * @return Betweeness centrality of each node */ public DoubleMatrix1D getNodeBetweenessCentrality() { if (nodeBetweenessCentrality == null) computeBetweenessCentrality(); return nodeBetweenessCentrality; } /** * Returns the node connectivity. The node connectivity is equal to the smallest * number of node-disjoint paths between each node pair. * * <p>Internally it makes use of the (modified) Edmonds-Karp algorithm to compute the maximum * flow between each node pair, assuming a link capacity equal to one for every link.</p> * * @return Node connectivity */ public int getNodeConnectivity() { if (E == 0) return 0; DirectedGraph<Node, Link> graph = getGraph_JGraphT(); int k = Integer.MAX_VALUE; for (Node originNode : nodes) { for (Node destinationNode : nodes) { if (originNode.equals(destinationNode)) continue; DirectedGraph<Node, Link> auxGraph = (DirectedGraph<Node, Link>) JGraphTUtils .buildAuxiliaryNodeDisjointGraph(graph, originNode, destinationNode); EdmondsKarpMaximumFlow<Node, Node> ek = new EdmondsKarpMaximumFlow(auxGraph); ek.calculateMaximumFlow(originNode, destinationNode); k = Math.min(k, ek.getMaximumFlowValue().intValue()); if (k == 0) break; } } return k == Integer.MAX_VALUE ? 0 : k; } /** * Returns the number of outgoing links for each node. * * @return Number of outgoing links per node */ public DoubleMatrix1D getOutNodeDegree() { if (outNodeDegree == null) { Graph<Node, Link> aux_graph = getGraph_JUNG(); outNodeDegree = DoubleFactory1D.dense.make(N); for (Node node : nodes) { outNodeDegree.set(node.getIndex(), aux_graph.outDegree(node)); } } return outNodeDegree; } /** * Returns the spectral radius of the network. The spectral radius is equal * to the largest eigenvalue of the adjacency matrix. * * @return Spectral radius */ public double getSpectralRadius() { if (E == 0) return 0; double[] eig = getAdjacencyMatrixEigenvalues(); return eig[eig.length - 1]; } /** * Returns the symmetry ratio. The symmetry ratio is equal to the number * of distinct eigenvalues of the adjacency matrix divided by the network * density plus one. * * @return Symmetry ratio */ public double getSymmetryRatio() { if (E == 0) return 0; double[] eig = getAdjacencyMatrixEigenvalues(); return (double) DoubleUtils.unique(eig).length / (getDensity() + 1); } }