graph_rewiring.hh 30.5 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

Tiago Peixoto's avatar
Tiago Peixoto committed
21
22
#include <unordered_set>
#include <tuple>
23

24
#include <boost/functional/hash.hpp>
25

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

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

33
34
#include "random.hh"

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

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

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

Tiago Peixoto's avatar
Tiago Peixoto committed
64

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

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

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

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,
100
         EdgeIndexMap, Graph& g)
101
    {
102
        // swap the target of the edge 'e' with the target of edge 'te', as
103
        // such:
104
        //
105
106
        //  (s)    -e--> (t)          (s)    -e--> (nt)
        //  (te_s) -te-> (nt)   =>    (te_s) -te-> (t)
107

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

111
112
        // new edges which will replace the old ones
        typename graph_traits<Graph>::edge_descriptor ne, nte;
113
114
115
116
117
118
        typename graph_traits<Graph>::vertex_descriptor
            s_e = source(edges[e], g),
            t_e = target(edges[e], g),
            s_te = source(te, edges, g),
            t_te = target(te, edges, g);
        remove_edge(edges[e], g);
119
120
        remove_edge(edges[te.first], g);

121
        ne = add_edge(s_e, t_te, g).first;
122
        edges[e] = ne;
123
124
125
126
127
        if (!te.second)
            nte = add_edge(s_te, t_e, g).first;
        else // keep invertedness (only for undirected graphs)
            nte = add_edge(t_e, s_te,  g).first;
        edges[te.first] = nte;
128
    }
129

130
131
};

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

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
//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,
177
                      const Graph&) const
178
179
180
181
182
183
184
    {
        return get(_p, v);
    }

private:
    PropertyMap _p;
};
Tiago Peixoto's avatar
Tiago Peixoto committed
185

186
// main rewire loop
187
188
template <template <class Graph, class EdgeIndexMap, class CorrProb,
                    class BlockDeg>
189
          class RewireStrategy>
190
191
struct graph_rewire
{
192
193
194

    template <class Graph, class EdgeIndexMap, class CorrProb,
              class BlockDeg>
195
    void operator()(Graph& g, EdgeIndexMap edge_index, CorrProb corr_prob,
196
                    bool self_loops, bool parallel_edges,
197
                    pair<size_t, bool> iter_sweep,
Tiago Peixoto's avatar
Tiago Peixoto committed
198
                    std::tuple<bool, bool, bool> cache_verbose,
199
                    size_t& pcount, rng_t& rng, BlockDeg bd)
200
        const
201
202
    {
        typedef typename graph_traits<Graph>::edge_descriptor edge_t;
Tiago Peixoto's avatar
Tiago Peixoto committed
203
204
205
        bool persist = std::get<0>(cache_verbose);
        bool cache = std::get<1>(cache_verbose);
        bool verbose = std::get<2>(cache_verbose);
206

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

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

220
        RewireStrategy<Graph, EdgeIndexMap, CorrProb, BlockDeg>
221
            rewire(g, edge_index, edges, corr_prob, bd, cache, rng);
222

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

            // for each edge rewire its source or target
237
            for (random_edge_iter ei = ei_begin; ei != ei_end; ++ei)
238
            {
239
                size_t e_pos = ei - ei_begin;
240
                if (verbose)
241
242
                    print_progress(i, niter, e_pos, no_sweep ? 1 : edges.size(),
                                   str);
243
244
245
246
247
248
249
                bool success = false;
                do
                {
                    success = rewire(*ei, self_loops, parallel_edges);
                }
                while(persist && !success);

250
251
252
253
254
255
                if (!success)
                    ++pcount;

                if (no_sweep)
                    break;
            }
256
        }
257
258
        if (verbose)
            cout << endl;
259
    }
260
261
262
263

    template <class Graph, class EdgeIndexMap, class CorrProb>
    void operator()(Graph& g, EdgeIndexMap edge_index, CorrProb corr_prob,
                    bool self_loops, bool parallel_edges,
264
                    pair<size_t, bool> iter_sweep,
Tiago Peixoto's avatar
Tiago Peixoto committed
265
                    std::tuple<bool, bool, bool> cache_verbose,
266
                    size_t& pcount, rng_t& rng)
