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
            edge_pos.push_back(edge_pos.size());
        }
199

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

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

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

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

    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());
    }
245
246
};

Tiago Peixoto's avatar
Tiago Peixoto committed
247

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

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

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

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

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

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

        return true;
299
300
301
302
303
    }

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


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

320
321
    typedef typename EdgeIndexMap::value_type index_t;

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

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

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

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

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

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

350
        if (ei != et.first)
351
        {
352
            self.update_edge(ei, false);
353
            self.update_edge(et.first, false);
Tiago Peixoto's avatar
Tiago Peixoto committed
354

355
            swap_edge::swap_target(ei, et, _edges, _edge_index, _g);
356

357
            self.update_edge(ei, true);
358
359
            self.update_edge(et.first, true);
        }
360
361
362
363
        else
        {
            return false;
        }
Tiago Peixoto's avatar
Tiago Peixoto committed
364

365
        return true;
366
367
    }

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

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

    typedef Graph graph_t;
    typedef EdgeIndexMap edge_index_t;

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

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

399
    RandomRewireStrategy(Graph& g, EdgeIndexMap edge_index,
400
401
402
                         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
403

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

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

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

423
424
425

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

    typedef Graph graph_t;

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

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

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

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

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

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

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

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

protected:
    const Graph& _g;
492
493
};

494
495
496

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

    typedef Graph graph_t;
    typedef EdgeIndexMap edge_index_t;

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

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

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

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

548
    pair<size_t, bool> get_target_edge(size_t ei)
549
    {
550
551
        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
552

553
554
555
556
557
558
559
560
561
        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);
        }
562
563
564
565

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

566
567
568
569
        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)]);
Tiago Peixoto's avatar
Tiago Peixoto committed
570

571
572
        if (pf >= pi)
            return ep;
573

574
575
        if (pf == 0)
            return make_pair(ei, false); // reject
Tiago Peixoto's avatar
Tiago Peixoto committed
576

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

579
580
581
582
583
        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
584
        else
585
            return ep;
586
587
    }

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

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

598
599
600
601
602
//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
603

604
605
606
607
608
609
610
    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));
    }
};
611

612
613
614
615
616
617
//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
618

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

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

632
633
634
} // graph_tool namespace

#endif // GRAPH_REWIRING_HH