Commit ab6bef60 authored by Tiago Peixoto's avatar Tiago Peixoto
Browse files

Refactor SBM code and implement overlapping model

parent 2905fe8c
......@@ -4,6 +4,7 @@
.. autofunction:: minimize_blockmodel_dl
.. autoclass:: BlockState
.. autoclass:: OverlapBlockState
.. autofunction:: mcmc_sweep
.. autoclass:: MinimizeState
.. autofunction:: multilevel_minimize
......@@ -22,5 +23,6 @@
.. autofunction:: nested_mcmc_sweep
.. autofunction:: nested_tree_sweep
.. autofunction:: get_hierarchy_tree
.. autofunction:: get_block_edge_gradient
.. autofunction:: community_structure
.. autofunction:: modularity
......@@ -16,14 +16,17 @@ libgraph_tool_community_la_LDFLAGS = $(MOD_LDFLAGS)
libgraph_tool_community_la_SOURCES = \
graph_blockmodel.cc \
graph_blockmodel_overlap.cc \
graph_community.cc \
graph_community_network.cc \
graph_community_network_edges.cc \
graph_community_network_vavg.cc \
graph_community_network_eavg.cc \
graph_community_network_eavg_imp1.cc
graph_community_network_eavg_imp1.cc \
spence.cc
libgraph_tool_community_la_include_HEADERS = \
graph_blockmodel.hh \
graph_blockmodel_overlap.hh \
graph_community.hh \
graph_community_network.hh
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
......@@ -121,6 +121,7 @@ void community_network_eavg(GraphInterface& gi, GraphInterface& cgi,
extern void export_blockmodel();
extern void export_blockmodel_overlap();
BOOST_PYTHON_MODULE(libgraph_tool_community)
{
......@@ -131,4 +132,5 @@ BOOST_PYTHON_MODULE(libgraph_tool_community)
def("community_network_eavg", &community_network_eavg);
export_blockmodel();
export_blockmodel_overlap();
}
......@@ -133,10 +133,9 @@ struct get_community_network_edges
typedef std::unordered_map<cvertex_t, cedge_t> ecomms_t;
#endif
unchecked_vector_property_map<ecomms_t,
typename property_map<CommunityGraph, vertex_index_t>::type>
comm_edges(get(vertex_index_t(), cg),
num_vertices(cg));
auto index_map = get(vertex_index_t(), cg);
unchecked_vector_property_map<ecomms_t, decltype(index_map)>
comm_edges(index_map, num_vertices(cg));
typename graph_traits<CommunityGraph>::vertex_iterator v, v_end;
for (tie(v, v_end) = vertices(cg); v != v_end; ++v)
......
// graph-tool -- a general graph modification and manipulation thingy
//
// Copyright (C) 2006-2014 Tiago de Paula Peixoto <tiago@skewed.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/>.
/*
* Cephes Math Library Release 2.1: January, 1989
* Copyright 1985, 1987, 1989 by Stephen L. Moshier
* Direct inquiries to 30 Frost Street, Cambridge, MA 02140
*/
#include <cmath>
#include <limits>
double MACHEP = 1.11022302462515654042E-16; /* 2**-53 */
static double A[8] = {
4.65128586073990045278E-5,
7.31589045238094711071E-3,
1.33847639578309018650E-1,
8.79691311754530315341E-1,
2.71149851196553469920E0,
4.25697156008121755724E0,
3.29771340985225106936E0,
1.00000000000000000126E0,
};
static double B[8] = {
6.90990488912553276999E-4,
2.54043763932544379113E-2,
2.82974860602568089943E-1,
1.41172597751831069617E0,
3.63800533345137075418E0,
5.03278880143316990390E0,
3.54771340985225096217E0,
9.99999999999999998740E-1,
};
double polevl(double x, double coef[], int N)
{
double ans;
int i;
double *p;
p = coef;
ans = *p++;
i = N;
do
ans = ans * x + *p++;
while (--i);
return (ans);
}
double spence(double x)
{
double w, y, z;
int flag;
if (x < 0.0)
return std::numeric_limits<double>::quiet_NaN();
if (x == 1.0)
return (0.0);
if (x == 0.0)
return (M_PI * M_PI / 6.0);
flag = 0;
if (x > 2.0) {
x = 1.0 / x;
flag |= 2;
}
if (x > 1.5) {
w = (1.0 / x) - 1.0;
flag |= 2;
}
else if (x < 0.5) {
w = -x;
flag |= 1;
}
else
w = x - 1.0;
y = -w * polevl(w, A, 7) / polevl(w, B, 7);
if (flag & 1)
y = (M_PI * M_PI) / 6.0 - std::log(x) * std::log(1.0 - x) - y;
if (flag & 2) {
z = std::log(x);
y = -0.5 * z * z - y;
}
return (y);
}
......@@ -92,6 +92,8 @@ public:
return _items[_alias[i]];
}
size_t size() const { return _items.size(); }
private:
typedef typename mpl::if_<KeepReference,
......@@ -104,6 +106,14 @@ private:
vector<size_t> _large;
};
// uniform sampling from containers
template <class Container, class RNG>
typename Container::value_type& uniform_sample(Container& v, RNG& rng)
{
std::uniform_int_distribution<size_t> i_rand(0, v.size() - 1);
return v[i_rand(rng)];
}
} // namespace graph_tool
......
......@@ -64,7 +64,8 @@ graph_tool_centralitydir = $(MOD_DIR)/centrality
graph_tool_community_PYTHON = \
community/__init__.py \
community/blockmodel.py \
community/nested_blockmodel.py
community/nested_blockmodel.py \
community/overlap_blockmodel.py
graph_tool_communitydir = $(MOD_DIR)/community
graph_tool_draw_PYTHON = \
......
......@@ -704,6 +704,8 @@ class PropertyMap(object):
def __getstate__(self):
g = self.get_graph()
if g is None:
raise ValueError("cannot pickle orphaned property map")
value_type = self.value_type()
key_type = self.key_type()
if not self.is_writable():
......
......@@ -40,6 +40,7 @@ Summary
minimize_blockmodel_dl
BlockState
OverlapBlockState
mcmc_sweep
MinimizeState
multilevel_minimize
......@@ -51,6 +52,7 @@ Summary
get_max_B
get_akc
condensation_graph
get_block_edge_gradient
Hierarchical models
===================
......@@ -112,6 +114,7 @@ __all__ = ["minimize_blockmodel_dl",
"get_max_B",
"get_akc",
"condensation_graph",
"OverlapBlockState",
"minimize_nested_blockmodel_dl",
"NestedBlockState",
"NestedMinimizeState",
......@@ -119,6 +122,7 @@ __all__ = ["minimize_blockmodel_dl",
"nested_mcmc_sweep",
"nested_tree_sweep",
"get_hierarchy_tree",
"get_block_edge_gradient",
"community_structure",
"modularity"]
......@@ -126,14 +130,30 @@ from . blockmodel import minimize_blockmodel_dl, BlockState, mcmc_sweep, \
multilevel_minimize, model_entropy, get_max_B, get_akc, condensation_graph, \
collect_edge_marginals, collect_vertex_marginals, bethe_entropy, mf_entropy, MinimizeState
from . overlap_blockmodel import OverlapBlockState, get_block_edge_gradient
from . nested_blockmodel import NestedBlockState, NestedMinimizeState, init_nested_state, \
nested_mcmc_sweep, nested_tree_sweep, minimize_nested_blockmodel_dl, get_hierarchy_tree
def community_structure(g, n_iter, n_spins, gamma=1.0, corr="erdos",
spins=None, weight=None, t_range=(100.0, 0.01),
verbose=False, history_file=None):
r"""
Obtain the community structure for the given graph, using a Potts model approach.
r"""Obtain the community structure for the given graph, using a Potts model
approach, which is a generalization of modularity maximization).
.. warning::
**The use of this function is discouraged.** Although community detection
based on modularity maximization is very common, it is very
problematic. It will find high-scoring partitions where there is none
[guimera-modularity-2004]_, and at the same will not find actual
structure in large graphs [fortunato-resolution-2007]_. Furthermore, in
many empirical networks, the partitions found in this way are largely
meaningless [good-performance-2010]_.
Instead, use the methods based on statistical inference
(i.e. :func:`~graph_tool.community.minimize_blockmodel_dl` and
:func:`~graph_tool.community.minimize_nested_blockmodel_dl`).
Parameters
----------
......@@ -226,8 +246,7 @@ def community_structure(g, n_iter, n_spins, gamma=1.0, corr="erdos",
>>> seed(42)
>>> g = gt.load_graph("community.xml")
>>> pos = g.vertex_properties["pos"]
>>> spins = gt.community_structure(g, 10000, 20, t_range=(5, 0.1),
... history_file="community-history1")
>>> spins = gt.community_structure(g, 10000, 20, t_range=(5, 0.1))
>>> gt.graph_draw(g, pos=pos, vertex_fill_color=spins, output_size=(420, 420), output="comm1.pdf")
<...>
......@@ -236,8 +255,7 @@ def community_structure(g, n_iter, n_spins, gamma=1.0, corr="erdos",
gt.graph_draw(g, pos=pos, vertex_fill_color=spins, output_size=(420, 420), output="comm1.png")
>>> spins = gt.community_structure(g, 10000, 40, t_range=(5, 0.1),
... gamma=2.5, history_file="community-history2")
>>> spins = gt.community_structure(g, 10000, 40, t_range=(5, 0.1), gamma=2.5)
>>> gt.graph_draw(g, pos=pos, vertex_fill_color=spins, output_size=(420, 420), output="comm2.pdf")
<...>
......@@ -246,51 +264,31 @@ def community_structure(g, n_iter, n_spins, gamma=1.0, corr="erdos",
gt.graph_draw(g, pos=pos, vertex_fill_color=spins, output_size=(420, 420), output="comm2.png")
>>> figure(figsize=(6, 4))
<...>
>>> xlabel("iterations")
<...>
>>> ylabel("number of communities")
<...>
>>> a = loadtxt("community-history1").transpose()
>>> plot(a[0], a[2])
[...]
>>> savefig("comm1-hist.pdf")
.. testcode::
:hide:
savefig("comm1-hist.png")
>>> clf()
>>> xlabel("iterations")
<...>
>>> ylabel("number of communities")
<...>
>>> a = loadtxt("community-history2").transpose()
>>> plot(a[0], a[2])
[...]
>>> savefig("comm2-hist.pdf")
.. figure:: comm1.*
:align: center
.. testcode::
:hide:
The community structure with :math:`\gamma=1`.
savefig("comm2-hist.png")
.. figure:: comm2.*
:align: center
The community structure with :math:`\gamma=1`:
.. image:: comm1.*
.. image:: comm1-hist.*
The community structure with :math:`\gamma=2.5`:
.. image:: comm2.*
.. image:: comm2-hist.*
The community structure with :math:`\gamma=2.5`.
References
----------
.. [guimera-modularity-2004] Roger Guimerà, Marta Sales-Pardo, and Luís
A. Nunes Amaral, "Modularity from fluctuations in random graphs and
complex networks", Phys. Rev. E 70, 025101(R) (2004),
:doi:`10.1103/PhysRevE.70.025101`
.. [fortunato-resolution-2007] Santo Fortunato and Marc Barthélemy,
"Resolution limit in community detection", Proc. Natl. Acad. Sci. USA
104(1): 36–41 (2007), :doi:`10.1073/pnas.0605965104`
.. [good-performance-2010] Benjamin H. Good, Yves-Alexandre de Montjoye, and
Aaron Clauset, "Performance of modularity maximization in practical
contexts", Phys. Rev. E 81, 046106 (2010),
:doi:`10.1103/PhysRevE.81.046106`
.. [reichard-statistical-2006] Joerg Reichardt and Stefan Bornholdt,
"Statistical Mechanics of Community Detection", Phys. Rev. E 74
016110 (2006), :doi:`10.1103/PhysRevE.74.016110`, :arxiv:`cond-mat/0603718`
......
This diff is collapsed.
This diff is collapsed.
......@@ -1156,7 +1156,8 @@ def get_hierarchy_control_points(g, t, tpos, beta=0.8):
>>> g = gt.collection.data["netscience"]
>>> g = gt.GraphView(g, vfilt=gt.label_largest_component(g))
>>> g.purge_vertices()
>>> bstack, mdl = gt.minimize_nested_blockmodel_dl(g, deg_corr=True)
>>> state = gt.minimize_nested_blockmodel_dl(g, deg_corr=True)
>>> bstack = state.get_bstack()
>>> t = gt.get_hierarchy_tree(bstack)[0]
>>> tpos = pos = gt.radial_tree_layout(t, t.vertex(t.num_vertices() - 1), weighted=True)
>>> cts = gt.get_hierarchy_control_points(g, t, tpos)
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment