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

21
22
23
24
25
26
27
#if (GCC_VERSION >= 40400)
#   include <tr1/unordered_set>
#   include <tr1/random>
#else
#   include <boost/tr1/unordered_set.hpp>
#   include <boost/tr1/random.hpp>
#endif
28
#include <boost/functional/hash.hpp>
29

Tiago Peixoto's avatar
Tiago Peixoto committed
30
31
32
33
34
35
#include <iostream>

#include <boost/multi_index_container.hpp>
#include <boost/multi_index/hashed_index.hpp>
#include <boost/multi_index/random_access_index.hpp>
#include <boost/multi_index/identity.hpp>
36
37
38

#include "graph.hh"
#include "graph_filtering.hh"
39
#include "graph_util.hh"
40
#include "sampler.hh"
41
42
43
44
45

namespace graph_tool
{
using namespace std;
using namespace boost;
Tiago Peixoto's avatar
Tiago Peixoto committed
46
using namespace boost::multi_index;
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62

// returns true if vertices u and v are adjacent. This is O(k(u)).
template <class Graph>
bool is_adjacent(typename graph_traits<Graph>::vertex_descriptor u,
                 typename graph_traits<Graph>::vertex_descriptor v,
                 const Graph& g )
{
    typename graph_traits<Graph>::out_edge_iterator e, e_end;
    for (tie(e, e_end) = out_edges(u, g); e != e_end; ++e)
    {
        if (target(*e,g) == v)
            return true;
    }
    return false;
}

63
64
template <class Graph>
typename graph_traits<Graph>::vertex_descriptor
65
66
source(const pair<size_t, bool>& e,
       const vector<typename graph_traits<Graph>::edge_descriptor>& edges,
67
68
69
       const Graph& g)
{
    if (e.second)
70
        return target(edges[e.first], g);
71
    else
72
        return source(edges[e.first], g);
73
74
75
76
}

template <class Graph>
typename graph_traits<Graph>::vertex_descriptor
77
78
target(const pair<size_t, bool>& e,
       const vector<typename graph_traits<Graph>::edge_descriptor>& edges,
79
80
81
       const Graph& g)
{
    if (e.second)
82
        return source(edges[e.first], g);
83
    else
84
        return target(edges[e.first], g);
85
86
}

Tiago Peixoto's avatar
Tiago Peixoto committed
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
// This will iterate over a random permutation of a random access sequence, by
// swapping the values of the sequence as it iterates
template <class RandomAccessIterator, class RNG,
          class RandomDist = tr1::uniform_int<size_t> >
class random_permutation_iterator : public
    std::iterator<input_iterator_tag, typename RandomAccessIterator::value_type>
{
public:
    random_permutation_iterator(RandomAccessIterator begin,
                                RandomAccessIterator end, RNG& rng)
        : _i(begin), _end(end), _rng(&rng)
    {
        if(_i != _end)
        {
            RandomDist random(0,  _end - _i - 1);
            std::iter_swap(_i, _i + random(*_rng));
        }
    }

    typename RandomAccessIterator::value_type operator*()
    {
        return *_i;
    }

    random_permutation_iterator& operator++()
    {
        ++_i;
        if(_i != _end)
        {
            RandomDist random(0,  _end - _i - 1);
            std::iter_swap(_i, _i + random(*_rng));
        }
        return *this;
    }

    bool operator==(const random_permutation_iterator& ri)
    {
        return _i == ri._i;
    }

    bool operator!=(const random_permutation_iterator& ri)
    {
        return _i != ri._i;
    }

    size_t operator-(const random_permutation_iterator& ri)
    {
        return _i - ri._i;
    }

private:
    RandomAccessIterator _i, _end;
    RNG* _rng;
};

142
143
// 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
144
struct swap_edge
145
{
146
    template <class Graph>
147
    static bool
148
149
150
    parallel_check_source (size_t e, pair<size_t, bool> se,
                           vector<typename graph_traits<Graph>::edge_descriptor>& edges,
                           const Graph &g)
151
152
    {
        // We want to check that if we swap the source of 'e' with the source of
Tiago Peixoto's avatar
Tiago Peixoto committed
153
        // 'se', as such
154
        //
Tiago Peixoto's avatar
Tiago Peixoto committed
155
        //  (s)    -e--> (t)          (ns)   -e--> (t)
156
157
        //  (ns)   -se-> (se_t)   =>  (s)    -se-> (se_t)
        //
158
        // no parallel edges are introduced.
159
160

        typename graph_traits<Graph>::vertex_descriptor
161
162
163
164
            s = source(edges[e], g),          // current source
            t = target(edges[e], g),          // current target
            ns = source(se, edges, g),        // new source
            se_t = target(se, edges, g);      // source edge target
165

166
        if (is_adjacent(ns,  t, g))
Tiago Peixoto's avatar
Tiago Peixoto committed
167
            return true; // e would clash with an existing (new) edge
168
        if (is_adjacent(s, se_t, g))
Tiago Peixoto's avatar
Tiago Peixoto committed
169
170
171
172
            return true; // se would clash with an existing (new) edge
        return false; // the coast is clear - hooray!
    }

173
    template <class Graph>
Tiago Peixoto's avatar
Tiago Peixoto committed
174
    static bool
175
176
177
    parallel_check_target (size_t e, const pair<size_t, bool>& te,
                           vector<typename graph_traits<Graph>::edge_descriptor>& edges,
                           const Graph &g)
Tiago Peixoto's avatar
Tiago Peixoto committed
178
179
180
181
182
183
184
    {
        // 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)
        //
185
        // no parallel edges are introduced.
Tiago Peixoto's avatar
Tiago Peixoto committed
186
187

        typename graph_traits<Graph>::vertex_descriptor
188
189
190
191
            s = source(edges[e], g),          // current source
            t = target(edges[e], g),          // current target
            nt = target(te, edges, g),        // new target
            te_s = source(te, edges, g);      // target edge source
192

193
        if (is_adjacent(s,  nt, g))
194
            return true; // e would clash with an existing (new) edge
195
        if (is_adjacent(te_s, t, g))
196
197
198
199
            return true; // te would clash with an existing (new) edge
        return false; // the coast is clear - hooray!
    }

200
201
202
203
    template <class Graph, class EdgeIndexMap>
    static void swap_source (size_t e, pair<size_t, bool>& se,
                             vector<typename graph_traits<Graph>::edge_descriptor>& edges,
                             EdgeIndexMap edge_index, Graph& g)
204
    {
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
        // swap the source of the edge 'e' with the source of edge 'se', as
        // such:
        //
        //  (s)    -e--> (t)          (ns)   -e--> (t)
        //  (ns)   -se-> (se_t)   =>  (s)    -se-> (se_t)

        if (e == se.first)
            return;

        // new edges which will replace the old ones
        typename graph_traits<Graph>::edge_descriptor ne, nse;
        ne = add_edge(source(se, edges, g), target(edges[e], g), g).first;
        if (!se.second)
            nse = add_edge(source(edges[e], g), target(se, edges, g), g).first;
        else // keep invertedness (only for undirected graphs)
            nse = add_edge(target(se, edges, g), source(edges[e], g),  g).first;
        edge_index[nse] = se.first;
        remove_edge(edges[se.first], g);
        edges[se.first] = nse;

        edge_index[ne] = e;
        remove_edge(edges[e], g);
        edges[e] = ne;
228
229
    }

230
231
232
233
234
    template <class Graph, class EdgeIndexMap>
    static void swap_target
        (size_t e, const pair<size_t, bool>& te,
         vector<typename graph_traits<Graph>::edge_descriptor>& edges,
         EdgeIndexMap edge_index, Graph& g)
235
    {
236
237
        // swap the source of the edge 'e' with the source of edge 'se', as
        // such:
238
        //
239
240
        //  (s)    -e--> (t)          (s)    -e--> (nt)
        //  (te_s) -te-> (nt)   =>    (te_s) -te-> (t)
241

242
243
        if (e == te.first)
            return;
244

245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
        // new edges which will replace the old ones
        typename graph_traits<Graph>::edge_descriptor ne, nte;
        ne = add_edge(source(edges[e], g), target(te, edges, g), g).first;
        if (!te.second)
            nte = add_edge(source(te, edges, g), target(edges[e], g), g).first;
        else // keep invertedness (only for undirected graphs)
            nte = add_edge(target(edges[e], g), source(te, edges, g),  g).first;

        edge_index[nte] = te.first;
        remove_edge(edges[te.first], g);
        edges[te.first] = nte;

        edge_index[ne] = e;
        remove_edge(edges[e], g);
        edges[e] = ne;
260
    }
261

262
263
};

Tiago Peixoto's avatar
Tiago Peixoto committed
264
265
266
// used for verbose display
void print_progress(size_t current, size_t total, stringstream& str);

267
// main rewire loop
268
269
template <template <class Graph, class EdgeIndexMap, class CorrProb>
          class RewireStrategy>
270
271
struct graph_rewire
{
272
273
274
275
    template <class Graph, class EdgeIndexMap, class CorrProb>
    void operator()(Graph& g, EdgeIndexMap edge_index, CorrProb corr_prob,
                    rng_t& rng, bool self_loops, bool parallel_edges,
                    bool verbose) const
276
277
278
279
280
    {
        typedef typename graph_traits<Graph>::vertex_descriptor vertex_t;
        typedef typename graph_traits<Graph>::edge_descriptor edge_t;

        vector<edge_t> edges(num_edges(g));
Tiago Peixoto's avatar
Tiago Peixoto committed
281
        vector<size_t> edges_sequence;
282
283
        typename graph_traits<Graph>::edge_iterator e, e_end;
        for (tie(e, e_end) = boost::edges(g); e != e_end; ++e)
284
        {
285
286
287
            if (edge_index[*e] >= edges.size())
                edges.resize(edge_index[*e] + 1);
            edges[edge_index[*e]] = *e;
Tiago Peixoto's avatar
Tiago Peixoto committed
288
            edges_sequence.push_back(edge_index[*e]);
289
290
        }

291
292
        RewireStrategy<Graph, EdgeIndexMap, CorrProb>
            rewire(g, edge_index, edges, corr_prob, rng);
293

Tiago Peixoto's avatar
Tiago Peixoto committed
294
295
296
        typedef random_permutation_iterator<typename vector<size_t>::iterator,
                                            rng_t>
            random_edge_iter;
297

Tiago Peixoto's avatar
Tiago Peixoto committed
298
299
300
        random_edge_iter
            ei(edges_sequence.begin(), edges_sequence.end(), rng),
            ei_end(edges_sequence.end(), edges_sequence.end(), rng);
301

302
303
304
        if (verbose)
            cout << "rewiring edges: ";

Tiago Peixoto's avatar
Tiago Peixoto committed
305
        stringstream str;
306
307
        size_t count = 0;
        // for each edge rewire its source or target
Tiago Peixoto's avatar
Tiago Peixoto committed
308
        for (; ei != ei_end; ++ei)
309
        {
310
311
312
            if (verbose)
                print_progress(count++, edges_sequence.size(), str);
            rewire(*ei, self_loops, parallel_edges);
313
        }
314
315
        if (verbose)
            cout << endl;
316
    }
317
318
};

Tiago Peixoto's avatar
Tiago Peixoto committed
319

320
321
// this will rewire the edges so that the resulting graph will be entirely
// random (i.e. Erdos-Renyi)
322
template <class Graph, class EdgeIndexMap, class CorrProb>
323
324
325
326
327
328
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
329
    ErdosRewireStrategy(Graph& g, EdgeIndexMap edge_index,
330
                        vector<edge_t>& edges, CorrProb corr, rng_t& rng)
Tiago Peixoto's avatar
Tiago Peixoto committed
331
332
        : _g(g), _edge_index(edge_index), _edges(edges),
          _vertices(HardNumVertices()(g)), _rng(rng)
333
334
335
336
337
338
339
    {
        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;
    }

340
    void operator()(size_t e, bool self_loops, bool parallel_edges)
341
342
343
344
345
346
347
    {
        //try randomly drawn pairs of vertices until one satisfies all the
        //consistency checks
        typedef random_permutation_iterator
            <typename graph_traits<Graph>::vertex_iterator, rng_t>
            random_vertex_iter;

348
        tr1::uniform_int<size_t> sample(0, _vertices.size()-1);
349
350
351
352
353
354
355
356
        typename graph_traits<Graph>::vertex_descriptor s, t;
        while (true)
        {
            s = sample(_rng);
            t = sample(_rng);

            if(s == t && !self_loops) // reject self-loops if not allowed
                continue;
357
            if (!parallel_edges && is_adjacent(s, t, _g))
358
359
360
361
                continue;  // reject parallel edges if not allowed
            break;
        }
        edge_t ne = add_edge(s, t, _g).first;
362
363
364
        _edge_index[ne] = e;
        remove_edge(_edges[e], _g);
        _edges[e] = ne;
365
366
367
368
369
    }

private:
    Graph& _g;
    EdgeIndexMap _edge_index;
Tiago Peixoto's avatar
Tiago Peixoto committed
370
    vector<edge_t>& _edges;
371
372
373
374
375
    vector<typename graph_traits<Graph>::vertex_descriptor> _vertices;
    rng_t& _rng;
};


