graph_rewiring.hh 39.9 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-2014 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
#ifdef HAVE_SPARSEHASH
36
#include SPARSEHASH_INCLUDE(dense_hash_map)
37
38
#endif

39
40
41
42
43
namespace graph_tool
{
using namespace std;
using namespace boost;

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

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

Tiago Peixoto's avatar
Tiago Peixoto committed
68

69
70
// 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
71
struct swap_edge
72
{
73
    template <class Graph>
Tiago Peixoto's avatar
Tiago Peixoto committed
74
    static bool
75
76
    parallel_check_target (const pair<size_t, bool>& e,
                           const pair<size_t, bool>& te,
77
78
                           vector<typename graph_traits<Graph>::edge_descriptor>& edges,
                           const Graph &g)
Tiago Peixoto's avatar
Tiago Peixoto committed
79
80
81
82
83
84
85
    {
        // 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)
        //
86
        // no parallel edges are introduced.
Tiago Peixoto's avatar
Tiago Peixoto committed
87
88

        typename graph_traits<Graph>::vertex_descriptor
89
90
            s = source(e, edges, g),          // current source
            t = target(e, edges, g),          // current target
91
92
            nt = target(te, edges, g),        // new target
            te_s = source(te, edges, g);      // target edge source
93

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

101
    template <class Graph>
102
    static void swap_target
103
104
        (const pair<size_t, bool>& e,
         const pair<size_t, bool>& te,
105
         vector<typename graph_traits<Graph>::edge_descriptor>& edges,
106
         Graph& g)
107
    {
108
        // swap the target of the edge 'e' with the target of edge 'te', as
109
        // such:
110
        //
111
112
        //  (s)    -e--> (t)          (s)    -e--> (nt)
        //  (te_s) -te-> (nt)   =>    (te_s) -te-> (t)
113

114
        if (e.first == te.first)
115
            return;
116

117
118
        // new edges which will replace the old ones
        typename graph_traits<Graph>::edge_descriptor ne, nte;
119
        typename graph_traits<Graph>::vertex_descriptor
120
121
            s_e = source(e, edges, g),
            t_e = target(e, edges, g),
122
123
            s_te = source(te, edges, g),
            t_te = target(te, edges, g);
124
        remove_edge(edges[e.first], g);
125
126
        remove_edge(edges[te.first], g);

127
128
129
130
131
132
        if (is_directed::apply<Graph>::type::value || !e.second)
            ne = add_edge(s_e, t_te, g).first;
        else // keep invertedness (only for undirected graphs)
            ne = add_edge(t_te, s_e, g).first;
        edges[e.first] = ne;
        if (is_directed::apply<Graph>::type::value || !te.second)
133
134
135
136
            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;
137
    }
138

139
140
};

Tiago Peixoto's avatar
Tiago Peixoto committed
141
// used for verbose display
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
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;
    }
}

161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
//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,
186
                      const Graph&) const
187
188
189
190
191
192
193
    {
        return get(_p, v);
    }

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

195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
// select an appropriate "null" key for densehash
template <class Type>
struct get_null_key
{
    Type operator()() const
    {
        return numeric_limits<Type>::max();
    }
};

template <>
struct get_null_key<string>
{
    string operator()() const
    {
        return lexical_cast<string>(get_null_key<size_t>()());
    }
};

template <>
struct get_null_key<boost::python::object>
{
    boost::python::object operator()() const
    {
        return boost::python::object();
    }
};

template <class Type>
struct get_null_key<vector<Type>>
{
    vector<Type> operator()() const
    {
        vector<Type> v(1);
        v[0] = get_null_key<Type>()();
        return v;
    }
};

template <class Type1, class Type2>
struct get_null_key<pair<Type1, Type2>>
{
    pair<Type1, Type2> operator()() const
    {
        return make_pair(get_null_key<Type1>()(),
                         get_null_key<Type2>()());
    }
};


245
// main rewire loop
246
247
template <template <class Graph, class EdgeIndexMap, class CorrProb,
                    class BlockDeg>
248
          class RewireStrategy>
249
250
struct graph_rewire
{
251
252
253

