graph_rewiring.hh 20.6 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-2011 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
180
                    bool self_loops, bool parallel_edges,
                    pair<size_t, bool> iter_sweep, bool verbose, size_t& pcount,
181
                    rng_t& rng, BlockDeg bd)
182
        const
183
184
185
186
    {
        typedef typename graph_traits<Graph>::vertex_descriptor vertex_t;
        typedef typename graph_traits<Graph>::edge_descriptor edge_t;

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

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

201
202
        RewireStrategy<Graph, EdgeIndexMap, CorrProb, BlockDeg>
            rewire(g, edge_index, edges, corr_prob, bd, rng);
203

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

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

                if (no_sweep)
                    break;
            }
231
        }
232
233
        if (verbose)
            cout << endl;
234
    }
235
236
237
238
239
240
241
242
243
244
245

    template <class Graph, class EdgeIndexMap, class CorrProb>
    void operator()(Graph& g, EdgeIndexMap edge_index, CorrProb corr_prob,
                    bool self_loops, bool parallel_edges,
                    pair<size_t, bool> iter_sweep, bool verbose, size_t& pcount,
                    rng_t& rng)
        const
    {
        operator()(g, edge_index, corr_prob, self_loops, parallel_edges,
                   iter_sweep, verbose, pcount, rng, DegreeBlock());
    }
246
247
};

Tiago Peixoto's avatar
Tiago Peixoto committed
248

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

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

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

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

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

294
        edge_t ne = add_edge(s, t, _g).first;
295
296
297
        _edge_index[ne] = _edge_index[_edges[ei]];
        remove_edge(_edges[ei], _g);
        _edges[ei] = ne;
298
299

        return true;
300
301
302
303
304
    }

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


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

321
322
    typedef typename EdgeIndexMap::value_type index_t;

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

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

331
332
        // try randomly drawn pairs of edges and check if they satisfy all the
        // consistency checks
333

334
        // rewire target
335
        pair<size_t, bool> et = self.get_target_edge(ei);
336
337

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

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

351
        if (ei != et.first)
352
        {
353
            self.update_edge(ei, false);
354
355
            self.update_edge(et.first, false);
            
356
            swap_edge::swap_target(ei, et, _edges, _edge_index, _g);
357

358
            self.update_edge(ei, true);
359
360
            self.update_edge(et.first, true);
        }
361
362
363
364
365
        else
        {
            return false;
        }
        
366
        return true;
367
368
    }

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

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

    typedef Graph graph_t;
    typedef EdgeIndexMap edge_index_t;

393
394
    typedef typename graph_traits<Graph>::vertex_descriptor vertex_t;
    typedef typename graph_traits<Graph>::edge_descriptor edge_t;
395
    typedef typename EdgeIndexMap::value_type index_t;
396

Tiago Peixoto's avatar
Tiago Peixoto committed
397
398
399
    struct hash_index {};
    struct random_index {};

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

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

417
    void update_edge(size_t e, bool insert) {}
418

419
private:
Tiago Peixoto's avatar
Tiago Peixoto committed
420
421
    Graph& _g;
    EdgeIndexMap _edge_index;
422
};
423

424
425
426

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

    typedef Graph graph_t;

441
442
443
    typedef typename graph_traits<Graph>::vertex_descriptor vertex_t;
    typedef typename graph_traits<Graph>::edge_descriptor edge_t;

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

Tiago Peixoto's avatar
Tiago Peixoto committed
456
457
458
            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));
459
460
461

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

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

478
        return elist[sample(base_t::_rng)];
Tiago Peixoto's avatar
Tiago Peixoto committed
479
    }
480

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

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

protected:
    const Graph& _g;
493
494
};

495
496
497

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

    typedef Graph graph_t;
    typedef EdgeIndexMap edge_index_t;

512
    typedef typename BlockDeg::block_t deg_t;
513
514
515
516
517
518

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

531
        // cache probabilities
532
533
534
535
536
        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)
            {
537
538
539
540
                double p = corr_prob(*s_iter, *t_iter);
                if (isnan(p) || isinf(p))
                    p = 0;
                _probs[make_pair(*s_iter, *t_iter)] = p;
541
542
543
544
545
            }
    }

    deg_t get_deg(vertex_t v, Graph& g)
    {
546
        return _blockdeg.get_block(v, g);
547
548
    }

549
    pair<size_t, bool> get_target_edge(size_t ei)
550
    {
551
552
        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);
553
     
554
555
556
557
558
559
560
561
562
        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);
        }
563
564
565
566

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

567
568
569
570
        double pi = (_probs[make_pair(s_deg, t_deg)] *
                     _probs[make_pair(ep_s_deg, ep_t_deg)]);
        double pf = (_probs[make_pair(s_deg, ep_t_deg)] *
                     _probs[make_pair(ep_s_deg, t_deg)]);
571
        
572
573
        if (pf >= pi)
            return ep;
574

575
576
577
578
579
580
581
582
583
584
        if (pf == 0)
            return make_pair(ei, false); // reject
   
        double a = exp(log(pf) - log(pi));
   
        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
585
        else
586
            return ep;
587
588
    }

589
590
    void update_edge(size_t e, bool insert) {}

591
592
593
private:
    Graph& _g;
    EdgeIndexMap _edge_index;
594
595
596
597
    BlockDeg _blockdeg;
    tr1::unordered_map<pair<deg_t, deg_t>, double, hash<pair<deg_t, deg_t> > >
        _probs;
};
598

599
600
601
602
603
604
605
606
607
608
609
610
611
//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));
    }
};
612

613
614
615
616
617
618
619
620
//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) {}
621

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

633
634
635
} // graph_tool namespace

#endif // GRAPH_REWIRING_HH