* correct graphml headers

* fix ticket [9]
* fix --for bug with multiple options with different parameters
* Add initial support for community detection (ticket [6])


git-svn-id: https://svn.forked.de/graph-tool/trunk@47 d4600afd-f417-0410-95de-beed9576f240

src/graph-tool

@@ -125,6 +125,10 @@ layout =  parser.add_option_group("Layout")
 layout.add_option("--compute-spring-block-layout", action="callback", callback=push_option, type="string", metavar="ITERATIONS[|SEED]", help="compute the spring block layout")
 layout.add_option("--compute-gursoy-atun-layout", action="callback", callback=push_option, type="string", metavar="ITERATIONS[|SEED]", help="compute the Gursoy-Atun layout")
 
+community =  parser.add_option_group("Community")
+clustering.add_option("--community-structure", action="callback", callback=push_option, type="string", metavar="PROPERTY|CORR_TYPE|GAMMA|N_ITER|SEED[|WEIGHT]", help="calculate the community structure and assign it to PROPERTY. The value of CORR_TYPE can be 'random', 'uncorrelated' or 'correlated'.")
+clustering.add_option("--get-modularity", action="callback", callback=push_option, type="string", metavar="PROPERTY[|WEIGHT]|FILE", help="calculate the modularity, given a community partition specified by PROPERTY")
+
 layout =  parser.add_option_group("History")
 layout.add_option("--for", action="callback", callback=push_option, type="string", metavar="INIT|CONDITION|STEP", help="simplified scripting")
 layout.add_option("--history", action="callback", callback=push_option, type="string", metavar="INIT|CONDITION|STEP", help="simplified scripting (does not overwrite previous results)")
@@ -285,15 +289,27 @@ def inv_corr_ceil(p,r,j,k):
 
 
 def parse_option(opt, just_file=False):
-    "this will execute an aption, and return either None, or a tuple with the result and the respective file name, if it exists"
+    "this will execute an option, and return either None, or a tuple with the result and the respective file name, if it exists"
     if opt.name == "load":
+	values = parse_values(opt.value)
+	if len(values) > 2 or len(values) < 1:
+	    raise OptionError(opt.name, "invalid value '$s'" % opt.value)
         if just_file:
             return None
-        graph.ReadFromFile(opt.value)
+	if len(values) == 1:
+	    graph.ReadFromFile(values[0])
+        else:
+            graph.ReadFromFile(values[0], values[1])
     elif opt.name == "save":
+	values = parse_values(opt.value)
+        if len(values) > 2 or len(values) < 1:
+	    raise OptionError(opt.name, "invalid value '$s'" % opt.value)
         if just_file:
             return None
-        graph.WriteToFile(opt.value)
+	if len(values) == 1:
+	    graph.WriteToFile(values[0])
+        else:
+            graph.WriteToFile(values[0], values[1])
     elif opt.name == "correlated-configurational-model":
         if just_file:
             return None
@@ -501,6 +517,46 @@ def parse_option(opt, just_file=False):
         if len(values) == 0 or len(values) > 2:
             raise OptionError(opt.name, "invalid value '%s'" % opt.value)
         graph.SetExtendedClusteringToProperty(values[0],int(values[1]))
+    elif opt.name == "community-structure":
+        if just_file:
+            return None
+        values = parse_values(opt.value)
+        if len(values) < 5 or len(values) > 6:
+            raise OptionError(opt.name, "invalid value '%s'" % opt.value)
+        offset = 0;
+        prop = values[0]
+        corr = None
+        if values[1].strip() == "random":
+            corr = CommCorr.ErdosReyni
+        elif values[1].strip() == "uncorrelated":
+            corr = CommCorr.Uncorrelated
+        elif values[1].strip() == "correlated":
+            corr = CommCorr.Correlated
+        else:
+            raise OptionError(opt.name, "invalid correlation type value '%s'" % values[1])
+        try:
+            g = float(values[2])
+            N = int(values[3])
+            seed = int(values[4])
+            if seed == 0:
+                seed = int(time())
+            weights = ""
+            if len(values) == 6:
+                weights = values[5]
+        except ValueError,e:
+            raise OptionError(opt.name, "invalid value '%s': %s" % (opt.value, e))
+        graph.GetCommunityStructure(g, corr, N, seed, weights, prop)
+    elif opt.name == "get-modularity":
+        values = parse_values(opt.value)
+        if len(values) < 2 or len(values) > 3:
+            raise OptionError(opt.name, "invalid value '%s'" % opt.value)
+        prop, weight, file_name = values[0], "", values[1]
+        if len(values) > 2:
+            weight = values[1]
+            file_name = values[2]
+        if just_file:
+            return file_name
+        return (graph.GetModularity(weight, prop), file_name)
     elif opt.name == "compute-spring-block-layout":
         if just_file:
             return None
