graph_rewiring.hh 22.7 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
// graph-tool -- a general graph modification and manipulation thingy
//
// Copyright (C) 2007  Tiago de Paula Peixoto <tiago@forked.de>
//
// 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_REWIRING_HH
#define GRAPH_REWIRING_HH

#include <tr1/unordered_set>
22
#include <tr1/random>
23
24
25
26
#include <boost/functional/hash.hpp>

#include "graph.hh"
#include "graph_filtering.hh"
27
#include "graph_util.hh"
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49

namespace graph_tool
{
using namespace std;
using namespace boost;


// returns true if vertices u and v are adjacent. This is O(k(u)).
template <class Graph>
bool is_adjacent(typename graph_traits<Graph>::vertex_descriptor u,
                 typename graph_traits<Graph>::vertex_descriptor v,
                 const Graph& g )
{
    typename graph_traits<Graph>::out_edge_iterator e, e_end;
    for (tie(e, e_end) = out_edges(u, g); e != e_end; ++e)
    {
        if (target(*e,g) == v)
            return true;
    }
    return false;
}

50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
template <class Graph>
typename graph_traits<Graph>::vertex_descriptor
source(const pair<typename graph_traits<Graph>::edge_descriptor, bool>& e,
       const Graph& g)
{
    if (e.second)
        return target(e.first, g);
    else
        return source(e.first, g);
}

template <class Graph>
typename graph_traits<Graph>::vertex_descriptor
target(const pair<typename graph_traits<Graph>::edge_descriptor, bool>& e,
       const Graph& g)
{
    if (e.second)
        return source(e.first, g);
    else
        return target(e.first, g);
}

72
73
74
75
76
// this functor will swap the source of the edge e with the source of edge se
// and the target of edge e with the target of te
struct swap_edge_triad
{
    template <class Graph, class NewEdgeMap>
77
78
79
80
81
82
83
    static bool
    parallel_check(const typename graph_traits<Graph>::edge_descriptor& e,
                   const pair<typename graph_traits<Graph>::edge_descriptor,
                              bool>& se,
                   const pair<typename graph_traits<Graph>::edge_descriptor,
                              bool>& te,
                   NewEdgeMap edge_is_new, const Graph &g)
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
    {
        // We want to check that if we swap the source of 'e' with the source of
        // 'se', and the target of 'te' with the target of 'e', as such
        //
        //  (s)    -e--> (t)          (ns)   -e--> (nt)
        //  (ns)   -se-> (se_t)   =>  (s)    -se-> (se_t)
        //  (te_s) -te-> (nt)         (te_s) -te-> (t),
        //
        // no parallel edges are introduced. We must considered only "new
        // edges", i.e., edges which were already sampled and swapped. "Old
        // edges" will have their chance of being swapped, and then they'll be
        // checked for parallelism.

        typename graph_traits<Graph>::vertex_descriptor
            s = source(e, g),          // current source
            t = target(e, g),          // current target
            ns = source(se, g),        // new source
101
102
            nt = target(te, g),        // new target
            te_s = source(te, g),      // target edge source
103
104
105
            se_t = target(se, g);      // source edge target


106
        if (edge_is_new[se.first] && (ns == s) && (nt == se_t))
107
            return true; // e is parallel to se after swap
108
        if (edge_is_new[te.first] && (te_s == ns) && (nt == t))
109
            return true; // e is parallel to te after swap
110
        if (edge_is_new[te.first] && edge_is_new[se.first] && (te != se) &&
111
112
             (s == te_s) && (t == se_t))
            return true; // se is parallel to te after swap
113
        if (is_adjacent_in_new(ns,  nt, edge_is_new, g))
114
            return true; // e would clash with an existing (new) edge
115
        if (edge_is_new[te.first] && is_adjacent_in_new(te_s, t, edge_is_new,g))
116
            return true; // te would clash with an existing (new) edge
117
        if (edge_is_new[se.first] && is_adjacent_in_new(s, se_t, edge_is_new,g))
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
            return true; // se would clash with an existing (new) edge
        return false; // the coast is clear - hooray!
    }

