graph_rewiring.hh 34.1 KB
Newer Older
1
2
// graph-tool -- a general graph modification and manipulation thingy
//
Tiago Peixoto's avatar
Tiago Peixoto committed
3
// Copyright (C) 2006-2017 Tiago de Paula Peixoto <tiago@skewed.de>
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
//
// This program is free software; you can redistribute it and/or
// modify it under the terms of the GNU General Public License
// as published by the Free Software Foundation; either version 3
// of the License, or (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <http://www.gnu.org/licenses/>.

#ifndef GRAPH_REWIRING_HH
#define GRAPH_REWIRING_HH

Tiago Peixoto's avatar
Tiago Peixoto committed
21
#include <tuple>
Tiago Peixoto's avatar
Tiago Peixoto committed
22
#include <iostream>
23
#include <boost/functional/hash.hpp>
Tiago Peixoto's avatar
Tiago Peixoto committed
24

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

30
31
#include "random.hh"

32
#include "hash_map_wrap.hh"
33

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

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

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

Tiago Peixoto's avatar
Tiago Peixoto committed
63

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

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

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

101
102
// this functor will swap the source of the edge e with the source of edge se
// and the target of edge e with the target of te
Tiago Peixoto's avatar
Tiago Peixoto committed
103
struct swap_edge
104
{
105
    template <class Nmap, class Graph>
Tiago Peixoto's avatar
Tiago Peixoto committed
106
    static bool
107
108
    parallel_check_target (const pair<size_t, bool>& e,
                           const pair<size_t, bool>& te,
109
                           vector<typename graph_traits<Graph>::edge_descriptor>& edges,
110
                           Nmap& nmap,
111
                           const Graph &g)
Tiago Peixoto's avatar
Tiago Peixoto committed
112
113
114
115
116
117
118
    {
        // We want to check that if we swap the target of 'e' with the target of
        // 'te', as such
        //
        //  (s)    -e--> (t)          (s)    -e--> (nt)
        //  (te_s) -te-> (nt)   =>    (te_s) -te-> (t)
        //
119
        // no parallel edges are introduced.
Tiago Peixoto's avatar
Tiago Peixoto committed
120
121

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

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

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

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

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

160
161
162
163
164
165
        if (is_directed::apply<Graph>::type::value || !e.second)
            ne = add_edge(s_e, t_te, g).first;
        else // keep invertedness (only for undirected graphs)
            ne = add_edge(t_te, s_e, g).first;
        edges[e.first] = ne;
        if (is_directed::apply<Graph>::type::value || !te.second)
166
167
168
169
            nte = add_edge(s_te, t_e, g).first;
        else // keep invertedness (only for undirected graphs)
            nte = add_edge(t_e, s_te,  g).first;
        edges[te.first] = nte;
170
    }
171

172
173
};

Tiago Peixoto's avatar
Tiago Peixoto committed
174
// used for verbose display
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
void print_progress(size_t i, size_t n_iter, size_t current, size_t total,
                    stringstream& str)
{
    size_t atom = (total > 200) ? total / 100 : 1;
    if ( ( (current+1) % atom == 0) || (current + 1) == total)
    {
        size_t size = str.str().length();
        for (size_t j = 0; j < str.str().length(); ++j)
            cout << "\b";
        str.str("");
        str << "(" << i + 1 << " / " << n_iter << ") "
            << current + 1 << " of " << total << " ("
            << (current + 1) * 100 / total << "%)";
        for (int j = 0; j < int(size - str.str().length()); ++j)
            str << " ";
        cout << str.str() << flush;
    }
}

194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
//select blocks based on in/out degrees
class DegreeBlock
{
public:
    typedef pair<size_t, size_t> block_t;

    template <class Graph>
    block_t get_block(typename graph_traits<Graph>::vertex_descriptor v,
                      const Graph& g) const
    {
        return make_pair(in_degreeS()(v, g), out_degree(v, g));
    }
};

//select blocks based on property map
template <class PropertyMap>
class PropertyBlock
{
public:
    typedef typename property_traits<PropertyMap>::value_type block_t;

    PropertyBlock(PropertyMap p): _p(p) {}