@@ -783,14 +839,15 @@ try:
                     prefix = str(count)
                     if variables.has_key("prefix"):
                         prefix = str(variables["prefix"])
-                    if not history_files.has_key(opt.name):
-                        history_files[opt.name] = HistoryFile(file_name, overwrite_history)
+                    hist_key = opt.name + opt.value
+                    if not history_files.has_key(hist_key):
+                        history_files[hist_key] = HistoryFile(file_name, overwrite_history)
                     else:
-                        if history_files[opt.name].file_name != os.path.expanduser(file_name):
-                            history_files[opt.name] = HistoryFile(file_name, overwrite_history)
+                        if history_files[hist_key].file_name != os.path.expanduser(file_name):
+                            history_files[hist_key] = HistoryFile(file_name, overwrite_history)
                     try:
-                        if history_files[opt.name].last_prefix != None:
-                            is_new = float(prefix) > float(history_files[opt.name].last_prefix)
+                        if history_files[hist_key].last_prefix != None:
+                            is_new = float(prefix) > float(history_files[hist_key].last_prefix)
                         else:
                             is_new = True
                     except ValueError:
@@ -798,9 +855,9 @@ try:
                     if is_new or overwrite_history:
                         data = parse_option(new_opt)[0]
                         if uses_prefix:
-                            history_files[opt.name].write(data, prefix)
+                            history_files[hist_key].write(data, prefix)
                         else:
-                            history_files[opt.name].write(data, "")
+                            history_files[hist_key].write(data, "")
                 else:
                     parse_option(new_opt)
             exec values[2] in variables # step

src/graph/Makefile.am

@@ -36,6 +36,7 @@ libgraph_tool_la_SOURCES = \
     graph_distance_sampled.cc\
     graph_reciprocity.cc\
     graph_minimum_spanning_tree.cc\
+    graph_community.cc\
     graph_io.cc\
     graph_bind.cc\
     graphml.hpp\

src/graph/graph.hh

@@ -96,6 +96,17 @@ public:
     double GetReciprocity() const;
     void   GetMinimumSpanningTree(std::string weight, std::string property);
 
+    // community structure
+    enum comm_corr_t
+    {
+	ERDOS_REYNI,
+	UNCORRELATED,
+	CORRELATED
+    };
+
+    void   GetCommunityStructure(double gamma, comm_corr_t corr, size_t n_iter, size_t seed, std::string weight, std::string property);
+    double GetModularity(std::string weight, std::string property);
+
     // filtering
     void SetDirected(bool directed) {_directed = directed;}
     bool GetDirected() const {return _directed;}
@@ -132,7 +143,9 @@ public:
 
     // i/o
     void WriteToFile(std::string s);
+    void WriteToFile(std::string s, std::string format);
     void ReadFromFile(std::string s);
+    void ReadFromFile(std::string s, std::string format);
 
     // signal handling
     void InitSignalHandling();

src/graph/graph_bind.cc

@@ -207,6 +207,13 @@ struct LibInfo
     string GetVersion()   const {return VERSION;}
 };
 
