graph_rewiring.hh 24 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>
Tiago Peixoto's avatar
Tiago Peixoto committed
147
    static bool
148
149
150
    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
151
152
153
154
155
156
157
    {
        // 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)
        //
158
        // no parallel edges are introduced.
Tiago Peixoto's avatar
Tiago Peixoto committed
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
            nt = target(te, edges, g),        // new target
            te_s = source(te, edges, g);      // target edge source
165

166
        if (is_adjacent(s,  nt, g))
167
            return true; // e would clash with an existing (new) edge
168
        if (is_adjacent(te_s, t, g))
169
170
171
172
            return true; // te would clash with an existing (new) edge
        return false; // the coast is clear - hooray!
    }

173
174
175
176
177
    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)
178
    {
179
180
        // swap the source of the edge 'e' with the source of edge 'se', as
        // such:
181
        //
182
183
        //  (s)    -e--> (t)          (s)    -e--> (nt)
        //  (te_s) -te-> (nt)   =>    (te_s) -te-> (t)
184

185
186
        if (e == te.first)
            return;
187

188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
        // 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;
203
    }
204

205
206
};

Tiago Peixoto's avatar
Tiago Peixoto committed
207
// used for verbose display
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
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;
    }
}


Tiago Peixoto's avatar
Tiago Peixoto committed
228

229
// main rewire loop
230
231
template <template <class Graph, class EdgeIndexMap, class CorrProb>
          class RewireStrategy>
232
233
struct graph_rewire
{
234
235
    template <class Graph, class EdgeIndexMap, class CorrProb>
    void operator()(Graph& g, EdgeIndexMap edge_index, CorrProb corr_prob,
236
237
238
239
                    bool self_loops, bool parallel_edges,
                    pair<size_t, bool> iter_sweep, bool verbose, size_t& pcount,
                    rng_t& rng)
        const
240
241
242
243
244
    {
        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
245
        vector<size_t> edges_sequence;
246
247
        typename graph_traits<Graph>::edge_iterator e, e_end;
        for (tie(e, e_end) = boost::edges(g); e != e_end; ++e)
248
        {
249
250
251
            if (edge_index[*e] >= edges.size())
                edges.resize(edge_index[*e] + 1);
            edges[edge_index[*e]] = *e;
Tiago Peixoto's avatar
Tiago Peixoto committed
252
            edges_sequence.push_back(edge_index[*e]);
253
254
        }

255
256
        RewireStrategy<Graph, EdgeIndexMap, CorrProb>
            rewire(g, edge_index, edges, corr_prob, rng);
257

Tiago Peixoto's avatar
Tiago Peixoto committed
258
259
260
        typedef random_permutation_iterator<typename vector<size_t>::iterator,
                                            rng_t>
            random_edge_iter;
261

262
263
264
265
        size_t niter;
        bool no_sweep;
        tie(niter, no_sweep) = iter_sweep;
        pcount = 0;
266
267
        if (verbose)
            cout << "rewiring edges: ";
Tiago Peixoto's avatar
Tiago Peixoto committed
268
        stringstream str;
269
        for (size_t i = 0; i < niter; ++i)
270
        {
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
            random_edge_iter
                ei(edges_sequence.begin(), edges_sequence.end(), rng),
                ei_end(edges_sequence.end(), edges_sequence.end(), rng);

            size_t count = 0;
            // for each edge rewire its source or target
            for (; ei != ei_end; ++ei)
            {
                if (verbose)
                    print_progress(i, niter, count++, 
                                   no_sweep ? 1 : edges_sequence.size(), str);
                bool success = rewire(*ei, self_loops, parallel_edges);
                if (!success)
                    ++pcount;

                if (no_sweep)
                    break;
            }
289
        }
290
291
        if (verbose)
            cout << endl;
292
    }
293
294
};

Tiago Peixoto's avatar
Tiago Peixoto committed
295

