graph_rewiring.hh 27.5 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
#include <boost/functional/hash.hpp>
Tiago Peixoto's avatar
Tiago Peixoto committed
24
25
26
27
28
29
#include <iostream>

#include <boost/multi_index_container.hpp>
#include <boost/multi_index/hashed_index.hpp>
#include <boost/multi_index/random_access_index.hpp>
#include <boost/multi_index/identity.hpp>
30
31
32

#include "graph.hh"
#include "graph_filtering.hh"
33
#include "graph_util.hh"
34
35
36
37
38

namespace graph_tool
{
using namespace std;
using namespace boost;
Tiago Peixoto's avatar
Tiago Peixoto committed
39
using namespace boost::multi_index;
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55

// 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;
}

56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
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);
}

Tiago Peixoto's avatar
Tiago Peixoto committed
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
// This will iterate over a random permutation of a random access sequence, by
// swapping the values of the sequence as it iterates
template <class RandomAccessIterator, class RNG,
          class RandomDist = tr1::uniform_int<size_t> >
class random_permutation_iterator : public
    std::iterator<input_iterator_tag, typename RandomAccessIterator::value_type>
{
public:
    random_permutation_iterator(RandomAccessIterator begin,
                                RandomAccessIterator end, RNG& rng)
        : _i(begin), _end(end), _rng(&rng)
    {
        if(_i != _end)
        {
            RandomDist random(0,  _end - _i - 1);
            std::iter_swap(_i, _i + random(*_rng));
        }
    }

    typename RandomAccessIterator::value_type operator*()
    {
        return *_i;
    }

    random_permutation_iterator& operator++()
    {
        ++_i;
        if(_i != _end)
        {
            RandomDist random(0,  _end - _i - 1);
            std::iter_swap(_i, _i + random(*_rng));
        }
        return *this;
    }

    bool operator==(const random_permutation_iterator& ri)
    {
        return _i == ri._i;
    }

    bool operator!=(const random_permutation_iterator& ri)
    {
        return _i != ri._i;
    }

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

private:
    RandomAccessIterator _i, _end;
    RNG* _rng;
};

133
134
// 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
Tiago Peixoto's avatar
Tiago Peixoto committed
135
struct swap_edge
136
137
{
    template <class Graph, class NewEdgeMap>
138
    static bool
Tiago Peixoto's avatar
Tiago Peixoto committed
139
140
141
142
    parallel_check_source
    (const typename graph_traits<Graph>::edge_descriptor& e,
     const pair<typename graph_traits<Graph>::edge_descriptor, bool>& se,
     NewEdgeMap edge_is_new, const Graph &g)
143
144
    {
        // We want to check that if we swap the source of 'e' with the source of
Tiago Peixoto's avatar
Tiago Peixoto committed
145
        // 'se', as such
146
        //
Tiago Peixoto's avatar
Tiago Peixoto committed
147
        //  (s)    -e--> (t)          (ns)   -e--> (t)
148
149
150
151
152
153
154
155
156
157
158
159
160
        //  (ns)   -se-> (se_t)   =>  (s)    -se-> (se_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
            se_t = target(se, g);      // source edge target

Tiago Peixoto's avatar
Tiago Peixoto committed
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
        if (is_adjacent_in_new(ns,  t, edge_is_new, g))
            return true; // e would clash with an existing (new) edge
        if (edge_is_new[se.first] && is_adjacent_in_new(s, se_t, edge_is_new,g))
            return true; // se would clash with an existing (new) edge
        return false; // the coast is clear - hooray!
    }

    template <class Graph, class NewEdgeMap>
    static bool
    parallel_check_target
    (const typename graph_traits<Graph>::edge_descriptor& e,
     const pair<typename graph_traits<Graph>::edge_descriptor, bool>& te,
     NewEdgeMap edge_is_new, const Graph &g)
    {
        // We want to check that if we swap the target of 'e' with the target of
        // 'te', as such
        //
        //  (s)    -e--> (t)          (s)    -e--> (nt)
        //  (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
            nt = target(te, g),        // new target
            te_s = source(te, g);      // target edge source
191

Tiago Peixoto's avatar
Tiago Peixoto committed
192
        if (is_adjacent_in_new(s,  nt, edge_is_new, g))
193
            return true; // e would clash with an existing (new) edge
194
        if (edge_is_new[te.first] && is_adjacent_in_new(te_s, t, edge_is_new,g))
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
            return true; // te 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)
        {
210
            if (edge_is_new[*e] && target(*e,g) == v)
211
212
213
214
215
216
                return true;
        }
        return false;
    }

    template <class Graph, class EdgeIndexMap, class EdgesType>
217
218
219
220
221
    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,
222
223
224
225
226
227
228
229
230
231
232
233
234
235
                    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
236
        if(se.first != te.first)
237
        {
238
            ne = add_edge(source(se, g), target(te, g), g).first;
239
            if(e != se.first)
240
            {
241
242
243
244
                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;
245
246
                edge_index[nse] = edge_index[se.first];
                remove_edge(se.first, g);
247
248
                edges[edge_index[nse]] = nse;
            }
249
            if(e != te.first)
250
            {
251
252
253
254
                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;
255
256
                edge_index[nte] = edge_index[te.first];
                remove_edge(te.first, g);
257
258
259
260
261
262
263
264
                edges[edge_index[nte]] = nte;
            }
            edge_index[ne] = edge_index[e];
            remove_edge(e, g);
            edges[edge_index[ne]] = ne;
        }
        else
        {
265
            if(e != se.first)
266
267
            {
                // se and te are the same. swapping indexes only.
268
269
                swap(edge_index[se.first], edge_index[e]);
                edges[edge_index[se.first]] = se.first;
270
271
272
273
274
275
                edges[edge_index[e]] = e;
            }
        }
    }
};