    template <class Graph, class EdgeIndexMap, class CorrProb,
              class BlockDeg>
254
    void operator()(Graph& g, EdgeIndexMap edge_index, CorrProb corr_prob,
255
                    bool self_loops, bool parallel_edges,
256
                    pair<size_t, bool> iter_sweep,
Tiago Peixoto's avatar
Tiago Peixoto committed
257
                    std::tuple<bool, bool, bool> cache_verbose,
258
                    size_t& pcount, rng_t& rng, BlockDeg bd)
259
        const
260
261
    {
        typedef typename graph_traits<Graph>::edge_descriptor edge_t;
Tiago Peixoto's avatar
Tiago Peixoto committed
262
263
264
        bool persist = std::get<0>(cache_verbose);
        bool cache = std::get<1>(cache_verbose);
        bool verbose = std::get<2>(cache_verbose);
265

266
        vector<edge_t> edges;
Tiago Peixoto's avatar
Tiago Peixoto committed
267
        vector<size_t> edge_pos;
268
269
        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
270
        {
271
            edges.push_back(*e);
Tiago Peixoto's avatar
Tiago Peixoto committed
272
273
            edge_pos.push_back(edge_pos.size());
        }
274

275
276
277
278
        typedef random_permutation_iterator<typename vector<size_t>::iterator,
                                            rng_t>
            random_edge_iter;

279
        RewireStrategy<Graph, EdgeIndexMap, CorrProb, BlockDeg>
280
            rewire(g, edge_index, edges, corr_prob, bd, cache, rng);
281

282
283
284
285
        size_t niter;
        bool no_sweep;
        tie(niter, no_sweep) = iter_sweep;
        pcount = 0;
286
287
        if (verbose)
            cout << "rewiring edges: ";
Tiago Peixoto's avatar
Tiago Peixoto committed
288
        stringstream str;
289
        for (size_t i = 0; i < niter; ++i)
290
        {
Tiago Peixoto's avatar
Tiago Peixoto committed
291
            random_edge_iter
Tiago Peixoto's avatar
Tiago Peixoto committed
292
293
                ei_begin(edge_pos.begin(), edge_pos.end(), rng),
                ei_end(edge_pos.end(), edge_pos.end(), rng);
294

295
            for (random_edge_iter ei = ei_begin; ei != ei_end; ++ei)
296
            {
297
                size_t e_pos = ei - ei_begin;
298
                if (verbose)
299
300
                    print_progress(i, niter, e_pos, no_sweep ? 1 : edges.size(),
                                   str);
301
302
303

                size_t e = *ei;

304
305
306
                bool success = false;
                do
                {
307
                    success = rewire(e, self_loops, parallel_edges);
308
309
310
                }
                while(persist && !success);

311
312
313
314
315
316
                if (!success)
                    ++pcount;

                if (no_sweep)
                    break;
            }
317
        }
318
319
        if (verbose)
            cout << endl;
320
    }
321
322
323
324

    template <class Graph, class EdgeIndexMap, class CorrProb>
    void operator()(Graph& g, EdgeIndexMap edge_index, CorrProb corr_prob,
                    bool self_loops, bool parallel_edges,
325
                    pair<size_t, bool> iter_sweep,
Tiago Peixoto's avatar
Tiago Peixoto committed
326
                    std::tuple<bool, bool, bool> cache_verbose,
327
                    size_t& pcount, rng_t& rng)
328
329
330
        const
    {
        operator()(g, edge_index, corr_prob, self_loops, parallel_edges,
331
                   iter_sweep, cache_verbose, pcount, rng, DegreeBlock());
332
    }
333
334
};

Tiago Peixoto's avatar
Tiago Peixoto committed
335

336
337
// this will rewire the edges so that the resulting graph will be entirely
// random (i.e. Erdos-Renyi)
338
template <class Graph, class EdgeIndexMap, class CorrProb, class BlockDeg>
339
340
341
342
343
344
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
345
    ErdosRewireStrategy(Graph& g, EdgeIndexMap edge_index,
346
                        vector<edge_t>& edges, CorrProb, BlockDeg,
347
                        bool, rng_t& rng)
Tiago Peixoto's avatar
Tiago Peixoto committed
348
349
        : _g(g), _edge_index(edge_index), _edges(edges),
          _vertices(HardNumVertices()(g)), _rng(rng)
