1// Copyright (C) 2006-2009 Dmitry Bufistov and Andrey Parfenov
2
3// Use, modification and distribution is subject to the Boost Software
4// License, Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
5// http://www.boost.org/LICENSE_1_0.txt)
6
7#ifndef BOOST_GRAPH_CYCLE_RATIO_HOWARD_HPP
8#define BOOST_GRAPH_CYCLE_RATIO_HOWARD_HPP
9
10#include <vector>
11#include <list>
12#include <algorithm>
13#include <limits>
14
15#include <boost/bind.hpp>
16#include <boost/type_traits/is_same.hpp>
17#include <boost/type_traits/remove_const.hpp>
18#include <boost/concept_check.hpp>
19#include <boost/pending/queue.hpp>
20#include <boost/property_map/property_map.hpp>
21#include <boost/graph/graph_traits.hpp>
22#include <boost/graph/graph_concepts.hpp>
23#include <boost/concept/assert.hpp>
24
25/** @file howard_cycle_ratio.hpp
26 * @brief The implementation of the maximum/minimum cycle ratio/mean algorithm.
27 * @author Dmitry Bufistov
28 * @author Andrey Parfenov
29 */
30
31namespace boost {
32
33 /**
34 * The mcr_float is like numeric_limits, but only for floating point types
35 * and only defines infinity() and epsilon(). This class is primarily used
36 * to encapsulate a less-precise epsilon than natively supported by the
37 * floating point type.
38 */
39 template <typename Float = double> struct mcr_float {
40 typedef Float value_type;
41
42 static Float infinity()
43 { return std::numeric_limits<value_type>::infinity(); }
44
45 static Float epsilon()
46 { return Float(-0.005); }
47 };
48
49 namespace detail {
50
51 template <typename FloatTraits> struct
52 min_comparator_props {
53 typedef std::greater<typename FloatTraits::value_type> comparator;
54 static const int multiplier = 1;
55 };
56
57 template <typename FloatTraits> struct
58 max_comparator_props {
59 typedef std::less<typename FloatTraits::value_type> comparator;
60 static const int multiplier = -1;
61 };
62
63 template <typename FloatTraits, typename ComparatorProps>
64 struct float_wrapper {
65 typedef typename FloatTraits::value_type value_type;
66 typedef ComparatorProps comparator_props_t;
67 typedef typename ComparatorProps::comparator comparator;
68
69 static value_type infinity()
70 { return FloatTraits::infinity() * ComparatorProps::multiplier; }
71
72 static value_type epsilon()
73 { return FloatTraits::epsilon() * ComparatorProps::multiplier; }
74
75 };
76
77 /*! @class mcr_howard
78 * @brief Calculates optimum (maximum/minimum) cycle ratio of a directed graph.
79 * Uses Howard's iteration policy algorithm. </br>(It is described in the paper
80 * "Experimental Analysis of the Fastest Optimum Cycle Ratio and Mean Algorithm"
81 * by Ali Dasdan).
82 */
83 template <typename FloatTraits,
84 typename Graph, typename VertexIndexMap,
85 typename EdgeWeight1, typename EdgeWeight2>
86 class mcr_howard
87 {
88 public:
89 typedef typename FloatTraits::value_type float_t;
90 typedef typename FloatTraits::comparator_props_t cmp_props_t;
91 typedef typename FloatTraits::comparator comparator_t;
92 typedef enum{ my_white = 0, my_black } my_color_type;
93 typedef typename graph_traits<Graph>::vertex_descriptor vertex_t;
94 typedef typename graph_traits<Graph>::edge_descriptor edge_t;
95 typedef typename graph_traits<Graph>::vertices_size_type vn_t;
96 typedef std::vector<float_t> vp_t;
97 typedef typename boost::iterator_property_map<
98 typename vp_t::iterator, VertexIndexMap
99 > distance_map_t; //V -> float_t
100
101 typedef typename std::vector<edge_t> ve_t;
102 typedef std::vector<my_color_type> vcol_t;
103 typedef typename ::boost::iterator_property_map<
104 typename ve_t::iterator, VertexIndexMap
105 > policy_t; //Vertex -> Edge
106 typedef typename ::boost::iterator_property_map<
107 typename vcol_t::iterator, VertexIndexMap
108 > color_map_t;
109
110 typedef typename std::list<vertex_t> pinel_t;// The in_edges list of the policy graph
111 typedef typename std::vector<pinel_t> inedges1_t;
112 typedef typename ::boost::iterator_property_map<
113 typename inedges1_t::iterator, VertexIndexMap
114 > inedges_t;
115 typedef typename std::vector<edge_t> critical_cycle_t;
116
117 //Bad vertex flag. If true, then the vertex is "bad".
118 // Vertex is "bad" if its out_degree is equal to zero.
119 typedef typename boost::iterator_property_map<
120 std::vector<int>::iterator, VertexIndexMap
121 > badv_t;
122
123 /*!
124 * Constructor
125 * \param g = (V, E) - a directed multigraph.
126 * \param vim Vertex Index Map. Read property Map: V -> [0, num_vertices(g)).
127 * \param ewm edge weight map. Read property map: E -> R
128 * \param ew2m edge weight map. Read property map: E -> R+
129 * \param infty A big enough value to guaranty that there exist a cycle with
130 * better ratio.
131 * \param cmp The compare operator for float_ts.
132 */
133 mcr_howard(const Graph &g, VertexIndexMap vim,
134 EdgeWeight1 ewm, EdgeWeight2 ew2m) :
135 m_g(g), m_vim(vim), m_ew1m(ewm), m_ew2m(ew2m),
136 m_bound(mcr_bound()),
137 m_cr(m_bound),
138 m_V(num_vertices(m_g)),
139 m_dis(m_V, 0), m_dm(m_dis.begin(), m_vim),
140 m_policyc(m_V), m_policy(m_policyc.begin(), m_vim),
141 m_inelc(m_V), m_inel(m_inelc.begin(), m_vim),
142 m_badvc(m_V, false), m_badv(m_badvc.begin(), m_vim),
143 m_colcv(m_V),
144 m_col_bfs(m_V)
145 { }
146
147 /*!
148 * \return maximum/minimum_{for all cycles C}
149 * [sum_{e in C} w1(e)] / [sum_{e in C} w2(e)],
150 * or FloatTraits::infinity() if graph has no cycles.
151 */
152 float_t ocr_howard()
153 {
154 construct_policy_graph();
155 int k = 0;
156 float_t mcr = 0;
157 do
158 {
159 mcr = policy_mcr();
160 ++k;
161 }
162 while (try_improve_policy(cr: mcr) && k < 100); //To avoid infinite loop
163
164 const float_t eps_ = -0.00000001 * cmp_props_t::multiplier;
165 if (m_cmp(mcr, m_bound + eps_))
166 {
167 return FloatTraits::infinity();
168 }
169 else
170 {
171 return mcr;
172 }
173 }
174 virtual ~mcr_howard() {}
175
176 protected:
177 virtual void store_critical_edge(edge_t, critical_cycle_t &) {}
178 virtual void store_critical_cycle(critical_cycle_t &) {}
179
180 private:
181 /*!
182 * \return lower/upper bound for the maximal/minimal cycle ratio
183 */
184 float_t mcr_bound()
185 {
186 typename graph_traits<Graph>::vertex_iterator vi, vie;
187 typename graph_traits<Graph>::out_edge_iterator oei, oeie;
188 float_t cz = (std::numeric_limits<float_t>::max)(); //Closest to zero value
189 float_t s = 0;
190 const float_t eps_ = std::numeric_limits<float_t>::epsilon();
191 for (boost::tie(vi, vie) = vertices(m_g); vi != vie; ++vi)
192 {
193 for (boost::tie(oei, oeie) = out_edges(*vi, m_g); oei != oeie; ++oei)
194 {
195 s += std::abs(m_ew1m[*oei]);
196 float_t a = std::abs(m_ew2m[*oei]);
197 if ( a > eps_ && a < cz)
198 {
199 cz = a;
200 }
201 }
202 }
203 return cmp_props_t::multiplier * (s / cz);
204 }
205
206
207 /*!
208 * Constructs an arbitrary policy graph.
209 */
210 void construct_policy_graph()
211 {
212 m_sink = graph_traits<Graph>().null_vertex();
213 typename graph_traits<Graph>::vertex_iterator vi, vie;
214 typename graph_traits<Graph>::out_edge_iterator oei, oeie;
215 for ( boost::tie(vi, vie) = vertices(m_g); vi != vie; ++vi )
216 {
217 boost::tie(oei, oeie) = out_edges(*vi, m_g);
218 typename graph_traits<Graph>::out_edge_iterator mei =
219 std::max_element(oei, oeie,
220 boost::bind(m_cmp,
221 boost::bind(&EdgeWeight1::operator[], m_ew1m, _1),
222 boost::bind(&EdgeWeight1::operator[], m_ew1m, _2)
223 )
224 );
225 if (mei == oeie)
226 {
227 if (m_sink == graph_traits<Graph>().null_vertex())
228 {
229 m_sink = *vi;
230 }
231 m_badv[*vi] = true;
232 m_inel[m_sink].push_back(*vi);
233 }
234 else
235 {
236 m_inel[target(*mei, m_g)].push_back(*vi);
237 m_policy[*vi] = *mei;
238 }
239 }
240 }
241 /*! Sets the distance value for all vertices "v" such that there is
242 * a path from "v" to "sv". It does "inverse" breadth first visit of the policy
243 * graph, starting from the vertex "sv".
244 */
245 void mcr_bfv(vertex_t sv, float_t cr, color_map_t c)
246 {
247 boost::queue<vertex_t> Q;
248 c[sv] = my_black;
249 Q.push(sv);
250 while (!Q.empty())
251 {
252 vertex_t v = Q.top(); Q.pop();
253 for (typename pinel_t::const_iterator itr = m_inel[v].begin();
254 itr != m_inel[v].end(); ++itr)
255 //For all in_edges of the policy graph
256 {
257 if (*itr != sv)
258 {
259 if (m_badv[*itr])
260 {
261 m_dm[*itr] = m_dm[v] + m_bound - cr;
262 }
263 else
264 {
265 m_dm[*itr] = m_dm[v] + m_ew1m[m_policy[*itr]] -
266 m_ew2m[m_policy[*itr]] * cr;
267 }
268 c[*itr] = my_black;
269 Q.push(*itr);
270 }
271 }
272 }
273 }
274
275 /*!
276 * \param sv an arbitrary (undiscovered) vertex of the policy graph.
277 * \return a vertex in the policy graph that belongs to a cycle.
278 * Performs a depth first visit until a cycle edge is found.
279 */
280 vertex_t find_cycle_vertex(vertex_t sv)
281 {
282 vertex_t gv = sv;
283 std::fill(m_colcv.begin(), m_colcv.end(), my_white);
284 color_map_t cm(m_colcv.begin(), m_vim);
285 do
286 {
287 cm[gv] = my_black;
288 if (! m_badv[gv])
289 {
290 gv = target(m_policy[gv], m_g);
291 }
292 else
293 {
294 gv = m_sink;
295 }
296 }
297 while (cm[gv] != my_black);
298 return gv;
299 }
300
301 /*!
302 * \param sv - vertex that belongs to a cycle in the policy graph.
303 */
304 float_t cycle_ratio(vertex_t sv)
305 {
306 if (sv == m_sink) return m_bound;
307 std::pair<float_t, float_t> sums_(float_t(0), float_t(0));
308 vertex_t v = sv;
309 critical_cycle_t cc;
310 do
311 {
312 store_critical_edge(m_policy[v], cc);
313 sums_.first += m_ew1m[m_policy[v]];
314 sums_.second += m_ew2m[m_policy[v]];
315 v = target(m_policy[v], m_g);
316 }
317 while (v != sv);
318 float_t cr = sums_.first / sums_.second;
319 if ( m_cmp(m_cr, cr) )
320 {
321 m_cr = cr;
322 store_critical_cycle(cc);
323 }
324 return cr;
325 }
326
327 /*!
328 * Finds the optimal cycle ratio of the policy graph
329 */
330 float_t policy_mcr()
331 {
332 std::fill(m_col_bfs.begin(), m_col_bfs.end(), my_white);
333 color_map_t vcm_ = color_map_t(m_col_bfs.begin(), m_vim);
334 typename graph_traits<Graph>::vertex_iterator uv_itr, vie;
335 boost::tie(uv_itr, vie) = vertices(m_g);
336 float_t mcr = m_bound;
337 while ( (uv_itr = std::find_if(uv_itr, vie,
338 boost::bind(std::equal_to<my_color_type>(),
339 my_white,
340 boost::bind(&color_map_t::operator[], vcm_, _1)
341 )
342 )
343 ) != vie )
344 ///While there are undiscovered vertices
345 {
346 vertex_t gv = find_cycle_vertex(sv: *uv_itr);
347 float_t cr = cycle_ratio(sv: gv) ;
348 mcr_bfv(sv: gv, cr, c: vcm_);
349 if ( m_cmp(mcr, cr) ) mcr = cr;
350 ++uv_itr;
351 }
352 return mcr;
353 }
354
355 /*!
356 * Changes the edge m_policy[s] to the new_edge.
357 */
358 void improve_policy(vertex_t s, edge_t new_edge)
359 {
360 vertex_t t = target(m_policy[s], m_g);
361 typename property_traits<VertexIndexMap>::value_type ti = m_vim[t];
362 m_inelc[ti].erase( std::find(m_inelc[ti].begin(), m_inelc[ti].end(), s));
363 m_policy[s] = new_edge;
364 t = target(new_edge, m_g);
365 m_inel[t].push_back(s); ///Maintain in_edge list
366 }
367
368 /*!
369 * A negative cycle detector.
370 */
371 bool try_improve_policy(float_t cr)
372 {
373 bool improved = false;
374 typename graph_traits<Graph>::vertex_iterator vi, vie;
375 typename graph_traits<Graph>::out_edge_iterator oei, oeie;
376 const float_t eps_ = FloatTraits::epsilon();
377 for (boost::tie(vi, vie) = vertices(m_g); vi != vie; ++vi)
378 {
379 if (!m_badv[*vi])
380 {
381 for (boost::tie(oei, oeie) = out_edges(*vi, m_g); oei != oeie; ++oei)
382 {
383 vertex_t t = target(*oei, m_g);
384 //Current distance from *vi to some vertex
385 float_t dis_ = m_ew1m[*oei] - m_ew2m[*oei] * cr + m_dm[t];
386 if ( m_cmp(m_dm[*vi] + eps_, dis_) )
387 {
388 improve_policy(s: *vi, new_edge: *oei);
389 m_dm[*vi] = dis_;
390 improved = true;
391 }
392 }
393 }
394 else
395 {
396 float_t dis_ = m_bound - cr + m_dm[m_sink];
397 if ( m_cmp(m_dm[*vi] + eps_, dis_) )
398 {
399 m_dm[*vi] = dis_;
400 }
401 }
402 }
403 return improved;
404 }
405 private:
406 const Graph &m_g;
407 VertexIndexMap m_vim;
408 EdgeWeight1 m_ew1m;
409 EdgeWeight2 m_ew2m;
410 comparator_t m_cmp;
411 float_t m_bound; //> The lower/upper bound to the maximal/minimal cycle ratio
412 float_t m_cr; //>The best cycle ratio that has been found so far
413
414 vn_t m_V; //>The number of the vertices in the graph
415 vp_t m_dis; //>Container for the distance map
416 distance_map_t m_dm; //>Distance map
417
418 ve_t m_policyc; //>Container for the policy graph
419 policy_t m_policy; //>The interface for the policy graph
420
421 inedges1_t m_inelc; //>Container fot in edges list
422 inedges_t m_inel; //>Policy graph, input edges list
423
424 std::vector<int> m_badvc;
425 badv_t m_badv; //Marks "bad" vertices
426
427 vcol_t m_colcv, m_col_bfs; //Color maps
428 vertex_t m_sink; //To convert any graph to "good"
429 };
430
431 /*! \class mcr_howard1
432 * \brief Finds optimum cycle raio and a critical cycle
433 */
434 template <typename FloatTraits,
435 typename Graph, typename VertexIndexMap,
436 typename EdgeWeight1, typename EdgeWeight2>
437 class mcr_howard1 : public
438 mcr_howard<FloatTraits, Graph, VertexIndexMap,
439 EdgeWeight1, EdgeWeight2>
440 {
441 public:
442 typedef mcr_howard<FloatTraits, Graph, VertexIndexMap,
443 EdgeWeight1, EdgeWeight2> inhr_t;
444 mcr_howard1(const Graph &g, VertexIndexMap vim,
445 EdgeWeight1 ewm, EdgeWeight2 ew2m) :
446 inhr_t(g, vim, ewm, ew2m)
447 { }
448
449 void get_critical_cycle(typename inhr_t::critical_cycle_t &cc)
450 { return cc.swap(m_cc); }
451
452 protected:
453 void store_critical_edge(typename inhr_t::edge_t ed,
454 typename inhr_t::critical_cycle_t &cc)
455 { cc.push_back(ed); }
456
457 void store_critical_cycle(typename inhr_t::critical_cycle_t &cc)
458 { m_cc.swap(cc); }
459
460 private:
461 typename inhr_t::critical_cycle_t m_cc; //Critical cycle
462 };
463
464 /*!
465 * \param g a directed multigraph.
466 * \param vim Vertex Index Map. A map V->[0, num_vertices(g))
467 * \param ewm Edge weight1 map.
468 * \param ew2m Edge weight2 map.
469 * \param pcc pointer to the critical edges list.
470 * \return Optimum cycle ratio of g or FloatTraits::infinity() if g has no cycles.
471 */
472 template <typename FT,
473 typename TG, typename TVIM,
474 typename TEW1, typename TEW2,
475 typename EV>
476 typename FT::value_type
477 optimum_cycle_ratio(const TG &g, TVIM vim, TEW1 ewm, TEW2 ew2m, EV* pcc)
478 {
479 typedef typename graph_traits<TG>::directed_category DirCat;
480 BOOST_STATIC_ASSERT((is_convertible<DirCat*, directed_tag*>::value == true));
481 BOOST_CONCEPT_ASSERT(( IncidenceGraphConcept<TG> ));
482 BOOST_CONCEPT_ASSERT(( VertexListGraphConcept<TG> ));
483 typedef typename graph_traits<TG>::vertex_descriptor Vertex;
484 BOOST_CONCEPT_ASSERT(( ReadablePropertyMapConcept<TVIM, Vertex> ));
485 typedef typename graph_traits<TG>::edge_descriptor Edge;
486 BOOST_CONCEPT_ASSERT(( ReadablePropertyMapConcept<TEW1, Edge> ));
487 BOOST_CONCEPT_ASSERT(( ReadablePropertyMapConcept<TEW2, Edge> ));
488
489 if(pcc == 0) {
490 return detail::mcr_howard<FT,TG, TVIM, TEW1, TEW2>(
491 g, vim, ewm, ew2m
492 ).ocr_howard();
493 }
494
495 detail::mcr_howard1<FT, TG, TVIM, TEW1, TEW2> obj(g, vim, ewm, ew2m);
496 double ocr = obj.ocr_howard();
497 obj.get_critical_cycle(*pcc);
498 return ocr;
499 }
500 } // namespace detail
501
502// Algorithms
503// Maximum Cycle Ratio
504
505template <
506 typename FloatTraits,
507 typename Graph,
508 typename VertexIndexMap,
509 typename EdgeWeight1Map,
510 typename EdgeWeight2Map>
511inline typename FloatTraits::value_type
512maximum_cycle_ratio(const Graph &g, VertexIndexMap vim, EdgeWeight1Map ew1m,
513 EdgeWeight2Map ew2m,
514 std::vector<typename graph_traits<Graph>::edge_descriptor>* pcc = 0,
515 FloatTraits = FloatTraits())
516{
517 typedef detail::float_wrapper<
518 FloatTraits, detail::max_comparator_props<FloatTraits>
519 > Traits;
520 return detail::optimum_cycle_ratio<Traits>(g, vim, ew1m, ew2m, pcc);
521}
522
523template <
524 typename Graph,
525 typename VertexIndexMap,
526 typename EdgeWeight1Map,
527 typename EdgeWeight2Map>
528inline double
529maximum_cycle_ratio(const Graph &g, VertexIndexMap vim,
530 EdgeWeight1Map ew1m, EdgeWeight2Map ew2m,
531 std::vector<typename graph_traits<Graph>::edge_descriptor>* pcc = 0)
532{ return maximum_cycle_ratio(g, vim, ew1m, ew2m, pcc, mcr_float<>()); }
533
534// Minimum Cycle Ratio
535
536template <
537 typename FloatTraits,
538 typename Graph,
539 typename VertexIndexMap,
540 typename EdgeWeight1Map,
541 typename EdgeWeight2Map>
542typename FloatTraits::value_type
543minimum_cycle_ratio(const Graph &g, VertexIndexMap vim,
544 EdgeWeight1Map ew1m, EdgeWeight2Map ew2m,
545 std::vector<typename graph_traits<Graph>::edge_descriptor> *pcc = 0,
546 FloatTraits = FloatTraits())
547{
548 typedef detail::float_wrapper<
549 FloatTraits, detail::min_comparator_props<FloatTraits>
550 > Traits;
551 return detail::optimum_cycle_ratio<Traits>(g, vim, ew1m, ew2m, pcc);
552}
553
554template <
555 typename Graph,
556 typename VertexIndexMap,
557 typename EdgeWeight1Map,
558 typename EdgeWeight2Map>
559inline double
560minimum_cycle_ratio(const Graph &g, VertexIndexMap vim,
561 EdgeWeight1Map ew1m, EdgeWeight2Map ew2m,
562 std::vector<typename graph_traits<Graph>::edge_descriptor>* pcc = 0)
563{ return minimum_cycle_ratio(g, vim, ew1m, ew2m, pcc, mcr_float<>()); }
564
565// Maximum Cycle Mean
566
567template <
568 typename FloatTraits,
569 typename Graph,
570 typename VertexIndexMap,
571 typename EdgeWeightMap,
572 typename EdgeIndexMap>
573inline typename FloatTraits::value_type
574maximum_cycle_mean(const Graph &g, VertexIndexMap vim,
575 EdgeWeightMap ewm, EdgeIndexMap eim,
576 std::vector<typename graph_traits<Graph>::edge_descriptor>* pcc = 0,
577 FloatTraits ft = FloatTraits())
578{
579 typedef typename remove_const<
580 typename property_traits<EdgeWeightMap>::value_type
581 >::type Weight;
582 typename std::vector<Weight> ed_w2(boost::num_edges(g), 1);
583 return maximum_cycle_ratio(g, vim, ewm,
584 make_iterator_property_map(ed_w2.begin(), eim),
585 pcc, ft);
586}
587
588template <
589 typename Graph,
590 typename VertexIndexMap,
591 typename EdgeWeightMap,
592 typename EdgeIndexMap>
593inline double
594maximum_cycle_mean(const Graph& g, VertexIndexMap vim,
595 EdgeWeightMap ewm, EdgeIndexMap eim,
596 std::vector<typename graph_traits<Graph>::edge_descriptor>* pcc = 0)
597{ return maximum_cycle_mean(g, vim, ewm, eim, pcc, mcr_float<>()); }
598
599// Minimum Cycle Mean
600
601template <
602 typename FloatTraits,
603 typename Graph,
604 typename VertexIndexMap,
605 typename EdgeWeightMap,
606 typename EdgeIndexMap>
607inline typename FloatTraits::value_type
608minimum_cycle_mean(const Graph &g, VertexIndexMap vim,
609 EdgeWeightMap ewm, EdgeIndexMap eim,
610 std::vector<typename graph_traits<Graph>::edge_descriptor>* pcc = 0,
611 FloatTraits ft = FloatTraits())
612{
613 typedef typename remove_const<
614 typename property_traits<EdgeWeightMap>::value_type
615 >::type Weight;
616 typename std::vector<Weight> ed_w2(boost::num_edges(g), 1);
617 return minimum_cycle_ratio(g, vim, ewm,
618 make_iterator_property_map(ed_w2.begin(), eim),
619 pcc, ft);
620}
621
622template <
623 typename Graph,
624 typename VertexIndexMap,
625 typename EdgeWeightMap,
626 typename EdgeIndexMap>
627inline double
628minimum_cycle_mean(const Graph &g, VertexIndexMap vim,
629 EdgeWeightMap ewm, EdgeIndexMap eim,
630 std::vector<typename graph_traits<Graph>::edge_descriptor>* pcc = 0)
631{ return minimum_cycle_mean(g, vim, ewm, eim, pcc, mcr_float<>()); }
632
633} //namespace boost
634
635#endif
636

source code of boost/boost/graph/howard_cycle_ratio.hpp