Include probabilistic random_rewire() strategy

This includes also some internal refactoring, and bug fixes.

src/graph/generation/graph_rewiring.cc

@@ -19,35 +19,58 @@
 #include "graph_filtering.hh"
 
 #include <tr1/random>
-typedef std::tr1::mt19937 rng_t;
 
 #include "graph_rewiring.hh"
 
 #include <boost/bind.hpp>
+#include <boost/python.hpp>
 
 using namespace graph_tool;
 using namespace boost;
 
+class PythonFuncWrap
+{
+public:
+    PythonFuncWrap(python::object o): _o(o) {}
+
+    double operator()(pair<size_t, size_t> deg, pair<size_t, size_t> degl) const
+    {
+        python::object ret = _o(python::make_tuple(deg.first, deg.second),
+                                python::make_tuple(degl.first, degl.second));
+        return python::extract<double>(ret);
+    }
+
+private:
+    python::object _o;
+};
+
 void random_rewire(GraphInterface& gi, string strat, bool self_loops,
-                   bool parallel_edges, size_t seed)
+                   bool parallel_edges, python::object corr_prob, size_t seed,
+                   bool verbose)
 {
     rng_t rng(static_cast<rng_t::result_type>(seed));
+    PythonFuncWrap corr(corr_prob);
 
     if (strat == "erdos")
         run_action<graph_tool::detail::never_reversed>()
             (gi, bind<void>(graph_rewire<ErdosRewireStrategy>(),
-                            _1, gi.GetEdgeIndex(), ref(rng), self_loops,
-                            parallel_edges))();
+                            _1, gi.GetEdgeIndex(), ref(corr), ref(rng),
+                            self_loops, parallel_edges, verbose))();
     else if (strat == "uncorrelated")
         run_action<graph_tool::detail::never_reversed>()
             (gi, bind<void>(graph_rewire<RandomRewireStrategy>(),
-                            _1, gi.GetEdgeIndex(), ref(rng), self_loops,
-                            parallel_edges))();
+                            _1, gi.GetEdgeIndex(), ref(corr), ref(rng),
+                            self_loops, parallel_edges, verbose))();
     else if (strat == "correlated")
         run_action<graph_tool::detail::never_reversed>()
             (gi, bind<void>(graph_rewire<CorrelatedRewireStrategy>(),
-                            _1, gi.GetEdgeIndex(), ref(rng), self_loops,
-                            parallel_edges))();
+                            _1, gi.GetEdgeIndex(), ref(corr), ref(rng),
+                            self_loops, parallel_edges, verbose))();
+    else if (strat == "probabilistic")
+        run_action<>()
+            (gi, bind<void>(graph_rewire<ProbabilisticRewireStrategy>(),
+                            _1, gi.GetEdgeIndex(), ref(corr), ref(rng),
+                            self_loops, parallel_edges, verbose))();
     else
         throw ValueException("invalid random rewire strategy: " + strat);
 }

src/graph/generation/sampler.hh