    template <class Graph>
    block_t get_block(typename graph_traits<Graph>::vertex_descriptor v,
219
                      const Graph&) const
220
221
222
223
224
225
226
    {
        return get(_p, v);
    }

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

228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
// select an appropriate "null" key for densehash
template <class Type>
struct get_null_key
{
    Type operator()() const
    {
        return numeric_limits<Type>::max();
    }
};

template <>
struct get_null_key<string>
{
    string operator()() const
    {
        return lexical_cast<string>(get_null_key<size_t>()());
    }
};

template <>
struct get_null_key<boost::python::object>
{
    boost::python::object operator()() const
    {
        return boost::python::object();
    }
};

template <class Type>
struct get_null_key<vector<Type>>
{
    vector<Type> operator()() const
    {
        vector<Type> v(1);
        v[0] = get_null_key<Type>()();
        return v;
    }
};

template <class Type1, class Type2>
struct get_null_key<pair<Type1, Type2>>
{
    pair<Type1, Type2> operator()() const
    {
        return make_pair(get_null_key<Type1>()(),
                         get_null_key<Type2>()());
    }
};


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

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

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

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

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

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

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

                size_t e = *ei;

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

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

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

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

Tiago Peixoto's avatar
Tiago Peixoto committed
370

371
372
// this will rewire the edges so that the resulting graph will be entirely
// random (i.e. Erdos-Renyi)
373
template <class Graph, class EdgeIndexMap, class CorrProb, class BlockDeg>
374
375
376
377
378
379
class ErdosRewireStrategy
{
public:
    typedef typename graph_traits<Graph>::edge_descriptor edge_t;
    typedef typename EdgeIndexMap::value_type index_t;

Tiago Peixoto's avatar
Tiago Peixoto committed
380
    ErdosRewireStrategy(Graph& g, EdgeIndexMap edge_index,
381
                        vector<edge_t>& edges, CorrProb, BlockDeg,
382
                        bool, rng_t& rng, bool, bool configuration)
Tiago Peixoto's avatar
Tiago Peixoto committed
383
        : _g(g), _edge_index(edge_index), _edges(edges),
384
385
386
          _vertices(HardNumVertices()(g)), _rng(rng),
          _configuration(configuration),
          _nmap(get(vertex_index, g), num_vertices(g))
387
    {
388
        decltype(_vertices.begin()) viter = _vertices.begin();
389
390
391
        typename graph_traits<Graph>::vertex_iterator v, v_end;
        for (tie(v, v_end) = vertices(_g); v != v_end; ++v)
            *(viter++) = *v;
392
393
394
395
396
397

        if (!configuration)
        {
            for (size_t i = 0; i < edges.size(); ++i)
                add_count(source(edges[i], g), target(edges[i], g), _nmap, g);
        }
398
399
    }

400
    bool operator()(size_t ei, bool self_loops, bool parallel_edges)
401
    {
402
403
404
405
406
407
        size_t e_s = source(_edges[ei], _g);
        size_t e_t = target(_edges[ei], _g);

        if (!is_directed::apply<Graph>::type::value && e_s > e_t)
            std::swap(e_s, e_t);

408
        //try randomly drawn pairs of vertices
Tiago Peixoto's avatar
Tiago Peixoto committed
409
        std::uniform_int_distribution<size_t> sample(0, _vertices.size() - 1);
410
411
412
413
414
        typename graph_traits<Graph>::vertex_descriptor s, t;
        while (true)
        {
            s = sample(_rng);
            t = sample(_rng);
415
416
417
418
            if (s == t)
            {
                if (!self_loops) // reject self-loops if not allowed
                    continue;
419

420
421
422
423
424
425
426
427
428
            }
            else if (!is_directed::apply<Graph>::type::value && self_loops)
            {
                // sample self-loops w/ correct probability for undirected
                // graphs
                std::bernoulli_distribution reject(.5);
                if (reject(_rng))
                    continue;
            }
429
430
            break;
        }
431

432
433
434
435
436
437
        if (!is_directed::apply<Graph>::type::value && s > t)
            std::swap(s, t);

        if (s == e_s && t == e_t)
            return false;

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

442
443
444
445
446
447
448
449
450
451
452
453
        if (!_configuration)
        {
            size_t m = get_count(s, t, _nmap, _g);
            size_t m_e = get_count(e_s, e_t, _nmap, _g);

            double a = (m + 1) / double(m_e);

            std::bernoulli_distribution accept(std::min(a, 1.));
            if (!accept(_rng))
                return false;
        }

454
        remove_edge(_edges[ei], _g);
455
        edge_t ne = add_edge(s, t, _g).first;
456
        _edges[ei] = ne;
457

458
459
460
461
462
463
        if (!_configuration)
        {
            remove_count(e_s, e_t, _nmap, _g);
            add_count(s, t, _nmap, _g);
        }

464
        return true;
465
466
467
468
469
    }

private:
    Graph& _g;
    EdgeIndexMap _edge_index;
Tiago Peixoto's avatar
Tiago Peixoto committed
470
    vector<edge_t>& _edges;
471
472
    vector<typename graph_traits<Graph>::vertex_descriptor> _vertices;
    rng_t& _rng;
473
474
475
476
477
478
    bool _configuration;
    typedef gt_hash_map<size_t, size_t> nmapv_t;
    typedef typename property_map_type::apply<nmapv_t,
                                              typename property_map<Graph, vertex_index_t>::type>
        ::type::unchecked_t nmap_t;
    nmap_t _nmap;
479
480
481
};


482

483
484
485
// 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.
486
487
template <class Graph, class EdgeIndexMap, class RewireStrategy>
class RewireStrategyBase
488
489
490
{
public:
    typedef typename graph_traits<Graph>::edge_descriptor edge_t;
491
    typedef typename graph_traits<Graph>::vertex_descriptor vertex_t;
492

493
494
    typedef typename EdgeIndexMap::value_type index_t;

Tiago Peixoto's avatar
Tiago Peixoto committed
495
    RewireStrategyBase(Graph& g, EdgeIndexMap edge_index, vector<edge_t>& edges,
496
                       rng_t& rng, bool parallel_edges, bool configuration)
497
        : _g(g), _edge_index(edge_index), _edges(edges), _rng(rng),
498
499
          _nmap(get(vertex_index, g), num_vertices(g)),
          _configuration(configuration)
500
    {
501
        if (!parallel_edges || !configuration)
502
503
504
505
506
        {
            for (size_t i = 0; i < edges.size(); ++i)
                add_count(source(edges[i], g), target(edges[i], g), _nmap, g);
        }
    }
507

508
    bool operator()(size_t ei, bool self_loops, bool parallel_edges)
509
    {
Tiago Peixoto's avatar
Tiago Peixoto committed
510
        RewireStrategy& self = *static_cast<RewireStrategy*>(this);
511

512
513
        // try randomly drawn pairs of edges and check if they satisfy all the
        // consistency checks
514

515
516
        pair<size_t, bool> e = make_pair(ei, false);

517
        // rewire target
518
        pair<size_t, bool> et = self.get_target_edge(e, parallel_edges);
519

520
521
522
523
524
525
526
527
        if (et.first == ei)
            return false;

        auto s = source(e,  _edges, _g);
        auto t = target(e,  _edges, _g);
        auto ts = source(et,  _edges, _g);
        auto tt = target(et,  _edges, _g);

528
        if (!self_loops) // reject self-loops if not allowed
529
        {
530
            if(s == tt || ts == t)
531
532
                return false;
        }
533

534
        // reject parallel edges if not allowed
535
        if (!parallel_edges && (et.first != e.first))
536
        {
537
            if (swap_edge::parallel_check_target(e, et, _edges, _nmap, _g))
538
                return false;
539
        }
540

541
542
543
544
545
546
547
548
549
        double a = 0;

        if (!is_directed::apply<Graph>::type::value)
        {
            a -= log(2 + (s == t) + (ts == tt));
            a += log(2 + (s == tt) + (ts == t));
        }

        if (!_configuration)
550
        {
551
552
553
554
555
556
            map<std::pair<size_t, size_t>, int> delta;

            delta[std::make_pair(s, t)] -= 1;
            delta[std::make_pair(ts, tt)] -= 1;
            delta[std::make_pair(s, tt)] += 1;
            delta[std::make_pair(ts, t)] += 1;
Tiago Peixoto's avatar
Tiago Peixoto committed
557

558
            for (auto& e_d : delta)
559
            {
560
561
562
563
564
565
566
                auto u = e_d.first.first;
                auto v = e_d.first.second;
                int d = e_d.second;
                size_t m = get_count(u,  v,  _nmap, _g);
                a -= lgamma(m + 1) - lgamma((m + 1) + d);
                if (!is_directed::apply<Graph>::type::value && u == v)
                    a += d * log(2);
567
568
            }

569
        }
570

571
572
573
        std::bernoulli_distribution accept(std::min(exp(a), 1.));
        if (!accept(_rng))
            return false;
574

575
576
577
578
579
580
581
        self.update_edge(e.first, false);
        self.update_edge(et.first, false);

        if (!parallel_edges || !_configuration)
        {
            remove_count(source(e, _edges, _g), target(e, _edges, _g), _nmap, _g);
            remove_count(source(et, _edges, _g), target(et, _edges, _g), _nmap, _g);
582
        }
583
584
585
586
587
588
589

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

        self.update_edge(e.first, true);
        self.update_edge(et.first, true);

        if (!parallel_edges || !_configuration)
590
        {
591
592
            add_count(source(e, _edges, _g), target(e, _edges, _g), _nmap, _g);
            add_count(source(et, _edges, _g), target(et, _edges, _g), _nmap, _g);
593
        }
Tiago Peixoto's avatar
Tiago Peixoto committed
594

595
        return true;
596
597
    }

Tiago Peixoto's avatar
Tiago Peixoto committed
598
protected:
599
600
    Graph& _g;
    EdgeIndexMap _edge_index;
Tiago Peixoto's avatar
Tiago Peixoto committed
601
    vector<edge_t>& _edges;
602
    rng_t& _rng;
603

604
    typedef gt_hash_map<size_t, size_t> nmapv_t;
605
606
607
608
    typedef typename property_map_type::apply<nmapv_t,
                                              typename property_map<Graph, vertex_index_t>::type>
        ::type::unchecked_t nmap_t;
    nmap_t _nmap;
609
    bool _configuration;
610
611
};

612
613
// this will rewire the edges so that the combined (in, out) degree distribution
// will be the same, but all the rest is random
614
template <class Graph, class EdgeIndexMap, class CorrProb, class BlockDeg>
615
616
class RandomRewireStrategy:
    public RewireStrategyBase<Graph, EdgeIndexMap,
617
                              RandomRewireStrategy<Graph, EdgeIndexMap,
618
                                                   CorrProb, BlockDeg> >
619
620
{
public:
621
    typedef RewireStrategyBase<Graph, EdgeIndexMap,
622
                               RandomRewireStrategy<Graph, EdgeIndexMap,
623
                                                    CorrProb, BlockDeg> >
624
625
626
627
628
        base_t;

    typedef Graph graph_t;
    typedef EdgeIndexMap edge_index_t;

629
630
    typedef typename graph_traits<Graph>::vertex_descriptor vertex_t;
    typedef typename graph_traits<Graph>::edge_descriptor edge_t;
631
    typedef typename EdgeIndexMap::value_type index_t;
632

Tiago Peixoto's avatar
Tiago Peixoto committed
633
634
635
    struct hash_index {};
    struct random_index {};

636
    RandomRewireStrategy(Graph& g, EdgeIndexMap edge_index,
637
                         vector<edge_t>& edges, CorrProb, BlockDeg,
638
639
640
641
                         bool, rng_t& rng, bool parallel_edges,
                         bool configuration)
        : base_t(g, edge_index, edges, rng, parallel_edges, configuration),
          _g(g) {}
Tiago Peixoto's avatar
Tiago Peixoto committed
642

643
    pair<size_t,bool> get_target_edge(pair<size_t,bool>& e, bool)
Tiago Peixoto's avatar
Tiago Peixoto committed
644
    {
Tiago Peixoto's avatar
Tiago Peixoto committed
645
        std::uniform_int_distribution<> sample(0, base_t::_edges.size() - 1);
646
647
        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
648
        {
Tiago Peixoto's avatar
Tiago Peixoto committed
649
            std::bernoulli_distribution coin(0.5);
650
            et.second = coin(base_t::_rng);
651
            e.second = coin(base_t::_rng);
Tiago Peixoto's avatar
Tiago Peixoto committed
652
        }
653
        return et;
Tiago Peixoto's avatar
Tiago Peixoto committed
654
655
    }

656
    void update_edge(size_t, bool) {}
657

658
private:
Tiago Peixoto's avatar
Tiago Peixoto committed
659
660
    Graph& _g;
    EdgeIndexMap _edge_index;
661
};
662

663
664
665

// 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
666
template <class Graph, class EdgeIndexMap, class CorrProb, class BlockDeg>
667
668
class CorrelatedRewireStrategy:
    public RewireStrategyBase<Graph, EdgeIndexMap,
669
                              CorrelatedRewireStrategy<Graph, EdgeIndexMap,
670
                                                       CorrProb, BlockDeg> >
671
672
{
public:
673
    typedef RewireStrategyBase<Graph, EdgeIndexMap,
674
                               CorrelatedRewireStrategy<Graph, EdgeIndexMap,
675
                                                        CorrProb, BlockDeg> >
676
677
678
679
        base_t;

    typedef Graph graph_t;

680
681
682
    typedef typename graph_traits<Graph>::vertex_descriptor vertex_t;
    typedef typename graph_traits<Graph>::edge_descriptor edge_t;

683
684
    typedef typename BlockDeg::block_t deg_t;

685
    CorrelatedRewireStrategy(Graph& g, EdgeIndexMap edge_index,
686
                             vector<edge_t>& edges, CorrProb, BlockDeg blockdeg,
687
688
689
690
                             bool, rng_t& rng, bool parallel_edges,
                             bool configuration)
        : base_t(g, edge_index, edges, rng, parallel_edges, configuration),
          _blockdeg(blockdeg), _g(g)
691
    {
692
        for (size_t ei = 0; ei < base_t::_edges.size(); ++ei)
693
        {
694
            // For undirected graphs, there is no difference between source and
Tiago Peixoto's avatar
Tiago Peixoto committed
695
            // target, and each edge will appear _twice_ in the list below,
696
            // once for each different ordering of source and target.
697
            edge_t& e = base_t::_edges[ei];
698

Tiago Peixoto's avatar
Tiago Peixoto committed
699
            vertex_t t = target(e, _g);
700
            deg_t tdeg = get_deg(t, _g);;
Tiago Peixoto's avatar
Tiago Peixoto committed
701
            _edges_by_target[tdeg].push_back(make_pair(ei, false));
702
703
704

            if (!is_directed::apply<Graph>::type::value)
            {
Tiago Peixoto's avatar
Tiago Peixoto committed
705
                t = source(e, _g);
706
                deg_t tdeg = get_deg(t, _g);
Tiago Peixoto's avatar
Tiago Peixoto committed
707
                _edges_by_target[tdeg].push_back(make_pair(ei, true));
708
709
710
            }
        }
    }
711

712
    pair<size_t,bool> get_target_edge(pair<size_t, bool>& e, bool)
713
    {
714
715
716
717
718
        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
719

720
721
722
723
724
725
726
727
        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
728
    }
729

730
    void update_edge(size_t, bool) {}
Tiago Peixoto's avatar
Tiago Peixoto committed
731

732
733
734
735
736
    deg_t get_deg(vertex_t v, const Graph& g)
    {
        return _blockdeg.get_block(v, g);
    }

737
private:
738
739
    BlockDeg _blockdeg;

740
741
    typedef std::unordered_map<deg_t,
                          vector<pair<size_t, bool>>>
742
743
        edges_by_end_deg_t;

744
    edges_by_end_deg_t _edges_by_target;
745
746
747

protected:
    const Graph& _g;
748
749
};

750
751

// general stochastic blockmodel
752
// this version is based on rejection sampling
753
template <class Graph, class EdgeIndexMap, class CorrProb, class BlockDeg>
754
755
756
class ProbabilisticRewireStrategy:
    public RewireStrategyBase<Graph, EdgeIndexMap,
                              ProbabilisticRewireStrategy<Graph, EdgeIndexMap,
757
                                                          CorrProb, BlockDeg> >
758
759
760
{
public:
    typedef RewireStrategyBase<Graph, EdgeIndexMap,
761
                               ProbabilisticRewireStrategy<Graph, EdgeIndexMap,
762
                                                           CorrProb, BlockDeg> >
763
764
765
766
767
        base_t;

    typedef Graph graph_t;
    typedef EdgeIndexMap edge_index_t;

768
    typedef typename BlockDeg::block_t deg_t;
769
770
771
772
773
774

    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,
775
                                vector<edge_t>& edges, CorrProb corr_prob,
776
                                BlockDeg blockdeg, bool cache, rng_t& rng,
777
778
779
                                bool parallel_edges, bool configuration)
        : base_t(g, edge_index, edges, rng, parallel_edges, configuration),
          _g(g), _corr_prob(corr_prob), _blockdeg(blockdeg)
780
    {
781
782
783
        if (cache)
        {
            // cache probabilities
784
            _corr_prob.get_probs(_probs);
785

786
787
            if (_probs.empty())
            {
788
                std::unordered_set<deg_t> deg_set;
789
                for (size_t ei = 0; ei < base_t::_edges.size(); ++ei)
790
                {
791
792
793
                    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));
794
                }
795
796
797
798
799
800
801
802
803
804
805
806

                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;
807
808
809
                // avoid zero probability to not get stuck in rejection step
                if (std::isnan(p) || std::isinf(p) || p <= 0)
                    p = numeric_limits<double>::min();
810
811
                p = log(p);
            }
812
813
814
815
816
817
818
819
        }
    }

    double get_prob(const deg_t& s_deg, const deg_t& t_deg)
    {
        if (_probs.empty())
        {
            double p = _corr_prob(s_deg, t_deg);
820
            // avoid zero probability to not get stuck in rejection step
821
            if (std::isnan(p) || std::isinf(p) || p <= 0)
822
                p = numeric_limits<double>::min();
823
            return log(p);
824
        }
825
826
827
828
829
        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;
830
831
832
833
    }

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

