COIN-OR::LEMON - Graph Library

source: lemon-1.2/lemon/cancel_and_tighten.h @ 814:0643a9c2c3ae

Last change on this file since 814:0643a9c2c3ae was 814:0643a9c2c3ae, checked in by Peter Kovacs <kpeter@…>, 14 years ago

Port cycle canceling algorithms from SVN -r3524 (#180)

File size: 26.3 KB
Line 
1/* -*- C++ -*-
2 *
3 * This file is a part of LEMON, a generic C++ optimization library
4 *
5 * Copyright (C) 2003-2008
6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
7 * (Egervary Research Group on Combinatorial Optimization, EGRES).
8 *
9 * Permission to use, modify and distribute this software is granted
10 * provided that this copyright notice appears in all copies. For
11 * precise terms see the accompanying LICENSE file.
12 *
13 * This software is provided "AS IS" with no warranty of any kind,
14 * express or implied, and with no claim as to its suitability for any
15 * purpose.
16 *
17 */
18
19#ifndef LEMON_CANCEL_AND_TIGHTEN_H
20#define LEMON_CANCEL_AND_TIGHTEN_H
21
22/// \ingroup min_cost_flow
23///
24/// \file
25/// \brief Cancel and Tighten algorithm for finding a minimum cost flow.
26
27#include <vector>
28
29#include <lemon/circulation.h>
30#include <lemon/bellman_ford.h>
31#include <lemon/howard.h>
32#include <lemon/adaptors.h>
33#include <lemon/tolerance.h>
34#include <lemon/math.h>
35
36#include <lemon/static_graph.h>
37
38namespace lemon {
39
40  /// \addtogroup min_cost_flow
41  /// @{
42
43  /// \brief Implementation of the Cancel and Tighten algorithm for
44  /// finding a minimum cost flow.
45  ///
46  /// \ref CancelAndTighten implements the Cancel and Tighten algorithm for
47  /// finding a minimum cost flow.
48  ///
49  /// \tparam Digraph The digraph type the algorithm runs on.
50  /// \tparam LowerMap The type of the lower bound map.
51  /// \tparam CapacityMap The type of the capacity (upper bound) map.
52  /// \tparam CostMap The type of the cost (length) map.
53  /// \tparam SupplyMap The type of the supply map.
54  ///
55  /// \warning
56  /// - Arc capacities and costs should be \e non-negative \e integers.
57  /// - Supply values should be \e signed \e integers.
58  /// - The value types of the maps should be convertible to each other.
59  /// - \c CostMap::Value must be signed type.
60  ///
61  /// \author Peter Kovacs
62  template < typename Digraph,
63             typename LowerMap = typename Digraph::template ArcMap<int>,
64             typename CapacityMap = typename Digraph::template ArcMap<int>,
65             typename CostMap = typename Digraph::template ArcMap<int>,
66             typename SupplyMap = typename Digraph::template NodeMap<int> >
67  class CancelAndTighten
68  {
69    TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
70
71    typedef typename CapacityMap::Value Capacity;
72    typedef typename CostMap::Value Cost;
73    typedef typename SupplyMap::Value Supply;
74    typedef typename Digraph::template ArcMap<Capacity> CapacityArcMap;
75    typedef typename Digraph::template NodeMap<Supply> SupplyNodeMap;
76
77    typedef ResidualDigraph< const Digraph,
78      CapacityArcMap, CapacityArcMap > ResDigraph;
79
80  public:
81
82    /// The type of the flow map.
83    typedef typename Digraph::template ArcMap<Capacity> FlowMap;
84    /// The type of the potential map.
85    typedef typename Digraph::template NodeMap<Cost> PotentialMap;
86
87  private:
88
89    /// \brief Map adaptor class for handling residual arc costs.
90    ///
91    /// Map adaptor class for handling residual arc costs.
92    class ResidualCostMap : public MapBase<typename ResDigraph::Arc, Cost>
93    {
94      typedef typename ResDigraph::Arc Arc;
95     
96    private:
97
98      const CostMap &_cost_map;
99
100    public:
101
102      ///\e
103      ResidualCostMap(const CostMap &cost_map) : _cost_map(cost_map) {}
104
105      ///\e
106      Cost operator[](const Arc &e) const {
107        return ResDigraph::forward(e) ? _cost_map[e] : -_cost_map[e];
108      }
109
110    }; //class ResidualCostMap
111
112    /// \brief Map adaptor class for handling reduced arc costs.
113    ///
114    /// Map adaptor class for handling reduced arc costs.
115    class ReducedCostMap : public MapBase<Arc, Cost>
116    {
117    private:
118
119      const Digraph &_gr;
120      const CostMap &_cost_map;
121      const PotentialMap &_pot_map;
122
123    public:
124
125      ///\e
126      ReducedCostMap( const Digraph &gr,
127                      const CostMap &cost_map,
128                      const PotentialMap &pot_map ) :
129        _gr(gr), _cost_map(cost_map), _pot_map(pot_map) {}
130
131      ///\e
132      inline Cost operator[](const Arc &e) const {
133        return _cost_map[e] + _pot_map[_gr.source(e)]
134                            - _pot_map[_gr.target(e)];
135      }
136
137    }; //class ReducedCostMap
138
139    struct BFOperationTraits {
140      static double zero() { return 0; }
141
142      static double infinity() {
143        return std::numeric_limits<double>::infinity();
144      }
145
146      static double plus(const double& left, const double& right) {
147        return left + right;
148      }
149
150      static bool less(const double& left, const double& right) {
151        return left + 1e-6 < right;
152      }
153    }; // class BFOperationTraits
154
155  private:
156
157    // The digraph the algorithm runs on
158    const Digraph &_graph;
159    // The original lower bound map
160    const LowerMap *_lower;
161    // The modified capacity map
162    CapacityArcMap _capacity;
163    // The original cost map
164    const CostMap &_cost;
165    // The modified supply map
166    SupplyNodeMap _supply;
167    bool _valid_supply;
168
169    // Arc map of the current flow
170    FlowMap *_flow;
171    bool _local_flow;
172    // Node map of the current potentials
173    PotentialMap *_potential;
174    bool _local_potential;
175
176    // The residual digraph
177    ResDigraph *_res_graph;
178    // The residual cost map
179    ResidualCostMap _res_cost;
180
181  public:
182
183    /// \brief General constructor (with lower bounds).
184    ///
185    /// General constructor (with lower bounds).
186    ///
187    /// \param digraph The digraph the algorithm runs on.
188    /// \param lower The lower bounds of the arcs.
189    /// \param capacity The capacities (upper bounds) of the arcs.
190    /// \param cost The cost (length) values of the arcs.
191    /// \param supply The supply values of the nodes (signed).
192    CancelAndTighten( const Digraph &digraph,
193                      const LowerMap &lower,
194                      const CapacityMap &capacity,
195                      const CostMap &cost,
196                      const SupplyMap &supply ) :
197      _graph(digraph), _lower(&lower), _capacity(digraph), _cost(cost),
198      _supply(digraph), _flow(NULL), _local_flow(false),
199      _potential(NULL), _local_potential(false),
200      _res_graph(NULL), _res_cost(_cost)
201    {
202      // Check the sum of supply values
203      Supply sum = 0;
204      for (NodeIt n(_graph); n != INVALID; ++n) {
205        _supply[n] = supply[n];
206        sum += _supply[n];
207      }
208      _valid_supply = sum == 0;
209
210      // Remove non-zero lower bounds
211      for (ArcIt e(_graph); e != INVALID; ++e) {
212        _capacity[e] = capacity[e];
213        if (lower[e] != 0) {
214          _capacity[e] -= lower[e];
215          _supply[_graph.source(e)] -= lower[e];
216          _supply[_graph.target(e)] += lower[e];
217        }
218      }
219    }
220/*
221    /// \brief General constructor (without lower bounds).
222    ///
223    /// General constructor (without lower bounds).
224    ///
225    /// \param digraph The digraph the algorithm runs on.
226    /// \param capacity The capacities (upper bounds) of the arcs.
227    /// \param cost The cost (length) values of the arcs.
228    /// \param supply The supply values of the nodes (signed).
229    CancelAndTighten( const Digraph &digraph,
230                      const CapacityMap &capacity,
231                      const CostMap &cost,
232                      const SupplyMap &supply ) :
233      _graph(digraph), _lower(NULL), _capacity(capacity), _cost(cost),
234      _supply(supply), _flow(NULL), _local_flow(false),
235      _potential(NULL), _local_potential(false),
236      _res_graph(NULL), _res_cost(_cost)
237    {
238      // Check the sum of supply values
239      Supply sum = 0;
240      for (NodeIt n(_graph); n != INVALID; ++n) sum += _supply[n];
241      _valid_supply = sum == 0;
242    }
243
244    /// \brief Simple constructor (with lower bounds).
245    ///
246    /// Simple constructor (with lower bounds).
247    ///
248    /// \param digraph The digraph the algorithm runs on.
249    /// \param lower The lower bounds of the arcs.
250    /// \param capacity The capacities (upper bounds) of the arcs.
251    /// \param cost The cost (length) values of the arcs.
252    /// \param s The source node.
253    /// \param t The target node.
254    /// \param flow_value The required amount of flow from node \c s
255    /// to node \c t (i.e. the supply of \c s and the demand of \c t).
256    CancelAndTighten( const Digraph &digraph,
257                      const LowerMap &lower,
258                      const CapacityMap &capacity,
259                      const CostMap &cost,
260                      Node s, Node t,
261                      Supply flow_value ) :
262      _graph(digraph), _lower(&lower), _capacity(capacity), _cost(cost),
263      _supply(digraph, 0), _flow(NULL), _local_flow(false),
264      _potential(NULL), _local_potential(false),
265      _res_graph(NULL), _res_cost(_cost)
266    {
267      // Remove non-zero lower bounds
268      _supply[s] =  flow_value;
269      _supply[t] = -flow_value;
270      for (ArcIt e(_graph); e != INVALID; ++e) {
271        if (lower[e] != 0) {
272          _capacity[e] -= lower[e];
273          _supply[_graph.source(e)] -= lower[e];
274          _supply[_graph.target(e)] += lower[e];
275        }
276      }
277      _valid_supply = true;
278    }
279
280    /// \brief Simple constructor (without lower bounds).
281    ///
282    /// Simple constructor (without lower bounds).
283    ///
284    /// \param digraph The digraph the algorithm runs on.
285    /// \param capacity The capacities (upper bounds) of the arcs.
286    /// \param cost The cost (length) values of the arcs.
287    /// \param s The source node.
288    /// \param t The target node.
289    /// \param flow_value The required amount of flow from node \c s
290    /// to node \c t (i.e. the supply of \c s and the demand of \c t).
291    CancelAndTighten( const Digraph &digraph,
292                      const CapacityMap &capacity,
293                      const CostMap &cost,
294                      Node s, Node t,
295                      Supply flow_value ) :
296      _graph(digraph), _lower(NULL), _capacity(capacity), _cost(cost),
297      _supply(digraph, 0), _flow(NULL), _local_flow(false),
298      _potential(NULL), _local_potential(false),
299      _res_graph(NULL), _res_cost(_cost)
300    {
301      _supply[s] =  flow_value;
302      _supply[t] = -flow_value;
303      _valid_supply = true;
304    }
305*/
306    /// Destructor.
307    ~CancelAndTighten() {
308      if (_local_flow) delete _flow;
309      if (_local_potential) delete _potential;
310      delete _res_graph;
311    }
312
313    /// \brief Set the flow map.
314    ///
315    /// Set the flow map.
316    ///
317    /// \return \c (*this)
318    CancelAndTighten& flowMap(FlowMap &map) {
319      if (_local_flow) {
320        delete _flow;
321        _local_flow = false;
322      }
323      _flow = &map;
324      return *this;
325    }
326
327    /// \brief Set the potential map.
328    ///
329    /// Set the potential map.
330    ///
331    /// \return \c (*this)
332    CancelAndTighten& potentialMap(PotentialMap &map) {
333      if (_local_potential) {
334        delete _potential;
335        _local_potential = false;
336      }
337      _potential = &map;
338      return *this;
339    }
340
341    /// \name Execution control
342
343    /// @{
344
345    /// \brief Run the algorithm.
346    ///
347    /// Run the algorithm.
348    ///
349    /// \return \c true if a feasible flow can be found.
350    bool run() {
351      return init() && start();
352    }
353
354    /// @}
355
356    /// \name Query Functions
357    /// The result of the algorithm can be obtained using these
358    /// functions.\n
359    /// \ref lemon::CancelAndTighten::run() "run()" must be called before
360    /// using them.
361
362    /// @{
363
364    /// \brief Return a const reference to the arc map storing the
365    /// found flow.
366    ///
367    /// Return a const reference to the arc map storing the found flow.
368    ///
369    /// \pre \ref run() must be called before using this function.
370    const FlowMap& flowMap() const {
371      return *_flow;
372    }
373
374    /// \brief Return a const reference to the node map storing the
375    /// found potentials (the dual solution).
376    ///
377    /// Return a const reference to the node map storing the found
378    /// potentials (the dual solution).
379    ///
380    /// \pre \ref run() must be called before using this function.
381    const PotentialMap& potentialMap() const {
382      return *_potential;
383    }
384
385    /// \brief Return the flow on the given arc.
386    ///
387    /// Return the flow on the given arc.
388    ///
389    /// \pre \ref run() must be called before using this function.
390    Capacity flow(const Arc& arc) const {
391      return (*_flow)[arc];
392    }
393
394    /// \brief Return the potential of the given node.
395    ///
396    /// Return the potential of the given node.
397    ///
398    /// \pre \ref run() must be called before using this function.
399    Cost potential(const Node& node) const {
400      return (*_potential)[node];
401    }
402
403    /// \brief Return the total cost of the found flow.
404    ///
405    /// Return the total cost of the found flow. The complexity of the
406    /// function is \f$ O(e) \f$.
407    ///
408    /// \pre \ref run() must be called before using this function.
409    Cost totalCost() const {
410      Cost c = 0;
411      for (ArcIt e(_graph); e != INVALID; ++e)
412        c += (*_flow)[e] * _cost[e];
413      return c;
414    }
415
416    /// @}
417
418  private:
419
420    /// Initialize the algorithm.
421    bool init() {
422      if (!_valid_supply) return false;
423
424      // Initialize flow and potential maps
425      if (!_flow) {
426        _flow = new FlowMap(_graph);
427        _local_flow = true;
428      }
429      if (!_potential) {
430        _potential = new PotentialMap(_graph);
431        _local_potential = true;
432      }
433
434      _res_graph = new ResDigraph(_graph, _capacity, *_flow);
435
436      // Find a feasible flow using Circulation
437      Circulation< Digraph, ConstMap<Arc, Capacity>,
438                   CapacityArcMap, SupplyMap >
439        circulation( _graph, constMap<Arc>(Capacity(0)),
440                     _capacity, _supply );
441      return circulation.flowMap(*_flow).run();
442    }
443
444    bool start() {
445      const double LIMIT_FACTOR = 0.01;
446      const int MIN_LIMIT = 3;
447
448      typedef typename Digraph::template NodeMap<double> FloatPotentialMap;
449      typedef typename Digraph::template NodeMap<int> LevelMap;
450      typedef typename Digraph::template NodeMap<bool> BoolNodeMap;
451      typedef typename Digraph::template NodeMap<Node> PredNodeMap;
452      typedef typename Digraph::template NodeMap<Arc> PredArcMap;
453      typedef typename ResDigraph::template ArcMap<double> ResShiftCostMap;
454      FloatPotentialMap pi(_graph);
455      LevelMap level(_graph);
456      BoolNodeMap reached(_graph);
457      BoolNodeMap processed(_graph);
458      PredNodeMap pred_node(_graph);
459      PredArcMap pred_arc(_graph);
460      int node_num = countNodes(_graph);
461      typedef std::pair<Arc, bool> pair;
462      std::vector<pair> stack(node_num);
463      std::vector<Node> proc_vector(node_num);
464      ResShiftCostMap shift_cost(*_res_graph);
465
466      Tolerance<double> tol;
467      tol.epsilon(1e-6);
468
469      Timer t1, t2, t3;
470      t1.reset();
471      t2.reset();
472      t3.reset();
473
474      // Initialize epsilon and the node potentials
475      double epsilon = 0;
476      for (ArcIt e(_graph); e != INVALID; ++e) {
477        if (_capacity[e] - (*_flow)[e] > 0 && _cost[e] < -epsilon)
478          epsilon = -_cost[e];
479        else if ((*_flow)[e] > 0 && _cost[e] > epsilon)
480          epsilon = _cost[e];
481      }
482      for (NodeIt v(_graph); v != INVALID; ++v) {
483        pi[v] = 0;
484      }
485
486      // Start phases
487      int limit = int(LIMIT_FACTOR * node_num);
488      if (limit < MIN_LIMIT) limit = MIN_LIMIT;
489      int iter = limit;
490      while (epsilon * node_num >= 1) {
491        t1.start();
492        // Find and cancel cycles in the admissible digraph using DFS
493        for (NodeIt n(_graph); n != INVALID; ++n) {
494          reached[n] = false;
495          processed[n] = false;
496        }
497        int stack_head = -1;
498        int proc_head = -1;
499
500        for (NodeIt start(_graph); start != INVALID; ++start) {
501          if (reached[start]) continue;
502
503          // New start node
504          reached[start] = true;
505          pred_arc[start] = INVALID;
506          pred_node[start] = INVALID;
507
508          // Find the first admissible residual outgoing arc
509          double p = pi[start];
510          Arc e;
511          _graph.firstOut(e, start);
512          while ( e != INVALID && (_capacity[e] - (*_flow)[e] == 0 ||
513                  !tol.negative(_cost[e] + p - pi[_graph.target(e)])) )
514            _graph.nextOut(e);
515          if (e != INVALID) {
516            stack[++stack_head] = pair(e, true);
517            goto next_step_1;
518          }
519          _graph.firstIn(e, start);
520          while ( e != INVALID && ((*_flow)[e] == 0 ||
521                  !tol.negative(-_cost[e] + p - pi[_graph.source(e)])) )
522            _graph.nextIn(e);
523          if (e != INVALID) {
524            stack[++stack_head] = pair(e, false);
525            goto next_step_1;
526          }
527          processed[start] = true;
528          proc_vector[++proc_head] = start;
529          continue;
530        next_step_1:
531
532          while (stack_head >= 0) {
533            Arc se = stack[stack_head].first;
534            bool sf = stack[stack_head].second;
535            Node u, v;
536            if (sf) {
537              u = _graph.source(se);
538              v = _graph.target(se);
539            } else {
540              u = _graph.target(se);
541              v = _graph.source(se);
542            }
543
544            if (!reached[v]) {
545              // A new node is reached
546              reached[v] = true;
547              pred_node[v] = u;
548              pred_arc[v] = se;
549              // Find the first admissible residual outgoing arc
550              double p = pi[v];
551              Arc e;
552              _graph.firstOut(e, v);
553              while ( e != INVALID && (_capacity[e] - (*_flow)[e] == 0 ||
554                      !tol.negative(_cost[e] + p - pi[_graph.target(e)])) )
555                _graph.nextOut(e);
556              if (e != INVALID) {
557                stack[++stack_head] = pair(e, true);
558                goto next_step_2;
559              }
560              _graph.firstIn(e, v);
561              while ( e != INVALID && ((*_flow)[e] == 0 ||
562                      !tol.negative(-_cost[e] + p - pi[_graph.source(e)])) )
563                _graph.nextIn(e);
564              stack[++stack_head] = pair(e, false);
565            next_step_2: ;
566            } else {
567              if (!processed[v]) {
568                // A cycle is found
569                Node n, w = u;
570                Capacity d, delta = sf ? _capacity[se] - (*_flow)[se] :
571                                         (*_flow)[se];
572                for (n = u; n != v; n = pred_node[n]) {
573                  d = _graph.target(pred_arc[n]) == n ?
574                      _capacity[pred_arc[n]] - (*_flow)[pred_arc[n]] :
575                      (*_flow)[pred_arc[n]];
576                  if (d <= delta) {
577                    delta = d;
578                    w = pred_node[n];
579                  }
580                }
581
582/*
583                std::cout << "CYCLE FOUND: ";
584                if (sf)
585                  std::cout << _cost[se] + pi[_graph.source(se)] - pi[_graph.target(se)];
586                else
587                  std::cout << _graph.id(se) << ":" << -(_cost[se] + pi[_graph.source(se)] - pi[_graph.target(se)]);
588                for (n = u; n != v; n = pred_node[n]) {
589                  if (_graph.target(pred_arc[n]) == n)
590                    std::cout << " " << _cost[pred_arc[n]] + pi[_graph.source(pred_arc[n])] - pi[_graph.target(pred_arc[n])];
591                  else
592                    std::cout << " " << -(_cost[pred_arc[n]] + pi[_graph.source(pred_arc[n])] - pi[_graph.target(pred_arc[n])]);
593                }
594                std::cout << "\n";
595*/
596                // Augment along the cycle
597                (*_flow)[se] = sf ? (*_flow)[se] + delta :
598                                    (*_flow)[se] - delta;
599                for (n = u; n != v; n = pred_node[n]) {
600                  if (_graph.target(pred_arc[n]) == n)
601                    (*_flow)[pred_arc[n]] += delta;
602                  else
603                    (*_flow)[pred_arc[n]] -= delta;
604                }
605                for (n = u; stack_head > 0 && n != w; n = pred_node[n]) {
606                  --stack_head;
607                  reached[n] = false;
608                }
609                u = w;
610              }
611              v = u;
612
613              // Find the next admissible residual outgoing arc
614              double p = pi[v];
615              Arc e = stack[stack_head].first;
616              if (!stack[stack_head].second) {
617                _graph.nextIn(e);
618                goto in_arc_3;
619              }
620              _graph.nextOut(e);
621              while ( e != INVALID && (_capacity[e] - (*_flow)[e] == 0 ||
622                      !tol.negative(_cost[e] + p - pi[_graph.target(e)])) )
623                _graph.nextOut(e);
624              if (e != INVALID) {
625                stack[stack_head] = pair(e, true);
626                goto next_step_3;
627              }
628              _graph.firstIn(e, v);
629            in_arc_3:
630              while ( e != INVALID && ((*_flow)[e] == 0 ||
631                      !tol.negative(-_cost[e] + p - pi[_graph.source(e)])) )
632                _graph.nextIn(e);
633              stack[stack_head] = pair(e, false);
634            next_step_3: ;
635            }
636
637            while (stack_head >= 0 && stack[stack_head].first == INVALID) {
638              processed[v] = true;
639              proc_vector[++proc_head] = v;
640              if (--stack_head >= 0) {
641                v = stack[stack_head].second ?
642                    _graph.source(stack[stack_head].first) :
643                    _graph.target(stack[stack_head].first);
644                // Find the next admissible residual outgoing arc
645                double p = pi[v];
646                Arc e = stack[stack_head].first;
647                if (!stack[stack_head].second) {
648                  _graph.nextIn(e);
649                  goto in_arc_4;
650                }
651                _graph.nextOut(e);
652                while ( e != INVALID && (_capacity[e] - (*_flow)[e] == 0 ||
653                        !tol.negative(_cost[e] + p - pi[_graph.target(e)])) )
654                  _graph.nextOut(e);
655                if (e != INVALID) {
656                  stack[stack_head] = pair(e, true);
657                  goto next_step_4;
658                }
659                _graph.firstIn(e, v);
660              in_arc_4:
661                while ( e != INVALID && ((*_flow)[e] == 0 ||
662                        !tol.negative(-_cost[e] + p - pi[_graph.source(e)])) )
663                  _graph.nextIn(e);
664                stack[stack_head] = pair(e, false);
665              next_step_4: ;
666              }
667            }
668          }
669        }
670        t1.stop();
671
672        // Tighten potentials and epsilon
673        if (--iter > 0) {
674          // Compute levels
675          t2.start();
676          for (int i = proc_head; i >= 0; --i) {
677            Node v = proc_vector[i];
678            double p = pi[v];
679            int l = 0;
680            for (InArcIt e(_graph, v); e != INVALID; ++e) {
681              Node u = _graph.source(e);
682              if ( _capacity[e] - (*_flow)[e] > 0 &&
683                   tol.negative(_cost[e] + pi[u] - p) &&
684                   level[u] + 1 > l ) l = level[u] + 1;
685            }
686            for (OutArcIt e(_graph, v); e != INVALID; ++e) {
687              Node u = _graph.target(e);
688              if ( (*_flow)[e] > 0 &&
689                   tol.negative(-_cost[e] + pi[u] - p) &&
690                   level[u] + 1 > l ) l = level[u] + 1;
691            }
692            level[v] = l;
693          }
694
695          // Modify potentials
696          double p, q = -1;
697          for (ArcIt e(_graph); e != INVALID; ++e) {
698            Node u = _graph.source(e);
699            Node v = _graph.target(e);
700            if (_capacity[e] - (*_flow)[e] > 0 && level[u] - level[v] > 0) {
701              p = (_cost[e] + pi[u] - pi[v] + epsilon) /
702                  (level[u] - level[v] + 1);
703              if (q < 0 || p < q) q = p;
704            }
705            else if ((*_flow)[e] > 0 && level[v] - level[u] > 0) {
706              p = (-_cost[e] - pi[u] + pi[v] + epsilon) /
707                  (level[v] - level[u] + 1);
708              if (q < 0 || p < q) q = p;
709            }
710          }
711          for (NodeIt v(_graph); v != INVALID; ++v) {
712            pi[v] -= q * level[v];
713          }
714
715          // Modify epsilon
716          epsilon = 0;
717          for (ArcIt e(_graph); e != INVALID; ++e) {
718            double curr = _cost[e] + pi[_graph.source(e)]
719                                   - pi[_graph.target(e)];
720            if (_capacity[e] - (*_flow)[e] > 0 && curr < -epsilon)
721              epsilon = -curr;
722            else if ((*_flow)[e] > 0 && curr > epsilon)
723              epsilon = curr;
724          }
725          t2.stop();
726        } else {
727          // Set epsilon to the minimum cycle mean
728          t3.start();
729
730/**/
731          StaticDigraph static_graph;
732          typename ResDigraph::template NodeMap<typename StaticDigraph::Node> node_ref(*_res_graph);
733          typename ResDigraph::template ArcMap<typename StaticDigraph::Arc> arc_ref(*_res_graph);
734          static_graph.build(*_res_graph, node_ref, arc_ref);
735          typename StaticDigraph::template NodeMap<double> static_pi(static_graph);
736          typename StaticDigraph::template ArcMap<double> static_cost(static_graph);
737
738          for (typename ResDigraph::ArcIt e(*_res_graph); e != INVALID; ++e)
739            static_cost[arc_ref[e]] = _res_cost[e];
740
741          Howard<StaticDigraph, typename StaticDigraph::template ArcMap<double> >
742            mmc(static_graph, static_cost);
743          mmc.findMinMean();
744          epsilon = -mmc.cycleMean();
745/**/
746
747/*
748          Howard<ResDigraph, ResidualCostMap> mmc(*_res_graph, _res_cost);
749          mmc.findMinMean();
750          epsilon = -mmc.cycleMean();
751*/
752
753          // Compute feasible potentials for the current epsilon
754          for (typename StaticDigraph::ArcIt e(static_graph); e != INVALID; ++e)
755            static_cost[e] += epsilon;
756          typename BellmanFord<StaticDigraph, typename StaticDigraph::template ArcMap<double> >::
757            template SetDistMap<typename StaticDigraph::template NodeMap<double> >::
758            template SetOperationTraits<BFOperationTraits>::Create
759              bf(static_graph, static_cost);
760          bf.distMap(static_pi).init(0);
761          bf.start();
762          for (NodeIt n(_graph); n != INVALID; ++n)
763            pi[n] = static_pi[node_ref[n]];
764         
765/*
766          for (typename ResDigraph::ArcIt e(*_res_graph); e != INVALID; ++e)
767            shift_cost[e] = _res_cost[e] + epsilon;
768          typename BellmanFord<ResDigraph, ResShiftCostMap>::
769            template SetDistMap<FloatPotentialMap>::
770            template SetOperationTraits<BFOperationTraits>::Create
771              bf(*_res_graph, shift_cost);
772          bf.distMap(pi).init(0);
773          bf.start();
774*/
775
776          iter = limit;
777          t3.stop();
778        }
779      }
780
781//      std::cout << t1.realTime() << " " << t2.realTime() << " " << t3.realTime() << "\n";
782
783      // Handle non-zero lower bounds
784      if (_lower) {
785        for (ArcIt e(_graph); e != INVALID; ++e)
786          (*_flow)[e] += (*_lower)[e];
787      }
788      return true;
789    }
790
791  }; //class CancelAndTighten
792
793  ///@}
794
795} //namespace lemon
796
797#endif //LEMON_CANCEL_AND_TIGHTEN_H
Note: See TracBrowser for help on using the repository browser.