Improve modularity() computation

src/graph/community/graph_community.cc

@@ -88,7 +88,7 @@ double modularity(GraphInterface& g, boost::any weight, boost::any property)
 
     run_action<graph_tool::detail::never_directed>()
         (g, bind<void>(get_modularity(), _1, _2, _3, ref(modularity)),
-         edge_props_t(), vertex_properties())
+         edge_props_t(), vertex_scalar_properties())
         (weight, property);
     return modularity;
 }

src/graph/community/graph_community.hh

@@ -489,35 +489,50 @@ struct get_communities_selector
 struct get_modularity
 {
     template <class Graph, class WeightMap, class CommunityMap>
-    void operator()(const Graph& g, WeightMap weights, CommunityMap s,
-                    double& modularity) const
+    void operator()(const Graph& g, WeightMap weights, CommunityMap b, double& Q) const
     {
         typedef typename property_traits<WeightMap>::key_type weight_key_t;
         typedef typename property_traits<WeightMap>::value_type weight_val_t;
         typedef typename property_traits<CommunityMap>::value_type s_val_t;
 
-        modularity = 0.0;
+        vector<double> er, err;
         double W = 0;
 
         typename graph_traits<Graph>::edge_iterator e, e_end;
-        for (tie(e,e_end) = edges(g); e != e_end; ++e)
-            if (target(*e,g) != source(*e,g))
+        for (tie(e, e_end) = edges(g); e != e_end; ++e)
+        {
+            size_t r = get(b, source(*e,g));
+            size_t s = get(b, target(*e,g));
+
+            double w = get(weights, *e);
+            W += 2 * w;
+
+            if (er.size() <= r)
+                er.resize(r + 1);
+            er[r] += w;
+
+            if (er.size() <= s)
+                er.resize(s + 1);
+            er[s] += w;
+
+            if (r == s)
             {
-                W += get(weights, *e);
-                if (get(s, target(*e,g)) == get(s, source(*e,g)))
-                    modularity += 2 * get(weights, weight_key_t(*e));
+                if (err.size() <= r)
+                    err.resize(r + 1);
+                err[r] += 2 * w;
             }
+        }
 
-        unordered_map<s_val_t, weight_val_t> Ks;
-
-        typename graph_traits<Graph>::vertex_iterator v, v_end;
-        for (tie(v,v_end) = vertices(g); v != v_end; ++v)
-            Ks[get(s, *v)] += out_degree_no_loops_weighted(*v, weights, g);
+        Q = 0;
+        for (size_t r = 0; r < er.size(); ++r)
+        {
+            if (err.size() <= r)
+                err.resize(r + 1);
 
-        for (typeof(Ks.begin()) iter = Ks.begin(); iter != Ks.end(); ++iter)
-            modularity -= (iter->second*iter->second)/double(2*W);
+            Q += err[r] - (er[r] * er[r]) / W;
+        }
+        Q /= W;
 
-        modularity /= 2*W;
     }
 };
 

src/graph_tool/community/__init__.py

@@ -305,12 +305,15 @@ def modularity(g, prop, weight=None):
 
     .. math::
 
-          Q = \sum_s e_{ss}-\left(\sum_r e_{rs}\right)^2
+          Q = \frac{1}{2E} \sum_r e_{rr}- \frac{e_r^2}{2E}
 
-    where :math:`e_{rs}` is the fraction of edges which fall between
-    vertices with spin s and r.
+    where :math:`e_{rs}` is the number of edges which fall between
+    vertices in communities s and r, or twice that number if :math:`r = s`, and
+    :math:`e_r = \sum_s e_{rs}`.
 
-    If enabled during compilation, this algorithm runs in parallel.
+    If weights are provided, the matrix :math:`e_{rs}` corresponds to the sum
+    of edge weights instead of number of edges, and the value of :math:`E`
+    becomes the total sum of edge weights.
 
     Examples
     --------
@@ -318,9 +321,9 @@ def modularity(g, prop, weight=None):
     >>> from numpy.random import seed
     >>> seed(42)
     >>> g = gt.load_graph("community.xml")
-    >>> spins = gt.community_structure(g, 10000, 10)
-    >>> gt.modularity(g, spins)
-    0.535314188562404
+    >>> b = gt.community_structure(g, 10000, 10)
+    >>> gt.modularity(g, b)
+    0.5353141885624041
 
     References
     ----------