graph_rewiring.hh 43.3 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
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
template <class Nmap, class Graph>
void add_count(size_t s, size_t t, Nmap& nvmap, Graph& g)
{
    if (!is_directed::apply<Graph>::type::value && s > t)
        std::swap(s, t);
    auto& nmap = nvmap[s];
    nmap[t]++;
}

template <class Nmap, class Graph>
void remove_count(size_t s, size_t t, Nmap& nvmap, Graph& g)
{
    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>
size_t get_count(size_t s, size_t t, Nmap& nvmap, Graph& g)
{
    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;
}

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

        typename graph_traits<Graph>::vertex_descriptor
127
128
            s = source(e, edges, g),          // current source
            t = target(e, edges, g),          // current target
129
130
            nt = target(te, edges, g),        // new target
            te_s = source(te, edges, g);      // target edge source
131

132
        if (get_count(s,  nt, nmap, g) > 0)
133
            return true; // e would clash with an existing edge
134
        if (get_count(te_s, t, nmap, g) > 0)
135
            return true; // te would clash with an existing edge
136
        return false; // no parallel edges
137
138
    }

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

152
        if (e.first == te.first)
153
            return;
154

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

165
166
167
168
169
170
        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)
171
172
173
174
            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;
175
    }
176

177
178
};

Tiago Peixoto's avatar
Tiago Peixoto committed
179
// used for verbose display
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
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;
    }
}

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
//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,
224
                      const Graph&) const
225
226
227
228
229
230
231
    {
        return get(_p, v);
    }

private:
    PropertyMap _p;
};
Tiago Peixoto's avatar
Tiago Peixoto committed
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
278
279
280
281
282
// 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>()());
    }
};


283
// main rewire loop
284
285
template <template <class Graph, class EdgeIndexMap, class CorrProb,
                    class BlockDeg>
286
          class RewireStrategy>
287
288
struct graph_rewire
{
289
290
291

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

304
        vector<edge_t> edges;
Tiago Peixoto's avatar
Tiago Peixoto committed
305
        vector<size_t> edge_pos;
306
307
        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
308
        {
309
            edges.push_back(*e);
Tiago Peixoto's avatar
Tiago Peixoto committed
310
311
            edge_pos.push_back(edge_pos.size());
        }
312

313
314
315
316
        typedef random_permutation_iterator<typename vector<size_t>::iterator,
                                            rng_t>
            random_edge_iter;

317
        RewireStrategy<Graph, EdgeIndexMap, CorrProb, BlockDeg>
318
            rewire(g, edge_index, edges, corr_prob, bd, cache, rng, parallel_edges);
319

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

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

                size_t e = *ei;

342
343
344
                bool success = false;
                do
                {
345
                    success = rewire(e, self_loops, parallel_edges);
346
347
348
                }
                while(persist && !success);

349
350
351
352
353
354
                if (!success)
                    ++pcount;

                if (no_sweep)
                    break;
            }
355
        }
356
357
        if (verbose)
            cout << endl;
358
    }
359
360
361
362

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

Tiago Peixoto's avatar
Tiago Peixoto committed
373

374
375
// this will rewire the edges so that the resulting graph will be entirely
// random (i.e. Erdos-Renyi)
376
template <class Graph, class EdgeIndexMap, class CorrProb, class BlockDeg>
377
378
379
380
381
382
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
383
    ErdosRewireStrategy(Graph& g, EdgeIndexMap edge_index,
384
                        vector<edge_t>& edges, CorrProb, BlockDeg,
385
                        bool, rng_t& rng, bool parallel_edges)
Tiago Peixoto's avatar
Tiago Peixoto committed
386
387
        : _g(g), _edge_index(edge_index), _edges(edges),
          _vertices(HardNumVertices()(g)), _rng(rng)
388
389
390
391
392
393
394
    {
        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;
    }

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

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

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

415
        remove_edge(_edges[ei], _g);
416
        edge_t ne = add_edge(s, t, _g).first;
417
        _edges[ei] = ne;
418
419

        return true;
420
421
422
423
424
    }

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


431

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

442
443
    typedef typename EdgeIndexMap::value_type index_t;

Tiago Peixoto's avatar
Tiago Peixoto committed
444
    RewireStrategyBase(Graph& g, EdgeIndexMap edge_index, vector<edge_t>& edges,
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
                       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)
        {
            typename graph_traits<Graph>::vertex_iterator v, v_end;
            for (tie(v, v_end) = boost::vertices(g); v != v_end; ++v)
            {
    #ifdef HAVE_SPARSEHASH
                _nmap[*v].set_empty_key(get_null_key<size_t>()());
                _nmap[*v].set_deleted_key(get_null_key<size_t>()() - 1);
    #endif
            }

            for (size_t i = 0; i < edges.size(); ++i)
                add_count(source(edges[i], g), target(edges[i], g), _nmap, g);
        }
    }