350
351
352
353
354
355
356
    {
        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;
    }

357
    bool operator()(size_t ei, bool self_loops, bool parallel_edges)
358
    {
359
        //try randomly drawn pairs of vertices
Tiago Peixoto's avatar
Tiago Peixoto committed
360
        std::uniform_int_distribution<size_t> sample(0, _vertices.size() - 1);
361
362
363
364
365
366
        typename graph_traits<Graph>::vertex_descriptor s, t;
        while (true)
        {
            s = sample(_rng);
            t = sample(_rng);

367
368
            // reject self-loops if not allowed
            if(s == t && !self_loops)
369
370
371
                continue;
            break;
        }
372
373
374
375
376

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

377
        remove_edge(_edges[ei], _g);
378
        edge_t ne = add_edge(s, t, _g).first;
379
        _edges[ei] = ne;
380
381

        return true;
382
383
384
385
386
    }

private:
    Graph& _g;
    EdgeIndexMap _edge_index;
Tiago Peixoto's avatar
Tiago Peixoto committed
387
    vector<edge_t>& _edges;
388
389
390
391
392
    vector<typename graph_traits<Graph>::vertex_descriptor> _vertices;
    rng_t& _rng;
};


393
394
395
// 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.
396
397
template <class Graph, class EdgeIndexMap, class RewireStrategy>
class RewireStrategyBase
398
399
400
{
public:
    typedef typename graph_traits<Graph>::edge_descriptor edge_t;
401
    typedef typename graph_traits<Graph>::vertex_descriptor vertex_t;
402

403
404
    typedef typename EdgeIndexMap::value_type index_t;

Tiago Peixoto's avatar
Tiago Peixoto committed
405
406
    RewireStrategyBase(Graph& g, EdgeIndexMap edge_index, vector<edge_t>& edges,
                       rng_t& rng)
407
        : _g(g), _edge_index(edge_index), _edges(edges),  _rng(rng) {}
408

409
    bool operator()(size_t ei, bool self_loops, bool parallel_edges)
410
    {
Tiago Peixoto's avatar
Tiago Peixoto committed
411
        RewireStrategy& self = *static_cast<RewireStrategy*>(this);
412

413
414
        // try randomly drawn pairs of edges and check if they satisfy all the
        // consistency checks
415

416
417
        pair<size_t, bool> e = make_pair(ei, false);

418
        // rewire target
419
        pair<size_t, bool> et = self.get_target_edge(e);
420
421

        if (!self_loops) // reject self-loops if not allowed
422
        {
423
424
            if((source(e, _edges, _g) == target(et, _edges, _g)) ||
               (target(e, _edges, _g) == source(et, _edges, _g)))
425
426
                return false;
        }
427

428
        // reject parallel edges if not allowed
429
        if (!parallel_edges && (et.first != e.first))
430
        {
431
            if (swap_edge::parallel_check_target(e, et, _edges, _g))
432
                return false;
433
        }
434

435
        if (e.first != et.first)
436
        {
437
            self.update_edge(e.first, false);
438
            self.update_edge(et.first, false);
Tiago Peixoto's avatar
Tiago Peixoto committed
439

440
            swap_edge::swap_target(e, et, _edges, _g);
441

442
            self.update_edge(e.first, true);
443
444
            self.update_edge(et.first, true);
        }
445
446
447
448
        else
        {
            return false;
        }
Tiago Peixoto's avatar
Tiago Peixoto committed
449

450
        return true;
451
452
    }

Tiago Peixoto's avatar
Tiago Peixoto committed
453
protected:
454
455
    Graph& _g;
    EdgeIndexMap _edge_index;
Tiago Peixoto's avatar
Tiago Peixoto committed
456
    vector<edge_t>& _edges;
457
    rng_t& _rng;
458
459
};

460
461
// this will rewire the edges so that the combined (in, out) degree distribution
// will be the same, but all the rest is random
462
template <class Graph, class EdgeIndexMap, class CorrProb, class BlockDeg>
463
464
class RandomRewireStrategy:
    public RewireStrategyBase<Graph, EdgeIndexMap,
465
                              RandomRewireStrategy<Graph, EdgeIndexMap,
466
                                                   CorrProb, BlockDeg> >