Tiago Peixoto's avatar
Tiago Peixoto committed
276
277
278
// used for verbose display
void print_progress(size_t current, size_t total, stringstream& str);

279
280
281
282
283
// main rewire loop
template <template <class Graph, class EdgeIndexMap> class RewireStrategy>
struct graph_rewire
{
    template <class Graph, class EdgeIndexMap>
284
    void operator()(Graph& g, EdgeIndexMap edge_index, rng_t& rng,
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
                    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
306
                throw ValueException("Self-loop detected. Can't rewire graph "
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
                                     "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
336
                throw ValueException("Parallel edge detected. Can't rewire "
337
338
339
340
341
                                     "graph without parallel edges if it "
                                     "already contains parallel edges!");
        }

        vector<edge_t> edges(num_edges(g));
Tiago Peixoto's avatar
Tiago Peixoto committed
342
        vector<size_t> edges_sequence;
343
344
        typename graph_traits<Graph>::edge_iterator e, e_end;
        for (tie(e, e_end) = boost::edges(g); e != e_end; ++e)
345
        {
346
347
348
            if (edge_index[*e] >= edges.size())
                edges.resize(edge_index[*e] + 1);
            edges[edge_index[*e]] = *e;
Tiago Peixoto's avatar
Tiago Peixoto committed
349
            edges_sequence.push_back(edge_index[*e]);
350
351
        }

Tiago Peixoto's avatar
Tiago Peixoto committed
352
        RewireStrategy<Graph, EdgeIndexMap> rewire(g, edge_index, edges, rng);
353

Tiago Peixoto's avatar
Tiago Peixoto committed
354
355
356
        typedef random_permutation_iterator<typename vector<size_t>::iterator,
                                            rng_t>
            random_edge_iter;
357

Tiago Peixoto's avatar
Tiago Peixoto committed
358
359
360
        random_edge_iter
            ei(edges_sequence.begin(), edges_sequence.end(), rng),
            ei_end(edges_sequence.end(), edges_sequence.end(), rng);
361

Tiago Peixoto's avatar
Tiago Peixoto committed
362
363
364
365
        stringstream str;
        //size_t count = 0;
        // for each edge simultaneously rewire its source or target
        for (; ei != ei_end; ++ei)
366
        {
Tiago Peixoto's avatar
Tiago Peixoto committed
367
368
            //print_progress(count++, edges_sequence.size(), str);
            rewire(edges[*ei], self_loops, parallel_edges);
369
        }
370
    }
371
372
};

Tiago Peixoto's avatar
Tiago Peixoto committed
373

374
375
376
377
378
379
380
381
382
// 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;

