graph_rewiring.hh 33.9 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-2020 Tiago de Paula Peixoto <tiago@skewed.de>
4
//
5
6
7
8
// This program is free software; you can redistribute it and/or modify it under
// the terms of the GNU Lesser General Public License as published by the Free
// Software Foundation; either version 3 of the License, or (at your option) any
// later version.
9
//
10
11
12
13
// 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 Lesser General Public License for more
// details.
14
//
15
// You should have received a copy of the GNU Lesser General Public License
16
17
18
19
20
// 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
23
#include <iostream>

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

29
30
#include "random.hh"

31
#include "hash_map_wrap.hh"
32

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

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

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

Tiago Peixoto's avatar
Tiago Peixoto committed
62

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

template <class Nmap, class Graph>
73
void remove_count(size_t s, size_t t, Nmap& nvmap, Graph& g)
74
{
75
    if (!graph_tool::is_directed(g) && s > t)
76
77
78
79
80
81
82
83
84
        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>
85
size_t get_count(size_t s, size_t t, Nmap& nvmap, Graph& g)
86
{
87
    if (!graph_tool::is_directed(g) && s > t)
88
89
90
91
92
93
94
95
96
97
98
99
        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;
}

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

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

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

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

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

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

159
        if (graph_tool::is_directed(g) || !e.second)
160
161
162
163
            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;
164
        if (graph_tool::is_directed(g) || !te.second)
165
166
167
168
            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;
169
    }
170

171
172
};

Tiago Peixoto's avatar
Tiago Peixoto committed
173
// used for verbose display
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
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;
    }
}

193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
//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,
218
                      const Graph&) const
219
220
221
222
223
224
225
    {
        return get(_p, v);
    }

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

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
// 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>()());
    }
};


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

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

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

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

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

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

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

                size_t e = *ei;

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

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

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

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

Tiago Peixoto's avatar
Tiago Peixoto committed
369

370
371
// this will rewire the edges so that the resulting graph will be entirely
// random (i.e. Erdos-Renyi)
372
template <class Graph, class EdgeIndexMap, class CorrProb, class BlockDeg>
373
374
375
376
377
378
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
379
    ErdosRewireStrategy(Graph& g, EdgeIndexMap edge_index,
380
                        vector<edge_t>& edges, CorrProb, BlockDeg,
381
                        bool, rng_t& rng, bool, bool configuration)
Tiago Peixoto's avatar
Tiago Peixoto committed
382
        : _g(g), _edge_index(edge_index), _edges(edges),
383
384
385
          _vertices(HardNumVertices()(g)), _rng(rng),
          _configuration(configuration),
          _nmap(get(vertex_index, g), num_vertices(g))
386
    {
387
        decltype(_vertices.begin()) viter = _vertices.begin();
388
389
390
        typename graph_traits<Graph>::vertex_iterator v, v_end;
        for (tie(v, v_end) = vertices(_g); v != v_end; ++v)
            *(viter++) = *v;
391
392
393
394
395
396

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

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

404
        if (!graph_tool::is_directed(_g) && e_s > e_t)
405
406
            std::swap(e_s, e_t);

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

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

431
        if (!graph_tool::is_directed(_g) && s > t)
432
433
434
435
436
            std::swap(s, t);

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

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

441
442
443
444
445
446
447
448
449
450
451
452
        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;
        }

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

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

463
        return true;
464
465
466
467
468
    }

private:
    Graph& _g;
    EdgeIndexMap _edge_index;
Tiago Peixoto's avatar
Tiago Peixoto committed
469
    vector<edge_t>& _edges;
470
471
    vector<typename graph_traits<Graph>::vertex_descriptor> _vertices;
    rng_t& _rng;
472
473
474
475
476
477
    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;
478
479
480
};


481

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

492
493
    typedef typename EdgeIndexMap::value_type index_t;

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

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

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

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

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

519
520
521
522
523
524
525
526
        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);

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

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

540
541
        double a = 0;

542
        if (!graph_tool::is_directed(_g))
543
544
545
546
547
548
        {
            a -= log(2 + (s == t) + (ts == tt));
            a += log(2 + (s == tt) + (ts == t));
        }

        if (!_configuration)
