graph_sfdp.hh 14.1 KB
Newer Older
Tiago Peixoto's avatar
Tiago Peixoto committed
1
2
// graph-tool -- a general graph modification and manipulation thingy
//
Tiago Peixoto's avatar
Tiago Peixoto committed
3
// Copyright (C) 2006-2016 Tiago de Paula Peixoto <tiago@skewed.de>
Tiago Peixoto's avatar
Tiago Peixoto committed
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
//
// This program is free software; you can redistribute it and/or
// modify it under the terms of the GNU General Public License
// as published by the Free Software Foundation; either version 3
// of the License, or (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <http://www.gnu.org/licenses/>.

#ifndef GRAPH_FDP_HH
#define GRAPH_FDP_HH

#include <limits>
#include <iostream>
Tiago Peixoto's avatar
Tiago Peixoto committed
23
24

#ifndef __clang__
25
26
#include <ext/numeric>
using __gnu_cxx::power;
Tiago Peixoto's avatar
Tiago Peixoto committed
27
28
29
30
31
32
33
#else
template <class Value>
Value power(Value value, int n)
{
    return pow(value, n);
}
#endif
Tiago Peixoto's avatar
Tiago Peixoto committed
34
35
36
37
38
39
40
41
42
43
44
45
46
47

namespace graph_tool
{
using namespace std;
using namespace boost;

template <class Pos, class Weight>
class QuadTree
{
public:
    QuadTree(const Pos& ll, const Pos& ur, int max_level)
        :_ll(ll), _ur(ur), _cm(2, 0), _count(0),
         _max_level(max_level)
    {
48
49
        _w = sqrt(power(_ur[0] - _ll[0], 2) +
                  power(_ur[1] - _ll[1], 2));
Tiago Peixoto's avatar
Tiago Peixoto committed
50
51
    }

52
    vector<QuadTree>& get_leafs()
Tiago Peixoto's avatar
Tiago Peixoto committed
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
    {
        if (_max_level > 0 && _leafs.empty())
        {
            _leafs.reserve(4);
            for (size_t i = 0; i < 4; ++i)
            {
                Pos lll = _ll, lur = _ur;
                if (i % 2)
                    lll[0] += (_ur[0] - _ll[0]) / 2;
                else
                    lur[0] -= (_ur[0] - _ll[0]) / 2;
                if (i / 2)
                    lll[1] += (_ur[1] - _ll[1]) / 2;
                else
                    lur[1] -= (_ur[1] - _ll[1]) / 2;
68
                _leafs.emplace_back(lll, lur, _max_level - 1);
Tiago Peixoto's avatar
Tiago Peixoto committed
69
70
71
            }
        }

72
        return _leafs;
Tiago Peixoto's avatar
Tiago Peixoto committed
73
74
    }

Tiago Peixoto's avatar
Tiago Peixoto committed
75
    vector<std::tuple<Pos,Weight> >& get_dense_leafs()
Tiago Peixoto's avatar
Tiago Peixoto committed
76
77
78
79
80
81
82
83
84
85
86
87
    {
        return _dense_leafs;
    }

    size_t put_pos(Pos& p, Weight w)
    {
        _count += w;
        _cm[0] += p[0] * w;
        _cm[1] += p[1] * w;

        if (_max_level == 0)
        {
Tiago Peixoto's avatar
Tiago Peixoto committed
88
            _dense_leafs.push_back(std::make_tuple(p, w));
Tiago Peixoto's avatar
Tiago Peixoto committed
89
90
91
92
93
94
95
            return 1;
        }
        else
        {
            int i = p[0] > (_ll[0] + (_ur[0] - _ll[0]) / 2);
            int j = p[1] > (_ll[1] + (_ur[1] - _ll[1]) / 2);
            size_t sc = (_max_level > 0 && _leafs.empty()) ? 4 : 0;
96
            return sc + get_leafs()[i + 2 * j].put_pos(p, w);
Tiago Peixoto's avatar
Tiago Peixoto committed
97
98
99
100
101
102
103
104
105
106
        }
        return 0;
    }

    void get_cm(Pos& cm)
    {
        for (size_t i = 0; i < 2; ++i)
            cm[i] = _cm[i] / _count;
    }

107
    double get_w() const
Tiago Peixoto's avatar
Tiago Peixoto committed
108
    {
109
        return _w;
Tiago Peixoto's avatar
Tiago Peixoto committed
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
    }

    Weight get_count()
    {
        return _count;
    }

