graph_rewiring.hh 20.5 KB
Newer Older
1
2
// graph-tool -- a general graph modification and manipulation thingy
//
Tiago Peixoto's avatar
Tiago Peixoto committed
3
// Copyright (C) 2006-2013 Tiago de Paula Peixoto <tiago@skewed.de>
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
//
// 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

21
22
#include "tr1_include.hh"
#include TR1_HEADER(unordered_set)
23
#include <boost/functional/hash.hpp>
24

Tiago Peixoto's avatar
Tiago Peixoto committed
25
26
#include <iostream>

27
28
#include "graph.hh"
#include "graph_filtering.hh"
29
#include "graph_util.hh"
30
#include "sampler.hh"
31

32
33
#include "random.hh"

34
35
36
37
38
namespace graph_tool
{
using namespace std;
using namespace boost;

39
40
template <class Graph>
typename graph_traits<Graph>::vertex_descriptor
41
42
source(const pair<size_t, bool>& e,
       const vector<typename graph_traits<Graph>::edge_descriptor>& edges,
43
44
45
       const Graph& g)
{
    if (e.second)
46
        return target(edges[e.first], g);
47
    else
48
        return source(edges[e.first], g);
49
50
51
52
}

template <class Graph>
typename graph_traits<Graph>::vertex_descriptor
53
54
target(const pair<size_t, bool>& e,
       const vector<typename graph_traits<Graph>::edge_descriptor>& edges,
55
56
57
       const Graph& g)
{
    if (e.second)
58
        return source(edges[e.first], g);
59
    else
60
        return target(edges[e.first], g);
61
62
}

Tiago Peixoto's avatar
Tiago Peixoto committed
63

64
65
// 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
66
struct swap_edge
67
{
68
    template <class Graph>
Tiago Peixoto's avatar
Tiago Peixoto committed
69
    static bool
70
71
72
    parallel_check_target (size_t e, const pair<size_t, bool>& te,
                           vector<typename graph_traits<Graph>::edge_descriptor>& edges,
                           const Graph &g)
Tiago Peixoto's avatar
Tiago Peixoto committed
73
74
75
76
77
78
79
    {
        // 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)
        //
80
        // no parallel edges are introduced.
Tiago Peixoto's avatar
Tiago Peixoto committed
81
82

        typename graph_traits<Graph>::vertex_descriptor
83
84
85
86
            s = source(edges[e], g),          // current source
            t = target(edges[e], g),          // current target
            nt = target(te, edges, g),        // new target
            te_s = source(te, edges, g);      // target edge source
87

88
        if (is_adjacent(s,  nt, g))
89
            return true; // e would clash with an existing edge
90
        if (is_adjacent(te_s, t, g))
91
            return true; // te would clash with an existing edge
92
93
94
        return false; // the coast is clear - hooray!
    }

95
96
97
98
    template <class Graph, class EdgeIndexMap>
    static void swap_target
        (size_t e, const pair<size_t, bool>& te,
         vector<typename graph_traits<Graph>::edge_descriptor>& edges,
99
         EdgeIndexMap, Graph& g)
100
    {
101
102
        // swap the source of the edge 'e' with the source of edge 'se', as
        // such:
103
        //
104
105
        //  (s)    -e--> (t)          (s)    -e--> (nt)
        //  (te_s) -te-> (nt)   =>    (te_s) -te-> (t)
106

107
108
        if (e == te.first)
            return;
109

110
111
        // new edges which will replace the old ones
        typename graph_traits<Graph>::edge_descriptor ne, nte;
112
113
114
115
116
117
        typename graph_traits<Graph>::vertex_descriptor
            s_e = source(edges[e], g),
            t_e = target(edges[e], g),
            s_te = source(te, edges, g),
            t_te = target(te, edges, g);
        remove_edge(edges[e], g);
118
119
        remove_edge(edges[te.first], g);

120
        ne = add_edge(s_e, t_te, g).first;
121
        edges[e] = ne;
122
123
124
125
126
        if (!te.second)
            nte = add_edge(s_te, t_e, g).first;
        else // keep invertedness (only for undirected graphs)
            nte = add_edge(t_e, s_te,  g).first;
        edges[te.first] = nte;
127
    }
128

129
130
};

Tiago Peixoto's avatar
Tiago Peixoto committed
131
// used for verbose display
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
void print_progress(size_t i, size_t n_iter, size_t current, size_t total,
                    stringstream& str)
{
    size_t atom = (total > 200) ? total / 100 : 1;
    if ( ( (current+1) % atom == 0) || (current + 1) == total)
    {
        size_t size = str.str().length();
        for (size_t j = 0; j < str.str().length(); ++j)
            cout << "\b";
        str.str("");
        str << "(" << i + 1 << " / " << n_iter << ") "
            << current + 1 << " of " << total << " ("
            << (current + 1) * 100 / total << "%)";
        for (int j = 0; j < int(size - str.str().length()); ++j)
            str << " ";
        cout << str.str() << flush;
    }
}

151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
//select blocks based on in/out degrees
class DegreeBlock
{
public:
    typedef pair<size_t, size_t> block_t;

