graph_rewiring.hh 29.5 KB
Newer Older
1
2
// graph-tool -- a general graph modification and manipulation thingy
//
3
// Copyright (C) 2007-2010 Tiago de Paula Peixoto <tiago@forked.de>
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
//
// 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

#include <tr1/unordered_set>
22
#include <tr1/random>
23
#include <boost/functional/hash.hpp>
Tiago Peixoto's avatar
Tiago Peixoto committed
24
25
26
27
28
29
#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>
30
31
32

#include "graph.hh"
#include "graph_filtering.hh"
33
#include "graph_util.hh"
34
#include "sampler.hh"
35
36
37
38
39

namespace graph_tool
{
using namespace std;
using namespace boost;
Tiago Peixoto's avatar
Tiago Peixoto committed
40
using namespace boost::multi_index;
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56

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

57
58
template <class Graph>
typename graph_traits<Graph>::vertex_descriptor
59
60
source(const pair<size_t, bool>& e,
       const vector<typename graph_traits<Graph>::edge_descriptor>& edges,
61
62
63
       const Graph& g)
{
    if (e.second)
64
        return target(edges[e.first], g);
65
    else
66
        return source(edges[e.first], g);
67
68
69
70
}

template <class Graph>
typename graph_traits<Graph>::vertex_descriptor
71
72
target(const pair<size_t, bool>& e,
       const vector<typename graph_traits<Graph>::edge_descriptor>& edges,
73
74
75
       const Graph& g)
{
    if (e.second)
76
        return source(edges[e.first], g);
77
    else
78
        return target(edges[e.first], g);
79
80
}

Tiago Peixoto's avatar
Tiago Peixoto committed
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
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
// 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;
};

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

        typename graph_traits<Graph>::vertex_descriptor
155
156
157
158
            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
159

160
        if (is_adjacent(ns,  t, g))
Tiago Peixoto's avatar
Tiago Peixoto committed
161
            return true; // e would clash with an existing (new) edge
162
        if (is_adjacent(s, se_t, g))
Tiago Peixoto's avatar
Tiago Peixoto committed
163
164
165
166
            return true; // se would clash with an existing (new) edge
        return false; // the coast is clear - hooray!
    }

167
    template <class Graph>
Tiago Peixoto's avatar
Tiago Peixoto committed
168
    static bool
169
170
171
    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
172
173
174
175
176
177
178
    {
        // 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)
        //
179
        // no parallel edges are introduced.
Tiago Peixoto's avatar
Tiago Peixoto committed
180
181

        typename graph_traits<Graph>::vertex_descriptor
182
183
184
185
            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
186

187
        if (is_adjacent(s,  nt, g))
188
            return true; // e would clash with an existing (new) edge
189
        if (is_adjacent(te_s, t, g))
190
191
192
193
            return true; // te would clash with an existing (new) edge
        return false; // the coast is clear - hooray!
    }

194
195
196
197
    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)
198
    {
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
        // 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;
222
223
    }

224
225
226
227
228
    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)
229
    {
230
231
        // swap the source of the edge 'e' with the source of edge 'se', as
        // such:
232
        //
233
234
        //  (s)    -e--> (t)          (s)    -e--> (nt)
        //  (te_s) -te-> (nt)   =>    (te_s) -te-> (t)
235

236
237
        if (e == te.first)
            return;
238

239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
        // 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;
254
    }
255

256
257
};

Tiago Peixoto's avatar
Tiago Peixoto committed
258
259
260
// used for verbose display
void print_progress(size_t current, size_t total, stringstream& str);

261
// main rewire loop
262
263
template <template <class Graph, class EdgeIndexMap, class CorrProb>
          class RewireStrategy>
