graph_rewiring.hh 43.5 KB
Newer Older
1
2
// graph-tool -- a general graph modification and manipulation thingy
//
Tiago Peixoto's avatar
Tiago Peixoto committed
3
// Copyright (C) 2006-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
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

    template <class Graph, class EdgeIndexMap, class CorrProb,
291
              class BlockDeg, class PinMap>
292
    void operator()(Graph& g, EdgeIndexMap edge_index, CorrProb corr_prob,
293
                    PinMap pin, 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
310
            if (pin[*e])
                continue;
311
            edges.push_back(*e);
Tiago Peixoto's avatar
Tiago Peixoto committed
312
313
            edge_pos.push_back(edge_pos.size());
        }
314

315
316
317
318
        typedef random_permutation_iterator<typename vector<size_t>::iterator,
                                            rng_t>
            random_edge_iter;

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

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

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

                size_t e = *ei;

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

351
352
353
354
355
356
                if (!success)
                    ++pcount;

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

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

Tiago Peixoto's avatar
Tiago Peixoto committed
375

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

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

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

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

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

        return true;
422
423
424
425
426
    }

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


433

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

444
445
    typedef typename EdgeIndexMap::value_type index_t;

Tiago Peixoto's avatar
Tiago Peixoto committed
446
    RewireStrategyBase(Graph& g, EdgeIndexMap edge_index, vector<edge_t>& edges,
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
                       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);
        }
    }
466

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

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

474
475
        pair<size_t, bool> e = make_pair(ei, false);

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

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

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

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

498
499
500
501
502
503
            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);
            }

504
            swap_edge::swap_target(e, et, _edges, _g);
505

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

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

520
        return true;
521
522
    }

Tiago Peixoto's avatar
Tiago Peixoto committed
523
protected:
524
525
    Graph& _g;
    EdgeIndexMap _edge_index;
Tiago Peixoto's avatar
Tiago Peixoto committed
526
    vector<edge_t>& _edges;
527
    rng_t& _rng;
528
529

#ifdef HAVE_SPARSEHASH
530
    typedef google::dense_hash_map<size_t, size_t, std::hash<size_t>> nmapv_t;
531
#else
532
    typedef unordered_map<size_t, size_t> nmapv_t;
533
534
535
536
537
#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;
538
539
};

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

    typedef Graph graph_t;
    typedef EdgeIndexMap edge_index_t;

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

Tiago Peixoto's avatar
Tiago Peixoto committed
561
562
563
    struct hash_index {};
    struct random_index {};

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

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

582
    void update_edge(size_t, bool) {}
583

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

589
590
591

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

    typedef Graph graph_t;

606
607
608
    typedef typename graph_traits<Graph>::vertex_descriptor vertex_t;
    typedef typename graph_traits<Graph>::edge_descriptor edge_t;

609
610
    typedef typename BlockDeg::block_t deg_t;

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

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

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

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

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

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

658
    void update_edge(size_t, bool) {}
Tiago Peixoto's avatar
Tiago Peixoto committed
659

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

665
private:
666
667
668
669
670
    BlockDeg _blockdeg;

#ifdef HAVE_SPARSEHASH
    typedef google::dense_hash_map<deg_t,
                                   vector<pair<size_t, bool>>,
671
                                   std::hash<deg_t>>
672
673
        edges_by_end_deg_t;
#else
Tiago Peixoto's avatar
Tiago Peixoto committed
674
    typedef std::unordered_map<deg_t,
675
                               vector<pair<size_t, bool>>>
676
        edges_by_end_deg_t;
677
678
#endif

679
    edges_by_end_deg_t _edges_by_target;
680
681
682

protected:
    const Graph& _g;
683
684
};

685
686

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

    typedef Graph graph_t;
    typedef EdgeIndexMap edge_index_t;

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

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

724
725
            if (_probs.empty())
            {
726
                std::unordered_set<deg_t> deg_set;
727
                for (size_t ei = 0; ei < base_t::_edges.size(); ++ei)
728
                {
729
730
731
                    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));
732
                }
733
734
735
736
737
738
739
740
741
742
743
744
745
746

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

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

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

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

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

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

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

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

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

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

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

817
    void update_edge(size_t, bool) {}
818

819
820
821
private:
    Graph& _g;
    EdgeIndexMap _edge_index;
822
    CorrProb _corr_prob;
823
    BlockDeg _blockdeg;
824
825
826
827


#ifdef HAVE_SPARSEHASH
    typedef google::dense_hash_map<pair<deg_t, deg_t>, double,
828
                                   std::hash<pair<deg_t, deg_t>>> prob_map_t;
829
#else
830
    typedef std::unordered_map<pair<deg_t, deg_t>, double> prob_map_t;
831
832
833
#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