467
468
{
public:
469
    typedef RewireStrategyBase<Graph, EdgeIndexMap,
470
                               RandomRewireStrategy<Graph, EdgeIndexMap,
471
                                                    CorrProb, BlockDeg> >
472
473
474
475
476
        base_t;

    typedef Graph graph_t;
    typedef EdgeIndexMap edge_index_t;

477
478
    typedef typename graph_traits<Graph>::vertex_descriptor vertex_t;
    typedef typename graph_traits<Graph>::edge_descriptor edge_t;
479
    typedef typename EdgeIndexMap::value_type index_t;
480

Tiago Peixoto's avatar
Tiago Peixoto committed
481
482
483
    struct hash_index {};
    struct random_index {};

484
    RandomRewireStrategy(Graph& g, EdgeIndexMap edge_index,
485
                         vector<edge_t>& edges, CorrProb, BlockDeg,
486
                         bool, rng_t& rng)
487
        : base_t(g, edge_index, edges, rng), _g(g) {}
Tiago Peixoto's avatar
Tiago Peixoto committed
488

489
    pair<size_t,bool> get_target_edge(pair<size_t,bool>& e)
Tiago Peixoto's avatar
Tiago Peixoto committed
490
    {
Tiago Peixoto's avatar
Tiago Peixoto committed
491
        std::uniform_int_distribution<> sample(0, base_t::_edges.size() - 1);
492
493
        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
494
        {
Tiago Peixoto's avatar
Tiago Peixoto committed
495
            std::bernoulli_distribution coin(0.5);
496
            et.second = coin(base_t::_rng);
497
            e.second = coin(base_t::_rng);
Tiago Peixoto's avatar
Tiago Peixoto committed
498
        }
499
        return et;
Tiago Peixoto's avatar
Tiago Peixoto committed
500
501
    }

502
    void update_edge(size_t, bool) {}
503

504
private:
Tiago Peixoto's avatar
Tiago Peixoto committed
505
506
    Graph& _g;
    EdgeIndexMap _edge_index;
507
};
508

509
510
511

// 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
512
template <class Graph, class EdgeIndexMap, class CorrProb, class BlockDeg>
513
514
class CorrelatedRewireStrategy:
    public RewireStrategyBase<Graph, EdgeIndexMap,
515
                              CorrelatedRewireStrategy<Graph, EdgeIndexMap,
516
                                                       CorrProb, BlockDeg> >