376
377
378
// 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.
379
380
template <class Graph, class EdgeIndexMap, class RewireStrategy>
class RewireStrategyBase
381
382
383
{
public:
    typedef typename graph_traits<Graph>::edge_descriptor edge_t;
384
385
    typedef typename graph_traits<Graph>::edge_descriptor vertex_t;

386
387
    typedef typename EdgeIndexMap::value_type index_t;

Tiago Peixoto's avatar
Tiago Peixoto committed
388
389
    RewireStrategyBase(Graph& g, EdgeIndexMap edge_index, vector<edge_t>& edges,
                       rng_t& rng)
390
        : _g(g), _edge_index(edge_index), _edges(edges),  _rng(rng) {}
391

392
    void operator()(size_t e, bool self_loops, bool parallel_edges)
393
    {
Tiago Peixoto's avatar
Tiago Peixoto committed
394
395
        RewireStrategy& self = *static_cast<RewireStrategy*>(this);
        tr1::bernoulli_distribution coin;
396
397
398

        //try randomly drawn pairs of edges until one satisfies all the
        //consistency checks
399
        pair<size_t, bool> es, et;
400

Tiago Peixoto's avatar
Tiago Peixoto committed
401
        while (true)
402
        {
Tiago Peixoto's avatar
Tiago Peixoto committed
403
404
            es.first = et.first = e;
            es.second = et.second = false;
405

Tiago Peixoto's avatar
Tiago Peixoto committed
406
            if (coin(_rng))
407
            {
Tiago Peixoto's avatar
Tiago Peixoto committed
408
409
                // rewire source
                es = self.get_source_edge(e);
410

Tiago Peixoto's avatar
Tiago Peixoto committed
411
412
                if(!self_loops) // reject self-loops if not allowed
                {
413
414
                    if((source(_edges[e], _g) == target(es, _edges, _g)) ||
                       (target(_edges[e], _g) == source(es, _edges, _g)))
Tiago Peixoto's avatar
Tiago Peixoto committed
415
416
                        continue;
                }
417
418
419

                // reject parallel edges if not allowed
                if (!parallel_edges && (es.first != e))
Tiago Peixoto's avatar
Tiago Peixoto committed
420
                {
421
                    if (swap_edge::parallel_check_source(e, es, _edges, _g))
Tiago Peixoto's avatar
Tiago Peixoto committed
422
423
                        continue;
                }
424
425
426
427
428
429
430
431
432
433
434

                if (e != es.first)
                {
                    self.update_edge(e, false, true);
                    self.update_edge(es.first, false, false);

                    swap_edge::swap_source(e, es, _edges, _edge_index, _g);

                    self.update_edge(e, true, true);
                    self.update_edge(es.first, true, false);
                }
Tiago Peixoto's avatar
Tiago Peixoto committed
435
436
            }
            else
437
            {
Tiago Peixoto's avatar
Tiago Peixoto committed
438
439
                // rewire target
                et = self.get_target_edge(e);
440

441
442
                if (!self_loops) // reject self-loops if not allowed
                {
443
444
                    if((source(_edges[e], _g) == target(et, _edges, _g)) ||
                       (target(_edges[e], _g) == source(et, _edges, _g)))
445
446
                        continue;
                }
447
448
449

                // reject parallel edges if not allowed
                if (!parallel_edges && (et.first != e))
450
                {
451
                    if (swap_edge::parallel_check_target(e, et, _edges, _g))
452
453
                        continue;
                }
Tiago Peixoto's avatar
Tiago Peixoto committed
454

455
456
457
458
                if (e != et.first)
                {
                    self.update_edge(e, false, true);
                    self.update_edge(et.first, false, false);
Tiago Peixoto's avatar
Tiago Peixoto committed
459

460
                    swap_edge::swap_target(e, et, _edges, _edge_index, _g);
Tiago Peixoto's avatar
Tiago Peixoto committed
461

462
463
464
465
466
467
                    self.update_edge(e, true, true);
                    self.update_edge(et.first, true, false);
                }
            }
            break;
        }
468
469
    }

Tiago Peixoto's avatar
Tiago Peixoto committed
470
protected:
471
472
    Graph& _g;
    EdgeIndexMap _edge_index;
Tiago Peixoto's avatar
Tiago Peixoto committed
473
    vector<edge_t>& _edges;
474
    rng_t& _rng;
475
476
};

