// 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> | |
/** @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)); | |
function_requires< IncidenceGraphConcept<TG> >(); | |
function_requires< VertexListGraphConcept<TG> >(); | |
typedef typename graph_traits<TG>::vertex_descriptor Vertex; | |
function_requires< ReadablePropertyMapConcept<TVIM, Vertex> >(); | |
typedef typename graph_traits<TG>::edge_descriptor Edge; | |
function_requires< ReadablePropertyMapConcept<TEW1, Edge> >(); | |
function_requires< 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 |