267
268
269
        const
    {
        operator()(g, edge_index, corr_prob, self_loops, parallel_edges,
270
                   iter_sweep, cache_verbose, pcount, rng, DegreeBlock());
271
    }
272
273
};

Tiago Peixoto's avatar
Tiago Peixoto committed
274

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

296
    bool operator()(size_t ei, bool self_loops, bool parallel_edges)
297
    {
298
        //try randomly drawn pairs of vertices
Tiago Peixoto's avatar
Tiago Peixoto committed
299
        std::uniform_int_distribution<size_t> sample(0, _vertices.size() - 1);
300
301
302
303
304
305
        typename graph_traits<Graph>::vertex_descriptor s, t;
        while (true)
        {
            s = sample(_rng);
            t = sample(_rng);

306
307
            // reject self-loops if not allowed
            if(s == t && !self_loops)
308
309
310
                continue;
            break;
        }
311
312
313
314
315

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

316
        remove_edge(_edges[ei], _g);
317
        edge_t ne = add_edge(s, t, _g).first;
318
        _edges[ei] = ne;
319
320

        return true;
321
322
323
324
325
    }

private:
    Graph& _g;
    EdgeIndexMap _edge_index;
Tiago Peixoto's avatar
Tiago Peixoto committed
326
    vector<edge_t>& _edges;
327
328
329
330
331
    vector<typename graph_traits<Graph>::vertex_descriptor> _vertices;
    rng_t& _rng;
};


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

342
343
    typedef typename EdgeIndexMap::value_type index_t;

Tiago Peixoto's avatar
Tiago Peixoto committed
344
345
    RewireStrategyBase(Graph& g, EdgeIndexMap edge_index, vector<edge_t>& edges,
                       rng_t& rng)
346
        : _g(g), _edge_index(edge_index), _edges(edges),  _rng(rng) {}
347

348
    bool operator()(size_t ei, bool self_loops, bool parallel_edges)
349
    {
Tiago Peixoto's avatar
Tiago Peixoto committed
350
        RewireStrategy& self = *static_cast<RewireStrategy*>(this);
351

352
353
        // try randomly drawn pairs of edges and check if they satisfy all the
        // consistency checks
354

355
        // rewire target
356
        pair<size_t, bool> et = self.get_target_edge(ei);
357
358

        if (!self_loops) // reject self-loops if not allowed
359
        {
360
361
            if((source(_edges[ei], _g) == target(et, _edges, _g)) ||
               (target(_edges[ei], _g) == source(et, _edges, _g)))
362
363
                return false;
        }
364

365
        // reject parallel edges if not allowed
366
        if (!parallel_edges && (et.first != ei))
367
        {
368
            if (swap_edge::parallel_check_target(ei, et, _edges, _g))
369
                return false;
370
        }
371

372
        if (ei != et.first)
373
        {
374
            self.update_edge(ei, false);
375
            self.update_edge(et.first, false);
Tiago Peixoto's avatar
Tiago Peixoto committed
376

377
            swap_edge::swap_target(ei, et, _edges, _edge_index, _g);
378

379
            self.update_edge(ei, true);
380
381
            self.update_edge(et.first, true);
        }
382
383
384
385
        else
        {
            return false;
        }
Tiago Peixoto's avatar
Tiago Peixoto committed
386

387
        return true;
388
389
    }

Tiago Peixoto's avatar
Tiago Peixoto committed
390
protected:
391
392
    Graph& _g;
    EdgeIndexMap _edge_index;
Tiago Peixoto's avatar
Tiago Peixoto committed
393
    vector<edge_t>& _edges;
394
    rng_t& _rng;
395
396
};

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

    typedef Graph graph_t;
    typedef EdgeIndexMap edge_index_t;

414
415
    typedef typename graph_traits<Graph>::vertex_descriptor vertex_t;
    typedef typename graph_traits<Graph>::edge_descriptor edge_t;
416
    typedef typename EdgeIndexMap::value_type index_t;
417

Tiago Peixoto's avatar
Tiago Peixoto committed
418
419
420
    struct hash_index {};
    struct random_index {};

421
    RandomRewireStrategy(Graph& g, EdgeIndexMap edge_index,
422
                         vector<edge_t>& edges, CorrProb, BlockDeg,
423
                         bool, rng_t& rng)