477
478
// this will rewire the edges so that the combined (in, out) degree distribution
// will be the same, but all the rest is random
479
template <class Graph, class EdgeIndexMap, class CorrProb>
480
481
class RandomRewireStrategy:
    public RewireStrategyBase<Graph, EdgeIndexMap,
482
483
                              RandomRewireStrategy<Graph, EdgeIndexMap,
                                                   CorrProb> >
484
485
{
public:
486
    typedef RewireStrategyBase<Graph, EdgeIndexMap,
487
488
                               RandomRewireStrategy<Graph, EdgeIndexMap,
                                                    CorrProb> >
489
490
491
492
493
        base_t;

    typedef Graph graph_t;
    typedef EdgeIndexMap edge_index_t;

494
495
    typedef typename graph_traits<Graph>::vertex_descriptor vertex_t;
    typedef typename graph_traits<Graph>::edge_descriptor edge_t;
496
    typedef typename EdgeIndexMap::value_type index_t;
497

Tiago Peixoto's avatar
Tiago Peixoto committed
498
499
500
501
502
503
504
505
506
507
    struct hash_index {};
    struct random_index {};

    typedef multi_index_container<
        pair<index_t,bool>,
        indexed_by<
            hashed_unique<tag<hash_index>, identity<pair<index_t,bool> > >,
            random_access<tag<random_index> > > >
        edge_set_t;

508
    RandomRewireStrategy(Graph& g, EdgeIndexMap edge_index,
509
                         vector<edge_t>& edges, CorrProb corr, rng_t& rng)
Tiago Peixoto's avatar
Tiago Peixoto committed
510
        : base_t(g, edge_index, edges, rng), _g(g)
511
    {
Tiago Peixoto's avatar
Tiago Peixoto committed
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
        size_t E = HardNumEdges()(_g);

        typename graph_traits<Graph>::vertex_iterator v_i, v_i_end;
        for (tie(v_i, v_i_end) = vertices(g); v_i != v_i_end; ++v_i)
        {
            if (out_degree(*v_i, _g) > E/3)
            {
                _hub_non_adjacent[*v_i] = edge_set_t();
                if(!is_directed::apply<Graph>::type::value)
                    _hub_non_incident[*v_i] = edge_set_t();
            }

            if (in_degreeS()(*v_i, _g) > E/3)
                _hub_non_incident[*v_i] = edge_set_t();
        }
527

528
        typename graph_traits<Graph>::edge_iterator e_i, e_i_end;
Tiago Peixoto's avatar
Tiago Peixoto committed
529
        for (tie(e_i, e_i_end) = boost::edges(g); e_i != e_i_end; ++e_i)
530
        {
Tiago Peixoto's avatar
Tiago Peixoto committed
531
532
533
            // For undirected graphs, there is no difference between source and
            // target, and each edge will appear _twice_ in the list below,
            // once for each different ordering of source and target.
534

Tiago Peixoto's avatar
Tiago Peixoto committed
535
            _all_edges.push_back(make_pair(edge_index[*e_i], false));
536
            if (!is_directed::apply<Graph>::type::value)
Tiago Peixoto's avatar
Tiago Peixoto committed
537
                _all_edges.push_back(make_pair(edge_index[*e_i], true));
538
            update_edge(edge_index[*e_i], true);
Tiago Peixoto's avatar
Tiago Peixoto committed
539
540
541
        }
    }

542
    pair<size_t,bool> get_edge(size_t e, bool src)
Tiago Peixoto's avatar
Tiago Peixoto committed
543
544
545
546
    {
        tr1::unordered_map<size_t, edge_set_t>& hub =
            (src || !is_directed::apply<Graph>::type::value) ?
            _hub_non_adjacent : _hub_non_incident;
547

Tiago Peixoto's avatar
Tiago Peixoto committed
548
        typeof(hub.begin()) iter =
549
550
            hub.find(src ? target(base_t::_edges[e], _g) :
                     source(base_t::_edges[e], _g));
Tiago Peixoto's avatar
Tiago Peixoto committed
551
552
553
554
555
556
557
558
559
        if (iter != hub.end())
        {
            edge_set_t& eset = iter->second;

            tr1::uniform_int<> sample(0, eset.size());
            size_t i = sample(base_t::_rng);
            if (i == eset.size()) // no rewire option
                return make_pair(e, false);
            else
560
                return *(eset.get<random_index>().begin() + i);
561
562
        }

Tiago Peixoto's avatar
Tiago Peixoto committed
563
        tr1::uniform_int<> sample(0, _all_edges.size()-1);
564
        return _all_edges[sample(base_t::_rng)];
565
    }
566

Tiago Peixoto's avatar
Tiago Peixoto committed
567

568
    pair<size_t,bool> get_source_edge(size_t e)
569
    {
Tiago Peixoto's avatar
Tiago Peixoto committed
570
571
572
        return get_edge(e, true);
    }

573
    pair<size_t,bool> get_target_edge(size_t e)
Tiago Peixoto's avatar
Tiago Peixoto committed
574
575
576
577
    {
        return get_edge(e, false);
    }

578
    void update_edge_list(index_t v, const pair<size_t, bool>& e,
Tiago Peixoto's avatar
Tiago Peixoto committed
579
580
581
                          edge_set_t& edge_set, bool src, bool insert)
    {
        pair<index_t,bool> ep;
582
583
        if ((src && source(e, base_t::_edges, _g) != v) ||
            (!src && target(e, base_t::_edges, _g) != v))
Tiago Peixoto's avatar
Tiago Peixoto committed
584
585
        {
            if (insert)
586
                edge_set.get<hash_index>().insert(e);
Tiago Peixoto's avatar
Tiago Peixoto committed
587
            else
588
                edge_set.get<hash_index>().erase(e);
Tiago Peixoto's avatar
Tiago Peixoto committed
589
590
591
        }
    }

592
    void update_edge(size_t e, bool insert, bool final = false)
Tiago Peixoto's avatar
Tiago Peixoto committed
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
    {
        for (typeof(_hub_non_adjacent.begin()) iter = _hub_non_adjacent.begin();
             iter != _hub_non_adjacent.end(); ++iter)
        {
            update_edge_list(iter->first, make_pair(e, false), iter->second,
                             true, insert);
            if (!is_directed::apply<Graph>::type::value)
                update_edge_list(iter->first, make_pair(e, true), iter->second,
                                 true, insert);
        }

        for (typeof(_hub_non_incident.begin()) iter = _hub_non_incident.begin();
             iter != _hub_non_incident.end(); ++iter)
        {
            update_edge_list(iter->first, make_pair(e, false), iter->second,
                             false, insert);
            if (!is_directed::apply<Graph>::type::value)
                update_edge_list(iter->first, make_pair(e, true), iter->second,
                                 false, insert);
        }
613
    }
614

615
private:
Tiago Peixoto's avatar
Tiago Peixoto committed
616
617
    Graph& _g;
    EdgeIndexMap _edge_index;
618
    vector<pair<index_t,bool> > _all_edges;
Tiago Peixoto's avatar
Tiago Peixoto committed
619
620
621

    tr1::unordered_map<size_t, edge_set_t>  _hub_non_incident;
    tr1::unordered_map<size_t, edge_set_t>  _hub_non_adjacent;
622
};
623