    template <class Graph>
    block_t get_block(typename graph_traits<Graph>::vertex_descriptor v,
                      const Graph& g) const
    {
        return make_pair(in_degreeS()(v, g), out_degree(v, g));
    }
};

//select blocks based on property map
template <class PropertyMap>
class PropertyBlock
{
public:
    typedef typename property_traits<PropertyMap>::value_type block_t;

    PropertyBlock(PropertyMap p): _p(p) {}

    template <class Graph>
    block_t get_block(typename graph_traits<Graph>::vertex_descriptor v,
176
                      const Graph&) const
177
178
179
180
181
182
183
    {
        return get(_p, v);
    }

private:
    PropertyMap _p;
};
Tiago Peixoto's avatar
Tiago Peixoto committed
184

185
// main rewire loop
186
187
template <template <class Graph, class EdgeIndexMap, class CorrProb,
                    class BlockDeg>
188
          class RewireStrategy>
189
190
struct graph_rewire
{
191
192
193

    template <class Graph, class EdgeIndexMap, class CorrProb,
              class BlockDeg>
194
    void operator()(Graph& g, EdgeIndexMap edge_index, CorrProb corr_prob,
195
                    bool self_loops, bool parallel_edges,
196
197
198
                    pair<size_t, bool> iter_sweep,
                    pair<bool, bool> cache_verbose,
                    size_t& pcount, rng_t& rng, BlockDeg bd)
199
        const
200
201
202
    {
        typedef typename graph_traits<Graph>::vertex_descriptor vertex_t;
        typedef typename graph_traits<Graph>::edge_descriptor edge_t;
203
204
        bool cache = cache_verbose.first;
        bool verbose = cache_verbose.second;
205

206
        vector<edge_t> edges;
Tiago Peixoto's avatar
Tiago Peixoto committed
207
208
        vector<size_t> edge_pos;
        typedef random_permutation_iterator<typename vector<size_t>::iterator,
209
210
211
                                            rng_t>
            random_edge_iter;

212
213
        typename graph_traits<Graph>::edge_iterator e, e_end;
        for (tie(e, e_end) = boost::edges(g); e != e_end; ++e)
Tiago Peixoto's avatar
Tiago Peixoto committed
214
        {
215
            edges.push_back(*e);
Tiago Peixoto's avatar
Tiago Peixoto committed
216
217
            edge_pos.push_back(edge_pos.size());
        }
218

219
        RewireStrategy<Graph, EdgeIndexMap, CorrProb, BlockDeg>
220
            rewire(g, edge_index, edges, corr_prob, bd, cache, rng);
221

222
223
224
225
        size_t niter;
        bool no_sweep;
        tie(niter, no_sweep) = iter_sweep;
        pcount = 0;
226
227
        if (verbose)
            cout << "rewiring edges: ";
Tiago Peixoto's avatar
Tiago Peixoto committed
228
        stringstream str;
229
        for (size_t i = 0; i < niter; ++i)
230
        {
Tiago Peixoto's avatar
Tiago Peixoto committed
231
            random_edge_iter
Tiago Peixoto's avatar
Tiago Peixoto committed
232
233
                ei_begin(edge_pos.begin(), edge_pos.end(), rng),
                ei_end(edge_pos.end(), edge_pos.end(), rng);
234
235

            // for each edge rewire its source or target
236
            for (random_edge_iter ei = ei_begin; ei != ei_end; ++ei)
237
            {
238
                size_t e_pos = ei - ei_begin;
239
                if (verbose)
240
241
                    print_progress(i, niter, e_pos, no_sweep ? 1 : edges.size(),
                                   str);
Tiago Peixoto's avatar
Tiago Peixoto committed
242
                bool success = rewire(*ei, self_loops, parallel_edges);
243
244
245
246
247
248
                if (!success)
                    ++pcount;

                if (no_sweep)
                    break;
            }
249
        }
250
251
        if (verbose)
            cout << endl;
252
    }
253
254
255
256