424
        : base_t(g, edge_index, edges, rng), _g(g) {}
Tiago Peixoto's avatar
Tiago Peixoto committed
425

426
    pair<size_t,bool> get_target_edge(size_t)
Tiago Peixoto's avatar
Tiago Peixoto committed
427
    {
Tiago Peixoto's avatar
Tiago Peixoto committed
428
        std::uniform_int_distribution<> sample(0, base_t::_edges.size() - 1);
429
430
        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
431
        {
Tiago Peixoto's avatar
Tiago Peixoto committed
432
            std::bernoulli_distribution coin(0.5);
433
            et.second = coin(base_t::_rng);
Tiago Peixoto's avatar
Tiago Peixoto committed
434
        }
435
        return et;
Tiago Peixoto's avatar
Tiago Peixoto committed
436
437
    }

438
    void update_edge(size_t, bool) {}
439

440
private:
Tiago Peixoto's avatar
Tiago Peixoto committed
441
442
    Graph& _g;
    EdgeIndexMap _edge_index;
443
};
444

445
446
447

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

    typedef Graph graph_t;

462
463
464
    typedef typename graph_traits<Graph>::vertex_descriptor vertex_t;
    typedef typename graph_traits<Graph>::edge_descriptor edge_t;

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

Tiago Peixoto's avatar
Tiago Peixoto committed
477
478
479
            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));
480
481
482

            if (!is_directed::apply<Graph>::type::value)
            {
Tiago Peixoto's avatar
Tiago Peixoto committed
483
484
485
                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));
486
487
488
            }
        }
    }
489

490
    pair<size_t,bool> get_target_edge(size_t ei)
491
    {
Tiago Peixoto's avatar
Tiago Peixoto committed
492
493
494
495
496
        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];
Tiago Peixoto's avatar
Tiago Peixoto committed
497
        std::uniform_int_distribution<> sample(0, elist.size() - 1);
Tiago Peixoto's avatar
Tiago Peixoto committed
498

499
        return elist[sample(base_t::_rng)];
Tiago Peixoto's avatar
Tiago Peixoto committed
500
    }
501

502
    void update_edge(size_t, bool) {}
Tiago Peixoto's avatar
Tiago Peixoto committed
503

504
private:
Tiago Peixoto's avatar
Tiago Peixoto committed
505
    typedef pair<size_t, size_t> deg_t;
Tiago Peixoto's avatar
Tiago Peixoto committed
506
    typedef std::unordered_map<deg_t,
Tiago Peixoto's avatar
Tiago Peixoto committed
507
                               vector<pair<size_t, bool> >,
Tiago Peixoto's avatar
Tiago Peixoto committed
508
                               boost::hash<deg_t> >
509
        edges_by_end_deg_t;
510
    edges_by_end_deg_t _edges_by_target;
511
512
513

protected:
    const Graph& _g;
514
515
};

516
517

// general stochastic blockmodel
518
// this version is based on rejection sampling
519
template <class Graph, class EdgeIndexMap, class CorrProb, class BlockDeg>
520
521
522
class ProbabilisticRewireStrategy:
    public RewireStrategyBase<Graph, EdgeIndexMap,
                              ProbabilisticRewireStrategy<Graph, EdgeIndexMap,
523
                                                          CorrProb, BlockDeg> >
524
525
526
{
public:
    typedef RewireStrategyBase<Graph, EdgeIndexMap,
527
                               ProbabilisticRewireStrategy<Graph, EdgeIndexMap,
528
                                                           CorrProb, BlockDeg> >
529
530
531
532
533
        base_t;

    typedef Graph graph_t;
    typedef EdgeIndexMap edge_index_t;

534
    typedef typename BlockDeg::block_t deg_t;
535
536
537
538
539
540

    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,
541
                                vector<edge_t>& edges, CorrProb corr_prob,
542
543
544
                                BlockDeg blockdeg, bool cache, rng_t& rng)
        : base_t(g, edge_index, edges, rng), _g(g), _corr_prob(corr_prob),
          _blockdeg(blockdeg)
