Fix rewiring bias bug for undirected graphs

This fixes a problem, where the underlying edge directionality would
cause correlations to arise when rewiring undirected graphs. Now each
edge is correctly considered in both possible orientations, which are
chosen randomly.

src/graph/generation/graph_rewiring.hh

@@ -47,15 +47,40 @@ bool is_adjacent(typename graph_traits<Graph>::vertex_descriptor u,
     return false;
 }
 
+template <class Graph>
+typename graph_traits<Graph>::vertex_descriptor
+source(const pair<typename graph_traits<Graph>::edge_descriptor, bool>& e,
+       const Graph& g)
+{
+    if (e.second)
+        return target(e.first, g);
+    else
+        return source(e.first, g);
+}
+
+template <class Graph>
+typename graph_traits<Graph>::vertex_descriptor
+target(const pair<typename graph_traits<Graph>::edge_descriptor, bool>& e,
+       const Graph& g)
+{
+    if (e.second)
+        return source(e.first, g);
+    else
+        return target(e.first, g);
+}
+
 // 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
 struct swap_edge_triad
 {
     template <class Graph, class NewEdgeMap>
-    static bool parallel_check(typename graph_traits<Graph>::edge_descriptor e,
-                               typename graph_traits<Graph>::edge_descriptor se,
-                               typename graph_traits<Graph>::edge_descriptor te,
-                               NewEdgeMap edge_is_new, const Graph &g)
+    static bool
+    parallel_check(const typename graph_traits<Graph>::edge_descriptor& e,
+                   const pair<typename graph_traits<Graph>::edge_descriptor,
+                              bool>& se,
+                   const pair<typename graph_traits<Graph>::edge_descriptor,
+                              bool>& te,
+                   NewEdgeMap edge_is_new, const Graph &g)
     {
         // We want to check that if we swap the source of 'e' with the source of
         // 'se', and the target of 'te' with the target of 'e', as such
@@ -78,18 +103,18 @@ struct swap_edge_triad
             se_t = target(se, g);      // source edge target
 
 
-        if (edge_is_new[se] && (ns == s) && (nt == se_t))
+        if (edge_is_new[se.first] && (ns == s) && (nt == se_t))
             return true; // e is parallel to se after swap
-        if (edge_is_new[te] && (te_s == ns) && (nt == t))
+        if (edge_is_new[te.first] && (te_s == ns) && (nt == t))
             return true; // e is parallel to te after swap
-        if (edge_is_new[te] && edge_is_new[se] && (te != se) &&
+        if (edge_is_new[te.first] && edge_is_new[se.first] && (te != se) &&
              (s == te_s) && (t == se_t))
             return true; // se is parallel to te after swap
         if (is_adjacent_in_new(ns,  nt, edge_is_new, g))
             return true; // e would clash with an existing (new) edge
-        if (edge_is_new[te] && is_adjacent_in_new(te_s, t, edge_is_new, g))
+        if (edge_is_new[te.first] && is_adjacent_in_new(te_s, t, edge_is_new, g))
             return true; // te would clash with an existing (new) edge
-        if (edge_is_new[se] && is_adjacent_in_new(s, se_t, edge_is_new, g))
+        if (edge_is_new[se.first] && is_adjacent_in_new(s, se_t, edge_is_new, g))
             return true; // se would clash with an existing (new) edge
         return false; // the coast is clear - hooray!
     }
@@ -112,9 +137,11 @@ struct swap_edge_triad
     }
 
     template <class Graph, class EdgeIndexMap, class EdgesType>
-    void operator()(typename graph_traits<Graph>::edge_descriptor e,
-                    typename graph_traits<Graph>::edge_descriptor se,
-                    typename graph_traits<Graph>::edge_descriptor te,
+    void operator()(const typename graph_traits<Graph>::edge_descriptor& e,
+                    const pair<typename graph_traits<Graph>::edge_descriptor,
+                               bool>& se,
+                    const pair<typename graph_traits<Graph>::edge_descriptor,
+                               bool>& te,
                     EdgesType& edges, EdgeIndexMap edge_index, Graph& g)
     {
         // swap the source of the edge 'e' with the source of edge 'se' and the
@@ -129,21 +156,21 @@ struct swap_edge_triad
 
         // split cases where different combinations of the three edges are
         // the same
-        if(se != te)
+        if(se.first != te.first)
         {
             ne = add_edge(source(se, g), target(te, g), g).first;
-            if(e != se)
+            if(e != se.first)
             {
                 nse = add_edge(source(e, g), target(se, g), g).first;
-                edge_index[nse] = edge_index[se];
-                remove_edge(se, g);
+                edge_index[nse] = edge_index[se.first];
+                remove_edge(se.first, g);
                 edges[edge_index[nse]] = nse;
             }
-            if(e != te)
+            if(e != te.first)
             {
                 nte = add_edge(source(te, g), target(e, g), g).first;
-                edge_index[nte] = edge_index[te];
-                remove_edge(te, g);
+                edge_index[nte] = edge_index[te.first];
+                remove_edge(te.first, g);
                 edges[edge_index[nte]] = nte;
             }
             edge_index[ne] = edge_index[e];
@@ -152,11 +179,11 @@ struct swap_edge_triad
         }
         else
         {
-            if(e != se)
+            if(e != se.first)
             {
                 // se and te are the same. swapping indexes only.
-                swap(edge_index[se], edge_index[e]);
-                edges[edge_index[se]] = se;
+                swap(edge_index[se.first], edge_index[e]);
+                edges[edge_index[se.first]] = se.first;
                 edges[edge_index[e]] = e;
             }
         }
@@ -382,6 +409,8 @@ class RewireStrategyBase
 {
 public:
     typedef typename graph_traits<Graph>::edge_descriptor edge_t;
+    typedef typename graph_traits<Graph>::edge_descriptor vertex_t;
+
     typedef typename EdgeIndexMap::value_type index_t;
 
     RewireStrategyBase(Graph& g, EdgeIndexMap edge_index, rng_t& rng)
@@ -392,16 +421,18 @@ public:
                     bool self_loops, bool parallel_edges)
     {
         // where should we sample the edges from
-        vector<index_t>* edges_source=0, *edges_target=0;
+        vector<pair<index_t,bool> >* edges_source=0, *edges_target=0;
         static_cast<RewireStrategy*>(this)->get_edges(e, edges_source,
                                                       edges_target);
 
+
         //try randomly drawn pairs of edges until one satisfies all the
         //consistency checks
         bool found = false;
-        edge_t es = e, et = e;
+        pair<edge_t,bool> es, et;
         typedef random_permutation_iterator
-            <typename vector<index_t>::iterator, rng_t> random_edge_iter;
+            <typename vector<pair<index_t,bool> >::iterator, rng_t>
+            random_edge_iter;
 
         random_edge_iter esi(edges_source->begin(), edges_source->end(),
                              _rng),
@@ -409,9 +440,9 @@ public:
                              _rng);
         for (; esi != esi_end && !found; ++esi)
         {
-            if (!is_edge[*esi])
+            if (!is_edge[(*esi).first])
                 continue;
-            es = edges[*esi];
+            es = make_pair(edges[(*esi).first], (*esi).second);
             static_cast<RewireStrategy*>(this)->check_source_edge(es, e);
 
             if(!self_loops) // reject self-loops if not allowed
@@ -426,9 +457,9 @@ public:
                                  _rng);
             for (; eti != eti_end && !found; ++eti)
             {
-                if (!is_edge[*eti])
+                if (!is_edge[(*eti).first])
                     continue;
-                et = edges[*eti];
+                et = make_pair(edges[(*eti).first], (*eti).second);
                 static_cast<RewireStrategy*>(this)->check_target_edge(et, e);
 
                 if (!self_loops) // reject self-loops if not allowed
@@ -483,25 +514,36 @@ public:
                          rng_t& rng)
         : base_t(g, edge_index, rng)
     {
+        //bernoulli_distribution coin(0.5);
         typename graph_traits<Graph>::edge_iterator e_i, e_i_end;
         for (tie(e_i, e_i_end) = edges(g); e_i != e_i_end; ++e_i)
-            _all_edges.push_back(edge_index[*e_i]);
+        {
+            // For undirected graphs, there is no difference between source
+            // and target, and each edge will appear _twice_ on the list
+            // once for each different ordering of source and target.
+
+            _all_edges.push_back(make_pair(edge_index[*e_i], false));
+            if (!is_directed::apply<Graph>::type::value)
+                _all_edges.push_back(make_pair(edge_index[*e_i], true));
+        }
         _all_edges2 = _all_edges;
+
     }
 
-    void get_edges(const edge_t& e, vector<index_t>*& edges_source,
-                   vector<index_t>*& edges_target)
+    void get_edges(const edge_t& e,
+                   vector<pair<index_t,bool> >*& edges_source,
+                   vector<pair<index_t,bool> >*& edges_target)
     {
         edges_source = &_all_edges;
         edges_target = &_all_edges2;
     }
 
-    void check_source_edge(edge_t& se, const edge_t& e) {}
-    void check_target_edge(edge_t& te, const edge_t& e) {}
+    void check_source_edge(pair<edge_t,bool>& se, const edge_t& e) {}
+    void check_target_edge(pair<edge_t,bool>& te, const edge_t& e) {}
 
 private:
-    vector<index_t> _all_edges;
-    vector<index_t> _all_edges2;
+    vector<pair<index_t,bool> > _all_edges;
+    vector<pair<index_t,bool> > _all_edges2;
 };
 
 
@@ -527,34 +569,39 @@ public:
     CorrelatedRewireStrategy (Graph& g, EdgeIndexMap edge_index,
                               rng_t& rng) : base_t(g, edge_index, rng), _g(g)
     {
-        int i, N = num_vertices(_g);
-        for (i = 0; i < N; ++i)
+        typename graph_traits<Graph>::edge_iterator e_i, e_i_end;
+        for (tie(e_i, e_i_end) = edges(g); e_i != e_i_end; ++e_i)
         {
-            vertex_t v = vertex(i, _g);
-            if (v == graph_traits<Graph>::null_vertex())
-                continue;
-            typename graph_traits<Graph>::out_edge_iterator e_i, e_i_end;
-            for (tie(e_i, e_i_end) = out_edges(v, _g); e_i != e_i_end; ++e_i)
-            {
-                // For undirected graphs, there is no difference between source
-                // and target, and each edge will appear _twice_ on the lists
-                // below, once for each different ordering of source and target.
+            // For undirected graphs, there is no difference between source and
+            // target, and each edge will appear _twice_ on the lists below,
+            // once for each different ordering of source and target.
 
-                _edges_by_source
-                    [make_pair(in_degreeS()(source(*e_i, _g), _g),
-                               out_degree(source(*e_i, _g), _g))]
-                    .push_back(edge_index[*e_i]);
+            _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));
 
-                _edges_by_target
+            _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
                     [make_pair(in_degreeS()(target(*e_i, _g), _g),
                                out_degree(target(*e_i, _g), _g))]
-                    .push_back(edge_index[*e_i]);
+                    .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));
             }
         }
     }
 
