graph_rewiring.hh 21.2 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
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56

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

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

template <class Graph>
typename graph_traits<Graph>::vertex_descriptor
71
72
target(const pair<size_t, bool>& e,
       const vector<typename graph_traits<Graph>::edge_descriptor>& edges,
73
74
75
       const Graph& g)
{
    if (e.second)
76
        return source(edges[e.first], g);
77
    else
78
        return target(edges[e.first], g);
79
80
}

Tiago Peixoto's avatar
Tiago Peixoto committed
81

82
83
// 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
84
struct swap_edge
85
{
86
    template <class Graph>
Tiago Peixoto's avatar
Tiago Peixoto committed
87
    static bool
88
89
90
    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
91
92
93
94
95
96
97
    {
        // 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)
        //
98
        // no parallel edges are introduced.
Tiago Peixoto's avatar
Tiago Peixoto committed
99
100

        typename graph_traits<Graph>::vertex_descriptor
101
102
103
104
            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
105

106
        if (is_adjacent(s,  nt, g))
107
            return true; // e would clash with an existing edge
108
        if (is_adjacent(te_s, t, g))
109
            return true; // te would clash with an existing edge
110
111
112
        return false; // the coast is clear - hooray!
    }

113
114
115
116
117
    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)
118
    {
119
120
        // swap the source of the edge 'e' with the source of edge 'se', as
        // such:
121
        //
122
123
        //  (s)    -e--> (t)          (s)    -e--> (nt)
        //  (te_s) -te-> (nt)   =>    (te_s) -te-> (t)
124

125
126
        if (e == te.first)
            return;
127

128
129
130
131
132
133
134
135
        // 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;

136
        edge_index[nte] = edge_index[edges[te.first]];
137
138
139
        remove_edge(edges[te.first], g);
        edges[te.first] = nte;

140
        edge_index[ne] = edge_index[edges[e]];
141
142
        remove_edge(edges[e], g);
        edges[e] = ne;
143
    }
144

145
146
};

Tiago Peixoto's avatar
Tiago Peixoto committed
147
// used for verbose display
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
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;
    }
}

167
class DegreeBlock;
Tiago Peixoto's avatar
Tiago Peixoto committed
168

169
// main rewire loop
170
171
template <template <class Graph, class EdgeIndexMap, class CorrProb,
                    class BlockDeg>
172
          class RewireStrategy>
173
174
struct graph_rewire
{
175
176
177

    template <class Graph, class EdgeIndexMap, class CorrProb,
              class BlockDeg>
178
    void operator()(Graph& g, EdgeIndexMap edge_index, CorrProb corr_prob,
179
                    bool self_loops, bool parallel_edges,
180
181
182
                    pair<size_t, bool> iter_sweep,
                    pair<bool, bool> cache_verbose,
                    size_t& pcount, rng_t& rng, BlockDeg bd)
183
        const
184
185
186
    {
        typedef typename graph_traits<Graph>::vertex_descriptor vertex_t;
        typedef typename graph_traits<Graph>::edge_descriptor edge_t;
187
188
        bool cache = cache_verbose.first;
        bool verbose = cache_verbose.second;
189

190
        vector<edge_t> edges;
Tiago Peixoto's avatar
Tiago Peixoto committed
191
192
        vector<size_t> edge_pos;
        typedef random_permutation_iterator<typename vector<size_t>::iterator,
193
194
195
                                            rng_t>
            random_edge_iter;

196
197
        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
198
        {
199
            edges.push_back(*e);
Tiago Peixoto's avatar
Tiago Peixoto committed
200
201
            edge_pos.push_back(edge_pos.size());
        }
202

203
        RewireStrategy<Graph, EdgeIndexMap, CorrProb, BlockDeg>
204
            rewire(g, edge_index, edges, corr_prob, bd, cache, rng);
205

206
207
208
209
        size_t niter;
        bool no_sweep;
        tie(niter, no_sweep) = iter_sweep;
        pcount = 0;
210
211
        if (verbose)
            cout << "rewiring edges: ";
Tiago Peixoto's avatar
Tiago Peixoto committed
212
        stringstream str;
213
        for (size_t i = 0; i < niter; ++i)
214
        {
Tiago Peixoto's avatar
Tiago Peixoto committed
215
            random_edge_iter
Tiago Peixoto's avatar
Tiago Peixoto committed
216
217
                ei_begin(edge_pos.begin(), edge_pos.end(), rng),
                ei_end(edge_pos.end(), edge_pos.end(), rng);
218
219

            // for each edge rewire its source or target
220
            for (random_edge_iter ei = ei_begin; ei != ei_end; ++ei)
221
            {
222
                size_t e_pos = ei - ei_begin;
223
                if (verbose)
224
225
                    print_progress(i, niter, e_pos, no_sweep ? 1 : edges.size(),
                                   str);
Tiago Peixoto's avatar
Tiago Peixoto committed
226
                bool success = rewire(*ei, self_loops, parallel_edges);
227
228
229
230
231
232
                if (!success)
                    ++pcount;

                if (no_sweep)
                    break;
            }
233
        }
234
235
        if (verbose)
            cout << endl;
236
    }
237
238
239
240