    int max_level()
    {
        return _max_level;
    }

private:
    Pos _ll, _ur;
    vector<QuadTree> _leafs;
Tiago Peixoto's avatar
Tiago Peixoto committed
125
    vector<std::tuple<Pos,Weight> > _dense_leafs;
Tiago Peixoto's avatar
Tiago Peixoto committed
126
127
128
    Pos _cm;
    Weight _count;
    int _max_level;
129
    double _w;
Tiago Peixoto's avatar
Tiago Peixoto committed
130
131
132
133
134
135
136
};

template <class Pos>
inline double dist(const Pos& p1, const Pos& p2)
{
    double r = 0;
    for (size_t i = 0; i < 2; ++i)
Tiago Peixoto's avatar
Tiago Peixoto committed
137
        r += power(double(p1[i] - p2[i]), 2);
Tiago Peixoto's avatar
Tiago Peixoto committed
138
139
140
141
142
143
144
145
146
    return sqrt(r);
}

template <class Pos>
inline double f_r(double C, double K, double p, const Pos& p1, const Pos& p2)
{
    double d = dist(p1, p2);
    if (d == 0)
        return 0;
147
148
149
150
    if (round(p) == p)
        return -C * power(K, int(1 + p)) / power(d, int(p));
    else
        return -C * pow(K, 1 + p) / pow(d, p);
Tiago Peixoto's avatar
Tiago Peixoto committed
151
152
153
154
155
}

template <class Pos>
inline double f_a(double K, const Pos& p1, const Pos& p2)
{
Tiago Peixoto's avatar
Tiago Peixoto committed
156
    return power(dist(p1, p2), 2) / K;
Tiago Peixoto's avatar
Tiago Peixoto committed
157
158
159
160
161
162
163
164
165
}

template <class Pos>
inline double get_diff(const Pos& p1, const Pos& p2, Pos& r)
{
    double abs = 0;
    for (size_t i = 0; i < 2; ++i)
    {
        r[i] = p1[i] - p2[i];
166
        abs += r[i] * r[i];
Tiago Peixoto's avatar
Tiago Peixoto committed
167
168
169
    }
    if (abs == 0)
        abs = 1;
170
    abs = sqrt(abs);
Tiago Peixoto's avatar
Tiago Peixoto committed
171
    for (size_t i = 0; i < 2; ++i)
172
173
        r[i] /= abs;
    return abs;
Tiago Peixoto's avatar
Tiago Peixoto committed
174
175
176
177
178
179
180
}

template <class Pos>
inline double norm(Pos& x)
{
    double abs = 0;
    for (size_t i = 0; i < 2; ++i)
Tiago Peixoto's avatar
Tiago Peixoto committed
181
        abs += power(x[i], 2);
Tiago Peixoto's avatar
Tiago Peixoto committed
182
183
184
185
186
187
188
    for (size_t i = 0; i < 2; ++i)
        x[i] /= sqrt(abs);
    return sqrt(abs);
}

struct get_sfdp_layout
{
189
    get_sfdp_layout(double C, double K, double p, double theta, double gamma,
190
191
192
193
194
195
196
197
198
                    double mu, double mu_p, double init_step,
                    double step_schedule, size_t max_level, double epsilon,
                    size_t max_iter, bool simple)
        : C(C), K(K), p(p), theta(theta), gamma(gamma), mu(mu), mu_p(mu_p),
          init_step(init_step), step_schedule(step_schedule),
          epsilon(epsilon), max_level(max_level), max_iter(max_iter),
          simple(simple) {}

