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

21
22
#include "tr1_include.hh"
#include TR1_HEADER(unordered_set)
23
24
#include TR1_HEADER(tuple)

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

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

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

34
35
#include "random.hh"

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

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

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

Tiago Peixoto's avatar
Tiago Peixoto committed
65

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

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

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

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

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

112
113
        // new edges which will replace the old ones
        typename graph_traits<Graph>::edge_descriptor ne, nte;
114
115
116
117
118
119
        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);
120
121
        remove_edge(edges[te.first], g);

122
        ne = add_edge(s_e, t_te, g).first;
123
        edges[e] = ne;
124
125
126
127
128
        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;
129
    }
130

131
132
};

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

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

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

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

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

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

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

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

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

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

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

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

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

Tiago Peixoto's avatar
Tiago Peixoto committed
275

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

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

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

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

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

        return true;
322
323
324
325
326
    }

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


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

343
344
    typedef typename EdgeIndexMap::value_type index_t;

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

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

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

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

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

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

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

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

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

388
        return true;
389
390
    }

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

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

    typedef Graph graph_t;
    typedef EdgeIndexMap edge_index_t;

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

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

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

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

439
    void update_edge(size_t, bool) {}
440

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

446
447
448

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

    typedef Graph graph_t;

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

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

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

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

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

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

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

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

protected:
    const Graph& _g;
515
516
};

517
518

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

    typedef Graph graph_t;
    typedef EdgeIndexMap edge_index_t;

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

    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,
542
                                vector<edge_t>& edges, CorrProb corr_prob,
543
544
545
                                BlockDeg blockdeg, bool cache, rng_t& rng)
        : base_t(g, edge_index, edges, rng), _g(g), _corr_prob(corr_prob),
          _blockdeg(blockdeg)
546
    {
547
548
549
        if (cache)
        {
            // cache probabilities
550
551
552
553
554
555
556
557
            tr1::unordered_set<deg_t, hash<deg_t> > deg_set;
            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));
            }

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

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

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

599
600
601
602
603
604
605
606
607
        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);
        }
608
609
610
611

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

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

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

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

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

625
626
627
628
629
        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
630
        else
631
            return ep;
632
633
    }

634
    void update_edge(size_t, bool) {}
635

636
637
638
private:
    Graph& _g;
    EdgeIndexMap _edge_index;
639
    CorrProb _corr_prob;
640
641
642
643
    BlockDeg _blockdeg;
    tr1::unordered_map<pair<deg_t, deg_t>, double, hash<pair<deg_t, deg_t> > >
        _probs;
};
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
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
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
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
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
834
835
836
837
838
839
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
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
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
929
930
931
932
933
934
935
936
937
938

// 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)
    {
        tr1::unordered_set<deg_t, hash<deg_t> > deg_set;
        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);
                if (isnan(p) || isinf(p) || p < 0)
                    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;
        tr1::bernoulli_distribution coin(_in_edges[nt].size() /
                                         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];

            tr1::uniform_int<> sample(0, ies.size() - 1);
            size_t epi = ies[sample(base_t::_rng)];

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

            tr1::uniform_int<> sample(0, oes.size() - 1);
            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));

        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
        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;
    tr1::unordered_map<deg_t, Sampler<deg_t>* > _sampler;
    tr1::unordered_map<pair<deg_t, deg_t>, double, hash<pair<deg_t, deg_t> > >
        _probs;

    tr1::unordered_map<deg_t, vector<size_t> > _in_edges;
    tr1::unordered_map<deg_t, vector<size_t> > _out_edges;
    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);
                if (isnan(p) || isinf(p) || p < 0)
                    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];
        tr1::uniform_int<size_t> s_sample(0, svs.size() - 1);
        tr1::uniform_int<size_t> t_sample(0, tvs.size() - 1);

        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;

    tr1::unordered_map<deg_t, vector<vertex_t> > _vertices;

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

};

939
940
941
} // graph_tool namespace

#endif // GRAPH_REWIRING_HH