+// overloads
+void  (GraphInterfaceWrap::*ReadFromFile1) (string)          = &GraphInterfaceWrap::ReadFromFile;
+void  (GraphInterfaceWrap::*ReadFromFile2) (string, string)  = &GraphInterfaceWrap::ReadFromFile;
+void  (GraphInterfaceWrap::*WriteToFile1)  (string)          = &GraphInterfaceWrap::WriteToFile;
+void  (GraphInterfaceWrap::*WriteToFile2)  (string, string)  = &GraphInterfaceWrap::WriteToFile;
+
+
 BOOST_PYTHON_MODULE(libgraph_tool)
 {
     class_<GraphInterfaceWrap>("GraphInterface")
@@ -230,6 +237,8 @@ BOOST_PYTHON_MODULE(libgraph_tool)
 	.def("GetSampledDistanceHistogram", &GraphInterfaceWrap::GetSampledDistanceHistogram)
 	.def("GetReciprocity", &GraphInterfaceWrap::GetReciprocity)
 	.def("GetMinimumSpanningTree", &GraphInterfaceWrap::GetMinimumSpanningTree)
+	.def("GetCommunityStructure", &GraphInterfaceWrap::GetCommunityStructure)
+	.def("GetModularity", &GraphInterfaceWrap::GetModularity)
 	.def("SetDirected", &GraphInterfaceWrap::SetDirected)
 	.def("GetDirected", &GraphInterfaceWrap::GetDirected)
 	.def("SetReversed", &GraphInterfaceWrap::SetReversed)
@@ -254,8 +263,10 @@ BOOST_PYTHON_MODULE(libgraph_tool)
 	.def("InsertVertexIndexProperty",  &GraphInterfaceWrap::InsertVertexIndexProperty)
 	.def("ComputeGraphLayoutGursoy", &GraphInterfaceWrap::ComputeGraphLayoutGursoy)
 	.def("ComputeGraphLayoutSpringBlock", &GraphInterfaceWrap::ComputeGraphLayoutSpringBlock)
-	.def("WriteToFile", &GraphInterfaceWrap::WriteToFile)
-	.def("ReadFromFile", &GraphInterfaceWrap::ReadFromFile)
+	.def("WriteToFile", WriteToFile1)
+	.def("WriteToFile", WriteToFile2)
+	.def("ReadFromFile", ReadFromFile1)
+	.def("ReadFromFile", ReadFromFile2)
 	.def("InitSignalHandling", &GraphInterfaceWrap::InitSignalHandling);
 	
     enum_<GraphInterfaceWrap::degree_t>("Degree")
@@ -263,6 +274,11 @@ BOOST_PYTHON_MODULE(libgraph_tool)
 	.value("Out", GraphInterfaceWrap::OUT_DEGREE)
 	.value("Total", GraphInterfaceWrap::TOTAL_DEGREE);
 
+    enum_<GraphInterfaceWrap::comm_corr_t>("CommCorr")
+	.value("ErdosReyni", GraphInterfaceWrap::ERDOS_REYNI)
+	.value("Uncorrelated", GraphInterfaceWrap::UNCORRELATED)
+	.value("Correlated", GraphInterfaceWrap::CORRELATED);
+
     variant_from_python<string>();
     variant_from_python<GraphInterfaceWrap::degree_t>();
     to_python_converter<pair<double,double>, pair_to_tuple<double,double> >();

src/graph/graph_community.cc

+// graph-tool -- a general graph modification and manipulation thingy
+//
+// Copyright (C) 2006  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 2
+// 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, write to the Free Software
+// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
+
+#include <boost/lambda/lambda.hpp>
+#include <boost/lambda/bind.hpp>
+#include <boost/random.hpp>
+#include <tr1/unordered_set>
+
+#include "graph.hh"
+#include "histogram.hh"
+#include "graph_filtering.hh"
+#include "graph_selectors.hh"
+#include "graph_properties.hh"
+
+using namespace std;
+using namespace boost;
+using namespace boost::lambda;
+using namespace graph_tool;
+
+using std::tr1::unordered_map;
+using std::tr1::unordered_set;
+
+typedef boost::mt19937 rng_t;
+
+template <class Graph>
+size_t out_degree_no_loops(typename graph_traits<Graph>::vertex_descriptor v, const Graph &g)
+{
+    size_t k = 0;
+    typename graph_traits<Graph>::adjacency_iterator a,a_end;
+    for (tie(a,a_end) = adjacent_vertices(v,g); a != a_end; ++a)
+	if (*a != v)
+	    k++;
+    return k;
+}
+
+//==============================================================================
+// GetCommunityStructure()
+// computes the community structure through a spin glass system with 
+// simulated annealing
+//==============================================================================
+
+template <template <class G> class NNKS>
+struct get_communities
+{
+    template <class Graph, class WeightMap, class CommunityMap>
+    void operator()(Graph& g, WeightMap weights, CommunityMap s, double gamma, size_t n_iter, size_t seed) const
+    {
+	rng_t rng(seed);
+	boost::uniform_real<double> uniform_p(0.0,1.0);
+
+	unordered_map<size_t, size_t> Nk; // degree histogram
+	unordered_map<size_t, size_t> Ns; // spin histogram
+	unordered_map<size_t, unordered_map<size_t, size_t> > kNs; // spin histogram per degree
+	unordered_map<size_t, size_t> Ks; // sum of degrees per spin
+	unordered_map<size_t, map<double, unordered_set<size_t> > > global_term; // global energy term
+
+	NNKS<Graph> Nnnks(g); // this will retrieve the expected number of neighbours with given spin, in funcion of degree
+
+	// init spins from [0,N-1] and global info	
+	size_t index = 0;
+	typename graph_traits<Graph>::vertex_iterator v,v_end;
+	for (tie(v,v_end) = vertices(g); v != v_end; ++v)
+	{
+	    s[*v] = index++; 
+	    Ns[s[*v]]++;
+	    size_t k = out_degree_no_loops(*v,g);
+	    Nk[k]++;
+	    kNs[k][s[*v]]++;
+	    Ks[s[*v]] += k;
+	}
+
+	for (tie(v,v_end) = vertices(g); v != v_end; ++v)
+	{
+	    size_t k = out_degree_no_loops(*v,g);
+	    global_term[k][gamma*Nnnks(k,s[*v],Ns,kNs,Ks)].insert(s[*v]);
+	}
+
+	// define cooling rate so that temperature starts at "infinity" (numeric_limits::max()) at temp_count == 0
+	// and reaches "zero" (numeric_limits::epsilon()) at temp_count == n_iter - 1
+	double cooling_rate = -(log(numeric_limits<double>::epsilon())-log(numeric_limits<double>::max()))/(n_iter-1);
+
+	for (size_t temp_count = 0; temp_count < n_iter; ++temp_count)
+	{
+	    double T = numeric_limits<double>::max()*exp(-cooling_rate*temp_count);
+	    
+            // sample a new spin for every vertex
+	    for (tie(v,v_end) = vertices(g); v != v_end; ++v)
+	    {
+		double E = -T*log(1-uniform_p(rng)); // sampled target jump energy
+		double local_term = 0;
+		
+		unordered_map<size_t, double> ns; // number of neighbours with spin 's' (weighted)
+		
+		// neighbourhood spins info
+		typename graph_traits<Graph>::out_edge_iterator e,e_end;
+		for (tie(e,e_end) = out_edges(*v,g); e != e_end; ++e)
+		{
+		    typename graph_traits<Graph>::vertex_descriptor t = target(*e,g);
+		    if (t != *v)
+			ns[s[t]] += get(weights, *e);
+		}
+		    
+		size_t k = out_degree_no_loops(*v,g);
+		
+		// local energy term
+		local_term = ns[s[*v]] + gamma*Nnnks(k,s[*v],Ns,kNs,Ks);
+		
+		// update global info with local info
+		for (typeof(ns.begin()) iter = ns.begin(); iter != ns.end(); ++iter)
+		{
+		    double old_E = gamma*Nnnks(k,iter->first,Ns,kNs,Ks);
+		    global_term[k][old_E].erase(iter->first);
+		    if (global_term[k][old_E].empty())
+			global_term[k].erase(old_E);
+		    global_term[k][old_E - ns[iter->first]].insert(iter->first);
+		}
+		    
+		// fetch spins with closest energy
+		E -= local_term;
+		
+		typeof(global_term[k].begin()) lower = global_term[k].lower_bound(E);
+		if (lower == global_term[k].end())
+		    lower = global_term[k].begin();
+		typeof(global_term[k].begin()) upper = global_term[k].upper_bound(E);
+		if (upper == global_term[k].end())
+		    --upper;
+		    
+		typeof(global_term[k].begin()) closest = (abs(E - lower->first) < abs(E-upper->first))?lower:upper;
+		    
+		//new spin (randomly chosen amongst those with equal energy)
+		uniform_int<size_t> sample_spin(0,closest->second.size()-1);
+		typeof(closest->second.begin()) iter = closest->second.begin();
+		advance(iter, sample_spin(rng));
+		size_t a = *iter;
+		
+		//cleanup global info
+		for (typeof(ns.begin()) iter = ns.begin(); iter != ns.end(); ++iter)
+		{
+		    double old_E = gamma*Nnnks(k,iter->first,Ns,kNs,Ks) - ns[iter->first];		    
+		    global_term[k][old_E].erase(iter->first);
+		    if (global_term[k][old_E].empty())
+			global_term[k].erase(old_E);
+		    global_term[k][old_E + ns[iter->first]].insert(iter->first);
+		}
+
+		//update global info
+		double old_E = gamma*Nnnks(k,s[*v],Ns,kNs,Ks);
+		global_term[k][old_E].erase(s[*v]);
+		if (global_term[k][old_E].empty())
+		    global_term[k].erase(old_E);
+		Ns[s[*v]]--;
+		kNs[k][s[*v]]--;
+		Ks[s[*v]] -= k;
+		Ns[a]++;
+		kNs[k][a]++;
+		Ks[a] += k;
+		global_term[k][gamma*Nnnks(k,s[*v],Ns,kNs,Ks)].insert(s[*v]);
+		global_term[k][gamma*Nnnks(k,a,Ns,kNs,Ks)].insert(a);
+		
+		// update spin
+		s[*v] = a;
+	    }		
+	}	
+    }
+};
+
+template <class Graph>
+class NNKSErdosReyni
+{
+public:
+    NNKSErdosReyni(Graph &g)
+    {
+	size_t N = 0;
+	double _avg_k = 0.0;
+	typename graph_traits<Graph>::vertex_iterator v,v_end;
+	for (tie(v,v_end) = vertices(g); v != v_end; ++v)
+	{
+	    size_t k = out_degree_no_loops(*v,g); 
+	    _avg_k += k;
+	    N++;
+	}
+	_p = _avg_k/(N*N);
+    }
+
+    double operator()(size_t k, size_t s, unordered_map<size_t,size_t>& Ns, unordered_map<size_t,unordered_map<size_t,size_t> >& kNs,  
+		      unordered_map<size_t, size_t>& Ks) const 
+    {
+	return _p*Ns[s];
+    }
+
+private:
+    double _p;
+};
+
+template <class Graph>
+class NNKSUncorr
+{
+public:
+    NNKSUncorr(Graph &g)
+	:_K(0)
+    {
+	typename graph_traits<Graph>::vertex_iterator v,v_end;
+	for (tie(v,v_end) = vertices(g); v != v_end; ++v)
+	    _K += out_degree_no_loops(*v,g); 
+    }
+
+    double operator()(size_t k, size_t s, unordered_map<size_t,size_t>& Ns, unordered_map<size_t,unordered_map<size_t,size_t> >& kNs,  
+		      unordered_map<size_t, size_t>& Ks) const 
+    {
+	return k*Ks[s]/double(_K);
+    }
+private:
+    size_t _K;
+};
+
+template <class Graph>
+class NNKSCorr
+{
+public:
+    NNKSCorr(Graph &g)
+    {
+	typename graph_traits<Graph>::vertex_iterator v,v_end;
+	for (tie(v,v_end) = vertices(g); v != v_end; ++v)
+	{
+	    size_t k = out_degree_no_loops(*v,g);
+	    _Nk[k]++;
+	    typename graph_traits<Graph>::adjacency_iterator a,a_end;
+	    for (tie(a,a_end) = adjacent_vertices(*v,g); a != a_end; ++a)
+		if (*a != *v)
+		    _Nkk[k][out_degree_no_loops(*a,g)]++;
+	}
+    }
+
+    double operator()(size_t k, size_t s, unordered_map<size_t,size_t>& Ns, unordered_map<size_t,unordered_map<size_t,size_t> >& kNs,  
+		      unordered_map<size_t, size_t>& Ks)
+    {
+	double N = 0;
+	unordered_map<size_t, size_t>& nkk = _Nkk[k];
+	for (typeof(nkk.begin()) iter = nkk.begin(); iter != nkk.end(); ++iter)
+	    N += iter->second * kNs[iter->first][s]/double(_Nk[iter->first]);
+	return N;
+    }
+private:
+    unordered_map<size_t, size_t> _Nk;
+    unordered_map<size_t, unordered_map<size_t, size_t> > _Nkk;
+};
+
+
+struct get_communities_selector
+{
+    get_communities_selector(GraphInterface::comm_corr_t corr):_corr(corr) {}
+    GraphInterface::comm_corr_t _corr;
+
+    template <class Graph, class WeightMap, class CommunityMap>
+    void operator()(Graph& g, WeightMap weights, CommunityMap s, double gamma, size_t n_iter, size_t seed) const
+    {
+	switch (_corr)
+	{
+	case GraphInterface::ERDOS_REYNI:
+	    get_communities<NNKSErdosReyni>()(g, weights, s, gamma, n_iter, seed);
+	    break;
+	case GraphInterface::UNCORRELATED:
+	    get_communities<NNKSUncorr>()(g, weights, s, gamma, n_iter, seed);
+	    break;
+	case GraphInterface::CORRELATED:
+	    get_communities<NNKSCorr>()(g, weights, s, gamma, n_iter, seed);
+	    break;
+	}
+    }
+};
+
+void GraphInterface::GetCommunityStructure(double gamma, comm_corr_t corr, size_t n_iter, size_t seed, string weight, string property)
+{
+    typedef HashedDescriptorMap<vertex_index_map_t,size_t> comm_map_t;
+    comm_map_t comm_map(_vertex_index);
+
+    bool directed = _directed;
+    _directed = false;
+
+    if(weight != "")
+    {
+	try 
+	{
+	    dynamic_property_map& weight_prop = find_property_map(_properties, weight, typeid(graph_traits<multigraph_t>::edge_descriptor));
+	    if (get_static_property_map<vector_property_map<double,edge_index_map_t> >(&weight_prop))
+	    {
+		vector_property_map<double,edge_index_map_t> weight_map = 
+		    get_static_property_map<vector_property_map<double,edge_index_map_t> >(weight_prop);
+		check_filter(*this, bind<void>(get_communities_selector(corr), _1, var(weight_map), var(comm_map), gamma, n_iter, seed),
+			     reverse_check(), always_undirected());
+	    }
+	    else
+	    {
+		DynamicPropertyMapWrap<double,graph_traits<multigraph_t>::edge_descriptor> weight_map(weight_prop);
+		check_filter(*this, bind<void>(get_communities_selector(corr), _1, var(weight_map), var(comm_map), gamma, n_iter, seed),
+			     reverse_check(), always_undirected());
+	    }
+	}
+	catch (property_not_found& e)
+	{
+	    throw GraphException("error getting scalar property: " + string(e.what()));
+	}
+    }
+    else
+    {
+	ConstantPropertyMap<double,graph_traits<multigraph_t>::edge_descriptor> weight_map(1.0);
+	check_filter(*this, bind<void>(get_communities_selector(corr), _1, var(weight_map), var(comm_map), gamma, n_iter, seed),
+		     reverse_check(), always_undirected());
+    }
+    _directed = directed;
+
+    try
+    {
+	find_property_map(_properties, property, typeid(graph_traits<multigraph_t>::edge_descriptor));
+	RemoveVertexProperty(property);
+    }
+    catch (property_not_found) {}
+
+    _properties.property(property, comm_map);
+}
+
+//==============================================================================
+// GetModularity()
+// get Newman's modularity of a given community partition
+//==============================================================================
+
+struct get_modularity
+{
+    template <class Graph, class WeightMap, class CommunityMap>
+    void operator()(Graph& g, WeightMap weights, CommunityMap s, double& modularity) const
+    {
+	modularity = 0.0;
+	
+	size_t E = 0;
+	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))
+	    {
+		W += get(weights, *e);
+		E++;
+	    }
+
+	typename graph_traits<Graph>::vertex_iterator v,v_end;
+	for (tie(v,v_end) = vertices(g); v != v_end; ++v)
+	{
+	    size_t k = out_degree_no_loops(*v,g);
+	    typename graph_traits<Graph>::out_edge_iterator e,e_end;
+	    for(tie(e, e_end) = out_edges(*v,g); e != e_end; ++e)
+	    {
+		typename graph_traits<Graph>::vertex_descriptor t = target(*e,g);
+		if (t != *v)
+		{
+		    if (get(s, t) == get(s, *v))
+			modularity += get(weights, *e) - k*out_degree_no_loops(t,g)/double(2*E);
+		}
+	    }
+	}
+	modularity /= 2*W;
+    }
+};
+
+double GraphInterface::GetModularity(string weight, string property)
+{
+    double modularity = 0;
+
+    bool directed = _directed;
+    _directed = false;
+    try
+    {
+	dynamic_property_map& comm_prop = find_property_map(_properties, property, typeid(graph_traits<multigraph_t>::vertex_descriptor));