graph_rewiring.hh 40.6 KB
Newer Older
1
2
// graph-tool -- a general graph modification and manipulation thingy
//
Tiago Peixoto's avatar
Tiago Peixoto committed
3
// Copyright (C) 2006-2017 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
#include <tuple>
Tiago Peixoto's avatar
Tiago Peixoto committed
22
#include <iostream>
23
#include <boost/functional/hash.hpp>
Tiago Peixoto's avatar
Tiago Peixoto committed
24

25
26
#include "graph.hh"
#include "graph_filtering.hh"
27
#include "graph_util.hh"
28
#include "sampler.hh"
29

30
31
#include "random.hh"

32
#include "hash_map_wrap.hh"
33

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

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

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

Tiago Peixoto's avatar
Tiago Peixoto committed
63

64
template <class Nmap, class Graph>
65
void add_count(size_t s, size_t t, Nmap& nvmap, Graph&)
66
67
68
69
70
71
72
73
{
    if (!is_directed::apply<Graph>::type::value && s > t)
        std::swap(s, t);
    auto& nmap = nvmap[s];
    nmap[t]++;
}

template <class Nmap, class Graph>
74
void remove_count(size_t s, size_t t, Nmap& nvmap, Graph&)
75
76
77
78
79
80
81
82
83
84
85
{
    if (!is_directed::apply<Graph>::type::value && s > t)
        std::swap(s, t);
    auto& nmap = nvmap[s];
    auto iter = nmap.find(t);
    iter->second--;
    if (iter->second == 0)
        nmap.erase(iter);
}

template <class Nmap, class Graph>
86
size_t get_count(size_t s, size_t t, Nmap& nvmap, Graph&)
87
88
89
90
91
92
93
94
95
96
97
98
99
100
{
    if (!is_directed::apply<Graph>::type::value && s > t)
        std::swap(s, t);
    auto& nmap = nvmap[s];
    auto iter = nmap.find(t);
    if (iter == nmap.end())
        return 0;
    return iter->second;
    // if (s != t)
    //     return iter->second;
    // else
    //     return iter->second / 2;
}

101
102
// 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
103
struct swap_edge
104
{
105
    template <class Nmap, class Graph>
Tiago Peixoto's avatar
Tiago Peixoto committed
106
    static bool
107
108
    parallel_check_target (const pair<size_t, bool>& e,
                           const pair<size_t, bool>& te,
109
                           vector<typename graph_traits<Graph>::edge_descriptor>& edges,
110
                           Nmap& nmap,
111
                           const Graph &g)
Tiago Peixoto's avatar
Tiago Peixoto committed
112
113
114
115
116
117
118
    {
        // 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)
        //
119
        // no parallel edges are introduced.
Tiago Peixoto's avatar
Tiago Peixoto committed
120
121

        typename graph_traits<Graph>::vertex_descriptor
122
123
            s = source(e, edges, g),          // current source
            t = target(e, edges, g),          // current target
124
125
            nt = target(te, edges, g),        // new target
            te_s = source(te, edges, g);      // target edge source
126

127
        if (get_count(s,  nt, nmap, g) > 0)
128
            return true; // e would clash with an existing edge
129
        if (get_count(te_s, t, nmap, g) > 0)
130
            return true; // te would clash with an existing edge
131
        return false; // no parallel edges
132
133
    }

134
    template <class Graph>
135
    static void swap_target
136
137
        (const pair<size_t, bool>& e,
         const pair<size_t, bool>& te,
138
         vector<typename graph_traits<Graph>::edge_descriptor>& edges,
139
         Graph& g)
140
    {
141
        // swap the target of the edge 'e' with the target of edge 'te', as
142
        // such:
143
        //
144
145
        //  (s)    -e--> (t)          (s)    -e--> (nt)
        //  (te_s) -te-> (nt)   =>    (te_s) -te-> (t)
146

147
        if (e.first == te.first)
148
            return;
149

150
151
        // new edges which will replace the old ones
        typename graph_traits<Graph>::edge_descriptor ne, nte;
152
        typename graph_traits<Graph>::vertex_descriptor
153
154
            s_e = source(e, edges, g),
            t_e = target(e, edges, g),
155
156
            s_te = source(te, edges, g),
            t_te = target(te, edges, g);
157
        remove_edge(edges[e.first], g);
158
159
        remove_edge(edges[te.first], g);

160
161
162
163
164
165
        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)
166
167
168
169
            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;
170
    }
171

172
173
};