517
518
{
public:
519
    typedef RewireStrategyBase<Graph, EdgeIndexMap,
520
                               CorrelatedRewireStrategy<Graph, EdgeIndexMap,
521
                                                        CorrProb, BlockDeg> >
522
523
524
525
        base_t;

    typedef Graph graph_t;

526
527
528
    typedef typename graph_traits<Graph>::vertex_descriptor vertex_t;
    typedef typename graph_traits<Graph>::edge_descriptor edge_t;

529
530
    typedef typename BlockDeg::block_t deg_t;

531
    CorrelatedRewireStrategy(Graph& g, EdgeIndexMap edge_index,
532
                             vector<edge_t>& edges, CorrProb, BlockDeg blockdeg,
533
                             bool, rng_t& rng)
534
        : base_t(g, edge_index, edges, rng), _blockdeg(blockdeg), _g(g)
535
    {
536
537
538
539
#ifdef HAVE_SPARSEHASH
        _edges_by_target.set_empty_key(get_null_key<deg_t>()());
#endif

540
        for (size_t ei = 0; ei < base_t::_edges.size(); ++ei)
541
        {
542
            // For undirected graphs, there is no difference between source and
Tiago Peixoto's avatar
Tiago Peixoto committed
543
            // target, and each edge will appear _twice_ in the list below,
544
            // once for each different ordering of source and target.
545
            edge_t& e = base_t::_edges[ei];
546

Tiago Peixoto's avatar
Tiago Peixoto committed
547
            vertex_t t = target(e, _g);
548
            deg_t tdeg = get_deg(t, _g);;
Tiago Peixoto's avatar
Tiago Peixoto committed
549
            _edges_by_target[tdeg].push_back(make_pair(ei, false));
550
551
552

            if (!is_directed::apply<Graph>::type::value)
            {
Tiago Peixoto's avatar
Tiago Peixoto committed
553
                t = source(e, _g);
554
                deg_t tdeg = get_deg(t, _g);
Tiago Peixoto's avatar
Tiago Peixoto committed
555
                _edges_by_target[tdeg].push_back(make_pair(ei, true));
556
557
558
            }
        }
    }
559

560
    pair<size_t,bool> get_target_edge(pair<size_t, bool>& e)
561
    {
562
563
564
565
566
        if (!is_directed::apply<Graph>::type::value)
        {
            std::bernoulli_distribution coin(0.5);
            e.second = coin(base_t::_rng);
        }
Tiago Peixoto's avatar
Tiago Peixoto committed
567

568
569
570
571
572
573
574
575
        vertex_t t = target(e, base_t::_edges, _g);
        deg_t tdeg = get_deg(t, _g);
        auto& elist = _edges_by_target[tdeg];
        std::uniform_int_distribution<> sample(0, elist.size() - 1);
        auto ep = elist[sample(base_t::_rng)];
        if (get_deg(target(ep, base_t::_edges, _g), _g) != tdeg)
            ep.second = not ep.second;
        return ep;
Tiago Peixoto's avatar
Tiago Peixoto committed
576
    }
577

578
    void update_edge(size_t, bool) {}
Tiago Peixoto's avatar
Tiago Peixoto committed
579

580
581
582
583
584
    deg_t get_deg(vertex_t v, const Graph& g)
    {
        return _blockdeg.get_block(v, g);
    }

585
private:
586
587
588
589
590
591
592
593
    BlockDeg _blockdeg;

#ifdef HAVE_SPARSEHASH
    typedef google::dense_hash_map<deg_t,
                                   vector<pair<size_t, bool>>,
                                   boost::hash<deg_t>>
        edges_by_end_deg_t;
#else
Tiago Peixoto's avatar
Tiago Peixoto committed
594
    typedef std::unordered_map<deg_t,
595
596
                               vector<pair<size_t, bool>>,
                               boost::hash<deg_t>>
597
        edges_by_end_deg_t;
598
599
#endif

600
    edges_by_end_deg_t _edges_by_target;
601
602
603

protected:
    const Graph& _g;
604
605
};

606
607

// general stochastic blockmodel
608
// this version is based on rejection sampling
609
template <class Graph, class EdgeIndexMap, class CorrProb, class BlockDeg>
610
611
612
class ProbabilisticRewireStrategy:
    public RewireStrategyBase<Graph, EdgeIndexMap,
                              ProbabilisticRewireStrategy<Graph, EdgeIndexMap,
613
                                                          CorrProb, BlockDeg> >
614
615
616
{
public:
    typedef RewireStrategyBase<Graph, EdgeIndexMap,
617
                               ProbabilisticRewireStrategy<Graph, EdgeIndexMap,
618
                                                           CorrProb, BlockDeg> >
619
620
621
622
623
        base_t;

    typedef Graph graph_t;
    typedef EdgeIndexMap edge_index_t;

624
    typedef typename BlockDeg::block_t deg_t;
625
626
627
628
629
630

    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,
631
                                vector<edge_t>& edges, CorrProb corr_prob,
632
633
634
                                BlockDeg blockdeg, bool cache, rng_t& rng)
        : base_t(g, edge_index, edges, rng), _g(g), _corr_prob(corr_prob),
          _blockdeg(blockdeg)
635
    {
636
637
638
#ifdef HAVE_SPARSEHASH
        _probs.set_empty_key(get_null_key<pair<deg_t, deg_t>>()());
#endif
639
640
641
        if (cache)
        {
            // cache probabilities
642
            _corr_prob.get_probs(_probs);
643

644
645
646
647
            if (_probs.empty())
            {
                std::unordered_set<deg_t, boost::hash<deg_t> > deg_set;
                for (size_t ei = 0; ei < base_t::_edges.size(); ++ei)
648
                {
649
650
651
                    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));
652
                }
653
654
655
656
657
658
659
660
661
662
663
664
665
666

                for (auto s_iter = deg_set.begin(); s_iter != deg_set.end(); ++s_iter)
                    for (auto t_iter = deg_set.begin(); t_iter != deg_set.end(); ++t_iter)
                    {
                        double p = _corr_prob(*s_iter, *t_iter);
                        _probs[make_pair(*s_iter, *t_iter)] = p;
                    }
            }

            for (auto iter = _probs.begin(); iter != _probs.end(); ++iter)
            {
                double& p = iter->second;
                p = log(p);
            }
