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) 2006-2013 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
#include "tr1_include.hh"
#include TR1_HEADER(unordered_set)
23
#include <boost/functional/hash.hpp>
24

Tiago Peixoto's avatar
Tiago Peixoto committed
25
26
#include <iostream>

27
28
#include "graph.hh"
#include "graph_filtering.hh"
29
#include "graph_util.hh"
30
#include "sampler.hh"
31

32
33
#include "random.hh"

34
35
36
37
38
namespace graph_tool
{
using namespace std;
using namespace boost;

39
40
template <class Graph>
typename graph_traits<Graph>::vertex_descriptor
41
42
source(const pair<size_t, bool>& e,
       const vector<typename graph_traits<Graph>::edge_descriptor>& edges,
43
44
45
       const Graph& g)
{
    if (e.second)
46
        return target(edges[e.first], g);
47
    else
48
        return source(edges[e.first], g);
49
50
51
52
}

template <class Graph>
typename graph_traits<Graph>::vertex_descriptor
53
54
target(const pair<size_t, bool>& e,
       const vector<typename graph_traits<Graph>::edge_descriptor>& edges,
55
56
57
       const Graph& g)
{
    if (e.second)
58
        return source(edges[e.first], g);
59
    else
60
        return target(edges[e.first], g);
61
62
}

Tiago Peixoto's avatar
Tiago Peixoto committed
63

64
65
// 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
66
struct swap_edge
67
{
68
    template <class Graph>
Tiago Peixoto's avatar
Tiago Peixoto committed
69
    static bool
70
71
72
    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
73
74
75
76
77
78
79
    {
        // 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)
        //
80
        // no parallel edges are introduced.
Tiago Peixoto's avatar
Tiago Peixoto committed
81
82

        typename graph_traits<Graph>::vertex_descriptor
83
84
85
86
            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
87

88
        if (is_adjacent(s,  nt, g))
89
            return true; // e would clash with an existing edge
90
        if (is_adjacent(te_s, t, g))
91
            return true; // te would clash with an existing edge
92
93
94
        return false; // the coast is clear - hooray!
    }

95
96
97
98
99
    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)
100
    {
101
102
        // swap the source of the edge 'e' with the source of edge 'se', as
        // such:
103
        //
104
105
        //  (s)    -e--> (t)          (s)    -e--> (nt)
        //  (te_s) -te-> (nt)   =>    (te_s) -te-> (t)
106

107
108
        if (e == te.first)
            return;
109

110
111
112
113
114
115
116
117
        // 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;

118
        edge_index[nte] = edge_index[edges[te.first]];
119
120
121
        remove_edge(edges[te.first], g);
        edges[te.first] = nte;

122
        edge_index[ne] = edge_index[edges[e]];
123
124
        remove_edge(edges[e], g);
        edges[e] = ne;
125
    }
126

127
128
};

Tiago Peixoto's avatar
Tiago Peixoto committed
129
// used for verbose display
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
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;
    }
}

149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
//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));
    }
};

//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) {}

    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;
};
Tiago Peixoto's avatar
Tiago Peixoto committed
182

183
// main rewire loop
184
185
template <template <class Graph, class EdgeIndexMap, class CorrProb,
                    class BlockDeg>
186
          class RewireStrategy>
187
188
struct graph_rewire
{
189
190
191

    template <class Graph, class EdgeIndexMap, class CorrProb,
              class BlockDeg>
192
    void operator()(Graph& g, EdgeIndexMap edge_index, CorrProb corr_prob,
193
                    bool self_loops, bool parallel_edges,
194
195
196
                    pair<size_t, bool> iter_sweep,
                    pair<bool, bool> cache_verbose,
                    size_t& pcount, rng_t& rng, BlockDeg bd)
197
        const
198
199
200
    {
        typedef typename graph_traits<Graph>::vertex_descriptor vertex_t;
        typedef typename graph_traits<Graph>::edge_descriptor edge_t;
201
202
        bool cache = cache_verbose.first;
        bool verbose = cache_verbose.second;
203

204
        vector<edge_t> edges;
Tiago Peixoto's avatar
Tiago Peixoto committed
205
206
        vector<size_t> edge_pos;
        typedef random_permutation_iterator<typename vector<size_t>::iterator,
207
208
209
                                            rng_t>
            random_edge_iter;

210
211
        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
212
        {
213
            edges.push_back(*e);
Tiago Peixoto's avatar
Tiago Peixoto committed
214
215
            edge_pos.push_back(edge_pos.size());
        }
216

217
        RewireStrategy<Graph, EdgeIndexMap, CorrProb, BlockDeg>
218
            rewire(g, edge_index, edges, corr_prob, bd, cache, rng);
219

220
221
222
223
        size_t niter;
        bool no_sweep;
        tie(niter, no_sweep) = iter_sweep;
        pcount = 0;
224
225
        if (verbose)
            cout << "rewiring edges: ";
Tiago Peixoto's avatar
Tiago Peixoto committed
226
        stringstream str;
227
        for (size_t i = 0; i < niter; ++i)
228
        {
Tiago Peixoto's avatar
Tiago Peixoto committed
229
            random_edge_iter
Tiago Peixoto's avatar
Tiago Peixoto committed
230
231
                ei_begin(edge_pos.begin(), edge_pos.end(), rng),
                ei_end(edge_pos.end(), edge_pos.end(), rng);
232
233

            // for each edge rewire its source or target
234
            for (random_edge_iter ei = ei_begin; ei != ei_end; ++ei)
235
            {
236
                size_t e_pos = ei - ei_begin;
237
                if (verbose)
238
239
                    print_progress(i, niter, e_pos, no_sweep ? 1 : edges.size(),
                                   str);
Tiago Peixoto's avatar
Tiago Peixoto committed
240
                bool success = rewire(*ei, self_loops, parallel_edges);
241
242
243
244
245
246
                if (!success)
                    ++pcount;

                if (no_sweep)
                    break;
            }
247
        }
248
249
        if (verbose)
            cout << endl;
250
    }
251
252
253
254