    template <class Graph, class EdgeIndexMap, class CorrProb>
    void operator()(Graph& g, EdgeIndexMap edge_index, CorrProb corr_prob,
                    bool self_loops, bool parallel_edges,
241
242
243
                    pair<size_t, bool> iter_sweep,
                    pair<bool, bool> cache_verbose,
                    size_t& pcount, rng_t& rng)
244
245
246
        const
    {
        operator()(g, edge_index, corr_prob, self_loops, parallel_edges,
247
                   iter_sweep, cache_verbose, pcount, rng, DegreeBlock());
248
    }
249
250
};

Tiago Peixoto's avatar
Tiago Peixoto committed
251

252
253
// this will rewire the edges so that the resulting graph will be entirely
// random (i.e. Erdos-Renyi)
254
template <class Graph, class EdgeIndexMap, class CorrProb, class BlockDeg>
255
256
257
258
259
260
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
261
    ErdosRewireStrategy(Graph& g, EdgeIndexMap edge_index,
262
                        vector<edge_t>& edges, CorrProb, BlockDeg,
263
                        bool, rng_t& rng)
Tiago Peixoto's avatar
Tiago Peixoto committed
264
265
        : _g(g), _edge_index(edge_index), _edges(edges),
          _vertices(HardNumVertices()(g)), _rng(rng)
266
267
268
269
270
271
272
    {
        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;
    }

273
    bool operator()(size_t ei, bool self_loops, bool parallel_edges)
274
    {
275
        //try randomly drawn pairs of vertices
276
277
278
279
        typedef random_permutation_iterator
            <typename graph_traits<Graph>::vertex_iterator, rng_t>
            random_vertex_iter;

280
        tr1::uniform_int<size_t> sample(0, _vertices.size() - 1);
281
282
283
284
285
286
        typename graph_traits<Graph>::vertex_descriptor s, t;
        while (true)
        {
            s = sample(_rng);
            t = sample(_rng);

287
288
            // reject self-loops if not allowed
            if(s == t && !self_loops)
289
290
291
                continue;
            break;
        }
292
293
294
295
296

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

297
        edge_t ne = add_edge(s, t, _g).first;
298
299
300
        _edge_index[ne] = _edge_index[_edges[ei]];
        remove_edge(_edges[ei], _g);
        _edges[ei] = ne;
301
302

        return true;
303
304
305
306
307
    }

private:
    Graph& _g;
    EdgeIndexMap _edge_index;
Tiago Peixoto's avatar
Tiago Peixoto committed
308
    vector<edge_t>& _edges;
309
310
311
312
313
    vector<typename graph_traits<Graph>::vertex_descriptor> _vertices;
    rng_t& _rng;
};


314
315
316
// 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.
317
318
template <class Graph, class EdgeIndexMap, class RewireStrategy>
class RewireStrategyBase
319
320
321
{
public:
    typedef typename graph_traits<Graph>::edge_descriptor edge_t;
322
323
    typedef typename graph_traits<Graph>::edge_descriptor vertex_t;

324
325
    typedef typename EdgeIndexMap::value_type index_t;

Tiago Peixoto's avatar
Tiago Peixoto committed
326
327
    RewireStrategyBase(Graph& g, EdgeIndexMap edge_index, vector<edge_t>& edges,
                       rng_t& rng)
328
        : _g(g), _edge_index(edge_index), _edges(edges),  _rng(rng) {}
329

330
    bool operator()(size_t ei, bool self_loops, bool parallel_edges)
331
    {
Tiago Peixoto's avatar
Tiago Peixoto committed
332
        RewireStrategy& self = *static_cast<RewireStrategy*>(this);
333

334
335
        // try randomly drawn pairs of edges and check if they satisfy all the
        // consistency checks
336

337
        // rewire target
338
        pair<size_t, bool> et = self.get_target_edge(ei);
339
340

        if (!self_loops) // reject self-loops if not allowed
341
        {
342
343
            if((source(_edges[ei], _g) == target(et, _edges, _g)) ||
               (target(_edges[ei], _g) == source(et, _edges, _g)))
344
345
                return false;
        }
346

347
        // reject parallel edges if not allowed
348
        if (!parallel_edges && (et.first != ei))
349
        {
350
            if (swap_edge::parallel_check_target(ei, et, _edges, _g))
351
                return false;
352
        }
353

354
        if (ei != et.first)
355
        {
356
            self.update_edge(ei, false);
357
            self.update_edge(et.first, false);
Tiago Peixoto's avatar
Tiago Peixoto committed
358

359
            swap_edge::swap_target(ei, et, _edges, _edge_index, _g);
360

361
            self.update_edge(ei, true);
362
363
            self.update_edge(et.first, true);
        }
364
365
366
367
        else
        {
            return false;
        }
Tiago Peixoto's avatar
Tiago Peixoto committed
368

369
        return true;
370
371
    }

Tiago Peixoto's avatar
Tiago Peixoto committed
372
protected:
373
374
    Graph& _g;
    EdgeIndexMap _edge_index;
Tiago Peixoto's avatar
Tiago Peixoto committed
375
    vector<edge_t>& _edges;
376
    rng_t& _rng;
377
378
};