+// graph-tool -- a general graph modification and manipulation thingy
+//
+// Copyright (C) 2007  Tiago de Paula Peixoto <tiago@forked.de>
+//
+// 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 SAMPLER_HH
+#define SAMPLER_HH
+
+#include <tr1/random>
+#include <iostream>
+
+namespace graph_tool
+{
+using namespace std;
+using namespace boost;
+
+typedef tr1::mt19937 rng_t;
+
+// utility class to sample uniformly from a collection of values
+template <class ValueType>
+class Sampler
+{
+public:
+    Sampler(bool biased=false): _biased(biased), _erased_prob(0) {}
+
+    template <class Iterator>
+    Sampler(Iterator iter, Iterator end):
+        _biased(false)
+    {
+        for(; iter != end; ++iter)
+        {
+            _candidates.push_back(*iter);
+            _candidates_set.insert(make_pair(*iter, _candidates.size()-1));
+        }
+        //assert(!_candidates.empty());
+    }
+
+    void Insert(const ValueType& v, double p = 0.0)
+    {
+        _candidates.push_back(v);
+        _candidates_set.insert(make_pair(v, _candidates.size()-1));
+        if (_biased)
+        {
+            if (_probs.size() > 0)
+                _probs.push_back(_probs.back()+p);
+            else
+                _probs.push_back(p);
+            _erased.push_back(false);
+        }
+    }
+
+    bool HasValue(const ValueType& v)
+    {
+        typeof(_candidates_set.begin()) iter, end;
+        tie(iter, end) = _candidates_set.equal_range(v);
+        return (iter != end);
+    }
+
+    void Remove(const ValueType& v)
+    {
+        typeof(_candidates_set.begin()) iter, end, temp;
+        tie(iter, end) = _candidates_set.equal_range(v);
+        //assert(iter != end);
+
+        if (_biased)
+        {
+            while(_erased[iter->second])
+            {
+                temp = iter++;
+                _candidates_set.erase(temp);
+            }
+
+            size_t index = iter->second;
+            _erased[index] = true;
+            _erased_prob += (index > 0) ?
+                _probs[index]-_probs[index-1] : _probs[index];
+        }
+        else
+        {
+            size_t index = iter->second;
+            temp = _candidates_set.find(_candidates.back());
+            swap(_candidates[index], _candidates.back());
+            _candidates.pop_back();
+            if (!_candidates.empty() && temp != iter)
+            {
+                _candidates_set.erase(temp);
+                _candidates_set.insert(make_pair(_candidates[index], index));
+            }
+        }
+        _candidates_set.erase(iter);
+
+        clean();
+    }
+
+    bool Empty()
+    {
+        return _candidates.empty();
+    }
+
+    size_t Size()
+    {
+        return _candidates.size();
+    }
+
+    ValueType operator()(rng_t& rng, bool remove = false)
+    {
+        //assert(!_candidates.empty());
+        if (!_biased)
+        {
+            tr1::uniform_int<> sample(0, _candidates.size() - 1);
+            int i = sample(rng);
+            if (remove)
+            {
+                swap(_candidates[i], _candidates.back());
+                ValueType ret = _candidates.back();
+                _candidates.pop_back();
+                return ret;
+            }
+            else
+            {
+                return _candidates[i];
+            }
+        }
+        else
+        {
+            size_t i = 0;
+            do
+            {
+                if (_probs.back() > 0)
+                {
+                    tr1::variate_generator<rng_t&, tr1::uniform_real<> >
+                        sample(rng, tr1::uniform_real<>(0.0, _probs.back()));
+                    double r = sample();
+                    i = upper_bound(_probs.begin(), _probs.end(), r) -
+                        _probs.begin();
+                }
+                else
+                {
+                    cout << "all probabilities are zero... sample randomly.\n";
+                    // all probabilities are zero... sample randomly.
+                    tr1::uniform_int<size_t>
+                        sample(0, _candidates_set.size()-1);
+                    size_t j = sample(rng), count = 0;
+                    for (typeof(_candidates_set.begin()) iter =
+                             _candidates_set.begin();
+                         iter != _candidates_set.end(); ++iter)
+                    {
+                        if (count == j)
+                        {
+                            i = iter->second;
+                            break;
+                        }
+                        count++;
+                    }
+                }
+            } while (_erased[i]);
+
+            if (remove)
+            {
+                _erased[i] = true;
+                _erased_prob += (i > 0) ? _probs[i] - _probs[i-1] : _probs[i];
+                clean();
+            }
+
+            return _candidates[i];
+        }
+    }
+
+    void clean()
+    {
+        // if too many elements were erased, we need to make things less sparse
+        if (_biased && !_candidates_set.empty() &&
+            _erased_prob >= _probs.back()/3)
+        {
+            for (int i = int(_probs.size()) - 1; i > 0; --i)
+                _probs[i] -= _probs[i-1];
+
+            for (size_t i = 0; i < _candidates.size(); ++i)
+            {
+                while (i < _erased.size() && _erased[i])
+                {
+                    swap(_candidates[i], _candidates.back());
+                    _candidates.pop_back();
+
+                    swap(_probs[i], _probs.back());
+                    _probs.pop_back();
+
+                    swap(_erased[i], _erased.back());
+                    _erased.pop_back();
+                }
+            }
+
+            for (size_t i = 1; i < _probs.size(); i++)
+                _probs[i] += _probs[i-1];
+
+            _candidates_set.clear();
+            for (size_t i = 0; i < _candidates.size(); i++)
+                _candidates_set.insert(make_pair(_candidates[i],i));
+            _erased_prob = 0.0;
+        }
+    }
+
+private:
+    bool _biased;
+    vector<ValueType> _candidates;
+    tr1::unordered_multimap<ValueType, size_t, hash<ValueType> >
+        _candidates_set;
+    vector<double> _probs;
+    vector<uint8_t> _erased;
+    double _erased_prob;
+};
+
+} // namespace graph_tool
+
+#endif // SAMPLER_HH

src/graph_tool/generation/__init__.py