    // returns true if vertices u and v are adjacent in the new graph. This is
    // O(k(u)).
    template <class Graph, class EdgeIsNew>
    static bool is_adjacent_in_new
        (typename graph_traits<Graph>::vertex_descriptor u,
         typename graph_traits<Graph>::vertex_descriptor v,
         EdgeIsNew edge_is_new, const Graph& g)
    {
        typename graph_traits<Graph>::out_edge_iterator e, e_end;
        for (tie(e, e_end) = out_edges(u, g); e != e_end; ++e)
        {
133
            if (edge_is_new[*e] && target(*e,g) == v)
134
135
136
137
138
139
                return true;
        }
        return false;
    }

    template <class Graph, class EdgeIndexMap, class EdgesType>
140
141
142
143
144
    void operator()(const typename graph_traits<Graph>::edge_descriptor& e,
                    const pair<typename graph_traits<Graph>::edge_descriptor,
                               bool>& se,
                    const pair<typename graph_traits<Graph>::edge_descriptor,
                               bool>& te,
145
146
147
148
149
150
151
152
153
154
155
156
157
158
                    EdgesType& edges, EdgeIndexMap edge_index, Graph& g)
    {
        // swap the source of the edge 'e' with the source of edge 'se' and the
        // target of edge 'e' with the target of 'te', as such:
        //
        //  (s)    -e--> (t)          (ns)   -e--> (nt)
        //  (ns)   -se-> (se_t)   =>  (s)    -se-> (se_t)
        //  (te_s) -te-> (nt)         (te_s) -te-> (t),

        // new edges which will replace the old ones
        typename graph_traits<Graph>::edge_descriptor ne, nse, nte;

        // split cases where different combinations of the three edges are
        // the same
159
        if(se.first != te.first)
160
        {
161
            ne = add_edge(source(se, g), target(te, g), g).first;
162
            if(e != se.first)
163
            {
164
165
166
167
                if (!se.second)
                    nse = add_edge(source(e, g), target(se, g), g).first;
                else // keep invertedness (only for undirected graphs)
                    nse = add_edge(target(se, g), source(e, g), g).first;
168
169
                edge_index[nse] = edge_index[se.first];
                remove_edge(se.first, g);
170
171
                edges[edge_index[nse]] = nse;
            }
172
            if(e != te.first)
173
            {
174
175
176
177
                if (!te.second)
                    nte = add_edge(source(te, g), target(e, g), g).first;
                else // keep invertedness (only for undirected graphs)
                    nte = add_edge(target(e, g), source(te, g), g).first;
178
179
                edge_index[nte] = edge_index[te.first];
                remove_edge(te.first, g);
180
181
182
183
184
185
186
187
                edges[edge_index[nte]] = nte;
            }
            edge_index[ne] = edge_index[e];
            remove_edge(e, g);
            edges[edge_index[ne]] = ne;
        }
        else
        {
188
            if(e != se.first)
189
190
            {
                // se and te are the same. swapping indexes only.
191
192
                swap(edge_index[se.first], edge_index[e]);
                edges[edge_index[se.first]] = se.first;
193
194
195
196
197
198
199
200
201
202
203
                edges[edge_index[e]] = e;
            }
        }
    }
};

// main rewire loop
template <template <class Graph, class EdgeIndexMap> class RewireStrategy>
struct graph_rewire
{
    template <class Graph, class EdgeIndexMap>
204
    void operator()(Graph& g, EdgeIndexMap edge_index, rng_t& rng,
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
                    bool self_loops, bool parallel_edges) const
    {
        typedef typename graph_traits<Graph>::vertex_descriptor vertex_t;
        typedef typename graph_traits<Graph>::edge_descriptor edge_t;

        if (!self_loops)
        {
            // check the existence of self-loops
            bool has_self_loops = false;
            int i, N = num_vertices(g);
            #pragma omp parallel for default(shared) private(i) \
                schedule(dynamic)
            for (i = 0; i < N; ++i)
            {
                vertex_t v = vertex(i, g);
                if (v == graph_traits<Graph>::null_vertex())
                    continue;
                if (is_adjacent(v, v, g))
                    has_self_loops = true;
            }
            if (has_self_loops)
Tiago Peixoto's avatar
Tiago Peixoto committed
226
                throw ValueException("Self-loop detected. Can't rewire graph "
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
                                     "without self-loops if it already contains"
                                     " self-loops!");
        }

        if (!parallel_edges)
        {
            // check the existence of parallel edges
            bool has_parallel_edges = false;
            int i, N = num_vertices(g);
            #pragma omp parallel for default(shared) private(i) \
                schedule(dynamic)
            for (i = 0; i < N; ++i)
            {
                vertex_t v = vertex(i, g);
                if (v == graph_traits<Graph>::null_vertex())
                    continue;

                tr1::unordered_set<vertex_t> targets;
                typename graph_traits<Graph>::out_edge_iterator e, e_end;
                for (tie(e, e_end) = out_edges(v, g); e != e_end; ++e)
                {
                    if (targets.find(target(*e, g)) != targets.end())
                        has_parallel_edges = true;
                    else
                        targets.insert(target(*e, g));
                }
            }

            if (has_parallel_edges)
Tiago Peixoto's avatar
Tiago Peixoto committed
256
                throw ValueException("Parallel edge detected. Can't rewire "
257
258
259
260
261
262
263
                                     "graph without parallel edges if it "
                                     "already contains parallel edges!");
        }

        RewireStrategy<Graph, EdgeIndexMap> rewire(g, edge_index, rng);

        vector<edge_t> edges(num_edges(g));
264
265
266
        vector<bool> is_edge(num_edges(g), false);
        typename graph_traits<Graph>::edge_iterator e, e_end;
        for (tie(e, e_end) = boost::edges(g); e != e_end; ++e)
267
        {
268
269
270
271
272
273
274
            if (edge_index[*e] >= edges.size())
            {
                edges.resize(edge_index[*e] + 1);
                is_edge.resize(edge_index[*e] + 1, false);
            }
            edges[edge_index[*e]] = *e;
            is_edge[edge_index[*e]] = true;
275
276
277
        }

        // for each edge simultaneously rewire its source and target
278
        for (size_t i = 0; i < edges.size(); ++i)
279
        {
280
281
            if (!is_edge[i])
                continue;
282
283
            typename graph_traits<Graph>::edge_descriptor e = edges[i];
            typename graph_traits<Graph>::edge_descriptor se, te;
284
            rewire(e, edges, is_edge, self_loops, parallel_edges);
285
286
287
288
289
290
        }
    }
};

// This will iterate over a random permutation of a random access sequence, by
// swapping the values of the sequence as it iterates
291
292
template <class RandomAccessIterator, class RNG,
          class RandomDist = tr1::uniform_int<size_t> >
293
294
class random_permutation_iterator : public
    std::iterator<input_iterator_tag, typename RandomAccessIterator::value_type>
295
296
{
public:
297
298
299
    random_permutation_iterator(RandomAccessIterator begin,
                                RandomAccessIterator end, RNG& rng)
        : _i(begin), _end(end), _rng(&rng)
300
    {
301
302
303
304
305
        if(_i != _end)
        {
            RandomDist random(0,  _end - _i - 1);
            std::iter_swap(_i, _i + random(*_rng));
        }
306
    }
307

308
309
310
311
    typename RandomAccessIterator::value_type operator*()
    {
        return *_i;
    }
312

313
314
315
    random_permutation_iterator& operator++()
    {
        ++_i;
316
        if(_i != _end)
317
        {
318
319
            RandomDist random(0,  _end - _i - 1);
            std::iter_swap(_i, _i + random(*_rng));
320
        }
321
322
        return *this;
    }
323

324
    bool operator==(const random_permutation_iterator& ri)
325
    {
326
        return _i == ri._i;
327
    }
328

329
    bool operator!=(const random_permutation_iterator& ri)
330
    {
331
        return _i != ri._i;
332
    }
333
334
335
336
337
338

    size_t operator-(const random_permutation_iterator& ri)
    {
        return _i - ri._i;
    }

339
private:
340
341
    RandomAccessIterator _i, _end;
    RNG* _rng;
342
343
};

344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
// this will rewire the edges so that the resulting graph will be entirely
// random (i.e. Erdos-Renyi)
template <class Graph, class EdgeIndexMap>
class ErdosRewireStrategy
{
public:
    typedef typename graph_traits<Graph>::edge_descriptor edge_t;
    typedef typename EdgeIndexMap::value_type index_t;