667
668
669
670
671
672
673
674
        }
    }

    double get_prob(const deg_t& s_deg, const deg_t& t_deg)
    {
        if (_probs.empty())
        {
            double p = _corr_prob(s_deg, t_deg);
675
            // avoid zero probability to not get stuck in rejection step
676
            if (std::isnan(p) || std::isinf(p) || p <= 0)
677
                p = numeric_limits<double>::min();
678
            return log(p);
679
        }
680
681
682
683
684
        auto k = make_pair(s_deg, t_deg);
        auto iter = _probs.find(k);
        if (iter == _probs.end())
            return log(numeric_limits<double>::min());
        return iter->second;
685
686
687
688
    }

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

692
    pair<size_t, bool> get_target_edge(pair<size_t, bool>& e)
693
    {
694
695
696
697
698
699
700
701
        if (!is_directed::apply<Graph>::type::value)
        {
            std::bernoulli_distribution coin(0.5);
            e.second = coin(base_t::_rng);
        }

        deg_t s_deg = get_deg(source(e, base_t::_edges, _g), _g);
        deg_t t_deg = get_deg(target(e, base_t::_edges, _g), _g);
Tiago Peixoto's avatar
Tiago Peixoto committed
702

Tiago Peixoto's avatar
Tiago Peixoto committed
703
        std::uniform_int_distribution<> sample(0, base_t::_edges.size() - 1);
704
705
706
707
708
        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
709
            std::bernoulli_distribution coin(0.5);
710
711
            ep.second = coin(base_t::_rng);
        }
712
713
714
715

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

716
717
        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
718

719
        if (pf >= pi)
720
            return ep;
721

722
        double a = exp(pf - pi);
Tiago Peixoto's avatar
Tiago Peixoto committed
723

Tiago Peixoto's avatar
Tiago Peixoto committed
724
725
        std::uniform_real_distribution<> rsample(0.0, 1.0);
        double r = rsample(base_t::_rng);
726
        if (r > a)
727
            return e; // reject
728
        else
729
            return ep;
730
731
    }

732
    void update_edge(size_t, bool) {}
733

734
735
736
private:
    Graph& _g;
    EdgeIndexMap _edge_index;
737
    CorrProb _corr_prob;
738
    BlockDeg _blockdeg;
739
740
741
742
743
744
745
746
747
748
749


#ifdef HAVE_SPARSEHASH
    typedef google::dense_hash_map<pair<deg_t, deg_t>, double,
                                   boost::hash<pair<deg_t, deg_t>>> prob_map_t;
#else
    typedef std::unordered_map<pair<deg_t, deg_t>, double,
                               boost::hash<pair<deg_t, deg_t>>> prob_map_t;
#endif

    prob_map_t _probs;
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

// 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)
    {
782
783
784
785
786
787
788
789
790
791
792
793
794

#ifdef HAVE_SPARSEHASH
        _probs.set_empty_key(get_null_key<pair<deg_t, deg_t>>()());
        _sampler.set_empty_key(get_null_key<deg_t>()());
        _in_edges.set_empty_key(get_null_key<deg_t>()());
        _out_edges.set_empty_key(get_null_key<deg_t>()());
        _sprob.set_empty_key(get_null_key<deg_t>()());
#endif

        _in_pos.resize(base_t::_edges.size());
        if (!is_directed::apply<Graph>::type::value)
            _out_pos.resize(base_t::_edges.size());

Tiago Peixoto's avatar
Tiago Peixoto committed
795
        std::unordered_set<deg_t> deg_set;
796
797
798
799
800
        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));
801
802
803
804
805
806
807
808
809
810
811
812
813

            vertex_t v = target(e, g);
            deg_t r = get_deg(v, g);
            _in_edges[r].push_back(ei);
            _in_pos[ei] = _in_edges[r].size() - 1;

            if (!is_directed::apply<Graph>::type::value)
            {
                v = source(e, g);
                deg_t s = get_deg(v, g);
                _out_edges[s].push_back(ei);
                _out_pos[ei] = _out_edges[s].size() - 1;
            }