545
    {
546
547
548
        if (cache)
        {
            // cache probabilities
Tiago Peixoto's avatar
Tiago Peixoto committed
549
            std::unordered_set<deg_t, boost::hash<deg_t> > deg_set;
550
551
552
553
554
555
556
            for (size_t ei = 0; ei < base_t::_edges.size(); ++ei)
            {
                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));
            }

557
558
559
560
561
562
            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);
Tiago Peixoto's avatar
Tiago Peixoto committed
563
                    if (std::isnan(p) || std::isinf(p) || p < 0)
564
                        p = 0;
565
566
567
                    // avoid zero probability to not get stuck in rejection step
                    if (p == 0)
                        p = numeric_limits<double>::min();
568
569
570
571
572
573
574
575
576
577
                    _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);
Tiago Peixoto's avatar
Tiago Peixoto committed
578
            if (std::isnan(p) || std::isinf(p) || p < 0)
579
                p = 0;
580
581
582
            // avoid zero probability to not get stuck in rejection step
            if (p == 0)
                p = numeric_limits<double>::min();
583
584
585
            return p;
        }
        return _probs[make_pair(s_deg, t_deg)];
586
587
588
589
    }

    deg_t get_deg(vertex_t v, Graph& g)
    {
590
        return _blockdeg.get_block(v, g);
591
592
    }

593
    pair<size_t, bool> get_target_edge(size_t ei)
594
    {
595
596
        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
597

Tiago Peixoto's avatar
Tiago Peixoto committed
598
        std::uniform_int_distribution<> sample(0, base_t::_edges.size() - 1);
599
600
601
602
603
        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
Tiago Peixoto's avatar
Tiago Peixoto committed
604
            std::bernoulli_distribution coin(0.5);
605
606
            ep.second = coin(base_t::_rng);
        }
607
608
609
610

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

611
612
613
614
        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
615

616
        if (pf >= pi)
617
            return ep;
618

619
620
        if (pf == 0)
            return make_pair(ei, false); // reject
Tiago Peixoto's avatar
Tiago Peixoto committed
621

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

Tiago Peixoto's avatar
Tiago Peixoto committed
624
625
        std::uniform_real_distribution<> rsample(0.0, 1.0);
        double r = rsample(base_t::_rng);
626
627
        if (r > a)
            return make_pair(ei, false); // reject
628
        else
629
            return ep;
630
631
    }

632
    void update_edge(size_t, bool) {}
633

634
635
636
private:
    Graph& _g;
    EdgeIndexMap _edge_index;
637
    CorrProb _corr_prob;
638
    BlockDeg _blockdeg;
Tiago Peixoto's avatar
Tiago Peixoto committed
639
    std::unordered_map<pair<deg_t, deg_t>, double, boost::hash<pair<deg_t, deg_t> > >
640
641
        _probs;
};
642

643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672

// general "degree-corrected" stochastic blockmodel
// this version is based on the alias method
template <class Graph, class EdgeIndexMap, class CorrProb, class BlockDeg>
class AliasProbabilisticRewireStrategy:
    public RewireStrategyBase<Graph, EdgeIndexMap,
                              AliasProbabilisticRewireStrategy<Graph, EdgeIndexMap,
                                                               CorrProb, BlockDeg> >
{
public:
    typedef RewireStrategyBase<Graph, EdgeIndexMap,
                               AliasProbabilisticRewireStrategy<Graph, EdgeIndexMap,
                                                                CorrProb, BlockDeg> >
        base_t;

    typedef Graph graph_t;
    typedef EdgeIndexMap edge_index_t;

    typedef typename BlockDeg::block_t deg_t;

    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;

    AliasProbabilisticRewireStrategy(Graph& g, EdgeIndexMap edge_index,
                                     vector<edge_t>& edges, CorrProb corr_prob,
                                     BlockDeg blockdeg, bool, rng_t& rng)
        : base_t(g, edge_index, edges, rng), _g(g), _corr_prob(corr_prob),
          _blockdeg(blockdeg)
    {
Tiago Peixoto's avatar
Tiago Peixoto committed
673
        std::unordered_set<deg_t> deg_set;
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
        for (size_t ei = 0; ei < base_t::_edges.size(); ++ei)
        {
            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));
        }

        for (typeof(deg_set.begin()) s_iter = deg_set.begin();
             s_iter != deg_set.end(); ++s_iter)
            _items.push_back(*s_iter);

        for (typeof(deg_set.begin()) s_iter = deg_set.begin();
             s_iter != deg_set.end(); ++s_iter)
        {
            vector<double> probs;
            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);
