Commit 4974404c authored by Tiago Peixoto's avatar Tiago Peixoto
Browse files

Slightly refactor random graph generation and add generate_sbm()

parent dd1e9dfa
......@@ -28,6 +28,7 @@ libgraph_tool_generation_la_SOURCES = \
graph_predecessor.cc \
graph_price.cc \
graph_rewiring.cc \
graph_sbm.cc \
graph_triangulation.cc \
graph_union.cc \
graph_union_eprop.cc \
......@@ -44,6 +45,7 @@ libgraph_tool_generation_la_include_HEADERS = \
graph_predecessor.hh \
graph_price.hh \
graph_rewiring.hh \
graph_sbm.hh \
graph_triangulation.hh \
graph_union.hh \
sampler.hh
......@@ -68,6 +68,10 @@ void generate_graph(GraphInterface& gi, size_t N,
std::ref(rng), verbose, verify))();
}
void generate_sbm(GraphInterface& gi, boost::any ab, boost::python::object ors,
boost::python::object oss, boost::python::object oprobs,
boost::any ain_deg, boost::any aout_deg, rng_t& rng);
size_t random_rewire(GraphInterface& gi, string strat, size_t niter,
bool no_sweep, bool self_loops, bool parallel_edges,
bool alias, bool traditional, bool persist,
......@@ -119,6 +123,7 @@ using namespace boost::python;
BOOST_PYTHON_MODULE(libgraph_tool_generation)
{
def("gen_graph", &generate_graph);
def("gen_sbm", &generate_sbm);
def("random_rewire", &random_rewire);
def("predecessor_graph", &predecessor_graph);
def("line_graph", &line_graph);
......
......@@ -139,7 +139,6 @@ size_t random_rewire(GraphInterface& gi, string strat, size_t niter,
emap_t::unchecked_t pin =
any_cast<emap_t>(apin).get_unchecked(gi.get_edge_index_range());
if (strat == "erdos")
{
run_action<graph_tool::detail::never_reversed>()
......@@ -151,7 +150,7 @@ size_t random_rewire(GraphInterface& gi, string strat, size_t niter,
std::make_tuple(persist, cache, verbose),
std::ref(pcount), std::ref(rng)))();
}
else if (strat == "uncorrelated")
else if (strat == "configuration")
{
run_action<graph_tool::detail::never_reversed>()
(gi, std::bind(graph_rewire<RandomRewireStrategy>(),
......@@ -161,7 +160,7 @@ size_t random_rewire(GraphInterface& gi, string strat, size_t niter,
std::make_tuple(persist, cache, verbose),
std::ref(pcount), std::ref(rng)))();
}
else if (strat == "correlated")
else if (strat == "constrained-configuration")
{
if (block.empty())
{
......@@ -186,7 +185,7 @@ size_t random_rewire(GraphInterface& gi, string strat, size_t niter,
vertex_properties())(block);
}
}
else if (strat == "probabilistic")
else if (strat == "probabilistic-configuration")
{
run_action<>()
(gi, std::bind(graph_rewire<ProbabilisticRewireStrategy>(),
......
......@@ -710,6 +710,9 @@ public:
for (auto iter = _probs.begin(); iter != _probs.end(); ++iter)
{
double& p = iter->second;
// avoid zero probability to not get stuck in rejection step
if (std::isnan(p) || std::isinf(p) || p <= 0)
p = numeric_limits<double>::min();
p = log(p);
}
}
......@@ -789,7 +792,6 @@ private:
CorrProb _corr_prob;
BlockDeg _blockdeg;
typedef std::unordered_map<pair<deg_t, deg_t>, double> prob_map_t;
prob_map_t _probs;
};
......@@ -1020,7 +1022,7 @@ public:
double pi = get_prob(s_deg, t_deg) + get_prob(ep_s_deg, ep_t_deg);
double pf = get_prob(s_deg, ep_t_deg) + get_prob(ep_s_deg, t_deg);
double a = exp(pf - pi);
double a = pf - pi;
if (is_directed::apply<Graph>::type::value)
{
......@@ -1032,10 +1034,10 @@ public:
e_t -= in_degreeS()(ep_t, _g);
}
a /= (get_sprob(s_deg, ep_t_deg) / e_ep_t +
get_sprob(ep_s_deg, t_deg) / e_t); // forwards
a *= (get_sprob(s_deg, t_deg) / e_t +
get_sprob(ep_s_deg, ep_t_deg) / e_ep_t); // backwards
a -= log(get_sprob(s_deg, ep_t_deg)) - log(e_ep_t) +
log(get_sprob(ep_s_deg, t_deg)) - log(e_t); // forwards
a += log(get_sprob(s_deg, t_deg)) - log(e_t) +
log(get_sprob(ep_s_deg, ep_t_deg)) - log(e_ep_t); // backwards
}
else
{
......@@ -1055,22 +1057,22 @@ public:
e_s -= out_degree(ep_s, _g);
}
a /= (get_sprob(s_deg, ep_t_deg) / e_ep_t +
get_sprob(ep_t_deg, s_deg) / e_s +
get_sprob(ep_s_deg, t_deg) / e_t +
get_sprob(t_deg, ep_s_deg) / e_ep_s); // forwards
a *= (get_sprob(s_deg, t_deg) / e_t +
get_sprob(t_deg, s_deg) / e_s +
get_sprob(ep_s_deg, ep_t_deg) / e_ep_t +
get_sprob(ep_t_deg, ep_s_deg) / e_ep_s); // backwards
a -= log(get_sprob(s_deg, ep_t_deg)) - log(e_ep_t) +
log(get_sprob(ep_t_deg, s_deg)) - log(e_s) +
log(get_sprob(ep_s_deg, t_deg)) - log(e_t) +
log(get_sprob(t_deg, ep_s_deg)) - log(e_ep_s); // forwards
a += log(get_sprob(s_deg, t_deg)) - log(e_t) +
log(get_sprob(t_deg, s_deg)) - log(e_s) +
log(get_sprob(ep_s_deg, ep_t_deg)) - log(e_ep_t) +
log(get_sprob(ep_t_deg, ep_s_deg)) - log(e_ep_s); // backwards
}
if (a > 1)
if (a > 0)
return ep;
std::uniform_real_distribution<> rsample(0.0, 1.0);
double r = rsample(base_t::_rng);
if (r > a)
if (r > exp(a))
return e; // reject
return ep;
}
......@@ -1163,11 +1165,10 @@ public:
: _g(g), _edge_index(edge_index), _edges(edges), _corr_prob(corr_prob),
_blockdeg(blockdeg), _rng(rng), _sampler(nullptr)
{
typename graph_traits<Graph>::vertex_iterator v, v_end;
for (tie(v, v_end) = vertices(_g); v != v_end; ++v)
for (auto v : vertices_range(_g))
{
deg_t d = _blockdeg.get_block(*v, g);
_vertices[d].push_back(*v);
deg_t d = _blockdeg.get_block(v, g);
_vertices[d].push_back(v);
}
std::unordered_map<pair<deg_t, deg_t>, double> probs;
......@@ -1176,31 +1177,31 @@ public:
vector<double> dprobs;
if (probs.empty())
{
for (auto s_iter = _vertices.begin(); s_iter != _vertices.end(); ++s_iter)
for (auto& s : _vertices)
{
for (auto t_iter = _vertices.begin(); t_iter != _vertices.end(); ++t_iter)
for (auto& t : _vertices)
{
double p = _corr_prob(s_iter->first, t_iter->first);
double p = _corr_prob(s.first, t.first);
if (std::isnan(p) || std::isinf(p) || p <= 0)
continue;
_items.push_back(make_pair(s_iter->first, t_iter->first));
dprobs.push_back(p);
_items.push_back(make_pair(s.first, t.first));
dprobs.push_back(p * s.second.size() * t.second.size());
}
}
}
else
{
for (auto iter = probs.begin(); iter != probs.end(); ++iter)
for (auto& stp : probs)
{
deg_t s = iter->first.first;
deg_t t = iter->first.second;
double p = iter->second;
deg_t s = stp.first.first;
deg_t t = stp.first.second;
double p = stp.second;
// avoid zero probability to not get stuck in rejection step
if (std::isnan(p) || std::isinf(p) || p <= 0)
continue;
_items.push_back(make_pair(s, t));
dprobs.push_back(p);
dprobs.push_back(p * _vertices[s].size() * _vertices[t].size());
}
}
......@@ -1230,11 +1231,8 @@ public:
if (svs.empty() || tvs.empty())
continue;
std::uniform_int_distribution<size_t> s_sample(0, svs.size() - 1);
std::uniform_int_distribution<size_t> t_sample(0, tvs.size() - 1);
s = svs[s_sample(_rng)];
t = tvs[t_sample(_rng)];
s = uniform_sample(svs, _rng);
t = uniform_sample(tvs, _rng);
break;
}
......
// graph-tool -- a general graph modification and manipulation thingy
//
// Copyright (C) 2006-2016 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/>.
#include "graph_sbm.hh"
#include "numpy_bind.hh"
using namespace std;
using namespace boost;
using namespace graph_tool;
void generate_sbm(GraphInterface& gi, boost::any ab, boost::python::object ors,
boost::python::object oss, boost::python::object oprobs,
boost::any ain_deg, boost::any aout_deg, rng_t& rng)
{
auto rs = get_array<int64_t, 1>(ors);
auto ss = get_array<int64_t, 1>(oss);
auto probs = get_array<double, 2>(oprobs);
typedef vprop_map_t<int32_t>::type bmap_t;
auto b = any_cast<bmap_t>(ab).get_unchecked();
typedef vprop_map_t<double>::type dmap_t;
auto in_deg = any_cast<dmap_t>(ain_deg).get_unchecked();
auto out_deg = any_cast<dmap_t>(aout_deg).get_unchecked();
run_action<graph_tool::detail::always_directed_never_reversed>()
(gi, [&](auto& g) { gen_sbm(g, b, rs, ss, probs, in_deg, out_deg, rng); })();
}
// graph-tool -- a general graph modification and manipulation thingy
//
// Copyright (C) 2006-2016 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/>.
#ifndef GRAPH_SBM_HH
#define GRAPH_SBM_HH
#include <tuple>
#include <iostream>
#include <boost/functional/hash.hpp>
#include "graph.hh"
#include "graph_filtering.hh"
#include "graph_util.hh"
#include "sampler.hh"
#include "random.hh"
#include "hash_map_wrap.hh"
namespace graph_tool
{
using namespace std;
using namespace boost;
template <class Graph, class VProp, class IVec, class FVec, class VDProp,
class RNG>
void gen_sbm(Graph& g, VProp b, IVec& rs, IVec& ss, FVec probs, VDProp in_deg,
VDProp out_deg, RNG& rng)
{
vector<vector<size_t>> rvs;
vector<vector<double>> v_in_probs, v_out_probs;
for (auto v : vertices_range(g))
{
size_t r = b[v];
if (r >= v_in_probs.size())
{
v_in_probs.resize(r+1);
v_out_probs.resize(r+1);
rvs.resize(r+1);
}
rvs[r].push_back(v);
v_in_probs[r].push_back(in_deg[v]);
v_out_probs[r].push_back(out_deg[v]);
}
vector<Sampler<size_t>> v_in_sampler, v_out_sampler;
for (size_t r = 0; r < rvs.size(); ++r)
{
v_in_sampler.emplace_back(rvs[r], v_in_probs[r]);
v_out_sampler.emplace_back(rvs[r], v_out_probs[r]);
}
for (size_t i = 0; i < rs.shape()[0]; ++i)
{
size_t r = rs[i];
size_t s = ss[i];
double p = probs[0][i];
if (r >= v_out_sampler.size() || v_out_sampler[r].prob_sum() == 0 ||
s >= v_in_sampler.size() || v_in_sampler[s].prob_sum() == 0)
throw GraphException("Inconsistent SBM parameters: edge probabilities given for empty groups");
std::poisson_distribution<> poi(p);
size_t ers = poi(rng);
for (size_t j = 0; j < ers; ++j)
{
size_t u = v_out_sampler[r].sample(rng);
size_t v = v_in_sampler[s].sample(rng);
add_edge(u, v, g);
}
}
}
} // graph_tool namespace
#endif // GRAPH_SBM_HH
......@@ -38,18 +38,18 @@ class Sampler
public:
Sampler(const vector<Value>& items,
const vector<double>& probs)
: _items(items), _probs(probs), _alias(items.size())
: _items(items), _probs(probs), _alias(items.size()),
_S(0)
{
double S = 0;
for (size_t i = 0; i < _probs.size(); ++i)
S += _probs[i];
_S += _probs[i];
vector<size_t> small;
vector<size_t> large;
for (size_t i = 0; i < _probs.size(); ++i)
{
_probs[i] *= _probs.size() / S;
_probs[i] *= _probs.size() / _S;
if (_probs[i] < 1)
small.push_back(i);
else
......@@ -95,6 +95,7 @@ public:
size_t size() const { return _items.size(); }
bool empty() const { return _items.empty(); }
double prob_sum() const { return _S; }
const Value& operator[](size_t i) const
{
......@@ -125,7 +126,7 @@ private:
vector<double> _probs;
vector<size_t> _alias;
uniform_int_distribution<size_t> _sample;
double _S;
};
// uniform sampling from containers
......
......@@ -252,8 +252,9 @@ def corr_hist(g, deg_source, deg_target, bins=[[0, 1], [0, 1]], weight=None,
... accept = np.random.random() < 1.0/k
... return k
...
>>> g = gt.random_graph(10000, lambda: sample_k(40), model="probabilistic",
... vertex_corr=lambda i, j: (sin(i / pi) * sin(j / pi) + 1) / 2,
>>> g = gt.random_graph(10000, lambda: sample_k(40),
... model="probabilistic-configuration",
... edge_probs=lambda i, j: (sin(i / pi) * sin(j / pi) + 1) / 2,
... directed=False, n_iter=100)
>>> h = gt.corr_hist(g, "out", "out")
>>> clf()
......@@ -453,8 +454,9 @@ def avg_neighbour_corr(g, deg_source, deg_target, bins=[0, 1], weight=None):
... accept = np.random.random() < 1.0 / k
... return k
...
>>> g = gt.random_graph(10000, lambda: sample_k(40), model="probabilistic",
... vertex_corr=lambda i, j: (sin(i / pi) * sin(j / pi) + 1) / 2,
>>> g = gt.random_graph(10000, lambda: sample_k(40),
... model="probabilistic-configuration",
... edge_probs=lambda i, j: (sin(i / pi) * sin(j / pi) + 1) / 2,
... directed=False, n_iter=100)
>>> h = gt.avg_neighbour_corr(g, "out", "out")
>>> clf()
......
This diff is collapsed.
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