Tiago Peixoto's avatar
Tiago Peixoto committed
174
// used for verbose display
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
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;
    }
}

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
//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,
219
                      const Graph&) const
220
221
222
223
224
225
226
    {
        return get(_p, v);
    }

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

228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
// 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>()());
    }
};


278
// main rewire loop
279
280
template <template <class Graph, class EdgeIndexMap, class CorrProb,
                    class BlockDeg>
281
          class RewireStrategy>
282
283
struct graph_rewire
{
284
285

    template <class Graph, class EdgeIndexMap, class CorrProb,
286
              class BlockDeg, class PinMap>
287
    void operator()(Graph& g, EdgeIndexMap edge_index, CorrProb corr_prob,
288
                    PinMap pin, bool self_loops, bool parallel_edges,
289
                    pair<size_t, bool> iter_sweep,
Tiago Peixoto's avatar
Tiago Peixoto committed
290
                    std::tuple<bool, bool, bool> cache_verbose,
291
                    size_t& pcount, rng_t& rng, BlockDeg bd)
292
        const
293
294
    {
        typedef typename graph_traits<Graph>::edge_descriptor edge_t;
Tiago Peixoto's avatar
Tiago Peixoto committed
295
296
297
        bool persist = std::get<0>(cache_verbose);
        bool cache = std::get<1>(cache_verbose);
        bool verbose = std::get<2>(cache_verbose);
298

299
        vector<edge_t> edges;
Tiago Peixoto's avatar
Tiago Peixoto committed
300
        vector<size_t> edge_pos;
301
302
        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
303
        {
304
305
            if (pin[*e])
                continue;
306
            edges.push_back(*e);
Tiago Peixoto's avatar
Tiago Peixoto committed
307
308
            edge_pos.push_back(edge_pos.size());
        }
309

310
311
312
313
        typedef random_permutation_iterator<typename vector<size_t>::iterator,
                                            rng_t>
            random_edge_iter;

314
        RewireStrategy<Graph, EdgeIndexMap, CorrProb, BlockDeg>
315
            rewire(g, edge_index, edges, corr_prob, bd, cache, rng, parallel_edges);
316

317
318
319
320
        size_t niter;
        bool no_sweep;
        tie(niter, no_sweep) = iter_sweep;
        pcount = 0;
321
322
        if (verbose)
            cout << "rewiring edges: ";
Tiago Peixoto's avatar
Tiago Peixoto committed
323
        stringstream str;
324
        for (size_t i = 0; i < niter; ++i)
325
        {
Tiago Peixoto's avatar
Tiago Peixoto committed
326
            random_edge_iter
Tiago Peixoto's avatar
Tiago Peixoto committed
327
328
                ei_begin(edge_pos.begin(), edge_pos.end(), rng),
                ei_end(edge_pos.end(), edge_pos.end(), rng);
329

330
            for (random_edge_iter ei = ei_begin; ei != ei_end; ++ei)
331
            {
332
                size_t e_pos = ei - ei_begin;
333
                if (verbose)
334
335
                    print_progress(i, niter, e_pos, no_sweep ? 1 : edges.size(),
                                   str);
336
337
338

                size_t e = *ei;

339
340
341
                bool success = false;
                do
                {
342
                    success = rewire(e, self_loops, parallel_edges);
343
344
345
                }
                while(persist && !success);

346
347
348
349
350
351
                if (!success)
                    ++pcount;

                if (no_sweep)
                    break;
            }
352
        }
353
354
        if (verbose)
            cout << endl;
355
    }
356

357
    template <class Graph, class EdgeIndexMap, class CorrProb, class PinMap>
358
    void operator()(Graph& g, EdgeIndexMap edge_index, CorrProb corr_prob,
359
                    PinMap pin, bool self_loops, bool parallel_edges,
360
                    pair<size_t, bool> iter_sweep,
Tiago Peixoto's avatar
Tiago Peixoto committed
361
                    std::tuple<bool, bool, bool> cache_verbose,
362
                    size_t& pcount, rng_t& rng)
363
364
        const
    {
365
        operator()(g, edge_index, corr_prob, pin, self_loops, parallel_edges,
366
                   iter_sweep, cache_verbose, pcount, rng, DegreeBlock());
367
    }
368
369
};

Tiago Peixoto's avatar
Tiago Peixoto committed
370