837
    pair<size_t, bool> get_target_edge(pair<size_t, bool>& e, bool)
838
    {
839
840
841
842
843
844
845
846
        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
847

Tiago Peixoto's avatar
Tiago Peixoto committed
848
        std::uniform_int_distribution<> sample(0, base_t::_edges.size() - 1);
849
850
851
852
853
        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
854
            std::bernoulli_distribution coin(0.5);
855
856
            ep.second = coin(base_t::_rng);
        }
857

858
859
860
861
        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

862
863
864
        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);

865
866
        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
867

868
        if (pf >= pi)
869
            return ep;
870

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

Tiago Peixoto's avatar
Tiago Peixoto committed
873
874
        std::uniform_real_distribution<> rsample(0.0, 1.0);
        double r = rsample(base_t::_rng);
875
        if (r > a)
876
            return e; // reject
877
        else
878
            return ep;
879
880
    }

881
    void update_edge(size_t, bool) {}
882

883
884
885
private:
    Graph& _g;
    EdgeIndexMap _edge_index;
886
    CorrProb _corr_prob;
887
    BlockDeg _blockdeg;
888

889
    typedef std::unordered_map<pair<deg_t, deg_t>, double> prob_map_t;
890
    prob_map_t _probs;
891
};
892

893
894
895