Tiago Peixoto's avatar
Tiago Peixoto committed
693
                if (std::isnan(p) || std::isinf(p) || p < 0)
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
                    p = 0;
                // avoid zero probability to not get stuck in rejection step
                if (p == 0)
                    p = numeric_limits<double>::min();
                probs.push_back(p);
                _probs[make_pair(*s_iter, *t_iter)] = p;
            }

            _sampler[*s_iter] = new Sampler<deg_t>(_items, probs);
        }

        _in_pos.resize(base_t::_edges.size());
        if (!is_directed::apply<Graph>::type::value)
            _out_pos.resize(base_t::_edges.size());
        for (size_t i = 0; i < base_t::_edges.size(); ++i)
        {
            vertex_t v = target(base_t::_edges[i], g);
            deg_t d = get_deg(v, g);
            _in_edges[d].push_back(i);
            _in_pos[i] = _in_edges[d].size() - 1;

            if (!is_directed::apply<Graph>::type::value)
            {
                vertex_t v = source(base_t::_edges[i], g);
                deg_t d = get_deg(v, g);
                _out_edges[d].push_back(i);
                _out_pos[i] = _out_edges[d].size() - 1;
            }
        }
    }

    ~AliasProbabilisticRewireStrategy()
    {
        for (typeof(_sampler.begin()) iter = _sampler.begin();
             iter != _sampler.end(); ++iter)
            delete iter->second;
    }

    double get_prob(const deg_t& s_deg, const deg_t& t_deg)
    {
        return _probs[make_pair(s_deg, t_deg)];
    }

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

    pair<size_t, bool> get_target_edge(size_t ei)
    {
        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);

        deg_t nt = _sampler[s_deg]->sample(base_t::_rng);

        pair<size_t, bool> ep;
Tiago Peixoto's avatar
Tiago Peixoto committed
750
        std::bernoulli_distribution coin(_in_edges[nt].size() /
751
752
753
754
755
756
                                         double(_in_edges[nt].size() +
                                                _out_edges[nt].size()));
        if (is_directed::apply<Graph>::type::value || coin(base_t::_rng))
        {
            vector<size_t>& ies = _in_edges[nt];

Tiago Peixoto's avatar
Tiago Peixoto committed
757
            std::uniform_int_distribution<> sample(0, ies.size() - 1);
758
759
760
761
762
763
764
765
            size_t epi = ies[sample(base_t::_rng)];

            ep = make_pair(epi, false);
        }
        else
        {
            vector<size_t>& oes = _out_edges[nt];

Tiago Peixoto's avatar
Tiago Peixoto committed
766
            std::uniform_int_distribution<> sample(0, oes.size() - 1);
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
            size_t epi = oes[sample(base_t::_rng)];

            ep = make_pair(epi, true);
        }

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

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

        if (pf >= pi)
            return ep;

        if (pf == 0)
            return make_pair(ei, false); // reject

        double a = exp(log(pf) - log(pi));

Tiago Peixoto's avatar
Tiago Peixoto committed
788
789
        std::uniform_real_distribution<> rsample(0.0, 1.0);
        double r = rsample(base_t::_rng);
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
        if (r > a)
            return make_pair(ei, false); // reject
        else
            return ep;
    }

    void update_edge(size_t ei, bool insert)
    {
        if (insert)
        {
            vertex_t v = target(base_t::_edges[ei], _g);
            deg_t d = get_deg(v, _g);
            _in_edges[d].push_back(ei);
            _in_pos[ei] = _in_edges[d].size() - 1;

            if (!is_directed::apply<Graph>::type::value)
            {
                v = source(base_t::_edges[ei], _g);
                deg_t d = get_deg(v, _g);
                _out_edges[d].push_back(ei);
                _out_pos[ei] = _out_edges[d].size() - 1;
            }
        }
        else
        {
            if (!is_directed::apply<Graph>::type::value)
            {
                vertex_t v = target(base_t::_edges[ei], _g);
                deg_t d = get_deg(v, _g);
                size_t j = _in_pos[ei];
                _in_pos[_in_edges[d].back()] = j;
                _in_edges[d][j] = _in_edges[d].back();
                _in_edges[d].pop_back();
            }
        }
    }

