194 class ResidualDijkstra |
194 class ResidualDijkstra |
195 { |
195 { |
196 private: |
196 private: |
197 |
197 |
198 int _node_num; |
198 int _node_num; |
|
199 bool _geq; |
199 const IntVector &_first_out; |
200 const IntVector &_first_out; |
200 const IntVector &_target; |
201 const IntVector &_target; |
201 const CostVector &_cost; |
202 const CostVector &_cost; |
202 const ValueVector &_res_cap; |
203 const ValueVector &_res_cap; |
203 const ValueVector &_excess; |
204 const ValueVector &_excess; |
208 CostVector _dist; |
209 CostVector _dist; |
209 |
210 |
210 public: |
211 public: |
211 |
212 |
212 ResidualDijkstra(CapacityScaling& cs) : |
213 ResidualDijkstra(CapacityScaling& cs) : |
213 _node_num(cs._node_num), _first_out(cs._first_out), |
214 _node_num(cs._node_num), _geq(cs._sum_supply < 0), |
214 _target(cs._target), _cost(cs._cost), _res_cap(cs._res_cap), |
215 _first_out(cs._first_out), _target(cs._target), _cost(cs._cost), |
215 _excess(cs._excess), _pi(cs._pi), _pred(cs._pred), |
216 _res_cap(cs._res_cap), _excess(cs._excess), _pi(cs._pi), |
216 _dist(cs._node_num) |
217 _pred(cs._pred), _dist(cs._node_num) |
217 {} |
218 {} |
218 |
219 |
219 int run(int s, Value delta = 1) { |
220 int run(int s, Value delta = 1) { |
220 RangeMap<int> heap_cross_ref(_node_num, Heap::PRE_HEAP); |
221 RangeMap<int> heap_cross_ref(_node_num, Heap::PRE_HEAP); |
221 Heap heap(heap_cross_ref); |
222 Heap heap(heap_cross_ref); |
230 _dist[u] = heap.prio(); |
231 _dist[u] = heap.prio(); |
231 _proc_nodes.push_back(u); |
232 _proc_nodes.push_back(u); |
232 heap.pop(); |
233 heap.pop(); |
233 |
234 |
234 // Traverse outgoing residual arcs |
235 // Traverse outgoing residual arcs |
235 for (int a = _first_out[u]; a != _first_out[u+1]; ++a) { |
236 int last_out = _geq ? _first_out[u+1] : _first_out[u+1] - 1; |
|
237 for (int a = _first_out[u]; a != last_out; ++a) { |
236 if (_res_cap[a] < delta) continue; |
238 if (_res_cap[a] < delta) continue; |
237 v = _target[a]; |
239 v = _target[a]; |
238 switch (heap.state(v)) { |
240 switch (heap.state(v)) { |
239 case Heap::PRE_HEAP: |
241 case Heap::PRE_HEAP: |
240 heap.push(v, d + _cost[a] - _pi[v]); |
242 heap.push(v, d + _cost[a] - _pi[v]); |
685 for (int i = 0; i != _root; ++i) { |
687 for (int i = 0; i != _root; ++i) { |
686 _sum_supply += _supply[i]; |
688 _sum_supply += _supply[i]; |
687 } |
689 } |
688 if (_sum_supply > 0) return INFEASIBLE; |
690 if (_sum_supply > 0) return INFEASIBLE; |
689 |
691 |
690 // Initialize maps |
692 // Initialize vectors |
691 for (int i = 0; i != _root; ++i) { |
693 for (int i = 0; i != _root; ++i) { |
692 _pi[i] = 0; |
694 _pi[i] = 0; |
693 _excess[i] = _supply[i]; |
695 _excess[i] = _supply[i]; |
694 } |
696 } |
695 |
697 |
696 // Remove non-zero lower bounds |
698 // Remove non-zero lower bounds |
|
699 const Value MAX = std::numeric_limits<Value>::max(); |
|
700 int last_out; |
697 if (_have_lower) { |
701 if (_have_lower) { |
698 for (int i = 0; i != _root; ++i) { |
702 for (int i = 0; i != _root; ++i) { |
699 for (int j = _first_out[i]; j != _first_out[i+1]; ++j) { |
703 last_out = _first_out[i+1]; |
|
704 for (int j = _first_out[i]; j != last_out; ++j) { |
700 if (_forward[j]) { |
705 if (_forward[j]) { |
701 Value c = _lower[j]; |
706 Value c = _lower[j]; |
702 if (c >= 0) { |
707 if (c >= 0) { |
703 _res_cap[j] = _upper[j] < INF ? _upper[j] - c : INF; |
708 _res_cap[j] = _upper[j] < MAX ? _upper[j] - c : INF; |
704 } else { |
709 } else { |
705 _res_cap[j] = _upper[j] < INF + c ? _upper[j] - c : INF; |
710 _res_cap[j] = _upper[j] < MAX + c ? _upper[j] - c : INF; |
706 } |
711 } |
707 _excess[i] -= c; |
712 _excess[i] -= c; |
708 _excess[_target[j]] += c; |
713 _excess[_target[j]] += c; |
709 } else { |
714 } else { |
710 _res_cap[j] = 0; |
715 _res_cap[j] = 0; |
716 _res_cap[j] = _forward[j] ? _upper[j] : 0; |
721 _res_cap[j] = _forward[j] ? _upper[j] : 0; |
717 } |
722 } |
718 } |
723 } |
719 |
724 |
720 // Handle negative costs |
725 // Handle negative costs |
721 for (int u = 0; u != _root; ++u) { |
726 for (int i = 0; i != _root; ++i) { |
722 for (int a = _first_out[u]; a != _first_out[u+1]; ++a) { |
727 last_out = _first_out[i+1] - 1; |
723 Value rc = _res_cap[a]; |
728 for (int j = _first_out[i]; j != last_out; ++j) { |
724 if (_cost[a] < 0 && rc > 0) { |
729 Value rc = _res_cap[j]; |
725 if (rc == INF) return UNBOUNDED; |
730 if (_cost[j] < 0 && rc > 0) { |
726 _excess[u] -= rc; |
731 if (rc >= MAX) return UNBOUNDED; |
727 _excess[_target[a]] += rc; |
732 _excess[i] -= rc; |
728 _res_cap[a] = 0; |
733 _excess[_target[j]] += rc; |
729 _res_cap[_reverse[a]] += rc; |
734 _res_cap[j] = 0; |
|
735 _res_cap[_reverse[j]] += rc; |
730 } |
736 } |
731 } |
737 } |
732 } |
738 } |
733 |
739 |
734 // Handle GEQ supply type |
740 // Handle GEQ supply type |
735 if (_sum_supply < 0) { |
741 if (_sum_supply < 0) { |
736 _pi[_root] = 0; |
742 _pi[_root] = 0; |
737 _excess[_root] = -_sum_supply; |
743 _excess[_root] = -_sum_supply; |
738 for (int a = _first_out[_root]; a != _res_arc_num; ++a) { |
744 for (int a = _first_out[_root]; a != _res_arc_num; ++a) { |
739 int u = _target[a]; |
745 int ra = _reverse[a]; |
740 if (_excess[u] < 0) { |
746 _res_cap[a] = -_sum_supply + 1; |
741 _res_cap[a] = -_excess[u] + 1; |
747 _res_cap[ra] = 0; |
742 } else { |
|
743 _res_cap[a] = 1; |
|
744 } |
|
745 _res_cap[_reverse[a]] = 0; |
|
746 _cost[a] = 0; |
748 _cost[a] = 0; |
747 _cost[_reverse[a]] = 0; |
749 _cost[ra] = 0; |
748 } |
750 } |
749 } else { |
751 } else { |
750 _pi[_root] = 0; |
752 _pi[_root] = 0; |
751 _excess[_root] = 0; |
753 _excess[_root] = 0; |
752 for (int a = _first_out[_root]; a != _res_arc_num; ++a) { |
754 for (int a = _first_out[_root]; a != _res_arc_num; ++a) { |
|
755 int ra = _reverse[a]; |
753 _res_cap[a] = 1; |
756 _res_cap[a] = 1; |
754 _res_cap[_reverse[a]] = 0; |
757 _res_cap[ra] = 0; |
755 _cost[a] = 0; |
758 _cost[a] = 0; |
756 _cost[_reverse[a]] = 0; |
759 _cost[ra] = 0; |
757 } |
760 } |
758 } |
761 } |
759 |
762 |
760 // Initialize delta value |
763 // Initialize delta value |
761 if (_factor > 1) { |
764 if (_factor > 1) { |
762 // With scaling |
765 // With scaling |
763 Value max_sup = 0, max_dem = 0; |
766 Value max_sup = 0, max_dem = 0; |
764 for (int i = 0; i != _node_num; ++i) { |
767 for (int i = 0; i != _node_num; ++i) { |
765 if ( _excess[i] > max_sup) max_sup = _excess[i]; |
768 Value ex = _excess[i]; |
766 if (-_excess[i] > max_dem) max_dem = -_excess[i]; |
769 if ( ex > max_sup) max_sup = ex; |
|
770 if (-ex > max_dem) max_dem = -ex; |
767 } |
771 } |
768 Value max_cap = 0; |
772 Value max_cap = 0; |
769 for (int j = 0; j != _res_arc_num; ++j) { |
773 for (int j = 0; j != _res_arc_num; ++j) { |
770 if (_res_cap[j] > max_cap) max_cap = _res_cap[j]; |
774 if (_res_cap[j] > max_cap) max_cap = _res_cap[j]; |
771 } |
775 } |
787 else |
791 else |
788 pt = startWithoutScaling(); |
792 pt = startWithoutScaling(); |
789 |
793 |
790 // Handle non-zero lower bounds |
794 // Handle non-zero lower bounds |
791 if (_have_lower) { |
795 if (_have_lower) { |
792 for (int j = 0; j != _res_arc_num - _node_num + 1; ++j) { |
796 int limit = _first_out[_root]; |
|
797 for (int j = 0; j != limit; ++j) { |
793 if (!_forward[j]) _res_cap[j] += _lower[j]; |
798 if (!_forward[j]) _res_cap[j] += _lower[j]; |
794 } |
799 } |
795 } |
800 } |
796 |
801 |
797 // Shift potentials if necessary |
802 // Shift potentials if necessary |
810 // Perform capacity scaling phases |
815 // Perform capacity scaling phases |
811 int s, t; |
816 int s, t; |
812 ResidualDijkstra _dijkstra(*this); |
817 ResidualDijkstra _dijkstra(*this); |
813 while (true) { |
818 while (true) { |
814 // Saturate all arcs not satisfying the optimality condition |
819 // Saturate all arcs not satisfying the optimality condition |
|
820 int last_out; |
815 for (int u = 0; u != _node_num; ++u) { |
821 for (int u = 0; u != _node_num; ++u) { |
816 for (int a = _first_out[u]; a != _first_out[u+1]; ++a) { |
822 last_out = _sum_supply < 0 ? |
|
823 _first_out[u+1] : _first_out[u+1] - 1; |
|
824 for (int a = _first_out[u]; a != last_out; ++a) { |
817 int v = _target[a]; |
825 int v = _target[a]; |
818 Cost c = _cost[a] + _pi[u] - _pi[v]; |
826 Cost c = _cost[a] + _pi[u] - _pi[v]; |
819 Value rc = _res_cap[a]; |
827 Value rc = _res_cap[a]; |
820 if (c < 0 && rc >= _delta) { |
828 if (c < 0 && rc >= _delta) { |
821 _excess[u] -= rc; |
829 _excess[u] -= rc; |
828 |
836 |
829 // Find excess nodes and deficit nodes |
837 // Find excess nodes and deficit nodes |
830 _excess_nodes.clear(); |
838 _excess_nodes.clear(); |
831 _deficit_nodes.clear(); |
839 _deficit_nodes.clear(); |
832 for (int u = 0; u != _node_num; ++u) { |
840 for (int u = 0; u != _node_num; ++u) { |
833 if (_excess[u] >= _delta) _excess_nodes.push_back(u); |
841 Value ex = _excess[u]; |
834 if (_excess[u] <= -_delta) _deficit_nodes.push_back(u); |
842 if (ex >= _delta) _excess_nodes.push_back(u); |
|
843 if (ex <= -_delta) _deficit_nodes.push_back(u); |
835 } |
844 } |
836 int next_node = 0, next_def_node = 0; |
845 int next_node = 0, next_def_node = 0; |
837 |
846 |
838 // Find augmenting shortest paths |
847 // Find augmenting shortest paths |
839 while (next_node < int(_excess_nodes.size())) { |
848 while (next_node < int(_excess_nodes.size())) { |