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_interface.h"
27 #include "igraph_random.h"
28
29 #include "core/math.h"
30 #include "core/interruption.h"
31
32 /**
33 * \ingroup layout
34 * \function igraph_layout_gem
35 *
36 * The GEM layout algorithm, as described in Arne Frick, Andreas Ludwig,
37 * Heiko Mehldau: A Fast Adaptive Layout Algorithm for Undirected Graphs,
38 * Proc. Graph Drawing 1994, LNCS 894, pp. 388-403, 1995.
39 * \param graph The input graph. Edge directions are ignored in
40 * directed graphs.
41 * \param res The result is stored here. If the \p use_seed argument
42 * is true (non-zero), then this matrix is also used as the
43 * starting point of the algorithm.
44 * \param use_seed Boolean, whether to use the supplied coordinates in
45 * \p res as the starting point. If false (zero), then a
46 * uniform random starting point is used.
47 * \param maxiter The maximum number of iterations to
48 * perform. Updating a single vertex counts as an iteration.
49 * A reasonable default is 40 * n * n, where n is the number of
50 * vertices. The original paper suggests 4 * n * n, but this
51 * usually only works if the other parameters are set up carefully.
52 * \param temp_max The maximum allowed local temperature. A reasonable
53 * default is the number of vertices.
54 * \param temp_min The global temperature at which the algorithm
55 * terminates (even before reaching \p maxiter iterations). A
56 * reasonable default is 1/10.
57 * \param temp_init Initial local temperature of all vertices. A
58 * reasonable default is the square root of the number of
59 * vertices.
60 * \return Error code.
61 *
62 * Time complexity: O(t * n * (n+e)), where n is the number of vertices,
63 * e is the number of edges and t is the number of time steps
64 * performed.
65 */
66
igraph_layout_gem(const igraph_t * graph,igraph_matrix_t * res,igraph_bool_t use_seed,igraph_integer_t maxiter,igraph_real_t temp_max,igraph_real_t temp_min,igraph_real_t temp_init)67 int igraph_layout_gem(const igraph_t *graph, igraph_matrix_t *res,
68 igraph_bool_t use_seed, igraph_integer_t maxiter,
69 igraph_real_t temp_max, igraph_real_t temp_min,
70 igraph_real_t temp_init) {
71
72 igraph_integer_t no_nodes = igraph_vcount(graph);
73 igraph_vector_int_t perm;
74 igraph_vector_float_t impulse_x, impulse_y, temp, skew_gauge;
75 igraph_integer_t i;
76 float temp_global;
77 igraph_integer_t perm_pointer = 0;
78 float barycenter_x = 0, barycenter_y = 0;
79 igraph_vector_t phi;
80 igraph_vector_t neis;
81 const float elen_des2 = 128 * 128;
82 const float gamma = 1 / 16.0f;
83 const float alpha_o = (float)M_PI;
84 const float alpha_r = (float)M_PI / 3.0f;
85 const float sigma_o = 1.0f / 3.0f;
86 const float sigma_r = 1.0f / 2.0f / no_nodes;
87
88 if (maxiter < 0) {
89 IGRAPH_ERROR("Number of iterations must be non-negative in GEM layout",
90 IGRAPH_EINVAL);
91 }
92 if (use_seed && (igraph_matrix_nrow(res) != no_nodes ||
93 igraph_matrix_ncol(res) != 2)) {
94 IGRAPH_ERROR("Invalid start position matrix size in GEM layout",
95 IGRAPH_EINVAL);
96 }
97 if (temp_max <= 0) {
98 IGRAPH_ERROR("Maximum temperature should be positive in GEM layout",
99 IGRAPH_EINVAL);
100 }
101 if (temp_min <= 0) {
102 IGRAPH_ERROR("Minimum temperature should be positive in GEM layout",
103 IGRAPH_EINVAL);
104 }
105 if (temp_init <= 0) {
106 IGRAPH_ERROR("Initial temperature should be positive in GEM layout",
107 IGRAPH_EINVAL);
108 }
109 if (temp_max < temp_init || temp_init < temp_min) {
110 IGRAPH_ERROR("Minimum <= Initial <= Maximum temperature is required "
111 "in GEM layout", IGRAPH_EINVAL);
112 }
113
114 if (no_nodes == 0) {
115 return 0;
116 }
117
118 IGRAPH_CHECK(igraph_vector_float_init(&impulse_x, no_nodes));
119 IGRAPH_FINALLY(igraph_vector_float_destroy, &impulse_x);
120 IGRAPH_CHECK(igraph_vector_float_init(&impulse_y, no_nodes));
121 IGRAPH_FINALLY(igraph_vector_float_destroy, &impulse_y);
122 IGRAPH_CHECK(igraph_vector_float_init(&temp, no_nodes));
123 IGRAPH_FINALLY(igraph_vector_float_destroy, &temp);
124 IGRAPH_CHECK(igraph_vector_float_init(&skew_gauge, no_nodes));
125 IGRAPH_FINALLY(igraph_vector_float_destroy, &skew_gauge);
126 IGRAPH_CHECK(igraph_vector_int_init_seq(&perm, 0, no_nodes - 1));
127 IGRAPH_FINALLY(igraph_vector_int_destroy, &perm);
128 IGRAPH_VECTOR_INIT_FINALLY(&phi, no_nodes);
129 IGRAPH_VECTOR_INIT_FINALLY(&neis, 10);
130
131 RNG_BEGIN();
132
133 /* Initialization */
134 igraph_degree(graph, &phi, igraph_vss_all(), IGRAPH_ALL, IGRAPH_LOOPS);
135 if (!use_seed) {
136 const igraph_real_t width_half = no_nodes * 100, height_half = width_half;
137 IGRAPH_CHECK(igraph_matrix_resize(res, no_nodes, 2));
138 for (i = 0; i < no_nodes; i++) {
139 MATRIX(*res, i, 0) = RNG_UNIF(-width_half, width_half);
140 MATRIX(*res, i, 1) = RNG_UNIF(-height_half, height_half);
141 barycenter_x += MATRIX(*res, i, 0);
142 barycenter_y += MATRIX(*res, i, 1);
143 VECTOR(phi)[i] *= (VECTOR(phi)[i] / 2.0 + 1.0);
144 }
145 } else {
146 for (i = 0; i < no_nodes; i++) {
147 barycenter_x += MATRIX(*res, i, 0);
148 barycenter_y += MATRIX(*res, i, 1);
149 VECTOR(phi)[i] *= (VECTOR(phi)[i] / 2.0 + 1.0);
150 }
151 }
152 igraph_vector_float_fill(&temp, temp_init);
153 temp_global = temp_init * no_nodes;
154
155 while (temp_global > temp_min * no_nodes && maxiter > 0) {
156 igraph_integer_t u, v, nlen, j;
157 float px, py, pvx, pvy;
158
159 IGRAPH_ALLOW_INTERRUPTION();
160
161 /* choose a vertex v to update */
162 if (!perm_pointer) {
163 igraph_vector_int_shuffle(&perm);
164 perm_pointer = no_nodes - 1;
165 }
166 v = VECTOR(perm)[perm_pointer--];
167
168 /* compute v's impulse */
169 px = (barycenter_x / no_nodes - MATRIX(*res, v, 0)) * gamma * VECTOR(phi)[v];
170 py = (barycenter_y / no_nodes - MATRIX(*res, v, 1)) * gamma * VECTOR(phi)[v];
171 px += RNG_UNIF(-32.0, 32.0);
172 py += RNG_UNIF(-32.0, 32.0);
173
174 for (u = 0; u < no_nodes; u++) {
175 float dx, dy, dist2;
176 if (u == v) {
177 continue;
178 }
179 dx = MATRIX(*res, v, 0) - MATRIX(*res, u, 0);
180 dy = MATRIX(*res, v, 1) - MATRIX(*res, u, 1);
181 dist2 = dx * dx + dy * dy;
182 if (dist2 != 0) {
183 px += dx * elen_des2 / dist2;
184 py += dy * elen_des2 / dist2;
185 }
186 }
187
188 IGRAPH_CHECK(igraph_neighbors(graph, &neis, v, IGRAPH_ALL));
189 nlen = igraph_vector_size(&neis);
190 for (j = 0; j < nlen; j++) {
191 igraph_integer_t u = VECTOR(neis)[j];
192 float dx = MATRIX(*res, v, 0) - MATRIX(*res, u, 0);
193 float dy = MATRIX(*res, v, 1) - MATRIX(*res, u, 1);
194 float dist2 = dx * dx + dy * dy;
195 px -= dx * dist2 / (elen_des2 * VECTOR(phi)[v]);
196 py -= dy * dist2 / (elen_des2 * VECTOR(phi)[v]);
197 }
198
199 /* update v's position and temperature */
200 if (px != 0 || py != 0) {
201 float plen = sqrtf(px * px + py * py);
202 px *= VECTOR(temp)[v] / plen;
203 py *= VECTOR(temp)[v] / plen;
204 MATRIX(*res, v, 0) += px;
205 MATRIX(*res, v, 1) += py;
206 barycenter_x += px;
207 barycenter_y += py;
208 }
209
210 pvx = VECTOR(impulse_x)[v]; pvy = VECTOR(impulse_y)[v];
211 if (pvx != 0 || pvy != 0) {
212 float beta = atan2f(pvy - py, pvx - px);
213 float sin_beta = sinf(beta);
214 float sign_sin_beta = (sin_beta > 0) ? 1 : ((sin_beta < 0) ? -1 : 0);
215 float cos_beta = cosf(beta);
216 float abs_cos_beta = fabsf(cos_beta);
217 float old_temp = VECTOR(temp)[v];
218 if (sin(beta) >= sin(M_PI_2 + alpha_r / 2.0)) {
219 VECTOR(skew_gauge)[v] += sigma_r * sign_sin_beta;
220 }
221 if (abs_cos_beta >= cosf(alpha_o / 2.0)) {
222 VECTOR(temp)[v] *= sigma_o * cos_beta;
223 }
224 VECTOR(temp)[v] *= (1 - fabsf(VECTOR(skew_gauge)[v]));
225 if (VECTOR(temp)[v] > temp_max) {
226 VECTOR(temp)[v] = temp_max;
227 }
228 VECTOR(impulse_x)[v] = px;
229 VECTOR(impulse_y)[v] = py;
230 temp_global += VECTOR(temp)[v] - old_temp;
231 }
232
233 maxiter--;
234
235 } /* while temp && iter */
236
237
238 RNG_END();
239
240 igraph_vector_destroy(&neis);
241 igraph_vector_destroy(&phi);
242 igraph_vector_int_destroy(&perm);
243 igraph_vector_float_destroy(&skew_gauge);
244 igraph_vector_float_destroy(&temp);
245 igraph_vector_float_destroy(&impulse_y);
246 igraph_vector_float_destroy(&impulse_x);
247 IGRAPH_FINALLY_CLEAN(7);
248
249 return 0;
250 }
251