296
297
// this will rewire the edges so that the resulting graph will be entirely
// random (i.e. Erdos-Renyi)
298
template <class Graph, class EdgeIndexMap, class CorrProb>
299
300
301
302
303
304
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
305
    ErdosRewireStrategy(Graph& g, EdgeIndexMap edge_index,
306
                        vector<edge_t>& edges, CorrProb corr, rng_t& rng)
Tiago Peixoto's avatar
Tiago Peixoto committed
307
308
        : _g(g), _edge_index(edge_index), _edges(edges),
          _vertices(HardNumVertices()(g)), _rng(rng)
309
310
311
312
313
314
315
    {
        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;
    }

316
    bool operator()(size_t e, bool self_loops, bool parallel_edges)
317
    {
318
        //try randomly drawn pairs of vertices
319
320
321
322
        typedef random_permutation_iterator
            <typename graph_traits<Graph>::vertex_iterator, rng_t>
            random_vertex_iter;

323
        tr1::uniform_int<size_t> sample(0, _vertices.size() - 1);
324
325
326
327
328
329
        typename graph_traits<Graph>::vertex_descriptor s, t;
        while (true)
        {
            s = sample(_rng);
            t = sample(_rng);

330
331
            // reject self-loops if not allowed
            if(s == t && !self_loops)
332
333
334
                continue;
            break;
        }
335
336
337
338
339

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

340
        edge_t ne = add_edge(s, t, _g).first;
341
342
343
        _edge_index[ne] = e;
        remove_edge(_edges[e], _g);
        _edges[e] = ne;
344
345

        return true;
346
347
348
349
350
    }

private:
    Graph& _g;
    EdgeIndexMap _edge_index;
Tiago Peixoto's avatar
Tiago Peixoto committed
351
    vector<edge_t>& _edges;
352
353
354
355
356
    vector<typename graph_traits<Graph>::vertex_descriptor> _vertices;
    rng_t& _rng;
};


357
358
359
// 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.
360
361
template <class Graph, class EdgeIndexMap, class RewireStrategy>
class RewireStrategyBase
362
363
364
{
public:
    typedef typename graph_traits<Graph>::edge_descriptor edge_t;
365
366
    typedef typename graph_traits<Graph>::edge_descriptor vertex_t;

367
368
    typedef typename EdgeIndexMap::value_type index_t;

Tiago Peixoto's avatar
Tiago Peixoto committed
369
370
    RewireStrategyBase(Graph& g, EdgeIndexMap edge_index, vector<edge_t>& edges,
                       rng_t& rng)
371
        : _g(g), _edge_index(edge_index), _edges(edges),  _rng(rng) {}
372

373
    bool operator()(size_t e, bool self_loops, bool parallel_edges)
374
    {
Tiago Peixoto's avatar
Tiago Peixoto committed
375
        RewireStrategy& self = *static_cast<RewireStrategy*>(this);
376

377
378
        // try randomly drawn pairs of edges and check if they satisfy all the
        // consistency checks
379

380
381
382
383
        // rewire target
        pair<size_t, bool> et = self.get_target_edge(e);

        if (!self_loops) // reject self-loops if not allowed
384
        {
385
386
387
388
            if((source(_edges[e], _g) == target(et, _edges, _g)) ||
               (target(_edges[e], _g) == source(et, _edges, _g)))
                return false;
        }
389

390
391
392
393
394
        // reject parallel edges if not allowed
        if (!parallel_edges && (et.first != e))
        {
            if (swap_edge::parallel_check_target(e, et, _edges, _g))
                return false;
395
        }
396
397
398
399
400
401
402
403
404
405
406
407

        if (e != et.first)
        {
            self.update_edge(e, false);
            self.update_edge(et.first, false);
            
            swap_edge::swap_target(e, et, _edges, _edge_index, _g);

            self.update_edge(e, true);
            self.update_edge(et.first, true);
        }
        return true;
408
409
    }

Tiago Peixoto's avatar
Tiago Peixoto committed
410
protected:
411
412
    Graph& _g;
    EdgeIndexMap _edge_index;
Tiago Peixoto's avatar
Tiago Peixoto committed
413
    vector<edge_t>& _edges;
414
    rng_t& _rng;
415
416
};

