Modify absolute_centrality implementation

This dumps the random-path approach, and implements a deterministic,
priority queue algorithm.

src/graph/centrality/graph_absolute_trust.cc

@@ -22,9 +22,6 @@
 #include "graph.hh"
 #include "graph_selectors.hh"
 
-#include <tr1/random>
-typedef std::tr1::mt19937 rng_t;
-
 #include "graph_absolute_trust.hh"
 
 using namespace std;
@@ -32,10 +29,8 @@ using namespace boost;
 using namespace graph_tool;
 
 void absolute_trust(GraphInterface& g, int64_t source, boost::any c,
-                    boost::any t, size_t iter, bool reversed, size_t seed)
+                    boost::any t, size_t n_paths, double epsilon, bool reversed)
 {
-    rng_t rng(static_cast<rng_t::result_type>(seed));
-
     if (!belongs<edge_floating_properties>()(c))
         throw ValueException("edge property must be of floating point value type");
     if (!belongs<vertex_floating_vector_properties>()(t))
@@ -43,8 +38,7 @@ void absolute_trust(GraphInterface& g, int64_t source, boost::any c,
 
     run_action<>()(g,
                    bind<void>(get_absolute_trust(), _1, g.GetVertexIndex(),
-                              source, _2, _3, iter, reversed,
-                              ref(rng)),
+                              source, _2, _3, n_paths, epsilon, reversed),
                    edge_floating_properties(),
                    vertex_floating_vector_properties())(c, t);
 }

src/graph/centrality/graph_absolute_trust.hh

@@ -23,8 +23,10 @@
 #include "graph_util.hh"
 
 #include <tr1/unordered_set>
-#include <tr1/random>
-#include <boost/functional/hash.hpp>
+#include <tr1/tuple>
+#include <algorithm>
+
+#include "minmax.hh"
 
 #include <iostream>
 
@@ -32,183 +34,151 @@ namespace graph_tool
 {
 using namespace std;
 using namespace boost;
+using std::tr1::get;
+using std::tr1::tuple;
+
+template <class Path>
+struct path_cmp
+{
+    path_cmp(vector<Path>& paths): _paths(paths) {}
+    vector<Path>& _paths;
+
+    typedef size_t first_argument_type;
+    typedef size_t second_argument_type;
+    typedef bool result_type;
+    inline bool operator()(size_t a, size_t b)
+    {
+        if (get<0>(_paths[a]).first == get<0>(_paths[b]).first)
+            return get<1>(_paths[a]).size() > get<1>(_paths[b]).size();
+        return get<0>(_paths[a]).first < get<0>(_paths[b]).first;
+    }
+};
 
 struct get_absolute_trust
 {
     template <class Graph, class VertexIndex, class TrustMap,
               class InferredTrustMap>
     void operator()(Graph& g, VertexIndex vertex_index, int64_t source,
-                    TrustMap c, InferredTrustMap t, size_t n_iter,
-                    bool reversed, rng_t& rng) const
+                    TrustMap c, InferredTrustMap t, size_t n_paths,
+                    double epsilon, bool reversed) const
     {
+        typedef typename graph_traits<Graph>::vertex_descriptor vertex_t;
         typedef typename property_traits<TrustMap>::value_type c_type;
         typedef typename property_traits<InferredTrustMap>::value_type
             ::value_type t_type;
 
-        unchecked_vector_property_map<vector<t_type>, VertexIndex>
-            t_count(vertex_index, num_vertices(g));
-        unchecked_vector_property_map<vector<size_t>, VertexIndex>
-            v_mark(vertex_index, num_vertices(g));
+        typedef tuple<pair<t_type, t_type>, tr1::unordered_set<vertex_t>,
+                      size_t > path_t;
 
-        // init inferred trust t
+        double delta = epsilon+1;
         int i, N = num_vertices(g);
         #pragma omp parallel for default(shared) private(i) schedule(dynamic)
         for (i = (source == -1) ? 0 : source;
              i < ((source == -1) ? N : source + 1); ++i)
         {
-            typename graph_traits<Graph>::vertex_descriptor v =
-                vertex(i, g);
+            vertex_t v = vertex(i, g);
             if (v == graph_traits<Graph>::null_vertex())
                 continue;
-            t[v].resize(N);
-            t_count[v].resize(N,0);
-            v_mark[v].resize(N,0);
-        }
 
-        #pragma omp parallel for default(shared) private(i) schedule(dynamic)
-        for (i = (source == -1) ? 0 : source;
-             i < ((source == -1) ? N : source + 1); ++i)
-        {
-            typename graph_traits<Graph>::vertex_descriptor v =
-                vertex(i, g);
-            if (v == graph_traits<Graph>::null_vertex())
-                continue;
+            // path priority queue
+            vector<path_t> paths(1);
+            vector<size_t> free_indexes;
+            typedef double_priority_queue<size_t, path_cmp<path_t> > queue_t;
+            queue_t queue = queue_t(path_cmp<path_t>(paths));
+
+            get<0>(paths.back()).first = get<0>(paths.back()).second = 1;
+            get<1>(paths.back()).insert(v);
+            get<2>(paths.back()) = vertex_index[v];
+            queue.push(0);
 
-            // walk hash set
-            tr1::unordered_set<size_t> path_set;
-            tr1::uniform_int<size_t>
-                random_salt(0, numeric_limits<size_t>::max());
-            size_t salt;
+            t[v].resize(num_vertices(g));
+            unchecked_vector_property_map<t_type, VertexIndex>
+                sum_weight(vertex_index, num_vertices(g));
+
+            size_t count = 1;
+            while (!queue.empty())
             {
-                #pragma omp critical
-                salt = random_salt(rng);
-            }
+                size_t pi = queue.top();
+                queue.pop_top();
+                vertex_t w = vertex(get<2>(paths[pi]), g);
 
-            for (size_t bias = 0; bias < 3; ++bias)
-                for (size_t iter = 0; iter < n_iter; ++iter)
+                typename graph_traits<Graph>::out_edge_iterator e, e_end;
+                for (tie(e, e_end) = out_edges(w, g); e != e_end; ++e)
                 {
-                    typename graph_traits<Graph>::vertex_descriptor pos = v;
-                    t_type pos_t = 1.0;
-                    t_type pweight = 1.0;
-                    v_mark[v][vertex_index[v]] = bias*n_iter + iter + 1;
-
-                    size_t path_hash = salt;
-                    hash_combine(path_hash, i);
-
-                    // start a self-avoiding walk from vertex v
-                    vector<typename graph_traits<Graph>::edge_descriptor>
-                        out_es;
-                    vector<t_type> out_prob;
-                    do
+                    vertex_t a = target(*e, g);
+                    // no loops
+                    if (get<1>(paths[pi]).find(a) == get<1>(paths[pi]).end())
                     {
-                        out_es.clear();
-                        out_prob.clear();
-
-                        // obtain list of candidate edges to follow
-                        typename graph_traits<Graph>::out_edge_iterator e,
-                            e_end;
-                        for (tie(e, e_end) = out_edges(pos, g); e != e_end; ++e)
+                        // new path;  only follow non-zero paths
+                        t_type old = (sum_weight[a] > 0) ?
+                            t[v][a]/sum_weight[a] : 0;
+                        if (c[*e] > 0 || (reversed &&
+                                          get<1>(paths[pi]).size() == 1))
                         {
-                            typename graph_traits<Graph>::vertex_descriptor t =
-                                target(*e,g);
-                            if (v_mark[v][vertex_index[t]] <
-                                bias*n_iter + iter + 1)
+                            size_t npi;
+                            if (free_indexes.empty())
                             {
-                                if (c[*e] < 0 || c[*e] > 1)
-                                    continue;
-                                out_es.push_back(*e);
-
-                                if (bias != 1)
-                                {
-                                    double prob = c[*e];
-                                    if (bias == 2)
-                                        prob = 1.0 - prob;
-
-                                    if (out_prob.empty())
-                                        out_prob.push_back(prob);
-                                    else
-                                        out_prob.push_back(out_prob.back() +
-                                                           prob);
-                                }
+                                paths.push_back(paths[pi]); // clone last path
+                                npi = paths.size()-1;
+                            }
+                            else
+                            {
+                                npi = free_indexes.back();
+                                free_indexes.pop_back();
+                                paths[npi] = paths[pi];
                             }
-                        }
 
-                        if (!out_es.empty())
-                        {
-                            // select edge according to its probability
-                            typename graph_traits<Graph>::edge_descriptor e;
-                            if (!out_prob.empty() && out_prob.back() > 0)
+                            path_t& np = paths[npi];
+                            if (!reversed)
                             {
-                                typedef tr1::uniform_real<t_type> dist_t;
-                                tr1::variate_generator<rng_t&, dist_t>
-                                    random(rng, dist_t(0, out_prob.back()));
-
-                                t_type u;
-                                {
-                                    #pragma omp critical
-                                    u = random();
-                                }
-                                e = out_es[upper_bound(out_prob.begin(),
-                                                       out_prob.end(), u) -
-                                           out_prob.begin()];
+                                get<0>(np).second = get<0>(np).first;
+                                get<0>(np).first *= c[*e];
                             }
                             else
                             {
-                                tr1::uniform_int<size_t>
-                                    random(0, out_es.size()-1);
-                                #pragma omp critical
-                                e = out_es[random(rng)];
+                                if (get<1>(np).size() > 1)
+                                    get<0>(np).second *= c[*e];
+                                get<0>(np).first *= c[*e];
                             }
+                            get<1>(np).insert(a);
+                            get<2>(np) = vertex_index[a];
 
-                            // new position in random walk
-                            pos = target(e,g);
-                            size_t posi = vertex_index[pos];
-
-                            // update current path trust
-                            pos_t *= c[e];
-
-                            if (reversed && boost::source(e, g) != v)
-                                pweight *= c[e];
-
-                            // get path hash
-                            hash_combine(path_hash, posi);
-                            if (path_set.find(path_hash) == path_set.end())
+                            t[v][a] += get<0>(np).first*get<0>(np).second;
+                            sum_weight[a] += get<0>(np).second;
+                            queue.push(npi);
+                            if (n_paths > 0 && queue.size() > n_paths)
                             {
-                                // if new path, modify vertex trust score
-                                path_set.insert(path_hash);
-
-                                t_type old = 0;
-                                if (t_count[v][posi] > 0)
-                                    old = t[v][posi]/t_count[v][posi];
-                                t[v][posi] += pos_t*pweight;
-                                t_count[v][posi] += pweight;
+                                size_t bi = queue.bottom();
+                                queue.pop_bottom();
+                                free_indexes.push_back(bi);
                             }
-
-                            if (!reversed)
-                                pweight *= c[e];
-
-                            // stop if weight drops to zero
-                            if (pweight == 0)
-                                break;
-
-                            // mark vertex
-                            v_mark[v][posi] = bias*n_iter + iter + 1;
                         }
+                        else
+                        {
+                            sum_weight[a] +=  get<0>(paths[pi]).first;
+                        }
+                        if (sum_weight[a] > 0)
+                            delta += abs(old-t[v][a]/sum_weight[a]);
                     }
-                    while (!out_es.empty());
                 }
-        }
+                free_indexes.push_back(pi);
 
-        #pragma omp parallel for default(shared) private(i) schedule(dynamic)
-        for (i = (source == -1) ? 0 : source;
-             i < ((source == -1) ? N : source + 1); ++i)
-        {
-            typename graph_traits<Graph>::vertex_descriptor v = vertex(i, g);
-            if (v == graph_traits<Graph>::null_vertex())
-                continue;
-            for (size_t j = 0; j < size_t(N); ++j)
-                if (t_count[v][j] > 0)
-                    t[v][j] /= t_count[v][j];
+                if ((count % N) == 0)
+                {
+                    if (delta < epsilon)
+                        break;
+                    else
+                        delta = 0;
+                }
+                count++;
+            }
+
+            typename graph_traits<Graph>::vertex_iterator w, w_end;
+            for (tie(w, w_end) = vertices(g); w != w_end; ++w)
+            {
+                if (sum_weight[*w] > 0)
+                    t[v][*w] /= sum_weight[*w];
+            }
         }
     }
 };

src/graph/centrality/minmax.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 MINMAX_HH
+#define MINMAX_HH
+
+#include <vector>
+#include <boost/mpl/bool.hpp>
+#include <iostream>
+using namespace  std;
+
+template <class Value, class Compare = std::less<Value> >
+class double_priority_queue
+{
+public:
+    double_priority_queue() {}
+    double_priority_queue(const Compare& cmp): _cmp(cmp) {}
+
+    size_t size() { return _queue.size(); }
+    bool empty() { return _queue.empty(); }
+
+    void push(const Value& v)
+    {
+        _queue.push_back(v);
+        bubble_up(_queue.size()-1);
+    }
+
+    void pop_top()
+    {
+        // for (size_t i = 0; i < _queue.size(); ++i)
+        //     if (_cmp(top(), _queue[i]))
+        //     {
+        //         cout << size_t(log2(i+1)) << endl;
+        //         cout << "ops top" << endl;
+        //         abort();
+        //     }
+        size_t i = &top() - &_queue[0];
+        swap(_queue[i], _queue.back());
+        _queue.pop_back();
+        if (!_queue.empty())
+            trickle_down(i);
+    }
+
+    void pop_bottom()
+    {
+        return;
+        swap(_queue.front(), _queue.back());
+        _queue.pop_back();
+        if (!_queue.empty())
+            trickle_down(0);
+    }
+
+    const Value& top()
+    {
+        if (_queue.size() < 2)
+            return _queue[0];
+        if (_queue.size() < 3)
+            return _queue[1];
+        if (_cmp(_queue[1], _queue[2]))
+            return _queue[2];
+        else
+            return _queue[1];
+    }
+
+    const Value& bottom()
+    {
+        // for (size_t i = 0; i < _queue.size(); ++i)
+        //     if (_cmp(_queue[i], _queue.front()))
+        //         cout << "ops bottom" << endl;
+        return _queue.front();
+    }
+
+    // general comparison: if IsMin is mpl::true_, then a < b, otherwise a > b
+    // is computed
+    template <class IsMin>
+    bool cmp(size_t a, size_t b, IsMin)
+    {
+        if (IsMin::value)
+            return _cmp(_queue[a], _queue[b]);
+        else
+            return _cmp(_queue[b], _queue[a]);
+    }
+
+    void trickle_down(size_t i)
+    {
+        if ((size_t(log2(i+1)) % 2) == 0)
+            trickle_down(i, boost::mpl::true_());
+        else
+            trickle_down(i, boost::mpl::false_());
+    }
+
+
+    template <class IsMin>
+    void trickle_down(size_t pos, IsMin is_min)
+    {
+        size_t i = pos;
+        while (true)
+        {
+            size_t m = 2*i+1;
+            if (m >= _queue.size())
+                break;
+
+            for (size_t j = 2*i+1; j < 2*i+3; ++j)
+            {
+                if (j >= _queue.size())
+                    break;
+                if (cmp(j, m, is_min))
+                    m = j;
+                for (size_t l = 2*j+1; l < 2*j+3; ++l)
+                {
+                    if (l >= _queue.size())
+                        break;
+                    if (cmp(l, m, is_min))
+                        m = l;
+                }
+            }
+
+            if (m > 2*i+2) // is grandchild
+            {
+                if (cmp(m, i, is_min))
+                {
+                    swap(_queue[m], _queue[i]);
+                    if (!cmp(m, (m-1)/2, is_min))
+                        swap(_queue[m], _queue[(m-1)/2]);
+                    i = m;
+                }
+                else
+                {
+                    break;
+                }
+            }
+            else
+            {
+                if (cmp(m, i, is_min))
+                    swap(_queue[m], _queue[i]);
+                break;
+            }
+        }
+    }
+
+    void bubble_up(size_t i)
+    {
+        if ((size_t(log2(i+1)) % 2) == 0)
+        {
+            if ((i > 0) && cmp(i, (i-1)/2, boost::mpl::false_()))
+            {
+                swap(_queue[i], _queue[(i-1)/2]);
+                bubble_up((i-1)/2,  boost::mpl::false_());
+            }
+            else
+            {
+                bubble_up(i,  boost::mpl::true_());
+            }
+        }
+        else
+        {
+            if ((i > 0) && cmp(i, (i-1)/2, boost::mpl::true_()))
+            {
+                swap(_queue[i], _queue[(i-1)/2]);
+                bubble_up((i-1)/2, boost::mpl::true_());
+            }
+            else
+            {
+                bubble_up(i, boost::mpl::false_());
+            }
+        }
+    }
+
+    template <class IsMin>
+    void bubble_up(size_t pos, IsMin is_min)
+    {
+        size_t i = pos;
+        while (true)
+        {
+            if (i > 2) // has grandparent
+            {
+                size_t m = ((i-1)/2 - 1)/2;
+                if (cmp(i, m, is_min))
+                {
+                    swap(_queue[m], _queue[i]);
+                    i = m;
+                }
+                else
+                {
+                    break;
+                }
+            }
+            else
+            {
+                break;
+            }
+        }
+    }
+
+private:
+    std::vector<Value> _queue;
+    Compare _cmp;
+};
+
+
+#endif // MINMAX_HH

src/graph_tool/centrality/__init__.py

@@ -390,10 +390,10 @@ def eigentrust(g, trust_map, vprop=None, norm=False, epslon=1e-6, max_iter=0,
     else:
         return vprop
 
-def absolute_trust(g, trust_map, source=None, vprop=None, n_iter=100,
-                   reversed=False, seed=None):
+def absolute_trust(g, trust_map, source=None, vprop=None, epsilon=1e-6,
+                   n_paths=100000, reversed=False):
     r"""
-    Samples the absolute trust centrality of each vertex in the graph, or only
+    Calculate the absolute trust centrality of each vertex in the graph, or only
     for a given source, if one is provided.
 
     Parameters
@@ -408,16 +408,15 @@ def absolute_trust(g, trust_map, source=None, vprop=None, n_iter=100,
         values, instead of all the vertices in the graph.
     vprop : PropertyMap, optional (default: None)
         Vertex property map where the values of eigentrust must be stored.
-    n_iter : int, optional (default: 100)
-        Number of iterations (independent self-avoiding walks) per source
-        vertex.
+    epsilon : float, optional (default: 1e-6)
+        Convergence criterion for the algorithm. If the summed difference of
+        trust values drops below this value, the algorithm stops.
+    n_paths : int, optimal (default: 100000)
+        Maximum number of paths to keep following at any given time.
     reversed : bool, optional (default: False)
         Calculates the "reversed" trust instead: The direction of the edges are
         inverted, but the path weighting is preserved in the original direction
         (see Notes below).
-    seed : int, optional (default: None)
-         The initializing seed for the random number generator. If not supplied
-         a different random value will be chosen each time.
 
     Returns
     -------
@@ -455,14 +454,11 @@ def absolute_trust(g, trust_map, source=None, vprop=None, n_iter=100,
 
     such that the direct trust of the last edge on the path is not considered.
 
-    The algorithm performs only an approximation of the above measure, by doing
-    several self-avoiding random walks per source vertex, and computing the
-    trust for all different paths found. Each complete walk is done by one of
-    three types of biased choices, which govern the probability of following a
-    specific edge: 1. With probability proportional to an edge's trust value;
-    2. With probability proportional to the complement of an edge's trust value
-    :math:`(1-t_{e})`; 3. All edges with equal probability. The parameter
-    "n_iter" controls how many walks of each type are performed.