-    void get_edges(const edge_t& e, vector<index_t>*& edges_source,
-                   vector<index_t>*& edges_target)
+    void get_edges(const edge_t& e, vector<pair<index_t,bool> >*& edges_source,
+                   vector<pair<index_t,bool> >*& edges_target)
     {
         pair<size_t, size_t> deg_source =
             make_pair(in_degreeS()(source(e, _g), _g),
@@ -570,59 +617,59 @@ public:
     }
 
 
-    void check_source_edge(edge_t& se, const edge_t& e)
+    void check_source_edge(pair<edge_t,bool>& se, const edge_t& e)
     {
         check_source_edge_if_undirected
             (se, e, typename is_directed::apply<Graph>::type());
     }
-    void check_target_edge(edge_t& te, const edge_t& e)
+    void check_target_edge(pair<edge_t,bool>& te, const edge_t& e)
     {
         check_target_edge_if_undirected
             (te, e, typename is_directed::apply<Graph>::type());
     }
 
-    void check_source_edge_if_undirected(edge_t& se, const edge_t& e,
+    void check_source_edge_if_undirected(pair<edge_t,bool>& se, const edge_t& e,
                                          boost::true_type) {}
-    void check_target_edge_if_undirected(edge_t& te, const edge_t& e,
+    void check_target_edge_if_undirected(pair<edge_t,bool>& te, const edge_t& e,
                                          boost::true_type) {}
 
-    void check_source_edge_if_undirected(edge_t& se, const edge_t& e,
+    void check_source_edge_if_undirected(pair<edge_t,bool>& se, const edge_t& e,
                                          boost::false_type)
     {
         // check if the edge direction is correct, otherwise invert it.
         pair<size_t, size_t> deg_source1 =
-            make_pair(in_degreeS()(source(e, _g), _g),
+            make_pair(out_degree(source(e, _g), _g),
                       out_degree(source(e, _g), _g));
 
         pair<size_t, size_t> deg_source2 =
-            make_pair(in_degreeS()(source(se, _g), _g),
+            make_pair(out_degree(source(se, _g), _g),
                       out_degree(source(se, _g), _g));
 
         if (deg_source1 != deg_source2)
-            se = edge_t(se, !se.IsInverted());
+            se.first = edge_t(se.first, !se.first.IsInverted());
     }
 
-    void check_target_edge_if_undirected(edge_t& te, const edge_t& e,
+    void check_target_edge_if_undirected(pair<edge_t,bool>& te, const edge_t& e,
                                          boost::false_type)
     {
         // check if the edge direction is correct, otherwise invert it.
         pair<size_t, size_t> deg_target1 =
-            make_pair(in_degreeS()(target(e, _g), _g),
+            make_pair(out_degree(target(e, _g), _g),
                       out_degree(target(e, _g), _g));
 
         pair<size_t, size_t> deg_target2 =
-            make_pair(in_degreeS()(target(te, _g), _g),
+            make_pair(out_degree(target(te, _g), _g),
                       out_degree(target(te, _g), _g));
 
         if (deg_target1 != deg_target2)
-            te = edge_t(te, !te.IsInverted());
+            te.first = edge_t(te.first, !te.first.IsInverted());
     }
 
 private:
-    typedef tr1::unordered_map<pair<size_t, size_t>, vector<index_t>,
+    typedef tr1::unordered_map<pair<size_t, size_t>,
+                               vector<pair<index_t, bool> >,
                                hash<pair<size_t, size_t> > > edges_by_end_deg_t;
     edges_by_end_deg_t _edges_by_source, _edges_by_target;
-    vector<index_t> _temp;
 
 protected:
     const Graph& _g;