371
372
// this will rewire the edges so that the resulting graph will be entirely
// random (i.e. Erdos-Renyi)
373
template <class Graph, class EdgeIndexMap, class CorrProb, class BlockDeg>
374
375
376
377
378
379
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
380
    ErdosRewireStrategy(Graph& g, EdgeIndexMap edge_index,
381
                        vector<edge_t>& edges, CorrProb, BlockDeg,
382
                        bool, rng_t& rng, bool)
Tiago Peixoto's avatar
Tiago Peixoto committed
383
384
        : _g(g), _edge_index(edge_index), _edges(edges),
          _vertices(HardNumVertices()(g)), _rng(rng)
385
    {
386
        decltype(_vertices.begin()) viter = _vertices.begin();
387
388
389
390
391
        typename graph_traits<Graph>::vertex_iterator v, v_end;
        for (tie(v, v_end) = vertices(_g); v != v_end; ++v)
            *(viter++) = *v;
    }

392
    bool operator()(size_t ei, bool self_loops, bool parallel_edges)
393
    {
394
        //try randomly drawn pairs of vertices
Tiago Peixoto's avatar
Tiago Peixoto committed
395
        std::uniform_int_distribution<size_t> sample(0, _vertices.size() - 1);
396
397
398
399
400
401
        typename graph_traits<Graph>::vertex_descriptor s, t;
        while (true)
        {
            s = sample(_rng);
            t = sample(_rng);

402
403
            // reject self-loops if not allowed
            if(s == t && !self_loops)
404
405
406
                continue;
            break;
        }
407
408
409
410
411

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

412
        remove_edge(_edges[ei], _g);
413
        edge_t ne = add_edge(s, t, _g).first;
414
        _edges[ei] = ne;
415
416

        return true;
417
418
419
420
421
    }

private:
    Graph& _g;
    EdgeIndexMap _edge_index;
Tiago Peixoto's avatar
Tiago Peixoto committed
422
    vector<edge_t>& _edges;
423
424
425
426
427
    vector<typename graph_traits<Graph>::vertex_descriptor> _vertices;
    rng_t& _rng;
};


428

429
430
431
// 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.
432
433
template <class Graph, class EdgeIndexMap, class RewireStrategy>
class RewireStrategyBase
434
435
436
{
public:
    typedef typename graph_traits<Graph>::edge_descriptor edge_t;
437
    typedef typename graph_traits<Graph>::vertex_descriptor vertex_t;
438

439
440
    typedef typename EdgeIndexMap::value_type index_t;

Tiago Peixoto's avatar
Tiago Peixoto committed
441
    RewireStrategyBase(Graph& g, EdgeIndexMap edge_index, vector<edge_t>& edges,
442
443
444
445
446
447
448
449
450
451
                       rng_t& rng, bool parallel_edges)
        : _g(g), _edge_index(edge_index), _edges(edges), _rng(rng),
          _nmap(get(vertex_index, g), num_vertices(g))
    {
        if (!parallel_edges)
        {
            for (size_t i = 0; i < edges.size(); ++i)
                add_count(source(edges[i], g), target(edges[i], g), _nmap, g);
        }
    }
452

453
    bool operator()(size_t ei, bool self_loops, bool parallel_edges)
454
    {
Tiago Peixoto's avatar
Tiago Peixoto committed
455
        RewireStrategy& self = *static_cast<RewireStrategy*>(this);
456

457
458
        // try randomly drawn pairs of edges and check if they satisfy all the
        // consistency checks
459

460
461
        pair<size_t, bool> e = make_pair(ei, false);

462
        // rewire target
463
        pair<size_t, bool> et = self.get_target_edge(e, parallel_edges);
464
465

        if (!self_loops) // reject self-loops if not allowed
466
        {
467
468
            if((source(e, _edges, _g) == target(et, _edges, _g)) ||
               (target(e, _edges, _g) == source(et, _edges, _g)))
469
470
                return false;
        }
471

472
        // reject parallel edges if not allowed
473
        if (!parallel_edges && (et.first != e.first))
474
        {
475
            if (swap_edge::parallel_check_target(e, et, _edges, _nmap, _g))
476
                return false;
477
        }
478

479
        if (e.first != et.first)
480
        {
481
            self.update_edge(e.first, false);
482
            self.update_edge(et.first, false);
Tiago Peixoto's avatar
Tiago Peixoto committed
483

484
485
486
487
488
489
            if (!parallel_edges)
            {
                remove_count(source(e, _edges, _g), target(e, _edges, _g), _nmap, _g);
                remove_count(source(et, _edges, _g), target(et, _edges, _g), _nmap, _g);
            }

490
            swap_edge::swap_target(e, et, _edges, _g);
491

492
            self.update_edge(e.first, true);
493
            self.update_edge(et.first, true);
494
495
496
497
498
499

            if (!parallel_edges)
            {
                add_count(source(e, _edges, _g), target(e, _edges, _g), _nmap, _g);
                add_count(source(et, _edges, _g), target(et, _edges, _g), _nmap, _g);
            }
500
        }
501
502
503
504
        else
        {
            return false;
        }
Tiago Peixoto's avatar
Tiago Peixoto committed
505

506
        return true;
507
508
    }

Tiago Peixoto's avatar
Tiago Peixoto committed
509
protected:
510
511
    Graph& _g;
    EdgeIndexMap _edge_index;
Tiago Peixoto's avatar
Tiago Peixoto committed
512
    vector<edge_t>& _edges;
513
    rng_t& _rng;
514

515
    typedef gt_hash_map<size_t, size_t> nmapv_t;
516
517
518
519
    typedef typename property_map_type::apply<nmapv_t,
                                              typename property_map<Graph, vertex_index_t>::type>
        ::type::unchecked_t nmap_t;
    nmap_t _nmap;
520
521
};

