graph_rewiring.hh 43.4 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
template <class Nmap, class Graph>
70
void add_count(size_t s, size_t t, Nmap& nvmap, Graph&)
71
72
73
74
75
76
77
78
{
    if (!is_directed::apply<Graph>::type::value && s > t)
        std::swap(s, t);
    auto& nmap = nvmap[s];
    nmap[t]++;
}

template <class Nmap, class Graph>
79
void remove_count(size_t s, size_t t, Nmap& nvmap, Graph&)
80
81
82
83
84
85
86
87
88
89
90
{
    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>
91
size_t get_count(size_t s, size_t t, Nmap& nvmap, Graph&)
92
93
94
95
96
97
98
99
100
101
102
103
104
105
{
    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)
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)
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)
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
    pair<size_t, bool> get_target_edge(pair<size_t, bool>& e, bool)
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


#ifdef HAVE_SPARSEHASH
    typedef google::dense_hash_map<pair<deg_t, deg_t>, double,
827
                                   std::hash<pair<deg_t, deg_t>>> prob_map_t;
828
#else
829
    typedef std::unordered_map<pair<deg_t, deg_t>, double> prob_map_t;
830
831
832
#endif

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

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

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

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

900
        _corr_prob.get_probs(_probs);
901

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

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

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

960
961

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

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

    deg_t get_deg(vertex_t v, Graph& g)
    {
        return