1 /* -*- mode: C -*- */
2 /* vim:set ts=4 sw=4 sts=4 et: */
3 /*
4 IGraph library.
5 Copyright (C) 2006-2012 Gabor Csardi <csardi.gabor@gmail.com>
6 334 Harvard street, Cambridge, MA 02139 USA
7
8 This program is free software; you can redistribute it and/or modify
9 it under the terms of the GNU General Public License as published by
10 the Free Software Foundation; either version 2 of the License, or
11 (at your option) any later version.
12
13 This program is distributed in the hope that it will be useful,
14 but WITHOUT ANY WARRANTY; without even the implied warranty of
15 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 GNU General Public License for more details.
17
18 You should have received a copy of the GNU General Public License
19 along with this program; if not, write to the Free Software
20 Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21 02110-1301 USA
22
23 */
24
25 #include "igraph_structural.h"
26 #include "igraph_interface.h"
27
28 #include <math.h>
29
igraph_i_weighted_laplacian(const igraph_t * graph,igraph_matrix_t * res,igraph_sparsemat_t * sparseres,igraph_bool_t normalized,const igraph_vector_t * weights)30 static int igraph_i_weighted_laplacian(const igraph_t *graph, igraph_matrix_t *res,
31 igraph_sparsemat_t *sparseres,
32 igraph_bool_t normalized,
33 const igraph_vector_t *weights) {
34
35 igraph_eit_t edgeit;
36 int no_of_nodes = (int) igraph_vcount(graph);
37 int no_of_edges = (int) igraph_ecount(graph);
38 igraph_bool_t directed = igraph_is_directed(graph);
39 igraph_vector_t degree;
40 long int i;
41
42 if (igraph_vector_size(weights) != no_of_edges) {
43 IGRAPH_ERROR("Invalid edge weight vector length", IGRAPH_EINVAL);
44 }
45
46 if (res) {
47 IGRAPH_CHECK(igraph_matrix_resize(res, no_of_nodes, no_of_nodes));
48 igraph_matrix_null(res);
49 }
50 if (sparseres) {
51 int nz = directed ? no_of_edges + no_of_nodes :
52 no_of_edges * 2 + no_of_nodes;
53 igraph_sparsemat_resize(sparseres, no_of_nodes, no_of_nodes, nz);
54 }
55
56 IGRAPH_CHECK(igraph_eit_create(graph, igraph_ess_all(0), &edgeit));
57 IGRAPH_FINALLY(igraph_eit_destroy, &edgeit);
58
59 IGRAPH_VECTOR_INIT_FINALLY(°ree, no_of_nodes);
60
61 if (directed) {
62
63 if (!normalized) {
64
65 while (!IGRAPH_EIT_END(edgeit)) {
66 long int edge = IGRAPH_EIT_GET(edgeit);
67 long int from = IGRAPH_FROM(graph, edge);
68 long int to = IGRAPH_TO (graph, edge);
69 igraph_real_t weight = VECTOR(*weights)[edge];
70 if (from != to) {
71 if (res) {
72 MATRIX(*res, from, to) -= weight;
73 }
74 if (sparseres) {
75 IGRAPH_CHECK(igraph_sparsemat_entry(sparseres, (int) from, (int)to,
76 -weight));
77 }
78 VECTOR(degree)[from] += weight;
79 }
80 IGRAPH_EIT_NEXT(edgeit);
81 }
82
83 /* And the diagonal */
84 for (i = 0; i < no_of_nodes; i++) {
85 if (res) {
86 MATRIX(*res, i, i) = VECTOR(degree)[i];
87 }
88 if (sparseres) {
89 IGRAPH_CHECK(igraph_sparsemat_entry(sparseres, (int) i, (int) i,
90 VECTOR(degree)[i]));
91 }
92 }
93
94 } else { /* normalized */
95
96 while (!IGRAPH_EIT_END(edgeit)) {
97 long int edge = IGRAPH_EIT_GET(edgeit);
98 long int from = IGRAPH_FROM(graph, edge);
99 long int to = IGRAPH_TO (graph, edge);
100 igraph_real_t weight = VECTOR(*weights)[edge];
101 if (from != to) {
102 VECTOR(degree)[from] += weight;
103 }
104 IGRAPH_EIT_NEXT(edgeit);
105 }
106
107 for (i = 0; i < no_of_nodes; i++) {
108 int t = VECTOR(degree)[i] > 0 ? 1 : 0;
109 if (res) {
110 MATRIX(*res, i, i) = t;
111 }
112 if (sparseres) {
113 IGRAPH_CHECK(igraph_sparsemat_entry(sparseres, (int) i, (int) i, t));
114 }
115 }
116
117 IGRAPH_EIT_RESET(edgeit);
118 while (!IGRAPH_EIT_END(edgeit)) {
119 long int edge = IGRAPH_EIT_GET(edgeit);
120 long int from = IGRAPH_FROM(graph, edge);
121 long int to = IGRAPH_TO (graph, edge);
122 igraph_real_t weight = VECTOR(*weights)[edge];
123 if (from != to) {
124 igraph_real_t t = weight / VECTOR(degree)[from];
125 if (res) {
126 MATRIX(*res, from, to) -= t;
127 }
128 if (sparseres) {
129 IGRAPH_CHECK(igraph_sparsemat_entry(sparseres, (int) from, (int) to,
130 -t));
131 }
132 }
133 IGRAPH_EIT_NEXT(edgeit);
134 }
135
136 }
137
138 } else { /* undirected */
139
140 if (!normalized) {
141
142 while (!IGRAPH_EIT_END(edgeit)) {
143 long int edge = IGRAPH_EIT_GET(edgeit);
144 long int from = IGRAPH_FROM(graph, edge);
145 long int to = IGRAPH_TO (graph, edge);
146 igraph_real_t weight = VECTOR(*weights)[edge];
147 if (from != to) {
148 if (res) {
149 MATRIX(*res, from, to) -= weight;
150 MATRIX(*res, to, from) -= weight;
151 }
152 if (sparseres) {
153 IGRAPH_CHECK(igraph_sparsemat_entry(sparseres, (int) from, (int) to,
154 -weight));
155 IGRAPH_CHECK(igraph_sparsemat_entry(sparseres, (int) to, (int) from,
156 -weight));
157 }
158 VECTOR(degree)[from] += weight;
159 VECTOR(degree)[to] += weight;
160 }
161 IGRAPH_EIT_NEXT(edgeit);
162 }
163
164 /* And the diagonal */
165 for (i = 0; i < no_of_nodes; i++) {
166 if (res) {
167 MATRIX(*res, i, i) = VECTOR(degree)[i];
168 }
169 if (sparseres) {
170 IGRAPH_CHECK(igraph_sparsemat_entry(sparseres, (int) i, (int) i,
171 VECTOR(degree)[i]));
172 }
173 }
174
175 } else { /* normalized */
176
177 while (!IGRAPH_EIT_END(edgeit)) {
178 long int edge = IGRAPH_EIT_GET(edgeit);
179 long int from = IGRAPH_FROM(graph, edge);
180 long int to = IGRAPH_TO (graph, edge);
181 igraph_real_t weight = VECTOR(*weights)[edge];
182 if (from != to) {
183 VECTOR(degree)[from] += weight;
184 VECTOR(degree)[to] += weight;
185 }
186 IGRAPH_EIT_NEXT(edgeit);
187 }
188
189 for (i = 0; i < no_of_nodes; i++) {
190 int t = VECTOR(degree)[i] > 0 ? 1 : 0;
191 if (res) {
192 MATRIX(*res, i, i) = t;
193 }
194 if (sparseres) {
195 IGRAPH_CHECK(igraph_sparsemat_entry(sparseres, (int) i, (int) i, t));
196 }
197 VECTOR(degree)[i] = sqrt(VECTOR(degree)[i]);
198 }
199
200 IGRAPH_EIT_RESET(edgeit);
201 while (!IGRAPH_EIT_END(edgeit)) {
202 long int edge = IGRAPH_EIT_GET(edgeit);
203 long int from = IGRAPH_FROM(graph, edge);
204 long int to = IGRAPH_TO (graph, edge);
205 igraph_real_t weight = VECTOR(*weights)[edge];
206 if (from != to) {
207 double diff = weight / (VECTOR(degree)[from] * VECTOR(degree)[to]);
208 if (res) {
209 MATRIX(*res, from, to) -= diff;
210 MATRIX(*res, to, from) -= diff;
211 }
212 if (sparseres) {
213 IGRAPH_CHECK(igraph_sparsemat_entry(sparseres, (int) from, (int) to,
214 -diff));
215 IGRAPH_CHECK(igraph_sparsemat_entry(sparseres, (int) to, (int) from,
216 -diff));
217 }
218 }
219 IGRAPH_EIT_NEXT(edgeit);
220 }
221
222 }
223
224 }
225
226 igraph_vector_destroy(°ree);
227 igraph_eit_destroy(&edgeit);
228 IGRAPH_FINALLY_CLEAN(2);
229
230 return 0;
231 }
232
233 /**
234 * \function igraph_laplacian
235 * \brief Returns the Laplacian matrix of a graph
236 *
237 * </para><para>
238 * The graph Laplacian matrix is similar to an adjacency matrix but
239 * contains -1's instead of 1's and the vertex degrees are included in
240 * the diagonal. So the result for edge i--j is -1 if i!=j and is equal
241 * to the degree of vertex i if i==j. igraph_laplacian will work on a
242 * directed graph; in this case, the diagonal will contain the out-degrees.
243 * Loop edges will be ignored.
244 *
245 * </para><para>
246 * The normalized version of the Laplacian matrix has 1 in the diagonal and
247 * -1/sqrt(d[i]d[j]) if there is an edge from i to j.
248 *
249 * </para><para>
250 * The first version of this function was written by Vincent Matossian.
251 * \param graph Pointer to the graph to convert.
252 * \param res Pointer to an initialized matrix object, the result is
253 * stored here. It will be resized if needed.
254 * If it is a null pointer, then it is ignored.
255 * At least one of \p res and \p sparseres must be a non-null pointer.
256 * \param sparseres Pointer to an initialized sparse matrix object, the
257 * result is stored here, if it is not a null pointer.
258 * At least one of \p res and \p sparseres must be a non-null pointer.
259 * \param normalized Whether to create a normalized Laplacian matrix.
260 * \param weights An optional vector containing edge weights, to calculate
261 * the weighted Laplacian matrix. Set it to a null pointer to
262 * calculate the unweighted Laplacian.
263 * \return Error code.
264 *
265 * Time complexity: O(|V||V|),
266 * |V| is the
267 * number of vertices in the graph.
268 *
269 * \example examples/simple/igraph_laplacian.c
270 */
271
igraph_laplacian(const igraph_t * graph,igraph_matrix_t * res,igraph_sparsemat_t * sparseres,igraph_bool_t normalized,const igraph_vector_t * weights)272 int igraph_laplacian(const igraph_t *graph, igraph_matrix_t *res,
273 igraph_sparsemat_t *sparseres,
274 igraph_bool_t normalized,
275 const igraph_vector_t *weights) {
276
277 igraph_eit_t edgeit;
278 int no_of_nodes = (int) igraph_vcount(graph);
279 int no_of_edges = (int) igraph_ecount(graph);
280 igraph_bool_t directed = igraph_is_directed(graph);
281 int from, to;
282 igraph_integer_t ffrom, fto;
283 igraph_vector_t degree;
284 int i;
285
286 if (!res && !sparseres) {
287 IGRAPH_ERROR("Laplacian: give at least one of `res' or `sparseres'",
288 IGRAPH_EINVAL);
289 }
290
291 if (weights) {
292 return igraph_i_weighted_laplacian(graph, res, sparseres, normalized,
293 weights);
294 }
295
296 if (res) {
297 IGRAPH_CHECK(igraph_matrix_resize(res, no_of_nodes, no_of_nodes));
298 igraph_matrix_null(res);
299 }
300 if (sparseres) {
301 int nz = directed ? no_of_edges + no_of_nodes :
302 no_of_edges * 2 + no_of_nodes;
303 IGRAPH_CHECK(igraph_sparsemat_resize(sparseres, no_of_nodes,
304 no_of_nodes, nz));
305 }
306 IGRAPH_CHECK(igraph_eit_create(graph, igraph_ess_all(0), &edgeit));
307 IGRAPH_FINALLY(igraph_eit_destroy, &edgeit);
308
309 IGRAPH_VECTOR_INIT_FINALLY(°ree, no_of_nodes);
310
311 IGRAPH_CHECK(igraph_degree(graph, °ree, igraph_vss_all(),
312 IGRAPH_OUT, IGRAPH_NO_LOOPS));
313
314 if (directed) {
315 if (!normalized) {
316 for (i = 0; i < no_of_nodes; i++) {
317 if (res) {
318 MATRIX(*res, i, i) = VECTOR(degree)[i];
319 }
320 if (sparseres) {
321 IGRAPH_CHECK(igraph_sparsemat_entry(sparseres, i, i,
322 VECTOR(degree)[i]));
323 }
324 }
325 while (!IGRAPH_EIT_END(edgeit)) {
326 igraph_edge(graph, IGRAPH_EIT_GET(edgeit), &ffrom, &fto);
327 from = ffrom;
328 to = fto;
329 if (from != to) {
330 if (res) {
331 MATRIX(*res, from, to) -= 1;
332 }
333 if (sparseres) {
334 IGRAPH_CHECK(igraph_sparsemat_entry(sparseres, from, to, -1.0));
335 }
336 }
337 IGRAPH_EIT_NEXT(edgeit);
338 }
339 } else {
340 for (i = 0; i < no_of_nodes; i++) {
341 int t = VECTOR(degree)[i] > 0 ? 1 : 0;
342 if (res) {
343 MATRIX(*res, i, i) = t;
344 }
345 if (sparseres) {
346 IGRAPH_CHECK(igraph_sparsemat_entry(sparseres, i, i, t));
347 }
348 if (VECTOR(degree)[i] > 0) {
349 VECTOR(degree)[i] = 1.0 / VECTOR(degree)[i];
350 }
351 }
352
353 while (!IGRAPH_EIT_END(edgeit)) {
354 igraph_edge(graph, IGRAPH_EIT_GET(edgeit), &ffrom, &fto);
355 from = ffrom; to = fto;
356 if (from != to) {
357 if (res) {
358 MATRIX(*res, from, to) -= VECTOR(degree)[from];
359 }
360 if (sparseres) {
361 IGRAPH_CHECK(igraph_sparsemat_entry(sparseres, from, to,
362 -VECTOR(degree)[from]));
363 }
364 }
365 IGRAPH_EIT_NEXT(edgeit);
366 }
367 }
368
369 } else {
370
371 if (!normalized) {
372 for (i = 0; i < no_of_nodes; i++) {
373 if (res) {
374 MATRIX(*res, i, i) = VECTOR(degree)[i];
375 }
376 if (sparseres) {
377 IGRAPH_CHECK(igraph_sparsemat_entry(sparseres, i, i,
378 VECTOR(degree)[i]));
379 }
380 }
381
382 while (!IGRAPH_EIT_END(edgeit)) {
383 igraph_edge(graph, IGRAPH_EIT_GET(edgeit), &ffrom, &fto);
384 from = ffrom;
385 to = fto;
386
387 if (from != to) {
388 if (res) {
389 MATRIX(*res, to, from) -= 1;
390 MATRIX(*res, from, to) -= 1;
391 }
392 if (sparseres) {
393 IGRAPH_CHECK(igraph_sparsemat_entry(sparseres, to, from, -1.0));
394 IGRAPH_CHECK(igraph_sparsemat_entry(sparseres, from, to, -1.0));
395 }
396 }
397
398 IGRAPH_EIT_NEXT(edgeit);
399 }
400 } else {
401 for (i = 0; i < no_of_nodes; i++) {
402 int t = VECTOR(degree)[i] > 0 ? 1 : 0;
403 if (res) {
404 MATRIX(*res, i, i) = t;
405 }
406 if (sparseres) {
407 IGRAPH_CHECK(igraph_sparsemat_entry(sparseres, i, i, t));
408 }
409 VECTOR(degree)[i] = sqrt(VECTOR(degree)[i]);
410 }
411
412 while (!IGRAPH_EIT_END(edgeit)) {
413 igraph_edge(graph, IGRAPH_EIT_GET(edgeit), &ffrom, &fto);
414 from = ffrom; to = fto;
415 if (from != to) {
416 double diff = 1.0 / (VECTOR(degree)[from] * VECTOR(degree)[to]);
417 if (res) {
418 MATRIX(*res, from, to) -= diff;
419 MATRIX(*res, to, from) -= diff;
420 }
421 if (sparseres) {
422 IGRAPH_CHECK(igraph_sparsemat_entry(sparseres, from, to, -diff));
423 IGRAPH_CHECK(igraph_sparsemat_entry(sparseres, to, from, -diff));
424 }
425 }
426 IGRAPH_EIT_NEXT(edgeit);
427 }
428 }
429
430 }
431
432 igraph_vector_destroy(°ree);
433 igraph_eit_destroy(&edgeit);
434 IGRAPH_FINALLY_CLEAN(2);
435 return 0;
436 }
437