    template <class Graph, class EdgeIndexMap, class CorrProb>
    void operator()(Graph& g, EdgeIndexMap edge_index, CorrProb corr_prob,
                    bool self_loops, bool parallel_edges,
255
256
257
                    pair<size_t, bool> iter_sweep,
                    pair<bool, bool> cache_verbose,
                    size_t& pcount, rng_t& rng)
258
259
260
        const
    {
        operator()(g, edge_index, corr_prob, self_loops, parallel_edges,
261
                   iter_sweep, cache_verbose, pcount, rng, DegreeBlock());
262
    }
263
264
};

Tiago Peixoto's avatar
Tiago Peixoto committed
265

266
267
// this will rewire the edges so that the resulting graph will be entirely
// random (i.e. Erdos-Renyi)
268
template <class Graph, class EdgeIndexMap, class CorrProb, class BlockDeg>
269
270
271
272
273
274
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
275
    ErdosRewireStrategy(Graph& g, EdgeIndexMap edge_index,
276
                        vector<edge_t>& edges, CorrProb, BlockDeg,
277
                        bool, rng_t& rng)
Tiago Peixoto's avatar
Tiago Peixoto committed
278
279
        : _g(g), _edge_index(edge_index), _edges(edges),
          _vertices(HardNumVertices()(g)), _rng(rng)
280
281
282
283
284
285
286
    {
        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;
    }

287
    bool operator()(size_t ei, bool self_loops, bool parallel_edges)
288
    {
289
        //try randomly drawn pairs of vertices
290
291
292
293
        typedef random_permutation_iterator
            <typename graph_traits<Graph>::vertex_iterator, rng_t>
            random_vertex_iter;

294
        tr1::uniform_int<size_t> sample(0, _vertices.size() - 1);
295
296
297
298
299
300
        typename graph_traits<Graph>::vertex_descriptor s, t;
        while (true)
        {
            s = sample(_rng);
            t = sample(_rng);

301
302
            // reject self-loops if not allowed
            if(s == t && !self_loops)
303
304
305
                continue;
            break;
        }
306
307
308
309
310

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

311
        edge_t ne = add_edge(s, t, _g).first;
312
313
314
        _edge_index[ne] = _edge_index[_edges[ei]];
        remove_edge(_edges[ei], _g);
        _edges[ei] = ne;
315
316

        return true;
317
318
319
320
321
    }

private:
    Graph& _g;
    EdgeIndexMap _edge_index;
Tiago Peixoto's avatar
Tiago Peixoto committed
322
    vector<edge_t>& _edges;
323
324
325
326
327
    vector<typename graph_traits<Graph>::vertex_descriptor> _vertices;
    rng_t& _rng;
};