624
625
626

// 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
627
template <class Graph, class EdgeIndexMap, class CorrProb>
628
629
class CorrelatedRewireStrategy:
    public RewireStrategyBase<Graph, EdgeIndexMap,
630
631
                              CorrelatedRewireStrategy<Graph, EdgeIndexMap,
                                                       CorrProb> >
632
633
{
public:
634
    typedef RewireStrategyBase<Graph, EdgeIndexMap,
635
636
                               CorrelatedRewireStrategy<Graph, EdgeIndexMap,
                                                        CorrProb> >
637
638
639
640
641
        base_t;

    typedef Graph graph_t;
    typedef EdgeIndexMap edge_index_t;

642
643
644
645
    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;

646
    CorrelatedRewireStrategy (Graph& g, EdgeIndexMap edge_index,
647
                              vector<edge_t>& edges, CorrProb corr, rng_t& rng)
Tiago Peixoto's avatar
Tiago Peixoto committed
648
        : base_t(g, edge_index, edges, rng), _g(g)
649
    {
650
        typename graph_traits<Graph>::edge_iterator e_i, e_i_end;
Tiago Peixoto's avatar
Tiago Peixoto committed
651
        for (tie(e_i, e_i_end) = boost::edges(g); e_i != e_i_end; ++e_i)
652
        {
653
            // For undirected graphs, there is no difference between source and
Tiago Peixoto's avatar
Tiago Peixoto committed
654
            // target, and each edge will appear _twice_ in the lists below,
655
            // once for each different ordering of source and target.
656

657
658
659
660
            _edges_by_source
                [make_pair(in_degreeS()(source(*e_i, _g), _g),
                           out_degree(source(*e_i, _g), _g))]
                .push_back(make_pair(edge_index[*e_i], false));
661

662
663
664
665
666
667
668
669
            _edges_by_target
                [make_pair(in_degreeS()(target(*e_i, _g), _g),
                           out_degree(target(*e_i, _g), _g))]
                .push_back(make_pair(edge_index[*e_i], false));

            if (!is_directed::apply<Graph>::type::value)
            {
                _edges_by_source
670
671
                    [make_pair(in_degreeS()(target(*e_i, _g), _g),
                               out_degree(target(*e_i, _g), _g))]
672
673
674
675
676
                    .push_back(make_pair(edge_index[*e_i], true));
                _edges_by_target
                    [make_pair(in_degreeS()(source(*e_i, _g), _g),
                               out_degree(source(*e_i, _g), _g))]
                    .push_back(make_pair(edge_index[*e_i], true));
677
678
679
            }
        }
    }
680

681
    pair<size_t,bool> get_edge(size_t e, bool src)
682
    {
Tiago Peixoto's avatar
Tiago Peixoto committed
683
        pair<size_t, size_t> deg =
684
685
686
687
            make_pair(in_degreeS()(src ? source(base_t::_edges[e], _g)
                                   : target(base_t::_edges[e], _g), _g),
                      out_degree(src ? source(base_t::_edges[e], _g) :
                                 target(base_t::_edges[e], _g), _g));
Tiago Peixoto's avatar
Tiago Peixoto committed
688
689
690
        edges_by_end_deg_t& edges = src ? _edges_by_source : _edges_by_target;
        typename edges_by_end_deg_t::mapped_type& elist = edges[deg];
        tr1::uniform_int<> sample(0, elist.size()-1);
691
        return elist[sample(base_t::_rng)];
Tiago Peixoto's avatar
Tiago Peixoto committed
692
    }
693

694
    pair<size_t,bool> get_source_edge(size_t e)
Tiago Peixoto's avatar
Tiago Peixoto committed
695
696
697
    {
        return get_edge(e, true);
    }
698

699
    pair<size_t,bool> get_target_edge(size_t e)
Tiago Peixoto's avatar
Tiago Peixoto committed
700
701
    {
        return get_edge(e, false);
702
703
    }

704
    void update_edge(size_t e, bool insert, bool) {}
Tiago Peixoto's avatar
Tiago Peixoto committed
705

706
private:
707
708
    typedef tr1::unordered_map<pair<size_t, size_t>,
                               vector<pair<index_t, bool> >,
709
710
                               hash<pair<size_t, size_t> > >
        edges_by_end_deg_t;
711
    edges_by_end_deg_t _edges_by_source, _edges_by_target;
712
713
714

protected:
    const Graph& _g;
715
716
};