417
418
// this will rewire the edges so that the combined (in, out) degree distribution
// will be the same, but all the rest is random
419
template <class Graph, class EdgeIndexMap, class CorrProb>
420
421
class RandomRewireStrategy:
    public RewireStrategyBase<Graph, EdgeIndexMap,
422
423
                              RandomRewireStrategy<Graph, EdgeIndexMap,
                                                   CorrProb> >
424
425
{
public:
426
    typedef RewireStrategyBase<Graph, EdgeIndexMap,
427
428
                               RandomRewireStrategy<Graph, EdgeIndexMap,
                                                    CorrProb> >
429
430
431
432
433
        base_t;

    typedef Graph graph_t;
    typedef EdgeIndexMap edge_index_t;

434
435
    typedef typename graph_traits<Graph>::vertex_descriptor vertex_t;
    typedef typename graph_traits<Graph>::edge_descriptor edge_t;
436
    typedef typename EdgeIndexMap::value_type index_t;
437

Tiago Peixoto's avatar
Tiago Peixoto committed
438
439
440
441
442
443
444
445
446
447
    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;

448
    RandomRewireStrategy(Graph& g, EdgeIndexMap edge_index,
449
                         vector<edge_t>& edges, CorrProb corr, rng_t& rng)
Tiago Peixoto's avatar
Tiago Peixoto committed
450
        : base_t(g, edge_index, edges, rng), _g(g)
451
    {
Tiago Peixoto's avatar
Tiago Peixoto committed
452
453
454
455
456
        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)
        {
457
            if (out_degree(*v_i, _g) > E / 3)
Tiago Peixoto's avatar
Tiago Peixoto committed
458
459
                _hub_non_adjacent[*v_i] = edge_set_t();
        }
460

461
        typename graph_traits<Graph>::edge_iterator e_i, e_i_end;
Tiago Peixoto's avatar
Tiago Peixoto committed
462
        for (tie(e_i, e_i_end) = boost::edges(g); e_i != e_i_end; ++e_i)
463
        {
Tiago Peixoto's avatar
Tiago Peixoto committed
464
465
466
            // 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.
467

Tiago Peixoto's avatar
Tiago Peixoto committed
468
            _all_edges.push_back(make_pair(edge_index[*e_i], false));
469
            if (!is_directed::apply<Graph>::type::value)
Tiago Peixoto's avatar
Tiago Peixoto committed
470
                _all_edges.push_back(make_pair(edge_index[*e_i], true));
471
            update_edge(edge_index[*e_i], true);
Tiago Peixoto's avatar
Tiago Peixoto committed
472
473
474
        }
    }

475
    pair<size_t,bool> get_target_edge(size_t e)
Tiago Peixoto's avatar
Tiago Peixoto committed
476
    {
477
        tr1::unordered_map<size_t, edge_set_t>& hub = _hub_non_adjacent;
478

Tiago Peixoto's avatar
Tiago Peixoto committed
479
        typeof(hub.begin()) iter =
480
            hub.find(source(base_t::_edges[e], _g));
Tiago Peixoto's avatar
Tiago Peixoto committed
481
482
483
484
485
486
487
488
489
        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
490
                return *(eset.get<random_index>().begin() + i);
491
492
        }

493
        tr1::uniform_int<> sample(0, _all_edges.size() - 1);
494
        return _all_edges[sample(base_t::_rng)];
495
    }
496

Tiago Peixoto's avatar
Tiago Peixoto committed
497