private:
    Graph& _g;
    EdgeIndexMap _edge_index;
    CorrProb _corr_prob;
    BlockDeg _blockdeg;

    vector<deg_t> _items;
Tiago Peixoto's avatar
Tiago Peixoto committed
834
835
    std::unordered_map<deg_t, Sampler<deg_t>* > _sampler;
    std::unordered_map<pair<deg_t, deg_t>, double, boost::hash<pair<deg_t, deg_t> > >
836
837
        _probs;

Tiago Peixoto's avatar
Tiago Peixoto committed
838
839
    std::unordered_map<deg_t, vector<size_t> > _in_edges;
    std::unordered_map<deg_t, vector<size_t> > _out_edges;
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
    vector<size_t> _in_pos;
    vector<size_t> _out_pos;
};


// general "traditional" stochastic blockmodel
// this version is based on the alias method, and does not keep the degrees fixed
template <class Graph, class EdgeIndexMap, class CorrProb, class BlockDeg>
class TradBlockRewireStrategy
{
public:
    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;
    typedef typename BlockDeg::block_t deg_t;

    TradBlockRewireStrategy(Graph& g, EdgeIndexMap edge_index,
                            vector<edge_t>& edges, CorrProb corr_prob,
                            BlockDeg blockdeg, bool, rng_t& rng)

        : _g(g), _edge_index(edge_index), _edges(edges), _corr_prob(corr_prob),
          _blockdeg(blockdeg), _rng(rng)
    {
        typename graph_traits<Graph>::vertex_iterator v, v_end;
        for (tie(v, v_end) = vertices(_g); v != v_end; ++v)
        {
            deg_t d = _blockdeg.get_block(*v, g);
            _vertices[d].push_back(*v);
        }

        vector<double> probs;
        for (typeof(_vertices.begin()) s_iter = _vertices.begin();
             s_iter != _vertices.end(); ++s_iter)
        {
            for (typeof(_vertices.begin()) t_iter = _vertices.begin();
                 t_iter != _vertices.end(); ++t_iter)
            {
                double p = _corr_prob(s_iter->first, t_iter->first);
Tiago Peixoto's avatar
Tiago Peixoto committed
878
                if (std::isnan(p) || std::isinf(p) || p < 0)
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
                    p = 0;

                _items.push_back(make_pair(s_iter->first, t_iter->first));
                probs.push_back(p);
            }
        }
        _sampler = new Sampler<pair<deg_t, deg_t> >(_items, probs);
    }

    ~TradBlockRewireStrategy()
    {
        delete _sampler;
    }

    bool operator()(size_t ei, bool self_loops, bool parallel_edges)
    {
        const pair<deg_t, deg_t>& deg = _sampler->sample(_rng);

        vector<vertex_t>& svs = _vertices[deg.first];
        vector<vertex_t>& tvs = _vertices[deg.second];
Tiago Peixoto's avatar
Tiago Peixoto committed
899
900
        std::uniform_int_distribution<size_t> s_sample(0, svs.size() - 1);
        std::uniform_int_distribution<size_t> t_sample(0, tvs.size() - 1);
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928

        typename graph_traits<Graph>::vertex_descriptor s, t;
        s = svs[s_sample(_rng)];
        t = tvs[t_sample(_rng)];

        // reject self-loops if not allowed
        if(!self_loops && s == t)
            return false;

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

        remove_edge(_edges[ei], _g);
        edge_t ne = add_edge(s, t, _g).first;
        _edges[ei] = ne;

        return true;
    }

private:
    Graph& _g;
    EdgeIndexMap _edge_index;
    vector<edge_t>& _edges;
    CorrProb _corr_prob;
    BlockDeg _blockdeg;
    rng_t& _rng;

Tiago Peixoto's avatar
Tiago Peixoto committed
929
    std::unordered_map<deg_t, vector<vertex_t> > _vertices;
930
931
932
933
934
935

    vector<pair<deg_t, deg_t> > _items;
    Sampler<pair<deg_t, deg_t> >* _sampler;

};

936
937
938
} // graph_tool namespace

#endif // GRAPH_REWIRING_HH