Tiago Peixoto's avatar
Tiago Peixoto committed
383
384
385
386
    ErdosRewireStrategy(Graph& g, EdgeIndexMap edge_index,
                        vector<edge_t>& edges, rng_t& rng)
        : _g(g), _edge_index(edge_index), _edges(edges),
          _vertices(HardNumVertices()(g)), _rng(rng)
387
388
389
390
391
392
393
    {
        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;
    }

Tiago Peixoto's avatar
Tiago Peixoto committed
394
    void operator()(const edge_t& e, bool self_loops, bool parallel_edges)
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
    {
        //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 &&
Tiago Peixoto's avatar
Tiago Peixoto committed
412
                swap_edge::is_adjacent_in_new(s, t, _edge_is_new, _g))
413
414
415
416
                continue;  // reject parallel edges if not allowed
            break;
        }
        edge_t ne = add_edge(s, t, _g).first;
Tiago Peixoto's avatar
Tiago Peixoto committed
417
        _edges[_edge_index[e]] = ne;
418
        remove_edge(e, _g);
Tiago Peixoto's avatar
Tiago Peixoto committed
419
420
421
        if (_edge_index[ne] >= _edges.size())
            _edges.resize(_edge_index[ne] + 1);
        _edges[_edge_index[ne]] = ne;
422
423
424
425
426
427
        _edge_is_new[ne] = true;
    }

private:
    Graph& _g;
    EdgeIndexMap _edge_index;
Tiago Peixoto's avatar
Tiago Peixoto committed
428
    vector<edge_t>& _edges;
429
    vector<typename graph_traits<Graph>::vertex_descriptor> _vertices;
430
    checked_vector_property_map<bool, EdgeIndexMap> _edge_is_new;
431
432
433
434
    rng_t& _rng;
};


435
436
437
// 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.
438
439
template <class Graph, class EdgeIndexMap, class RewireStrategy>
class RewireStrategyBase
440
441
442
{
public:
    typedef typename graph_traits<Graph>::edge_descriptor edge_t;
443
444
    typedef typename graph_traits<Graph>::edge_descriptor vertex_t;

445
446
    typedef typename EdgeIndexMap::value_type index_t;

Tiago Peixoto's avatar
Tiago Peixoto committed
447
448
449
450
451
452
453
454
455
    RewireStrategyBase(Graph& g, EdgeIndexMap edge_index, vector<edge_t>& edges,
                       rng_t& rng)
        : _g(g), _edge_index(edge_index), _edges(edges),
          _edge_is_new(edge_index, 1), _rng(rng)
    {
        typename graph_traits<Graph>::edge_iterator e_i, e_i_end;
        for (tie(e_i, e_i_end) = boost::edges(g); e_i != e_i_end; ++e_i)
            _edge_is_new.get_checked()[*e_i] == true;
    }
456

Tiago Peixoto's avatar
Tiago Peixoto committed
457
    void operator()(const edge_t& e, bool self_loops, bool parallel_edges)
458
    {
Tiago Peixoto's avatar
Tiago Peixoto committed
459
460
        RewireStrategy& self = *static_cast<RewireStrategy*>(this);
        tr1::bernoulli_distribution coin;
461
462
463

        //try randomly drawn pairs of edges until one satisfies all the
        //consistency checks
464
        pair<edge_t,bool> es, et;
465

Tiago Peixoto's avatar
Tiago Peixoto committed
466
467
        _edge_is_new[e] = false;
        while (true)
468
        {
Tiago Peixoto's avatar
Tiago Peixoto committed
469
470
            es.first = et.first = e;
            es.second = et.second = false;
471

Tiago Peixoto's avatar
Tiago Peixoto committed
472
            if (coin(_rng))
473
            {
Tiago Peixoto's avatar
Tiago Peixoto committed
474
475
                // rewire source
                es = self.get_source_edge(e);
476

Tiago Peixoto's avatar
Tiago Peixoto committed
477
478
479
480
481
482
483
484
485
486
487
488
489
490
                if(!self_loops) // reject self-loops if not allowed
                {
                    if((source(e, _g) == target(es, _g)) ||
                       (target(e, _g) == source(es, _g)))
                        continue;
                }
                if (!parallel_edges) // reject parallel edges if not allowed
                {
                    if (swap_edge::parallel_check_source(e, es, _edge_is_new,
                                                         _g))
                        continue;
                }
            }
            else
491
            {
Tiago Peixoto's avatar
Tiago Peixoto committed
492
493
                // rewire target
                et = self.get_target_edge(e);
494

495
496
                if (!self_loops) // reject self-loops if not allowed
                {
Tiago Peixoto's avatar
Tiago Peixoto committed
497
498
                    if((source(e, _g) == target(et, _g)) ||
                       (target(e, _g) == source(et, _g)))
499
500
501
502
                        continue;
                }
                if (!parallel_edges) // reject parallel edges if not allowed
                {
Tiago Peixoto's avatar
Tiago Peixoto committed
503
504
                    if (swap_edge::parallel_check_target(e, et, _edge_is_new,
                                                         _g))
505
506
507
                        continue;
                }
            }
Tiago Peixoto's avatar
Tiago Peixoto committed
508
            break;
509
        }
Tiago Peixoto's avatar
Tiago Peixoto committed
510
511

        swap_edge()(e, es, et, _edges, _edge_index, _g);
512
        _edge_is_new[e] = true;
Tiago Peixoto's avatar
Tiago Peixoto committed
513
514
515
516
517
518
519
520

        const edge_t* elist[3] = {&e, &es.first, &et.first};
        for (size_t i = 0; i < 3; ++i)
        {
            self.update_edge(*elist[i], false);
            self.update_edge(*elist[i], true);
        }

521
522
    }

Tiago Peixoto's avatar
Tiago Peixoto committed
523
protected:
524
525
    Graph& _g;
    EdgeIndexMap _edge_index;
Tiago Peixoto's avatar
Tiago Peixoto committed
526
527
    vector<edge_t>& _edges;
    unchecked_vector_property_map<bool, EdgeIndexMap> _edge_is_new;
528
    rng_t& _rng;
529
530
};