    template <class Graph, class EdgeIndexMap, class CorrProb>
    void operator()(Graph& g, EdgeIndexMap edge_index, CorrProb corr_prob,
                    bool self_loops, bool parallel_edges,
257
258
259
                    pair<size_t, bool> iter_sweep,
                    pair<bool, bool> cache_verbose,
                    size_t& pcount, rng_t& rng)
260
261
262
        const
    {
        operator()(g, edge_index, corr_prob, self_loops, parallel_edges,
263
                   iter_sweep, cache_verbose, pcount, rng, DegreeBlock());
264
    }
265
266
};

Tiago Peixoto's avatar
Tiago Peixoto committed
267

268
269
// this will rewire the edges so that the resulting graph will be entirely
// random (i.e. Erdos-Renyi)
270
template <class Graph, class EdgeIndexMap, class CorrProb, class BlockDeg>
271
272
273
274
275
276
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
277
    ErdosRewireStrategy(Graph& g, EdgeIndexMap edge_index,
278
                        vector<edge_t>& edges, CorrProb, BlockDeg,
279
                        bool, rng_t& rng)
Tiago Peixoto's avatar
Tiago Peixoto committed
280
281
        : _g(g), _edge_index(edge_index), _edges(edges),
          _vertices(HardNumVertices()(g)), _rng(rng)
282
283
284
285
286
287
288
    {
        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;
    }

289
    bool operator()(size_t ei, bool self_loops, bool parallel_edges)
290
    {
291
        //try randomly drawn pairs of vertices
292
293
294
295
        typedef random_permutation_iterator
            <typename graph_traits<Graph>::vertex_iterator, rng_t>
            random_vertex_iter;

296
        tr1::uniform_int<size_t> sample(0, _vertices.size() - 1);
297
298
299
300
301
302
        typename graph_traits<Graph>::vertex_descriptor s, t;
        while (true)
        {
            s = sample(_rng);
            t = sample(_rng);

303
304
            // reject self-loops if not allowed
            if(s == t && !self_loops)
305
306
307
                continue;
            break;
        }
308
309
310
311
312

        // reject parallel edges if not allowed
        if (!parallel_edges && is_adjacent(s, t, _g))
            return false;

313
        remove_edge(_edges[ei], _g);
314
        edge_t ne = add_edge(s, t, _g).first;
315
        _edges[ei] = ne;
316
317

        return true;
318
319
320
321
322
    }

private:
    Graph& _g;
    EdgeIndexMap _edge_index;
Tiago Peixoto's avatar
Tiago Peixoto committed
323
    vector<edge_t>& _edges;
324
325
326
327
328
    vector<typename graph_traits<Graph>::vertex_descriptor> _vertices;
    rng_t& _rng;
};