522
523
// this will rewire the edges so that the combined (in, out) degree distribution
// will be the same, but all the rest is random
524
template <class Graph, class EdgeIndexMap, class CorrProb, class BlockDeg>
525
526
class RandomRewireStrategy:
    public RewireStrategyBase<Graph, EdgeIndexMap,
527
                              RandomRewireStrategy<Graph, EdgeIndexMap,
528
                                                   CorrProb, BlockDeg> >
529
530
{
public:
531
    typedef RewireStrategyBase<Graph, EdgeIndexMap,
532
                               RandomRewireStrategy<Graph, EdgeIndexMap,
533
                                                    CorrProb, BlockDeg> >
534
535
536
537
538
        base_t;

    typedef Graph graph_t;
    typedef EdgeIndexMap edge_index_t;

539
540
    typedef typename graph_traits<Graph>::vertex_descriptor vertex_t;
    typedef typename graph_traits<Graph>::edge_descriptor edge_t;
541
    typedef typename EdgeIndexMap::value_type index_t;
542

Tiago Peixoto's avatar
Tiago Peixoto committed
543
544
545
    struct hash_index {};
    struct random_index {};

546
    RandomRewireStrategy(Graph& g, EdgeIndexMap edge_index,
547
                         vector<edge_t>& edges, CorrProb, BlockDeg,
548
549
                         bool, rng_t& rng, bool parallel_edges)
        : base_t(g, edge_index, edges, rng, parallel_edges), _g(g) {}
Tiago Peixoto's avatar
Tiago Peixoto committed
550

551
    pair<size_t,bool> get_target_edge(pair<size_t,bool>& e, bool)
Tiago Peixoto's avatar
Tiago Peixoto committed
552
    {
Tiago Peixoto's avatar
Tiago Peixoto committed
553
        std::uniform_int_distribution<> sample(0, base_t::_edges.size() - 1);
554
555
        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
556
        {
Tiago Peixoto's avatar
Tiago Peixoto committed
557
            std::bernoulli_distribution coin(0.5);
558
            et.second = coin(base_t::_rng);
559
            e.second = coin(base_t::_rng);
Tiago Peixoto's avatar
Tiago Peixoto committed
560
        }
561
        return et;
Tiago Peixoto's avatar
Tiago Peixoto committed
562
563
    }

564
    void update_edge(size_t, bool) {}
565

566
private:
Tiago Peixoto's avatar
Tiago Peixoto committed
567
568
    Graph& _g;
    EdgeIndexMap _edge_index;
569
};
570

571
572
573

// 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
574
template <class Graph, class EdgeIndexMap, class CorrProb, class BlockDeg>
575
576
class CorrelatedRewireStrategy:
    public RewireStrategyBase<Graph, EdgeIndexMap,
577
                              CorrelatedRewireStrategy<Graph, EdgeIndexMap,
578
                                                       CorrProb, BlockDeg> >
