graph_rewiring.hh 40.7 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-2015 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
387
388
389
390
391
    {
        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;
    }

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
713
714

                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);
            }
715
716
717
718
719
720
721
722
        }
    }

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

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

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

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

761
762
763
764
        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

765
766
767
        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);

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

771
        if (pf >= pi)
772
            return ep;
773

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

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

784
    void update_edge(size_t, bool) {}
785

786
787
788
private:
    Graph& _g;
    EdgeIndexMap _edge_index;
789
    CorrProb _corr_prob;
790
    BlockDeg _blockdeg;
791
792


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

// 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,
823
824
825
826
                                     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)
827
    {
828
829
830
831
        _in_pos.resize(base_t::_edges.size());
        if (!is_directed::apply<Graph>::type::value)
            _out_pos.resize(base_t::_edges.size());

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

            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;
            }
851
852
        }

853
        _corr_prob.get_probs(_probs);
854

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

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

879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
                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
895
        {
896
897
            std::unordered_map<deg_t, vector<double>> sprobs;
            std::unordered_map<deg_t, vector<deg_t>> sitems;
898
899
900
901
902
903
904
905
906
907
908
            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);
            }
909

910
911

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

                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));
                }
930
931
932
933
934
935
936
937
938
939
940
941
942
            }
        }
    }

    ~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)
    {
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
        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;
958
959
960
961
962
963
964
    }

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

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

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

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

988
        pair<size_t, bool> ep;
Tiago Peixoto's avatar
Tiago Peixoto committed
989
        std::bernoulli_distribution coin(_in_edges[nt].size() /
990
991
992
993
994
995
                                         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
996
            std::uniform_int_distribution<> sample(0, ies.size() - 1);
997
998
999
1000
1001
1002
1003
1004
            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
1005