329
330
331
// 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.
332
333
template <class Graph, class EdgeIndexMap, class RewireStrategy>
class RewireStrategyBase
334
335
336
{
public:
    typedef typename graph_traits<Graph>::edge_descriptor edge_t;
337
338
    typedef typename graph_traits<Graph>::edge_descriptor vertex_t;

339
340
    typedef typename EdgeIndexMap::value_type index_t;

Tiago Peixoto's avatar
Tiago Peixoto committed
341
342
    RewireStrategyBase(Graph& g, EdgeIndexMap edge_index, vector<edge_t>& edges,
                       rng_t& rng)
343
        : _g(g), _edge_index(edge_index), _edges(edges),  _rng(rng) {}
344

345
    bool operator()(size_t ei, bool self_loops, bool parallel_edges)
346
    {
Tiago Peixoto's avatar
Tiago Peixoto committed
347
        RewireStrategy& self = *static_cast<RewireStrategy*>(this);
348

349
350
        // try randomly drawn pairs of edges and check if they satisfy all the
        // consistency checks
351

352
        // rewire target
353
        pair<size_t, bool> et = self.get_target_edge(ei);
354
355

        if (!self_loops) // reject self-loops if not allowed
356
        {
357
358
            if((source(_edges[ei], _g) == target(et, _edges, _g)) ||
               (target(_edges[ei], _g) == source(et, _edges, _g)))
359
360
                return false;
        }
361

362
        // reject parallel edges if not allowed
363
        if (!parallel_edges && (et.first != ei))
364
        {
365
            if (swap_edge::parallel_check_target(ei, et, _edges, _g))
366
                return false;
367
        }
368

369
        if (ei != et.first)
370
        {
371
            self.update_edge(ei, false);
372
            self.update_edge(et.first, false);
Tiago Peixoto's avatar
Tiago Peixoto committed
373

374
            swap_edge::swap_target(ei, et, _edges, _edge_index, _g);
375

376
            self.update_edge(ei, true);
377
378
            self.update_edge(et.first, true);
        }
379
380
381
382
        else
        {
            return false;
        }
Tiago Peixoto's avatar
Tiago Peixoto committed
383

384
        return true;
385
386
    }

Tiago Peixoto's avatar
Tiago Peixoto committed
387
protected:
388
389
    Graph& _g;
    EdgeIndexMap _edge_index;
Tiago Peixoto's avatar
Tiago Peixoto committed
390
    vector<edge_t>& _edges;
391
    rng_t& _rng;
392
393
};

394
395
// this will rewire the edges so that the combined (in, out) degree distribution
// will be the same, but all the rest is random
396
template <class Graph, class EdgeIndexMap, class CorrProb, class BlockDeg>
397
398
class RandomRewireStrategy:
    public RewireStrategyBase<Graph, EdgeIndexMap,
399
                              RandomRewireStrategy<Graph, EdgeIndexMap,
400
                                                   CorrProb, BlockDeg> >
401
402
{
public:
403
    typedef RewireStrategyBase<Graph, EdgeIndexMap,
404
                               RandomRewireStrategy<Graph, EdgeIndexMap,
405
                                                    CorrProb, BlockDeg> >
406
407
408
409
410
        base_t;

    typedef Graph graph_t;
    typedef EdgeIndexMap edge_index_t;

411
412
    typedef typename graph_traits<Graph>::vertex_descriptor vertex_t;
    typedef typename graph_traits<Graph>::edge_descriptor edge_t;
413
    typedef typename EdgeIndexMap::value_type index_t;
414

Tiago Peixoto's avatar
Tiago Peixoto committed
415
416
417
    struct hash_index {};
    struct random_index {};

418
    RandomRewireStrategy(Graph& g, EdgeIndexMap edge_index,
419
                         vector<edge_t>& edges, CorrProb, BlockDeg,
420
                         bool, rng_t& rng)
421
        : base_t(g, edge_index, edges, rng), _g(g) {}
Tiago Peixoto's avatar
Tiago Peixoto committed
422

423
    pair<size_t,bool> get_target_edge(size_t)
Tiago Peixoto's avatar
Tiago Peixoto committed
424
    {
425
426
427
        tr1::uniform_int<> sample(0, base_t::_edges.size() - 1);
        pair<size_t, bool> et = make_pair(sample(base_t::_rng), false);
        if (!is_directed::apply<Graph>::type::value)
Tiago Peixoto's avatar
Tiago Peixoto committed
428
        {
429
430
            tr1::bernoulli_distribution coin(0.5);
            et.second = coin(base_t::_rng);
Tiago Peixoto's avatar
Tiago Peixoto committed
431
        }
432
        return et;
Tiago Peixoto's avatar
Tiago Peixoto committed
433
434
    }

435
    void update_edge(size_t, bool) {}
436

437
private:
Tiago Peixoto's avatar
Tiago Peixoto committed
438
439
    Graph& _g;
    EdgeIndexMap _edge_index;
440
};
441

442
443
444

// 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
445
template <class Graph, class EdgeIndexMap, class CorrProb, class BlockDeg>
446
447
class CorrelatedRewireStrategy:
    public RewireStrategyBase<Graph, EdgeIndexMap,
