Mercurial > hg > vamp-build-and-test
diff DEPENDENCIES/generic/include/boost/graph/howard_cycle_ratio.hpp @ 16:2665513ce2d3
Add boost headers
author | Chris Cannam |
---|---|
date | Tue, 05 Aug 2014 11:11:38 +0100 |
parents | |
children |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DEPENDENCIES/generic/include/boost/graph/howard_cycle_ratio.hpp Tue Aug 05 11:11:38 2014 +0100 @@ -0,0 +1,635 @@ +// Copyright (C) 2006-2009 Dmitry Bufistov and Andrey Parfenov + +// 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) + +#ifndef BOOST_GRAPH_CYCLE_RATIO_HOWARD_HPP +#define BOOST_GRAPH_CYCLE_RATIO_HOWARD_HPP + +#include <vector> +#include <list> +#include <algorithm> +#include <limits> + +#include <boost/bind.hpp> +#include <boost/type_traits/is_same.hpp> +#include <boost/type_traits/remove_const.hpp> +#include <boost/concept_check.hpp> +#include <boost/pending/queue.hpp> +#include <boost/property_map/property_map.hpp> +#include <boost/graph/graph_traits.hpp> +#include <boost/graph/graph_concepts.hpp> +#include <boost/concept/assert.hpp> + +/** @file howard_cycle_ratio.hpp + * @brief The implementation of the maximum/minimum cycle ratio/mean algorithm. + * @author Dmitry Bufistov + * @author Andrey Parfenov + */ + +namespace boost { + + /** + * The mcr_float is like numeric_limits, but only for floating point types + * and only defines infinity() and epsilon(). This class is primarily used + * to encapsulate a less-precise epsilon than natively supported by the + * floating point type. + */ + template <typename Float = double> struct mcr_float { + typedef Float value_type; + + static Float infinity() + { return std::numeric_limits<value_type>::infinity(); } + + static Float epsilon() + { return Float(-0.005); } + }; + + namespace detail { + + template <typename FloatTraits> struct + min_comparator_props { + typedef std::greater<typename FloatTraits::value_type> comparator; + static const int multiplier = 1; + }; + + template <typename FloatTraits> struct + max_comparator_props { + typedef std::less<typename FloatTraits::value_type> comparator; + static const int multiplier = -1; + }; + + template <typename FloatTraits, typename ComparatorProps> + struct float_wrapper { + typedef typename FloatTraits::value_type value_type; + typedef ComparatorProps comparator_props_t; + typedef typename ComparatorProps::comparator comparator; + + static value_type infinity() + { return FloatTraits::infinity() * ComparatorProps::multiplier; } + + static value_type epsilon() + { return FloatTraits::epsilon() * ComparatorProps::multiplier; } + + }; + + /*! @class mcr_howard + * @brief Calculates optimum (maximum/minimum) cycle ratio of a directed graph. + * Uses Howard's iteration policy algorithm. </br>(It is described in the paper + * "Experimental Analysis of the Fastest Optimum Cycle Ratio and Mean Algorithm" + * by Ali Dasdan). + */ + template <typename FloatTraits, + typename Graph, typename VertexIndexMap, + typename EdgeWeight1, typename EdgeWeight2> + class mcr_howard + { + public: + typedef typename FloatTraits::value_type float_t; + typedef typename FloatTraits::comparator_props_t cmp_props_t; + typedef typename FloatTraits::comparator comparator_t; + typedef enum{ my_white = 0, my_black } my_color_type; + typedef typename graph_traits<Graph>::vertex_descriptor vertex_t; + typedef typename graph_traits<Graph>::edge_descriptor edge_t; + typedef typename graph_traits<Graph>::vertices_size_type vn_t; + typedef std::vector<float_t> vp_t; + typedef typename boost::iterator_property_map< + typename vp_t::iterator, VertexIndexMap + > distance_map_t; //V -> float_t + + typedef typename std::vector<edge_t> ve_t; + typedef std::vector<my_color_type> vcol_t; + typedef typename ::boost::iterator_property_map< + typename ve_t::iterator, VertexIndexMap + > policy_t; //Vertex -> Edge + typedef typename ::boost::iterator_property_map< + typename vcol_t::iterator, VertexIndexMap + > color_map_t; + + typedef typename std::list<vertex_t> pinel_t;// The in_edges list of the policy graph + typedef typename std::vector<pinel_t> inedges1_t; + typedef typename ::boost::iterator_property_map< + typename inedges1_t::iterator, VertexIndexMap + > inedges_t; + typedef typename std::vector<edge_t> critical_cycle_t; + + //Bad vertex flag. If true, then the vertex is "bad". + // Vertex is "bad" if its out_degree is equal to zero. + typedef typename boost::iterator_property_map< + std::vector<int>::iterator, VertexIndexMap + > badv_t; + + /*! + * Constructor + * \param g = (V, E) - a directed multigraph. + * \param vim Vertex Index Map. Read property Map: V -> [0, num_vertices(g)). + * \param ewm edge weight map. Read property map: E -> R + * \param ew2m edge weight map. Read property map: E -> R+ + * \param infty A big enough value to guaranty that there exist a cycle with + * better ratio. + * \param cmp The compare operator for float_ts. + */ + mcr_howard(const Graph &g, VertexIndexMap vim, + EdgeWeight1 ewm, EdgeWeight2 ew2m) : + m_g(g), m_vim(vim), m_ew1m(ewm), m_ew2m(ew2m), + m_bound(mcr_bound()), + m_cr(m_bound), + m_V(num_vertices(m_g)), + m_dis(m_V, 0), m_dm(m_dis.begin(), m_vim), + m_policyc(m_V), m_policy(m_policyc.begin(), m_vim), + m_inelc(m_V), m_inel(m_inelc.begin(), m_vim), + m_badvc(m_V, false), m_badv(m_badvc.begin(), m_vim), + m_colcv(m_V), + m_col_bfs(m_V) + { } + + /*! + * \return maximum/minimum_{for all cycles C} + * [sum_{e in C} w1(e)] / [sum_{e in C} w2(e)], + * or FloatTraits::infinity() if graph has no cycles. + */ + float_t ocr_howard() + { + construct_policy_graph(); + int k = 0; + float_t mcr = 0; + do + { + mcr = policy_mcr(); + ++k; + } + while (try_improve_policy(mcr) && k < 100); //To avoid infinite loop + + const float_t eps_ = -0.00000001 * cmp_props_t::multiplier; + if (m_cmp(mcr, m_bound + eps_)) + { + return FloatTraits::infinity(); + } + else + { + return mcr; + } + } + virtual ~mcr_howard() {} + + protected: + virtual void store_critical_edge(edge_t, critical_cycle_t &) {} + virtual void store_critical_cycle(critical_cycle_t &) {} + + private: + /*! + * \return lower/upper bound for the maximal/minimal cycle ratio + */ + float_t mcr_bound() + { + typename graph_traits<Graph>::vertex_iterator vi, vie; + typename graph_traits<Graph>::out_edge_iterator oei, oeie; + float_t cz = (std::numeric_limits<float_t>::max)(); //Closest to zero value + float_t s = 0; + const float_t eps_ = std::numeric_limits<float_t>::epsilon(); + for (boost::tie(vi, vie) = vertices(m_g); vi != vie; ++vi) + { + for (boost::tie(oei, oeie) = out_edges(*vi, m_g); oei != oeie; ++oei) + { + s += std::abs(m_ew1m[*oei]); + float_t a = std::abs(m_ew2m[*oei]); + if ( a > eps_ && a < cz) + { + cz = a; + } + } + } + return cmp_props_t::multiplier * (s / cz); + } + + + /*! + * Constructs an arbitrary policy graph. + */ + void construct_policy_graph() + { + m_sink = graph_traits<Graph>().null_vertex(); + typename graph_traits<Graph>::vertex_iterator vi, vie; + typename graph_traits<Graph>::out_edge_iterator oei, oeie; + for ( boost::tie(vi, vie) = vertices(m_g); vi != vie; ++vi ) + { + boost::tie(oei, oeie) = out_edges(*vi, m_g); + typename graph_traits<Graph>::out_edge_iterator mei = + std::max_element(oei, oeie, + boost::bind(m_cmp, + boost::bind(&EdgeWeight1::operator[], m_ew1m, _1), + boost::bind(&EdgeWeight1::operator[], m_ew1m, _2) + ) + ); + if (mei == oeie) + { + if (m_sink == graph_traits<Graph>().null_vertex()) + { + m_sink = *vi; + } + m_badv[*vi] = true; + m_inel[m_sink].push_back(*vi); + } + else + { + m_inel[target(*mei, m_g)].push_back(*vi); + m_policy[*vi] = *mei; + } + } + } + /*! Sets the distance value for all vertices "v" such that there is + * a path from "v" to "sv". It does "inverse" breadth first visit of the policy + * graph, starting from the vertex "sv". + */ + void mcr_bfv(vertex_t sv, float_t cr, color_map_t c) + { + boost::queue<vertex_t> Q; + c[sv] = my_black; + Q.push(sv); + while (!Q.empty()) + { + vertex_t v = Q.top(); Q.pop(); + for (typename pinel_t::const_iterator itr = m_inel[v].begin(); + itr != m_inel[v].end(); ++itr) + //For all in_edges of the policy graph + { + if (*itr != sv) + { + if (m_badv[*itr]) + { + m_dm[*itr] = m_dm[v] + m_bound - cr; + } + else + { + m_dm[*itr] = m_dm[v] + m_ew1m[m_policy[*itr]] - + m_ew2m[m_policy[*itr]] * cr; + } + c[*itr] = my_black; + Q.push(*itr); + } + } + } + } + + /*! + * \param sv an arbitrary (undiscovered) vertex of the policy graph. + * \return a vertex in the policy graph that belongs to a cycle. + * Performs a depth first visit until a cycle edge is found. + */ + vertex_t find_cycle_vertex(vertex_t sv) + { + vertex_t gv = sv; + std::fill(m_colcv.begin(), m_colcv.end(), my_white); + color_map_t cm(m_colcv.begin(), m_vim); + do + { + cm[gv] = my_black; + if (! m_badv[gv]) + { + gv = target(m_policy[gv], m_g); + } + else + { + gv = m_sink; + } + } + while (cm[gv] != my_black); + return gv; + } + + /*! + * \param sv - vertex that belongs to a cycle in the policy graph. + */ + float_t cycle_ratio(vertex_t sv) + { + if (sv == m_sink) return m_bound; + std::pair<float_t, float_t> sums_(float_t(0), float_t(0)); + vertex_t v = sv; + critical_cycle_t cc; + do + { + store_critical_edge(m_policy[v], cc); + sums_.first += m_ew1m[m_policy[v]]; + sums_.second += m_ew2m[m_policy[v]]; + v = target(m_policy[v], m_g); + } + while (v != sv); + float_t cr = sums_.first / sums_.second; + if ( m_cmp(m_cr, cr) ) + { + m_cr = cr; + store_critical_cycle(cc); + } + return cr; + } + + /*! + * Finds the optimal cycle ratio of the policy graph + */ + float_t policy_mcr() + { + std::fill(m_col_bfs.begin(), m_col_bfs.end(), my_white); + color_map_t vcm_ = color_map_t(m_col_bfs.begin(), m_vim); + typename graph_traits<Graph>::vertex_iterator uv_itr, vie; + boost::tie(uv_itr, vie) = vertices(m_g); + float_t mcr = m_bound; + while ( (uv_itr = std::find_if(uv_itr, vie, + boost::bind(std::equal_to<my_color_type>(), + my_white, + boost::bind(&color_map_t::operator[], vcm_, _1) + ) + ) + ) != vie ) + ///While there are undiscovered vertices + { + vertex_t gv = find_cycle_vertex(*uv_itr); + float_t cr = cycle_ratio(gv) ; + mcr_bfv(gv, cr, vcm_); + if ( m_cmp(mcr, cr) ) mcr = cr; + ++uv_itr; + } + return mcr; + } + + /*! + * Changes the edge m_policy[s] to the new_edge. + */ + void improve_policy(vertex_t s, edge_t new_edge) + { + vertex_t t = target(m_policy[s], m_g); + typename property_traits<VertexIndexMap>::value_type ti = m_vim[t]; + m_inelc[ti].erase( std::find(m_inelc[ti].begin(), m_inelc[ti].end(), s)); + m_policy[s] = new_edge; + t = target(new_edge, m_g); + m_inel[t].push_back(s); ///Maintain in_edge list + } + + /*! + * A negative cycle detector. + */ + bool try_improve_policy(float_t cr) + { + bool improved = false; + typename graph_traits<Graph>::vertex_iterator vi, vie; + typename graph_traits<Graph>::out_edge_iterator oei, oeie; + const float_t eps_ = FloatTraits::epsilon(); + for (boost::tie(vi, vie) = vertices(m_g); vi != vie; ++vi) + { + if (!m_badv[*vi]) + { + for (boost::tie(oei, oeie) = out_edges(*vi, m_g); oei != oeie; ++oei) + { + vertex_t t = target(*oei, m_g); + //Current distance from *vi to some vertex + float_t dis_ = m_ew1m[*oei] - m_ew2m[*oei] * cr + m_dm[t]; + if ( m_cmp(m_dm[*vi] + eps_, dis_) ) + { + improve_policy(*vi, *oei); + m_dm[*vi] = dis_; + improved = true; + } + } + } + else + { + float_t dis_ = m_bound - cr + m_dm[m_sink]; + if ( m_cmp(m_dm[*vi] + eps_, dis_) ) + { + m_dm[*vi] = dis_; + } + } + } + return improved; + } + private: + const Graph &m_g; + VertexIndexMap m_vim; + EdgeWeight1 m_ew1m; + EdgeWeight2 m_ew2m; + comparator_t m_cmp; + float_t m_bound; //> The lower/upper bound to the maximal/minimal cycle ratio + float_t m_cr; //>The best cycle ratio that has been found so far + + vn_t m_V; //>The number of the vertices in the graph + vp_t m_dis; //>Container for the distance map + distance_map_t m_dm; //>Distance map + + ve_t m_policyc; //>Container for the policy graph + policy_t m_policy; //>The interface for the policy graph + + inedges1_t m_inelc; //>Container fot in edges list + inedges_t m_inel; //>Policy graph, input edges list + + std::vector<int> m_badvc; + badv_t m_badv; //Marks "bad" vertices + + vcol_t m_colcv, m_col_bfs; //Color maps + vertex_t m_sink; //To convert any graph to "good" + }; + + /*! \class mcr_howard1 + * \brief Finds optimum cycle raio and a critical cycle + */ + template <typename FloatTraits, + typename Graph, typename VertexIndexMap, + typename EdgeWeight1, typename EdgeWeight2> + class mcr_howard1 : public + mcr_howard<FloatTraits, Graph, VertexIndexMap, + EdgeWeight1, EdgeWeight2> + { + public: + typedef mcr_howard<FloatTraits, Graph, VertexIndexMap, + EdgeWeight1, EdgeWeight2> inhr_t; + mcr_howard1(const Graph &g, VertexIndexMap vim, + EdgeWeight1 ewm, EdgeWeight2 ew2m) : + inhr_t(g, vim, ewm, ew2m) + { } + + void get_critical_cycle(typename inhr_t::critical_cycle_t &cc) + { return cc.swap(m_cc); } + + protected: + void store_critical_edge(typename inhr_t::edge_t ed, + typename inhr_t::critical_cycle_t &cc) + { cc.push_back(ed); } + + void store_critical_cycle(typename inhr_t::critical_cycle_t &cc) + { m_cc.swap(cc); } + + private: + typename inhr_t::critical_cycle_t m_cc; //Critical cycle + }; + + /*! + * \param g a directed multigraph. + * \param vim Vertex Index Map. A map V->[0, num_vertices(g)) + * \param ewm Edge weight1 map. + * \param ew2m Edge weight2 map. + * \param pcc pointer to the critical edges list. + * \return Optimum cycle ratio of g or FloatTraits::infinity() if g has no cycles. + */ + template <typename FT, + typename TG, typename TVIM, + typename TEW1, typename TEW2, + typename EV> + typename FT::value_type + optimum_cycle_ratio(const TG &g, TVIM vim, TEW1 ewm, TEW2 ew2m, EV* pcc) + { + typedef typename graph_traits<TG>::directed_category DirCat; + BOOST_STATIC_ASSERT((is_convertible<DirCat*, directed_tag*>::value == true)); + BOOST_CONCEPT_ASSERT(( IncidenceGraphConcept<TG> )); + BOOST_CONCEPT_ASSERT(( VertexListGraphConcept<TG> )); + typedef typename graph_traits<TG>::vertex_descriptor Vertex; + BOOST_CONCEPT_ASSERT(( ReadablePropertyMapConcept<TVIM, Vertex> )); + typedef typename graph_traits<TG>::edge_descriptor Edge; + BOOST_CONCEPT_ASSERT(( ReadablePropertyMapConcept<TEW1, Edge> )); + BOOST_CONCEPT_ASSERT(( ReadablePropertyMapConcept<TEW2, Edge> )); + + if(pcc == 0) { + return detail::mcr_howard<FT,TG, TVIM, TEW1, TEW2>( + g, vim, ewm, ew2m + ).ocr_howard(); + } + + detail::mcr_howard1<FT, TG, TVIM, TEW1, TEW2> obj(g, vim, ewm, ew2m); + double ocr = obj.ocr_howard(); + obj.get_critical_cycle(*pcc); + return ocr; + } + } // namespace detail + +// Algorithms +// Maximum Cycle Ratio + +template < + typename FloatTraits, + typename Graph, + typename VertexIndexMap, + typename EdgeWeight1Map, + typename EdgeWeight2Map> +inline typename FloatTraits::value_type +maximum_cycle_ratio(const Graph &g, VertexIndexMap vim, EdgeWeight1Map ew1m, + EdgeWeight2Map ew2m, + std::vector<typename graph_traits<Graph>::edge_descriptor>* pcc = 0, + FloatTraits = FloatTraits()) +{ + typedef detail::float_wrapper< + FloatTraits, detail::max_comparator_props<FloatTraits> + > Traits; + return detail::optimum_cycle_ratio<Traits>(g, vim, ew1m, ew2m, pcc); +} + +template < + typename Graph, + typename VertexIndexMap, + typename EdgeWeight1Map, + typename EdgeWeight2Map> +inline double +maximum_cycle_ratio(const Graph &g, VertexIndexMap vim, + EdgeWeight1Map ew1m, EdgeWeight2Map ew2m, + std::vector<typename graph_traits<Graph>::edge_descriptor>* pcc = 0) +{ return maximum_cycle_ratio(g, vim, ew1m, ew2m, pcc, mcr_float<>()); } + +// Minimum Cycle Ratio + +template < + typename FloatTraits, + typename Graph, + typename VertexIndexMap, + typename EdgeWeight1Map, + typename EdgeWeight2Map> +typename FloatTraits::value_type +minimum_cycle_ratio(const Graph &g, VertexIndexMap vim, + EdgeWeight1Map ew1m, EdgeWeight2Map ew2m, + std::vector<typename graph_traits<Graph>::edge_descriptor> *pcc = 0, + FloatTraits = FloatTraits()) +{ + typedef detail::float_wrapper< + FloatTraits, detail::min_comparator_props<FloatTraits> + > Traits; + return detail::optimum_cycle_ratio<Traits>(g, vim, ew1m, ew2m, pcc); +} + +template < + typename Graph, + typename VertexIndexMap, + typename EdgeWeight1Map, + typename EdgeWeight2Map> +inline double +minimum_cycle_ratio(const Graph &g, VertexIndexMap vim, + EdgeWeight1Map ew1m, EdgeWeight2Map ew2m, + std::vector<typename graph_traits<Graph>::edge_descriptor>* pcc = 0) +{ return minimum_cycle_ratio(g, vim, ew1m, ew2m, pcc, mcr_float<>()); } + +// Maximum Cycle Mean + +template < + typename FloatTraits, + typename Graph, + typename VertexIndexMap, + typename EdgeWeightMap, + typename EdgeIndexMap> +inline typename FloatTraits::value_type +maximum_cycle_mean(const Graph &g, VertexIndexMap vim, + EdgeWeightMap ewm, EdgeIndexMap eim, + std::vector<typename graph_traits<Graph>::edge_descriptor>* pcc = 0, + FloatTraits ft = FloatTraits()) +{ + typedef typename remove_const< + typename property_traits<EdgeWeightMap>::value_type + >::type Weight; + typename std::vector<Weight> ed_w2(boost::num_edges(g), 1); + return maximum_cycle_ratio(g, vim, ewm, + make_iterator_property_map(ed_w2.begin(), eim), + pcc, ft); +} + +template < + typename Graph, + typename VertexIndexMap, + typename EdgeWeightMap, + typename EdgeIndexMap> +inline double +maximum_cycle_mean(const Graph& g, VertexIndexMap vim, + EdgeWeightMap ewm, EdgeIndexMap eim, + std::vector<typename graph_traits<Graph>::edge_descriptor>* pcc = 0) +{ return maximum_cycle_mean(g, vim, ewm, eim, pcc, mcr_float<>()); } + +// Minimum Cycle Mean + +template < + typename FloatTraits, + typename Graph, + typename VertexIndexMap, + typename EdgeWeightMap, + typename EdgeIndexMap> +inline typename FloatTraits::value_type +minimum_cycle_mean(const Graph &g, VertexIndexMap vim, + EdgeWeightMap ewm, EdgeIndexMap eim, + std::vector<typename graph_traits<Graph>::edge_descriptor>* pcc = 0, + FloatTraits ft = FloatTraits()) +{ + typedef typename remove_const< + typename property_traits<EdgeWeightMap>::value_type + >::type Weight; + typename std::vector<Weight> ed_w2(boost::num_edges(g), 1); + return minimum_cycle_ratio(g, vim, ewm, + make_iterator_property_map(ed_w2.begin(), eim), + pcc, ft); +} + +template < + typename Graph, + typename VertexIndexMap, + typename EdgeWeightMap, + typename EdgeIndexMap> +inline double +minimum_cycle_mean(const Graph &g, VertexIndexMap vim, + EdgeWeightMap ewm, EdgeIndexMap eim, + std::vector<typename graph_traits<Graph>::edge_descriptor>* pcc = 0) +{ return minimum_cycle_mean(g, vim, ewm, eim, pcc, mcr_float<>()); } + +} //namespace boost + +#endif