forked from jonstewart/boost-svn
-
Notifications
You must be signed in to change notification settings - Fork 0
/
howard_cycle_ratio.hpp
634 lines (575 loc) · 21.5 KB
/
howard_cycle_ratio.hpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
// 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 (tie(vi, vie) = vertices(m_g); vi != vie; ++vi)
{
for (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 ( tie(vi, vie) = vertices(m_g); vi != vie; ++vi )
{
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;
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 (tie(vi, vie) = vertices(m_g); vi != vie; ++vi)
{
if (!m_badv[*vi])
{
for (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