448
                              CorrelatedRewireStrategy<Graph, EdgeIndexMap,
449
                                                       CorrProb, BlockDeg> >
450
451
{
public:
452
    typedef RewireStrategyBase<Graph, EdgeIndexMap,
453
                               CorrelatedRewireStrategy<Graph, EdgeIndexMap,
454
                                                        CorrProb, BlockDeg> >
455
456
457
458
        base_t;

    typedef Graph graph_t;

459
460
461
    typedef typename graph_traits<Graph>::vertex_descriptor vertex_t;
    typedef typename graph_traits<Graph>::edge_descriptor edge_t;

462
463
    CorrelatedRewireStrategy(Graph& g, EdgeIndexMap edge_index,
                             vector<edge_t>& edges, CorrProb, BlockDeg,
464
                             bool, rng_t& rng)
Tiago Peixoto's avatar
Tiago Peixoto committed
465
        : base_t(g, edge_index, edges, rng), _g(g)
466
    {
467
        for (size_t ei = 0; ei < base_t::_edges.size(); ++ei)
468
        {
469
            // For undirected graphs, there is no difference between source and
Tiago Peixoto's avatar
Tiago Peixoto committed
470
            // target, and each edge will appear _twice_ in the list below,
471
            // once for each different ordering of source and target.
472
            edge_t& e = base_t::_edges[ei];
473

Tiago Peixoto's avatar
Tiago Peixoto committed
474
475
476
            vertex_t t = target(e, _g);
            deg_t tdeg = make_pair(in_degreeS()(t, _g), out_degree(t, _g));
            _edges_by_target[tdeg].push_back(make_pair(ei, false));
477
478
479

            if (!is_directed::apply<Graph>::type::value)
            {
Tiago Peixoto's avatar
Tiago Peixoto committed
480
481
482
                t = source(e, _g);
                tdeg = make_pair(in_degreeS()(t, _g), out_degree(t, _g));
                _edges_by_target[tdeg].push_back(make_pair(ei, true));
483
484
485
            }
        }
    }
486

487
    pair<size_t,bool> get_target_edge(size_t ei)
488
    {
Tiago Peixoto's avatar
Tiago Peixoto committed
489
490
491
492
493
        edge_t& e = base_t::_edges[ei];
        vertex_t t = target(e, _g);
        deg_t tdeg = make_pair(in_degreeS()(t, _g), out_degree(t, _g));
        typename edges_by_end_deg_t::mapped_type& elist =
            _edges_by_target[tdeg];
494
        tr1::uniform_int<> sample(0, elist.size() - 1);
Tiago Peixoto's avatar
Tiago Peixoto committed
495

496
        return elist[sample(base_t::_rng)];
Tiago Peixoto's avatar
Tiago Peixoto committed
497
    }
498

499
    void update_edge(size_t, bool) {}
Tiago Peixoto's avatar
Tiago Peixoto committed
500

501
private:
Tiago Peixoto's avatar
Tiago Peixoto committed
502
503
504
505
    typedef pair<size_t, size_t> deg_t;
    typedef tr1::unordered_map<deg_t,
                               vector<pair<size_t, bool> >,
                               hash<deg_t> >
506
        edges_by_end_deg_t;
507
    edges_by_end_deg_t _edges_by_target;
508
509
510

protected:
    const Graph& _g;
511
512
};

513
514
515

// general stochastic blockmodel
template <class Graph, class EdgeIndexMap, class CorrProb, class BlockDeg>
516
517
518
class ProbabilisticRewireStrategy:
    public RewireStrategyBase<Graph, EdgeIndexMap,
                              ProbabilisticRewireStrategy<Graph, EdgeIndexMap,
519
                                                          CorrProb, BlockDeg> >