549
        {
550
551
552
553
554
555
            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
556

557
            for (auto& e_d : delta)
558
            {
559
560
561
562
563
                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);
564
                if (!graph_tool::is_directed(_g) && u == v)
565
                    a += d * log(2);
566
567
            }

568
        }
569

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

574
575
576
577
578
579
580
        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);
581
        }
582
583
584
585
586
587
588

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

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

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

594
        return true;
595
596
    }

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

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

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

    typedef Graph graph_t;
    typedef EdgeIndexMap edge_index_t;

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

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

635
    RandomRewireStrategy(Graph& g, EdgeIndexMap edge_index,
636
                         vector<edge_t>& edges, CorrProb, BlockDeg,
637
638
639
640
                         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
641

642
    pair<size_t,bool> get_target_edge(pair<size_t,bool>& e, bool)
Tiago Peixoto's avatar
Tiago Peixoto committed
643
    {
Tiago Peixoto's avatar
Tiago Peixoto committed
644
        std::uniform_int_distribution<> sample(0, base_t::_edges.size() - 1);
645
        pair<size_t, bool> et = make_pair(sample(base_t::_rng), false);
646
        if (!graph_tool::is_directed(_g))
Tiago Peixoto's avatar
Tiago Peixoto committed
647
        {
Tiago Peixoto's avatar
Tiago Peixoto committed
648
            std::bernoulli_distribution coin(0.5);
649
            et.second = coin(base_t::_rng);
650
            e.second = coin(base_t::_rng);
Tiago Peixoto's avatar
Tiago Peixoto committed
651
        }
652
        return et;
Tiago Peixoto's avatar
Tiago Peixoto committed
653
654
    }

655
    void update_edge(size_t, bool) {}
656

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

662
663
664

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

    typedef Graph graph_t;

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

682
683
    typedef typename BlockDeg::block_t deg_t;

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

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

702
            if (!graph_tool::is_directed(_g))
703
            {
Tiago Peixoto's avatar
Tiago Peixoto committed
704
                t = source(e, _g);
705
                deg_t tdeg = get_deg(t, _g);
Tiago Peixoto's avatar
Tiago Peixoto committed
706
                _edges_by_target[tdeg].push_back(make_pair(ei, true));
707
708
709
            }
        }
    }
710

711
    pair<size_t,bool> get_target_edge(pair<size_t, bool>& e, bool)
712
    {
713
        if (!graph_tool::is_directed(_g))
714
715
716
717
        {
            std::bernoulli_distribution coin(0.5);
            e.second = coin(base_t::_rng);
        }
Tiago Peixoto's avatar
Tiago Peixoto committed
718

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

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

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

736
private:
737
738
    BlockDeg _blockdeg;

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

743
    edges_by_end_deg_t _edges_by_target;
744
745
746

protected:
    const Graph& _g;
747
748
};

749
750

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

    typedef Graph graph_t;
    typedef EdgeIndexMap edge_index_t;

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

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

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

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

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

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

836
    pair<size_t, bool> get_target_edge(pair<size_t, bool>& e, bool)
837
    {
838
        if (!graph_tool::is_directed(_g))
839
840
841
842
843
844
845
        {
            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
846

Tiago Peixoto's avatar
Tiago Peixoto committed
847
        std::uniform_int_distribution<> sample(0, base_t::_edges.size() - 1);
848
849
        size_t epi = sample(base_t::_rng);
        pair<size_t, bool> ep = make_pair(epi, false);
850
        if (!graph_tool::is_directed(_g))
851
852
        {
            // for undirected graphs we must select a random direction
Tiago Peixoto's avatar
Tiago Peixoto committed
853
            std::bernoulli_distribution coin(0.5);
854
855
            ep.second = coin(base_t::_rng);
        }
856

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

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

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

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

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

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

880
    void update_edge(size_t, bool) {}
881

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

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

892
893
894

// general "traditional" stochastic blockmodel
// this version is based on the alias method, and does not keep the degrees fixed
895
896
template <class Graph, class EdgeIndexMap, class CorrProb, class BlockDeg,
          bool micro>
897
898
899
900
901
902
903
904
905
906
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,
907
                            BlockDeg blockdeg, bool, rng_t& rng,
908
                            bool parallel_edges, bool configuration)
909
910

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

Tiago Peixoto's avatar