264
265
struct graph_rewire
{
266
267
268
269
    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
270
271
272
273
274
    {
        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
275
        vector<size_t> edges_sequence;
276
277
        typename graph_traits<Graph>::edge_iterator e, e_end;
        for (tie(e, e_end) = boost::edges(g); e != e_end; ++e)
278
        {
279
280
281
            if (edge_index[*e] >= edges.size())
                edges.resize(edge_index[*e] + 1);
            edges[edge_index[*e]] = *e;
Tiago Peixoto's avatar
Tiago Peixoto committed
282
            edges_sequence.push_back(edge_index[*e]);
283
284
        }

285
286
        RewireStrategy<Graph, EdgeIndexMap, CorrProb>
            rewire(g, edge_index, edges, corr_prob, rng);
287

Tiago Peixoto's avatar
Tiago Peixoto committed
288
289
290
        typedef random_permutation_iterator<typename vector<size_t>::iterator,
                                            rng_t>
            random_edge_iter;
291

Tiago Peixoto's avatar
Tiago Peixoto committed
292
293
294
        random_edge_iter
            ei(edges_sequence.begin(), edges_sequence.end(), rng),
            ei_end(edges_sequence.end(), edges_sequence.end(), rng);
295

296
297
298
        if (verbose)
            cout << "rewiring edges: ";

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

Tiago Peixoto's avatar
Tiago Peixoto committed
313

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

334
    void operator()(size_t e, bool self_loops, bool parallel_edges)
335
336
337
338
339
340
341
    {
        //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;

342
        tr1::uniform_int<size_t> sample(0, _vertices.size()-1);
343
344
345
346
347
348
349
350
        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;
351
            if (!parallel_edges && is_adjacent(s, t, _g))
352
353
354
355
                continue;  // reject parallel edges if not allowed
            break;
        }
        edge_t ne = add_edge(s, t, _g).first;
356
357
358
        _edge_index[ne] = e;
        remove_edge(_edges[e], _g);
        _edges[e] = ne;
359
360
361
362
363
    }

private:
    Graph& _g;
    EdgeIndexMap _edge_index;
Tiago Peixoto's avatar
Tiago Peixoto committed
364
    vector<edge_t>& _edges;
365
366
367
368
369
    vector<typename graph_traits<Graph>::vertex_descriptor> _vertices;
    rng_t& _rng;
};


