graph_rewiring.hh 22.3 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
32
33
34
35
#include <iostream>

#include <boost/multi_index_container.hpp>
#include <boost/multi_index/hashed_index.hpp>
#include <boost/multi_index/random_access_index.hpp>
#include <boost/multi_index/identity.hpp>
36
37
38

#include "graph.hh"
#include "graph_filtering.hh"
39
#include "graph_util.hh"
40
#include "sampler.hh"
41
42
43
44
45

namespace graph_tool
{
using namespace std;
using namespace boost;
Tiago Peixoto's avatar
Tiago Peixoto committed
46
using namespace boost::multi_index;
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62

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

63
64
template <class Graph>
typename graph_traits<Graph>::vertex_descriptor
65
66
source(const pair<size_t, bool>& e,
       const vector<typename graph_traits<Graph>::edge_descriptor>& edges,
67
68
69
       const Graph& g)
{
    if (e.second)
70
        return target(edges[e.first], g);
71
    else
72
        return source(edges[e.first], g);
73
74
75
76
}

template <class Graph>
typename graph_traits<Graph>::vertex_descriptor
77
78
target(const pair<size_t, bool>& e,
       const vector<typename graph_traits<Graph>::edge_descriptor>& edges,
79
80
81
       const Graph& g)
{
    if (e.second)
82
        return source(edges[e.first], g);
83
    else
84
        return target(edges[e.first], g);
85
86
}

Tiago Peixoto's avatar
Tiago Peixoto committed
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
// This will iterate over a random permutation of a random access sequence, by
// swapping the values of the sequence as it iterates
template <class RandomAccessIterator, class RNG,
          class RandomDist = tr1::uniform_int<size_t> >
class random_permutation_iterator : public
    std::iterator<input_iterator_tag, typename RandomAccessIterator::value_type>
{
public:
    random_permutation_iterator(RandomAccessIterator begin,
                                RandomAccessIterator end, RNG& rng)
        : _i(begin), _end(end), _rng(&rng)
    {
        if(_i != _end)
        {
            RandomDist random(0,  _end - _i - 1);
            std::iter_swap(_i, _i + random(*_rng));
        }
    }

    typename RandomAccessIterator::value_type operator*()
    {
        return *_i;
    }

    random_permutation_iterator& operator++()
    {
        ++_i;
        if(_i != _end)
        {
            RandomDist random(0,  _end - _i - 1);
            std::iter_swap(_i, _i + random(*_rng));
        }
        return *this;
    }

    bool operator==(const random_permutation_iterator& ri)
    {
        return _i == ri._i;
    }

    bool operator!=(const random_permutation_iterator& ri)
    {
        return _i != ri._i;
    }

    size_t operator-(const random_permutation_iterator& ri)
    {
        return _i - ri._i;
    }

private:
    RandomAccessIterator _i, _end;
    RNG* _rng;
};

142
143
// 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
144
struct swap_edge
145
{
146
    template <class Graph>
Tiago Peixoto's avatar
Tiago Peixoto committed
147
    static bool
148
149
150
    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
151
152
153
154
155
156
157
    {
        // 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)
        //
158
        // no parallel edges are introduced.
Tiago Peixoto's avatar
Tiago Peixoto committed
159
160

        typename graph_traits<Graph>::vertex_descriptor
161
162
163
164
            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
165

166
        if (is_adjacent(s,  nt, g))
167
            return true; // e would clash with an existing edge
168
        if (is_adjacent(te_s, t, g))
169
            return true; // te would clash with an existing edge
170
171
172
        return false; // the coast is clear - hooray!
    }

173
174
175
176
177
    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)
178
    {
179
180
        // swap the source of the edge 'e' with the source of edge 'se', as
        // such:
181
        //
182
183
        //  (s)    -e--> (t)          (s)    -e--> (nt)
        //  (te_s) -te-> (nt)   =>    (te_s) -te-> (t)
184

185
186
        if (e == te.first)
            return;
187

188
189
190
191
192
193
194
195
        // 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;

196
        edge_index[nte] = edge_index[edges[te.first]];
197
198
199
        remove_edge(edges[te.first], g);
        edges[te.first] = nte;

200
        edge_index[ne] = edge_index[edges[e]];
201
202
        remove_edge(edges[e], g);
        edges[e] = ne;
203
    }
204

205
206
};