328
329
330
// 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.
331
332
template <class Graph, class EdgeIndexMap, class RewireStrategy>
class RewireStrategyBase
333
334
335
{
public:
    typedef typename graph_traits<Graph>::edge_descriptor edge_t;
336
337
    typedef typename graph_traits<Graph>::edge_descriptor vertex_t;

338
339
    typedef typename EdgeIndexMap::value_type index_t;

Tiago Peixoto's avatar
Tiago Peixoto committed
340
341
    RewireStrategyBase(Graph& g, EdgeIndexMap edge_index, vector<edge_t>& edges,
                       rng_t& rng)
342
        : _g(g), _edge_index(edge_index), _edges(edges),  _rng(rng) {}
343

344
    bool operator()(size_t ei, bool self_loops, bool parallel_edges)
345
    {
Tiago Peixoto's avatar
Tiago Peixoto committed
346
        RewireStrategy& self = *static_cast<RewireStrategy*>(this);
347

348
349
        // try randomly drawn pairs of edges and check if they satisfy all the
        // consistency checks
350

351
        // rewire target
352
        pair<size_t, bool> et = self.get_target_edge(ei);
353
354

        if (!self_loops) // reject self-loops if not allowed
355
        {
356
357
            if((source(_edges[ei], _g) == target(et, _edges, _g)) ||
               (target(_edges[ei], _g) == source(et, _edges, _g)))
358
359
                return false;
        }
360

361
        // reject parallel edges if not allowed
362
        if (!parallel_edges && (et.first != ei))
363
        {
364
            if (swap_edge::parallel_check_target(ei, et, _edges, _g))
365
                return false;
366
        }
367

368
        if (ei != et.first)
369
        {
370
            self.update_edge(ei, false);
371
            self.update_edge(et.first, false);
Tiago Peixoto's avatar
Tiago Peixoto committed
372

373
            swap_edge::swap_target(ei, et, _edges, _edge_index, _g);
374

375
            self.update_edge(ei, true);
376
377
            self.update_edge(et.first, true);
        }
378
379
380
381
        else
        {
            return false;
        }
Tiago Peixoto's avatar
Tiago Peixoto committed
382

383
        return true;
384
385
    }

Tiago Peixoto's avatar
Tiago Peixoto committed
386
protected:
387
388
    Graph& _g;
    EdgeIndexMap _edge_index;
Tiago Peixoto's avatar
Tiago Peixoto committed
389
    vector<edge_t>& _edges;
390
    rng_t& _rng;
391
392
};

393
394
// this will rewire the edges so that the combined (in, out) degree distribution
// will be the same, but all the rest is random
395
template <class Graph, class EdgeIndexMap, class CorrProb, class BlockDeg>
396
397
class RandomRewireStrategy:
    public RewireStrategyBase<Graph, EdgeIndexMap,
398
                              RandomRewireStrategy<Graph, EdgeIndexMap,
399
                                                   CorrProb, BlockDeg> >
400
401
{
public:
402
    typedef RewireStrategyBase<Graph, EdgeIndexMap,
403
                               RandomRewireStrategy<Graph, EdgeIndexMap,
404
                                                    CorrProb, BlockDeg> >
405
406
407
408
409
        base_t;

    typedef Graph graph_t;
    typedef EdgeIndexMap edge_index_t;

410
411
    typedef typename graph_traits<Graph>::vertex_descriptor vertex_t;
    typedef typename graph_traits<Graph>::edge_descriptor edge_t;
412
    typedef typename EdgeIndexMap::value_type index_t;
413

Tiago Peixoto's avatar
Tiago Peixoto committed
414
415
416
    struct hash_index {};
    struct random_index {};

417
    RandomRewireStrategy(Graph& g, EdgeIndexMap edge_index,
418
                         vector<edge_t>& edges, CorrProb, BlockDeg,
419
                         bool, rng_t& rng)
420
        : base_t(g, edge_index, edges, rng), _g(g) {}
Tiago Peixoto's avatar
Tiago Peixoto committed
421

422
    pair<size_t,bool> get_target_edge(size_t e)
Tiago Peixoto's avatar
Tiago Peixoto committed
423
    {
424
425
426
        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
427
        {
428
429
            tr1::bernoulli_distribution coin(0.5);
            et.second = coin(base_t::_rng);
Tiago Peixoto's avatar
Tiago Peixoto committed
430
        }
431
        return et;
Tiago Peixoto's avatar
Tiago Peixoto committed
432
433
    }

434
    void update_edge(size_t e, bool insert) {}
435

436
private:
Tiago Peixoto's avatar
Tiago Peixoto committed
437
438
    Graph& _g;
    EdgeIndexMap _edge_index;
439
};
440

441
442
443

// 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
444
template <class Graph, class EdgeIndexMap, class CorrProb, class BlockDeg>
445
446
class CorrelatedRewireStrategy:
    public RewireStrategyBase<Graph, EdgeIndexMap,
447
                              CorrelatedRewireStrategy<Graph, EdgeIndexMap,
448
                                                       CorrProb, BlockDeg> >