579
580
{
public:
581
    typedef RewireStrategyBase<Graph, EdgeIndexMap,
582
                               CorrelatedRewireStrategy<Graph, EdgeIndexMap,
583
                                                        CorrProb, BlockDeg> >
584
585
586
587
        base_t;

    typedef Graph graph_t;

588
589
590
    typedef typename graph_traits<Graph>::vertex_descriptor vertex_t;
    typedef typename graph_traits<Graph>::edge_descriptor edge_t;

591
592
    typedef typename BlockDeg::block_t deg_t;

593
    CorrelatedRewireStrategy(Graph& g, EdgeIndexMap edge_index,
594
                             vector<edge_t>& edges, CorrProb, BlockDeg blockdeg,
595
596
                             bool, rng_t& rng, bool parallel_edges)
        : base_t(g, edge_index, edges, rng, parallel_edges), _blockdeg(blockdeg), _g(g)
597
    {
598
        for (size_t ei = 0; ei < base_t::_edges.size(); ++ei)
599
        {
600
            // For undirected graphs, there is no difference between source and
Tiago Peixoto's avatar
Tiago Peixoto committed
601
            // target, and each edge will appear _twice_ in the list below,
602
            // once for each different ordering of source and target.
603
            edge_t& e = base_t::_edges[ei];
604

Tiago Peixoto's avatar
Tiago Peixoto committed
605
            vertex_t t = target(e, _g);
606
            deg_t tdeg = get_deg(t, _g);;
Tiago Peixoto's avatar
Tiago Peixoto committed
607
            _edges_by_target[tdeg].push_back(make_pair(ei, false));
608
609
610

            if (!is_directed::apply<Graph>::type::value)
            {
Tiago Peixoto's avatar
Tiago Peixoto committed
611
                t = source(e, _g);
612
                deg_t tdeg = get_deg(t, _g);
Tiago Peixoto's avatar
Tiago Peixoto committed
613
                _edges_by_target[tdeg].push_back(make_pair(ei, true));
614
615
616
            }
        }
    }
617

618
    pair<size_t,bool> get_target_edge(pair<size_t, bool>& e, bool)
619
    {
620
621
622
623
624
        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
625

626
627
628
629
630
631
632
633
        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
634
    }
635

636
    void update_edge(size_t, bool) {}
Tiago Peixoto's avatar
Tiago Peixoto committed
637

638
639
640
641
642
    deg_t get_deg(vertex_t v, const Graph& g)
    {
        return _blockdeg.get_block(v, g);
    }

643
private:
644
645
    BlockDeg _blockdeg;

646
647
    typedef std::unordered_map<deg_t,
                          vector<pair<size_t, bool>>>
648
649
        edges_by_end_deg_t;

650
    edges_by_end_deg_t _edges_by_target;
651
652
653

protected:
    const Graph& _g;
654
655
};

656
657

// general stochastic blockmodel
658
// this version is based on rejection sampling
659
template <class Graph, class EdgeIndexMap, class CorrProb, class BlockDeg>
660
661
662
class ProbabilisticRewireStrategy:
    public RewireStrategyBase<Graph, EdgeIndexMap,
                              ProbabilisticRewireStrategy<Graph, EdgeIndexMap,
663
                                                          CorrProb, BlockDeg> >
664
665
666
{
public:
    typedef RewireStrategyBase<Graph, EdgeIndexMap,
667
                               ProbabilisticRewireStrategy<Graph, EdgeIndexMap,
668
                                                           CorrProb, BlockDeg> >
669
670
671
672
673
        base_t;

    typedef Graph graph_t;
    typedef EdgeIndexMap edge_index_t;

674
    typedef typename BlockDeg::block_t deg_t;
675
676
677
678
679
680

    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,
681
                                vector<edge_t>& edges, CorrProb corr_prob,
682
683
684
                                BlockDeg blockdeg, bool cache, rng_t& rng,
                                bool parallel_edges)
        : base_t(g, edge_index, edges, rng, parallel_edges), _g(g), _corr_prob(corr_prob),
685
          _blockdeg(blockdeg)
686
    {
687
688
689
        if (cache)
        {
            // cache probabilities
690
            _corr_prob.get_probs(_probs);
691

692
693
            if (_probs.empty())
            {
694
                std::unordered_set<deg_t> deg_set;
695
                for (size_t ei = 0; ei < base_t::_edges.size(); ++ei)
696
                {
697
698
699
                    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));
700
                }
701
702
703
704
705
706
707
708
709
710
711
712

                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;
713
714
715
                // avoid zero probability to not get stuck in rejection step
                if (std::isnan(p) || std::isinf(p) || p <= 0)
                    p = numeric_limits<double>::min();
716
717
                p = log(p);
            }
