...
 
Commits (12)
......@@ -104,100 +104,103 @@ for i in range(3 * 8):
g.add_edge(s, t)
for nested in [False, True]:
for directed in [False, True]:
g.set_directed(directed)
hists = {}
for directed in [False, True]:
g.set_directed(directed)
for nested in [False, True]:
if nested:
state = minimize_nested_blockmodel_dl(g, deg_corr=True)
state = state.copy(bs=[g.get_vertices()] + [zeros(1)] * 6, sampling=True)
else:
state = minimize_blockmodel_dl(g, deg_corr=True)
state = state.copy(B=g.num_vertices())
cs = list(reversed([numpy.inf, 1, 0.1, 0.01, "gibbs", -numpy.inf, -1]))
hists = {}
for i, c in enumerate(cs):
if c != "gibbs":
mcmc_args=dict(beta=.3, c=abs(c), niter=40)
else:
mcmc_args=dict(beta=.3, niter=40)
cs = list(reversed([numpy.inf, 1, 0.1, 0.01, "gibbs", -numpy.inf, -1]))
inits = [1, "N"]
for init in inits:
if nested:
get_B = lambda s: [sl.get_nonempty_B() for sl in s.levels]
bs = [g.get_vertices()] if init == "N" else [zeros(1)]
state = state.copy(bs=bs + [zeros(1)] * 6, sampling=True)
else:
get_B = lambda s: [s.get_nonempty_B()]
if i == 0:
mcmc_equilibrate(state,
mcmc_args=mcmc_args,
gibbs=c=="gibbs",
multiflip = c != "gibbs" and c < 0,
nbreaks=5,
wait=1000,
callback=get_B,
verbose=(1, "c = %s (t) " % str(c)) if verbose else False)
hists[c] = mcmc_equilibrate(state,
mcmc_args=mcmc_args,
gibbs=c=="gibbs",
multiflip = c != "gibbs" and c < 0,
wait=6000,
nbreaks=10,
callback=get_B,
verbose=(1, "c = %s " % str(c)) if verbose else False,
history=True)
state = state.copy(B=g.num_vertices() if init == "N" else 1)
for i, c in enumerate(cs):
if c != "gibbs":
mcmc_args=dict(beta=.3, c=abs(c), niter=40)
else:
mcmc_args=dict(beta=.3, niter=40)
if nested:
get_B = lambda s: [sl.get_nonempty_B() for sl in s.levels]
else:
get_B = lambda s: [s.get_nonempty_B()]
hists[(c, init)] = mcmc_equilibrate(state,
mcmc_args=mcmc_args,
gibbs=c=="gibbs",
multiflip = c != "gibbs" and c < 0,
wait=6000,
nbreaks=6,
callback=get_B,
verbose=(1, "c = %s " % str(c)) if verbose else False,
history=True)
output = open("test_mcmc_nested%s_directed%s-output" % (nested, directed), "w")
for c1 in cs:
for c2 in cs:
try:
if c2 < c1:
continue
except TypeError:
pass
Ss1 = array(list(zip(*hists[c1]))[2])
Ss2 = array(list(zip(*hists[c2]))[2])
# add very small normal noise, to solve discreteness issue
Ss1 += numpy.random.normal(0, 1e-2, len(Ss1))
Ss2 += numpy.random.normal(0, 1e-2, len(Ss2))
D, p = scipy.stats.ks_2samp(Ss1, Ss2)
D_c = 1.63 * sqrt((len(Ss1) + len(Ss2)) / (len(Ss1) * len(Ss2)))
print("nested:", nested, "directed:", directed, "c1:", c1,
"c2:", c2, "D:", D, "D_c:", D_c, "p-value:", p,
file=output)
if p < .001:
print(("Warning, distributions for directed=%s (c1, c2) = " +
"(%s, %s) are not the same, with a p-value: %g (D=%g, D_c=%g)") %
(str(directed), str(c1), str(c2), p, D, D_c),
file=output)
output.flush()
for init1 in inits:
for c2 in cs:
for init2 in inits:
try:
if (c2, init2) < (c1, init1):
continue
except TypeError:
pass
Ss1 = array(list(zip(*hists[(c1, init1)]))[2])
Ss2 = array(list(zip(*hists[(c2, init2)]))[2])
# add very small normal noise, to solve discreteness issue
Ss1 += numpy.random.normal(0, 1e-2, len(Ss1))
Ss2 += numpy.random.normal(0, 1e-2, len(Ss2))
D, p = scipy.stats.ks_2samp(Ss1, Ss2)
D_c = 1.63 * sqrt((len(Ss1) + len(Ss2)) / (len(Ss1) * len(Ss2)))
print("nested:", nested, "directed:", directed, "c1:", c1,
"init1:", init1, "c2:", c2, "init2:", init2, "D:",
D, "D_c:", D_c, "p-value:", p,
file=output)
if p < .001:
print(("Warning, distributions for directed=%s (c1, c2) = " +
"(%s, %s) are not the same, with a p-value: %g (D=%g, D_c=%g)") %
(str(directed), str((c1, init1)),
str((c2, init2)), p, D, D_c),
file=output)
output.flush()
for cum in [True, False]:
clf()
bins = None
for c in cs:
hist = hists[c]
if cum:
h = histogram(list(zip(*hist))[2], 1000000, density=True)
y = numpy.cumsum(h[0])
y /= y[-1]
plot(h[-1][:-1], y, "-", label="c=%s" % str(c))
if c != numpy.inf:
hist = hists[numpy.inf]
h2 = histogram(list(zip(*hist))[2], bins=h[-1], density=True)
y2 = numpy.cumsum(h2[0])
y2 /= y2[-1]
res = abs(y - y2)
i = res.argmax()
axvline(h[-1][i], color="grey")
else:
h = histogram(list(zip(*hist))[2], 40, density=True)
plot(h[-1][:-1], h[0], "s-", label="c=%s" % str(c))
gca().set_yscale("log")
for init in inits:
hist = hists[(c,init)]
if cum:
h = histogram(list(zip(*hist))[2], 1000000,
density=True)
y = numpy.cumsum(h[0])
y /= y[-1]
plot(h[-1][:-1], y, "-", label="c=%s, init=%s" % (str(c), init))
if c != numpy.inf:
hist = hists[(numpy.inf, 1)]
h2 = histogram(list(zip(*hist))[2], bins=h[-1],
density=True)
y2 = numpy.cumsum(h2[0])
y2 /= y2[-1]
res = abs(y - y2)
i = res.argmax()
axvline(h[-1][i], color="grey")
else:
h = histogram(list(zip(*hist))[2], 40, density=True)
plot(h[-1][:-1], h[0], "s-", label="c=%s, init=%s" % (str(c), init))
gca().set_yscale("log")
if cum:
legend(loc="lower right")
else:
......
......@@ -872,7 +872,12 @@ public:
void set_partition(VMap&& b)
{
for (auto v : vertices_range(_g))
move_vertex(v, b[v]);
{
size_t r = b[v];
while (r >= num_vertices(_bg))
add_block();
move_vertex(v, r);
}
}
void set_partition(boost::any& ab)
......@@ -1351,7 +1356,7 @@ public:
template <class MEntries>
double virtual_move(size_t v, size_t r, size_t nr, entropy_args_t ea,
double virtual_move(size_t v, size_t r, size_t nr, const entropy_args_t& ea,
MEntries& m_entries)
{
assert(size_t(_b[v]) == r || r == null_group);
......@@ -1478,7 +1483,7 @@ public:
double propagate_entries_dS(size_t u, size_t v, int du, int dv,
std::vector<std::tuple<size_t, size_t, GraphInterface::edge_t, int,
std::vector<double>>>& entries,
entropy_args_t& ea, std::vector<double>& dBdx,
const entropy_args_t& ea, std::vector<double>& dBdx,
int dL)
{
openmp_scoped_lock lock(_lock);
......@@ -1646,13 +1651,13 @@ public:
return dS;
}
double virtual_move(size_t v, size_t r, size_t nr, entropy_args_t ea)
double virtual_move(size_t v, size_t r, size_t nr, const entropy_args_t& ea)
{
return virtual_move(v, r, nr, ea, _m_entries);
}
double get_delta_partition_dl(size_t v, size_t r, size_t nr,
entropy_args_t& ea)
const entropy_args_t& ea)
{
if (r == nr)
return 0;
......@@ -2024,7 +2029,7 @@ public:
return S;
}
double entropy(entropy_args_t ea, bool propagate=false)
double entropy(const entropy_args_t& ea, bool propagate=false)
{
double S = 0, S_dl = 0;
......@@ -2157,7 +2162,7 @@ public:
}
template <bool Add>
double edge_entropy_term(size_t u, size_t v, entropy_args_t ea)
double edge_entropy_term(size_t u, size_t v, const entropy_args_t& ea)
{
double S = 0, S_dl = 0;
size_t r = _b[u];
......@@ -2343,14 +2348,14 @@ public:
return S + S_dl * ea.beta_dl;
}
double edge_entropy_term(size_t u, size_t v, entropy_args_t ea)
double edge_entropy_term(size_t u, size_t v, const entropy_args_t& ea)
{
return edge_entropy_term<true>(u, v, ea);
}
template <bool Add>
double modify_edge_dS(size_t u, size_t v, GraphInterface::edge_t& e,
const std::vector<double>& recs, entropy_args_t ea)
const std::vector<double>& recs, const entropy_args_t& ea)
{
double dS = 0;
dS -= edge_entropy_term<Add>(u, v, ea);
......@@ -2436,7 +2441,7 @@ public:
disable_partition_stats();
}
void couple_state(BlockStateVirtualBase& s, entropy_args_t ea)
void couple_state(BlockStateVirtualBase& s, const entropy_args_t& ea)
{
_coupled_state = &s;
_coupled_entropy_args = ea;
......@@ -2603,7 +2608,7 @@ public:
EGroups<g_t, is_weighted_t> _egroups;
bool _egroups_enabled = true;
typedef NeighborSampler<g_t, is_weighted_t, is_weighted_t>
typedef NeighborSampler<g_t, is_weighted_t, mpl::false_>
neighbor_sampler_t;
neighbor_sampler_t _neighbor_sampler;
......
......@@ -157,6 +157,7 @@ public:
const pair<size_t, size_t>& get_move() { return _rnr; }
template <bool First, bool Source>
inline __attribute__((always_inline)) __attribute__((flatten))
size_t& get_field_rnr(size_t s, size_t t)
{
auto& out_field = First ? _r_out_field : _nr_out_field;
......@@ -171,6 +172,7 @@ public:
}
}
inline __attribute__((always_inline)) __attribute__((flatten))
size_t& get_field(size_t s, size_t t)
{
if (s == _rnr.first)
......@@ -185,6 +187,7 @@ public:
}
template <bool Add, class... DVals>
inline __attribute__((always_inline)) __attribute__((flatten))
void insert_delta_dispatch(size_t s, size_t t, size_t& f, int d, DVals&&... delta)
{
if (f == _null)
......@@ -196,7 +199,7 @@ public:
_edelta.emplace_back();
}
if (Add)
if constexpr (Add)
{
_delta[f] += d;
tuple_op(_edelta[f], [&](auto&& r, auto&& v){ r += v; },
......@@ -211,7 +214,7 @@ public:
}
template <bool First, bool Source, bool Add, class... DVals>
__attribute__((flatten))
inline __attribute__((always_inline)) __attribute__((flatten))
void insert_delta_rnr(size_t s, size_t t, int d, DVals&&... delta)
{
auto& f = get_field_rnr<First, Source>(s, t);
......@@ -226,6 +229,7 @@ public:
insert_delta_dispatch<Add>(s, t, f, d, std::forward<DVals>(delta)...);
}
inline __attribute__((always_inline)) __attribute__((flatten))
int get_delta(size_t r, size_t s)
{
size_t f = get_field(r, s);
......
......@@ -49,7 +49,7 @@ void export_sbm_state()
void (state_t::*move_vertices)(python::object, python::object) =
&state_t::move_vertices;
double (state_t::*virtual_move)(size_t, size_t, size_t,
entropy_args_t) =
const entropy_args_t&) =
&state_t::virtual_move;
size_t (state_t::*sample_block)(size_t, double, double, rng_t&) =
&state_t::sample_block;
......
......@@ -44,7 +44,7 @@ using namespace std;
((allow_vacate,, bool, 0)) \
((sequential,, bool, 0)) \
((deterministic,, bool, 0)) \
((verbose,, bool, 0)) \
((verbose,, int, 0)) \
((niter,, size_t, 0))
......
......@@ -41,7 +41,7 @@ using namespace std;
((S_max, , double, 0)) \
((f, , double, 0)) \
((S, , double, 0)) \
((verbose,, bool, 0))
((verbose,, int, 0))
template <class State>
......
......@@ -74,7 +74,7 @@ public:
std::vector<size_t>& bmap)
: _bmap(bmap), _N(0), _E(E), _total_B(B)
{
if (!use_rmap)
if constexpr (!use_rmap)
{
_hist.resize(num_vertices(g));
_total.resize(num_vertices(g));
......@@ -110,7 +110,7 @@ public:
size_t get_r(size_t r)
{
if (use_rmap)
if constexpr (use_rmap)
{
constexpr size_t null =
std::numeric_limits<size_t>::max();
......
......@@ -33,13 +33,13 @@ namespace graph_tool
class BlockStateVirtualBase {
public:
virtual double entropy(entropy_args_t eargs, bool propagate) = 0;
virtual double entropy(const entropy_args_t& eargs, bool propagate) = 0;
virtual void add_partition_node(size_t v, size_t r) = 0;
virtual void remove_partition_node(size_t v, size_t r) = 0;
virtual void set_vertex_weight(size_t v, int w) = 0;
virtual void coupled_resize_vertex(size_t v) = 0;
virtual double virtual_move(size_t v, size_t r, size_t nr,
entropy_args_t eargs) = 0;
const entropy_args_t& eargs) = 0;
virtual void sample_branch(size_t v, size_t u, rng_t& rng) = 0;
virtual size_t sample_block(size_t v, double c, double d, rng_t& rng) = 0;
virtual double get_move_prob(size_t v, size_t r, size_t s, double c, double d,
......@@ -55,7 +55,7 @@ public:
const std::vector<double>& rec) = 0;
virtual void remove_edge(size_t u, size_t v, GraphInterface::edge_t& e,
const std::vector<double>& rec) = 0;
virtual double edge_entropy_term(size_t u, size_t v, entropy_args_t ea) = 0;
virtual double edge_entropy_term(size_t u, size_t v, const entropy_args_t& ea) = 0;
virtual void propagate_delta(size_t u, size_t v,
std::vector<std::tuple<size_t, size_t,
GraphInterface::edge_t, int,
......@@ -64,10 +64,10 @@ public:
std::vector<std::tuple<size_t, size_t,
GraphInterface::edge_t, int,
std::vector<double>>>& entries,
entropy_args_t& ea,
const entropy_args_t& ea,
std::vector<double>& dBdx, int dL) = 0;
virtual double get_delta_partition_dl(size_t v, size_t r, size_t nr,
entropy_args_t& ea) = 0;
const entropy_args_t& ea) = 0;
virtual vprop_map_t<int32_t>::type::unchecked_t& get_b() = 0;
virtual bool check_edge_counts(bool emat=true) = 0;
virtual bool allow_move(size_t r, size_t nr) = 0;
......
......@@ -116,7 +116,7 @@ double poisson_w_log_P(DT N, double x, double alpha, double beta)
template <class State>
std::tuple<double,double> rec_entropy(State& state, entropy_args_t& ea)
std::tuple<double,double> rec_entropy(State& state, const entropy_args_t& ea)
{
double S = 0, S_dl = 0;
for (size_t i = 0; i < state._rec_types.size(); ++i)
......@@ -200,7 +200,7 @@ std::tuple<double,double> rec_entropy(State& state, entropy_args_t& ea)
template <class State, class MEntries>
std::tuple<double, double> rec_entries_dS(State& state, MEntries& m_entries,
entropy_args_t& ea,
const entropy_args_t& ea,
std::vector<double>& dBdx, int& dL)
{
int ddL = 0;
......
......@@ -496,8 +496,8 @@ struct Layers
}
template <class MEntries>
double virtual_move(size_t v, size_t r, size_t s, entropy_args_t ea,
MEntries& m_entries)
double virtual_move(size_t v, size_t r, size_t s,
const entropy_args_t& ea, MEntries& m_entries)
{
if (s == r)
{
......@@ -582,7 +582,8 @@ struct Layers
return dS;
}
double virtual_move(size_t v, size_t r, size_t s, entropy_args_t ea)
double virtual_move(size_t v, size_t r, size_t s,
const entropy_args_t& ea)
{
return virtual_move(v, r, s, ea, _m_entries);
}
......@@ -679,7 +680,7 @@ struct Layers
// assert(check_edge_counts());
}
double entropy(entropy_args_t ea, bool propagate=false)
double entropy(const entropy_args_t& ea, bool propagate=false)
{
double S = 0, S_dl = 0;
if (_master)
......@@ -799,7 +800,7 @@ struct Layers
}
}
double edge_entropy_term(size_t, size_t, entropy_args_t) { return 0; }
double edge_entropy_term(size_t, size_t, const entropy_args_t&) { return 0; }
void enable_partition_stats()
{
......@@ -849,7 +850,7 @@ struct Layers
}
void couple_state(LayeredBlockStateVirtualBase& s,
entropy_args_t ea)
const entropy_args_t& ea)
{
_lcoupled_state = &s;
......@@ -880,7 +881,7 @@ struct Layers
}
void couple_state(BlockStateVirtualBase& s,
entropy_args_t ea)
const entropy_args_t& ea)
{
BaseState::couple_state(s, ea);
}
......@@ -1022,7 +1023,7 @@ struct Layers
double propagate_entries_dS(size_t u, size_t v, int du, int dv,
std::vector<std::tuple<size_t, size_t, GraphInterface::edge_t, int,
std::vector<double>>>& entries,
entropy_args_t& ea, std::vector<double>& dBdx,
const entropy_args_t& ea, std::vector<double>& dBdx,
int dL)
{
double dS = BaseState::propagate_entries_dS(u, v, du, dv, entries, ea, dBdx, dL);
......@@ -1044,7 +1045,7 @@ struct Layers
}
double get_delta_partition_dl(size_t v, size_t r, size_t nr,
entropy_args_t& ea)
const entropy_args_t& ea)
{
return BaseState::get_delta_partition_dl(v, r, nr, ea);
}
......
......@@ -56,7 +56,7 @@ void export_lsbm()
{
typedef typename std::remove_reference<decltype(*s)>::type state_t;
double (state_t::*virtual_move)(size_t, size_t, size_t, entropy_args_t) =
double (state_t::*virtual_move)(size_t, size_t, size_t, const entropy_args_t&) =
&state_t::virtual_move;
size_t (state_t::*sample_block)(size_t, double, double, rng_t&)
= &state_t::sample_block;
......@@ -74,7 +74,7 @@ void export_lsbm()
void (state_t::*add_vertices)(python::object, python::object) =
&state_t::add_vertices;
void (state_t::*couple_state)(LayeredBlockStateVirtualBase&,
entropy_args_t) =
const entropy_args_t&) =
&state_t::couple_state;
class_<state_t, bases<LayeredBlockStateVirtualBase>>
......
......@@ -76,7 +76,7 @@ void export_layered_overlap_blockmodel_state()
state_t;
double (state_t::*virtual_move)(size_t, size_t, size_t,
entropy_args_t) =
const entropy_args_t&) =
&state_t::virtual_move;
size_t (state_t::*sample_block)(size_t, double, double,
rng_t&)
......@@ -88,7 +88,7 @@ void export_layered_overlap_blockmodel_state()
python::object) =
&state_t::move_vertices;
void (state_t::*couple_state)(LayeredBlockStateVirtualBase&,
entropy_args_t) =
const entropy_args_t&) =
&state_t::couple_state;
class_<state_t> c(name_demangle(typeid(state_t).name()).c_str(),
......
......@@ -95,7 +95,7 @@ auto mcmc_sweep(MCMCState state, RNG& rng)
if (state.skip_node(v))
continue;
auto r = (state._verbose) ? state.node_state(v)
auto r = (state._verbose > 1) ? state.node_state(v)
: decltype(state.node_state(v))();
move_t s;
......@@ -125,7 +125,7 @@ auto mcmc_sweep(MCMCState state, RNG& rng)
state.step(v, s);
if (state._verbose)
if (state._verbose > 1)
cout << v << ": " << r << " -> " << s << " " << accept << " " << dS << " " << mP << " " << -dS * beta + mP << " " << S << endl;
}
......
......@@ -156,7 +156,7 @@ void export_overlap_blockmodel_state()
typedef typename std::remove_reference<decltype(*s)>::type state_t;
double (state_t::*virtual_move)(size_t, size_t, size_t,
entropy_args_t) =
const entropy_args_t&) =
&state_t::virtual_move;
size_t (state_t::*sample_block)(size_t, double, double, rng_t&)
= &state_t::sample_block;
......
......@@ -56,6 +56,8 @@ typedef mpl::vector2<std::true_type, std::false_type> use_hash_tr;
((candidate_pos,, vmap_t, 0)) \
((bclabel,, vmap_t, 0)) \
((pclabel,, vmap_t, 0)) \
((bfield,, vprop_map_t<std::vector<double>>::type, 0)) \
((Bfield, &, std::vector<double>&, 0)) \
((deg_corr,, bool, 0)) \
((rec_types,, std::vector<int>, 0)) \
((rec,, std::vector<eprop_map_t<double>::type>, 0)) \
......@@ -412,7 +414,7 @@ public:
}
template <class MEntries>
double virtual_move(size_t v, size_t r, size_t nr, entropy_args_t ea,
double virtual_move(size_t v, size_t r, size_t nr, const entropy_args_t& ea,
MEntries& m_entries)
{
if (r == nr)
......@@ -503,12 +505,13 @@ public:
return dS + ea.beta_dl * dS_dl;
}
double virtual_move(size_t v, size_t r, size_t nr, entropy_args_t ea)
double virtual_move(size_t v, size_t r, size_t nr, const entropy_args_t& ea)
{
return virtual_move(v, r, nr, ea, _m_entries);
}
double get_delta_partition_dl(size_t v, size_t r, size_t nr, entropy_args_t& ea)
double get_delta_partition_dl(size_t v, size_t r, size_t nr,
const entropy_args_t& ea)
{
if (r == nr)
return 0;
......@@ -833,7 +836,7 @@ public:
throw GraphException("Dense entropy for overlapping model not implemented!");
}
double entropy(entropy_args_t ea, bool propagate=false)
double entropy(const entropy_args_t& ea, bool propagate=false)
{
double S = 0, S_dl = 0;
if (ea.adjacency)
......@@ -926,7 +929,7 @@ public:
return S;
}
double edge_entropy_term(size_t, size_t, entropy_args_t)
double edge_entropy_term(size_t, size_t, const entropy_args_t&)
{
return 0;
}
......@@ -934,7 +937,7 @@ public:
std::vector<std::tuple<size_t, size_t,
GraphInterface::edge_t, int,
std::vector<double>>>&,
entropy_args_t&, std::vector<double>&, int)
const entropy_args_t&, std::vector<double>&, int)
{
return 0;
}
......@@ -992,7 +995,7 @@ public:
return _partition_stats[_pclabel[v]];
}
void couple_state(BlockStateVirtualBase& s, entropy_args_t ea)
void couple_state(BlockStateVirtualBase& s, const entropy_args_t& ea)
{
_coupled_state = &s;
_coupled_entropy_args = ea;
......
......@@ -44,7 +44,7 @@ using namespace std;
((allow_vacate,, bool, 0)) \
((sequential,, bool, 0)) \
((deterministic,, bool, 0)) \
((verbose,, bool, 0)) \
((verbose,, int, 0)) \
((niter,, size_t, 0))
......
......@@ -87,7 +87,7 @@ inline double lgamma_fast(T x)
{
if (size_t(x) >= __lgamma_cache.size())
{
if (Init)
if constexpr (Init)
init_lgamma(x);
else
return lgamma(x);
......
......@@ -39,7 +39,6 @@ public:
NeighborSampler(Graph& g, Eprop& eweight, bool self_loops=false)
: _g(g),
_sampler(get(vertex_index_t(), g), num_vertices(g)),
_sampler_pos(get(edge_index_t(), g)),
_self_loops(self_loops)
{
init(eweight);
......@@ -48,7 +47,14 @@ public:
template <class Eprop>
void init_vertex(size_t v, Eprop& eweight)
{
_sampler[v].clear();
init_vertex(v, eweight, Weighted());
}
template <class Eprop>
void init_vertex(size_t v, Eprop& eweight, std::true_type)
{
std::vector<item_t> us;
std::vector<double> ps;
for (auto e : out_edges_range(v, _g))
{
......@@ -67,17 +73,17 @@ public:
w /= 2;
}
insert(v, u, w, e);
us.push_back(u);
ps.push_back(w);
}
if constexpr (is_directed_::apply<Graph>::type::value)
{
for (auto e : in_edges_range(v, _g))
{
auto u = source(e, _g);
if (!_self_loops && u == v)
if (u == v)
continue;
auto w = eweight[e];
......@@ -85,7 +91,50 @@ public:
if (w == 0)
continue;
insert(v, u, w, e);
us.push_back(u);
ps.push_back(w);
}
}
_sampler[v] = sampler_t(us, ps);
}
template <class Eprop>
void init_vertex(size_t v, Eprop&, std::false_type)
{
auto& sampler = _sampler[v];
sampler.clear();
[[maybe_unused]] gt_hash_set<size_t> sl_set;
[[maybe_unused]] auto eindex = get(edge_index_t(), _g);
for (auto e : out_edges_range(v, _g))
{
auto u = target(e, _g);
if (u == v)
{
if (!_self_loops)
continue;
if constexpr (!is_directed_::apply<Graph>::type::value)
{
if (sl_set.find(eindex[e]) != sl_set.end())
continue;
sl_set.insert(eindex[e]);
}
}
sampler.push_back(u);
}
if constexpr (is_directed_::apply<Graph>::type::value)
{
for (auto e : in_edges_range(v, _g))
{
auto u = source(e, _g);
if (u == v)
continue;
sampler.push_back(u);
}
}
}
......@@ -101,8 +150,7 @@ public:
vertex_t sample(vertex_t v, RNG& rng)
{
auto& sampler = _sampler[v];
auto& item = sample_item(sampler, rng);
return get_u(item);
return sample_item(sampler, rng);
}
bool empty(vertex_t v)
......@@ -115,54 +163,8 @@ public:
_sampler.resize(n);
}
template <class Edge>
void remove(vertex_t v, vertex_t u, Edge&& e)
{
if (v == u && !_self_loops)
return;
auto& sampler = _sampler[v];
auto& pos = _sampler_pos[e];
bool is_src = (get_src(e) == u);
remove_item({is_src, e}, sampler, pos);
}
template <class Weight, class Edge>
void insert(vertex_t v, vertex_t u, Weight w, Edge&& e)
{
if (v == u && !_self_loops)
return;
auto& sampler = _sampler[v];
auto& pos = _sampler_pos[e];
bool is_src = (get_src(e) == u);
insert_item({is_src, e}, w, sampler, pos);
}
private:
typedef std::tuple<bool, edge_t> item_t;
vertex_t get_src(const edge_t& e)
{
if constexpr (is_directed_::apply<Graph>::type::value)
return source(e, _g);
else
return std::min(source(e, _g), target(e, _g));
}
vertex_t get_tgt(const edge_t& e)
{
if constexpr (is_directed_::apply<Graph>::type::value)
return target(e, _g);
else
return std::max(source(e, _g), target(e, _g));
}
vertex_t get_u(const item_t& item)
{
if (get<0>(item))
return get_src(get<1>(item));
else
return get_tgt(get<1>(item));
}
typedef vertex_t item_t;
template <class RNG>
const item_t& sample_item(std::vector<item_t>& sampler, RNG& rng)
......@@ -176,55 +178,6 @@ private:
return sampler.sample(rng);
}
size_t& get_pos(const item_t& u, std::tuple<size_t, size_t>& pos)
{
if (get<0>(u))
return get<0>(pos);
else
return get<1>(pos);
}
void remove_item(const item_t& u, std::vector<item_t>& sampler,
std::tuple<size_t, size_t>& pos)
{
auto u_pos = get_pos(u, pos);
if (u_pos >= sampler.size() || sampler[u_pos] != u)
return;
auto& back = sampler.back();
auto& e = get<1>(back);
auto& bpos = _sampler_pos[e];
get_pos(back, bpos) = u_pos;
sampler[u_pos] = back;
sampler.pop_back();
}
template <class Sampler>
void remove_item(const item_t& u, Sampler& sampler,
std::tuple<size_t, size_t>& pos)
{
auto i = get_pos(u, pos);
if (!sampler.is_valid(i) || sampler[i] != u)
return;
sampler.remove(i);
}
template <class Weight>
void insert_item(const item_t& u, Weight, std::vector<item_t>& sampler,
std::tuple<size_t, size_t>& pos)
{
get_pos(u, pos) = sampler.size();
sampler.push_back(u);
}
template <class Weight>
void insert_item(const item_t& u, Weight w, DynamicSampler<item_t>& sampler,
std::tuple<size_t, size_t>& pos)
{
get_pos(u, pos) = sampler.insert(u, w);
}
Graph& _g;
typedef typename std::conditional<Weighted::value,
......@@ -232,15 +185,12 @@ private:
DynamicSampler<item_t>,
Sampler<item_t,
boost::mpl::false_>>::type,
vector<item_t>>::type
std::vector<item_t>>::type
sampler_t;
typedef typename vprop_map_t<sampler_t>::type::unchecked_t vsampler_t;
vsampler_t _sampler;
typedef typename eprop_map_t<std::tuple<size_t, size_t>>::type pos_map_t;
pos_map_t _sampler_pos;
bool _self_loops;
};
......
......@@ -36,7 +36,7 @@ double log_q(T n, T k)
return 0;
if (k > n)
k = n;
if (n < T(__q_cache.shape()[0]))
if (size_t(n) < __q_cache.shape()[0])
return __q_cache[n][k];
return log_q_approx(n, k);
}
......
......@@ -34,12 +34,12 @@ using namespace std;
typedef std::vector<size_t> vlist_t;
#define MCMC_EPIDEMICS_STATE_params(State) \
#define MCMC_EPIDEMICS_STATE_params(State) \
((__class__,&, mpl::vector<python::object>, 1)) \
((state, &, State&, 0)) \
((beta,, double, 0)) \
((hstep,, double, 0)) \
((verbose,, bool, 0)) \
((verbose,, int, 0) ) \
((niter,, size_t, 0))
......
......@@ -66,7 +66,7 @@ typedef std::vector<size_t> vlist_t;
((xlog,, bool, 0)) \
((xdefault,, double, 0)) \
((entropy_args,, uentropy_args_t, 0)) \
((verbose,, bool, 0)) \
((verbose,, int, 0)) \
((niter,, size_t, 0))
......
......@@ -39,7 +39,7 @@ typedef std::vector<size_t> vlist_t;
((beta,, double, 0)) \
((n,, size_t, 0)) \
((hstep,, double, 0)) \
((verbose,, bool, 0)) \
((verbose,, int, 0)) \
((niter,, size_t, 0))
......
......@@ -40,7 +40,7 @@ typedef std::vector<size_t> vlist_t;
((beta,, double, 0)) \
((entropy_args,, uentropy_args_t, 0)) \
((edges_only,, bool, 0)) \
((verbose,, bool, 0)) \
((verbose,, int, 0)) \
((niter,, size_t, 0))
......
......@@ -24,6 +24,7 @@
#include "hash_map_wrap.hh"
#include "coroutine.hh"
#include "graph_python_interface.hh"
#include "../generation/sampler.hh"
#include <boost/graph/breadth_first_search.hpp>
#include <boost/graph/dijkstra_shortest_paths_no_color_map.hpp>
......@@ -595,6 +596,75 @@ python::object do_get_all_shortest_paths(GraphInterface& gi, size_t s, size_t t,
#endif // HAVE_BOOST_COROUTINE
}
void get_weighted_succs(size_t t, boost::any apred, boost::any asucc,
boost::any acount, boost::any avisited)
{
typedef vprop_map_t<vector<int64_t>>::type pred_map_t;
typedef vprop_map_t<int64_t>::type count_map_t;
typedef vprop_map_t<uint8_t>::type visited_map_t;
pred_map_t pred = any_cast<pred_map_t>(apred);
pred_map_t succ = any_cast<pred_map_t>(asucc);
count_map_t count = any_cast<count_map_t>(acount);
visited_map_t visited = any_cast<visited_map_t>(avisited);
count[t] = 1;
visited[t] = true;
deque<size_t> queue = {t};
while (!queue.empty())
{
size_t v = queue.front();
queue.pop_front();
for (auto w : pred[v])
{
count[w] += count[v];
succ[w].push_back(v);
if (!visited[w])
{
visited[w] = true;
queue.push_back(w);
}
}
}
};
void get_random_shortest_path(size_t s, size_t t, boost::any asucc,
boost::any acount, vector<size_t>& path,
rng_t& rng)
{
typedef vprop_map_t<vector<int64_t>>::type pred_map_t;
typedef vprop_map_t<int64_t>::type count_map_t;
pred_map_t succ = any_cast<pred_map_t>(asucc);
count_map_t count = any_cast<count_map_t>(acount);
vector<double> probs;
path.clear();
if (succ[s].empty())
return;
path.push_back(s);
size_t v = s;
while (v != t)
{
if (succ[v].size() == 1)
{
v = succ[v].front();
}
else
{
probs.clear();
for (auto w : succ[v])
probs.push_back(count[w]);
Sampler<int64_t> next(succ[v], probs);
v = next.sample(rng);
}
path.push_back(v);
}
};
template <bool edges, class Graph, class Yield, class VMap>
void get_all_paths(size_t s, size_t t, size_t cutoff, VMap visited,
......@@ -691,4 +761,7 @@ void export_dists()
python::def("get_all_preds", &do_get_all_preds);
python::def("get_all_shortest_paths", &do_get_all_shortest_paths);
python::def("get_all_paths", &do_get_all_paths);
python::def("get_weighted_succs", &get_weighted_succs);
python::def("get_random_shortest_path", &get_random_shortest_path);
};
......@@ -1499,13 +1499,13 @@ class PropertyDict(object):
def __len__(self):
count = 0
for k in self.properties.iterkeys():
for k in self.properties.keys():
if k[0] == self.t:
count += 1
return count
def __iter__(self):
return self.iterkeys()
return self.keys()
def iterkeys(self):
for k in self.properties.iterkeys():
......@@ -1531,23 +1531,23 @@ class PropertyDict(object):
if k[0] == self.t:
yield v
def keys(self):
return [k[1] for k in self.properties.keys() if k[0] == self.t]
if sys.version_info < (3,):
def keys(self):
return list(self.iterkeys())
def values(self):
return [v for k, v in self.properties.iteritems() if k[0] == self.t]
return list(self.itervalues())
def __repr__(self):
temp = dict([(k[1], v) for k, v in self.properties.iteritems() if k[0] == self.t])
return repr(temp)
else:
def keys(self):
return self.iterkeys()
def values(self):
return [v for k, v in self.properties.items() if k[0] == self.t]
return self.itervalues()
def __repr__(self):
temp = dict([(k[1], v) for k, v in self.properties.items() if k[0] == self.t])
return repr(temp)
def __getattr__(self, attr):
return self.__getitem__(attr)
......@@ -1761,6 +1761,19 @@ class Graph(object):
# provide more useful information
d = "directed" if self.is_directed() else "undirected"
fr = ", reversed" if self.is_reversed() and self.is_directed() else ""
p = []
if len(self.vp) == 1:
p.append("%d internal vertex %s" % (len(self.vp),
"property" if len(self.vp) == 1 else "properties"))
if len(self.ep) == 1:
p.append("%d internal edge %s" % (len(self.ep),
"property" if len(self.ep) == 1 else "properties"))
if len(self.gp) == 1:
p.append("%d internal graph %s" % (len(self.gp),
"property" if len(self.gp) == 1 else "properties"))
p = ", ".join(p)
if len(p) > 0:
p = ", " + p
f = ""
if self.get_edge_filter()[0] is not None:
f += ", edges filtered by %s" % (str(self.get_edge_filter()))
......@@ -1768,10 +1781,10 @@ class Graph(object):
f += ", vertices filtered by %s" % (str(self.get_vertex_filter()))
n = self.num_vertices()
e = self.num_edges()
return "<%s object, %s%s, with %d %s and %d edge%s%s at 0x%x>"\
return "<%s object, %s%s, with %d %s and %d edge%s%s%s, at 0x%x>"\
% (type(self).__name__, d, fr, n,
"vertex" if n == 1 else "vertices", e, "" if e == 1 else "s",
f, id(self))
p, f, id(self))
# Graph access
# ============
......
......@@ -367,6 +367,33 @@ class GraphWidget(Gtk.DrawingArea):
self.is_rotating = False
self.is_drag_gesture = False
def update(self, pos=None, vprops=None, eprops=None, vorder=None, eorder=None,
nodesfirst=None, display_props=None,
fit_view=True, bg_color=None, **kwargs):
props, kwargs = parse_props("vertex", kwargs)
vprops.update(props)
props, kwargs = parse_props("edge", kwargs)
eprops.update(props)
if pos is not None:
self.pos = pos
self.vprops.update(vprops)
self.eprops.update(eprops)
if vorder is not None:
self.vorder = vorder
if eorder is not None:
self.eorder = eorder
if nodesfirst is not None:
self.nodesfirst = nodesfirst
self.fit_view = fit_view
self.display_prop = self.g.vertex_index if display_props is None \
else display_props
self.bg_color = bg_color if bg_color is not None else [1, 1, 1, 1]
def cleanup(self):
"""Cleanup callbacks."""
if gobject is None:
......@@ -1187,9 +1214,9 @@ _window_list = []
def interactive_window(g, pos=None, vprops=None, eprops=None, vorder=None,
eorder=None, nodesfirst=False, geometry=(500, 400),
update_layout=True, sync=True, main=True, **kwargs):
r"""
Display an interactive GTK+ window containing the given graph.
update_layout=True, sync=True, main=True,
window=None, return_window=False, **kwargs):
r"""Display an interactive GTK+ window containing the given graph.
Parameters
----------
......@@ -1216,10 +1243,13 @@ def interactive_window(g, pos=None, vprops=None, eprops=None, vorder=None,
Window geometry.
update_layout : bool (optional, default: ``True``)
If ``True``, the layout will be updated dynamically.
sync : bool (optional, default: ``True``)
If ``False``, run asynchronously. (Requires :mod:`IPython`)
main : bool (optional, default: ``True``)
If ``False``, the GTK+ main loop will not be called.
window : :class:`~graph_tool.draw.GraphWindow` (optional, default: ``None``)
If provided, specifies the window where the drawing will
occur. Otherwise a new one will be created.
return_window : bool (optional, default: ``False``)
If ``True``, the GTK+ window will be returned.
**kwargs
Any extra parameters are passed to :class:`~graph_tool.draw.GraphWindow`,
:class:`~graph_tool.draw.GraphWidget` and :func:`~graph_tool.draw.cairo_draw`.
......@@ -1239,28 +1269,35 @@ def interactive_window(g, pos=None, vprops=None, eprops=None, vorder=None,
information.
"""
if pos is None:
if update_layout:
pos = random_layout(g, [1, 1])
else:
pos = sfdp_layout(g)
win = GraphWindow(g, pos, geometry, vprops, eprops, vorder, eorder,
nodesfirst, update_layout, **kwargs)
win.show_all()
_window_list.append(win)
if window is None:
if pos is None:
if update_layout:
pos = random_layout(g, [1, 1])
else:
pos = sfdp_layout(g)
win = GraphWindow(g, pos, geometry, vprops, eprops, vorder, eorder,
nodesfirst, update_layout, **kwargs)
win.show_all()
_window_list.append(win)
else:
win = window
win.graph.update(pos, vprops, eprops, vorder, eorder, nodesfirst,
**kwargs)
win.graph.regenerate_surface(complete=True)
win.graph.queue_draw()
if main:
if not sync:
# just a placeholder for a proper main loop integration with gtk3 when
# ipython implements it
import IPython.lib.inputhook
f = lambda: Gtk.main_iteration_do(False)
IPython.lib.inputhook.set_inputhook(f)
else:
def destroy_callback(*args, **kwargs):
global _window_list
for w in _window_list:
w.destroy()
Gtk.main_quit()
win.connect("delete_event", destroy_callback)
Gtk.main()
def destroy_callback(*args, **kwargs):
global _window_list
for w in _window_list:
w.destroy()
Gtk.main_quit()
win.connect("delete_event", destroy_callback)
Gtk.main()
else:
while Gtk.events_pending():
Gtk.main_iteration()
if return_window:
return win
return pos, win.graph.selected.copy()
......@@ -23,9 +23,9 @@ import sys
if sys.version_info < (3,):
range = xrange
from .. import _degree, _prop, Graph, GraphView, libcore, _get_rng, PropertyMap, \
conv_pickle_state, Vector_size_t, Vector_double, group_vector_property, \
perfect_prop_hash
from .. import _degree, _prop, Graph, GraphView, libcore, _get_rng, \
PropertyMap, VertexPropertyMap, conv_pickle_state, Vector_size_t, \
Vector_double, group_vector_property, perfect_prop_hash
from .. generation import condensation_graph, random_rewire, generate_sbm, \
solve_sbm_fugacities, generate_maxent_sbm
from .. stats import label_self_loops, remove_parallel_edges, remove_self_loops
......@@ -874,8 +874,10 @@ class BlockState(object):
r"""Returns the property map which contains the block labels for each vertex."""
return self.b
def set_blocks(self, b):
def set_state(self, b):
r"""Sets the internal partition of the state."""
if not isinstance(b, VertexPropertyMap):
b = self.g.new_vp("int32_t", vals=b)
if b.value_type() != "int32_t":
b = b.copy("int32_t")
self._state.set_partition(_prop("v", self.g, b))
......@@ -1493,7 +1495,7 @@ class BlockState(object):
_get_rng())
def multiflip_mcmc_sweep(self, beta=1., c=1., a1=.95, psplit=.5,
def multiflip_mcmc_sweep(self, beta=1., c=1., a1=.95, psplit=.5, prec=.3,
d=0.01, gibbs_sweeps=10, niter=1, entropy_args={},
verbose=False, **kwargs):
r"""Perform ``niter`` sweeps of a Metropolis-Hastings acceptance-rejection
......@@ -1515,6 +1517,8 @@ class BlockState(object):
psplit : ``float`` (optional, default: ``.5``)
Probability of proposing a group split. A group merge will be
proposed with probability ``1-psplit``.
prec : ``float`` (optional, default: ``.3``)
Probability of proposing a group recombination.
d : ``float`` (optional, default: ``.01``)
Probability of selecting a new (i.e. empty) group for a given
single-node move.
......
......@@ -613,6 +613,8 @@ class TemperingState(object):
if idx is None:
self.idx = list(range(len(betas)))
self.beta_dl = beta_dl
self.swap_attempts = 0
self.swaps = numpy.zeros(len(states), dtype="int")
def entropy(self, **kwargs):
"""Returns the sum of the entropy of the parallel states. All keyword
......@@ -638,7 +640,7 @@ class TemperingState(object):
return [s.entropy(**kwargs) * beta for s, beta in
zip(self.states, self.betas)]
def states_swap(self, adjacent=True, **kwargs):
def states_swap(self, **kwargs):
"""Perform a full sweep of the parallel states, where swaps are attempted. All
relevant keyword arguments are propagated to the individual states'
`entropy()` method."""
......@@ -646,20 +648,12 @@ class TemperingState(object):
verbose = kwargs.get("verbose", False)
eargs = kwargs.get("entropy_args", {})
if adjacent:
idx = numpy.arange(len(self.states) - 1)
else:
idx = numpy.arange(len(self.states))
idx = numpy.arange(len(self.states) - 1)
numpy.random.shuffle(idx)
nswaps = 0
dS = 0
for i in idx:
if adjacent:
j = i + 1
else:
j = i
while j == i:
j = numpy.random.randint(0, len(idx))
j = i + 1
s1 = self.states[i]
s2 = self.states[j]
......@@ -673,18 +667,22 @@ class TemperingState(object):
P1_b = -s1.entropy(beta_dl=b1, **eargs)
P2_b = -s2.entropy(beta_dl=b2, **eargs)
else:
P1_f = -s1.entropy(**eargs) * b2
P2_f = -s2.entropy(**eargs) * b1
S1 = s1.entropy(**eargs)
S2 = s2.entropy(**eargs)
P1_f = -S1 * b2
P2_f = -S2 * b1
P1_b = -s1.entropy(**eargs) * b1
P2_b = -s2.entropy(**eargs) * b2
P1_b = -S1 * b1
P2_b = -S2 * b2
ddS = -(P1_f + P2_f - P1_b - P2_b)
self.swap_attempts += 1
if ddS < 0 or numpy.random.random() < exp(-ddS):
self.states[j], self.states[i], self.idx[j], self.idx[i] = \
self.states[i], self.states[j], self.idx[i], self.idx[j]
nswaps += 1
self.swaps[i] += 1
dS += ddS
if check_verbose(verbose):
print(verbose_pad(verbose)
......
......@@ -215,6 +215,11 @@ class NestedBlockState(object):
"""
return [s.b.fa for s in self.levels]
def set_state(self, bs):
r"""Sets the internal nested partition of the state."""
for i in range(len(bs)):
self.levels[i].set_state(bs[i])
def get_levels(self):
"""Get hierarchy levels as a list of :class:`~graph_tool.inference.blockmodel.BlockState`
instances."""
......