Tiago Peixoto's avatar
Tiago Peixoto committed
207
// used for verbose display
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
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;
    }
}

227
class DegreeBlock;
Tiago Peixoto's avatar
Tiago Peixoto committed
228

229
// main rewire loop
230
231
template <template <class Graph, class EdgeIndexMap, class CorrProb,
                    class BlockDeg>
232
          class RewireStrategy>
233
234
struct graph_rewire
{
235
236
237

    template <class Graph, class EdgeIndexMap, class CorrProb,
              class BlockDeg>
238
    void operator()(Graph& g, EdgeIndexMap edge_index, CorrProb corr_prob,
239
240
                    bool self_loops, bool parallel_edges,
                    pair<size_t, bool> iter_sweep, bool verbose, size_t& pcount,
241
                    rng_t& rng, BlockDeg bd)
242
        const
243
244
245
246
    {
        typedef typename graph_traits<Graph>::vertex_descriptor vertex_t;
        typedef typename graph_traits<Graph>::edge_descriptor edge_t;

247
248
249
250
251
        vector<edge_t> edges;
        typedef random_permutation_iterator<typename vector<edge_t>::iterator,
                                            rng_t>
            random_edge_iter;

252
253
        typename graph_traits<Graph>::edge_iterator e, e_end;
        for (tie(e, e_end) = boost::edges(g); e != e_end; ++e)
254
            edges.push_back(*e);
255

256
257
        RewireStrategy<Graph, EdgeIndexMap, CorrProb, BlockDeg>
            rewire(g, edge_index, edges, corr_prob, bd, rng);
258

259
260
261
262
        size_t niter;
        bool no_sweep;
        tie(niter, no_sweep) = iter_sweep;
        pcount = 0;
263
264
        if (verbose)
            cout << "rewiring edges: ";
Tiago Peixoto's avatar
Tiago Peixoto committed
265
        stringstream str;
266
        for (size_t i = 0; i < niter; ++i)
267
        {
268
269
270
            random_edge_iter 
                ei_begin(edges.begin(), edges.end(), rng),
                ei_end(edges.end(), edges.end(), rng);
271
272

            // for each edge rewire its source or target
273
            for (random_edge_iter ei = ei_begin; ei != ei_end; ++ei)
274
            {
275
                size_t e_pos = ei - ei_begin;
276
                if (verbose)
277
278
279
                    print_progress(i, niter, e_pos, no_sweep ? 1 : edges.size(),
                                   str);
                bool success = rewire(e_pos, self_loops, parallel_edges);
280
281
282
283
284
285
                if (!success)
                    ++pcount;

                if (no_sweep)
                    break;
            }
286
        }
287
288
        if (verbose)
            cout << endl;
289
    }
290
291
292
293
294
295
296
297
298
299
300

    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());
    }
301
302
};

Tiago Peixoto's avatar
Tiago Peixoto committed
303

304
305
// this will rewire the edges so that the resulting graph will be entirely
// random (i.e. Erdos-Renyi)
306
template <class Graph, class EdgeIndexMap, class CorrProb, class BlockDeg>
307
308
309
310
311
312
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
313
    ErdosRewireStrategy(Graph& g, EdgeIndexMap edge_index,
314
315
                        vector<edge_t>& edges, CorrProb, BlockDeg,
                        rng_t& rng)
Tiago Peixoto's avatar
Tiago Peixoto committed
316
317
        : _g(g), _edge_index(edge_index), _edges(edges),
          _vertices(HardNumVertices()(g)), _rng(rng)