379
380
// this will rewire the edges so that the combined (in, out) degree distribution
// will be the same, but all the rest is random
381
template <class Graph, class EdgeIndexMap, class CorrProb, class BlockDeg>
382
383
class RandomRewireStrategy:
    public RewireStrategyBase<Graph, EdgeIndexMap,
384
                              RandomRewireStrategy<Graph, EdgeIndexMap,
385
                                                   CorrProb, BlockDeg> >
386
387
{
public:
388
    typedef RewireStrategyBase<Graph, EdgeIndexMap,
389
                               RandomRewireStrategy<Graph, EdgeIndexMap,
390
                                                    CorrProb, BlockDeg> >
391
392
393
394
395
        base_t;

    typedef Graph graph_t;
    typedef EdgeIndexMap edge_index_t;

396
397
    typedef typename graph_traits<Graph>::vertex_descriptor vertex_t;
    typedef typename graph_traits<Graph>::edge_descriptor edge_t;
398
    typedef typename EdgeIndexMap::value_type index_t;
399

Tiago Peixoto's avatar
Tiago Peixoto committed
400
401
402
    struct hash_index {};
    struct random_index {};

403
    RandomRewireStrategy(Graph& g, EdgeIndexMap edge_index,
404
                         vector<edge_t>& edges, CorrProb, BlockDeg,
405
                         bool, rng_t& rng)
406
        : base_t(g, edge_index, edges, rng), _g(g) {}
Tiago Peixoto's avatar
Tiago Peixoto committed
407

408
    pair<size_t,bool> get_target_edge(size_t e)
Tiago Peixoto's avatar
Tiago Peixoto committed
409
    {
410
411
412
        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
413
        {
414
415
            tr1::bernoulli_distribution coin(0.5);
            et.second = coin(base_t::_rng);
Tiago Peixoto's avatar
Tiago Peixoto committed
416
        }
417
        return et;
Tiago Peixoto's avatar
Tiago Peixoto committed
418
419
    }

420
    void update_edge(size_t e, bool insert) {}
421

422
private:
Tiago Peixoto's avatar
Tiago Peixoto committed
423
424
    Graph& _g;
    EdgeIndexMap _edge_index;
425
};
426

427
428
429

// 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
430
template <class Graph, class EdgeIndexMap, class CorrProb, class BlockDeg>
431
432
class CorrelatedRewireStrategy:
    public RewireStrategyBase<Graph, EdgeIndexMap,
433
                              CorrelatedRewireStrategy<Graph, EdgeIndexMap,
434
                                                       CorrProb, BlockDeg> >
435
436
{
public:
437
    typedef RewireStrategyBase<Graph, EdgeIndexMap,
438
                               CorrelatedRewireStrategy<Graph, EdgeIndexMap,
439
                                                        CorrProb, BlockDeg> >
440
441
442
443
        base_t;

    typedef Graph graph_t;

444
445
446
    typedef typename graph_traits<Graph>::vertex_descriptor vertex_t;
    typedef typename graph_traits<Graph>::edge_descriptor edge_t;

447
448
    CorrelatedRewireStrategy(Graph& g, EdgeIndexMap edge_index,
                             vector<edge_t>& edges, CorrProb, BlockDeg,
449
                             bool, rng_t& rng)
Tiago Peixoto's avatar
Tiago Peixoto committed
450
        : base_t(g, edge_index, edges, rng), _g(g)
451
    {
452
        for (size_t ei = 0; ei < base_t::_edges.size(); ++ei)
453
        {
454
            // For undirected graphs, there is no difference between source and
Tiago Peixoto's avatar
Tiago Peixoto committed
455
            // target, and each edge will appear _twice_ in the list below,
456
            // once for each different ordering of source and target.
457
            edge_t& e = base_t::_edges[ei];
458

Tiago Peixoto's avatar
Tiago Peixoto committed
459
460
461
            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));
462
463
464

            if (!is_directed::apply<Graph>::type::value)
            {
Tiago Peixoto's avatar
Tiago Peixoto committed
465
466
467
                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));