    double C, K, p, theta, gamma, mu, mu_p, init_step, step_schedule, epsilon;
Tiago Peixoto's avatar
Tiago Peixoto committed
199
200
201
    size_t max_level, max_iter;
    bool simple;

202
    template <class Graph, class PosMap, class VertexWeightMap,
203
              class EdgeWeightMap, class PinMap, class GroupMap, class RNG>
204
205
206
    void operator()(Graph& g, PosMap pos, VertexWeightMap vweight,
                    EdgeWeightMap eweight, PinMap pin, GroupMap group,
                    bool verbose, RNG& rng) const
Tiago Peixoto's avatar
Tiago Peixoto committed
207
208
209
210
211
212
    {
        typedef typename property_traits<PosMap>::value_type pos_t;
        typedef typename property_traits<PosMap>::value_type::value_type val_t;

        typedef typename property_traits<VertexWeightMap>::value_type vweight_t;

213
        vector<pos_t> group_cm;
214
        vector<vweight_t> group_size;
215
        vector<size_t> vertices;
216

217
218
        int  HN=0;
        for (auto v : vertices_range(g))
Tiago Peixoto's avatar
Tiago Peixoto committed
219
        {
220
            if (pin[v] == 0)
221
                vertices.push_back(v);
222
            pos[v].resize(2, 0);
223
            if (gamma != 0 || mu != 0)
224
            {
225
                size_t s = group[v];
Tiago Peixoto's avatar
Tiago Peixoto committed
226

227
228
229
230
231
232
233
234
235
236
                if (s >= group_cm.size())
                {
                    group_cm.resize(s + 1);
                    group_size.resize(s + 1, 0);
                }
                group_cm[s].resize(2, 0);
                group_size[s] += get(vweight, v);

                for (size_t j = 0; j < 2; ++j)
                    group_cm[s][j] += pos[v][j] * get(vweight, v);
Tiago Peixoto's avatar
Tiago Peixoto committed
237
            }
238
239
240
241
242
            HN++;
        }

        for (size_t s = 0; s < group_size.size(); ++s)
        {
243
244
            if (group_size[s] == 0)
                continue;
245
246
247
            group_cm[s].resize(2, 0);
            for (size_t j = 0; j < 2; ++j)
                group_cm[s][j] /= group_size[s];
Tiago Peixoto's avatar
Tiago Peixoto committed
248
249
250
251
252
253
254
255
256
257
258
259
260
261
        }

        val_t delta = epsilon * K + 1, E = 0, E0;
        E0 = numeric_limits<val_t>::max();
        size_t n_iter = 0;
        val_t step = init_step;
        size_t progress = 0;

        while (delta > epsilon * K && (max_iter == 0 || n_iter < max_iter))
        {
            delta = 0;
            E0 = E;
            E = 0;

262
263
264
265
266
267
268
269
270
271
            pos_t ll(2, numeric_limits<val_t>::max()),
                ur(2, -numeric_limits<val_t>::max());
            for (auto v : vertices_range(g))
            {
                for (size_t j = 0; j < 2; ++j)
                {
                    ll[j] = min(pos[v][j], ll[j]);
                    ur[j] = max(pos[v][j], ur[j]);
                }
            }
Tiago Peixoto's avatar
Tiago Peixoto committed
272

273
            if (gamma != 0 || mu != 0)
Tiago Peixoto's avatar
Tiago Peixoto committed
274
            {
275
276
277
278
279
280
281
282
283
284
285
286
287
288
                for (size_t s = 0; s < group_size.size(); ++s)
                {
                    if (group_size[s] == 0)
                        continue;
                    group_cm[s] = {0, 0};
                }

                for (auto v : vertices_range(g))
                {
                    size_t s = group[v];
                    for (size_t j = 0; j < 2; ++j)
                        group_cm[s][j] += pos[v][j] * get(vweight, v) /
                            group_size[s];
                }
Tiago Peixoto's avatar
Tiago Peixoto committed
289
290
            }

291
292
293
294
            QuadTree<pos_t, vweight_t> qt(ll, ur, max_level);
            for (auto v : vertices_range(g))
                qt.put_pos(pos[v], vweight[v]);

295
296
            std::shuffle(vertices.begin(), vertices.end(), rng);

Tiago Peixoto's avatar
Tiago Peixoto committed
297
            size_t nmoves = 0;
298
299
            int i = 0, N = vertices.size();
            vector<QuadTree<pos_t, vweight_t>*> Q;
300
            #pragma omp parallel for default(shared) private(i, Q) \
301
                reduction(+:E, delta, nmoves) schedule(runtime) if (N > 100)
Tiago Peixoto's avatar
Tiago Peixoto committed
302
303
            for (i = 0; i < N; ++i)
            {
304
                auto v = vertex(vertices[i], g);
Tiago Peixoto's avatar
Tiago Peixoto committed
305

306
                pos_t diff(2, 0), pos_u(2, 0), ftot(2, 0), cm(2, 0);
307
308

                // global repulsive forces
309
                Q.push_back(&qt);
Tiago Peixoto's avatar
Tiago Peixoto committed
310
311
                while (!Q.empty())
                {
312
                    auto& q = *Q.back();
Tiago Peixoto's avatar
Tiago Peixoto committed
313
314
315
316
                    Q.pop_back();

                    if (q.max_level() == 0)
                    {
317
                        auto& dleafs = q.get_dense_leafs();
318
                        for (auto& dleaf : dleafs)
Tiago Peixoto's avatar
Tiago Peixoto committed
319
                        {
320
                            val_t d = get_diff(get<0>(dleaf), pos[v], diff);
Tiago Peixoto's avatar
Tiago Peixoto committed
321
322
                            if (d == 0)
                                continue;
323
324
                            val_t f = f_r(C, K, p, pos[v], get<0>(dleaf));
                            f *= get<1>(dleaf) * get(vweight, v);
Tiago Peixoto's avatar
Tiago Peixoto committed
325
326
327
328
329
330
331
332
333
                            for (size_t l = 0; l < 2; ++l)
                                ftot[l] += f * diff[l];
                        }
                    }
                    else
                    {
                        double w = q.get_w();
                        q.get_cm(cm);
                        double d = get_diff(cm, pos[v], diff);
334
                        if (w > theta * d)
Tiago Peixoto's avatar
Tiago Peixoto committed
335
                        {
336
                            for (auto& leaf : q.get_leafs())
Tiago Peixoto's avatar
Tiago Peixoto committed
337
338
                            {
                                if (leaf.get_count() > 0)
339
                                    Q.push_back(&leaf);
Tiago Peixoto's avatar
Tiago Peixoto committed
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
                            }
                        }
                        else
                        {
                            if (d > 0)
                            {
                                val_t f = f_r(C, K, p, cm, pos[v]);
                                f *= q.get_count() * get(vweight, v);
                                for (size_t l = 0; l < 2; ++l)
                                    ftot[l] += f * diff[l];
                            }
                        }
                    }
                }

355
                // local attractive forces
356
                auto& pos_v = pos[v];
357
                for (auto e : out_edges_range(v, g))
Tiago Peixoto's avatar
Tiago Peixoto committed
358
                {
359
                    auto u = target(e, g);
Tiago Peixoto's avatar
Tiago Peixoto committed
360
361
                    if (u == v)
                        continue;
362
363
364
                    pos_u = pos[u];
                    get_diff(pos_u, pos_v, diff);
                    val_t f = f_a(K, pos_u, pos_v);
365
                    f *= get(eweight, e) * get(vweight, u) * get(vweight, v);
Tiago Peixoto's avatar
Tiago Peixoto committed
366
367
368
369
                    for (size_t l = 0; l < 2; ++l)
                        ftot[l] += f * diff[l];
                }

370
                // inter-group attractive forces
371
                if (gamma > 0)
372
                {
373
374
375
376
377
378
379
380
381
                    for (size_t s = 0; s < group_cm.size(); ++s)
                    {
                        if (group_size[s] == 0)
                            continue;
                        if (s == size_t(group[v]))
                            continue;
                        val_t d = get_diff(group_cm[s], pos[v], diff);
                        if (d == 0)
                            continue;
Tiago Peixoto's avatar
Tiago Peixoto committed
382
                        double Kp = K * power(HN, 2);
383
384
385
386
387
                        val_t f = f_a(Kp, group_cm[s], pos[v]) * gamma * \
                            group_size[s] * get(vweight, v);
                        for (size_t l = 0; l < 2; ++l)
                            ftot[l] += f * diff[l];
                    }
388
389
                }

390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
                // inter-group repulsive forces
                if (gamma < 0)
                {
                    for (size_t s = 0; s < group_cm.size(); ++s)
                    {
                        if (group_size[s] == 0)
                            continue;
                        if (s == size_t(group[v]))
                            continue;
                        val_t d = get_diff(group_cm[s], pos[v], diff);
                        if (d == 0)
                            continue;
                        val_t f = f_r(C, K, p, cm, pos[v]);
                        f *= group_size[s] * get(vweight, v) * abs(gamma);
                        for (size_t l = 0; l < 2; ++l)
                            ftot[l] += f * diff[l];
                    }
                }

                // intra-group attractive forces
410
                if (mu > 0 && group_size[group[v]] > 1)
411
412
413
414
                {
                    val_t d = get_diff(group_cm[group[v]], pos[v], diff);
                    if (d > 0)
                    {
415
                        double Kp = K * pow(double(group_size[group[v]]), mu_p);
416
417
418
419
420
421
                        val_t f = f_a(Kp, group_cm[group[v]], pos[v]) * mu * \
                            group_size[group[v]] * get(vweight, v);
                        for (size_t l = 0; l < 2; ++l)
                            ftot[l] += f * diff[l];
                    }
                }
422

Tiago Peixoto's avatar
Tiago Peixoto committed
423
                E += power(norm(ftot), 2);
Tiago Peixoto's avatar
Tiago Peixoto committed
424

425
                for (size_t l = 0; l < 2; ++l)
Tiago Peixoto's avatar
Tiago Peixoto committed
426
                {
427
428
                    ftot[l] *= step;
                    pos[v][l] += ftot[l];
Tiago Peixoto's avatar
Tiago Peixoto committed
429
                }
430

Tiago Peixoto's avatar
Tiago Peixoto committed
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
                delta += norm(ftot);
                nmoves++;
            }
            n_iter++;
            delta /= nmoves;

            if (verbose)
                cout << n_iter << " " << E << " " << step << " "
                     << delta << " " << max_level << endl;

            if (simple)
            {
                step *= step_schedule;
            }
            else
            {
                if (E < E0)
                {
                    ++progress;
                    if (progress >= 5)
                    {
                        progress = 0;
                        step /= step_schedule;
                    }
                }
                else
                {
                    progress = 0;
                    step *= step_schedule;
                }
            }
        }
    }
};

} // namespace graph_tool


#endif // GRAPH_FDP_HH