370
371
372
// 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.
373
374
template <class Graph, class EdgeIndexMap, class RewireStrategy>
class RewireStrategyBase
375
376
377
{
public:
    typedef typename graph_traits<Graph>::edge_descriptor edge_t;
378
379
    typedef typename graph_traits<Graph>::edge_descriptor vertex_t;

380
381
    typedef typename EdgeIndexMap::value_type index_t;

Tiago Peixoto's avatar
Tiago Peixoto committed
382
383
    RewireStrategyBase(Graph& g, EdgeIndexMap edge_index, vector<edge_t>& edges,
                       rng_t& rng)
384
        : _g(g), _edge_index(edge_index), _edges(edges),  _rng(rng) {}
385

386
    void operator()(size_t e, bool self_loops, bool parallel_edges)
387
    {
Tiago Peixoto's avatar
Tiago Peixoto committed
388
389
        RewireStrategy& self = *static_cast<RewireStrategy*>(this);
        tr1::bernoulli_distribution coin;
390
391
392

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

Tiago Peixoto's avatar
Tiago Peixoto committed
395
        while (true)
396
        {
Tiago Peixoto's avatar
Tiago Peixoto committed
397
398
            es.first = et.first = e;
            es.second = et.second = false;
399

Tiago Peixoto's avatar
Tiago Peixoto committed
400
            if (coin(_rng))
401
            {
Tiago Peixoto's avatar
Tiago Peixoto committed
402
403
                // rewire source
                es = self.get_source_edge(e);
404

Tiago Peixoto's avatar
Tiago Peixoto committed
405
406
                if(!self_loops) // reject self-loops if not allowed
                {
407
408
                    if((source(_edges[e], _g) == target(es, _edges, _g)) ||
                       (target(_edges[e], _g) == source(es, _edges, _g)))
Tiago Peixoto's avatar
Tiago Peixoto committed
409
410
                        continue;
                }
411
412
413

                // reject parallel edges if not allowed
                if (!parallel_edges && (es.first != e))
Tiago Peixoto's avatar
Tiago Peixoto committed
414
                {
415
                    if (swap_edge::parallel_check_source(e, es, _edges, _g))
Tiago Peixoto's avatar
Tiago Peixoto committed
416
417
                        continue;
                }
418
419
420
421
422
423
424
425
426
427
428

                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
429
430
            }
            else
431
            {
Tiago Peixoto's avatar
Tiago Peixoto committed
432
433
                // rewire target
                et = self.get_target_edge(e);
434

435
436
                if (!self_loops) // reject self-loops if not allowed
                {
437
438
                    if((source(_edges[e], _g) == target(et, _edges, _g)) ||
                       (target(_edges[e], _g) == source(et, _edges, _g)))
439
440
                        continue;
                }
441
442
443

                // reject parallel edges if not allowed
                if (!parallel_edges && (et.first != e))
444
                {
445
                    if (swap_edge::parallel_check_target(e, et, _edges, _g))
446
447
                        continue;
                }
Tiago Peixoto's avatar
Tiago Peixoto committed
448

449
450
451
452
                if (e != et.first)
                {
                    self.update_edge(e, false, true);
                    self.update_edge(et.first, false, false);
Tiago Peixoto's avatar
Tiago Peixoto committed
453

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

456
457
458
459
460
461
                    self.update_edge(e, true, true);
                    self.update_edge(et.first, true, false);
                }
            }
            break;
        }
462
463
    }

Tiago Peixoto's avatar
Tiago Peixoto committed
464
protected:
465
466
    Graph& _g;
    EdgeIndexMap _edge_index;
Tiago Peixoto's avatar
Tiago Peixoto committed
467
    vector<edge_t>& _edges;
468
    rng_t& _rng;
469
470
};

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

    typedef Graph graph_t;
    typedef EdgeIndexMap edge_index_t;

488
489
    typedef typename graph_traits<Graph>::vertex_descriptor vertex_t;
    typedef typename graph_traits<Graph>::edge_descriptor edge_t;
490
    typedef typename EdgeIndexMap::value_type index_t;
491

Tiago Peixoto's avatar
Tiago Peixoto committed
492
493
494
495
496
497
498
499
500
501
    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;

502
    RandomRewireStrategy(Graph& g, EdgeIndexMap edge_index,
503
                         vector<edge_t>& edges, CorrProb corr, rng_t& rng)
Tiago Peixoto's avatar
Tiago Peixoto committed
504
        : base_t(g, edge_index, edges, rng), _g(g)
505
    {
Tiago Peixoto's avatar
Tiago Peixoto committed
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
        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();
        }
521

522
        typename graph_traits<Graph>::edge_iterator e_i, e_i_end;
Tiago Peixoto's avatar
Tiago Peixoto committed
523
        for (tie(e_i, e_i_end) = boost::edges(g); e_i != e_i_end; ++e_i)
524
        {
Tiago Peixoto's avatar
Tiago Peixoto committed
525
526
527
            // 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.
528

Tiago Peixoto's avatar
Tiago Peixoto committed
529
            _all_edges.push_back(make_pair(edge_index[*e_i], false));
530
            if (!is_directed::apply<Graph>::type::value)
Tiago Peixoto's avatar
Tiago Peixoto committed
531
                _all_edges.push_back(make_pair(edge_index[*e_i], true));
532
            update_edge(edge_index[*e_i], true);
Tiago Peixoto's avatar
Tiago Peixoto committed
533
534
535
        }
    }

536
    pair<size_t,bool> get_edge(size_t e, bool src)