318
319
320
321
322
323
324
    {
        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;
    }

325
    bool operator()(size_t ei, bool self_loops, bool parallel_edges)
326
    {
327
        //try randomly drawn pairs of vertices
328
329
330
331
        typedef random_permutation_iterator
            <typename graph_traits<Graph>::vertex_iterator, rng_t>
            random_vertex_iter;

332
        tr1::uniform_int<size_t> sample(0, _vertices.size() - 1);
333
334
335
336
337
338
        typename graph_traits<Graph>::vertex_descriptor s, t;
        while (true)
        {
            s = sample(_rng);
            t = sample(_rng);

339
340
            // reject self-loops if not allowed
            if(s == t && !self_loops)
341
342
343
                continue;
            break;
        }
344
345
346
347
348

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

349
        edge_t ne = add_edge(s, t, _g).first;
350
351
352
        _edge_index[ne] = _edge_index[_edges[ei]];
        remove_edge(_edges[ei], _g);
        _edges[ei] = ne;
353
354

        return true;
355
356
357
358
359
    }

private:
    Graph& _g;
    EdgeIndexMap _edge_index;
Tiago Peixoto's avatar
Tiago Peixoto committed
360
    vector<edge_t>& _edges;
361
362
363
364
365
    vector<typename graph_traits<Graph>::vertex_descriptor> _vertices;
    rng_t& _rng;
};


366
367
368
// 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.
369
370
template <class Graph, class EdgeIndexMap, class RewireStrategy>
class RewireStrategyBase
371
372
373
{
public:
    typedef typename graph_traits<Graph>::edge_descriptor edge_t;
374
375
    typedef typename graph_traits<Graph>::edge_descriptor vertex_t;

376
377
    typedef typename EdgeIndexMap::value_type index_t;

Tiago Peixoto's avatar
Tiago Peixoto committed
378
379
    RewireStrategyBase(Graph& g, EdgeIndexMap edge_index, vector<edge_t>& edges,
                       rng_t& rng)
380
        : _g(g), _edge_index(edge_index), _edges(edges),  _rng(rng) {}
381

382
    bool operator()(size_t ei, bool self_loops, bool parallel_edges)
383
    {
Tiago Peixoto's avatar
Tiago Peixoto committed
384
        RewireStrategy& self = *static_cast<RewireStrategy*>(this);
385

386
387
        // try randomly drawn pairs of edges and check if they satisfy all the
        // consistency checks
388

389
        // rewire target
390
        pair<size_t, bool> et = self.get_target_edge(ei);
391
392

        if (!self_loops) // reject self-loops if not allowed
393
        {
394
395
            if((source(_edges[ei], _g) == target(et, _edges, _g)) ||
               (target(_edges[ei], _g) == source(et, _edges, _g)))
396
397
                return false;
        }
398

399
        // reject parallel edges if not allowed
400
        if (!parallel_edges && (et.first != ei))
401
        {
402
            if (swap_edge::parallel_check_target(ei, et, _edges, _g))
403
                return false;
404
        }
405

406
        if (ei != et.first)
407
        {
408
            self.update_edge(ei, false);
409
410
            self.update_edge(et.first, false);
            
411
            swap_edge::swap_target(ei, et, _edges, _edge_index, _g);
412

413
            self.update_edge(ei, true);
414
415
            self.update_edge(et.first, true);
        }
416
417
418
419
420
        else
        {
            return false;
        }
        
421
        return true;
422
423
    }

Tiago Peixoto's avatar
Tiago Peixoto committed
424
protected:
425
426
    Graph& _g;
    EdgeIndexMap _edge_index;
Tiago Peixoto's avatar
Tiago Peixoto committed
427
    vector<edge_t>& _edges;
428
    rng_t& _rng;
429
430
};

