Commit 324b5495 authored by Tiago Peixoto's avatar Tiago Peixoto
Browse files

inference: implement normalized cut

parent 7822a48f
......@@ -372,6 +372,11 @@ libgraph_tool_inference_la_SOURCES = \
src/graph/inference/modularity/graph_modularity_mcmc.cc \
src/graph/inference/modularity/graph_modularity_multiflip_mcmc.cc \
src/graph/inference/modularity/graph_modularity_multilevel_mcmc.cc \
src/graph/inference/norm_cut/graph_norm_cut.cc \
src/graph/inference/norm_cut/graph_norm_cut_gibbs.cc \
src/graph/inference/norm_cut/graph_norm_cut_mcmc.cc \
src/graph/inference/norm_cut/graph_norm_cut_multiflip_mcmc.cc \
src/graph/inference/norm_cut/graph_norm_cut_multilevel_mcmc.cc \
src/graph/inference/partition_centroid/graph_partition_centroid.cc \
src/graph/inference/partition_centroid/graph_partition_centroid_mcmc.cc \
src/graph/inference/partition_centroid/graph_partition_centroid_multiflip_mcmc.cc \
......@@ -477,6 +482,10 @@ graph_inference_modularitydir = $(MOD_DIR)/include/inference/modularity
dist_graph_inference_modularity_HEADERS = \
src/graph/inference/modularity/graph_modularity.hh
graph_inference_norm_cutdir = $(MOD_DIR)/include/inference/norm_cut
dist_graph_inference_norm_cut_HEADERS = \
src/graph/inference/norm_cut/graph_norm_cut.hh
graph_inference_partition_centroiddir = $(MOD_DIR)/include/inference/partition_centroid
dist_graph_inference_partition_centroid_HEADERS = \
src/graph/inference/partition_centroid/graph_partition_centroid.hh \
......@@ -773,6 +782,7 @@ graph_tool_inference_PYTHON = \
src/graph_tool/inference/mcmc.py \
src/graph_tool/inference/minimize.py \
src/graph_tool/inference/modularity.py \
src/graph_tool/inference/norm_cut.py \
src/graph_tool/inference/partition_centroid.py \
src/graph_tool/inference/partition_modes.py \
src/graph_tool/inference/planted_partition.py \
......
......@@ -149,6 +149,11 @@ extern void export_modularity_gibbs();
extern void export_modularity_mcmc();
extern void export_modularity_multiflip_mcmc();
extern void export_modularity_multilevel_mcmc();
extern void export_norm_cut_state();
extern void export_norm_cut_gibbs();
extern void export_norm_cut_mcmc();
extern void export_norm_cut_multiflip_mcmc();
extern void export_norm_cut_multilevel_mcmc();
extern void export_latent_closure_state();
extern void export_latent_closure_mcmc();
extern void export_hist_state();
......@@ -239,6 +244,11 @@ BOOST_PYTHON_MODULE(libgraph_tool_inference)
export_modularity_mcmc();
export_modularity_multiflip_mcmc();
export_modularity_multilevel_mcmc();
export_norm_cut_state();
export_norm_cut_gibbs();
export_norm_cut_mcmc();
export_norm_cut_multiflip_mcmc();
export_norm_cut_multilevel_mcmc();
export_latent_closure_state();
export_latent_closure_mcmc();
export_hist_state();
......
// graph-tool -- a general graph modification and manipulation thingy
//
// Copyright (C) 2006-2022 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 Lesser 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 Lesser General Public License for more
// details.
//
// You should have received a copy of the GNU Lesser General Public License
// along with this program. If not, see <http://www.gnu.org/licenses/>.
#include "graph_tool.hh"
#include <boost/python.hpp>
#include "graph_norm_cut.hh"
#include "../support/graph_state.hh"
using namespace boost;
using namespace graph_tool;
GEN_DISPATCH(block_state, NormCutState, BLOCK_STATE_params)
python::object make_norm_cut_state(boost::python::object ostate)
{
python::object state;
block_state::make_dispatch(ostate,
[&](auto& s){state = python::object(s);});
return state;
}
void export_norm_cut_state()
{
using namespace boost::python;
def("make_norm_cut_state", &make_norm_cut_state);
block_state::dispatch
([&](auto* s)
{
typedef typename std::remove_reference<decltype(*s)>::type state_t;
void (state_t::*move_vertex)(size_t, size_t) =
&state_t::move_vertex;
double (state_t::*virtual_move)(size_t, size_t, size_t,
const norm_cut_entropy_args_t& ea) =
&state_t::virtual_move;
class_<state_t, bases<>, std::shared_ptr<state_t>>
c(name_demangle(typeid(state_t).name()).c_str(),
no_init);
c.def("move_vertex", move_vertex)
.def("virtual_move", virtual_move)
.def("entropy", &state_t::entropy);
});
class_<norm_cut_entropy_args_t>("norm_cut_entropy_args");
}
// graph-tool -- a general graph modification and manipulation thingy
//
// Copyright (C) 2006-2022 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 Lesser 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 Lesser General Public License for more
// details.
//
// You should have received a copy of the GNU Lesser General Public License
// along with this program. If not, see <http://www.gnu.org/licenses/>.
#ifndef GRAPH_NORM_CUT_HH
#define GRAPH_NORM_CUT_HH
#include "config.h"
#include <vector>
#include "idx_map.hh"
#include "../blockmodel/graph_blockmodel_util.hh"
#include "../support/graph_state.hh"
namespace graph_tool
{
using namespace boost;
using namespace std;
struct norm_cut_entropy_args_t
{
};
typedef vprop_map_t<int32_t>::type vmap_t;
#define BLOCK_STATE_params \
((g, &, never_directed, 1)) \
((_abg, &, boost::any&, 0)) \
((b,, vmap_t, 0)) \
((er, &, vector<size_t>&, 0)) \
((err, &, vector<size_t>&, 0))
GEN_STATE_BASE(NormCutStateBase, BLOCK_STATE_params)
template <class... Ts>
class NormCutState
: public NormCutStateBase<Ts...>
{
public:
GET_PARAMS_USING(NormCutStateBase<Ts...>, BLOCK_STATE_params)
GET_PARAMS_TYPEDEF(Ts, BLOCK_STATE_params)
typedef partition_stats<false> partition_stats_t;
template <class... ATs,
typename std::enable_if_t<sizeof...(ATs) == sizeof...(Ts)>* = nullptr>
NormCutState(ATs&&... args)
: NormCutStateBase<Ts...>(std::forward<ATs>(args)...),
_bg(boost::any_cast<std::reference_wrapper<bg_t>>(__abg)),
_N(HardNumVertices()(_g)),
_E(HardNumEdges()(_g)),
_bclabel(_N),
_pclabel(_N),
_wr(_N)
{
_wr.resize(num_vertices(_g), 0);
_er.resize(num_vertices(_g), 0);
_err.resize(num_vertices(_g), 0);
for (auto v : vertices_range(_g))
{
auto r = _b[v];
_er[r] += out_degree(v, _g);
_wr[r]++;
}
for (size_t r = 0; r < _N; ++r)
{
if (_wr[r] == 0)
_empty_groups.insert(r);
else
_candidate_groups.insert(r);
}
for (auto e : edges_range(_g))
{
auto r = _b[source(e, _g)];
auto s = _b[target(e, _g)];
if (r == s)
_err[r] += 2;
}
}
typedef typename
std::conditional<is_directed_::apply<g_t>::type::value,
GraphInterface::multigraph_t,
undirected_adaptor<GraphInterface::multigraph_t>>::type
bg_t;
bg_t& _bg;
size_t _N;
size_t _E;
idx_set<size_t> _empty_groups;
idx_set<size_t> _candidate_groups;
std::vector<size_t> _bclabel;
std::vector<size_t> _pclabel;
std::vector<size_t> _wr;
constexpr static BlockStateVirtualBase* _coupled_state = nullptr;
typedef int m_entries_t;
UnityPropertyMap<int,GraphInterface::vertex_t> _vweight;
UnityPropertyMap<int,GraphInterface::edge_t> _eweight;
simple_degs_t _degs;
typedef norm_cut_entropy_args_t _entropy_args_t;
// =========================================================================
// State modification
// =========================================================================
void move_vertex(size_t v, size_t nr)
{
size_t r = _b[v];
if (nr == r)
return;
size_t k = 0;
size_t m = 0;
for (auto e : out_edges_range(v, _g))
{
++k;
auto u = target(e, _g);
if (u == v)
{
++m;
continue;
}
size_t s = _b[u];
if (s == r)
_err[r] -= 2;
else if (s == nr)
_err[nr] += 2;
}
_err[r] -= m;
_err[nr] += m;
_er[r] -= k;
_er[nr] += k;
_wr[r]--;
_wr[nr]++;
if (_wr[r] == 0)
{
_empty_groups.insert(r);
_candidate_groups.erase(r);
}
if (_wr[nr] == 1)
{
_empty_groups.erase(nr);
_candidate_groups.insert(nr);
}
_b[v] = nr;
}
size_t virtual_remove_size(size_t v)
{
return _wr[_b[v]] - 1;
}
constexpr void add_block(size_t)
{
}
double virtual_move(size_t v, size_t r, size_t nr,
const norm_cut_entropy_args_t&)
{
if (r == nr)
return 0;
std::array<int, 2> derr({0,0});
size_t k = 0;
size_t m = 0;
for (auto e : out_edges_range(v, _g))
{
++k;
auto u = target(e, _g);
if (u == v)
{
++m;
continue;
}
size_t s = _b[u];
if (s == r)
derr[0] -= 2;
else if (s == nr)
derr[1] += 2;
}
derr[0] -= m;
derr[1] += m;
double Cb = 0;
double Ca = 0;
Cb -= (_er[r] > 0) ? _err[r]/double(_er[r]) : 0;
Cb -= (_er[nr] > 0) ? _err[nr]/double(_er[nr]) : 0;
Ca -= (_er[r] - k > 0) ? (_err[r] + derr[0])/double(_er[r] - k) : 0;
Ca -= (_er[nr] + k > 0) ? (_err[nr] + derr[1])/double(_er[nr] + k) : 0;
int dB = 0;
if (_wr[r] == 1)
dB--;
if (_wr[nr] == 0)
dB++;
Cb += _candidate_groups.size();
Ca += _candidate_groups.size() + dB;
double dS = (Ca - Cb);
return dS;
}
size_t get_empty_block(size_t, bool=false)
{
return *(_empty_groups.end() - 1);
}
size_t sample_block(size_t v, double c, double d, rng_t& rng)
{
std::bernoulli_distribution new_r(d);
if (d > 0 && !_empty_groups.empty() && new_r(rng))
return uniform_sample(_empty_groups, rng);
c = std::max(std::min(c, 1.), 0.);
std::bernoulli_distribution adj(1.-c);
auto iter = out_neighbors(v, _g);
if (iter.first != iter.second && adj(rng))
{
auto w = uniform_sample(iter.first, iter.second, rng);
return _b[w];
}
return uniform_sample(_candidate_groups, rng);
}
size_t sample_block_local(size_t v, rng_t& rng)
{
return sample_block(v, 0, 0, rng);
}
// Computes the move proposal probability
double get_move_prob(size_t v, size_t r, size_t s, double c, double d,
bool reverse)
{
size_t B = _candidate_groups.size();
if (reverse)
{
if (_wr[s] == 1)
return log(d);
if (_wr[r] == 0)
B++;
}
else
{
if (_wr[s] == 0)
return log(d);
}
size_t k_s = 0;
size_t k = 0;
for (auto w : out_neighbors_range(v, _g))
{
if (size_t(_b[w]) == s)
k_s++;
k++;
}
if (B == _N)
d = 0;
if (k > 0)
{
double p = k_s / double(k);
c = 1 - std::max(std::min(c, 1.), 0.);
return log1p(-d) + log(c * p + (1. - c)/B);
}
else
{
return log1p(-d) - log(B);
}
}
template <class MEntries>
double get_move_prob(size_t v, size_t r, size_t s, double c, double d,
bool reverse, MEntries&&)
{
return get_move_prob(v, r, s, c, d, reverse);
}
template <class EArgs, class MEntries>
double virtual_move(size_t v, size_t r, size_t nr, EArgs&& ea, MEntries&&)
{
return virtual_move(v, r, nr, ea);
}
double entropy(const norm_cut_entropy_args_t&)
{
size_t B = _candidate_groups.size();
double C = B;
for (auto r : _candidate_groups)
C -= (_er[r] > 0) ? _err[r] / double(_er[r]) : 0;
return C;
}
template <class MCMCState>
void init_mcmc(MCMCState&)
{
}
constexpr size_t node_weight(size_t)
{
return 1;
}
bool is_last(size_t v)
{
return _wr[_b[v]] == 1;
}
constexpr bool allow_move(size_t, size_t)
{
return true;
}
template <class V>
void push_state(V&) {}
void pop_state() {}
void store_next_state(size_t) {}
void clear_next_state() {}
void relax_update(bool) {}
};
} // graph_tool namespace
#endif //GRAPH_NORM_CUT_HH
// graph-tool -- a general graph modification and manipulation thingy
//
// Copyright (C) 2006-2022 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 Lesser 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 Lesser General Public License for more
// details.
//
// You should have received a copy of the GNU Lesser General Public License
// along with this program. If not, see <http://www.gnu.org/licenses/>.
#include "graph_tool.hh"
#include "random.hh"
#include <boost/python.hpp>
#include "graph_norm_cut.hh"
#include "../blockmodel/graph_blockmodel_gibbs.hh"
#include "../loops/gibbs_loop.hh"
using namespace boost;
using namespace graph_tool;
GEN_DISPATCH(block_state, NormCutState, BLOCK_STATE_params)
template <class State>
GEN_DISPATCH(gibbs_block_state, Gibbs<State>::template GibbsBlockState,
GIBBS_BLOCK_STATE_params(State))
python::object norm_cut_gibbs_sweep(python::object ogibbs_state,
python::object oblock_state,
rng_t& rng)
{
python::object ret;
auto dispatch = [&](auto& block_state)
{
typedef typename std::remove_reference<decltype(block_state)>::type
state_t;
gibbs_block_state<state_t>::make_dispatch
(ogibbs_state,
[&](auto& s)
{
auto ret_ = gibbs_sweep(*s, rng);
ret = tuple_apply([&](auto&... args){ return python::make_tuple(args...); }, ret_);
});
};
block_state::dispatch(oblock_state, dispatch);
return ret;
}
void export_norm_cut_gibbs()
{
using namespace boost::python;
def("norm_cut_gibbs_sweep", &norm_cut_gibbs_sweep);
}
// graph-tool -- a general graph modification and manipulation thingy
//
// Copyright (C) 2006-2022 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 Lesser 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 Lesser General Public License for more
// details.
//
// You should have received a copy of the GNU Lesser General Public License
// along with this program. If not, see <http://www.gnu.org/licenses/>.
#include "graph_tool.hh"
#include "random.hh"
#include <boost/python.hpp>
#include "graph_norm_cut.hh"
#include "../blockmodel/graph_blockmodel_mcmc.hh"
#include "../loops/mcmc_loop.hh"
using namespace boost;
using namespace graph_tool;
GEN_DISPATCH(block_state, NormCutState, BLOCK_STATE_params)
template <class State>
GEN_DISPATCH(mcmc_block_state, MCMC<State>::template MCMCBlockState,
MCMC_BLOCK_STATE_params(State))
python::object norm_cut_mcmc_sweep(python::object omcmc_state,
python::object oblock_state,
rng_t& rng)
{
python::object ret;
auto dispatch = [&](auto& block_state)
{
typedef typename std::remove_reference<decltype(block_state)>::type
state_t;
mcmc_block_state<state_t>::make_dispatch
(omcmc_state,
[&](auto& s)
{
auto ret_ = mcmc_sweep(*s, rng);
ret = tuple_apply([&](auto&... args){ return python::make_tuple(args...); }, ret_);
});
};
block_state::dispatch(oblock_state, dispatch);
return ret;
}
void export_norm_cut_mcmc()
{
using namespace boost::python;
def("norm_cut_mcmc_sweep", &norm_cut_mcmc_sweep);
}