717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
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
861
862
863
864
865
866
867
868
869
template <class Graph, class EdgeIndexMap, class CorrProb>
class ProbabilisticRewireStrategy:
    public RewireStrategyBase<Graph, EdgeIndexMap,
                              ProbabilisticRewireStrategy<Graph, EdgeIndexMap,
                                                          CorrProb> >
{
public:
    typedef RewireStrategyBase<Graph, EdgeIndexMap,
                                ProbabilisticRewireStrategy<Graph, EdgeIndexMap,
                                                            CorrProb> >
        base_t;

    typedef Graph graph_t;
    typedef EdgeIndexMap edge_index_t;

    typedef pair<size_t,size_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;

    ProbabilisticRewireStrategy(Graph& g, EdgeIndexMap edge_index,
                                vector<edge_t>& edges,
                                CorrProb corr_prob, rng_t& rng)
        : base_t(g, edge_index, edges, rng), _g(g)
    {
        typename graph_traits<Graph>::vertex_iterator v_i, v_i_end;
        for (tie(v_i, v_i_end) = boost::vertices(g); v_i != v_i_end; ++v_i)
            _deg_count[get_deg(*v_i, g)]++;


        typename graph_traits<Graph>::edge_iterator e_i, e_i_end;
        set<pair<size_t, size_t> > deg_set;
        for (tie(e_i, e_i_end) = boost::edges(g); e_i != e_i_end; ++e_i)
        {
            update_edge(_edge_index[*e_i], true);
            deg_set.insert(get_deg(source(*e_i, g), g));
            deg_set.insert(get_deg(target(*e_i, g), g));
        }

        for (typeof(deg_set.begin()) s_iter = deg_set.begin();
             s_iter != deg_set.end(); ++s_iter)
        {
             _sample_source_deg[*s_iter] = Sampler<deg_t>(true);
             _sample_target_deg[*s_iter] = Sampler<deg_t>(true);
        }

        for (typeof(deg_set.begin()) s_iter = deg_set.begin();
             s_iter != deg_set.end(); ++s_iter)
            for (typeof(deg_set.begin()) t_iter = deg_set.begin();
                 t_iter != deg_set.end(); ++t_iter)
        {
            double p = corr_prob(*s_iter, *t_iter);
            if (!is_directed::apply<Graph>::type::value ||
                (s_iter->second > 0 && t_iter->first > 0))
            {
                _sample_source_deg[*t_iter].Insert(*s_iter, p);
                _sample_target_deg[*s_iter].Insert(*t_iter, p);
            }
        }
    }

    deg_t get_deg(vertex_t v, Graph& g)
    {
        return make_pair(in_degreeS()(v, g), out_degree(v,g));
    }

    pair<size_t,bool> get_source_edge(size_t e)
    {
        deg_t t_deg = get_deg(target(base_t::_edges[e],_g),_g), deg;
        while (true)
        {
            deg = _sample_source_deg[t_deg](base_t::_rng);
            if (_sample_edge_by_source_deg[deg].Empty())
                _sample_source_deg[t_deg].Remove(deg);
            else
                break;
        }
        pair<size_t, bool> ep = _sample_edge_by_source_deg[deg](base_t::_rng);
        //assert(get_deg(source(ep, base_t::_edges, _g), _g) == deg);
        return ep;
    }

    pair<size_t, bool> get_target_edge(size_t e)
    {
        deg_t s_deg = get_deg(source(base_t::_edges[e],_g),_g), deg;
        while (true)
        {
            deg = _sample_target_deg[s_deg](base_t::_rng);
            if (_sample_edge_by_target_deg[deg].Empty())
                _sample_target_deg[s_deg].Remove(deg);
            else
                break;
        }
        pair<size_t, bool> ep = _sample_edge_by_target_deg[deg](base_t::_rng);
        //assert(get_deg(target(ep, base_t::_edges, _g), _g) == deg);
        return ep;
    }

    void update_edge(size_t e, bool insert, bool final = false)
    {
        deg_t s_deg = get_deg(source(base_t::_edges[e], _g), _g),
            t_deg = get_deg(target(base_t::_edges[e], _g), _g);

        pair<size_t, bool> ep = make_pair(e, false);

        if (insert && !final)
        {
            _sample_edge_by_source_deg[s_deg].Insert(ep);
            _sample_edge_by_target_deg[t_deg].Insert(ep);
            if (!is_directed::apply<Graph>::type::value)
            {
                ep.second = true;
                _sample_edge_by_source_deg[t_deg].Insert(ep);
                _sample_edge_by_target_deg[s_deg].Insert(ep);
            }
        }

        if (!insert)
        {
            _sample_edge_by_source_deg[s_deg].Remove(ep);
            _sample_edge_by_target_deg[t_deg].Remove(ep);
            if (!is_directed::apply<Graph>::type::value)
            {
                ep.second = true;
                _sample_edge_by_source_deg[t_deg].Remove(ep);
                _sample_edge_by_target_deg[s_deg].Remove(ep);
            }
        }
    }

private:
    Graph& _g;
    EdgeIndexMap _edge_index;

    typedef tr1::unordered_map<pair<deg_t, deg_t>, double,
                               hash<pair<deg_t, deg_t> > > corr_prob_t;
    typedef tr1::unordered_map<deg_t, Sampler<deg_t>, hash<deg_t> >
        deg_sampler_t;

    typedef tr1::unordered_map<deg_t, Sampler<pair<size_t,bool> >,
                               hash<deg_t> >
        edge_sampler_t;

    edge_sampler_t _sample_edge_by_source_deg;
    edge_sampler_t _sample_edge_by_target_deg;

    deg_sampler_t _sample_source_deg;
    deg_sampler_t _sample_target_deg;

    tr1::unordered_map<deg_t, size_t, hash<deg_t> > _deg_count;
};

870
871
872
} // graph_tool namespace

#endif // GRAPH_REWIRING_HH