464

465
    bool operator()(size_t ei, bool self_loops, bool parallel_edges)
466
    {
Tiago Peixoto's avatar
Tiago Peixoto committed
467
        RewireStrategy& self = *static_cast<RewireStrategy*>(this);
468

469
470
        // try randomly drawn pairs of edges and check if they satisfy all the
        // consistency checks
471

472
473
        pair<size_t, bool> e = make_pair(ei, false);

474
        // rewire target
475
        pair<size_t, bool> et = self.get_target_edge(e, parallel_edges);
476
477

        if (!self_loops) // reject self-loops if not allowed
478
        {
479
480
            if((source(e, _edges, _g) == target(et, _edges, _g)) ||
               (target(e, _edges, _g) == source(et, _edges, _g)))
481
482
                return false;
        }
483

484
        // reject parallel edges if not allowed
485
        if (!parallel_edges && (et.first != e.first))
486
        {
487
            if (swap_edge::parallel_check_target(e, et, _edges, _nmap, _g))
488
                return false;
489
        }
490

491
        if (e.first != et.first)
492
        {
493
            self.update_edge(e.first, false);
494
            self.update_edge(et.first, false);
Tiago Peixoto's avatar
Tiago Peixoto committed
495

496
497
498
499
500
501
            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);
            }

502
            swap_edge::swap_target(e, et, _edges, _g);
503

504
            self.update_edge(e.first, true);
505
            self.update_edge(et.first, true);
506
507
508
509
510
511

            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);
            }
512
        }
513
514
515
516
        else
        {
            return false;
        }
Tiago Peixoto's avatar
Tiago Peixoto committed
517

518
        return true;
519
520
    }

Tiago Peixoto's avatar
Tiago Peixoto committed
521
protected:
522
523
    Graph& _g;
    EdgeIndexMap _edge_index;
Tiago Peixoto's avatar
Tiago Peixoto committed
524
    vector<edge_t>& _edges;
525
    rng_t& _rng;
526
527
528
529
530
531
532
533
534
535

#ifdef HAVE_SPARSEHASH
        typedef google::dense_hash_map<size_t, size_t> nmapv_t;
#else
        typedef unordered_map<size_t, size_t> nmapv_t;
#endif
    typedef typename property_map_type::apply<nmapv_t,
                                              typename property_map<Graph, vertex_index_t>::type>
        ::type::unchecked_t nmap_t;
    nmap_t _nmap;
536
537
};

538
539
// this will rewire the edges so that the combined (in, out) degree distribution
// will be the same, but all the rest is random
540
template <class Graph, class EdgeIndexMap, class CorrProb, class BlockDeg>
541
542
class RandomRewireStrategy:
    public RewireStrategyBase<Graph, EdgeIndexMap,
543
                              RandomRewireStrategy<Graph, EdgeIndexMap,
544
                                                   CorrProb, BlockDeg> >