// general "traditional" stochastic blockmodel
// this version is based on the alias method, and does not keep the degrees fixed
896
897
template <class Graph, class EdgeIndexMap, class CorrProb, class BlockDeg,
          bool micro>
898
899
900
901
902
903
904
905
906
907
class TradBlockRewireStrategy
{
public:
    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;
    typedef typename BlockDeg::block_t deg_t;

    TradBlockRewireStrategy(Graph& g, EdgeIndexMap edge_index,
                            vector<edge_t>& edges, CorrProb corr_prob,
908
                            BlockDeg blockdeg, bool, rng_t& rng,
909
                            bool parallel_edges, bool configuration)
910
911

        : _g(g), _edge_index(edge_index), _edges(edges), _corr_prob(corr_prob),
912
913
914
          _blockdeg(blockdeg), _rng(rng), _sampler(nullptr),
          _configuration(configuration),
          _nmap(get(vertex_index, g), num_vertices(g))
915
    {
916
        for (auto v : vertices_range(_g))
917
        {
918
919
            deg_t d = _blockdeg.get_block(v, g);
            _vertices[d].push_back(v);
920
921
        }

922
        if (!micro)
923
        {
924
925
926
927
928
            std::unordered_map<pair<deg_t, deg_t>, double> probs;
            _corr_prob.get_probs(probs);

            vector<double> dprobs;
            if (probs.empty())
929
            {
930
                for (auto& s : _vertices)
931
                {
932
933
934
935
936
                    for (auto& t : _vertices)
                    {
                        double p = _corr_prob(s.first, t.first);
                        if (std::isnan(p) || std::isinf(p) || p <= 0)
                            continue;
937

938
939
940
                        _items.push_back(make_pair(s.first, t.first));
                        dprobs.push_back(p * s.second.size() * t.second.size());
                    }
941
                }
942
            }
943
            else
944
            {
945
946
947
948
949
950
951
952
953
954
955
                for (auto& stp : probs)
                {
                    deg_t s = stp.first.first;
                    deg_t t = stp.first.second;
                    double p = stp.second;
                    // avoid zero probability to not get stuck in rejection step
                    if (std::isnan(p) || std::isinf(p) || p <= 0)
                        continue;
                    _items.push_back(make_pair(s, t));
                    dprobs.push_back(p * _vertices[s].size() * _vertices[t].size());
                }
956
957
            }

958
959
            if (_items.empty())
                throw GraphException(