814
815
        }

816
        _corr_prob.get_probs(_probs);
817

818
        if (_probs.empty())
819
        {
820
821
822
823
824
825
            vector<deg_t> items;

            for (auto s_iter = deg_set.begin(); s_iter != deg_set.end(); ++s_iter)
                items.push_back(*s_iter);

            for (auto s_iter = deg_set.begin(); s_iter != deg_set.end(); ++s_iter)
826
            {
827
828
829
830
831
832
833
834
835
836
837
838
                vector<double> probs;
                double sum = 0;
                for (auto t_iter = deg_set.begin(); t_iter != deg_set.end(); ++t_iter)
                {
                    double p = _corr_prob(*s_iter, *t_iter);
                    // avoid zero probability to not get stuck in rejection step
                    if (std::isnan(p) || std::isinf(p) || p <= 0)
                        continue;
                    probs.push_back(p);
                    _probs[make_pair(*s_iter, *t_iter)] = log(p);
                    sum += p;
                }
839

840
                _sampler[*s_iter] = new Sampler<deg_t, boost::mpl::false_>(items, probs);
841

842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
                auto& ps = _sprob[*s_iter];

#ifdef HAVE_SPARSEHASH
                ps.set_empty_key(get_null_key<deg_t>()());
#endif
                for (size_t i = 0; i < items.size(); ++i)
                {
                    double er = 0;
                    if (!is_directed::apply<Graph>::type::value)
                        er = max(_in_edges[items[i]].size() + _out_edges[items[i]].size(),
                                 numeric_limits<double>::min());
                    else
                        er = max(_in_edges[items[i]].size(),
                                 numeric_limits<double>::min());
                    //ps[items[i]] = exp(log(probs[i]) - log(sum) - log(er));
                    ps[items[i]] = exp(log(probs[i]) - log(sum) - log(er));
                }
            }
        }
        else
862
        {
863
864
865
866
867
868
869
870
871
872
873
874
875
            std::unordered_map<deg_t, vector<double>> sprobs;
            std::unordered_map<deg_t, vector<deg_t>> sitems;
            for (auto iter = _probs.begin(); iter != _probs.end(); ++iter)
            {
                deg_t s = iter->first.first;
                deg_t t = iter->first.second;
                double& p = iter->second;
                if (std::isnan(p) || std::isinf(p) || p <= 0)
                    p = numeric_limits<double>::min();
                sitems[s].push_back(t);
                sprobs[s].push_back(p);
                p = log(p);
            }
876

877
878

            for (auto iter = sitems.begin(); iter != sitems.end(); ++iter)
879
            {
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
                deg_t s = iter->first;
                _sampler[s] = new Sampler<deg_t, boost::mpl::false_>(iter->second, sprobs[s]);

                double sum = 0;
                for (size_t i = 0; i < iter->second.size(); ++i)
                {
                    deg_t t = iter->second[i];
                    sum += sprobs[s][i];
                }

                auto& pr = _sprob[s];

#ifdef HAVE_SPARSEHASH
                pr.set_empty_key(get_null_key<deg_t>()());
#endif

                for (size_t i = 0; i < iter->second.size(); ++i)
                {
                    deg_t t = iter->second[i];
                    pr[t] = exp(log(sprobs[s][i]) - log(sum));
                }
901
902
903
904
905
906
907
908
909
910
911
912
913
            }
        }
    }

    ~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)
    {
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
        static const double zero = log(numeric_limits<double>::min());
        auto k = make_pair(s_deg, t_deg);
        auto iter = _probs.find(k);
        if (iter == _probs.end())
            return zero;
        return iter->second;
    }

    double get_sprob(const deg_t& s_deg, const deg_t& t_deg)
    {
        auto& pr = _sprob[s_deg];
        auto iter = pr.find(t_deg);
        if (iter == pr.end())
            return numeric_limits<double>::min();
        return iter->second;
929
930
931
932
933
934
935
    }

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

936
    pair<size_t, bool> get_target_edge(pair<size_t, bool>& e)
