graph_rewiring.hh 20.7 KB
Newer Older
1
2
// graph-tool -- a general graph modification and manipulation thingy
//
Tiago Peixoto's avatar
Tiago Peixoto committed
3
// Copyright (C) 2007-2012 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
23
24
25
26
27
#if (GCC_VERSION >= 40400)
#   include <tr1/unordered_set>
#   include <tr1/random>
#else
#   include <boost/tr1/unordered_set.hpp>
#   include <boost/tr1/random.hpp>
#endif
28
#include <boost/functional/hash.hpp>
29

Tiago Peixoto's avatar
Tiago Peixoto committed
30
31
#include <iostream>

32
33
#include "graph.hh"
#include "graph_filtering.hh"
34
#include "graph_util.hh"
35
#include "sampler.hh"
36
37
38
39
40
41

namespace graph_tool
{
using namespace std;
using namespace boost;

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

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

Tiago Peixoto's avatar
Tiago Peixoto committed
66

67
68
// 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
69
struct swap_edge
70
{
71
    template <class Graph>
Tiago Peixoto's avatar
Tiago Peixoto committed
72
    static bool
73
74
75
    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
76
77
78
79
80
81
82
    {
        // 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)
        //
83
        // no parallel edges are introduced.
Tiago Peixoto's avatar
Tiago Peixoto committed
84
85

        typename graph_traits<Graph>::vertex_descriptor
86
87
88
89
            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
90

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

98
99
100
101
102
    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,
         EdgeIndexMap edge_index, Graph& g)
103
    {
104
105
        // swap the source of the edge 'e' with the source of edge 'se', as
        // such:
106
        //
107
108
        //  (s)    -e--> (t)          (s)    -e--> (nt)
        //  (te_s) -te-> (nt)   =>    (te_s) -te-> (t)
109

110
111
        if (e == te.first)
            return;
112

113
114
115
116
117
118
119
120
        // new edges which will replace the old ones
        typename graph_traits<Graph>::edge_descriptor ne, nte;
        ne = add_edge(source(edges[e], g), target(te, edges, g), g).first;
        if (!te.second)
            nte = add_edge(source(te, edges, g), target(edges[e], g), g).first;
        else // keep invertedness (only for undirected graphs)
            nte = add_edge(target(edges[e], g), source(te, edges, g),  g).first;

121
        edge_index[nte] = edge_index[edges[te.first]];
122
123
124
        remove_edge(edges[te.first], g);
        edges[te.first] = nte;

125
        edge_index[ne] = edge_index[edges[e]];
126
127
        remove_edge(edges[e], g);
        edges[e] = ne;
128
    }
129

130
131
};

Tiago Peixoto's avatar
Tiago Peixoto committed
132
// used for verbose display
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
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;
    }
}

152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
//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,
                      const Graph& g) const
    {
        return get(_p, v);
    }

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

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

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

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

213
214
        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
215
        {
216
            edges.push_back(*e);
Tiago Peixoto's avatar
Tiago Peixoto committed
217
218
            edge_pos.push_back(edge_pos.size());
        }
219

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

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

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

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

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

Tiago Peixoto's avatar
Tiago Peixoto committed
268

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

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

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

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

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

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

        return true;
320
321
322
323
324
    }

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


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

341
342
    typedef typename EdgeIndexMap::value_type index_t;

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

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

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

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

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

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

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

376
            swap_edge::swap_target(ei, et, _edges, _edge_index, _g);
377

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

386
        return true;
387
388
    }

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

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

    typedef Graph graph_t;
    typedef EdgeIndexMap edge_index_t;

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

Tiago Peixoto's avatar
Tiago Peixoto committed
417
418
419
    struct hash_index {};
    struct random_index {};

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

425
    pair<size_t,bool> get_target_edge(size_t e)
Tiago Peixoto's avatar
Tiago Peixoto committed
426
    {
427
428
429
        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
430
        {
431
432
            tr1::bernoulli_distribution coin(0.5);
            et.second = coin(base_t::_rng);
Tiago Peixoto's avatar
Tiago Peixoto committed
433
        }
434
        return et;
Tiago Peixoto's avatar
Tiago Peixoto committed
435
436
    }

437
    void update_edge(size_t e, bool insert) {}
438

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

444
445
446

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

    typedef Graph graph_t;

461
462
463
    typedef typename graph_traits<Graph>::vertex_descriptor vertex_t;
    typedef typename graph_traits<Graph>::edge_descriptor edge_t;

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

Tiago Peixoto's avatar
Tiago Peixoto committed
476
477
478
            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));
479
480
481

            if (!is_directed::apply<Graph>::type::value)
            {
Tiago Peixoto's avatar
Tiago Peixoto committed
482
483
484
                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));
485
486
487
            }
        }
    }
488

489
    pair<size_t,bool> get_target_edge(size_t ei)
490
    {
Tiago Peixoto's avatar
Tiago Peixoto committed
491
492
493
494
495
        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];
496
        tr1::uniform_int<> sample(0, elist.size() - 1);
Tiago Peixoto's avatar
Tiago Peixoto committed
497

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

501
    void update_edge(size_t e, bool insert) {}
Tiago Peixoto's avatar
Tiago Peixoto committed
502

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

protected:
    const Graph& _g;
513
514
};

515
516
517

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

    typedef Graph graph_t;
    typedef EdgeIndexMap edge_index_t;

532
    typedef typename BlockDeg::block_t deg_t;
533
534
535
536
537
538

    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,
539
                                vector<edge_t>& edges, CorrProb corr_prob,
540
541
542
                                BlockDeg blockdeg, bool cache, rng_t& rng)
        : base_t(g, edge_index, edges, rng), _g(g), _corr_prob(corr_prob),
          _blockdeg(blockdeg)
543
    {
544
        tr1::unordered_set<deg_t, hash<deg_t> > deg_set;
545
        for (size_t ei = 0; ei < base_t::_edges.size(); ++ei)
546
        {
547
548
549
            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));
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
576
577
        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)];
578
579
580
581
    }

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

585
    pair<size_t, bool> get_target_edge(size_t ei)
586
    {
587
588
        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
589

590
591
592
593
594
595
596
597
598
        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);
        }
599
600
601
602

        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);

603
604
605
606
        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
607

608
        if (pf >= pi || pi == 0)
609
            return ep;
610

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

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

616
617
618
619
620
        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
621
        else
622
            return ep;
623
624
    }

625
626
    void update_edge(size_t e, bool insert) {}

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

636
637
638
} // graph_tool namespace

#endif // GRAPH_REWIRING_HH