498
    void update_edge_list(index_t v, const pair<size_t, bool>& e,
Tiago Peixoto's avatar
Tiago Peixoto committed
499
500
501
                          edge_set_t& edge_set, bool src, bool insert)
    {
        pair<index_t,bool> ep;
502
503
        if ((src && source(e, base_t::_edges, _g) != v) ||
            (!src && target(e, base_t::_edges, _g) != v))
Tiago Peixoto's avatar
Tiago Peixoto committed
504
505
        {
            if (insert)
506
                edge_set.get<hash_index>().insert(e);
Tiago Peixoto's avatar
Tiago Peixoto committed
507
            else
508
                edge_set.get<hash_index>().erase(e);
Tiago Peixoto's avatar
Tiago Peixoto committed
509
510
511
        }
    }

512
    void update_edge(size_t e, bool insert)
Tiago Peixoto's avatar
Tiago Peixoto committed
513
514
515
516
517
518
519
520
521
522
    {
        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);
        }
523
    }
524

525
private:
Tiago Peixoto's avatar
Tiago Peixoto committed
526
527
    Graph& _g;
    EdgeIndexMap _edge_index;
528
    vector<pair<index_t,bool> > _all_edges;
Tiago Peixoto's avatar
Tiago Peixoto committed
529

530
    tr1::unordered_map<size_t, edge_set_t> _hub_non_adjacent;
531
};
532

533
534
535

// 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
536
template <class Graph, class EdgeIndexMap, class CorrProb>
537
538
class CorrelatedRewireStrategy:
    public RewireStrategyBase<Graph, EdgeIndexMap,
539
540
                              CorrelatedRewireStrategy<Graph, EdgeIndexMap,
                                                       CorrProb> >
541
542
{
public:
543
    typedef RewireStrategyBase<Graph, EdgeIndexMap,
544
545
                               CorrelatedRewireStrategy<Graph, EdgeIndexMap,
                                                        CorrProb> >
546
547
548
549
550
        base_t;

    typedef Graph graph_t;
    typedef EdgeIndexMap edge_index_t;

551
552
553
554
    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;

555
    CorrelatedRewireStrategy (Graph& g, EdgeIndexMap edge_index,
556
                              vector<edge_t>& edges, CorrProb corr, rng_t& rng)
Tiago Peixoto's avatar
Tiago Peixoto committed
557
        : base_t(g, edge_index, edges, rng), _g(g)
558
    {
559
        typename graph_traits<Graph>::edge_iterator e_i, e_i_end;
Tiago Peixoto's avatar
Tiago Peixoto committed
560
        for (tie(e_i, e_i_end) = boost::edges(g); e_i != e_i_end; ++e_i)
561
        {
562
            // For undirected graphs, there is no difference between source and
Tiago Peixoto's avatar
Tiago Peixoto committed
563
            // target, and each edge will appear _twice_ in the lists below,
564
            // once for each different ordering of source and target.
565

566
567
568
569
570
571
572
573
574
575
576
            _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_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));
577
578
579
            }
        }
    }
580

581
    pair<size_t,bool> get_target_edge(size_t e)
582
    {
Tiago Peixoto's avatar
Tiago Peixoto committed
583
        pair<size_t, size_t> deg =
584
585
586
            make_pair(in_degreeS()(target(base_t::_edges[e], _g), _g),
                      out_degree(target(base_t::_edges[e], _g), _g));
        edges_by_end_deg_t& edges = _edges_by_target;
Tiago Peixoto's avatar
Tiago Peixoto committed
587
        typename edges_by_end_deg_t::mapped_type& elist = edges[deg];
588
        tr1::uniform_int<> sample(0, elist.size() - 1);
589
        return elist[sample(base_t::_rng)];
Tiago Peixoto's avatar
Tiago Peixoto committed
590
    }
591

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

594
private:
595
596
    typedef tr1::unordered_map<pair<size_t, size_t>,
                               vector<pair<index_t, bool> >,
597
598
                               hash<pair<size_t, size_t> > >
        edges_by_end_deg_t;
599
    edges_by_end_deg_t _edges_by_target;
600
601
602

protected:
    const Graph& _g;