468
469
470
            }
        }
    }
471

472
    pair<size_t,bool> get_target_edge(size_t ei)
473
    {
Tiago Peixoto's avatar
Tiago Peixoto committed
474
475
476
477
478
        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];
479
        tr1::uniform_int<> sample(0, elist.size() - 1);
Tiago Peixoto's avatar
Tiago Peixoto committed
480

481
        return elist[sample(base_t::_rng)];
Tiago Peixoto's avatar
Tiago Peixoto committed
482
    }
483

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

486
private:
Tiago Peixoto's avatar
Tiago Peixoto committed
487
488
489
490
    typedef pair<size_t, size_t> deg_t;
    typedef tr1::unordered_map<deg_t,
                               vector<pair<size_t, bool> >,
                               hash<deg_t> >
491
        edges_by_end_deg_t;
492
    edges_by_end_deg_t _edges_by_target;
493
494
495

protected:
    const Graph& _g;
496
497
};

498
499
500

// general stochastic blockmodel
template <class Graph, class EdgeIndexMap, class CorrProb, class BlockDeg>
501
502
503
class ProbabilisticRewireStrategy:
    public RewireStrategyBase<Graph, EdgeIndexMap,
                              ProbabilisticRewireStrategy<Graph, EdgeIndexMap,
504
                                                          CorrProb, BlockDeg> >
505
506
507
{
public:
    typedef RewireStrategyBase<Graph, EdgeIndexMap,
508
                               ProbabilisticRewireStrategy<Graph, EdgeIndexMap,
509
                                                           CorrProb, BlockDeg> >
510
511
512
513
514
        base_t;

    typedef Graph graph_t;
    typedef EdgeIndexMap edge_index_t;

515
    typedef typename BlockDeg::block_t deg_t;
516
517
518
519
520
521

    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,
522
                                vector<edge_t>& edges, CorrProb corr_prob,
523
524
525
                                BlockDeg blockdeg, bool cache, rng_t& rng)
        : base_t(g, edge_index, edges, rng), _g(g), _corr_prob(corr_prob),
          _blockdeg(blockdeg)
526
    {
527
        tr1::unordered_set<deg_t, hash<deg_t> > deg_set;
528
        for (size_t ei = 0; ei < base_t::_edges.size(); ++ei)
529
        {
530
531
532
            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));
533
534
        }

535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
        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)];
561
562
563
564
    }

    deg_t get_deg(vertex_t v, Graph& g)
    {
565
        return _blockdeg.get_block(v, g);
566
567
    }

568
    pair<size_t, bool> get_target_edge(size_t ei)
569
    {
570
571
        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
572

573
574
575
576
577
578
579
580
581
        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);
        }
582
583
584
585

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

586
587
588
589
        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
590

591
        if (pf >= pi || pi == 0)
592
            return ep;
593

594
595
        if (pf == 0)
            return make_pair(ei, false); // reject
Tiago Peixoto's avatar
Tiago Peixoto committed
596

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

599
600
601
602
603
        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
604
        else
605
            return ep;
606
607
    }

608
609
    void update_edge(size_t e, bool insert) {}

610
611
612
private:
    Graph& _g;
    EdgeIndexMap _edge_index;
613
    CorrProb _corr_prob;
614
615
616
617
    BlockDeg _blockdeg;
    tr1::unordered_map<pair<deg_t, deg_t>, double, hash<pair<deg_t, deg_t> > >
        _probs;
};
618

619
620
621
622
623
//select blocks based on in/out degrees
class DegreeBlock
{
public:
    typedef pair<size_t, size_t> block_t;
Tiago Peixoto's avatar
Tiago Peixoto committed
624

625
626
627
628
629
630
631
    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));
    }
};
632

633
634
635
636
637
638
//select blocks based on property map
template <class PropertyMap>
class PropertyBlock
{
public:
    typedef typename property_traits<PropertyMap>::value_type block_t;
Tiago Peixoto's avatar
Tiago Peixoto committed
639

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

642
643
644
645
646
647
648
649
650
    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;
651
652
};

653
654
655
} // graph_tool namespace

#endif // GRAPH_REWIRING_HH