changeset 603 | 425cc8328c0e |
parent 601 | e8349c6f12ca |
child 604 | 8c3112a66878 |
0:fcef725c1143 | 1:7f8258550ebc |
---|---|
26 |
26 |
27 #include <vector> |
27 #include <vector> |
28 #include <limits> |
28 #include <limits> |
29 #include <algorithm> |
29 #include <algorithm> |
30 |
30 |
31 #include <lemon/core.h> |
|
31 #include <lemon/math.h> |
32 #include <lemon/math.h> |
32 |
33 |
33 namespace lemon { |
34 namespace lemon { |
34 |
35 |
35 /// \addtogroup min_cost_flow |
36 /// \addtogroup min_cost_flow |
118 }; |
119 }; |
119 |
120 |
120 private: |
121 private: |
121 |
122 |
122 // References for the original data |
123 // References for the original data |
123 const Digraph &_orig_graph; |
124 const Digraph &_graph; |
124 const LowerMap *_orig_lower; |
125 const LowerMap *_orig_lower; |
125 const CapacityMap &_orig_cap; |
126 const CapacityMap &_orig_cap; |
126 const CostMap &_orig_cost; |
127 const CostMap &_orig_cost; |
127 const SupplyMap *_orig_supply; |
128 const SupplyMap *_orig_supply; |
128 Node _orig_source; |
129 Node _orig_source; |
129 Node _orig_target; |
130 Node _orig_target; |
130 Capacity _orig_flow_value; |
131 Capacity _orig_flow_value; |
131 |
132 |
132 // Result maps |
133 // Result maps |
133 FlowMap *_flow_result; |
134 FlowMap *_flow_map; |
134 PotentialMap *_potential_result; |
135 PotentialMap *_potential_map; |
135 bool _local_flow; |
136 bool _local_flow; |
136 bool _local_potential; |
137 bool _local_potential; |
137 |
|
138 // Data structures for storing the graph |
|
139 ArcVector _arc; |
|
140 NodeVector _node; |
|
141 IntNodeMap _node_id; |
|
142 IntVector _source; |
|
143 IntVector _target; |
|
144 |
138 |
145 // The number of nodes and arcs in the original graph |
139 // The number of nodes and arcs in the original graph |
146 int _node_num; |
140 int _node_num; |
147 int _arc_num; |
141 int _arc_num; |
142 |
|
143 // Data structures for storing the graph |
|
144 IntNodeMap _node_id; |
|
145 ArcVector _arc_ref; |
|
146 IntVector _source; |
|
147 IntVector _target; |
|
148 |
148 |
149 // Node and arc maps |
149 // Node and arc maps |
150 CapacityVector _cap; |
150 CapacityVector _cap; |
151 CostVector _cost; |
151 CostVector _cost; |
152 CostVector _supply; |
152 CostVector _supply; |
153 CapacityVector _flow; |
153 CapacityVector _flow; |
154 CostVector _pi; |
154 CostVector _pi; |
155 |
155 |
156 // Node and arc maps for the spanning tree structure |
156 // Data for storing the spanning tree structure |
157 IntVector _depth; |
157 IntVector _depth; |
158 IntVector _parent; |
158 IntVector _parent; |
159 IntVector _pred; |
159 IntVector _pred; |
160 IntVector _thread; |
160 IntVector _thread; |
161 BoolVector _forward; |
161 BoolVector _forward; |
162 IntVector _state; |
162 IntVector _state; |
163 |
|
164 // The root node |
|
165 int _root; |
163 int _root; |
166 |
164 |
167 // The entering arc in the current pivot iteration |
|
168 int _in_arc; |
|
169 |
|
170 // Temporary data used in the current pivot iteration |
165 // Temporary data used in the current pivot iteration |
171 int join, u_in, v_in, u_out, v_out; |
166 int in_arc, join, u_in, v_in, u_out, v_out; |
172 int right, first, second, last; |
167 int first, second, right, last; |
173 int stem, par_stem, new_stem; |
168 int stem, par_stem, new_stem; |
174 Capacity delta; |
169 Capacity delta; |
175 |
170 |
176 private: |
171 private: |
177 |
172 |
185 class FirstEligiblePivotRule |
180 class FirstEligiblePivotRule |
186 { |
181 { |
187 private: |
182 private: |
188 |
183 |
189 // References to the NetworkSimplex class |
184 // References to the NetworkSimplex class |
190 const ArcVector &_arc; |
|
191 const IntVector &_source; |
185 const IntVector &_source; |
192 const IntVector &_target; |
186 const IntVector &_target; |
193 const CostVector &_cost; |
187 const CostVector &_cost; |
194 const IntVector &_state; |
188 const IntVector &_state; |
195 const CostVector &_pi; |
189 const CostVector &_pi; |
201 |
195 |
202 public: |
196 public: |
203 |
197 |
204 /// Constructor |
198 /// Constructor |
205 FirstEligiblePivotRule(NetworkSimplex &ns) : |
199 FirstEligiblePivotRule(NetworkSimplex &ns) : |
206 _arc(ns._arc), _source(ns._source), _target(ns._target), |
200 _source(ns._source), _target(ns._target), |
207 _cost(ns._cost), _state(ns._state), _pi(ns._pi), |
201 _cost(ns._cost), _state(ns._state), _pi(ns._pi), |
208 _in_arc(ns._in_arc), _arc_num(ns._arc_num), _next_arc(0) |
202 _in_arc(ns.in_arc), _arc_num(ns._arc_num), _next_arc(0) |
209 {} |
203 {} |
210 |
204 |
211 /// Find next entering arc |
205 /// Find next entering arc |
212 bool findEnteringArc() { |
206 bool findEnteringArc() { |
213 Cost c; |
207 Cost c; |
243 class BestEligiblePivotRule |
237 class BestEligiblePivotRule |
244 { |
238 { |
245 private: |
239 private: |
246 |
240 |
247 // References to the NetworkSimplex class |
241 // References to the NetworkSimplex class |
248 const ArcVector &_arc; |
|
249 const IntVector &_source; |
242 const IntVector &_source; |
250 const IntVector &_target; |
243 const IntVector &_target; |
251 const CostVector &_cost; |
244 const CostVector &_cost; |
252 const IntVector &_state; |
245 const IntVector &_state; |
253 const CostVector &_pi; |
246 const CostVector &_pi; |
256 |
249 |
257 public: |
250 public: |
258 |
251 |
259 /// Constructor |
252 /// Constructor |
260 BestEligiblePivotRule(NetworkSimplex &ns) : |
253 BestEligiblePivotRule(NetworkSimplex &ns) : |
261 _arc(ns._arc), _source(ns._source), _target(ns._target), |
254 _source(ns._source), _target(ns._target), |
262 _cost(ns._cost), _state(ns._state), _pi(ns._pi), |
255 _cost(ns._cost), _state(ns._state), _pi(ns._pi), |
263 _in_arc(ns._in_arc), _arc_num(ns._arc_num) |
256 _in_arc(ns.in_arc), _arc_num(ns._arc_num) |
264 {} |
257 {} |
265 |
258 |
266 /// Find next entering arc |
259 /// Find next entering arc |
267 bool findEnteringArc() { |
260 bool findEnteringArc() { |
268 Cost c, min = 0; |
261 Cost c, min = 0; |
289 class BlockSearchPivotRule |
282 class BlockSearchPivotRule |
290 { |
283 { |
291 private: |
284 private: |
292 |
285 |
293 // References to the NetworkSimplex class |
286 // References to the NetworkSimplex class |
294 const ArcVector &_arc; |
|
295 const IntVector &_source; |
287 const IntVector &_source; |
296 const IntVector &_target; |
288 const IntVector &_target; |
297 const CostVector &_cost; |
289 const CostVector &_cost; |
298 const IntVector &_state; |
290 const IntVector &_state; |
299 const CostVector &_pi; |
291 const CostVector &_pi; |
306 |
298 |
307 public: |
299 public: |
308 |
300 |
309 /// Constructor |
301 /// Constructor |
310 BlockSearchPivotRule(NetworkSimplex &ns) : |
302 BlockSearchPivotRule(NetworkSimplex &ns) : |
311 _arc(ns._arc), _source(ns._source), _target(ns._target), |
303 _source(ns._source), _target(ns._target), |
312 _cost(ns._cost), _state(ns._state), _pi(ns._pi), |
304 _cost(ns._cost), _state(ns._state), _pi(ns._pi), |
313 _in_arc(ns._in_arc), _arc_num(ns._arc_num), _next_arc(0) |
305 _in_arc(ns.in_arc), _arc_num(ns._arc_num), _next_arc(0) |
314 { |
306 { |
315 // The main parameters of the pivot rule |
307 // The main parameters of the pivot rule |
316 const double BLOCK_SIZE_FACTOR = 2.0; |
308 const double BLOCK_SIZE_FACTOR = 2.0; |
317 const int MIN_BLOCK_SIZE = 10; |
309 const int MIN_BLOCK_SIZE = 10; |
318 |
310 |
368 class CandidateListPivotRule |
360 class CandidateListPivotRule |
369 { |
361 { |
370 private: |
362 private: |
371 |
363 |
372 // References to the NetworkSimplex class |
364 // References to the NetworkSimplex class |
373 const ArcVector &_arc; |
|
374 const IntVector &_source; |
365 const IntVector &_source; |
375 const IntVector &_target; |
366 const IntVector &_target; |
376 const CostVector &_cost; |
367 const CostVector &_cost; |
377 const IntVector &_state; |
368 const IntVector &_state; |
378 const CostVector &_pi; |
369 const CostVector &_pi; |
387 |
378 |
388 public: |
379 public: |
389 |
380 |
390 /// Constructor |
381 /// Constructor |
391 CandidateListPivotRule(NetworkSimplex &ns) : |
382 CandidateListPivotRule(NetworkSimplex &ns) : |
392 _arc(ns._arc), _source(ns._source), _target(ns._target), |
383 _source(ns._source), _target(ns._target), |
393 _cost(ns._cost), _state(ns._state), _pi(ns._pi), |
384 _cost(ns._cost), _state(ns._state), _pi(ns._pi), |
394 _in_arc(ns._in_arc), _arc_num(ns._arc_num), _next_arc(0) |
385 _in_arc(ns.in_arc), _arc_num(ns._arc_num), _next_arc(0) |
395 { |
386 { |
396 // The main parameters of the pivot rule |
387 // The main parameters of the pivot rule |
397 const double LIST_LENGTH_FACTOR = 1.0; |
388 const double LIST_LENGTH_FACTOR = 1.0; |
398 const int MIN_LIST_LENGTH = 10; |
389 const int MIN_LIST_LENGTH = 10; |
399 const double MINOR_LIMIT_FACTOR = 0.1; |
390 const double MINOR_LIMIT_FACTOR = 0.1; |
480 class AlteringListPivotRule |
471 class AlteringListPivotRule |
481 { |
472 { |
482 private: |
473 private: |
483 |
474 |
484 // References to the NetworkSimplex class |
475 // References to the NetworkSimplex class |
485 const ArcVector &_arc; |
|
486 const IntVector &_source; |
476 const IntVector &_source; |
487 const IntVector &_target; |
477 const IntVector &_target; |
488 const CostVector &_cost; |
478 const CostVector &_cost; |
489 const IntVector &_state; |
479 const IntVector &_state; |
490 const CostVector &_pi; |
480 const CostVector &_pi; |
513 |
503 |
514 public: |
504 public: |
515 |
505 |
516 /// Constructor |
506 /// Constructor |
517 AlteringListPivotRule(NetworkSimplex &ns) : |
507 AlteringListPivotRule(NetworkSimplex &ns) : |
518 _arc(ns._arc), _source(ns._source), _target(ns._target), |
508 _source(ns._source), _target(ns._target), |
519 _cost(ns._cost), _state(ns._state), _pi(ns._pi), |
509 _cost(ns._cost), _state(ns._state), _pi(ns._pi), |
520 _in_arc(ns._in_arc), _arc_num(ns._arc_num), |
510 _in_arc(ns.in_arc), _arc_num(ns._arc_num), |
521 _next_arc(0), _cand_cost(ns._arc_num), _sort_func(_cand_cost) |
511 _next_arc(0), _cand_cost(ns._arc_num), _sort_func(_cand_cost) |
522 { |
512 { |
523 // The main parameters of the pivot rule |
513 // The main parameters of the pivot rule |
524 const double BLOCK_SIZE_FACTOR = 1.5; |
514 const double BLOCK_SIZE_FACTOR = 1.5; |
525 const int MIN_BLOCK_SIZE = 10; |
515 const int MIN_BLOCK_SIZE = 10; |
547 } |
537 } |
548 } |
538 } |
549 |
539 |
550 // Extend the list |
540 // Extend the list |
551 int cnt = _block_size; |
541 int cnt = _block_size; |
552 int last_edge = 0; |
542 int last_arc = 0; |
553 int limit = _head_length; |
543 int limit = _head_length; |
554 |
544 |
555 for (int e = _next_arc; e < _arc_num; ++e) { |
545 for (int e = _next_arc; e < _arc_num; ++e) { |
556 _cand_cost[e] = _state[e] * |
546 _cand_cost[e] = _state[e] * |
557 (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
547 (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
558 if (_cand_cost[e] < 0) { |
548 if (_cand_cost[e] < 0) { |
559 _candidates[_curr_length++] = e; |
549 _candidates[_curr_length++] = e; |
560 last_edge = e; |
550 last_arc = e; |
561 } |
551 } |
562 if (--cnt == 0) { |
552 if (--cnt == 0) { |
563 if (_curr_length > limit) break; |
553 if (_curr_length > limit) break; |
564 limit = 0; |
554 limit = 0; |
565 cnt = _block_size; |
555 cnt = _block_size; |
569 for (int e = 0; e < _next_arc; ++e) { |
559 for (int e = 0; e < _next_arc; ++e) { |
570 _cand_cost[e] = _state[e] * |
560 _cand_cost[e] = _state[e] * |
571 (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
561 (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
572 if (_cand_cost[e] < 0) { |
562 if (_cand_cost[e] < 0) { |
573 _candidates[_curr_length++] = e; |
563 _candidates[_curr_length++] = e; |
574 last_edge = e; |
564 last_arc = e; |
575 } |
565 } |
576 if (--cnt == 0) { |
566 if (--cnt == 0) { |
577 if (_curr_length > limit) break; |
567 if (_curr_length > limit) break; |
578 limit = 0; |
568 limit = 0; |
579 cnt = _block_size; |
569 cnt = _block_size; |
580 } |
570 } |
581 } |
571 } |
582 } |
572 } |
583 if (_curr_length == 0) return false; |
573 if (_curr_length == 0) return false; |
584 _next_arc = last_edge + 1; |
574 _next_arc = last_arc + 1; |
585 |
575 |
586 // Make heap of the candidate list (approximating a partial sort) |
576 // Make heap of the candidate list (approximating a partial sort) |
587 make_heap( _candidates.begin(), _candidates.begin() + _curr_length, |
577 make_heap( _candidates.begin(), _candidates.begin() + _curr_length, |
588 _sort_func ); |
578 _sort_func ); |
589 |
579 |
601 |
591 |
602 /// \brief General constructor (with lower bounds). |
592 /// \brief General constructor (with lower bounds). |
603 /// |
593 /// |
604 /// General constructor (with lower bounds). |
594 /// General constructor (with lower bounds). |
605 /// |
595 /// |
606 /// \param digraph The digraph the algorithm runs on. |
596 /// \param graph The digraph the algorithm runs on. |
607 /// \param lower The lower bounds of the arcs. |
597 /// \param lower The lower bounds of the arcs. |
608 /// \param capacity The capacities (upper bounds) of the arcs. |
598 /// \param capacity The capacities (upper bounds) of the arcs. |
609 /// \param cost The cost (length) values of the arcs. |
599 /// \param cost The cost (length) values of the arcs. |
610 /// \param supply The supply values of the nodes (signed). |
600 /// \param supply The supply values of the nodes (signed). |
611 NetworkSimplex( const Digraph &digraph, |
601 NetworkSimplex( const Digraph &graph, |
612 const LowerMap &lower, |
602 const LowerMap &lower, |
613 const CapacityMap &capacity, |
603 const CapacityMap &capacity, |
614 const CostMap &cost, |
604 const CostMap &cost, |
615 const SupplyMap &supply ) : |
605 const SupplyMap &supply ) : |
616 _orig_graph(digraph), _orig_lower(&lower), _orig_cap(capacity), |
606 _graph(graph), _orig_lower(&lower), _orig_cap(capacity), |
617 _orig_cost(cost), _orig_supply(&supply), |
607 _orig_cost(cost), _orig_supply(&supply), |
618 _flow_result(NULL), _potential_result(NULL), |
608 _flow_map(NULL), _potential_map(NULL), |
619 _local_flow(false), _local_potential(false), |
609 _local_flow(false), _local_potential(false), |
620 _node_id(digraph) |
610 _node_id(graph) |
621 {} |
611 {} |
622 |
612 |
623 /// \brief General constructor (without lower bounds). |
613 /// \brief General constructor (without lower bounds). |
624 /// |
614 /// |
625 /// General constructor (without lower bounds). |
615 /// General constructor (without lower bounds). |
626 /// |
616 /// |
627 /// \param digraph The digraph the algorithm runs on. |
617 /// \param graph The digraph the algorithm runs on. |
628 /// \param capacity The capacities (upper bounds) of the arcs. |
618 /// \param capacity The capacities (upper bounds) of the arcs. |
629 /// \param cost The cost (length) values of the arcs. |
619 /// \param cost The cost (length) values of the arcs. |
630 /// \param supply The supply values of the nodes (signed). |
620 /// \param supply The supply values of the nodes (signed). |
631 NetworkSimplex( const Digraph &digraph, |
621 NetworkSimplex( const Digraph &graph, |
632 const CapacityMap &capacity, |
622 const CapacityMap &capacity, |
633 const CostMap &cost, |
623 const CostMap &cost, |
634 const SupplyMap &supply ) : |
624 const SupplyMap &supply ) : |
635 _orig_graph(digraph), _orig_lower(NULL), _orig_cap(capacity), |
625 _graph(graph), _orig_lower(NULL), _orig_cap(capacity), |
636 _orig_cost(cost), _orig_supply(&supply), |
626 _orig_cost(cost), _orig_supply(&supply), |
637 _flow_result(NULL), _potential_result(NULL), |
627 _flow_map(NULL), _potential_map(NULL), |
638 _local_flow(false), _local_potential(false), |
628 _local_flow(false), _local_potential(false), |
639 _node_id(digraph) |
629 _node_id(graph) |
640 {} |
630 {} |
641 |
631 |
642 /// \brief Simple constructor (with lower bounds). |
632 /// \brief Simple constructor (with lower bounds). |
643 /// |
633 /// |
644 /// Simple constructor (with lower bounds). |
634 /// Simple constructor (with lower bounds). |
645 /// |
635 /// |
646 /// \param digraph The digraph the algorithm runs on. |
636 /// \param graph The digraph the algorithm runs on. |
647 /// \param lower The lower bounds of the arcs. |
637 /// \param lower The lower bounds of the arcs. |
648 /// \param capacity The capacities (upper bounds) of the arcs. |
638 /// \param capacity The capacities (upper bounds) of the arcs. |
649 /// \param cost The cost (length) values of the arcs. |
639 /// \param cost The cost (length) values of the arcs. |
650 /// \param s The source node. |
640 /// \param s The source node. |
651 /// \param t The target node. |
641 /// \param t The target node. |
652 /// \param flow_value The required amount of flow from node \c s |
642 /// \param flow_value The required amount of flow from node \c s |
653 /// to node \c t (i.e. the supply of \c s and the demand of \c t). |
643 /// to node \c t (i.e. the supply of \c s and the demand of \c t). |
654 NetworkSimplex( const Digraph &digraph, |
644 NetworkSimplex( const Digraph &graph, |
655 const LowerMap &lower, |
645 const LowerMap &lower, |
656 const CapacityMap &capacity, |
646 const CapacityMap &capacity, |
657 const CostMap &cost, |
647 const CostMap &cost, |
658 Node s, Node t, |
648 Node s, Node t, |
659 Capacity flow_value ) : |
649 Capacity flow_value ) : |
660 _orig_graph(digraph), _orig_lower(&lower), _orig_cap(capacity), |
650 _graph(graph), _orig_lower(&lower), _orig_cap(capacity), |
661 _orig_cost(cost), _orig_supply(NULL), |
651 _orig_cost(cost), _orig_supply(NULL), |
662 _orig_source(s), _orig_target(t), _orig_flow_value(flow_value), |
652 _orig_source(s), _orig_target(t), _orig_flow_value(flow_value), |
663 _flow_result(NULL), _potential_result(NULL), |
653 _flow_map(NULL), _potential_map(NULL), |
664 _local_flow(false), _local_potential(false), |
654 _local_flow(false), _local_potential(false), |
665 _node_id(digraph) |
655 _node_id(graph) |
666 {} |
656 {} |
667 |
657 |
668 /// \brief Simple constructor (without lower bounds). |
658 /// \brief Simple constructor (without lower bounds). |
669 /// |
659 /// |
670 /// Simple constructor (without lower bounds). |
660 /// Simple constructor (without lower bounds). |
671 /// |
661 /// |
672 /// \param digraph The digraph the algorithm runs on. |
662 /// \param graph The digraph the algorithm runs on. |
673 /// \param capacity The capacities (upper bounds) of the arcs. |
663 /// \param capacity The capacities (upper bounds) of the arcs. |
674 /// \param cost The cost (length) values of the arcs. |
664 /// \param cost The cost (length) values of the arcs. |
675 /// \param s The source node. |
665 /// \param s The source node. |
676 /// \param t The target node. |
666 /// \param t The target node. |
677 /// \param flow_value The required amount of flow from node \c s |
667 /// \param flow_value The required amount of flow from node \c s |
678 /// to node \c t (i.e. the supply of \c s and the demand of \c t). |
668 /// to node \c t (i.e. the supply of \c s and the demand of \c t). |
679 NetworkSimplex( const Digraph &digraph, |
669 NetworkSimplex( const Digraph &graph, |
680 const CapacityMap &capacity, |
670 const CapacityMap &capacity, |
681 const CostMap &cost, |
671 const CostMap &cost, |
682 Node s, Node t, |
672 Node s, Node t, |
683 Capacity flow_value ) : |
673 Capacity flow_value ) : |
684 _orig_graph(digraph), _orig_lower(NULL), _orig_cap(capacity), |
674 _graph(graph), _orig_lower(NULL), _orig_cap(capacity), |
685 _orig_cost(cost), _orig_supply(NULL), |
675 _orig_cost(cost), _orig_supply(NULL), |
686 _orig_source(s), _orig_target(t), _orig_flow_value(flow_value), |
676 _orig_source(s), _orig_target(t), _orig_flow_value(flow_value), |
687 _flow_result(NULL), _potential_result(NULL), |
677 _flow_map(NULL), _potential_map(NULL), |
688 _local_flow(false), _local_potential(false), |
678 _local_flow(false), _local_potential(false), |
689 _node_id(digraph) |
679 _node_id(graph) |
690 {} |
680 {} |
691 |
681 |
692 /// Destructor. |
682 /// Destructor. |
693 ~NetworkSimplex() { |
683 ~NetworkSimplex() { |
694 if (_local_flow) delete _flow_result; |
684 if (_local_flow) delete _flow_map; |
695 if (_local_potential) delete _potential_result; |
685 if (_local_potential) delete _potential_map; |
696 } |
686 } |
697 |
687 |
698 /// \brief Set the flow map. |
688 /// \brief Set the flow map. |
699 /// |
689 /// |
700 /// This function sets the flow map. |
690 /// This function sets the flow map. |
701 /// |
691 /// |
702 /// \return <tt>(*this)</tt> |
692 /// \return <tt>(*this)</tt> |
703 NetworkSimplex& flowMap(FlowMap &map) { |
693 NetworkSimplex& flowMap(FlowMap &map) { |
704 if (_local_flow) { |
694 if (_local_flow) { |
705 delete _flow_result; |
695 delete _flow_map; |
706 _local_flow = false; |
696 _local_flow = false; |
707 } |
697 } |
708 _flow_result = ↦ |
698 _flow_map = ↦ |
709 return *this; |
699 return *this; |
710 } |
700 } |
711 |
701 |
712 /// \brief Set the potential map. |
702 /// \brief Set the potential map. |
713 /// |
703 /// |
714 /// This function sets the potential map. |
704 /// This function sets the potential map. |
715 /// |
705 /// |
716 /// \return <tt>(*this)</tt> |
706 /// \return <tt>(*this)</tt> |
717 NetworkSimplex& potentialMap(PotentialMap &map) { |
707 NetworkSimplex& potentialMap(PotentialMap &map) { |
718 if (_local_potential) { |
708 if (_local_potential) { |
719 delete _potential_result; |
709 delete _potential_map; |
720 _local_potential = false; |
710 _local_potential = false; |
721 } |
711 } |
722 _potential_result = ↦ |
712 _potential_map = ↦ |
723 return *this; |
713 return *this; |
724 } |
714 } |
725 |
715 |
726 /// \name Execution control |
716 /// \name Execution control |
727 /// The algorithm can be executed using the |
717 /// The algorithm can be executed using the |
781 /// This function returns a const reference to an arc map storing |
771 /// This function returns a const reference to an arc map storing |
782 /// the found flow. |
772 /// the found flow. |
783 /// |
773 /// |
784 /// \pre \ref run() must be called before using this function. |
774 /// \pre \ref run() must be called before using this function. |
785 const FlowMap& flowMap() const { |
775 const FlowMap& flowMap() const { |
786 return *_flow_result; |
776 return *_flow_map; |
787 } |
777 } |
788 |
778 |
789 /// \brief Return a const reference to the potential map |
779 /// \brief Return a const reference to the potential map |
790 /// (the dual solution). |
780 /// (the dual solution). |
791 /// |
781 /// |
792 /// This function returns a const reference to a node map storing |
782 /// This function returns a const reference to a node map storing |
793 /// the found potentials (the dual solution). |
783 /// the found potentials (the dual solution). |
794 /// |
784 /// |
795 /// \pre \ref run() must be called before using this function. |
785 /// \pre \ref run() must be called before using this function. |
796 const PotentialMap& potentialMap() const { |
786 const PotentialMap& potentialMap() const { |
797 return *_potential_result; |
787 return *_potential_map; |
798 } |
788 } |
799 |
789 |
800 /// \brief Return the flow on the given arc. |
790 /// \brief Return the flow on the given arc. |
801 /// |
791 /// |
802 /// This function returns the flow on the given arc. |
792 /// This function returns the flow on the given arc. |
803 /// |
793 /// |
804 /// \pre \ref run() must be called before using this function. |
794 /// \pre \ref run() must be called before using this function. |
805 Capacity flow(const Arc& arc) const { |
795 Capacity flow(const Arc& arc) const { |
806 return (*_flow_result)[arc]; |
796 return (*_flow_map)[arc]; |
807 } |
797 } |
808 |
798 |
809 /// \brief Return the potential of the given node. |
799 /// \brief Return the potential of the given node. |
810 /// |
800 /// |
811 /// This function returns the potential of the given node. |
801 /// This function returns the potential of the given node. |
812 /// |
802 /// |
813 /// \pre \ref run() must be called before using this function. |
803 /// \pre \ref run() must be called before using this function. |
814 Cost potential(const Node& node) const { |
804 Cost potential(const Node& node) const { |
815 return (*_potential_result)[node]; |
805 return (*_potential_map)[node]; |
816 } |
806 } |
817 |
807 |
818 /// \brief Return the total cost of the found flow. |
808 /// \brief Return the total cost of the found flow. |
819 /// |
809 /// |
820 /// This function returns the total cost of the found flow. |
810 /// This function returns the total cost of the found flow. |
821 /// The complexity of the function is \f$ O(e) \f$. |
811 /// The complexity of the function is \f$ O(e) \f$. |
822 /// |
812 /// |
823 /// \pre \ref run() must be called before using this function. |
813 /// \pre \ref run() must be called before using this function. |
824 Cost totalCost() const { |
814 Cost totalCost() const { |
825 Cost c = 0; |
815 Cost c = 0; |
826 for (ArcIt e(_orig_graph); e != INVALID; ++e) |
816 for (ArcIt e(_graph); e != INVALID; ++e) |
827 c += (*_flow_result)[e] * _orig_cost[e]; |
817 c += (*_flow_map)[e] * _orig_cost[e]; |
828 return c; |
818 return c; |
829 } |
819 } |
830 |
820 |
831 /// @} |
821 /// @} |
832 |
822 |
833 private: |
823 private: |
834 |
824 |
835 // Initialize internal data structures |
825 // Initialize internal data structures |
836 bool init() { |
826 bool init() { |
837 // Initialize result maps |
827 // Initialize result maps |
838 if (!_flow_result) { |
828 if (!_flow_map) { |
839 _flow_result = new FlowMap(_orig_graph); |
829 _flow_map = new FlowMap(_graph); |
840 _local_flow = true; |
830 _local_flow = true; |
841 } |
831 } |
842 if (!_potential_result) { |
832 if (!_potential_map) { |
843 _potential_result = new PotentialMap(_orig_graph); |
833 _potential_map = new PotentialMap(_graph); |
844 _local_potential = true; |
834 _local_potential = true; |
845 } |
835 } |
846 |
836 |
847 // Initialize vectors |
837 // Initialize vectors |
848 _node_num = countNodes(_orig_graph); |
838 _node_num = countNodes(_graph); |
849 _arc_num = countArcs(_orig_graph); |
839 _arc_num = countArcs(_graph); |
850 int all_node_num = _node_num + 1; |
840 int all_node_num = _node_num + 1; |
851 int all_edge_num = _arc_num + _node_num; |
841 int all_arc_num = _arc_num + _node_num; |
852 |
842 |
853 _arc.resize(_arc_num); |
843 _arc_ref.resize(_arc_num); |
854 _node.reserve(_node_num); |
844 _source.resize(all_arc_num); |
855 _source.resize(all_edge_num); |
845 _target.resize(all_arc_num); |
856 _target.resize(all_edge_num); |
846 |
857 |
847 _cap.resize(all_arc_num); |
858 _cap.resize(all_edge_num); |
848 _cost.resize(all_arc_num); |
859 _cost.resize(all_edge_num); |
|
860 _supply.resize(all_node_num); |
849 _supply.resize(all_node_num); |
861 _flow.resize(all_edge_num, 0); |
850 _flow.resize(all_arc_num, 0); |
862 _pi.resize(all_node_num, 0); |
851 _pi.resize(all_node_num, 0); |
863 |
852 |
864 _depth.resize(all_node_num); |
853 _depth.resize(all_node_num); |
865 _parent.resize(all_node_num); |
854 _parent.resize(all_node_num); |
866 _pred.resize(all_node_num); |
855 _pred.resize(all_node_num); |
856 _forward.resize(all_node_num); |
|
867 _thread.resize(all_node_num); |
857 _thread.resize(all_node_num); |
868 _forward.resize(all_node_num); |
858 _state.resize(all_arc_num, STATE_LOWER); |
869 _state.resize(all_edge_num, STATE_LOWER); |
|
870 |
859 |
871 // Initialize node related data |
860 // Initialize node related data |
872 bool valid_supply = true; |
861 bool valid_supply = true; |
873 if (_orig_supply) { |
862 if (_orig_supply) { |
874 Supply sum = 0; |
863 Supply sum = 0; |
875 int i = 0; |
864 int i = 0; |
876 for (NodeIt n(_orig_graph); n != INVALID; ++n, ++i) { |
865 for (NodeIt n(_graph); n != INVALID; ++n, ++i) { |
877 _node.push_back(n); |
|
878 _node_id[n] = i; |
866 _node_id[n] = i; |
879 _supply[i] = (*_orig_supply)[n]; |
867 _supply[i] = (*_orig_supply)[n]; |
880 sum += _supply[i]; |
868 sum += _supply[i]; |
881 } |
869 } |
882 valid_supply = (sum == 0); |
870 valid_supply = (sum == 0); |
883 } else { |
871 } else { |
884 int i = 0; |
872 int i = 0; |
885 for (NodeIt n(_orig_graph); n != INVALID; ++n, ++i) { |
873 for (NodeIt n(_graph); n != INVALID; ++n, ++i) { |
886 _node.push_back(n); |
|
887 _node_id[n] = i; |
874 _node_id[n] = i; |
888 _supply[i] = 0; |
875 _supply[i] = 0; |
889 } |
876 } |
890 _supply[_node_id[_orig_source]] = _orig_flow_value; |
877 _supply[_node_id[_orig_source]] = _orig_flow_value; |
891 _supply[_node_id[_orig_target]] = -_orig_flow_value; |
878 _supply[_node_id[_orig_target]] = -_orig_flow_value; |
902 _pi[_root] = 0; |
889 _pi[_root] = 0; |
903 |
890 |
904 // Store the arcs in a mixed order |
891 // Store the arcs in a mixed order |
905 int k = std::max(int(sqrt(_arc_num)), 10); |
892 int k = std::max(int(sqrt(_arc_num)), 10); |
906 int i = 0; |
893 int i = 0; |
907 for (ArcIt e(_orig_graph); e != INVALID; ++e) { |
894 for (ArcIt e(_graph); e != INVALID; ++e) { |
908 _arc[i] = e; |
895 _arc_ref[i] = e; |
909 if ((i += k) >= _arc_num) i = (i % k) + 1; |
896 if ((i += k) >= _arc_num) i = (i % k) + 1; |
910 } |
897 } |
911 |
898 |
912 // Initialize arc maps |
899 // Initialize arc maps |
913 for (int i = 0; i != _arc_num; ++i) { |
900 for (int i = 0; i != _arc_num; ++i) { |
914 Arc e = _arc[i]; |
901 Arc e = _arc_ref[i]; |
915 _source[i] = _node_id[_orig_graph.source(e)]; |
902 _source[i] = _node_id[_graph.source(e)]; |
916 _target[i] = _node_id[_orig_graph.target(e)]; |
903 _target[i] = _node_id[_graph.target(e)]; |
917 _cost[i] = _orig_cost[e]; |
904 _cost[i] = _orig_cost[e]; |
918 _cap[i] = _orig_cap[e]; |
905 _cap[i] = _orig_cap[e]; |
919 } |
906 } |
920 |
907 |
921 // Remove non-zero lower bounds |
908 // Remove non-zero lower bounds |
922 if (_orig_lower) { |
909 if (_orig_lower) { |
923 for (int i = 0; i != _arc_num; ++i) { |
910 for (int i = 0; i != _arc_num; ++i) { |
924 Capacity c = (*_orig_lower)[_arc[i]]; |
911 Capacity c = (*_orig_lower)[_arc_ref[i]]; |
925 if (c != 0) { |
912 if (c != 0) { |
926 _cap[i] -= c; |
913 _cap[i] -= c; |
927 _supply[_source[i]] -= c; |
914 _supply[_source[i]] -= c; |
928 _supply[_target[i]] += c; |
915 _supply[_target[i]] += c; |
929 } |
916 } |
955 return true; |
942 return true; |
956 } |
943 } |
957 |
944 |
958 // Find the join node |
945 // Find the join node |
959 void findJoinNode() { |
946 void findJoinNode() { |
960 int u = _source[_in_arc]; |
947 int u = _source[in_arc]; |
961 int v = _target[_in_arc]; |
948 int v = _target[in_arc]; |
962 while (_depth[u] > _depth[v]) u = _parent[u]; |
949 while (_depth[u] > _depth[v]) u = _parent[u]; |
963 while (_depth[v] > _depth[u]) v = _parent[v]; |
950 while (_depth[v] > _depth[u]) v = _parent[v]; |
964 while (u != v) { |
951 while (u != v) { |
965 u = _parent[u]; |
952 u = _parent[u]; |
966 v = _parent[v]; |
953 v = _parent[v]; |
971 // Find the leaving arc of the cycle and returns true if the |
958 // Find the leaving arc of the cycle and returns true if the |
972 // leaving arc is not the same as the entering arc |
959 // leaving arc is not the same as the entering arc |
973 bool findLeavingArc() { |
960 bool findLeavingArc() { |
974 // Initialize first and second nodes according to the direction |
961 // Initialize first and second nodes according to the direction |
975 // of the cycle |
962 // of the cycle |
976 if (_state[_in_arc] == STATE_LOWER) { |
963 if (_state[in_arc] == STATE_LOWER) { |
977 first = _source[_in_arc]; |
964 first = _source[in_arc]; |
978 second = _target[_in_arc]; |
965 second = _target[in_arc]; |
979 } else { |
966 } else { |
980 first = _target[_in_arc]; |
967 first = _target[in_arc]; |
981 second = _source[_in_arc]; |
968 second = _source[in_arc]; |
982 } |
969 } |
983 delta = _cap[_in_arc]; |
970 delta = _cap[in_arc]; |
984 int result = 0; |
971 int result = 0; |
985 Capacity d; |
972 Capacity d; |
986 int e; |
973 int e; |
987 |
974 |
988 // Search the cycle along the path form the first node to the root |
975 // Search the cycle along the path form the first node to the root |
1018 |
1005 |
1019 // Change _flow and _state vectors |
1006 // Change _flow and _state vectors |
1020 void changeFlow(bool change) { |
1007 void changeFlow(bool change) { |
1021 // Augment along the cycle |
1008 // Augment along the cycle |
1022 if (delta > 0) { |
1009 if (delta > 0) { |
1023 Capacity val = _state[_in_arc] * delta; |
1010 Capacity val = _state[in_arc] * delta; |
1024 _flow[_in_arc] += val; |
1011 _flow[in_arc] += val; |
1025 for (int u = _source[_in_arc]; u != join; u = _parent[u]) { |
1012 for (int u = _source[in_arc]; u != join; u = _parent[u]) { |
1026 _flow[_pred[u]] += _forward[u] ? -val : val; |
1013 _flow[_pred[u]] += _forward[u] ? -val : val; |
1027 } |
1014 } |
1028 for (int u = _target[_in_arc]; u != join; u = _parent[u]) { |
1015 for (int u = _target[in_arc]; u != join; u = _parent[u]) { |
1029 _flow[_pred[u]] += _forward[u] ? val : -val; |
1016 _flow[_pred[u]] += _forward[u] ? val : -val; |
1030 } |
1017 } |
1031 } |
1018 } |
1032 // Update the state of the entering and leaving arcs |
1019 // Update the state of the entering and leaving arcs |
1033 if (change) { |
1020 if (change) { |
1034 _state[_in_arc] = STATE_TREE; |
1021 _state[in_arc] = STATE_TREE; |
1035 _state[_pred[u_out]] = |
1022 _state[_pred[u_out]] = |
1036 (_flow[_pred[u_out]] == 0) ? STATE_LOWER : STATE_UPPER; |
1023 (_flow[_pred[u_out]] == 0) ? STATE_LOWER : STATE_UPPER; |
1037 } else { |
1024 } else { |
1038 _state[_in_arc] = -_state[_in_arc]; |
1025 _state[in_arc] = -_state[in_arc]; |
1039 } |
1026 } |
1040 } |
1027 } |
1041 |
1028 |
1042 // Update _thread and _parent vectors |
1029 // Update _thread and _parent vectors |
1043 void updateThreadParent() { |
1030 void updateThreadParent() { |
1104 v = _parent[u]; |
1091 v = _parent[u]; |
1105 _pred[u] = _pred[v]; |
1092 _pred[u] = _pred[v]; |
1106 _forward[u] = !_forward[v]; |
1093 _forward[u] = !_forward[v]; |
1107 u = v; |
1094 u = v; |
1108 } |
1095 } |
1109 _pred[u_in] = _in_arc; |
1096 _pred[u_in] = in_arc; |
1110 _forward[u_in] = (u_in == _source[_in_arc]); |
1097 _forward[u_in] = (u_in == _source[in_arc]); |
1111 } |
1098 } |
1112 |
1099 |
1113 // Update _depth and _potential vectors |
1100 // Update _depth and _potential vectors |
1114 void updateDepthPotential() { |
1101 void updateDepthPotential() { |
1115 _depth[u_in] = _depth[v_in] + 1; |
1102 _depth[u_in] = _depth[v_in] + 1; |
1161 // Check if the flow amount equals zero on all the artificial arcs |
1148 // Check if the flow amount equals zero on all the artificial arcs |
1162 for (int e = _arc_num; e != _arc_num + _node_num; ++e) { |
1149 for (int e = _arc_num; e != _arc_num + _node_num; ++e) { |
1163 if (_flow[e] > 0) return false; |
1150 if (_flow[e] > 0) return false; |
1164 } |
1151 } |
1165 |
1152 |
1166 // Copy flow values to _flow_result |
1153 // Copy flow values to _flow_map |
1167 if (_orig_lower) { |
1154 if (_orig_lower) { |
1168 for (int i = 0; i != _arc_num; ++i) { |
1155 for (int i = 0; i != _arc_num; ++i) { |
1169 Arc e = _arc[i]; |
1156 Arc e = _arc_ref[i]; |
1170 (*_flow_result)[e] = (*_orig_lower)[e] + _flow[i]; |
1157 _flow_map->set(e, (*_orig_lower)[e] + _flow[i]); |
1171 } |
1158 } |
1172 } else { |
1159 } else { |
1173 for (int i = 0; i != _arc_num; ++i) { |
1160 for (int i = 0; i != _arc_num; ++i) { |
1174 (*_flow_result)[_arc[i]] = _flow[i]; |
1161 _flow_map->set(_arc_ref[i], _flow[i]); |
1175 } |
1162 } |
1176 } |
1163 } |
1177 // Copy potential values to _potential_result |
1164 // Copy potential values to _potential_map |
1178 for (int i = 0; i != _node_num; ++i) { |
1165 for (NodeIt n(_graph); n != INVALID; ++n) { |
1179 (*_potential_result)[_node[i]] = _pi[i]; |
1166 _potential_map->set(n, _pi[_node_id[n]]); |
1180 } |
1167 } |
1181 |
1168 |
1182 return true; |
1169 return true; |
1183 } |
1170 } |
1184 |
1171 |