520
521
522
{
public:
    typedef RewireStrategyBase<Graph, EdgeIndexMap,
523
                               ProbabilisticRewireStrategy<Graph, EdgeIndexMap,
524
                                                           CorrProb, BlockDeg> >
525
526
527
528
529
        base_t;

    typedef Graph graph_t;
    typedef EdgeIndexMap edge_index_t;

530
    typedef typename BlockDeg::block_t deg_t;
531
532
533
534
535
536

    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;

    ProbabilisticRewireStrategy(Graph& g, EdgeIndexMap edge_index,
537
                                vector<edge_t>& edges, CorrProb corr_prob,
538
539
540
                                BlockDeg blockdeg, bool cache, rng_t& rng)
        : base_t(g, edge_index, edges, rng), _g(g), _corr_prob(corr_prob),
          _blockdeg(blockdeg)
541
    {
542
        tr1::unordered_set<deg_t, hash<deg_t> > deg_set;
543
        for (size_t ei = 0; ei < base_t::_edges.size(); ++ei)
544
        {
545
546
547
            edge_t& e = base_t::_edges[ei];
            deg_set.insert(get_deg(source(e, g), g));
            deg_set.insert(get_deg(target(e, g), g));
548
549
        }

550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
        if (cache)
        {
            // cache probabilities
            for (typeof(deg_set.begin()) s_iter = deg_set.begin();
                 s_iter != deg_set.end(); ++s_iter)
                for (typeof(deg_set.begin()) t_iter = deg_set.begin();
                     t_iter != deg_set.end(); ++t_iter)
                {
                    double p = _corr_prob(*s_iter, *t_iter);
                    if (isnan(p) || isinf(p) || p < 0)
                        p = 0;
                    _probs[make_pair(*s_iter, *t_iter)] = p;
                }
        }
    }

    double get_prob(const deg_t& s_deg, const deg_t& t_deg)
    {
        if (_probs.empty())
        {
            double p = _corr_prob(s_deg, t_deg);
            if (isnan(p) || isinf(p) || p < 0)
                p = 0;
            return p;
        }
        return _probs[make_pair(s_deg, t_deg)];
576
577
578
579
    }

    deg_t get_deg(vertex_t v, Graph& g)
    {
580
        return _blockdeg.get_block(v, g);
581
582
    }

583
    pair<size_t, bool> get_target_edge(size_t ei)
584
    {
585
586
        deg_t s_deg = get_deg(source(base_t::_edges[ei], _g), _g);
        deg_t t_deg = get_deg(target(base_t::_edges[ei], _g), _g);
Tiago Peixoto's avatar
Tiago Peixoto committed
587

588
589
590
591
592
593
594
595
596
        tr1::uniform_int<> sample(0, base_t::_edges.size() - 1);
        size_t epi = sample(base_t::_rng);
        pair<size_t, bool> ep = make_pair(epi, false);
        if (!is_directed::apply<Graph>::type::value)
        {
            // for undirected graphs we must select a random direction
            tr1::bernoulli_distribution coin(0.5);
            ep.second = coin(base_t::_rng);
        }
597
598
599
600

        deg_t ep_s_deg = get_deg(source(ep, base_t::_edges, _g), _g);
        deg_t ep_t_deg = get_deg(target(ep, base_t::_edges, _g), _g);

601
602
603
604
        double pi = (get_prob(s_deg, t_deg) *
                     get_prob(ep_s_deg, ep_t_deg));
        double pf = (get_prob(s_deg, ep_t_deg) *
                     get_prob(ep_s_deg, t_deg));
Tiago Peixoto's avatar
Tiago Peixoto committed
605

606
        if (pf >= pi || pi == 0)
607
            return ep;
608

609
610
        if (pf == 0)
            return make_pair(ei, false); // reject
Tiago Peixoto's avatar
Tiago Peixoto committed
611

612
        double a = exp(log(pf) - log(pi));
Tiago Peixoto's avatar
Tiago Peixoto committed
613

614
615
616
617
618
        tr1::variate_generator<rng_t&, tr1::uniform_real<> >
            rsample(base_t::_rng, tr1::uniform_real<>(0.0, 1.0));
        double r = rsample();
        if (r > a)
            return make_pair(ei, false); // reject
619
        else
620
            return ep;
621
622
    }

623
    void update_edge(size_t, bool) {}
624

625
626
627
private:
    Graph& _g;
    EdgeIndexMap _edge_index;
628
    CorrProb _corr_prob;
629
630
631
632
    BlockDeg _blockdeg;
    tr1::unordered_map<pair<deg_t, deg_t>, double, hash<pair<deg_t, deg_t> > >
        _probs;
};
633

634
635
636
} // graph_tool namespace

#endif // GRAPH_REWIRING_HH