From 6e6c4f1a7cc822b787bc4a742a3cd56238479f43 Mon Sep 17 00:00:00 2001
From: Kilian Verhetsel <kilian.verhetsel@student.uclouvain.be>
Date: Tue, 29 Nov 2016 08:36:35 +0000
Subject: [PATCH] Added a solver for the Maximum Weight Independent Set problem
---
Mesh/mwis.hpp | 1037 +++++++++++++++++++++++++++++++++++++++++++++++++
1 file changed, 1037 insertions(+)
create mode 100644 Mesh/mwis.hpp
diff --git a/Mesh/mwis.hpp b/Mesh/mwis.hpp
new file mode 100644
index 0000000000..16277dbf89
--- /dev/null
+++ b/Mesh/mwis.hpp
@@ -0,0 +1,1037 @@
+#include "search.hpp"
+
+#include <boost/iterator/counting_iterator.hpp>
+
+#include <boost/graph/random.hpp>
+#include <boost/graph/graph_traits.hpp>
+#include <boost/graph/properties.hpp>
+#include <boost/graph/graph_concepts.hpp>
+#include <boost/graph/adjacency_list.hpp>
+
+#include <random>
+#include <vector>
+#include <functional>
+#include <limits>
+#include <numeric>
+#include <algorithm>
+#include <iostream>
+
+namespace mwis {
+
+/**
+ * A functor to calculate an upper bound to the maximum weight independent set
+ * in a graph.
+ *
+ * A Lagrangian relaxation of the problem is used:
+ *
+ * UB = min_{λ >= 0} max_{x_i \in {0, 1}}
+ * (\sum_{i \in V} x_iw_i) -
+ * \sum_{c \in C} λ_c(\sum_{i \in c} - 1)
+ *
+ * Where C is a set of cliques within the graph, such that if there is an edge
+ * between two nodes in the graph, a clique that contains both nodes is present
+ * in C.
+ *
+ * The set {{u, v} | (u, v) \in E} can always be used as the set of
+ * cliques. However, using larger cliques can lead to better upper bounds.
+ */
+template<typename Graph, typename WeightMap>
+class lagrangian_bound {
+ typedef typename boost::graph_traits<Graph>::vertex_descriptor vertex;
+ typedef typename boost::property_traits<WeightMap>::value_type weight;
+
+ typedef typename boost::property_map<
+ Graph, boost::vertex_index_t>::type index_map;
+
+ typedef boost::iterator_property_map<
+ std::vector<bool>::iterator, index_map> selected_map_type;
+ typedef boost::iterator_property_map<
+ std::vector<std::vector<size_t>>::iterator, index_map> clique_map_type;
+ typedef boost::iterator_property_map<
+ typename std::vector<weight>::iterator, index_map> weight_map_type;
+
+ const Graph &_graph;
+ WeightMap _weight;
+ size_t _max_size;
+
+ std::vector<std::vector<vertex>> _clique_contents;
+
+ std::vector<weight> _lambda, _gradient;
+
+ std::vector<std::vector<size_t>> _clique_storage;
+ std::vector<bool> _selected_storage;
+ std::vector<weight> _effective_weight_storage;
+
+ clique_map_type _cliques;
+ selected_map_type _selected;
+ weight_map_type _effective_weight;
+public:
+ typedef weight result_type;
+
+ template<typename CliqueIterator>
+ lagrangian_bound(const Graph &graph, WeightMap weight_map,
+ CliqueIterator clique_begin,
+ CliqueIterator clique_end,
+ size_t max_size = std::numeric_limits<size_t>::max()):
+ _graph(graph), _weight(weight_map), _max_size(max_size) {
+ auto vs = vertices(graph);
+ size_t vertex_count = vs.second - vs.first;
+
+ _clique_storage.resize(vertex_count);
+ _selected_storage.resize(vertex_count);
+ _effective_weight_storage.resize(vertex_count);
+
+ index_map ids = boost::get(boost::vertex_index, graph);
+
+ _cliques = make_iterator_property_map(_clique_storage.begin(), ids);
+ _selected = make_iterator_property_map(_selected_storage.begin(), ids);
+ _effective_weight = make_iterator_property_map(
+ _effective_weight_storage.begin(), ids);
+
+ size_t i = 0;
+ for (CliqueIterator it = clique_begin; it != clique_end; it++, i++) {
+ std::vector<vertex> to_add;
+ for (const vertex &v : *it) {
+ to_add.push_back(v);
+ get(_cliques, v).push_back(i);
+ }
+ _clique_contents.push_back(to_add);
+ }
+
+ _lambda.resize(_clique_contents.size());
+ _gradient.resize(_clique_contents.size());
+
+ for (weight &v : _lambda) v = weight(0);
+ }
+
+ /**
+ * Calculates an upper bound. This bound is always admissible, but may not be
+ * optimal.
+ *
+ * Vertices in [solution_begin, solution_end) are all forced to be part of the
+ * solution.
+ *
+ * Vertices in [rest_begin, rest_end) can be selected in any way to maximize
+ * the objective functions.
+ *
+ * Vertices in neither range are assumed not to be part of the solution.
+ */
+ template<typename SolutionIterator,
+ typename RestIterator>
+ weight operator()(SolutionIterator solution_begin,
+ SolutionIterator solution_end,
+ RestIterator rest_begin,
+ RestIterator rest_end) {
+ weight best(std::numeric_limits<weight>::max());
+
+ for (size_t i = 0; i < _selected_storage.size(); i++)
+ _selected_storage[i] = false;
+
+ for (SolutionIterator it = solution_begin; it != solution_end; it++)
+ put(_selected, *it, true);
+
+ /**
+ * We can get away with very few steps, because the starting point is the
+ * optimum that was found the last time the algorithm was used, and the
+ * input does not usually change much between two calls (therefore, one
+ * would expect the optima to be close to each other).
+ */
+ static const size_t n = 2;
+ for (size_t i = 0; i < n; i++) {
+ compute_weights();
+ if (_max_size == std::numeric_limits<size_t>::max())
+ select(rest_begin, rest_end);
+ else
+ select(rest_begin, rest_end,
+ _max_size - (size_t)std::distance(solution_begin,
+ solution_end));
+ best = std::min(best, evaluate());
+
+ if (i != n - 1) {
+ compute_gradient();
+ follow_gradient(weight(-1e-4 / 1));
+ }
+ }
+
+ return best;
+ }
+
+ /**
+ * For debugging purposes, displays any inconsistencies between the graph and
+ * the sequence of cliques.
+ *
+ * (u, v) \in E <=> there is a clique which contains both u and v.
+ */
+ void check_consistency() {
+ auto es = boost::edges(_graph);
+ for (auto eit = es.first; eit != es.second; eit++) {
+ bool found = false;
+
+ vertex v = boost::source(*eit, _graph);
+ vertex w = boost::target(*eit, _graph);
+
+ const std::vector<size_t> &cliques_v = get(_cliques, v);
+ const std::vector<size_t> &cliques_w = get(_cliques, w);
+
+ for (size_t clique : cliques_v) {
+ if (std::find(cliques_w.begin(), cliques_w.end(), clique) !=
+ cliques_w.end()) {
+ found = true;
+ break;
+ }
+ }
+
+ if (!found) {
+ std::cout << "no shared clique between connected vertices: " <<
+ v << " and " << w << "\n";
+ }
+ }
+
+ for (const std::vector<vertex> &clique : _clique_contents) {
+ for (vertex v : clique) {
+ for (vertex w : clique) {
+ if (v == w) continue;
+
+ if (!edge(v, w, _graph).second) {
+ std::cout << "no edge between " << v << " and " <<
+ w << "\n";
+ }
+ }
+ }
+ }
+ }
+private:
+ /**
+ * Calculates the gradient according to which nodes are currently selected.
+ */
+ void compute_gradient() {
+ std::transform(_clique_contents.begin(), _clique_contents.end(),
+ _gradient.begin(),
+ [&](const std::vector<vertex> &vertices) {
+ weight sum{1};
+ for (const vertex &v : vertices) {
+ sum -= weight(get(_selected, v) ? 1 : 0);
+ }
+
+ return sum;
+ });
+ }
+
+ /**
+ * Calculates the effective weight of each node (w_i - \sum_{c, x \in c} λ_c).
+ */
+ void compute_weights() {
+ for (size_t i = 0; i < _effective_weight_storage.size(); i++) {
+ weight result(_weight[i]);
+
+ for (size_t clique : _clique_storage[i])
+ result -= _lambda[clique];
+
+ _effective_weight_storage[i] = result;
+ }
+ }
+
+ /**
+ * Update the Lagrange multipliers in accordance with the last computed
+ * gradient, multiplied by scale.
+ */
+ void follow_gradient(weight scale) {
+ std::transform(_lambda.begin(), _lambda.end(), _gradient.begin(),
+ _lambda.begin(), [&](const weight &cur, const weight &grad) {
+ return std::max(weight(0), cur + scale * grad);
+ });
+ }
+
+ /**
+ * Iterates over the given range, and selects the subset of its nodes which have
+ * a non-negative weight, so as to maximize the objective function.
+ */
+ template<typename Iterator>
+ void select(Iterator begin, Iterator end) {
+ for (Iterator it = begin; it != end; it++) {
+ vertex v = *it;
+ put(_selected, v, get(_effective_weight, v) >= weight(0));
+ }
+ }
+
+ /**
+ * Iterates over the given range, and selects subset of its nodes which have
+ * a non-negative weight, selecting at most n nodes.
+ */
+ template<typename Iterator>
+ void select(Iterator begin, Iterator end, size_t n) {
+ n = std::min((size_t)std::distance(begin, end), n);
+ std::vector<vertex> largest(n);
+
+ std::partial_sort_copy(
+ begin, end,
+ largest.begin(), largest.end(),
+ [&](vertex a ,vertex b) {
+ return get(_effective_weight, a) > get(_effective_weight, b);
+ });
+
+ for (Iterator it = begin; it != end; it++) {
+ vertex v = *it;
+ put(_selected, v, false);
+ }
+
+ for (vertex v : largest) {
+ if (get(_effective_weight, v) <= weight(0))
+ break;
+ put(_selected, v, true);
+ }
+ }
+
+
+ /**
+ * Calculates the upper bound for the current values of the Lagrange
+ * multipliers and the current selection of nodes.
+ */
+ weight evaluate() {
+ weight result(std::accumulate(_lambda.begin(), _lambda.end(),
+ weight(0), std::plus<weight>()));
+
+ for (size_t i = 0; i < _selected_storage.size(); i++) {
+ if (_selected_storage[i])
+ result += _effective_weight_storage[i];
+ }
+
+ return result;
+ }
+};
+
+template<typename T>
+struct max_bound {
+ typedef T result_type;
+
+ template<typename SolutionIter, typename RestIter>
+ T operator()(SolutionIter solution_begin, SolutionIter solution_end,
+ RestIter rest_begin, RestIter rest_end) {
+ return std::numeric_limits<T>::max();
+ }
+};
+
+template<typename Graph, typename WeightMap>
+class state_change;
+
+template<typename Graph, typename WeightMap>
+class clique_assignment;
+
+template<typename Graph, typename WeightMap>
+struct state {
+ typedef typename boost::graph_traits<Graph>::vertex_descriptor vertex;
+ typedef typename boost::property_traits<WeightMap>::value_type weight;
+
+ typedef typename boost::property_map<
+ Graph, boost::vertex_index_t>::type index_map;
+ typedef boost::iterator_property_map<
+ std::vector<std::vector<size_t>>::iterator, index_map> clique_map_type;
+
+ const Graph &graph;
+ WeightMap weight_map;
+
+ std::vector<std::vector<vertex>> cliques;
+ std::vector<size_t> clique_sizes;
+ std::set<vertex> invalid;
+ std::set<vertex> selectable;
+ std::vector<bool> assigned;
+
+ std::vector<std::vector<size_t>> clique_map_storage;
+ clique_map_type clique_map;
+
+ std::vector<vertex> solution;
+ weight solution_value;
+
+ template<typename CliqueIterator>
+ state(const Graph &graph,
+ WeightMap weight_map,
+ CliqueIterator clique_begin,
+ CliqueIterator clique_end):
+ graph(graph), weight_map(weight_map),
+ cliques(clique_begin, clique_end),
+ clique_sizes(cliques.size()),
+ selectable(vertices(graph).first, vertices(graph).second),
+ assigned(cliques.size(), false),
+ solution_value(0)
+ {
+ clique_map_storage.resize(num_vertices(graph));
+ index_map id_map = get(boost::vertex_index, graph);
+ clique_map = make_iterator_property_map(
+ clique_map_storage.begin(), id_map);
+
+ for (size_t i = 0; i < cliques.size(); i++)
+ clique_sizes[i] = cliques[i].size();
+
+ size_t id = 0;
+ for (std::vector<vertex> &clique : cliques) {
+ for (vertex v : clique)
+ get(clique_map, v).push_back(id);
+
+ std::sort(clique.begin(), clique.end(), [&](vertex a, vertex b) {
+ return get(weight_map, a) > get(weight_map, b);
+ });
+
+ id++;
+ }
+ }
+};
+
+template<typename Graph, typename WeightMap>
+class clique_assignment {
+ typedef typename boost::graph_traits<Graph>::vertex_descriptor vertex;
+ typedef typename boost::property_traits<WeightMap>::value_type weight;
+
+ size_t _id;
+ vertex _vertex;
+ bool _selected;
+
+public:
+ typedef state_change<Graph, WeightMap> result_type;
+ typedef state<Graph, WeightMap> argument_type;
+
+ clique_assignment(size_t id):
+ _id(id), _selected(false)
+ {}
+
+ clique_assignment(size_t id, vertex vertex_id):
+ _id(id), _vertex(vertex_id), _selected(true)
+ {}
+
+ state_change<Graph, WeightMap> operator()(
+ const state<Graph, WeightMap> &state) const {
+ std::vector<vertex> selected_vertex, removed;
+ std::vector<size_t> assigned;
+
+ weight new_weight(state.solution_value);
+
+ if (_selected) {
+ selected_vertex.push_back(_vertex);
+ new_weight += get(state.weight_map, _vertex);
+
+ auto edges = out_edges(_vertex, state.graph);
+ for (auto eit = edges.first; eit != edges.second; eit++) {
+ const vertex &other = target(*eit, state.graph);
+ if (state.invalid.find(other) == state.invalid.end()) {
+ removed.push_back(other);
+ }
+ }
+
+ for (size_t i : get(state.clique_map, _vertex)) {
+ if (!state.assigned[i]) assigned.push_back(i);
+ }
+ }
+ else {
+ for (vertex v : state.cliques[_id]) {
+ if (state.invalid.find(v) == state.invalid.end())
+ removed.push_back(v);
+ }
+
+ assigned.push_back(_id);
+ }
+
+ std::map<size_t, size_t> size_delta;
+ for (vertex v : removed) {
+ for (size_t clique : get(state.clique_map, v)) {
+ size_delta[clique]++;
+ if (size_delta[clique] == state.clique_sizes[clique]) {
+ if (!state.assigned[clique])
+ assigned.push_back(clique);
+ }
+ }
+ }
+
+ return state_change<Graph, WeightMap>(
+ assigned,
+ size_delta,
+ selected_vertex, removed,
+ state.solution_value, new_weight);
+ }
+};
+
+template<typename Graph, typename WeightMap>
+class state_change {
+ typedef typename boost::graph_traits<Graph>::vertex_descriptor vertex;
+ typedef typename boost::property_traits<WeightMap>::value_type weight;
+
+ std::vector<size_t> _assigned_id;
+ std::map<size_t, size_t> _size_delta;
+ std::vector<vertex> _selected;
+ std::vector<vertex> _invalidated;
+ weight _old_weight, _new_weight;
+public:
+ state_change() {}
+
+ state_change(const std::vector<size_t> &assigned_id,
+ const std::map<size_t, size_t> &size_delta,
+ const std::vector<vertex> &selected,
+ const std::vector<vertex> &invalidated,
+ weight old_weight, weight new_weight):
+ _assigned_id(assigned_id),
+ _size_delta(size_delta),
+ _selected(selected), _invalidated(invalidated),
+ _old_weight(old_weight), _new_weight(new_weight)
+ {}
+
+ void apply(state<Graph, WeightMap> &state) const {
+ for (size_t i : _assigned_id)
+ state.assigned[i] = true;
+
+ size_t old_size = state.solution.size();
+ state.solution.resize(state.solution.size() + _selected.size());
+ std::copy(_selected.begin(), _selected.end(),
+ state.solution.begin() + old_size);
+
+ for (vertex v : _selected)
+ state.selectable.erase(v);
+
+ state.solution_value = _new_weight;
+
+ for (vertex v : _invalidated) {
+ state.invalid.insert(v);
+ state.selectable.erase(v);
+ }
+
+ for (auto pair : _size_delta)
+ state.clique_sizes[pair.first] -= pair.second;
+ }
+
+ void reverse(state<Graph, WeightMap> &state) const {
+ for (size_t i : _assigned_id)
+ state.assigned[i] = false;
+
+ state.solution.resize(state.solution.size() - _selected.size());
+
+ for (vertex v : _selected)
+ state.selectable.insert(v);
+
+ state.solution_value = _old_weight;
+
+ for (vertex v : _invalidated) {
+ state.invalid.erase(v);
+ state.selectable.insert(v);
+ }
+
+ for (auto pair : _size_delta)
+ state.clique_sizes[pair.first] += pair.second;
+ }
+};
+
+template<typename Graph, typename WeightMap>
+struct visit_state {
+ typedef typename boost::graph_traits<Graph>::vertex_descriptor vertex;
+ typedef typename boost::property_traits<WeightMap>::value_type weight;
+
+ bool quiet;
+
+ weight best_value;
+ std::vector<vertex> best_solution;
+
+ visit_state(bool quiet = true): quiet(quiet), best_value(0) {}
+
+ void operator()(state<Graph, WeightMap> &state) {
+ if (best_value < state.solution_value) {
+ if (!quiet) {
+ std::cout << "new best: " << state.solution_value <<
+ " (" << state.solution.size() << ")\n";
+ }
+
+ best_value = state.solution_value;
+ best_solution = state.solution;
+ }
+ }
+};
+
+template<typename Graph, typename WeightMap, typename UpperBound>
+class successor {
+ typedef typename boost::graph_traits<Graph>::vertex_descriptor vertex;
+ typedef typename std::vector<vertex>::iterator viterator;
+ typedef typename boost::property_traits<WeightMap>::value_type weight;
+
+ UpperBound &_bound;
+
+ size_t _visited;
+ weight *_best_weight;
+
+ size_t _limit;
+public:
+ typedef clique_assignment<Graph, WeightMap> action_type;
+
+ successor(UpperBound &bound, weight &best_weight,
+ size_t limit = std::numeric_limits<size_t>::max()):
+ _bound(bound), _visited(0), _best_weight(&best_weight), _limit(limit)
+ {}
+
+ template<typename OutputIterator>
+ void operator()(state<Graph, WeightMap> &state, OutputIterator it) {
+ if (_visited > 1000 * 1000 * 10)
+ return;
+ else
+ _visited++;
+
+ if (state.solution.size() == _limit)
+ return;
+
+ auto clique_begin = boost::make_counting_iterator<size_t>(0);
+ auto clique_end = boost::make_counting_iterator<size_t>(
+ state.cliques.size());
+
+ auto clique_it = std::min_element(
+ clique_begin, clique_end, [&](size_t a, size_t b) {
+ return clique_key(state, a) < clique_key(state, b);
+ });
+
+ if (clique_it == clique_end || state.assigned[*clique_it])
+ return;
+
+ std::vector<vertex> vertices;
+ size_t id = *clique_it;
+ for (vertex v : state.cliques[id]) {
+ if (state.invalid.find(v) == state.invalid.end()) {
+ vertices.push_back(v);
+ }
+ }
+
+ if (vertices.size() > 0) {
+ weight max(_bound(state.solution.begin(), state.solution.end(),
+ state.selectable.begin(), state.selectable.end()));
+ if (max <= *_best_weight)
+ return;
+ }
+
+ std::sort(vertices.begin(), vertices.end(),
+ [&](vertex a, vertex b) {
+ return vertex_regret(state, a) < vertex_regret(state, b);
+ });
+
+ for (vertex v : vertices)
+ *it++ = clique_assignment<Graph, WeightMap>(id, v);
+
+ *it++ = clique_assignment<Graph, WeightMap>(id);
+ }
+
+ std::tuple<int, size_t, weight>
+ clique_key(state<Graph, WeightMap> &state,
+ size_t id) {
+ if (state.assigned[id]) {
+ return std::make_tuple(1,
+ std::numeric_limits<size_t>::max(),
+ weight(0));
+ }
+ else if (state.clique_sizes[id] == 0) {
+ return std::make_tuple(0, 0, weight(0));
+ }
+ else {
+ vertex v = *std::find_if(
+ state.cliques[id].begin(), state.cliques[id].end(), [&](vertex v) {
+ return state.invalid.find(v) == state.invalid.end();
+ });
+
+ return std::make_tuple(id < 1361 ? 0 : 1,
+ state.clique_sizes[id],
+ -get(state.weight_map, v));
+ }
+ }
+
+ weight vertex_regret(const state<Graph, WeightMap> &state, vertex v) {
+ weight result(-get(state.weight_map, v));
+
+ auto es = out_edges(v, state.graph);
+ for (auto eit = es.first; eit != es.second; eit++) {
+ vertex u = target(*eit, state.graph);
+ if (state.invalid.find(u) != state.invalid.end())
+ result += get(state.weight_map, u);
+ }
+
+ return result;
+ }
+};
+
+template<typename Graph, typename WeightMap>
+class evaluator {
+ typedef typename boost::graph_traits<Graph>::vertex_descriptor vertex;
+ typedef typename boost::property_traits<WeightMap>::value_type weight;
+
+ std::vector<std::vector<vertex>> _cliques;
+ size_t _limit;
+public:
+ typedef weight result_type;
+ typedef state<Graph, WeightMap> argument_type;
+
+ template<typename Iterator>
+ evaluator(Iterator begin, Iterator end, size_t limit):
+ _cliques(begin, end), _limit(limit)
+ {}
+
+ weight operator()(state<Graph, WeightMap> &cur) const {
+ visit_state<Graph, WeightMap> visitor;
+
+ // max_bound<weight> bound;
+ lagrangian_bound<Graph, WeightMap>
+ bound(cur.graph, cur.weight_map, _cliques.begin(), _cliques.end(),
+ cur.solution.size() + _limit);
+ successor<Graph, WeightMap, decltype(bound)>
+ successor(bound, visitor.best_value, cur.solution.size() + _limit);
+
+ visitor(cur);
+ search::depth_limited_search(cur, visitor, successor, _limit);
+
+ return visitor.best_value;
+ }
+};
+
+template<typename Graph, typename WeightMap>
+struct lns_state {
+ typedef typename boost::graph_traits<Graph>::vertex_descriptor vertex;
+
+ const Graph &graph;
+ WeightMap weight_map;
+
+ std::vector<vertex> solution;
+
+ lns_state(const Graph &graph, WeightMap weight_map,
+ const std::vector<vertex> &solution):
+ graph(graph), weight_map(weight_map), solution(solution)
+ {}
+};
+
+template<typename Graph>
+struct lns_fragment {
+ typedef typename boost::graph_traits<Graph>::vertex_descriptor vertex;
+
+ std::set<vertex> vertices;
+
+ lns_fragment(const std::set<vertex> &vertices):
+ vertices(vertices)
+ {}
+};
+
+/**
+ * Behaves like a queue, but elements can only be inserted into the queue a
+ * limited amount of times.
+ */
+template<typename T>
+class limited_queue {
+ std::queue<T> _queue;
+ size_t _size, _max_size;
+public:
+ typedef T value_type;
+ typedef typename std::queue<T>::size_type size_type;
+
+ limited_queue(size_t max_size):
+ _size(0), _max_size(max_size) {}
+
+ void push(const T &t) {
+ if (_size != _max_size) {
+ _size++;
+ _queue.push(t);
+ }
+ }
+
+ void pop() { _queue.pop(); }
+
+ T &top() { return _queue.front(); }
+ const T &top() const { return _queue.front(); }
+
+ size_type size() const { return _queue.size(); }
+
+ bool empty() { return _queue.empty(); }
+};
+
+/**
+ * A visitor for graph traversal algoritms which records the target of every
+ * visited edge into a set.
+ */
+template<typename Graph, typename Tag>
+class set_recorder : public boost::base_visitor<set_recorder<Graph, Tag>> {
+ typedef typename boost::graph_traits<Graph>::vertex_descriptor vertex;
+ typedef typename boost::graph_traits<Graph>::edge_descriptor edge;
+
+ std::set<vertex> &_set;
+public:
+ typedef Tag event_filter;
+
+ set_recorder(std::set<vertex> &set): _set(set) {}
+
+ void operator()(edge e, const Graph &g) {
+ _set.insert(boost::target(e, g));
+ }
+};
+
+template<typename Graph, typename WeightMap>
+lns_fragment<Graph> fragment_selector(const lns_state<Graph, WeightMap> &state) {
+ typedef typename boost::graph_traits<Graph>::vertex_descriptor vertex;
+
+ std::random_device dev;
+ std::mt19937 gen(dev());
+
+ vertex v = boost::random_vertex(state.graph, gen);
+
+ size_t neighborhood_size = std::uniform_int_distribution<size_t>(50, 80)(gen);
+
+ std::set<vertex> set;
+ set.insert(v);
+ set_recorder<Graph, boost::on_black_target> recorder(set);
+
+ limited_queue<vertex> queue(neighborhood_size);
+
+ boost::breadth_first_search(
+ state.graph, v,
+ boost::buffer(queue).
+ visitor(boost::make_bfs_visitor(recorder)));
+
+ return lns_fragment<Graph>(set);
+}
+
+template<typename Graph>
+class lns_assignment {
+ typedef typename boost::graph_traits<Graph>::vertex_descriptor vertex;
+ std::vector<vertex> _vertices;
+public:
+ lns_assignment(const std::vector<vertex> &vertices): _vertices(vertices) {}
+
+ template<typename WeightMap>
+ void operator()(lns_state<Graph, WeightMap> &state,
+ const lns_fragment<Graph> &fragment) const {
+ static int i = 0;
+ std::cout << "\rsize: " << state.solution.size() << " " << i++;
+ std::cout.flush();
+
+ state.solution.erase(
+ std::remove_if(
+ state.solution.begin(), state.solution.end(),
+ [&](vertex v) {
+ return fragment.vertices.find(v) != fragment.vertices.end();
+ }),
+ state.solution.end());
+
+ std::copy(_vertices.begin(), _vertices.end(),
+ std::back_inserter(state.solution));
+ }
+};
+
+template<typename Graph, typename WeightMap>
+class lns_search {
+ typedef typename boost::graph_traits<Graph>::vertex_descriptor vertex;
+ typedef typename boost::property_traits<WeightMap>::value_type weight;
+
+ std::vector<std::vector<vertex>> _cliques;
+public:
+ typedef lns_assignment<Graph> result_type;
+ typedef lns_state<Graph, WeightMap> first_argument_type;
+ typedef lns_fragment<Graph> second_argument_type;
+
+ template<typename Iterator>
+ lns_search(Iterator begin, Iterator end): _cliques(begin, end) {}
+
+ lns_assignment<Graph> operator()(const lns_state<Graph, WeightMap> &l_state,
+ const lns_fragment<Graph> &fragment) const {
+ state<Graph, WeightMap> search_state(l_state.graph, l_state.weight_map,
+ _cliques.begin(), _cliques.end());
+
+ std::vector<vertex> selected_vertex;
+ std::set<vertex> removed;
+ std::set<size_t> assigned;
+
+ weight new_weight(0);
+
+ for (vertex v : l_state.solution) {
+ if (fragment.vertices.find(v) != fragment.vertices.end())
+ continue;
+
+ selected_vertex.push_back(v);
+
+ new_weight += get(l_state.weight_map, v);
+
+ auto edges = out_edges(v, l_state.graph);
+ for (auto eit = edges.first; eit != edges.second; eit++) {
+ const vertex &other = target(*eit, l_state.graph);
+ removed.insert(other);
+ }
+
+ for (size_t i : get(search_state.clique_map, v))
+ assigned.insert(i);
+ }
+
+ auto vs = boost::vertices(l_state.graph);
+ for (auto it = vs.first; it != vs.second; it++) {
+ vertex v = *it;
+ if (fragment.vertices.find(v) == fragment.vertices.end())
+ removed.insert(v);
+ }
+
+ std::map<size_t, size_t> size_delta;
+ for (vertex v : removed) {
+ for (size_t clique : get(search_state.clique_map, v)) {
+ size_delta[clique]++;
+ if (size_delta[clique] == search_state.clique_sizes[clique]) {
+ assigned.insert(clique);
+ }
+ }
+ }
+
+ std::vector<size_t> assigned_v(assigned.begin(), assigned.end());
+ std::vector<vertex> removed_v(removed.begin(), removed.end());
+
+ state_change<Graph, WeightMap> delta(
+ assigned_v, size_delta,
+ selected_vertex, removed_v,
+ weight(0), new_weight);
+
+ delta.apply(search_state);
+
+ lagrangian_bound<Graph, WeightMap>
+ bound(l_state.graph, l_state.weight_map,
+ _cliques.begin(), _cliques.end());
+
+ visit_state<Graph, WeightMap> visitor;
+ visitor.best_solution = l_state.solution;
+ visitor.best_value = weight(0);
+ for (vertex v : visitor.best_solution)
+ visitor.best_value += get(l_state.weight_map, v);
+
+ successor<Graph, WeightMap, decltype(bound)> successor(
+ bound, visitor.best_value);
+
+ search::depth_first_search(search_state, visitor, successor);
+
+ std::vector<vertex> newly_selected;
+ for (vertex v : visitor.best_solution) {
+ if (fragment.vertices.find(v) != fragment.vertices.end())
+ newly_selected.push_back(v);
+ }
+
+ return lns_assignment<Graph>(newly_selected);
+ }
+};
+
+template<typename Graph, typename WeightMap>
+class direct_evaluator {
+public:
+ typedef typename boost::property_traits<WeightMap>::value_type result_type;
+ typedef state<Graph, WeightMap> argument_type;
+
+private:
+ result_type _bound;
+
+public:
+ direct_evaluator(result_type bound): _bound(bound) {}
+
+ result_type operator()(const argument_type &state) const {
+ return state.solution_value / _bound;
+ }
+};
+
+template<typename Graph, typename OutputIterator>
+void find_cliques(const Graph &graph, OutputIterator out) {
+ typedef typename boost::graph_traits<Graph>::vertex_descriptor vertex;
+
+ std::set<std::pair<vertex, vertex>> visited_edges;
+
+ auto es = edges(graph);
+ for (auto eit = es.first; eit != es.second; eit++) {
+ auto u = source(*eit, graph);
+ auto v = target(*eit, graph);
+
+ if (visited_edges.find(std::make_pair(u, v)) == visited_edges.end()) {
+ std::vector<vertex> contents, candidates;
+
+ auto connected_es = out_edges(u, graph);
+ for (auto eit = connected_es.first; eit != connected_es.second; eit++) {
+ auto new_v = target(*eit, graph);
+ if (new_v != v && visited_edges.find(std::make_pair(u, new_v)) ==
+ visited_edges.end()) {
+ candidates.push_back(new_v);
+ }
+ }
+
+ candidates.push_back(v);
+
+ contents.push_back(u);
+
+ while (!candidates.empty()) {
+ vertex v = candidates.back();
+ candidates.pop_back();
+
+ contents.push_back(v);
+
+ candidates.erase(
+ std::remove_if(candidates.begin(), candidates.end(),
+ [&](vertex w) {
+ return visited_edges.find(std::make_pair(v, w)) ==
+ visited_edges.end();
+ }),
+ candidates.end());
+ }
+
+ for (size_t i = 0; i < contents.size(); i++) {
+ for (size_t j = 0; j < i; j++) {
+ visited_edges.insert(std::make_pair(contents[i], contents[j]));
+ visited_edges.insert(std::make_pair(contents[j], contents[i]));
+ }
+ }
+
+ *out++ = contents;
+ }
+ }
+}
+
+}
+
+template<typename Graph, typename WeightMap, typename OutputIterator>
+void maximum_weight_independent_set(const Graph &graph, WeightMap weight_map,
+ OutputIterator out) {
+ typedef typename boost::graph_traits<Graph>::vertex_descriptor vertex;
+ typedef typename boost::property_traits<WeightMap>::value_type weight;
+
+ std::vector<std::vector<vertex>> cliques;
+ mwis::find_cliques(graph, std::back_inserter(cliques));
+
+ mwis::lagrangian_bound<Graph, WeightMap> bound(
+ graph, weight_map,
+ cliques.begin(), cliques.end());
+
+ bound.check_consistency();
+
+ mwis::visit_state<Graph, WeightMap> visitor;
+ mwis::successor<Graph, WeightMap, decltype(bound)> successor(
+ bound, visitor.best_value);
+
+ mwis::state<Graph, WeightMap> state(
+ graph, weight_map,
+ cliques.begin(), cliques.end());
+
+ // search::depth_first_search(state, visitor, successor);
+ // for (vertex v : visitor.best_solution)
+ // *out++ = v;
+
+ auto vs = vertices(graph);
+ weight root_bound(bound(vs.second, vs.second, vs.first, vs.second));
+
+ mwis::direct_evaluator<Graph, WeightMap> direct_eval(root_bound);
+
+ mwis::max_bound<weight> max_bound;
+ mwis::successor<Graph, WeightMap, decltype(max_bound)>
+ max_successor(max_bound, visitor.best_value);
+
+ search::monte_carlo_tree_search(state, max_successor, direct_eval);
+ std::cout << "\n";
+
+ for (vertex v : state.solution)
+ *out++ = v;
+
+ // mwis::evaluator<Graph, WeightMap> eval(cliques.begin(), cliques.end(), 5);
+
+ // visitor(state);
+ // search::greedy_search(state, successor, eval);
+
+ // mwis::lns_state<Graph, WeightMap> l_state(graph, weight_map, state.solution);
+ // mwis::lns_fragment<Graph> (*l_selector)(const decltype(l_state)&) =
+ // mwis::fragment_selector<Graph, WeightMap>;
+ // mwis::lns_search<Graph, WeightMap> l_search(cliques.begin(), cliques.end());
+
+ // search::large_neighborhood_search(l_state, l_selector, l_search, 300);
+ // std::cout << "\n";
+
+ // for (vertex v : l_state.solution)
+ // *out++ = v;
+}
--
GitLab