Commit 1c48aec8 authored by Tiago Peixoto's avatar Tiago Peixoto

inference: Return number of attempted moves for mcmc sweeps

parent 1c155fb7
......@@ -144,8 +144,8 @@ for directed in [True, False]:
continue
except TypeError:
pass
Ss1 = array(list(zip(*hists[c1]))[0])
Ss2 = array(list(zip(*hists[c2]))[0])
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))
......
......@@ -49,8 +49,7 @@ python::object do_gibbs_sweep(python::object ogibbs_state,
[&](auto& s)
{
auto ret_ = gibbs_sweep(s, rng);
ret = python::make_tuple(get<0>(ret_), get<1>(ret_),
get<2>(ret_));
ret = tuple_apply([&](auto&... args){ return python::make_tuple(args...); }, ret_);
});
};
block_state::dispatch(oblock_state, dispatch);
......
......@@ -51,12 +51,12 @@ python::object do_mcmc_sweep(python::object omcmc_state,
if (s._parallel)
{
auto ret_ = mcmc_sweep_parallel(s, rng);
ret = python::make_tuple(ret_.first, ret_.second);
ret = tuple_apply([&](auto&... args){ return python::make_tuple(args...); }, ret_);
}
else
{
auto ret_ = mcmc_sweep(s, rng);
ret = python::make_tuple(ret_.first, ret_.second);
ret = tuple_apply([&](auto&... args){ return python::make_tuple(args...); }, ret_);
}
});
};
......
......@@ -49,7 +49,7 @@ python::object do_merge_sweep(python::object omerge_state,
[&](auto& s)
{
auto ret_ = merge_sweep(s, rng);
ret = python::make_tuple(ret_.first, ret_.second);
ret = tuple_apply([&](auto&... args){ return python::make_tuple(args...); }, ret_);
});
};
block_state::dispatch(oblock_state, dispatch);
......
......@@ -64,7 +64,7 @@ python::object do_multicanonical_sweep(python::object omulticanonical_state,
[&](auto& mc_state)
{
auto ret_ = mcmc_sweep(mc_state, rng);
ret = python::make_tuple(ret_.first, ret_.second);
ret = tuple_apply([&](auto&... args){ return python::make_tuple(args...); }, ret_);
});
});
};
......
......@@ -64,7 +64,7 @@ python::object do_multicanonical_multiflip_sweep(python::object omulticanonical_
[&](auto& mc_state)
{
auto ret_ = mcmc_sweep(mc_state, rng);
ret = python::make_tuple(ret_.first, ret_.second);
ret = tuple_apply([&](auto&... args){ return python::make_tuple(args...); }, ret_);
});
});
};
......
......@@ -49,7 +49,7 @@ python::object do_multiflip_mcmc_sweep(python::object omcmc_state,
[&](auto& s)
{
auto ret_ = mcmc_sweep(s, rng);
ret = python::make_tuple(ret_.first, ret_.second);
ret = tuple_apply([&](auto&... args){ return python::make_tuple(args...); }, ret_);
});
};
block_state::dispatch(oblock_state, dispatch);
......
......@@ -74,9 +74,7 @@ struct MCMC
_vpos(get(vertex_index_t(), _state._g),
num_vertices(_state._g)),
_rpos(get(vertex_index_t(), _state._bg),
num_vertices(_state._bg)),
_sequential(false),
_deterministic(false)
num_vertices(_state._bg))
{
_state.init_mcmc(_c,
(_entropy_args.partition_dl ||
......@@ -102,8 +100,7 @@ struct MCMC
std::vector<size_t> _vlist;
std::vector<size_t> _vs;
bool _sequential;
bool _deterministic;
size_t _null_move = null_group;
size_t node_state(size_t r)
{
......@@ -267,12 +264,12 @@ struct MCMC
bool is_deterministic()
{
return _deterministic;
return false;
}
bool is_sequential()
{
return _sequential;
return false;
}
auto& get_vlist()
......
......@@ -62,8 +62,7 @@ python::object gibbs_layered_sweep(python::object ogibbs_state,
[&](auto& s)
{
auto ret_ = gibbs_sweep(s, rng);
ret = python::make_tuple(get<0>(ret_), get<1>(ret_),
get<2>(ret_));
ret = tuple_apply([&](auto&... args){ return python::make_tuple(args...); }, ret_);
});
},
false);
......
......@@ -64,12 +64,12 @@ python::object mcmc_layered_sweep(python::object omcmc_state,
if (s._parallel)
{
auto ret_ = mcmc_sweep_parallel(s, rng);
ret = python::make_tuple(ret_.first, ret_.second);
ret = tuple_apply([&](auto&... args){ return python::make_tuple(args...); }, ret_);
}
else
{
auto ret_ = mcmc_sweep(s, rng);
ret = python::make_tuple(ret_.first, ret_.second);
ret = tuple_apply([&](auto&... args){ return python::make_tuple(args...); }, ret_);
}
});
},
......
......@@ -62,7 +62,7 @@ python::object merge_layered_sweep(python::object omerge_state,
[&](auto& s)
{
auto ret_ = merge_sweep(s, rng);
ret = python::make_tuple(ret_.first, ret_.second);
ret = tuple_apply([&](auto&... args){ return python::make_tuple(args...); }, ret_);
});
},
false);
......
......@@ -78,7 +78,7 @@ python::object multicanonical_layered_sweep(python::object omulticanonical_state
[&](auto& mc_state)
{
auto ret_ = mcmc_sweep(mc_state, rng);
ret = python::make_tuple(ret_.first, ret_.second);
ret = tuple_apply([&](auto&... args){ return python::make_tuple(args...); }, ret_);
});
});
},
......
......@@ -78,7 +78,7 @@ python::object multicanonical_layered_multiflip_sweep(python::object omulticanon
[&](auto& mc_state)
{
auto ret_ = mcmc_sweep(mc_state, rng);
ret = python::make_tuple(ret_.first, ret_.second);
ret = tuple_apply([&](auto&... args){ return python::make_tuple(args...); }, ret_);
});
});
},
......
......@@ -62,7 +62,7 @@ python::object multiflip_mcmc_layered_sweep(python::object omcmc_state,
[&](auto& s)
{
auto ret_ = mcmc_sweep(s, rng);
ret = python::make_tuple(ret_.first, ret_.second);
ret = tuple_apply([&](auto&... args){ return python::make_tuple(args...); }, ret_);
});
},
false);
......
......@@ -62,8 +62,7 @@ python::object gibbs_layered_overlap_sweep(python::object ogibbs_state,
[&](auto& s)
{
auto ret_ = gibbs_sweep(s, rng);
ret = python::make_tuple(get<0>(ret_), get<1>(ret_),
get<2>(ret_));
ret = tuple_apply([&](auto&... args){ return python::make_tuple(args...); }, ret_);
});
},
false);
......
......@@ -64,12 +64,12 @@ python::object mcmc_layered_overlap_sweep(python::object omcmc_state,
if (s._parallel)
{
auto ret_ = mcmc_sweep_parallel(s, rng);
ret = python::make_tuple(ret_.first, ret_.second);
ret = tuple_apply([&](auto&... args){ return python::make_tuple(args...); }, ret_);
}
else
{
auto ret_ = mcmc_sweep(s, rng);
ret = python::make_tuple(ret_.first, ret_.second);
ret = tuple_apply([&](auto&... args){ return python::make_tuple(args...); }, ret_);
}
});
},
......
......@@ -63,7 +63,7 @@ python::object mcmc_layered_overlap_bundled_sweep(python::object omcmc_state,
[&](auto& s)
{
auto ret_ = mcmc_sweep(s, rng);
ret = python::make_tuple(ret_.first, ret_.second);
ret = tuple_apply([&](auto&... args){ return python::make_tuple(args...); }, ret_);
});
},
false);
......
......@@ -78,7 +78,7 @@ multicanonical_layered_overlap_sweep(python::object omulticanonical_state,
[&](auto& mc_state)
{
auto ret_ = mcmc_sweep(mc_state, rng);
ret = python::make_tuple(ret_.first, ret_.second);
ret = tuple_apply([&](auto&... args){ return python::make_tuple(args...); }, ret_);
});
});
},
......
......@@ -78,7 +78,7 @@ multicanonical_layered_overlap_multiflip_sweep(python::object omulticanonical_st
[&](auto& mc_state)
{
auto ret_ = mcmc_sweep(mc_state, rng);
ret = python::make_tuple(ret_.first, ret_.second);
ret = tuple_apply([&](auto&... args){ return python::make_tuple(args...); }, ret_);
});
});
},
......
......@@ -62,7 +62,7 @@ python::object multiflip_mcmc_layered_overlap_sweep(python::object omcmc_state,
[&](auto& s)
{
auto ret_ = mcmc_sweep(s, rng);
ret = python::make_tuple(ret_.first, ret_.second);
ret = tuple_apply([&](auto&... args){ return python::make_tuple(args...); }, ret_);
});
},
false);
......
......@@ -63,7 +63,7 @@ python::object vacate_layered_overlap_sweep(python::object ovacate_state,
[&](auto& s)
{
auto ret_ = bundled_vacate_sweep(s, rng);
ret = python::make_tuple(ret_.first, ret_.second);
ret = tuple_apply([&](auto&... args){ return python::make_tuple(args...); }, ret_);
});
},
false);
......
......@@ -37,7 +37,9 @@ template <class MergeState, class RNG>
auto bundled_vacate_sweep(MergeState& state, RNG& rng)
{
if (state._nmerges == 0)
return make_pair(double(0), size_t(0));
return std::make_tuple(double(0), size_t(0), size_t(0));
size_t nattempts = 0;
// individual bundles can move in different directions
auto get_best_move = [&] (auto& bundle, auto& past_moves)
......@@ -64,6 +66,7 @@ auto bundled_vacate_sweep(MergeState& state, RNG& rng)
get<1>(best_move) = dS;
}
}
nattempts += state._niter;
};
find_candidates(false);
......@@ -115,6 +118,7 @@ auto bundled_vacate_sweep(MergeState& state, RNG& rng)
get<1>(best_move) = dS;
}
}
nattempts += state._niter;
};
find_candidates(false);
......@@ -244,7 +248,7 @@ auto bundled_vacate_sweep(MergeState& state, RNG& rng)
nmerges++;
}
return make_pair(S, nmerges);
return std::make_tuple(S, nattempts, nmerges);
}
} // graph_tool namespace
......
......@@ -175,7 +175,7 @@ auto gibbs_sweep(GibbsState state, RNG& rng_)
std::reverse(vlist.begin(), vlist.end());
}
return std::make_tuple(S, nmoves, nattempts);
return std::make_tuple(S, nattempts, nmoves);
}
} // graph_tool namespace
......
......@@ -65,6 +65,7 @@ auto mcmc_sweep(MCMCState state, RNG& rng)
auto beta = state.get_beta();
double S = 0;
size_t nattempts = 0;
size_t nmoves = 0;
for (size_t iter = 0; iter < state.get_niter(); ++iter)
......@@ -92,6 +93,8 @@ auto mcmc_sweep(MCMCState state, RNG& rng)
double dS, mP;
std::tie(dS, mP) = state.virtual_move_dS(v, s);
nattempts += state.node_weight(v);
if (metropolis_accept(dS, mP, beta, rng))
{
state.perform_move(v, s);
......@@ -108,7 +111,7 @@ auto mcmc_sweep(MCMCState state, RNG& rng)
if (state.is_sequential() && state.is_deterministic())
std::reverse(vlist.begin(), vlist.end());
}
return make_pair(S, nmoves);
return make_tuple(S, nattempts, nmoves);
}
......@@ -127,9 +130,9 @@ auto mcmc_sweep_parallel(MCMCState state, RNG& rng_)
auto& vlist = state._vlist;
auto& beta = state._beta;
double S = 0;
size_t nmoves = 0;
size_t nattempts = 0;
for (size_t iter = 0; iter < state._niter; ++iter)
{
......@@ -175,8 +178,10 @@ auto mcmc_sweep_parallel(MCMCState state, RNG& rng_)
cout << v << ": " << r << " -> " << s << " " << S << endl;
});
size_t nattempts = 0;
for (auto v : vlist)
{
nattempts++;
auto s = best_move[v].first;
double dS = best_move[v].second;
if (dS != numeric_limits<double>::max())
......@@ -192,7 +197,7 @@ auto mcmc_sweep_parallel(MCMCState state, RNG& rng_)
}
}
}
return make_pair(S, nmoves);
return std::make_tuple(S, nattempts, nmoves);
}
......
......@@ -52,7 +52,10 @@ auto merge_sweep(MergeState state, RNG& rng_)
make_tuple(size_t(0), size_t(0),
numeric_limits<double>::max()));
#pragma omp parallel firstprivate(state) if (state._parallel)
size_t nattempts = 0;
#pragma omp parallel firstprivate(state) reduction(+:nattempts) \
if (state._parallel)
parallel_loop_no_spawn
(state._available,
[&](size_t, auto v)
......@@ -78,6 +81,7 @@ auto merge_sweep(MergeState state, RNG& rng_)
if (dS < get<2>(best_merge[v]))
best_merge[v] = make_tuple(v, s, dS);
}
nattempts += state._niter;
};
find_candidates(false);
......@@ -129,7 +133,7 @@ auto merge_sweep(MergeState state, RNG& rng_)
// collapse merge tree
state.finalize();
return make_pair(S, nmerges);
return std::make_tuple(S, nattempts, nmerges);
}
} // graph_tool namespace
......
......@@ -53,8 +53,7 @@ python::object gibbs_overlap_sweep(python::object ogibbs_state,
[&](auto& s)
{
auto ret_ = gibbs_sweep(s, rng);
ret = python::make_tuple(get<0>(ret_), get<1>(ret_),
get<2>(ret_));
ret = tuple_apply([&](auto&... args){ return python::make_tuple(args...); }, ret_);
});
};
overlap_block_state::dispatch(oblock_state, dispatch);
......
......@@ -53,7 +53,7 @@ python::object overlap_mcmc_sweep(python::object omcmc_state,
[&](auto& s)
{
auto ret_ = mcmc_sweep(s, rng);
ret = python::make_tuple(ret_.first, ret_.second);
ret = tuple_apply([&](auto&... args){ return python::make_tuple(args...); }, ret_);
});
};
overlap_block_state::dispatch(oblock_state, dispatch);
......
......@@ -50,7 +50,7 @@ python::object do_overlap_mcmc_bundled_sweep(python::object omcmc_state,
[&](auto& s)
{
auto ret_ = mcmc_sweep(s, rng);
ret = python::make_tuple(ret_.first, ret_.second);
ret = tuple_apply([&](auto&... args){ return python::make_tuple(args...); }, ret_);
});
};
overlap_block_state::dispatch(oblock_state, dispatch);
......
......@@ -64,7 +64,7 @@ python::object multicanonical_overlap_sweep(python::object omulticanonical_state
[&](auto& mc_state)
{
auto ret_ = mcmc_sweep(mc_state, rng);
ret = python::make_tuple(ret_.first, ret_.second);
ret = tuple_apply([&](auto&... args){ return python::make_tuple(args...); }, ret_);
});
});
};
......
......@@ -64,7 +64,7 @@ python::object multicanonical_overlap_multiflip_sweep(python::object omulticanon
[&](auto& mc_state)
{
auto ret_ = mcmc_sweep(mc_state, rng);
ret = python::make_tuple(ret_.first, ret_.second);
ret = tuple_apply([&](auto&... args){ return python::make_tuple(args...); }, ret_);
});
});
};
......
......@@ -53,7 +53,7 @@ python::object overlap_multiflip_mcmc_sweep(python::object omcmc_state,
[&](auto& s)
{
auto ret_ = mcmc_sweep(s, rng);
ret = python::make_tuple(ret_.first, ret_.second);
ret = tuple_apply([&](auto&... args){ return python::make_tuple(args...); }, ret_);
});
};
overlap_block_state::dispatch(oblock_state, dispatch);
......
......@@ -50,7 +50,7 @@ python::object vacate_overlap_sweep(python::object omerge_state,
[&](auto& s)
{
auto ret_ = bundled_vacate_sweep(s, rng);
ret = python::make_tuple(ret_.first, ret_.second);
ret = tuple_apply([&](auto&... args){ return python::make_tuple(args...); }, ret_);
});
};
overlap_block_state::dispatch(oblock_state, dispatch);
......
......@@ -1447,6 +1447,8 @@ class BlockState(object):
-------
dS : ``float``
Entropy difference after the sweeps.
nattempts : ``int``
Number of vertex moves attempted.
nmoves : ``int``
Number of vertices moved.
......@@ -1490,7 +1492,7 @@ class BlockState(object):
assert self._check_clabel(), "invalid clabel before sweep"
Si = self.entropy(**entropy_args)
try:
dS, nmoves = self._mcmc_sweep_dispatch(mcmc_state)
dS, nattempts, nmoves = self._mcmc_sweep_dispatch(mcmc_state)
finally:
self.B = max(int(self.b.fa.max()) + 1, self.B)
......@@ -1509,7 +1511,7 @@ class BlockState(object):
if not dispatch:
return mcmc_state
return dS, nmoves
return dS, nattempts, nmoves
def _multiflip_mcmc_sweep_dispatch(self, mcmc_state):
return libinference.multiflip_mcmc_sweep(mcmc_state, self._state,
......@@ -1555,6 +1557,8 @@ class BlockState(object):
-------
dS : ``float``
Entropy difference after the sweeps.
nattempts : ``int``
Number of vertex moves attempted.
nmoves : ``int``
Number of vertices moved.
......@@ -1592,15 +1596,11 @@ class BlockState(object):
assert self._check_clabel(), "invalid clabel before sweep"
Si = self.entropy(**entropy_args)
nmoves = -(mcmc_state.maccept.a * arange(len(mcmc_state.maccept.a))).sum()
try:
dS, rnmoves = self._multiflip_mcmc_sweep_dispatch(mcmc_state)
dS, nattempts, nmoves = self._multiflip_mcmc_sweep_dispatch(mcmc_state)
finally:
self.B = max(int(self.b.fa.max()) + 1, self.B)
nmoves += (mcmc_state.maccept.a * arange(len(mcmc_state.maccept.a))).sum()
if "mproposals" in kwargs:
kwargs["mproposals"][:M] = mcmc_state.mproposals.a
del kwargs["mproposals"]
......@@ -1622,7 +1622,7 @@ class BlockState(object):
if not dispatch:
return mcmc_state
return dS, nmoves
return dS, nattempts, nmoves
def _gibbs_sweep_dispatch(self, gibbs_state):
return libinference.gibbs_sweep(gibbs_state, self._state,
......@@ -1674,10 +1674,10 @@ class BlockState(object):
-------
dS : ``float``
Entropy difference after the sweeps.
nattempts : ``int``
Number of vertex moves attempted.
nmoves : ``int``
Number of vertices moved.
nattempts : ``int``
Number of attempted moves.
Notes
-----
......@@ -1707,7 +1707,7 @@ class BlockState(object):
Si = self.entropy(**entropy_args)
try:
dS, nmoves, nattempts = self._gibbs_sweep_dispatch(gibbs_state)
dS, nattempts, nmoves = self._gibbs_sweep_dispatch(gibbs_state)
finally:
self.B = max(int(self.b.fa.max()) + 1, self.B)
......@@ -1722,7 +1722,7 @@ class BlockState(object):
raise ValueError("unrecognized keyword arguments: " +
str(list(kwargs.keys())))
return dS, nmoves, nattempts
return dS, nattempts, nmoves
def _multicanonical_sweep_dispatch(self, multicanonical_state):
if multicanonical_state.multiflip:
......@@ -1753,6 +1753,8 @@ class BlockState(object):
-------
dS : ``float``
Entropy difference after the sweeps.
nattempts : ``int``
Number of vertex moves attempted.
nmoves : ``int``
Number of vertices moved.
......@@ -1801,7 +1803,7 @@ class BlockState(object):
multi_state.S_max))
try:
S, nmoves = self._multicanonical_sweep_dispatch(multi_state)
S, nattempts, nmoves = self._multicanonical_sweep_dispatch(multi_state)
finally:
self.B = max(int(self.b.fa.max()) + 1, self.B)
......@@ -1812,7 +1814,7 @@ class BlockState(object):
"inconsistent entropy after sweep %g (%g): %s" % \
(S, Sf, str(entropy_args))
return S, nmoves
return S, nattempts, nmoves
def _exhaustive_sweep_dispatch(self, exhaustive_state, callback, hist):
if callback is not None:
......@@ -1961,6 +1963,8 @@ class BlockState(object):
-------
dS : ``float``
Entropy difference after the sweeps.
nattempts : ``int``
Number of attempted merges.
nmoves : ``int``
Number of vertices merged.
......@@ -1982,7 +1986,7 @@ class BlockState(object):
assert self._check_clabel(), "invalid clabel before sweep"
Si = self.entropy(**entropy_args)
dS, nmoves = self._merge_sweep_dispatch(merge_state)
dS, nattempts, nmoves = self._merge_sweep_dispatch(merge_state)
if _bm_test():
assert self._check_clabel(), "invalid clabel after sweep"
......@@ -1991,7 +1995,7 @@ class BlockState(object):
"inconsistent entropy delta %g (%g): %s" % (dS, Sf - Si,
str(entropy_args))
return dS, nmoves
return dS, nattempts, nmoves