    ErdosRewireStrategy(Graph& g, EdgeIndexMap edge_index, rng_t& rng)
        : _g(g), _edge_index(edge_index), _vertices(HardNumVertices()(g)),
          _rng(rng)
    {
        typeof(_vertices.begin()) viter = _vertices.begin();
        typename graph_traits<Graph>::vertex_iterator v, v_end;
        for (tie(v, v_end) = vertices(_g); v != v_end; ++v)
            *(viter++) = *v;
    }

    template<class EdgesType>
    void operator()(const edge_t& e, EdgesType& edges, vector<bool>& is_edge,
                    bool self_loops, bool parallel_edges)
    {
        //try randomly drawn pairs of vertices until one satisfies all the
        //consistency checks
        typedef random_permutation_iterator
            <typename graph_traits<Graph>::vertex_iterator, rng_t>
            random_vertex_iter;

        tr1::uniform_int<size_t> sample(0, _vertices.size());
        typename graph_traits<Graph>::vertex_descriptor s, t;
        while (true)
        {
            s = sample(_rng);
            t = sample(_rng);

            if(s == t && !self_loops) // reject self-loops if not allowed
                continue;
            if (!parallel_edges &&
                swap_edge_triad::is_adjacent_in_new(s, t, _edge_is_new, _g))
                continue;  // reject parallel edges if not allowed
            break;
        }
        edge_t ne = add_edge(s, t, _g).first;
        edges[_edge_index[e]] = ne;
        remove_edge(e, _g);
        if (_edge_index[ne] >= edges.size())
        {
            edges.resize(_edge_index[ne] + 1);
            is_edge.resize(_edge_index[ne] + 1, false);
        }
        edges[_edge_index[ne]] = ne;
        is_edge[_edge_index[ne]] = true;

        _edge_is_new[ne] = true;
    }

private:
    Graph& _g;
    EdgeIndexMap _edge_index;
    vector<typename graph_traits<Graph>::vertex_descriptor> _vertices;
405
    checked_vector_property_map<bool, EdgeIndexMap> _edge_is_new;
406
407
408
409
    rng_t& _rng;
};


410
411
412
// this is the mother class for edge-based rewire strategies
// it contains the common loop for finding edges to swap, so different
// strategies need only to specify where to sample the edges from.
413
414
template <class Graph, class EdgeIndexMap, class RewireStrategy>
class RewireStrategyBase
415
416
417
{
public:
    typedef typename graph_traits<Graph>::edge_descriptor edge_t;
418
419
    typedef typename graph_traits<Graph>::edge_descriptor vertex_t;

420
421
    typedef typename EdgeIndexMap::value_type index_t;

422
423
    RewireStrategyBase(Graph& g, EdgeIndexMap edge_index, rng_t& rng)
        : _g(g), _edge_index(edge_index), _edge_is_new(edge_index), _rng(rng) {}
424
425

    template<class EdgesType>
426
427
    void operator()(const edge_t& e, EdgesType& edges, vector<bool>& is_edge,
                    bool self_loops, bool parallel_edges)
428
    {
429
        // where should we sample the edges from
430
        vector<pair<index_t,bool> >* edges_source=0, *edges_target=0;
431
432
        static_cast<RewireStrategy*>(this)->get_edges(e, edges_source,
                                                      edges_target);
433
434
435
436

        //try randomly drawn pairs of edges until one satisfies all the
        //consistency checks
        bool found = false;
437
        pair<edge_t,bool> es, et;
438
        typedef random_permutation_iterator
439
440
            <typename vector<pair<index_t,bool> >::iterator, rng_t>
            random_edge_iter;
441

442
443
444
445
446
        random_edge_iter esi(edges_source->begin(), edges_source->end(),
                             _rng),
                         esi_end(edges_source->end(), edges_source->end(),
                             _rng);
        for (; esi != esi_end && !found; ++esi)
447
        {
448
            if (!is_edge[(*esi).first])
449
                continue;
450
            es = make_pair(edges[(*esi).first], (*esi).second);
451

452
453
            if(!self_loops) // reject self-loops if not allowed
            {
454
                if((source(e, _g) == target(es, _g)))
455
456
457
                    continue;
            }

458
            random_edge_iter eti(edges_target->begin(), edges_target->end(),
459
460
                                 _rng),
                             eti_end(edges_target->end(), edges_target->end(),
461
                                 _rng);
462
            for (; eti != eti_end && !found; ++eti)
463
            {
464
                if (!is_edge[(*eti).first])
465
                    continue;
466
                et = make_pair(edges[(*eti).first], (*eti).second);
467

468
469
                if (!self_loops) // reject self-loops if not allowed
                {
470
471
                    if ((source(es, _g) == target(et, _g)) ||
                        (source(et, _g) == target(e, _g)))
472
473
474
475
                        continue;
                }
                if (!parallel_edges) // reject parallel edges if not allowed
                {
476
477
                    if (swap_edge_triad::parallel_check(e, es, et, _edge_is_new,
                                                        _g))
478
479
480
481
482
483
484
485
                        continue;
                }
                found = true;
            }
        }
        if (!found)
            throw GraphException("Couldn't find random pair of edges to swap"
                                 "... This is a bug.");
486
        _edge_is_new[e] = true;
487
        swap_edge_triad()(e, es, et, edges, _edge_index, _g);
488
489
490
    }

private:
491
492
    Graph& _g;
    EdgeIndexMap _edge_index;
493
    checked_vector_property_map<bool, EdgeIndexMap> _edge_is_new;
494
    rng_t& _rng;
495
496
};

497
498
// this will rewire the edges so that the combined (in, out) degree distribution
// will be the same, but all the rest is random
499
template <class Graph, class EdgeIndexMap>
500
501
502
class RandomRewireStrategy:
    public RewireStrategyBase<Graph, EdgeIndexMap,
                              RandomRewireStrategy<Graph, EdgeIndexMap> >
503
504
{
public:
505
506
507
508
509
510
511
    typedef RewireStrategyBase<Graph, EdgeIndexMap,
                               RandomRewireStrategy<Graph, EdgeIndexMap> >
        base_t;

    typedef Graph graph_t;
    typedef EdgeIndexMap edge_index_t;

512
513
    typedef typename graph_traits<Graph>::vertex_descriptor vertex_t;
    typedef typename graph_traits<Graph>::edge_descriptor edge_t;
514
    typedef typename EdgeIndexMap::value_type index_t;
515

516
    RandomRewireStrategy(Graph& g, EdgeIndexMap edge_index,
517
518
                         rng_t& rng)
        : base_t(g, edge_index, rng)
519
    {
520
521
522
        tr1::bernoulli_distribution coin;
        bool invert = false;

523
524
        typename graph_traits<Graph>::edge_iterator e_i, e_i_end;
        for (tie(e_i, e_i_end) = edges(g); e_i != e_i_end; ++e_i)
525
        {
526
527
528
            // For undirected graphs, a random "direction" needs to be chosen,
            // and it needs to be uncorrelated with the implicit direction it
            // comes with.
529
530

            if (!is_directed::apply<Graph>::type::value)
531
532
533
                invert = coin(rng);

            _all_edges.push_back(make_pair(edge_index[*e_i], invert));
534
        }
535
        _all_edges2 = _all_edges;
536

537
    }
538

539
540
541
    void get_edges(const edge_t& e,
                   vector<pair<index_t,bool> >*& edges_source,
                   vector<pair<index_t,bool> >*& edges_target)
542
    {
543
        edges_source = &_all_edges;
544
        edges_target = &_all_edges2;
545
    }
546

547
private:
548
549
    vector<pair<index_t,bool> > _all_edges;
    vector<pair<index_t,bool> > _all_edges2;
550
};
551

552
553
554
555

// this will rewire the edges so that the (in,out) degree distributions and the
// (in,out)->(in,out) correlations will be the same, but all the rest is random
template <class Graph, class EdgeIndexMap>
556
557
558
class CorrelatedRewireStrategy:
    public RewireStrategyBase<Graph, EdgeIndexMap,
                              CorrelatedRewireStrategy<Graph, EdgeIndexMap> >
559
560
{
public:
561
562
563
564
565
566
567
    typedef RewireStrategyBase<Graph, EdgeIndexMap,
                               CorrelatedRewireStrategy<Graph, EdgeIndexMap> >
        base_t;

    typedef Graph graph_t;
    typedef EdgeIndexMap edge_index_t;

568
569
570
571
    typedef typename graph_traits<Graph>::vertex_descriptor vertex_t;
    typedef typename graph_traits<Graph>::edge_descriptor edge_t;
    typedef typename EdgeIndexMap::value_type index_t;

572
    CorrelatedRewireStrategy (Graph& g, EdgeIndexMap edge_index,
573
574
                              rng_t& rng)
        : base_t(g, edge_index, rng), _g(g)
575
    {
576
577
        typename graph_traits<Graph>::edge_iterator e_i, e_i_end;
        for (tie(e_i, e_i_end) = edges(g); e_i != e_i_end; ++e_i)
578
        {
579
580
581
            // For undirected graphs, there is no difference between source and
            // target, and each edge will appear _twice_ on the lists below,
            // once for each different ordering of source and target.
582

583
584
585
586
            _edges_by_source
                [make_pair(in_degreeS()(source(*e_i, _g), _g),
                           out_degree(source(*e_i, _g), _g))]
                .push_back(make_pair(edge_index[*e_i], false));
587

588
589
590
591
592
593
594
595
            _edges_by_target
                [make_pair(in_degreeS()(target(*e_i, _g), _g),
                           out_degree(target(*e_i, _g), _g))]
                .push_back(make_pair(edge_index[*e_i], false));

            if (!is_directed::apply<Graph>::type::value)
            {
                _edges_by_source
596
597
                    [make_pair(in_degreeS()(target(*e_i, _g), _g),
                               out_degree(target(*e_i, _g), _g))]
598
599
600
601
602
                    .push_back(make_pair(edge_index[*e_i], true));
                _edges_by_target
                    [make_pair(in_degreeS()(source(*e_i, _g), _g),
                               out_degree(source(*e_i, _g), _g))]
                    .push_back(make_pair(edge_index[*e_i], true));
603
604
605
            }
        }
    }
606

607
608
    void get_edges(const edge_t& e, vector<pair<index_t,bool> >*& edges_source,
                   vector<pair<index_t,bool> >*& edges_target)
609
    {
610
611
612
        pair<size_t, size_t> deg_source =
            make_pair(in_degreeS()(source(e, _g), _g),
                      out_degree(source(e, _g), _g));
613
614
        edges_source = &_edges_by_source[deg_source];

615
616

        pair<size_t, size_t> deg_target =
617
618
            make_pair(in_degreeS()(target(e, _g), _g),
                      out_degree(target(e, _g), _g));
619

620
621
622
        edges_target = &_edges_by_target[deg_target];
    }

623
private:
624
625
    typedef tr1::unordered_map<pair<size_t, size_t>,
                               vector<pair<index_t, bool> >,
626
                               hash<pair<size_t, size_t> > > edges_by_end_deg_t;
627
    edges_by_end_deg_t _edges_by_source, _edges_by_target;
628
629
630

protected:
    const Graph& _g;
631
632
633
634
635
};

} // graph_tool namespace

#endif // GRAPH_REWIRING_HH