545
546
{
public:
547
    typedef RewireStrategyBase<Graph, EdgeIndexMap,
548
                               RandomRewireStrategy<Graph, EdgeIndexMap,
549
                                                    CorrProb, BlockDeg> >
550
551
552
553
554
        base_t;

    typedef Graph graph_t;
    typedef EdgeIndexMap edge_index_t;

555
556
    typedef typename graph_traits<Graph>::vertex_descriptor vertex_t;
    typedef typename graph_traits<Graph>::edge_descriptor edge_t;
557
    typedef typename EdgeIndexMap::value_type index_t;
558

Tiago Peixoto's avatar
Tiago Peixoto committed
559
560
561
    struct hash_index {};
    struct random_index {};

562
    RandomRewireStrategy(Graph& g, EdgeIndexMap edge_index,
563
                         vector<edge_t>& edges, CorrProb, BlockDeg,
564
565
                         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
566

567
    pair<size_t,bool> get_target_edge(pair<size_t,bool>& e, bool parallel_edges)
Tiago Peixoto's avatar
Tiago Peixoto committed
568
    {
Tiago Peixoto's avatar
Tiago Peixoto committed
569
        std::uniform_int_distribution<> sample(0, base_t::_edges.size() - 1);
570
571
        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
572
        {
Tiago Peixoto's avatar
Tiago Peixoto committed
573
            std::bernoulli_distribution coin(0.5);
574
            et.second = coin(base_t::_rng);
575
            e.second = coin(base_t::_rng);
Tiago Peixoto's avatar
Tiago Peixoto committed
576
        }
577
        return et;
Tiago Peixoto's avatar
Tiago Peixoto committed
578
579
    }

580
    void update_edge(size_t, bool) {}
581

582
private:
Tiago Peixoto's avatar
Tiago Peixoto committed
583
584
    Graph& _g;
    EdgeIndexMap _edge_index;
585
};
586

587
588
589

// 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
590
template <class Graph, class EdgeIndexMap, class CorrProb, class BlockDeg>
591
592
class CorrelatedRewireStrategy:
    public RewireStrategyBase<Graph, EdgeIndexMap,
593
                              CorrelatedRewireStrategy<Graph, EdgeIndexMap,
594
                                                       CorrProb, BlockDeg> >
595
596
{
public:
597
    typedef RewireStrategyBase<Graph, EdgeIndexMap,
598
                               CorrelatedRewireStrategy<Graph, EdgeIndexMap,
599
                                                        CorrProb, BlockDeg> >
600
601
602
603
        base_t;

    typedef Graph graph_t;

604
605
606
    typedef typename graph_traits<Graph>::vertex_descriptor vertex_t;
    typedef typename graph_traits<Graph>::edge_descriptor edge_t;

607
608
    typedef typename BlockDeg::block_t deg_t;

609
    CorrelatedRewireStrategy(Graph& g, EdgeIndexMap edge_index,
610
                             vector<edge_t>& edges, CorrProb, BlockDeg blockdeg,
611
612
                             bool, rng_t& rng, bool parallel_edges)
        : base_t(g, edge_index, edges, rng, parallel_edges), _blockdeg(blockdeg), _g(g)
613
    {
614
615
616
617
#ifdef HAVE_SPARSEHASH
        _edges_by_target.set_empty_key(get_null_key<deg_t>()());
#endif

618
        for (size_t ei = 0; ei < base_t::_edges.size(); ++ei)
619
        {
620
            // For undirected graphs, there is no difference between source and
Tiago Peixoto's avatar
Tiago Peixoto committed
621
            // target, and each edge will appear _twice_ in the list below,
622
            // once for each different ordering of source and target.
623
            edge_t& e = base_t::_edges[ei];
624

Tiago Peixoto's avatar
Tiago Peixoto committed
625
            vertex_t t = target(e, _g);
626
            deg_t tdeg = get_deg(t, _g);;
Tiago Peixoto's avatar
Tiago Peixoto committed
627
            _edges_by_target[tdeg].push_back(make_pair(ei, false));
628
629
630

            if (!is_directed::apply<Graph>::type::value)
            {
Tiago Peixoto's avatar
Tiago Peixoto committed
631
                t = source(e, _g);
632
                deg_t tdeg = get_deg(t, _g);
Tiago Peixoto's avatar
Tiago Peixoto committed
633
                _edges_by_target[tdeg].push_back(make_pair(ei, true));
634
635
636
            }
        }
    }
637

638
    pair<size_t,bool> get_target_edge(pair<size_t, bool>& e, bool parallel_edge)
639
    {
640
641
642
643
644
        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
645

646
647
648
649
650
651
652
653
        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
654
    }
655

656
    void update_edge(size_t, bool) {}
Tiago Peixoto's avatar
Tiago Peixoto committed
657

658
659
660
661
662
    deg_t get_deg(vertex_t v, const Graph& g)
    {
        return _blockdeg.get_block(v, g);
    }

663
private:
664
665
666
667
668
669
670
671
    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
672
    typedef std::unordered_map<deg_t,
673
674
                               vector<pair<size_t, bool>>,
                               boost::hash<deg_t>>
675
        edges_by_end_deg_t;
676
677
#endif

678
    edges_by_end_deg_t _edges_by_target;
679
680
681

protected:
    const Graph& _g;
682
683
};

684
685

// general stochastic blockmodel
686
// this version is based on rejection sampling
687
template <class Graph, class EdgeIndexMap, class CorrProb, class BlockDeg>
688
689
690
class ProbabilisticRewireStrategy:
    public RewireStrategyBase<Graph, EdgeIndexMap,
                              ProbabilisticRewireStrategy<Graph, EdgeIndexMap,
691
                                                          CorrProb, BlockDeg> >
692
693
694
{
public:
    typedef RewireStrategyBase<Graph, EdgeIndexMap,
695
                               ProbabilisticRewireStrategy<Graph, EdgeIndexMap,
696
                                                           CorrProb, BlockDeg> >
697
698
699
700
701
        base_t;

    typedef Graph graph_t;
    typedef EdgeIndexMap edge_index_t;

702
    typedef typename BlockDeg::block_t deg_t;
703
704
705
706
707
708

    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,
709
                                vector<edge_t>& edges, CorrProb corr_prob,
710
711
712
                                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),
713
          _blockdeg(blockdeg)
714
    {
715
716
717
#ifdef HAVE_SPARSEHASH
        _probs.set_empty_key(get_null_key<pair<deg_t, deg_t>>()());
#endif
718
719
720
        if (cache)
        {
            // cache probabilities
721
            _corr_prob.get_probs(_probs);
722

723
724
725
726
            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)
727
                {
728
729
730
                    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));
731
                }
