All Classes Namespaces Files Functions Variables Typedefs Enumerations Friends Macros Pages
zelikovsky_11_per_6.hpp
Go to the documentation of this file.
1 //=======================================================================
2 // Copyright (c) 2013 Piotr Wygocki
3 //
4 // Distributed under the Boost Software License, Version 1.0. (See
5 // accompanying file LICENSE_1_0.txt or copy at
6 // http://www.boost.org/LICENSE_1_0.txt)
7 //=======================================================================
16 #ifndef PAAL_ZELIKOVSKY_11_PER_6_HPP
17 #define PAAL_ZELIKOVSKY_11_PER_6_HPP
18 
19 #include <boost/config.hpp>
20 
28 #include "paal/utils/functors.hpp"
29 #include "paal/utils/irange.hpp"
30 
31 #include <boost/graph/adjacency_list.hpp>
32 #include <boost/graph/connected_components.hpp>
33 #include <boost/graph/subgraph.hpp>
34 #include <boost/graph/prim_minimum_spanning_tree.hpp>
35 
36 #include <boost/range/combine.hpp>
37 #include <boost/range/join.hpp>
38 #include <boost/range/algorithm/copy.hpp>
39 
40 #include <boost/iterator/zip_iterator.hpp>
41 #include <boost/iterator/transform_iterator.hpp>
42 
43 #include <boost/functional/hash.hpp>
44 
45 #include <unordered_map>
46 #include <stack>
47 
48 namespace paal {
49 
50 namespace detail {
66 template <typename Metric, typename Voronoi> class steiner_tree {
67  using Idx = int;
68 
69  public:
71  using Dist = typename MT::DistanceType;
72  using VertexType = typename MT::VertexType;
73  static const int SUBSET_SIZE = 3;
74 
75  using ThreeTuple = k_tuple_t<Idx, SUBSET_SIZE>;
76  using Move = boost::tuple<ThreeTuple, Dist>;
77  using ResultSteinerVertices = std::vector<VertexType>;
78 
85  steiner_tree(const Metric &m, const Voronoi &vor)
86  : m_metric(m), m_voronoi(vor),
87  N(boost::distance(vor.get_generators())),
88  m_save(N) {}
89 
90  template <typename OutputIterator>
91  void get_result_steiner_vertices(OutputIterator out) {
92  ResultSteinerVertices res;
93 
94  if (m_voronoi.get_vertices().empty()) {
95  return;
96  }
97 
98  auto ti = irange(N);
99  auto subsets = make_three_subset_range(ti.begin(), ti.end());
100 
101  auto get_moves = [&](const AMatrix &) {
102  return boost::combine(subsets, m_subs_dists | boost::adaptors::transformed(
104  };
105 
106  auto obj_fun = [&](const AMatrix & m, const Move & t) {
107  return gain(t);
108  };
109 
110  auto commit_move = [&](AMatrix & m, const Move & t) {
111  contract(m, get<0>(t));
112  res.push_back(m_nearest_vertex[get<0>(t)]);
113  return true;
114  };
115 
117  get_moves, obj_fun, commit_move);
118 
120  m_metric, m_voronoi.get_generators(), m_t_idx);
121 
122  fill_sub_dists();
123 
124  find_save(ls_solution);
125  local_search::local_search(ls_solution,
126  local_search::best_improving_strategy{},
127  [ = ](AMatrix & a) {
128  find_save(a);
129  return true;
130  },
132 
133  unique_res(res);
134  boost::copy(res, out);
135  }
136 
137  private:
138 
139  // Spanning tree types
140  using SpanningTreeEdgeProp = boost::property<boost::edge_index_t, int,
141  boost::property<boost::edge_weight_t, Dist>>;
142  using SpanningTree = boost::subgraph<boost::adjacency_list<
143  boost::listS, boost::vecS, boost::undirectedS, boost::no_property,
144  SpanningTreeEdgeProp>>;
145  using GTraits = boost::graph_traits<SpanningTree>;
146  using SEdge = typename GTraits::edge_descriptor;
147 
148  // Adjacency Matrix types
149  using AMatrix = typename data_structures::adjacency_matrix<Metric>::type;
150  using MTraits = boost::graph_traits<AMatrix>;
151  using MEdge = typename MTraits::edge_descriptor;
152 
153  // other types
154  using ThreeSubsetsDists = std::vector<Dist>;
155  using NearstByThreeSubsets = std::unordered_map<ThreeTuple, VertexType,
156  boost::hash<ThreeTuple>>;
157 
158  template <typename Iter>
159  boost::iterator_range<
160  data_structures::subsets_iterator<SUBSET_SIZE, Iter>>
161  make_three_subset_range(Iter b, Iter e) {
162  return data_structures::make_subsets_iterator_range<SUBSET_SIZE>(b, e);
163  }
164 
165  void unique_res(ResultSteinerVertices &res) {
166  std::sort(res.begin(), res.end());
167  auto new_end = std::unique(res.begin(), res.end());
168  res.resize(std::distance(res.begin(), new_end));
169  }
170 
171  void contract(AMatrix &am, const ThreeTuple &t) {
172  utils::contract(am, get<0>(t), std::get<1>(t));
173  utils::contract(am, get<1>(t), std::get<2>(t));
174  }
175 
176  Dist gain(const Move &t) {
177  auto const &m = m_save;
178  Idx a, b, c;
179  std::tie(a, b, c) = get<0>(t);
180 
181  assert(m(a, b) == m(b, c) || m(b, c) == m(c, a) || m(c, a) == m(a, b));
182  return max3(m(a, b), m(b, c), m(c, a)) +
183  min3(m(a, b), m(b, c), m(c, a)) - get<1>(t);
184  }
185 
186  void fill_sub_dists() {
187  auto ti = irange(N);
188 
189  auto sub_range = make_three_subset_range(ti.begin(), ti.end());
190  m_subs_dists.reserve(boost::distance(sub_range));
191 
192  // finding nearest vertex to subset
193  for (const ThreeTuple &subset : sub_range) {
194  // TODO awful coding, need to be changed to loop (using fold
195  // form utils/fusion.hpp)
196  // TODO There is possible problem, one point could belong to two
197  // voronoi regions
198  // In our implementation the point will be in exactly one region and
199  // there
200  // it will not be contained in the range
201  auto v_range1 = m_voronoi.get_vertices_for_generator(
202  m_t_idx.get_val(std::get<0>(subset)));
203  auto v_range2 = m_voronoi.get_vertices_for_generator(
204  m_t_idx.get_val(std::get<1>(subset)));
205  auto v_range3 = m_voronoi.get_vertices_for_generator(
206  m_t_idx.get_val(std::get<2>(subset)));
207  auto range = boost::join(boost::join(v_range1, v_range2), v_range3);
208 
209  if (boost::empty(range)) {
210  m_nearest_vertex[subset] = *m_voronoi.get_vertices().begin();
211  } else {
212  m_nearest_vertex[subset] = *std::min_element(
213  std::begin(range), std::end(range),
214  utils::make_functor_to_comparator([&](VertexType v) {
215  return dist(v, subset);
216  }));
217  }
218  m_subs_dists.push_back(dist(m_nearest_vertex[subset], subset));
219  }
220  }
221 
222  Dist max3(Dist a, Dist b, Dist c) { return std::max(std::max(a, b), c); }
223 
224  Dist min3(Dist a, Dist b, Dist c) { return std::min(std::min(a, b), c); }
225 
226  Dist dist(VertexType steiner_point, Idx terminal_idx) {
227  return m_metric(steiner_point, m_t_idx.get_val(terminal_idx));
228  }
229 
230  // minor TODO could by more general somewhere
231  Dist dist(VertexType steiner_point, const ThreeTuple &tup) {
232  return dist(steiner_point, std::get<0>(tup)) +
233  dist(steiner_point, std::get<1>(tup)) +
234  dist(steiner_point, std::get<2>(tup));
235  }
236 
242  SpanningTree get_spanning_tree(const AMatrix &am) {
243  // compute spanning tree and write it to vector
244  std::vector<Idx> pm(N);
245  boost::prim_minimum_spanning_tree(am, &pm[0]);
246 
247  // transform vector into SpanningTree object
248  auto const &weight_map = get(boost::edge_weight, am);
249  SpanningTree spanning_tree(N);
250  for (Idx from = 0; from < N; ++from) {
251  if (from != pm[from]) {
252  bool succ = add_edge(
253  from, pm[from],
254  SpanningTreeEdgeProp(
255  from, get(weight_map, edge(from, pm[from], am).first)),
256  spanning_tree).second;
257  assert(succ);
258  }
259  }
260  return spanning_tree;
261  }
262 
263  template <typename WeightMap, typename EdgeRange>
264  SEdge max_edge(EdgeRange range, const WeightMap &weight_map) const {
265  assert(range.first != range.second);
266  return *std::max_element(range.first, range.second,
267  [&](SEdge e, SEdge f) {
268  return get(weight_map, e) < get(weight_map, f);
269  });
270  }
271 
272  void create_subgraphs(SpanningTree &g, SpanningTree &g1, SpanningTree &g2) {
273  int n = num_vertices(g);
274  std::vector<Idx> comps(n);
275  boost::connected_components(g, &comps[0]);
276  int c1 = comps[0];
277  int c2 = -1;
278 
279  for (auto i : irange(n)) {
280  if (comps[i] == c1) {
281  add_vertex(g.local_to_global(i), g1);
282  } else {
283  assert(c2 == -1 || comps[i] == c2);
284  c2 = comps[i];
285  add_vertex(g.local_to_global(i), g2);
286  }
287  }
288  }
289 
290  // setting m_save(v,w) = max_dist, for each v in g1 and w in g2
291  void move_save(const SpanningTree &g1, const SpanningTree &g2,
292  Dist max_dist) {
293  auto v1 = vertices(g1);
294  auto v2 = vertices(g2);
295  for (auto v : boost::make_iterator_range(v1)) {
296  for (auto w : boost::make_iterator_range(v2)) {
297  auto vg = g1.local_to_global(v);
298  auto wg = g2.local_to_global(w);
299  m_save(vg, wg) = max_dist;
300  m_save(wg, vg) = max_dist;
301  }
302  }
303  }
304 
305  // finds the longest edge between each pair of vertices
306  // in the spanning tree
307  // preforms recursive procedure
308  void find_save(const AMatrix &am) {
309  auto spanning_tree = get_spanning_tree(am);
310 
311  std::stack<SpanningTree *> s;
312  s.push(&spanning_tree);
313 
314  while (!s.empty()) {
315  // TODO delete children at once
316  SpanningTree &g = *s.top();
317  s.pop();
318  int n = num_vertices(g);
319  if (n == 1) {
320  continue;
321  }
322  auto e_range = edges(g);
323  assert(e_range.first != e_range.second);
324  auto const &weight_map = get(boost::edge_weight, g);
325  SEdge max_el = max_edge(e_range, weight_map);
326  Dist max_dist = get(weight_map, max_el);
327  remove_edge(max_el, g);
328  SpanningTree &g1 = g.create_subgraph();
329  SpanningTree &g2 = g.create_subgraph();
330  create_subgraphs(g, g1, g2);
331 
332  move_save(g1, g2, max_dist);
333 
334  s.push(&g1);
335  s.push(&g2);
336  }
337  }
338 
339  const Metric &m_metric;
340  const Voronoi &m_voronoi;
341  ThreeSubsetsDists m_subs_dists;
342  NearstByThreeSubsets m_nearest_vertex;
343  int N;
344  data_structures::array_metric<Dist> m_save;
345  data_structures::bimap<VertexType> m_t_idx;
346 };
347 } // !detail
348 
359 template <typename Metric, typename Voronoi, typename OutputIterator>
361  const Voronoi &v,
362  OutputIterator out) {
364  st.get_result_steiner_vertices(out);
365 }
366 
367 } // paal
368 
369 #endif // PAAL_ZELIKOVSKY_11_PER_6_HPP
adjacency_matrix< Metric >::type metric_to_bgl_with_index(const Metric &m, Vertices &&vertices, bimap< typename boost::range_value< Vertices >::type > &idx)
produces graph from metric with index
This is Alexander Zelikovsky 11/6 approximation algorithm for steiner tree.
functor return false
Definition: functors.hpp:222
auto make_search_components(Args &&...args)
make function for search components
removes reference
Definition: functors.hpp:257
const T & get_val(Idx i) const
get value for index i
Definition: bimap.hpp:166
steiner_tree(const Metric &m, const Voronoi &vor)
void steiner_tree_zelikovsky11per6approximation(const Metric &m, const Voronoi &v, OutputIterator out)
11/6 approximation for steiner_tree problem
typename components::type< Args...> search_components
search_components template alias
typename detail::k_tuple< T, k >::type k_tuple_t
returns tuple consisting of k times type T
auto make_functor_to_comparator(Functor functor, Compare compare=Compare())
make for functor to comparator
Definition: functors.hpp:654
auto irange(T begin, T end)
irange
Definition: irange.hpp:22
This file contains set of simple useful functors or functor adapters.
void contract(boost::adjacency_matrix< GraphArgs...> &amatrix, Idx v, Idx w)
contracts to vertices in adjacency_matrix
bool local_search(Solution &solution, SearchStrategy searchStrategy, ContinueOnSuccess succ, ContinueOnFail fail, components...comps)
detail
puretype(std::declval< Metric >()(std::declval< VertexType >(), std::declval< VertexType >())) DistanceType
Distance type.