431
432
// this will rewire the edges so that the combined (in, out) degree distribution
// will be the same, but all the rest is random
433
template <class Graph, class EdgeIndexMap, class CorrProb, class BlockDeg>
434
435
class RandomRewireStrategy:
    public RewireStrategyBase<Graph, EdgeIndexMap,
436
                              RandomRewireStrategy<Graph, EdgeIndexMap,
437
                                                   CorrProb, BlockDeg> >
438
439
{
public:
440
    typedef RewireStrategyBase<Graph, EdgeIndexMap,
441
                               RandomRewireStrategy<Graph, EdgeIndexMap,
442
                                                    CorrProb, BlockDeg> >
443
444
445
446
447
        base_t;

    typedef Graph graph_t;
    typedef EdgeIndexMap edge_index_t;

448
449
    typedef typename graph_traits<Graph>::vertex_descriptor vertex_t;
    typedef typename graph_traits<Graph>::edge_descriptor edge_t;
450
    typedef typename EdgeIndexMap::value_type index_t;
451

Tiago Peixoto's avatar
Tiago Peixoto committed
452
453
454
    struct hash_index {};
    struct random_index {};

455
    RandomRewireStrategy(Graph& g, EdgeIndexMap edge_index,
456
457
458
                         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
459

460
    pair<size_t,bool> get_target_edge(size_t e)
Tiago Peixoto's avatar
Tiago Peixoto committed
461
    {
462
463
464
        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
465
        {
466
467
            tr1::bernoulli_distribution coin(0.5);
            et.second = coin(base_t::_rng);
Tiago Peixoto's avatar
Tiago Peixoto committed
468
        }
469
        return et;
Tiago Peixoto's avatar
Tiago Peixoto committed
470
471
    }

472
    void update_edge(size_t e, bool insert) {}
473

474
private:
Tiago Peixoto's avatar
Tiago Peixoto committed
475
476
    Graph& _g;
    EdgeIndexMap _edge_index;
477
};
478

479
480
481

// 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
482
template <class Graph, class EdgeIndexMap, class CorrProb, class BlockDeg>
483
484
class CorrelatedRewireStrategy:
    public RewireStrategyBase<Graph, EdgeIndexMap,
485
                              CorrelatedRewireStrategy<Graph, EdgeIndexMap,
486
                                                       CorrProb, BlockDeg> >
487
488
{
public:
489
    typedef RewireStrategyBase<Graph, EdgeIndexMap,
490
                               CorrelatedRewireStrategy<Graph, EdgeIndexMap,
491
                                                        CorrProb, BlockDeg> >
492
493
494
495
496
        base_t;

    typedef Graph graph_t;
    typedef EdgeIndexMap edge_index_t;

497
498
499
500
    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;

501
502
503
    CorrelatedRewireStrategy(Graph& g, EdgeIndexMap edge_index,
                             vector<edge_t>& edges, CorrProb, BlockDeg,
                             rng_t& rng)
Tiago Peixoto's avatar
Tiago Peixoto committed
504
        : base_t(g, edge_index, edges, rng), _g(g)
505
    {
506
        for (size_t ei = 0; ei < base_t::_edges.size(); ++ei)
507
        {
508
            // For undirected graphs, there is no difference between source and
Tiago Peixoto's avatar
Tiago Peixoto committed
509
            // target, and each edge will appear _twice_ in the lists below,
510
            // once for each different ordering of source and target.
511
            edge_t& e = base_t::_edges[ei];
512

513
514
515
            _edges_by_target[make_pair(in_degreeS()(target(e, _g), _g),
                                       out_degree(target(e, _g), _g))]
                .push_back(make_pair(ei, false));
516
517
518

            if (!is_directed::apply<Graph>::type::value)
            {
519
520
521
                _edges_by_target[make_pair(in_degreeS()(source(e, _g), _g),
                                           out_degree(source(e, _g), _g))]
                    .push_back(make_pair(ei, true));
522
523
524
            }
        }
    }
525

526
    pair<size_t,bool> get_target_edge(size_t ei)
527
    {
Tiago Peixoto's avatar
Tiago Peixoto committed
528
        pair<size_t, size_t> deg =
529
530
            make_pair(in_degreeS()(target(base_t::_edges[ei], _g), _g),
                      out_degree(target(base_t::_edges[ei], _g), _g));
531
        edges_by_end_deg_t& edges = _edges_by_target;
Tiago Peixoto's avatar
Tiago Peixoto committed
532
        typename edges_by_end_deg_t::mapped_type& elist = edges[deg];
533
        tr1::uniform_int<> sample(0, elist.size() - 1);
534
        return elist[sample(base_t::_rng)];
Tiago Peixoto's avatar
Tiago Peixoto committed
535
    }
536

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

539
private:
540
541
    typedef tr1::unordered_map<pair<size_t, size_t>,
                               vector<pair<index_t, bool> >,
542
543
                               hash<pair<size_t, size_t> > >
        edges_by_end_deg_t;
544
    edges_by_end_deg_t _edges_by_target;
545
546
547

protected:
    const Graph& _g;
548
549
};

550
551
552

// general stochastic blockmodel
template <class Graph, class EdgeIndexMap, class CorrProb, class BlockDeg>
553
554
555
class ProbabilisticRewireStrategy:
    public RewireStrategyBase<Graph, EdgeIndexMap,
                              ProbabilisticRewireStrategy<Graph, EdgeIndexMap,
556
                                                          CorrProb, BlockDeg> >
557
558
559
{
public:
    typedef RewireStrategyBase<Graph, EdgeIndexMap,
560
                               ProbabilisticRewireStrategy<Graph, EdgeIndexMap,
561
                                                           CorrProb, BlockDeg> >
562
563
564
565
566
        base_t;

    typedef Graph graph_t;
    typedef EdgeIndexMap edge_index_t;

567
    typedef typename BlockDeg::block_t deg_t;
568
569
570
571
572
573

    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,
574
575
576
                                vector<edge_t>& edges, CorrProb corr_prob,
                                BlockDeg blockdeg, rng_t& rng)
        : base_t(g, edge_index, edges, rng), _g(g), _blockdeg(blockdeg)
577
    {
578
579
        unordered_set<deg_t, hash<deg_t> > deg_set;
        for (size_t ei = 0; ei < base_t::_edges.size(); ++ei)
580
        {
581
582
583
            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));
584
585
        }