603
604
};

605
606
607
608
609
610
611
612
template <class Graph, class EdgeIndexMap, class CorrProb>
class ProbabilisticRewireStrategy:
    public RewireStrategyBase<Graph, EdgeIndexMap,
                              ProbabilisticRewireStrategy<Graph, EdgeIndexMap,
                                                          CorrProb> >
{
public:
    typedef RewireStrategyBase<Graph, EdgeIndexMap,
613
614
                               ProbabilisticRewireStrategy<Graph, EdgeIndexMap,
                                                           CorrProb> >
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
        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)
    {
        set<pair<size_t, size_t> > deg_set;
632
        typename graph_traits<Graph>::edge_iterator e_i, e_i_end;
633
634
635
636
637
638
639
640
641
642
643
644
        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)
            for (typeof(deg_set.begin()) t_iter = deg_set.begin();
                 t_iter != deg_set.end(); ++t_iter)
            {
645
646
647
648
                double p = corr_prob(*s_iter, *t_iter);
                if (isnan(p) || isinf(p))
                    p = 0;
                _probs[make_pair(*s_iter, *t_iter)] = p;
649
650
651
652
653
654
655
656
657
658
659
660
                _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_target_edge(size_t e)
    {
661
662
663
664
665
666
667
668
669
670
671
672
673
        deg_t s_deg = get_deg(source(base_t::_edges[e], _g), _g);
        deg_t t_deg = get_deg(target(base_t::_edges[e], _g), _g);
     
        deg_t deg = _sample_target_deg[s_deg](base_t::_rng);
        pair<size_t, bool> ep = _sample_edge_by_target_deg[deg](base_t::_rng);

        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);

        double pi = _probs[make_pair(ep_s_deg, ep_t_deg)];
        double pf = _probs[make_pair(ep_s_deg, t_deg)];
   
        if (pi > pf && pf > 0)
674
        {
675
676
677
678
679
680
681
            double a = exp(log(pf) - log(pi));
   
            tr1::variate_generator<rng_t&, tr1::uniform_real<> >
                sample(base_t::_rng, tr1::uniform_real<>(0.0, 1.0));
            double r = sample();
            if (r > a)
                return make_pair(e, false); // reject
682
        }
683
        
684
685
686
        return ep;
    }

687
    void update_edge(size_t e, bool insert)
688
    {
689
690
        deg_t s_deg = get_deg(source(base_t::_edges[e], _g), _g);
        deg_t t_deg = get_deg(target(base_t::_edges[e], _g), _g);
691
692
693

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

694
        if (insert)
695
696
697
698
699
700
701
702
        {
            _sample_edge_by_target_deg[t_deg].Insert(ep);
            if (!is_directed::apply<Graph>::type::value)
            {
                ep.second = true;
                _sample_edge_by_target_deg[s_deg].Insert(ep);
            }
        }
703
        else
704
705
706
707
708
709
710
711
712
713
714
715
716
717
        {
            _sample_edge_by_target_deg[t_deg].Remove(ep);
            if (!is_directed::apply<Graph>::type::value)
            {
                ep.second = true;
                _sample_edge_by_target_deg[s_deg].Remove(ep);
            }
        }
    }

private:
    Graph& _g;
    EdgeIndexMap _edge_index;

718
    typedef tr1::unordered_map<deg_t, WeightedSampler<deg_t>, hash<deg_t> >
719
720
        deg_sampler_t;

721
    typedef tr1::unordered_map<deg_t, Sampler<pair<size_t, bool> >, hash<deg_t> >
722
723
724
725
        edge_sampler_t;

    edge_sampler_t _sample_edge_by_target_deg;
    deg_sampler_t _sample_target_deg;
726
727
    tr1::unordered_map<pair<deg_t, deg_t>, double, hash<pair<deg_t, deg_t> > >
        _probs;
728
729
};

730
731
732
} // graph_tool namespace

#endif // GRAPH_REWIRING_HH