531
532
// this will rewire the edges so that the combined (in, out) degree distribution
// will be the same, but all the rest is random
533
template <class Graph, class EdgeIndexMap>
534
535
536
class RandomRewireStrategy:
    public RewireStrategyBase<Graph, EdgeIndexMap,
                              RandomRewireStrategy<Graph, EdgeIndexMap> >
537
538
{
public:
539
540
541
542
543
544
545
    typedef RewireStrategyBase<Graph, EdgeIndexMap,
                               RandomRewireStrategy<Graph, EdgeIndexMap> >
        base_t;

    typedef Graph graph_t;
    typedef EdgeIndexMap edge_index_t;

546
547
    typedef typename graph_traits<Graph>::vertex_descriptor vertex_t;
    typedef typename graph_traits<Graph>::edge_descriptor edge_t;
548
    typedef typename EdgeIndexMap::value_type index_t;
549

Tiago Peixoto's avatar
Tiago Peixoto committed
550
551
552
553
554
555
556
557
558
559
    struct hash_index {};
    struct random_index {};

    typedef multi_index_container<
        pair<index_t,bool>,
        indexed_by<
            hashed_unique<tag<hash_index>, identity<pair<index_t,bool> > >,
            random_access<tag<random_index> > > >
        edge_set_t;

560
    RandomRewireStrategy(Graph& g, EdgeIndexMap edge_index,
Tiago Peixoto's avatar
Tiago Peixoto committed
561
562
                         vector<edge_t>& edges, rng_t& rng)
        : base_t(g, edge_index, edges, rng), _g(g)
563
    {
Tiago Peixoto's avatar
Tiago Peixoto committed
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
        size_t E = HardNumEdges()(_g);

        typename graph_traits<Graph>::vertex_iterator v_i, v_i_end;
        for (tie(v_i, v_i_end) = vertices(g); v_i != v_i_end; ++v_i)
        {
            if (out_degree(*v_i, _g) > E/3)
            {
                _hub_non_adjacent[*v_i] = edge_set_t();
                if(!is_directed::apply<Graph>::type::value)
                    _hub_non_incident[*v_i] = edge_set_t();
            }

            if (in_degreeS()(*v_i, _g) > E/3)
                _hub_non_incident[*v_i] = edge_set_t();
        }
579

580
        typename graph_traits<Graph>::edge_iterator e_i, e_i_end;
Tiago Peixoto's avatar
Tiago Peixoto committed
581
        for (tie(e_i, e_i_end) = boost::edges(g); e_i != e_i_end; ++e_i)
582
        {
Tiago Peixoto's avatar
Tiago Peixoto committed
583
584
585
            // For undirected graphs, there is no difference between source and
            // target, and each edge will appear _twice_ in the list below,
            // once for each different ordering of source and target.
586

Tiago Peixoto's avatar
Tiago Peixoto committed
587
            _all_edges.push_back(make_pair(edge_index[*e_i], false));
588
            if (!is_directed::apply<Graph>::type::value)
Tiago Peixoto's avatar
Tiago Peixoto committed
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
                _all_edges.push_back(make_pair(edge_index[*e_i], true));
            update_edge(*e_i, true);
        }

        for (typeof(_hub_non_adjacent.begin()) iter = _hub_non_adjacent.begin();
             iter != _hub_non_adjacent.end(); ++iter)
            cout << iter->first << " " << iter->second.size() << endl;
        for (typeof(_hub_non_incident.begin()) iter = _hub_non_incident.begin();
             iter != _hub_non_incident.end(); ++iter)
            cout << iter->first << " " << iter->second.size() << endl;
    }

    pair<edge_t,bool> get_edge(const edge_t& e, bool src)
    {
        tr1::unordered_map<size_t, edge_set_t>& hub =
            (src || !is_directed::apply<Graph>::type::value) ?
            _hub_non_adjacent : _hub_non_incident;
606

Tiago Peixoto's avatar
Tiago Peixoto committed
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
        typeof(hub.begin()) iter =
            hub.find(src ? target(e,_g) : source(e,_g));
        if (iter != hub.end())
        {
            edge_set_t& eset = iter->second;

            tr1::uniform_int<> sample(0, eset.size());
            size_t i = sample(base_t::_rng);
            if (i == eset.size()) // no rewire option
            {
                return make_pair(e, false);
            }
            else
            {
                pair<index_t,bool> ep =
                    *(eset.get<random_index>().begin() + i);
                return  make_pair(base_t::_edges[ep.first], ep.second);
            }
625
626
        }

Tiago Peixoto's avatar
Tiago Peixoto committed
627
628
629
        tr1::uniform_int<> sample(0, _all_edges.size()-1);
        pair<index_t,bool> ep = _all_edges[sample(base_t::_rng)];
        return make_pair(base_t::_edges[ep.first], ep.second);
630
    }
631

Tiago Peixoto's avatar
Tiago Peixoto committed
632
633

    pair<edge_t,bool> get_source_edge(const edge_t& e)
634
    {
Tiago Peixoto's avatar
Tiago Peixoto committed
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
        return get_edge(e, true);
    }

    pair<edge_t,bool> get_target_edge(const edge_t& e)
    {
        return get_edge(e, false);
    }

    void update_edge_list(index_t v, const pair<edge_t, bool>& e,
                          edge_set_t& edge_set, bool src, bool insert)
    {
        pair<index_t,bool> ep;
        if (!base_t::_edge_is_new[e.first] || !insert ||
            (src && source(e.first, _g) != v) ||
            (!src && target(e.first, _g) != v))
        {
            ep = make_pair(_edge_index[e.first], e.second);
            if (insert)
                edge_set.get<hash_index>().insert(ep);
            else
                edge_set.get<hash_index>().erase(ep);
        }
    }

    void update_edge(const edge_t& e, bool insert)
    {
        for (typeof(_hub_non_adjacent.begin()) iter = _hub_non_adjacent.begin();
             iter != _hub_non_adjacent.end(); ++iter)
        {
            update_edge_list(iter->first, make_pair(e, false), iter->second,
                             true, insert);
            if (!is_directed::apply<Graph>::type::value)
                update_edge_list(iter->first, make_pair(e, true), iter->second,
                                 true, insert);
        }

        for (typeof(_hub_non_incident.begin()) iter = _hub_non_incident.begin();
             iter != _hub_non_incident.end(); ++iter)
        {
            update_edge_list(iter->first, make_pair(e, false), iter->second,
                             false, insert);
            if (!is_directed::apply<Graph>::type::value)
                update_edge_list(iter->first, make_pair(e, true), iter->second,
                                 false, insert);
        }
680
    }
681

682
private:
Tiago Peixoto's avatar
Tiago Peixoto committed
683
684
    Graph& _g;
    EdgeIndexMap _edge_index;
685
    vector<pair<index_t,bool> > _all_edges;
Tiago Peixoto's avatar
Tiago Peixoto committed
686
687
688

    tr1::unordered_map<size_t, edge_set_t>  _hub_non_incident;
    tr1::unordered_map<size_t, edge_set_t>  _hub_non_adjacent;
689
};
690