586
        // cache probabilities
587
588
589
590
591
        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)
            {
592
593
594
595
                double p = corr_prob(*s_iter, *t_iter);
                if (isnan(p) || isinf(p))
                    p = 0;
                _probs[make_pair(*s_iter, *t_iter)] = p;
596
597
598
599
600
            }
    }

    deg_t get_deg(vertex_t v, Graph& g)
    {
601
        return _blockdeg.get_block(v, g);
602
603
    }

604
    pair<size_t, bool> get_target_edge(size_t ei)
605
    {
606
607
        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);
608
     
609
610
611
612
613
614
615
616
617
        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);
        }
618
619
620
621

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

622
623
624
625
        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)]);
626
        
627
628
        if (pf >= pi)
            return ep;
629

630
631
632
633
634
635
636
637
638
639
        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
640
        else
641
            return ep;
642
643
    }

644
645
    void update_edge(size_t e, bool insert) {}

646
647
648
private:
    Graph& _g;
    EdgeIndexMap _edge_index;
649
650
651
652
    BlockDeg _blockdeg;
    tr1::unordered_map<pair<deg_t, deg_t>, double, hash<pair<deg_t, deg_t> > >
        _probs;
};
653

654
655
656
657
658
659
660
661
662
663
664
665
666
//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));
    }
};
667

668
669
670
671
672
673
674
675
//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) {}
676

677
678
679
680
681
682
683
684
685
    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;
686
687
};

688
689
690
} // graph_tool namespace

#endif // GRAPH_REWIRING_HH