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