691
692
693
694

// 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>
695
696
697
class CorrelatedRewireStrategy:
    public RewireStrategyBase<Graph, EdgeIndexMap,
                              CorrelatedRewireStrategy<Graph, EdgeIndexMap> >
698
699
{
public:
700
701
702
703
704
705
706
    typedef RewireStrategyBase<Graph, EdgeIndexMap,
                               CorrelatedRewireStrategy<Graph, EdgeIndexMap> >
        base_t;

    typedef Graph graph_t;
    typedef EdgeIndexMap edge_index_t;

707
708
709
710
    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;

711
    CorrelatedRewireStrategy (Graph& g, EdgeIndexMap edge_index,
Tiago Peixoto's avatar
Tiago Peixoto committed
712
713
                              vector<edge_t>& edges, rng_t& rng)
        : base_t(g, edge_index, edges, rng), _g(g)
714
    {
715
        typename graph_traits<Graph>::edge_iterator e_i, e_i_end;
Tiago Peixoto's avatar
Tiago Peixoto committed
716
        for (tie(e_i, e_i_end) = boost::edges(g); e_i != e_i_end; ++e_i)
717
        {
718
            // For undirected graphs, there is no difference between source and
Tiago Peixoto's avatar
Tiago Peixoto committed
719
            // target, and each edge will appear _twice_ in the lists below,
720
            // once for each different ordering of source and target.
721

722
723
724
725
            _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));
726

727
728
729
730
731
732
733
734
            _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
735
736
                    [make_pair(in_degreeS()(target(*e_i, _g), _g),
                               out_degree(target(*e_i, _g), _g))]
737
738
739
740
741
                    .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));