@@ -50,7 +50,7 @@ def _corr_wrap(i, j, corr):
     return corr(i[1], j[1])
 
 def random_graph(N, deg_sampler, deg_corr=None, directed=True,
-                 parallel=False, self_loops=False, verbose=False):
+                 parallel_edges=False, self_loops=False, verbose=False):
     r"""
     Generate a random graph, with a given degree distribution and correlation.
 
@@ -73,10 +73,12 @@ def random_graph(N, deg_sampler, deg_corr=None, directed=True,
         an edge existing in the generated graph.
     directed : bool (optional, default: True)
         Whether the generated graph should be directed.
-    parallel : bool (optional, default: False)
+    parallel_edges : bool (optional, default: False)
         If True, parallel edges are allowed.
     self_loops : bool (optional, default: False)
         If True, self-loops are allowed.
+    verbose : bool (optional, default: False)
+        If True, verbose information is displayed.
 
     Returns
     -------
@@ -224,25 +226,26 @@ def random_graph(N, deg_sampler, deg_corr=None, directed=True,
         uncorrelated = True
     else:
         uncorrelated = False
-    if not directed and deg_corr != None:
-        corr = lambda i,j: _corr_wrap(i, j, deg_corr)
-    else:
-        corr = deg_corr
-    libgraph_tool_generation.gen_random_graph(g._Graph__graph, N,
-                                              deg_sampler, corr,
-                                              uncorrelated, not parallel,
+    libgraph_tool_generation.gen_random_graph(g._Graph__graph, N, deg_sampler,
+                                              uncorrelated, not parallel_edges,
                                               not self_loops, not directed,
                                               seed, verbose)
     g.set_directed(directed)
+    random_rewire(g, parallel_edges = parallel_edges, self_loops = self_loops,
+                  verbose = verbose)
+    if deg_corr != None:
+        random_rewire(g, strat = "probabilistic",
+                      parallel_edges = parallel_edges, deg_corr = deg_corr,
+                      self_loops = self_loops, verbose = verbose)
     return g
 
-@_limit_args({"strat":["erdos", "correlated", "uncorrelated"]})
+@_limit_args({"strat":["erdos", "correlated", "uncorrelated", "probabilistic"]})
 def random_rewire(g, strat= "uncorrelated", parallel_edges = False,
-                  self_loops = False):
+                  self_loops = False, deg_corr = None, verbose = False):
     r"""
     Shuffle the graph in-place. If `strat` != "erdos", the degrees (either in or
     out) of each vertex are always the same, but otherwise the edges are
-    randomly placed. If `strat` == "correlated", the degree correlations are
+    randomly placed. If `strat` = "correlated", the degree correlations are
     also maintained: The new source and target of each edge both have the same
     in and out-degree. If `strat` = "probabilistic", than edges are rewired
     according to the degree correlation given by the parameter `deg_corr`.
@@ -256,10 +259,22 @@ def random_rewire(g, strat= "uncorrelated", parallel_edges = False,
         `strat` == "uncorrelated" only the degrees of the vertices will be
         maintained, nothing else. If `strat` == "correlated", additionally the
         new source and target of each edge both have the same in and out-degree.
+        If `strat` == "probabilistic", than edges are rewired according to the
+        degree correlation given by the parameter `deg_corr`.
     parallel : bool (optional, default: False)
         If True, parallel edges are allowed.
     self_loops : bool (optional, default: False)
         If True, self-loops are allowed.
+    deg_corr : function (optional, default: None)
+        A function which gives the degree correlation of the graph. It should be
+        callable with two parameters: the in,out-degree pair of the source
+        vertex an edge, and the in,out-degree pair of the target of the same
+        edge (for undirected graphs, both parameters are single values). The
+        function should return a number proportional to the probability of such
+        an edge existing in the generated graph. This parameter is ignored,
+        unless `strat` = "probabilistic".
+    verbose : bool (optional, default: False)
+        If True, verbose information is displayed.
 
     See Also
     --------
@@ -405,10 +420,16 @@ def random_rewire(g, strat= "uncorrelated", parallel_edges = False,
 
     seed = numpy.random.randint(0, sys.maxint)
 
+    if not g.is_directed() and deg_corr != None:
+        corr = lambda i,j: _corr_wrap(i, j, deg_corr)
+    else:
+        corr = deg_corr
+
     g.stash_filter(reversed=True)
     try:
         libgraph_tool_generation.random_rewire(g._Graph__graph, strat,
-                                               self_loops, parallel_edges, seed)
+                                               self_loops, parallel_edges,
+                                               corr, seed, verbose)
     finally:
         g.pop_filter(reversed=True)