1 /* -*- mode: C -*- */
2 /* vim:set ts=4 sw=4 sts=4 et: */
3 /*
4 IGraph library.
5 Copyright (C) 2003-2020 The igraph development team
6
7 This program is free software; you can redistribute it and/or modify
8 it under the terms of the GNU General Public License as published by
9 the Free Software Foundation; either version 2 of the License, or
10 (at your option) any later version.
11
12 This program is distributed in the hope that it will be useful,
13 but WITHOUT ANY WARRANTY; without even the implied warranty of
14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 GNU General Public License for more details.
16
17 You should have received a copy of the GNU General Public License
18 along with this program; if not, write to the Free Software
19 Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
20 02110-1301 USA
21
22 */
23
24 #include "igraph_layout.h"
25
26 #include "igraph_blas.h"
27 #include "igraph_components.h"
28 #include "igraph_eigen.h"
29 #include "igraph_interface.h"
30 #include "igraph_memory.h"
31 #include "igraph_operators.h"
32 #include "igraph_paths.h"
33 #include "igraph_random.h"
34 #include "igraph_structural.h"
35
36 static int igraph_i_layout_mds_step(igraph_real_t *to, const igraph_real_t *from,
37 int n, void *extra);
38
39 static int igraph_i_layout_mds_single(const igraph_t* graph, igraph_matrix_t *res,
40 igraph_matrix_t *dist, long int dim);
41
igraph_i_layout_mds_step(igraph_real_t * to,const igraph_real_t * from,int n,void * extra)42 static int igraph_i_layout_mds_step(igraph_real_t *to, const igraph_real_t *from,
43 int n, void *extra) {
44 igraph_matrix_t* matrix = (igraph_matrix_t*)extra;
45 IGRAPH_UNUSED(n);
46 igraph_blas_dgemv_array(0, 1, matrix, from, 0, to);
47 return 0;
48 }
49
50 /* MDS layout for a connected graph, with no error checking on the
51 * input parameters. The distance matrix will be modified in-place. */
igraph_i_layout_mds_single(const igraph_t * graph,igraph_matrix_t * res,igraph_matrix_t * dist,long int dim)52 int igraph_i_layout_mds_single(const igraph_t* graph, igraph_matrix_t *res,
53 igraph_matrix_t *dist, long int dim) {
54
55 long int no_of_nodes = igraph_vcount(graph);
56 long int nev = dim;
57 igraph_matrix_t vectors;
58 igraph_vector_t values, row_means;
59 igraph_real_t grand_mean;
60 long int i, j, k;
61 igraph_eigen_which_t which;
62
63 /* Handle the trivial cases */
64 if (no_of_nodes == 1) {
65 IGRAPH_CHECK(igraph_matrix_resize(res, 1, dim));
66 igraph_matrix_fill(res, 0);
67 return IGRAPH_SUCCESS;
68 }
69 if (no_of_nodes == 2) {
70 IGRAPH_CHECK(igraph_matrix_resize(res, 2, dim));
71 igraph_matrix_fill(res, 0);
72 for (j = 0; j < dim; j++) {
73 MATRIX(*res, 1, j) = 1;
74 }
75 return IGRAPH_SUCCESS;
76 }
77
78 /* Initialize some stuff */
79 IGRAPH_VECTOR_INIT_FINALLY(&values, no_of_nodes);
80 IGRAPH_CHECK(igraph_matrix_init(&vectors, no_of_nodes, dim));
81 IGRAPH_FINALLY(igraph_matrix_destroy, &vectors);
82
83 /* Take the square of the distance matrix */
84 for (i = 0; i < no_of_nodes; i++) {
85 for (j = 0; j < no_of_nodes; j++) {
86 MATRIX(*dist, i, j) *= MATRIX(*dist, i, j);
87 }
88 }
89
90 /* Double centering of the distance matrix */
91 IGRAPH_VECTOR_INIT_FINALLY(&row_means, no_of_nodes);
92 igraph_vector_fill(&values, 1.0 / no_of_nodes);
93 igraph_blas_dgemv(0, 1, dist, &values, 0, &row_means);
94 grand_mean = igraph_vector_sum(&row_means) / no_of_nodes;
95 igraph_matrix_add_constant(dist, grand_mean);
96 for (i = 0; i < no_of_nodes; i++) {
97 for (j = 0; j < no_of_nodes; j++) {
98 MATRIX(*dist, i, j) -= VECTOR(row_means)[i] + VECTOR(row_means)[j];
99 MATRIX(*dist, i, j) *= -0.5;
100 }
101 }
102 igraph_vector_destroy(&row_means);
103 IGRAPH_FINALLY_CLEAN(1);
104
105 /* Calculate the top `dim` eigenvectors. */
106 which.pos = IGRAPH_EIGEN_LA;
107 which.howmany = (int) nev;
108 IGRAPH_CHECK(igraph_eigen_matrix_symmetric(/*A=*/ 0, /*sA=*/ 0,
109 /*fun=*/ igraph_i_layout_mds_step,
110 /*n=*/ (int) no_of_nodes, /*extra=*/ dist,
111 /*algorithm=*/ IGRAPH_EIGEN_LAPACK,
112 &which, /*options=*/ 0, /*storage=*/ 0,
113 &values, &vectors));
114
115 /* Calculate and normalize the final coordinates */
116 for (j = 0; j < nev; j++) {
117 VECTOR(values)[j] = sqrt(fabs(VECTOR(values)[j]));
118 }
119 IGRAPH_CHECK(igraph_matrix_resize(res, no_of_nodes, dim));
120 for (i = 0; i < no_of_nodes; i++) {
121 for (j = 0, k = nev - 1; j < nev; j++, k--) {
122 MATRIX(*res, i, k) = VECTOR(values)[j] * MATRIX(vectors, i, j);
123 }
124 }
125
126 igraph_matrix_destroy(&vectors);
127 igraph_vector_destroy(&values);
128 IGRAPH_FINALLY_CLEAN(2);
129
130 return IGRAPH_SUCCESS;
131 }
132
133 /**
134 * \function igraph_layout_mds
135 * \brief Place the vertices on a plane using multidimensional scaling.
136 *
137 * </para><para>
138 * This layout requires a distance matrix, where the intersection of
139 * row i and column j specifies the desired distance between vertex i
140 * and vertex j. The algorithm will try to place the vertices in a
141 * space having a given number of dimensions in a way that approximates
142 * the distance relations prescribed in the distance matrix. igraph
143 * uses the classical multidimensional scaling by Torgerson; for more
144 * details, see Cox & Cox: Multidimensional Scaling (1994), Chapman
145 * and Hall, London.
146 *
147 * </para><para>
148 * If the input graph is disconnected, igraph will decompose it
149 * first into its subgraphs, lay out the subgraphs one by one
150 * using the appropriate submatrices of the distance matrix, and
151 * then merge the layouts using \ref igraph_layout_merge_dla.
152 * Since \ref igraph_layout_merge_dla works for 2D layouts only,
153 * you cannot run the MDS layout on disconnected graphs for
154 * more than two dimensions.
155 *
156 * </para><para>
157 * Warning: if the graph is symmetric to the exchange of two vertices
158 * (as is the case with leaves of a tree connecting to the same parent),
159 * classical multidimensional scaling may assign the same coordinates to
160 * these vertices.
161 *
162 * \param graph A graph object.
163 * \param res Pointer to an initialized matrix object. This will
164 * contain the result and will be resized if needed.
165 * \param dist The distance matrix. It must be symmetric and this
166 * function does not check whether the matrix is indeed
167 * symmetric. Results are unspecified if you pass a non-symmetric
168 * matrix here. You can set this parameter to null; in this
169 * case, the shortest path lengths between vertices will be
170 * used as distances.
171 * \param dim The number of dimensions in the embedding space. For
172 * 2D layouts, supply 2 here.
173 * \return Error code.
174 *
175 * Added in version 0.6.
176 *
177 * </para><para>
178 * Time complexity: usually around O(|V|^2 dim).
179 */
180
igraph_layout_mds(const igraph_t * graph,igraph_matrix_t * res,const igraph_matrix_t * dist,long int dim)181 int igraph_layout_mds(const igraph_t* graph, igraph_matrix_t *res,
182 const igraph_matrix_t *dist, long int dim) {
183 long int i, no_of_nodes = igraph_vcount(graph);
184 igraph_matrix_t m;
185 igraph_bool_t conn;
186
187 RNG_BEGIN();
188
189 /* Check the distance matrix */
190 if (dist && (igraph_matrix_nrow(dist) != no_of_nodes ||
191 igraph_matrix_ncol(dist) != no_of_nodes)) {
192 IGRAPH_ERROR("invalid distance matrix size", IGRAPH_EINVAL);
193 }
194
195 /* Check the number of dimensions */
196 if (dim <= 1) {
197 IGRAPH_ERROR("dim must be positive", IGRAPH_EINVAL);
198 }
199 if (dim > no_of_nodes) {
200 IGRAPH_ERROR("dim must be less than the number of nodes", IGRAPH_EINVAL);
201 }
202
203 /* Copy or obtain the distance matrix */
204 if (dist == 0) {
205 IGRAPH_CHECK(igraph_matrix_init(&m, no_of_nodes, no_of_nodes));
206 IGRAPH_FINALLY(igraph_matrix_destroy, &m);
207 IGRAPH_CHECK(igraph_shortest_paths(graph, &m,
208 igraph_vss_all(), igraph_vss_all(), IGRAPH_ALL));
209 } else {
210 IGRAPH_CHECK(igraph_matrix_copy(&m, dist));
211 IGRAPH_FINALLY(igraph_matrix_destroy, &m);
212 /* Make sure that the diagonal contains zeroes only */
213 for (i = 0; i < no_of_nodes; i++) {
214 MATRIX(m, i, i) = 0.0;
215 }
216 }
217
218 /* Check whether the graph is connected */
219 IGRAPH_CHECK(igraph_is_connected(graph, &conn, IGRAPH_WEAK));
220 if (conn) {
221 /* Yes, it is, just do the MDS */
222 IGRAPH_CHECK(igraph_i_layout_mds_single(graph, res, &m, dim));
223 } else {
224 /* The graph is not connected, lay out the components one by one */
225 igraph_vector_ptr_t layouts;
226 igraph_vector_t comp, vertex_order;
227 igraph_t subgraph;
228 igraph_matrix_t *layout;
229 igraph_matrix_t dist_submatrix;
230 igraph_bool_t *seen_vertices;
231 long int j, n, processed_vertex_count = 0;
232
233 IGRAPH_VECTOR_INIT_FINALLY(&comp, 0);
234 IGRAPH_VECTOR_INIT_FINALLY(&vertex_order, no_of_nodes);
235
236 IGRAPH_CHECK(igraph_vector_ptr_init(&layouts, 0));
237 IGRAPH_FINALLY(igraph_vector_ptr_destroy_all, &layouts);
238 igraph_vector_ptr_set_item_destructor(&layouts, (igraph_finally_func_t*)igraph_matrix_destroy);
239
240 IGRAPH_CHECK(igraph_matrix_init(&dist_submatrix, 0, 0));
241 IGRAPH_FINALLY(igraph_matrix_destroy, &dist_submatrix);
242
243 seen_vertices = IGRAPH_CALLOC(no_of_nodes, igraph_bool_t);
244 if (seen_vertices == 0) {
245 IGRAPH_ERROR("cannot calculate MDS layout", IGRAPH_ENOMEM);
246 }
247 IGRAPH_FINALLY(igraph_free, seen_vertices);
248
249 for (i = 0; i < no_of_nodes; i++) {
250 if (seen_vertices[i]) {
251 continue;
252 }
253
254 /* This is a vertex whose component we did not lay out so far */
255 IGRAPH_CHECK(igraph_subcomponent(graph, &comp, i, IGRAPH_ALL));
256 /* Take the subgraph */
257 IGRAPH_CHECK(igraph_induced_subgraph(graph, &subgraph, igraph_vss_vector(&comp),
258 IGRAPH_SUBGRAPH_AUTO));
259 IGRAPH_FINALLY(igraph_destroy, &subgraph);
260 /* Calculate the submatrix of the distances */
261 IGRAPH_CHECK(igraph_matrix_select_rows_cols(&m, &dist_submatrix,
262 &comp, &comp));
263 /* Allocate a new matrix for storing the layout */
264 layout = IGRAPH_CALLOC(1, igraph_matrix_t);
265 if (layout == 0) {
266 IGRAPH_ERROR("cannot calculate MDS layout", IGRAPH_ENOMEM);
267 }
268 IGRAPH_FINALLY(igraph_free, layout);
269 IGRAPH_CHECK(igraph_matrix_init(layout, 0, 0));
270 IGRAPH_FINALLY(igraph_matrix_destroy, layout);
271 /* Lay out the subgraph */
272 IGRAPH_CHECK(igraph_i_layout_mds_single(&subgraph, layout, &dist_submatrix, dim));
273 /* Store the layout */
274 IGRAPH_CHECK(igraph_vector_ptr_push_back(&layouts, layout));
275 IGRAPH_FINALLY_CLEAN(2); /* ownership of layout taken by layouts */
276 /* Free the newly created subgraph */
277 igraph_destroy(&subgraph);
278 IGRAPH_FINALLY_CLEAN(1);
279 /* Mark all the vertices in the component as visited */
280 n = igraph_vector_size(&comp);
281 for (j = 0; j < n; j++) {
282 seen_vertices[(long int)VECTOR(comp)[j]] = 1;
283 VECTOR(vertex_order)[(long int)VECTOR(comp)[j]] = processed_vertex_count++;
284 }
285 }
286 /* Merge the layouts - reusing dist_submatrix here */
287 IGRAPH_CHECK(igraph_layout_merge_dla(0, &layouts, &dist_submatrix));
288 /* Reordering the rows of res to match the original graph */
289 IGRAPH_CHECK(igraph_matrix_select_rows(&dist_submatrix, res, &vertex_order));
290
291 igraph_free(seen_vertices);
292 igraph_matrix_destroy(&dist_submatrix);
293 igraph_vector_ptr_destroy_all(&layouts);
294 igraph_vector_destroy(&vertex_order);
295 igraph_vector_destroy(&comp);
296 IGRAPH_FINALLY_CLEAN(5);
297 }
298
299 RNG_END();
300
301 igraph_matrix_destroy(&m);
302 IGRAPH_FINALLY_CLEAN(1);
303
304 return IGRAPH_SUCCESS;
305 }
306