732
733
734
735
736
737
738
739
740
741
742
743
744
745

                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);
            }
746
747
748
749
750
751
752
753
        }
    }

    double get_prob(const deg_t& s_deg, const deg_t& t_deg)
    {
        if (_probs.empty())
        {
            double p = _corr_prob(s_deg, t_deg);
754
            // avoid zero probability to not get stuck in rejection step
755
            if (std::isnan(p) || std::isinf(p) || p <= 0)
756
                p = numeric_limits<double>::min();
757
            return log(p);
758
        }
759
760
761
762
763
        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;
764
765
766
767
    }

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

771
772
    pair<size_t, bool> get_target_edge(pair<size_t, bool>& e,
                                       bool parallel_edges)
773
    {
774
775
776
777
778
779
780
781
        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
782

Tiago Peixoto's avatar
Tiago Peixoto committed
783
        std::uniform_int_distribution<> sample(0, base_t::_edges.size() - 1);
784
785
786
787
788
        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
789
            std::bernoulli_distribution coin(0.5);
790
791
            ep.second = coin(base_t::_rng);
        }
792

793
794
795
796
        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

797
798
799
        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);

800
801
        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
802

803
        if (pf >= pi)
804
            return ep;
805

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

Tiago Peixoto's avatar
Tiago Peixoto committed
808
809
        std::uniform_real_distribution<> rsample(0.0, 1.0);
        double r = rsample(base_t::_rng);
810
        if (r > a)
811
            return e; // reject
812
        else
813
            return ep;
814
815
    }

816
    void update_edge(size_t, bool) {}
817

818
819
820
private:
    Graph& _g;
    EdgeIndexMap _edge_index;
821
    CorrProb _corr_prob;
822
    BlockDeg _blockdeg;
823
824
825
826
827
828
829
830
831
832
833


#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;
834
};
835

836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861

// 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,
862
863
864
865
                                     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)
866
    {
867
868
869
870
871
872
873
874
875
876
877
878
879

#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
880
        std::unordered_set<deg_t> deg_set;
881
882
883
884
885
        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));
886
887
888
889
890
891
892
893
894
895
896
897
898

            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;
            }
899
900
        }

901
        _corr_prob.get_probs(_probs);
902

903
        if (_probs.empty())
904
        {
905
906
907
908
909
910
            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)
911
            {
912
913
914
915
916
917
918
919
920
921
922
923
                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;
                }
924

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

927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
                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));
                }
            }
        }
        else
946
        {
947
948
949
950
951
952
953
954
955
956
957
958
959
            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);
            }
960

961
962

            for (auto iter = sitems.begin(); iter != sitems.end(); ++iter)
963
            {
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
                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));
                }
985
986
987
988
989
990
991
992
993
994
995
996
997
            }
        }
    }

    ~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)
    {
998
999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
        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;
1013
1014
1015
1016
1017
1018
1019
    }

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

1020
1021
    pair<size_t, bool> get_target_edge(pair<size_t, bool>& e,
                                       bool parallel_edges)
1022
    {
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
        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);
1034

1035
1036
1037
1038
1039
        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);
1040

1041
1042
1043
        if (_in_edges[nt].empty() && _out_edges[nt].empty())
            return e; // reject

1044
        pair<size_t, bool> ep;
Tiago Peixoto's avatar
Tiago Peixoto committed
1045
        std::bernoulli_distribution coin(_in_edges[nt].size() /
1046
1047
1048
1049
1050
1051
                                         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
1052
            std::uniform_int_distribution<> sample(0, ies.size() - 1);
1053
1054
1055
1056
1057
1058
1059
1060
            size_t epi = ies[sample(base_t::_rng)];

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