boost/graph/howard_cycle_ratio.hpp
/*! * Copyright 2007 Technical University of Catalonia * * Use, modification and distribution is subject to the Boost Software * License, Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at * http://www.boost.org/LICENSE_1_0.txt) * * Authors: Dmitry Bufistov * Andrey Parfenov */ #ifndef BOOST_GRAPH_HOWARD_CYCLE_RATIO_HOWARD_HPP #define BOOST_GRAPH_HOWARD_CYCLE_RATIO_HOWARD_HPP /*! * \file Maximum cycle ratio algorithm (Jean Cochet-Terrasson, Guy * Cochen and others) */ #include <exception> #include <set> #include <boost/bind.hpp> #include <boost/lexical_cast.hpp> #include <boost/type_traits/is_convertible.hpp> #include <boost/type_traits/remove_const.hpp> #include <boost/type_traits/is_signed.hpp> #include <boost/concept_check.hpp> #include <boost/graph/adjacency_list.hpp> #include <boost/graph/reverse_graph.hpp> #include <boost/graph/breadth_first_search.hpp> #include <boost/graph/iteration_macros.hpp> namespace boost { namespace detail { /// To avoid round error. static const double mcr_howard_ltolerance = 0.00001; /*! * Calculate maximum cycle ratio of "good" directed multigraph * g. Use Howard's iteration policy algorithm ("Numerical * Computation of Spectral Elements in MAX-PLUS algebra" by Jean * Cochet-Terrasson, Guy Cochen and others). * * \param g = (V, E) - a "good" directed multigraph (out_degree of * each vertex is greater then 0). If graph is strongly connected * then it is "good". * * \param vim - Vertex Index, read property Map: V -> [0, * num_vertices(g)). * * \param ewm - edge weight read property map: E -> R * * \param ewm2 - edge weight2 read property map: E -> R+ * * \return maximum_{for all cycles C}CR(C), or * -(std::numeric_limits<double>)::max() if g is not "good". */ template <typename TGraph, typename TVertexIndexMap, typename TWeight1EdgeMap, typename TWeight2EdgeMap > class Cmcr_Howard { public: Cmcr_Howard(const TGraph& g, TVertexIndexMap vim, TWeight1EdgeMap ewm, TWeight2EdgeMap ew2m) : m_g(g), m_vim(vim), m_ew1m(ewm), m_ew2m(ew2m), m_g2pi_g_vm(std::vector<pi_vertex_t>().end(), m_vim), /// Stupid dummy initialization m_minus_infinity(-(std::numeric_limits<double>::max)()) { typedef typename boost::graph_traits<TGraph>::directed_category DirCat; BOOST_STATIC_ASSERT((boost::is_convertible<DirCat*, boost::directed_tag*>::value == true)); m_cr = m_minus_infinity; } double operator()() { return maximum_cycle_ratio_Howard(); } virtual ~Cmcr_Howard() { } protected: typedef typename boost::graph_traits<TGraph>::vertex_descriptor mcr_vertex_t; typedef typename boost::graph_traits<TGraph>::edge_descriptor mcr_edge_t; const TGraph& m_g; typedef std::vector<double> eigenmode_t; eigenmode_t m_eigen_value; eigenmode_t m_eigen_vector; TVertexIndexMap m_vim; TWeight1EdgeMap m_ew1m; TWeight2EdgeMap m_ew2m; typedef typename boost::remove_const<typename boost::property_traits<TWeight1EdgeMap>::value_type>::type mcr_edge_weight1_t; typedef typename boost::remove_const<typename boost::property_traits<TWeight2EdgeMap>::value_type>::type mcr_edge_weight2_t; typedef typename boost::adjacency_list< boost::listS, boost::vecS, boost::bidirectionalS, boost::no_property, boost::property<boost::edge_weight_t, mcr_edge_weight1_t, boost::property<boost::edge_weight2_t, mcr_edge_weight2_t> > > pi_graph_t; typedef typename boost::property_map<pi_graph_t, boost::vertex_index_t>::type TPiGraphVertexIndexMap; typedef typename boost::property_map<pi_graph_t, boost::edge_weight_t>::type TPiGraphEdgeWeight1Map; typedef typename boost::property_map<pi_graph_t, boost::edge_weight2_t>::type TPiGraphEdgeWeight2Map; typedef typename boost::property_traits<TPiGraphVertexIndexMap>::value_type pigraph_vertex_index_t; pi_graph_t m_pi_g; typedef typename boost::graph_traits<pi_graph_t>::vertex_descriptor pi_vertex_t; typedef typename boost::graph_traits<pi_graph_t>::edge_descriptor pi_edge_t; typedef typename boost::iterator_property_map<typename std::vector<pi_vertex_t>::iterator, TVertexIndexMap> g2pi_g_vm_t; g2pi_g_vm_t m_g2pi_g_vm; ///Graph to Pi graph vertex map std::vector<pi_vertex_t> m_g2pig; int m_step_number; const double m_minus_infinity; typedef typename std::vector<mcr_edge_t> critical_cycle_t; double m_cr; ///Cycle ratio that already has been found class bad_graph { public: typedef typename boost::property_traits<TVertexIndexMap>::value_type v_index_t; bad_graph(v_index_t bvi) : bad_vertex_index(bvi) {} v_index_t what() const throw() { return bad_vertex_index; } private: v_index_t bad_vertex_index; }; double maximum_cycle_ratio_Howard() { try { construct_pi_graph(); } catch (const bad_graph& a) { return m_minus_infinity; } std::vector<double> max_eigen_val(boost::num_vertices(m_g)); m_eigen_value.resize(boost::num_vertices(m_g)); m_eigen_vector.resize(boost::num_vertices(m_g)); m_step_number = 0; do { pi_eingen_value(get(vertex_index, m_pi_g), get(boost::edge_weight, m_pi_g), get(boost::edge_weight2, m_pi_g)); ++m_step_number; } while (improve_policy_try1(max_eigen_val) || improve_policy_try2(max_eigen_val)); return *(std::max_element(m_eigen_value.begin(), m_eigen_value.end())); } /*! * Construct an arbitrary policy m_pi_g. */ void construct_pi_graph() { m_g2pig.resize(boost::num_vertices(m_g)); m_g2pi_g_vm = boost::make_iterator_property_map(m_g2pig.begin(), m_vim); BGL_FORALL_VERTICES_T(vd, m_g, TGraph) { m_g2pi_g_vm[vd] = boost::add_vertex(m_pi_g); store_pivertex(m_g2pi_g_vm[vd], vd); } BGL_FORALL_VERTICES_T(vd1, m_g, TGraph) { if (boost::out_edges(vd1, m_g).first == boost::out_edges(vd1, m_g).second) throw bad_graph(m_vim[vd1]); mcr_edge_t ed = *boost::out_edges(vd1, m_g).first; pi_edge_t pied = boost::add_edge(m_g2pi_g_vm[source(ed, m_g)], m_g2pi_g_vm[target(ed, m_g)], m_pi_g).first; boost::put(boost::edge_weight, m_pi_g, pied, m_ew1m[ed]); boost::put(boost::edge_weight2, m_pi_g, pied, m_ew2m[ed]); } } class bfs_eingmode_visitor : public boost::default_bfs_visitor { public: bfs_eingmode_visitor(TPiGraphVertexIndexMap vi_m, TPiGraphEdgeWeight1Map w_m, TPiGraphEdgeWeight2Map& d_m, eigenmode_t& e_val, eigenmode_t& e_vec, double ev) : m_index_map(vi_m), m_weight_map(w_m), m_delay_map(d_m), m_eig_value(&e_val), m_eig_vec(&e_vec), m_eigen_value(ev) { } template < typename Edge, typename g_t> void examine_edge(Edge e, const g_t & g) const { typedef typename boost::graph_traits<g_t>::vertex_descriptor Vertex; Vertex u = boost::target(e, g), v = boost::source(e, g); pigraph_vertex_index_t ind = m_index_map[u]; (*m_eig_value)[ind] = m_eigen_value; (*m_eig_vec)[ind] = m_weight_map[e] - m_eigen_value * m_delay_map[e] + (*m_eig_vec)[m_index_map[v]]; } private: TPiGraphVertexIndexMap m_index_map; TPiGraphEdgeWeight1Map m_weight_map; TPiGraphEdgeWeight2Map m_delay_map; eigenmode_t* m_eig_value; eigenmode_t* m_eig_vec; double m_eigen_value; }; /*! * Find a vertex in the Pi Graph which belongs to cycle, just a DFV until back edge found */ pi_vertex_t find_good_source(const pi_vertex_t start_vertex) { pi_vertex_t good_vertex = start_vertex; typename std::set<pi_vertex_t> s; s.insert(start_vertex); do { good_vertex = boost::target(*boost::out_edges(good_vertex, m_pi_g).first, m_pi_g); } while (s.insert(good_vertex).second); return good_vertex; } virtual void store_pivertex(pi_vertex_t pivd, mcr_vertex_t vd) {} virtual void store_critical_edge(pi_edge_t ed, critical_cycle_t& cc) {} virtual void store_critical_cycle(critical_cycle_t& cc) {} /*! * \param startV - vertex that belongs to a cycle in policy graph m_pi_g */ double calculate_eigen_value(pi_vertex_t startV) { std::pair<double, double> accum_sums(0., 0.); pi_vertex_t vd = startV; critical_cycle_t cc; do { pi_edge_t tmp_ed = *(boost::out_edges(vd, m_pi_g).first); store_critical_edge(tmp_ed, cc); accum_sums.first += boost::get(boost::edge_weight, m_pi_g, tmp_ed); accum_sums.second += boost::get(boost::edge_weight2, m_pi_g, tmp_ed); vd = boost::target(tmp_ed, m_pi_g); } while (vd != startV); //assert((std::abs<double>(accum_sums.first) <= 0.00000001) && "Division by zerro!"); double cr = accum_sums.first / accum_sums.second; if (cr > m_cr) { m_cr = cr; store_critical_cycle(cc); } else { } return cr; } /*! * Value determination. Find a generalized eigenmode (n^{k+1}, x^{k+1}) of A^{�_{k+1}} of the pi graph (Algorithm IV.1). */ void pi_eingen_value( TPiGraphVertexIndexMap index_map, TPiGraphEdgeWeight1Map weight_map, TPiGraphEdgeWeight2Map weigh2_map) { using namespace boost; typedef std::vector<default_color_type> color_map_t; color_map_t vcm(num_vertices(m_pi_g), white_color);//Vertex color map color_map_t::iterator uv_itr = vcm.begin(); //Undiscovered vertex reverse_graph<pi_graph_t> rev_g(m_pi_g); //For backward breadth visit while ((uv_itr = std::find_if(uv_itr, vcm.end(), boost::bind(std::equal_to<default_color_type>(), boost::white_color, _1))) != vcm.end()) ///While there are undiscovered vertices { pi_vertex_t gv = find_good_source(pi_vertex_t(uv_itr - vcm.begin())); pigraph_vertex_index_t gv_ind = index_map[gv]; m_eigen_value[gv_ind] = calculate_eigen_value(gv) ; bfs_eingmode_visitor bfs_vis(index_map, weight_map, weigh2_map, m_eigen_value, m_eigen_vector, m_eigen_value[gv_ind]); typename boost::queue<pi_vertex_t> Q; breadth_first_visit(rev_g, gv, Q, bfs_vis, make_iterator_property_map(vcm.begin(), index_map)); } } void improve_policy(mcr_vertex_t vd, mcr_edge_t new_edge) { remove_edge(*(out_edges(m_g2pi_g_vm[vd], m_pi_g).first), m_pi_g); pi_edge_t ned = add_edge(m_g2pi_g_vm[vd], m_g2pi_g_vm[target(new_edge, m_g)], m_pi_g).first; put(edge_weight, m_pi_g, ned, m_ew1m[new_edge]); put(edge_weight2, m_pi_g, ned, m_ew2m[new_edge]); } /*! * Policy Improvement. Improve the policy graph. The new policy graph has greater cycle ratio. * \return false if nothing can be improved. */ bool improve_policy_try1(std::vector<double>& max_eing_vals) { bool improved = false; BGL_FORALL_VERTICES_T(vd, m_g, TGraph) { double max_ev = m_minus_infinity;/// Maximum eigen value for vertex mcr_edge_t cr_ed;///Critical edge BGL_FORALL_OUTEDGES_T(vd, outed, m_g, TGraph) { if (m_eigen_value[m_vim[target(outed, m_g)]] > max_ev) { max_ev = m_eigen_value[m_vim[boost::target(outed, m_g)]]; cr_ed = outed; } } if (max_ev > m_eigen_value[get(m_vim,vd)]) { improve_policy(vd, cr_ed); improved = true; } max_eing_vals[get(m_vim,vd)] = max_ev; } return improved; } /*! * \param max_eigen_values[u] = max_(for all adjacent vertices (u,v)) m_eigen_value[v] */ bool improve_policy_try2(const std::vector<double>& max_eigen_values) { bool improved = false; BGL_FORALL_VERTICES_T(vd, m_g, TGraph) { mcr_edge_t impr_edge; double max_val = m_minus_infinity; BGL_FORALL_OUTEDGES_T(vd, outed, m_g, TGraph) { ///If vertex vd is in the K(vd) set if (max_eigen_values[get(m_vim, vd)] <= m_eigen_value[get(m_vim, target(outed, m_g))]) { double c_val = m_ew1m[outed] - m_ew2m[outed] * m_eigen_value[m_vim[boost::target(outed, m_g)]] + m_eigen_vector[m_vim[boost::target(outed, m_g)]]; if (c_val > max_val) { max_val = c_val; impr_edge = outed; } } } if ((max_val - m_eigen_vector[get(m_vim, vd)]) > mcr_howard_ltolerance) ///If m_eigen_vector[vd] == max_val { improve_policy(vd, impr_edge); improved = true; } } return improved; } };///Cmcr_Howard /*! * \return maximum cycle ratio and one critical cycle. */ template <typename TGraph, typename TVertexIndexMap, typename TWeight1EdgeMap, typename TWeight2EdgeMap> class Cmcr_Howard1 : public Cmcr_Howard<TGraph, TVertexIndexMap, TWeight1EdgeMap, TWeight2EdgeMap> { public: typedef Cmcr_Howard<TGraph, TVertexIndexMap, TWeight1EdgeMap, TWeight2EdgeMap> inhr_t; Cmcr_Howard1(const TGraph& g, TVertexIndexMap vim, TWeight1EdgeMap ewm, TWeight2EdgeMap ew2m) : inhr_t(g, vim, ewm, ew2m) { m_pi_g2g.resize(boost::num_vertices(g)); m_pi_g2g_vm = boost::make_iterator_property_map(m_pi_g2g.begin(), boost::get(boost::vertex_index, this->m_pi_g)); } void get_critical_cycle(typename inhr_t::critical_cycle_t& cc) { return cc.swap(m_critical_cycle); } protected: void store_pivertex(typename inhr_t::pi_vertex_t pivd, typename inhr_t::mcr_vertex_t vd) { m_pi_g2g_vm[pivd] = vd; } void store_critical_edge(typename inhr_t::pi_edge_t ed, typename inhr_t::critical_cycle_t& cc) { typename inhr_t::pi_vertex_t s = boost::source(ed, this->m_pi_g); typename inhr_t::pi_vertex_t t = boost::target(ed, this->m_pi_g); assert(boost::edge(m_pi_g2g_vm[s], m_pi_g2g_vm[t], this->m_g).second); cc.push_back(boost::edge(m_pi_g2g_vm[s], m_pi_g2g_vm[t], this->m_g).first); ///Store corresponding edge of the m_g } void store_critical_cycle(typename inhr_t::critical_cycle_t& cc) { m_critical_cycle.swap(cc); } private: typename inhr_t::critical_cycle_t m_critical_cycle; typedef typename boost::iterator_property_map<typename std::vector<typename inhr_t::mcr_vertex_t>::iterator, typename inhr_t::TPiGraphVertexIndexMap> pi_g2g_vm_t; pi_g2g_vm_t m_pi_g2g_vm; ///Maps policy graph vertices to input graph vertices typename std::vector<typename inhr_t::mcr_vertex_t> m_pi_g2g; }; /*! * Add sink vertex - this will make any graph good, the selfloop will have ratio equal to infinity * Properties must be "self increasing" */ template <typename TGraph, typename TWeight1EdgeMap, typename TWeight2EdgeMap> typename boost::graph_traits<TGraph>::vertex_descriptor make_graph_good(TGraph& g, TWeight1EdgeMap ewm, TWeight2EdgeMap ew2m, typename boost::property_traits<TWeight1EdgeMap>::value_type infinity) { typedef typename boost::graph_traits<TGraph>::edge_descriptor Edge; typename boost::graph_traits<TGraph>::vertex_descriptor sink = boost::add_vertex(g); BGL_FORALL_VERTICES_T(vd, g, TGraph) { Edge newed = boost::add_edge(vd, sink, g).first; boost::put(ewm, newed, 0); boost::put(ew2m, newed, 1); } Edge selfed = boost::edge(sink, sink, g).first; boost::put(ewm, selfed, infinity); return sink; } /*! * Construct from input graph g "safe" (suitable for maximum_cycle_ratio1() call) version - safeg */ template <typename TG, typename TIndVertexMap, typename TW1EdgeMap, typename TW2EdgeMap, typename TSafeG, typename SafeG2GEdgeMap> void construct_safe_graph(const TG& g, TIndVertexMap vim, TW1EdgeMap ew1m, TW2EdgeMap ew2m, TSafeG& safeg, SafeG2GEdgeMap& sg2gm) { assert(num_vertices(g) == num_vertices(safeg)); typedef typename graph_traits<TSafeG>::edge_descriptor tmp_edge_t; typedef typename graph_traits<TG>::edge_descriptor edge_t; typename graph_traits<TG>::edge_iterator ei, ei_end; for (tie(ei, ei_end) = edges(g); ei != ei_end; ++ei) { tmp_edge_t tmped = add_edge(vim[source(*ei, g)], vim[target(*ei, g)], safeg).first; sg2gm[tmped] = *ei; put(edge_weight, safeg, tmped, get(ew1m, *ei)); put(edge_weight2, safeg, tmped, get(ew2m, *ei)); } } template <typename TGraph, typename TVertexIndexMap, typename TWeight1EdgeMap, typename TWeight2EdgeMap> double maximum_cycle_ratio_good_graph(const TGraph& g, TVertexIndexMap vim, TWeight1EdgeMap ewm, TWeight2EdgeMap ew2m, typename std::vector<typename boost::graph_traits<TGraph>::edge_descriptor>* pcc = 0) { if (pcc == 0) { return detail::Cmcr_Howard<TGraph, TVertexIndexMap, TWeight1EdgeMap, TWeight2EdgeMap>(g, vim, ewm, ew2m)(); } else { detail::Cmcr_Howard1<TGraph, TVertexIndexMap, TWeight1EdgeMap, TWeight2EdgeMap> obj(g, vim, ewm, ew2m); double maxcr = obj(); obj.get_critical_cycle(*pcc); return maxcr; } } template <typename TGraph, typename TVertexIndexMap, typename TWeight1EdgeMap, typename TWeight2EdgeMap, typename TEdgeIndexMap> double minimum_cycle_ratio_good_graph(const TGraph& g, TVertexIndexMap vim, TWeight1EdgeMap ewm, TWeight2EdgeMap ew2m, TEdgeIndexMap eim, typename std::vector<typename boost::graph_traits<TGraph>::edge_descriptor>* pcc = 0) { typedef typename boost::remove_const<typename boost::property_traits<TWeight1EdgeMap>::value_type>::type weight_value_t; BOOST_STATIC_ASSERT(!is_integral<weight_value_t>::value || is_signed<weight_value_t>::value); typename std::vector<weight_value_t> ne_w(boost::num_edges(g)); BGL_FORALL_EDGES_T(ed, g, TGraph) ne_w[boost::get(eim, ed)] = -ewm[ed]; return -maximum_cycle_ratio_good_graph(g, vim, boost::make_iterator_property_map(ne_w.begin(), eim), ew2m, pcc); } /*! * \param g directed multigraph. * \param pcc - pointer to the critical edges list. * \param minus_infinity must be small enough to garanty that g has at least one cycle with greater ratio. * \return minus_infinity if there're no cycles in the graph */ template <typename TGraph, typename TWeight1EdgeMap, typename TWeight2EdgeMap> double maximum_cycle_ratio1(const TGraph& g, TWeight1EdgeMap ewm, TWeight2EdgeMap ew2m, typename std::vector<typename boost::graph_traits<TGraph>::edge_descriptor>* pcc = 0, typename boost::property_traits<TWeight1EdgeMap>::value_type minus_infinity = -(std::numeric_limits<int>::max)()) { typedef typename boost::graph_traits<TGraph>::vertex_descriptor Vertex; typedef typename boost::graph_traits<TGraph>::edge_descriptor Edge; boost::function_requires< boost::ReadWritePropertyMapConcept<TWeight1EdgeMap, Edge> >(); boost::function_requires< boost::ReadWritePropertyMapConcept<TWeight2EdgeMap, Edge> >(); TGraph& ncg = const_cast<TGraph&>(g); Vertex sink = detail::make_graph_good(ncg, ewm, ew2m, minus_infinity ); double res = maximum_cycle_ratio_good_graph(ncg, boost::get(boost::vertex_index, g), ewm, ew2m, pcc); boost::clear_vertex(sink, ncg); boost::remove_vertex(sink, ncg); return res; } /*! * Edge index MUST be in diapazon [0,..., num_edges(g)-1] * \return plus_infinity if g has no cycles. */ template <typename TGraph, typename TWeight1EdgeMap, typename TWeight2EdgeMap, typename TEdgeIndexMap> double minimum_cycle_ratio1(const TGraph& g, TWeight1EdgeMap ewm, TWeight2EdgeMap ew2m, TEdgeIndexMap eim, typename std::vector<typename boost::graph_traits<TGraph>::edge_descriptor>* pcc = 0, typename boost::property_traits<TWeight1EdgeMap>::value_type plus_infinity = (std::numeric_limits<int>::max)() ) { typedef typename boost::property_traits<TEdgeIndexMap>::value_type ei_t; typedef typename boost::graph_traits<TGraph>::vertex_descriptor Vertex; typedef typename boost::graph_traits<TGraph>::edge_descriptor Edge; boost::function_requires< boost::ReadWritePropertyMapConcept<TWeight1EdgeMap, Edge> >(); boost::function_requires< boost::ReadWritePropertyMapConcept<TWeight2EdgeMap, Edge> >(); boost::function_requires< boost::ReadWritePropertyMapConcept<TEdgeIndexMap, Edge> >(); TGraph& ncg = const_cast<TGraph&>(g); ei_t nei = ei_t(boost::num_edges(g)); Vertex sink = detail::make_graph_good(ncg, ewm, ew2m, plus_infinity ); ///Maintain edge index invariant BGL_FORALL_VERTICES_T(vd, ncg, TGraph) { typename boost::graph_traits<TGraph>::edge_descriptor ed = boost::edge(vd, sink, ncg).first; boost::put(eim, ed, nei++); } double res = minimum_cycle_ratio_good_graph(ncg, boost::get(boost::vertex_index, ncg), ewm, ew2m, eim, pcc); boost::clear_vertex(sink, ncg); boost::remove_vertex(sink, ncg); return res; } struct edge_less_than { template <typename TEdgeDescriptor> bool operator()(const TEdgeDescriptor& x, const TEdgeDescriptor& y) const { return x.get_property() < y.get_property(); } }; }///namespace detail namespace { template <typename TW1, typename TW2> struct safe_graph { typedef typename boost::adjacency_list<boost::vecS, boost::vecS, boost::directedS, boost::no_property, typename boost::property<boost::edge_weight_t, TW1, typename boost::property<boost::edge_weight2_t, TW2> > > type; }; } /*! * Calculate the maximum cycle ratio (mcr) of the directed multigraph g. * \param g directed multigraph * \param pcc - If provided then a critical cycle will be written to corresponding vector. * \param minus_infinity small enough value to garanty that g has at least one cycle with greater ratio. * \return mcr or minus_infinity if g has no cycles. */ template <typename TGraph, typename TVertexIndexMap, typename TW1EdgeMap, typename TW2EdgeMap> double maximum_cycle_ratio(const TGraph& g, TVertexIndexMap vim, TW1EdgeMap ew1m, TW2EdgeMap ew2m, typename std::vector<typename boost::graph_traits<TGraph>::edge_descriptor>* pcc = 0, typename boost::property_traits<TW1EdgeMap>::value_type minus_infinity = -(std::numeric_limits<int>::max)()) { typedef typename remove_const<typename property_traits<TW1EdgeMap>::value_type>::type w1_t; typedef typename remove_const<typename property_traits<TW2EdgeMap>::value_type>::type w2_t; typedef typename safe_graph<w1_t, w2_t>::type safe_graph_t; typedef typename graph_traits<safe_graph_t>::edge_descriptor tmp_edge_t; typedef typename graph_traits<TGraph>::edge_descriptor edge_t; typename std::map<tmp_edge_t, edge_t, detail::edge_less_than> tmpg2g; std::vector<tmp_edge_t> cc; safe_graph_t sg(num_vertices(g)); detail::construct_safe_graph(g, vim, ew1m, ew2m, sg, tmpg2g); double mcr = maximum_cycle_ratio1(sg, get(edge_weight, sg), get(edge_weight2, sg), pcc ? &cc : 0, minus_infinity); if (pcc && (mcr > minus_infinity)) { pcc->clear(); for (typename std::vector<tmp_edge_t>::iterator it = cc.begin(); it != cc.end(); ++it) pcc->push_back(tmpg2g[*it]); } return mcr; } template <typename TGraph, typename TVertexIndexMap, typename TW1EdgeMap, typename TW2EdgeMap, typename TIndEdgeMap> double minimum_cycle_ratio(const TGraph& g, TVertexIndexMap vim, TW1EdgeMap ew1m, TW2EdgeMap ew2m, TIndEdgeMap eim, typename std::vector<typename boost::graph_traits<TGraph>::edge_descriptor>* pcc = 0, typename boost::property_traits<TW1EdgeMap>::value_type plus_infinity = (std::numeric_limits<int>::max)()) { typedef typename boost::remove_const<typename boost::property_traits<TW1EdgeMap>::value_type>::type weight_value_t; BOOST_STATIC_ASSERT(!is_integral<weight_value_t>::value || is_signed<weight_value_t>::value); typename std::vector<weight_value_t> ne_w(boost::num_edges(g)); BGL_FORALL_EDGES_T(ed, g, TGraph) ne_w[boost::get(eim, ed)] = -ew1m[ed]; return -maximum_cycle_ratio(g, vim, boost::make_iterator_property_map(ne_w.begin(), eim), ew2m, pcc, -plus_infinity); } /*! * Calculate maximum mean cycle of directed weighted multigraph. * \param g directed multigraph * \return maximum mean cycle of g or minus_infinity if g has no cycles. */ template <typename TGraph, typename TVertexIndexMap, typename TWeightEdgeMap, typename TIndEdgeMap> double maximum_mean_cycle(const TGraph& g, TVertexIndexMap vim, TWeightEdgeMap ewm, TIndEdgeMap eim, typename std::vector<typename boost::graph_traits<TGraph>::edge_descriptor>* pcc = 0, typename boost::property_traits<TWeightEdgeMap>::value_type minus_infinity = -(std::numeric_limits<int>::max)()) { typedef typename boost::remove_const<typename boost::property_traits<TWeightEdgeMap>::value_type>::type weight_value_t; typedef typename boost::graph_traits<TGraph>::edge_descriptor Edge; typename std::vector<weight_value_t> ed_w2(boost::num_edges(g), 1); return maximum_cycle_ratio(g, vim, ewm, boost::make_iterator_property_map(ed_w2.begin(), eim), pcc, minus_infinity); } template <typename TGraph, typename TVertexIndexMap, typename TWeightEdgeMap, typename TIndEdgeMap> double minimum_mean_cycle(const TGraph& g, TVertexIndexMap vim, TWeightEdgeMap ewm, TIndEdgeMap eim, typename std::vector<typename boost::graph_traits<TGraph>::edge_descriptor>* pcc = 0, typename boost::property_traits<TWeightEdgeMap>::value_type plus_infinity = (std::numeric_limits<int>::max)()) { typedef typename boost::remove_const<typename boost::property_traits<TWeightEdgeMap>::value_type>::type weight_value_t; typedef typename boost::graph_traits<TGraph>::edge_descriptor Edge; typename std::vector<weight_value_t> ed_w2(boost::num_edges(g), 1); return minimum_cycle_ratio(g, vim, ewm, boost::make_iterator_property_map(ed_w2.begin(), eim), eim, pcc, plus_infinity); } } //namespace boost #endif