718
719
720
721
722
723
724
725
        }
    }

    double get_prob(const deg_t& s_deg, const deg_t& t_deg)
    {
        if (_probs.empty())
        {
            double p = _corr_prob(s_deg, t_deg);
726
            // avoid zero probability to not get stuck in rejection step
727
            if (std::isnan(p) || std::isinf(p) || p <= 0)
728
                p = numeric_limits<double>::min();
729
            return log(p);
730
        }
731
732
733
734
735
        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;
736
737
738
739
    }

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

743
    pair<size_t, bool> get_target_edge(pair<size_t, bool>& e, bool)
744
    {
745
746
747
748
749
750
751
752
        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
753

Tiago Peixoto's avatar
Tiago Peixoto committed
754
        std::uniform_int_distribution<> sample(0, base_t::_edges.size() - 1);
755
756
757
758
759
        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
760
            std::bernoulli_distribution coin(0.5);
761
762
            ep.second = coin(base_t::_rng);
        }
763

764
765
766
767
        if (source(e, base_t::_edges, _g) == source(ep, base_t::_edges, _g) ||
            target(e, base_t::_edges, _g) == target(ep, base_t::_edges, _g))
            return ep; // rewiring is a no-op

768
769
770
        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);

771
772
        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
773

774
        if (pf >= pi)
775
            return ep;
776

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

Tiago Peixoto's avatar
Tiago Peixoto committed
779
780
        std::uniform_real_distribution<> rsample(0.0, 1.0);
        double r = rsample(base_t::_rng);
781
        if (r > a)
782
            return e; // reject
783
        else
784
            return ep;
785
786
    }

787
    void update_edge(size_t, bool) {}
788

789
790
791
private:
    Graph& _g;
    EdgeIndexMap _edge_index;
792
    CorrProb _corr_prob;
793
    BlockDeg _blockdeg;
794

795
    typedef std::unordered_map<pair<deg_t, deg_t>, double> prob_map_t;
796
    prob_map_t _probs;
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

// 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,
825
826
827
828
                                     BlockDeg blockdeg, bool, rng_t& rng,
                                     bool parallel_edges)
        : base_t(g, edge_index, edges, rng, parallel_edges), _g(g),
          _corr_prob(corr_prob), _blockdeg(blockdeg)
829
    {
830
831
832
833
        _in_pos.resize(base_t::_edges.size());
        if (!is_directed::apply<Graph>::type::value)
            _out_pos.resize(base_t::_edges.size());

834
        std::unordered_set<deg_t> deg_set;
835
836
837
838
839
        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));
840
841
842
843
844
845
846
847
848
849
850
851
852

            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;
            }
853
854
        }

855
        _corr_prob.get_probs(_probs);
856

857
        if (_probs.empty())
858
        {
859
860
861
862
863
864
            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)
865
            {
866
867
868
869
870
871
872
873
874
875
876
877
                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;
                }
878

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

881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
                auto& ps = _sprob[*s_iter];

                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));
                }
            }
        }
        else
897
        {
898
899
            std::unordered_map<deg_t, vector<double>> sprobs;
            std::unordered_map<deg_t, vector<deg_t>> sitems;
900
901
902
903
904
905
906
907
908
909
910
            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);
            }
911

912
913

            for (auto iter = sitems.begin(); iter != sitems.end(); ++iter)
914
            {
915
916
917
918
919
920
921
922
923
924
925
926
927
928
                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)
                    sum += sprobs[s][i];

                auto& pr = _sprob[s];

                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));
                }
929
930
931
932
933
934
            }
        }
    }

    ~AliasProbabilisticRewireStrategy()
    {
935
        for (decltype(_sampler.begin()) iter = _sampler.begin();
936
937
938
939
940
941
             iter != _sampler.end(); ++iter)
            delete iter->second;
    }

    double get_prob(const deg_t& s_deg, const deg_t& t_deg)
    {
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
        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;
957
958
959
960
961
962
963
    }

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

964
    pair<size_t, bool> get_target_edge(pair<size_t, bool>& e, bool)
965
    {
966
967
968
969
970
971
972
973
974
975
976
        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);
977

978
979
980
981
982
        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);
983

984
985
986
        if (_in_edges[nt].empty() && _out_edges[nt].empty())
            return e; // reject

987
        pair<size_t, bool> ep;
Tiago Peixoto's avatar
Tiago Peixoto committed
988
        std::bernoulli_distribution coin(_in_edges[nt].size() /
989
990
991
992
993
994
                                         double(_in_edges[nt].