1 /*
2 * Copyright (c) 2003, the JUNG Project and the Regents of the University of
3 * California All rights reserved. This software is open-source under the BSD
4 * license; see either "license.txt" or http://jung.sourceforge.net/license.txt
5 * for a description.
6 */
7 package edu.uci.ics.jung.algorithms.matrix;
8
9 import java.util.ArrayList;
10 import java.util.Collection;
11 import java.util.List;
12 import java.util.Map;
13
14 import org.apache.commons.collections15.BidiMap;
15 import org.apache.commons.collections15.Factory;
16
17 import cern.colt.matrix.DoubleMatrix1D;
18 import cern.colt.matrix.DoubleMatrix2D;
19 import cern.colt.matrix.impl.DenseDoubleMatrix1D;
20 import cern.colt.matrix.impl.SparseDoubleMatrix2D;
21 import cern.colt.matrix.linalg.Algebra;
22 import edu.uci.ics.jung.algorithms.util.ConstantMap;
23 import edu.uci.ics.jung.algorithms.util.Indexer;
24 import edu.uci.ics.jung.graph.Graph;
25 import edu.uci.ics.jung.graph.UndirectedGraph;
26
27
28 /**
29 * Contains methods for performing the analogues of certain matrix operations on
30 * graphs.
31 * <p>
32 * These implementations are efficient on sparse graphs, but may not be the best
33 * implementations for very dense graphs.
34 *
35 * @author Joshua O'Madadhain
36 * @see MatrixElementOperations
37 */
38 public class GraphMatrixOperations
39 {
40 /**
41 * Returns the graph that corresponds to the square of the (weighted)
42 * adjacency matrix that the specified graph <code>g</code> encodes. The
43 * implementation of MatrixElementOperations that is furnished to the
44 * constructor specifies the implementation of the dot product, which is an
45 * integral part of matrix multiplication.
46 *
47 * @param g
48 * the graph to be squared
49 * @return the result of squaring g
50 */
51 @SuppressWarnings("unchecked")
52 public static <V,E> Graph<V,E> square(Graph<V,E> g,
53 Factory<E> edgeFactory, MatrixElementOperations<E> meo)
54 {
55 // create new graph of same type
56 Graph<V, E> squaredGraph = null;
57 try {
58 squaredGraph = g.getClass().newInstance();
59 } catch (InstantiationException e3) {
60 e3.printStackTrace();
61 } catch (IllegalAccessException e3) {
62 e3.printStackTrace();
63 }
64
65 Collection<V> vertices = g.getVertices();
66 for (V v : vertices)
67 {
68 squaredGraph.addVertex(v);
69 }
70 for (V v : vertices)
71 {
72 for (V src : g.getPredecessors(v))
73 {
74 // get the edge connecting src to v in G
75 E e1 = g.findEdge(src,v);
76 for (V dest : g.getSuccessors(v))
77 {
78 // get edge connecting v to dest in G
79 E e2 = g.findEdge(v,dest);
80 // collect data on path composed of e1 and e2
81 Number pathData = meo.computePathData(e1, e2);
82 E e = squaredGraph.findEdge(src,dest);
83 // if no edge from src to dest exists in G2, create one
84 if (e == null) {
85 e = edgeFactory.create();
86 squaredGraph.addEdge(e, src, dest);
87 }
88 meo.mergePaths(e, pathData);
89 }
90 }
91 }
92 return squaredGraph;
93 }
94
95 /**
96 * Creates a graph from a square (weighted) adjacency matrix. If
97 * <code>nev</code> is non-null then it will be used to store the edge weights.
98 *
99 * <p>Notes on implementation:
100 * <ul>
101 * <li>The matrix indices will be mapped onto vertices in the order in which the
102 * vertex factory generates the vertices. This means the user is responsible
103 * <li>The type of edges created (directed or undirected) depends
104 * entirely on the graph factory supplied, regardless of whether the
105 * matrix is symmetric or not. The Colt {@code Property.isSymmetric}
106 * method may be used to find out whether the matrix
107 * is symmetric prior to making this call.
108 * <li>The matrix supplied need not be square. If it is not square, then
109 * the
110 *
111 * @return a representation of <code>matrix</code> as a JUNG
112 * <code>Graph</code>
113 */
114 public static <V,E> Graph<V,E> matrixToGraph(DoubleMatrix2D matrix,
115 Factory<? extends Graph<V,E>> graphFactory,
116 Factory<V> vertexFactory, Factory<E> edgeFactory,
117 Map<E,Number> nev)
118 {
119 if (matrix.rows() != matrix.columns())
120 {
121 throw new IllegalArgumentException("Matrix must be square.");
122 }
123 int size = matrix.rows();
124
125 Graph<V,E> graph = graphFactory.create();
126
127 for(int i = 0; i < size; i++)
128 {
129 V vertex = vertexFactory.create();
130 graph.addVertex(vertex);
131 }
132
133 List<V> vertices = new ArrayList<V>(graph.getVertices());
134 for (int i = 0; i < size; i++)
135 {
136 for (int j = 0; j < size; j++)
137 {
138 double value = matrix.getQuick(i, j);
139 if (value != 0)
140 {
141 E e = edgeFactory.create();
142 if (graph.addEdge(e, vertices.get(i), vertices.get(j)))
143 {
144 if (e != null && nev != null)
145 nev.put(e, value);
146 }
147 }
148 }
149 }
150
151
152 return graph;
153 }
154
155
156 /**
157 * Creates a graph from a square (weighted) adjacency matrix.
158 *
159 * @return a representation of <code>matrix</code> as a JUNG <code>Graph</code>
160 */
161 public static <V,E> Graph<V,E> matrixToGraph(DoubleMatrix2D matrix,
162 Factory<? extends Graph<V,E>> graphFactory,
163 Factory<V> vertexFactory, Factory<E> edgeFactory)
164 {
165 return GraphMatrixOperations.<V,E>matrixToGraph(matrix,
166 graphFactory, vertexFactory, edgeFactory, null);
167 }
168
169 /**
170 * Returns an unweighted (0-1) adjacency matrix based on the specified graph.
171 * @param <V> the vertex type
172 * @param <E> the edge type
173 * @param g the graph to convert to a matrix
174 */
175 public static <V,E> SparseDoubleMatrix2D graphToSparseMatrix(Graph<V,E> g)
176 {
177 return graphToSparseMatrix(g, null);
178 }
179
180 /**
181 * Returns a SparseDoubleMatrix2D whose entries represent the edge weights for the
182 * edges in <code>g</code>, as specified by <code>nev</code>.
183 *
184 * <p>The <code>(i,j)</code> entry of the matrix returned will be equal to the sum
185 * of the weights of the edges connecting the vertex with index <code>i</code> to
186 * <code>j</code>.
187 *
188 * <p>If <code>nev</code> is <code>null</code>, then a constant edge weight of 1 is used.
189 *
190 * @param g
191 * @param nev
192 */
193 public static <V,E> SparseDoubleMatrix2D graphToSparseMatrix(Graph<V,E> g, Map<E,Number> nev)
194 {
195 if (nev == null)
196 nev = new ConstantMap<E,Number>(1);
197 int numVertices = g.getVertexCount();
198 SparseDoubleMatrix2D matrix = new SparseDoubleMatrix2D(numVertices,
199 numVertices);
200
201 BidiMap<V,Integer> indexer = Indexer.<V>create(g.getVertices());
202 int i=0;
203
204 for(V v : g.getVertices())
205 {
206 for (E e : g.getOutEdges(v))
207 {
208 V w = g.getOpposite(v,e);
209 int j = indexer.get(w);
210 matrix.set(i, j, matrix.getQuick(i,j) + nev.get(e).doubleValue());
211 }
212 i++;
213 }
214 return matrix;
215 }
216
217 /**
218 * Returns a diagonal matrix whose diagonal entries contain the degree for
219 * the corresponding node.
220 *
221 * <p>NOTE: the vertices will be traversed in the order given by the graph's vertex
222 * collection. If you want to be assured of a particular ordering, use a graph
223 * implementation that guarantees such an ordering (see the implementations with {@code Ordered}
224 * or {@code Sorted} in their name).
225 *
226 * @return SparseDoubleMatrix2D
227 */
228 public static <V,E> SparseDoubleMatrix2D createVertexDegreeDiagonalMatrix(Graph<V,E> graph)
229 {
230 int numVertices = graph.getVertexCount();
231 SparseDoubleMatrix2D matrix = new SparseDoubleMatrix2D(numVertices,
232 numVertices);
233 int i = 0;
234 for (V v : graph.getVertices())
235 {
236 matrix.set(i, i, graph.degree(v));
237 i++;
238 }
239 return matrix;
240 }
241
242 /**
243 * The idea here is based on the metaphor of an electric circuit. We assume
244 * that an undirected graph represents the structure of an electrical
245 * circuit where each edge has unit resistance. One unit of current is
246 * injected into any arbitrary vertex s and one unit of current is extracted
247 * from any arbitrary vertex t. The voltage at some vertex i for source
248 * vertex s and target vertex t can then be measured according to the
249 * equation: V_i^(s,t) = T_is - T-it where T is the voltage potential matrix
250 * returned by this method. *
251 *
252 * @param graph
253 * an undirected graph representing an electrical circuit
254 * @return the voltage potential matrix
255 * @see "P. Doyle and J. Snell, 'Random walks and electric networks,', 1989"
256 * @see "M. Newman, 'A measure of betweenness centrality based on random walks', pp. 5-7, 2003"
257 */
258 public static <V,E> DoubleMatrix2D computeVoltagePotentialMatrix(
259 UndirectedGraph<V,E> graph)
260 {
261 int numVertices = graph.getVertexCount();
262 //create adjacency matrix from graph
263 DoubleMatrix2D A = GraphMatrixOperations.graphToSparseMatrix(graph,
264 null);
265 //create diagonal matrix of vertex degrees
266 DoubleMatrix2D D = GraphMatrixOperations
267 .createVertexDegreeDiagonalMatrix(graph);
268 DoubleMatrix2D temp = new SparseDoubleMatrix2D(numVertices - 1,
269 numVertices - 1);
270 //compute D - A except for last row and column
271 for (int i = 0; i < numVertices - 1; i++)
272 {
273 for (int j = 0; j < numVertices - 1; j++)
274 {
275 temp.set(i, j, D.get(i, j) - A.get(i, j));
276 }
277 }
278 Algebra algebra = new Algebra();
279 DoubleMatrix2D tempInverse = algebra.inverse(temp);
280 DoubleMatrix2D T = new SparseDoubleMatrix2D(numVertices, numVertices);
281 //compute "voltage" matrix
282 for (int i = 0; i < numVertices - 1; i++)
283 {
284 for (int j = 0; j < numVertices - 1; j++)
285 {
286 T.set(i, j, tempInverse.get(i, j));
287 }
288 }
289 return T;
290 }
291
292 /**
293 * Converts a Map of (Vertex, Double) pairs to a DoubleMatrix1D.
294 *
295 * <p>Note: the vertices will appear in the output array in the order given
296 * by {@code map}'s iterator. If you want a particular ordering, use a {@code Map}
297 * implementation that provides such an ordering ({@code SortedMap, LinkedHashMap}, etc.).
298 */
299 public static <V> DoubleMatrix1D mapTo1DMatrix(Map<V,Number> map)
300 {
301 int numVertices = map.size();
302 DoubleMatrix1D vector = new DenseDoubleMatrix1D(numVertices);
303 int i = 0;
304 for (V v : map.keySet())
305 {
306 vector.set(i, map.get(v).doubleValue());
307 i++;
308 }
309 return vector;
310 }
311
312 /**
313 * Computes the all-pairs mean first passage time for the specified graph,
314 * given an existing stationary probability distribution.
315 * <P>
316 * The mean first passage time from vertex v to vertex w is defined, for a
317 * Markov network (in which the vertices represent states and the edge
318 * weights represent state->state transition probabilities), as the expected
319 * number of steps required to travel from v to w if the steps occur
320 * according to the transition probabilities.
321 * <P>
322 * The stationary distribution is the fraction of time, in the limit as the
323 * number of state transitions approaches infinity, that a given state will
324 * have been visited. Equivalently, it is the probability that a given state
325 * will be the current state after an arbitrarily large number of state
326 * transitions.
327 *
328 * @param G
329 * the graph on which the MFPT will be calculated
330 * @param edgeWeights
331 * the edge weights
332 * @param stationaryDistribution
333 * the asymptotic state probabilities
334 * @return the mean first passage time matrix
335 */
336 public static <V,E> DoubleMatrix2D computeMeanFirstPassageMatrix(Graph<V,E> G,
337 Map<E,Number> edgeWeights, DoubleMatrix1D stationaryDistribution)
338 {
339 DoubleMatrix2D temp = GraphMatrixOperations.graphToSparseMatrix(G,
340 edgeWeights);
341 for (int i = 0; i < temp.rows(); i++)
342 {
343 for (int j = 0; j < temp.columns(); j++)
344 {
345 double value = -1 * temp.get(i, j)
346 + stationaryDistribution.get(j);
347 if (i == j)
348 value += 1;
349 if (value != 0)
350 temp.set(i, j, value);
351 }
352 }
353 Algebra algebra = new Algebra();
354 DoubleMatrix2D fundamentalMatrix = algebra.inverse(temp);
355 temp = new SparseDoubleMatrix2D(temp.rows(), temp.columns());
356 for (int i = 0; i < temp.rows(); i++)
357 {
358 for (int j = 0; j < temp.columns(); j++)
359 {
360 double value = -1.0 * fundamentalMatrix.get(i, j);
361 value += fundamentalMatrix.get(j, j);
362 if (i == j)
363 value += 1;
364 if (value != 0)
365 temp.set(i, j, value);
366 }
367 }
368 DoubleMatrix2D stationaryMatrixDiagonal = new SparseDoubleMatrix2D(temp
369 .rows(), temp.columns());
370 int numVertices = stationaryDistribution.size();
371 for (int i = 0; i < numVertices; i++)
372 stationaryMatrixDiagonal.set(i, i, 1.0 / stationaryDistribution
373 .get(i));
374 DoubleMatrix2D meanFirstPassageMatrix = algebra.mult(temp,
375 stationaryMatrixDiagonal);
376 return meanFirstPassageMatrix;
377 }
378 }