Tiago Peixoto's avatar
Tiago Peixoto committed
537
538
539
540
    {
        tr1::unordered_map<size_t, edge_set_t>& hub =
            (src || !is_directed::apply<Graph>::type::value) ?
            _hub_non_adjacent : _hub_non_incident;
541

Tiago Peixoto's avatar
Tiago Peixoto committed
542
        typeof(hub.begin()) iter =
543
544
            hub.find(src ? target(base_t::_edges[e], _g) :
                     source(base_t::_edges[e], _g));
Tiago Peixoto's avatar
Tiago Peixoto committed
545
546
547
548
549
550
551
552
553
        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
554
                return *(eset.get<random_index>().begin() + i);
555
556
        }

Tiago Peixoto's avatar
Tiago Peixoto committed
557
        tr1::uniform_int<> sample(0, _all_edges.size()-1);
558
        return _all_edges[sample(base_t::_rng)];
559
    }
560

Tiago Peixoto's avatar
Tiago Peixoto committed
561

562
    pair<size_t,bool> get_source_edge(size_t e)
563
    {
Tiago Peixoto's avatar
Tiago Peixoto committed
564
565
566
        return get_edge(e, true);
    }

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

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

586
    void update_edge(size_t e, bool insert, bool final = false)
Tiago Peixoto's avatar
Tiago Peixoto committed
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
    {
        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);
        }
607
    }
608

609
private:
Tiago Peixoto's avatar
Tiago Peixoto committed
610
611
    Graph& _g;
    EdgeIndexMap _edge_index;
612
    vector<pair<index_t,bool> > _all_edges;
Tiago Peixoto's avatar
Tiago Peixoto committed
613
614
615

    tr1::unordered_map<size_t, edge_set_t>  _hub_non_incident;
    tr1::unordered_map<size_t, edge_set_t>  _hub_non_adjacent;
616
};
617

618
619
620

// 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
621
template <class Graph, class EdgeIndexMap, class CorrProb>
622
623
class CorrelatedRewireStrategy:
    public RewireStrategyBase<Graph, EdgeIndexMap,
624
625
                              CorrelatedRewireStrategy<Graph, EdgeIndexMap,
                                                       CorrProb> >
626
627
{
public:
628
    typedef RewireStrategyBase<Graph, EdgeIndexMap,
629
630
                               CorrelatedRewireStrategy<Graph, EdgeIndexMap,
                                                        CorrProb> >
631
632
633
634
635
        base_t;

    typedef Graph graph_t;
    typedef EdgeIndexMap edge_index_t;

636
637
638
639
    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;

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

651
652
653
654
            _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));
655

656
657
658
659
660
661
662
663
            _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
664
665
                    [make_pair(in_degreeS()(target(*e_i, _g), _g),
                               out_degree(target(*e_i, _g), _g))]
666
667
668
669
670
                    .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));
671
672
673
            }
        }
    }
674

675
    pair<size_t,bool> get_edge(size_t e, bool src)
676
    {
Tiago Peixoto's avatar
Tiago Peixoto committed
677
        pair<size_t, size_t> deg =
678
679
680
681
            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
682
683
684
        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);
685
        return elist[sample(base_t::_rng)];
Tiago Peixoto's avatar
Tiago Peixoto committed
686
    }
687

688
    pair<size_t,bool> get_source_edge(size_t e)
Tiago Peixoto's avatar
Tiago Peixoto committed
689
690
691
    {
        return get_edge(e, true);
    }
692

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

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

700
private:
701
702
    typedef tr1::unordered_map<pair<size_t, size_t>,
                               vector<pair<index_t, bool> >,
703
                               hash<pair<size_t, size_t> > > edges_by_end_deg_t;
704
    edges_by_end_deg_t _edges_by_source, _edges_by_target;
705
706
707

protected:
    const Graph& _g;
708
709
};

710
711
712
713
714
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
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;
};

863
864
865
} // graph_tool namespace

#endif // GRAPH_REWIRING_HH