742
743
744
            }
        }
    }
745

Tiago Peixoto's avatar
Tiago Peixoto committed
746
    pair<edge_t,bool> get_edge(const edge_t& e, bool src)
747
    {
Tiago Peixoto's avatar
Tiago Peixoto committed
748
749
750
751
752
753
754
755
756
        pair<size_t, size_t> deg =
            make_pair(in_degreeS()(src ? source(e, _g) : target(e, _g), _g),
                      out_degree(src ? source(e, _g) : target(e, _g), _g));
        edges_by_end_deg_t& edges = src ? _edges_by_source : _edges_by_target;
        typename edges_by_end_deg_t::mapped_type& elist = edges[deg];
        tr1::uniform_int<> sample(0, elist.size()-1);
        pair<index_t,bool> ep = elist[sample(base_t::_rng)];
        return make_pair(base_t::_edges[ep.first], ep.second);
    }
757

Tiago Peixoto's avatar
Tiago Peixoto committed
758
759
760
761
    pair<edge_t,bool> get_source_edge(const edge_t& e)
    {
        return get_edge(e, true);
    }
762

Tiago Peixoto's avatar
Tiago Peixoto committed
763
764
765
    pair<edge_t,bool> get_target_edge(const edge_t& e)
    {
        return get_edge(e, false);
766
767
    }

Tiago Peixoto's avatar
Tiago Peixoto committed
768
769
    void update_edge(const edge_t& e, bool insert) {}

770
private:
771
772
    typedef tr1::unordered_map<pair<size_t, size_t>,
                               vector<pair<index_t, bool> >,
773
                               hash<pair<size_t, size_t> > > edges_by_end_deg_t;
774
    edges_by_end_deg_t _edges_by_source, _edges_by_target;
775
776
777

protected:
    const Graph& _g;
778
779
780
781
782
};

} // graph_tool namespace

#endif // GRAPH_REWIRING_HH