449
450
{
public:
451
    typedef RewireStrategyBase<Graph, EdgeIndexMap,
452
                               CorrelatedRewireStrategy<Graph, EdgeIndexMap,
453
                                                        CorrProb, BlockDeg> >
454
455
456
457
        base_t;

    typedef Graph graph_t;

458
459
460
    typedef typename graph_traits<Graph>::vertex_descriptor vertex_t;
    typedef typename graph_traits<Graph>::edge_descriptor edge_t;

461
462
    CorrelatedRewireStrategy(Graph& g, EdgeIndexMap edge_index,
                             vector<edge_t>& edges, CorrProb, BlockDeg,
463
                             bool, rng_t& rng)
Tiago Peixoto's avatar
Tiago Peixoto committed
464
        : base_t(g, edge_index, edges, rng), _g(g)
465
    {
466
        for (size_t ei = 0; ei < base_t::_edges.size(); ++ei)
467
        {
468
            // For undirected graphs, there is no difference between source and
Tiago Peixoto's avatar
Tiago Peixoto committed
469
            // target, and each edge will appear _twice_ in the list below,
470
            // once for each different ordering of source and target.
471
            edge_t& e = base_t::_edges[ei];
472

Tiago Peixoto's avatar
Tiago Peixoto committed
473
474
475
            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));
476
477
478

            if (!is_directed::apply<Graph>::type::value)
            {
Tiago Peixoto's avatar
Tiago Peixoto committed
479
480
481
                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));
482
483
484
            }
        }
    }
485

486
    pair<size_t,bool> get_target_edge(size_t ei)
487
    {
Tiago Peixoto's avatar
Tiago Peixoto committed
488
489
490
491
492
        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];
493
        tr1::uniform_int<> sample(0, elist.size() - 1);
Tiago Peixoto's avatar
Tiago Peixoto committed
494

495
        return elist[sample(base_t::_rng)];
Tiago Peixoto's avatar
Tiago Peixoto committed
496
    }
497

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

500
private:
Tiago Peixoto's avatar
Tiago Peixoto committed
501
502
503
504
    typedef pair<size_t, size_t> deg_t;
    typedef tr1::unordered_map<deg_t,
                               vector<pair<size_t, bool> >,
                               hash<deg_t> >
505
        edges_by_end_deg_t;
506
    edges_by_end_deg_t _edges_by_target;
507
508
509

protected:
    const Graph& _g;
510
511
};

512
513
514

// general stochastic blockmodel
template <class Graph, class EdgeIndexMap, class CorrProb, class BlockDeg>
515
516
517
class ProbabilisticRewireStrategy:
    public RewireStrategyBase<Graph, EdgeIndexMap,
                              ProbabilisticRewireStrategy<Graph, EdgeIndexMap,
518
                                                          CorrProb, BlockDeg> >
519
520
521
{
public:
    typedef RewireStrategyBase<Graph, EdgeIndexMap,
522
                               ProbabilisticRewireStrategy<Graph, EdgeIndexMap,
523
                                                           CorrProb, BlockDeg> >
524
525
526
527
528
        base_t;

    typedef Graph graph_t;
    typedef EdgeIndexMap edge_index_t;

529
    typedef typename BlockDeg::block_t deg_t;
530
531
532
533
534
535

    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,
536
                                vector<edge_t>& edges, CorrProb corr_prob,
537
538
539
                                BlockDeg blockdeg, bool cache, rng_t& rng)
        : base_t(g, edge_index, edges, rng), _g(g), _corr_prob(corr_prob),
          _blockdeg(blockdeg)
540
    {
541
        tr1::unordered_set<deg_t, hash<deg_t> > deg_set;
542
        for (size_t ei = 0; ei < base_t::_edges.size(); ++ei)
543
        {
544
545
546
            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));
547
548
        }

549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
        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)];
575
576
577
578
    }

    deg_t get_deg(vertex_t v, Graph& g)
    {
579
        return _blockdeg.get_block(v, g);
580
581
    }

582
    pair<size_t, bool> get_target_edge(size_t ei)
583
    {
584
585
        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
586

587
588
589
590
591
592
593
594
595
        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);
        }
596
597
598
599

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

600
601
602
603
        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
604

605
        if (pf >= pi || pi == 0)
606
            return ep;
607

608
609
        if (pf == 0)
            return make_pair(ei, false); // reject
Tiago Peixoto's avatar
Tiago Peixoto committed
610

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

613
614
615
616
617
        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
618
        else
619
            return ep;
620
621
    }

622
623
    void update_edge(size_t e, bool insert) {}

624
625
626
private:
    Graph& _g;
    EdgeIndexMap _edge_index;
627
    CorrProb _corr_prob;
628
629
630
631
    BlockDeg _blockdeg;
    tr1::unordered_map<pair<deg_t, deg_t>, double, hash<pair<deg_t, deg_t> > >
        _probs;
};
632

633
634
635
} // graph_tool namespace

#endif // GRAPH_REWIRING_HH