937
    {
938
939
940
941
942
943
944
945
946
947
948
        if (!is_directed::apply<Graph>::type::value)
        {
            std::bernoulli_distribution coin(0.5);
            e.second = coin(base_t::_rng);
        }

        vertex_t s = source(e, base_t::_edges, _g);
        vertex_t t = target(e, base_t::_edges, _g);

        deg_t s_deg = get_deg(s, _g);
        deg_t t_deg = get_deg(t, _g);
949

950
951
952
953
954
        auto iter = _sampler.find(s_deg);
        if (iter == _sampler.end())
            throw GraphException("Block label without defined connection probability!");

        deg_t nt = iter->second->sample(base_t::_rng);
955
956

        pair<size_t, bool> ep;
Tiago Peixoto's avatar
Tiago Peixoto committed
957
        std::bernoulli_distribution coin(_in_edges[nt].size() /
958
959
960
961
962
963
                                         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
964
            std::uniform_int_distribution<> sample(0, ies.size() - 1);
965
966
967
968
969
970
971
972
            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
973
            std::uniform_int_distribution<> sample(0, oes.size() - 1);
974
975
976
977
978
            size_t epi = oes[sample(base_t::_rng)];

            ep = make_pair(epi, true);
        }

979
980
        vertex_t ep_s = source(ep, base_t::_edges, _g);
        vertex_t ep_t = target(ep, base_t::_edges, _g);
981

982
983
        if (t == ep_t)
            return e; // reject
984

985
986
        deg_t ep_s_deg = get_deg(ep_s, _g);
        deg_t ep_t_deg = get_deg(ep_t, _g);
987

988
        assert(ep_t_deg == nt);
989

990
991
992
993
994
995
996
997
998
999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
        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);

        double a = exp(pf - pi);

        if (is_directed::apply<Graph>::type::value)
        {
            double e_ep_t = _in_edges[ep_t_deg].size();
            double e_t = _in_edges[t_deg].size();
            if (t == ep_t)
            {
                e_ep_t -= in_degreeS()(t, _g);
                e_t -= in_degreeS()(ep_t, _g);
            }

            a /= (get_sprob(s_deg, ep_t_deg) / e_ep_t +
                  get_sprob(ep_s_deg, t_deg) / e_t);        // forwards
            a *= (get_sprob(s_deg, t_deg) / e_t +
                  get_sprob(ep_s_deg, ep_t_deg) / e_ep_t);  // backwards
        }
        else
        {
            double e_ep_t = _in_edges[ep_t_deg].size() + _out_edges[ep_t_deg].size();
            double e_t = _in_edges[t_deg].size()  + _out_edges[t_deg].size();
            if (t == ep_t)
            {
                e_ep_t -= out_degree(t, _g);
                e_t -= out_degree(ep_t, _g);
            }

            double e_ep_s = _in_edges[ep_s_deg].size() + _out_edges[ep_s_deg].size();
            double e_s = _in_edges[s_deg].size() + _out_edges[s_deg].size();
            if (s == ep_s)
            {
                e_ep_s -= out_degree(s, _g);
                e_s -= out_degree(ep_s, _g);
            }

            a /= (get_sprob(s_deg, ep_t_deg) / e_ep_t +
                  get_sprob(ep_t_deg, s_deg) / e_s +
                  get_sprob(ep_s_deg, t_deg) / e_t +
                  get_sprob(t_deg, ep_s_deg) / e_ep_s);      // forwards
            a *= (get_sprob(s_deg, t_deg) / e_t +
                  get_sprob(t_deg, s_deg) / e_s +
                  get_sprob(ep_s_deg, ep_t_deg) / e_ep_t +
                  get_sprob(ep_t_deg, ep_s_deg) / e_ep_s);   // backwards
        }

        if (a > 1)
            return ep;
1040

Tiago Peixoto's avatar
Tiago Peixoto committed
1041
1042
        std::uniform_real_distribution<> rsample(0.0, 1.0);
        double r = rsample(base_t::_rng);
1043
        if (r > a)
1044
1045
            return e; // reject
        return ep;
1046
1047
1048
1049
1050
1051
1052
1053
    }

    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);
1054
1055
1056
            auto& in_vec = _in_edges[d];
            in_vec.push_back(ei);
            _in_pos[ei] = in_vec.size() - 1;