Changeset 710:8b0df68370a4 in lemon for lemon
 Timestamp:
 05/12/09 12:06:40 (13 years ago)
 Branch:
 default
 Phase:
 public
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

lemon/network_simplex.h
r690 r710 20 20 #define LEMON_NETWORK_SIMPLEX_H 21 21 22 /// \ingroup min_cost_flow 22 /// \ingroup min_cost_flow_algs 23 23 /// 24 24 /// \file … … 34 34 namespace lemon { 35 35 36 /// \addtogroup min_cost_flow 36 /// \addtogroup min_cost_flow_algs 37 37 /// @{ 38 38 … … 103 103 /// constraints of the \ref min_cost_flow "minimum cost flow problem". 104 104 /// 105 /// The default supply type is \c GEQ, since this form is supported 106 /// by other minimum cost flow algorithms and the \ref Circulation 107 /// algorithm, as well. 108 /// The \c LEQ problem type can be selected using the \ref supplyType() 109 /// function. 110 /// 111 /// Note that the equality form is a special case of both supply types. 105 /// The default supply type is \c GEQ, the \c LEQ type can be 106 /// selected using \ref supplyType(). 107 /// The equality form is a special case of both supply types. 112 108 enum SupplyType { 113 114 109 /// This option means that there are <em>"greater or equal"</em> 115 /// supply/demand constraints in the definition, i.e. the exact 116 /// formulation of the problem is the following. 117 /** 118 \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f] 119 \f[ \sum_{uv\in A} f(uv)  \sum_{vu\in A} f(vu) \geq 120 sup(u) \quad \forall u\in V \f] 121 \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f] 122 */ 123 /// It means that the total demand must be greater or equal to the 124 /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or 125 /// negative) and all the supplies have to be carried out from 126 /// the supply nodes, but there could be demands that are not 127 /// satisfied. 110 /// supply/demand constraints in the definition of the problem. 128 111 GEQ, 129 /// It is just an alias for the \c GEQ option.130 CARRY_SUPPLIES = GEQ,131 132 112 /// This option means that there are <em>"less or equal"</em> 133 /// supply/demand constraints in the definition, i.e. the exact 134 /// formulation of the problem is the following. 135 /** 136 \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f] 137 \f[ \sum_{uv\in A} f(uv)  \sum_{vu\in A} f(vu) \leq 138 sup(u) \quad \forall u\in V \f] 139 \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f] 140 */ 141 /// It means that the total demand must be less or equal to the 142 /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or 143 /// positive) and all the demands have to be satisfied, but there 144 /// could be supplies that are not carried out from the supply 145 /// nodes. 146 LEQ, 147 /// It is just an alias for the \c LEQ option. 148 SATISFY_DEMANDS = LEQ 113 /// supply/demand constraints in the definition of the problem. 114 LEQ 149 115 }; 150 116 … … 216 182 int _node_num; 217 183 int _arc_num; 184 int _all_arc_num; 185 int _search_arc_num; 218 186 219 187 // Parameters of the problem … … 278 246 const CostVector &_pi; 279 247 int &_in_arc; 280 int _ arc_num;248 int _search_arc_num; 281 249 282 250 // Pivot rule data … … 289 257 _source(ns._source), _target(ns._target), 290 258 _cost(ns._cost), _state(ns._state), _pi(ns._pi), 291 _in_arc(ns.in_arc), _arc_num(ns._arc_num), _next_arc(0) 259 _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num), 260 _next_arc(0) 292 261 {} 293 262 … … 295 264 bool findEnteringArc() { 296 265 Cost c; 297 for (int e = _next_arc; e < _ arc_num; ++e) {266 for (int e = _next_arc; e < _search_arc_num; ++e) { 298 267 c = _state[e] * (_cost[e] + _pi[_source[e]]  _pi[_target[e]]); 299 268 if (c < 0) { … … 329 298 const CostVector &_pi; 330 299 int &_in_arc; 331 int _ arc_num;300 int _search_arc_num; 332 301 333 302 public: … … 337 306 _source(ns._source), _target(ns._target), 338 307 _cost(ns._cost), _state(ns._state), _pi(ns._pi), 339 _in_arc(ns.in_arc), _ arc_num(ns._arc_num)308 _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num) 340 309 {} 341 310 … … 343 312 bool findEnteringArc() { 344 313 Cost c, min = 0; 345 for (int e = 0; e < _ arc_num; ++e) {314 for (int e = 0; e < _search_arc_num; ++e) { 346 315 c = _state[e] * (_cost[e] + _pi[_source[e]]  _pi[_target[e]]); 347 316 if (c < min) { … … 368 337 const CostVector &_pi; 369 338 int &_in_arc; 370 int _ arc_num;339 int _search_arc_num; 371 340 372 341 // Pivot rule data … … 380 349 _source(ns._source), _target(ns._target), 381 350 _cost(ns._cost), _state(ns._state), _pi(ns._pi), 382 _in_arc(ns.in_arc), _arc_num(ns._arc_num), _next_arc(0) 351 _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num), 352 _next_arc(0) 383 353 { 384 354 // The main parameters of the pivot rule 385 const double BLOCK_SIZE_FACTOR = 2.0;355 const double BLOCK_SIZE_FACTOR = 0.5; 386 356 const int MIN_BLOCK_SIZE = 10; 387 357 388 358 _block_size = std::max( int(BLOCK_SIZE_FACTOR * 389 std::sqrt(double(_ arc_num))),359 std::sqrt(double(_search_arc_num))), 390 360 MIN_BLOCK_SIZE ); 391 361 } … … 396 366 int cnt = _block_size; 397 367 int e, min_arc = _next_arc; 398 for (e = _next_arc; e < _ arc_num; ++e) {368 for (e = _next_arc; e < _search_arc_num; ++e) { 399 369 c = _state[e] * (_cost[e] + _pi[_source[e]]  _pi[_target[e]]); 400 370 if (c < min) { … … 441 411 const CostVector &_pi; 442 412 int &_in_arc; 443 int _ arc_num;413 int _search_arc_num; 444 414 445 415 // Pivot rule data … … 455 425 _source(ns._source), _target(ns._target), 456 426 _cost(ns._cost), _state(ns._state), _pi(ns._pi), 457 _in_arc(ns.in_arc), _arc_num(ns._arc_num), _next_arc(0) 427 _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num), 428 _next_arc(0) 458 429 { 459 430 // The main parameters of the pivot rule … … 464 435 465 436 _list_length = std::max( int(LIST_LENGTH_FACTOR * 466 std::sqrt(double(_ arc_num))),437 std::sqrt(double(_search_arc_num))), 467 438 MIN_LIST_LENGTH ); 468 439 _minor_limit = std::max( int(MINOR_LIMIT_FACTOR * _list_length), … … 501 472 min = 0; 502 473 _curr_length = 0; 503 for (e = _next_arc; e < _ arc_num; ++e) {474 for (e = _next_arc; e < _search_arc_num; ++e) { 504 475 c = _state[e] * (_cost[e] + _pi[_source[e]]  _pi[_target[e]]); 505 476 if (c < 0) { … … 547 518 const CostVector &_pi; 548 519 int &_in_arc; 549 int _ arc_num;520 int _search_arc_num; 550 521 551 522 // Pivot rule data … … 575 546 _source(ns._source), _target(ns._target), 576 547 _cost(ns._cost), _state(ns._state), _pi(ns._pi), 577 _in_arc(ns.in_arc), _ arc_num(ns._arc_num),578 _next_arc(0), _cand_cost(ns._ arc_num), _sort_func(_cand_cost)548 _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num), 549 _next_arc(0), _cand_cost(ns._search_arc_num), _sort_func(_cand_cost) 579 550 { 580 551 // The main parameters of the pivot rule … … 585 556 586 557 _block_size = std::max( int(BLOCK_SIZE_FACTOR * 587 std::sqrt(double(_ arc_num))),558 std::sqrt(double(_search_arc_num))), 588 559 MIN_BLOCK_SIZE ); 589 560 _head_length = std::max( int(HEAD_LENGTH_FACTOR * _block_size), … … 611 582 int limit = _head_length; 612 583 613 for (int e = _next_arc; e < _ arc_num; ++e) {584 for (int e = _next_arc; e < _search_arc_num; ++e) { 614 585 _cand_cost[e] = _state[e] * 615 586 (_cost[e] + _pi[_source[e]]  _pi[_target[e]]); … … 679 650 _arc_num = countArcs(_graph); 680 651 int all_node_num = _node_num + 1; 681 int all_arc_num = _arc_num +_node_num;682 683 _source.resize( all_arc_num);684 _target.resize( all_arc_num);685 686 _lower.resize( all_arc_num);687 _upper.resize( all_arc_num);688 _cap.resize( all_arc_num);689 _cost.resize( all_arc_num);652 int max_arc_num = _arc_num + 2 * _node_num; 653 654 _source.resize(max_arc_num); 655 _target.resize(max_arc_num); 656 657 _lower.resize(_arc_num); 658 _upper.resize(_arc_num); 659 _cap.resize(max_arc_num); 660 _cost.resize(max_arc_num); 690 661 _supply.resize(all_node_num); 691 _flow.resize( all_arc_num);662 _flow.resize(max_arc_num); 692 663 _pi.resize(all_node_num); 693 664 … … 699 670 _succ_num.resize(all_node_num); 700 671 _last_succ.resize(all_node_num); 701 _state.resize( all_arc_num);672 _state.resize(max_arc_num); 702 673 703 674 // Copy the graph (store the arcs in a mixed order) … … 1070 1041 Cost ART_COST; 1071 1042 if (std::numeric_limits<Cost>::is_exact) { 1072 ART_COST = std::numeric_limits<Cost>::max() / 4+ 1;1043 ART_COST = std::numeric_limits<Cost>::max() / 2 + 1; 1073 1044 } else { 1074 1045 ART_COST = std::numeric_limits<Cost>::min(); … … 1094 1065 _last_succ[_root] = _root  1; 1095 1066 _supply[_root] = _sum_supply; 1096 _pi[_root] = _sum_supply < 0 ? ART_COST : ART_COST;1067 _pi[_root] = 0; 1097 1068 1098 1069 // Add artificial arcs and initialize the spanning tree data structure 1099 for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) { 1100 _parent[u] = _root; 1101 _pred[u] = e; 1102 _thread[u] = u + 1; 1103 _rev_thread[u + 1] = u; 1104 _succ_num[u] = 1; 1105 _last_succ[u] = u; 1106 _cost[e] = ART_COST; 1107 _cap[e] = INF; 1108 _state[e] = STATE_TREE; 1109 if (_supply[u] > 0  (_supply[u] == 0 && _sum_supply <= 0)) { 1110 _flow[e] = _supply[u]; 1111 _forward[u] = true; 1112 _pi[u] = ART_COST + _pi[_root]; 1113 } else { 1114 _flow[e] = _supply[u]; 1115 _forward[u] = false; 1116 _pi[u] = ART_COST + _pi[_root]; 1117 } 1070 if (_sum_supply == 0) { 1071 // EQ supply constraints 1072 _search_arc_num = _arc_num; 1073 _all_arc_num = _arc_num + _node_num; 1074 for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) { 1075 _parent[u] = _root; 1076 _pred[u] = e; 1077 _thread[u] = u + 1; 1078 _rev_thread[u + 1] = u; 1079 _succ_num[u] = 1; 1080 _last_succ[u] = u; 1081 _cap[e] = INF; 1082 _state[e] = STATE_TREE; 1083 if (_supply[u] >= 0) { 1084 _forward[u] = true; 1085 _pi[u] = 0; 1086 _source[e] = u; 1087 _target[e] = _root; 1088 _flow[e] = _supply[u]; 1089 _cost[e] = 0; 1090 } else { 1091 _forward[u] = false; 1092 _pi[u] = ART_COST; 1093 _source[e] = _root; 1094 _target[e] = u; 1095 _flow[e] = _supply[u]; 1096 _cost[e] = ART_COST; 1097 } 1098 } 1099 } 1100 else if (_sum_supply > 0) { 1101 // LEQ supply constraints 1102 _search_arc_num = _arc_num + _node_num; 1103 int f = _arc_num + _node_num; 1104 for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) { 1105 _parent[u] = _root; 1106 _thread[u] = u + 1; 1107 _rev_thread[u + 1] = u; 1108 _succ_num[u] = 1; 1109 _last_succ[u] = u; 1110 if (_supply[u] >= 0) { 1111 _forward[u] = true; 1112 _pi[u] = 0; 1113 _pred[u] = e; 1114 _source[e] = u; 1115 _target[e] = _root; 1116 _cap[e] = INF; 1117 _flow[e] = _supply[u]; 1118 _cost[e] = 0; 1119 _state[e] = STATE_TREE; 1120 } else { 1121 _forward[u] = false; 1122 _pi[u] = ART_COST; 1123 _pred[u] = f; 1124 _source[f] = _root; 1125 _target[f] = u; 1126 _cap[f] = INF; 1127 _flow[f] = _supply[u]; 1128 _cost[f] = ART_COST; 1129 _state[f] = STATE_TREE; 1130 _source[e] = u; 1131 _target[e] = _root; 1132 _cap[e] = INF; 1133 _flow[e] = 0; 1134 _cost[e] = 0; 1135 _state[e] = STATE_LOWER; 1136 ++f; 1137 } 1138 } 1139 _all_arc_num = f; 1140 } 1141 else { 1142 // GEQ supply constraints 1143 _search_arc_num = _arc_num + _node_num; 1144 int f = _arc_num + _node_num; 1145 for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) { 1146 _parent[u] = _root; 1147 _thread[u] = u + 1; 1148 _rev_thread[u + 1] = u; 1149 _succ_num[u] = 1; 1150 _last_succ[u] = u; 1151 if (_supply[u] <= 0) { 1152 _forward[u] = false; 1153 _pi[u] = 0; 1154 _pred[u] = e; 1155 _source[e] = _root; 1156 _target[e] = u; 1157 _cap[e] = INF; 1158 _flow[e] = _supply[u]; 1159 _cost[e] = 0; 1160 _state[e] = STATE_TREE; 1161 } else { 1162 _forward[u] = true; 1163 _pi[u] = ART_COST; 1164 _pred[u] = f; 1165 _source[f] = u; 1166 _target[f] = _root; 1167 _cap[f] = INF; 1168 _flow[f] = _supply[u]; 1169 _state[f] = STATE_TREE; 1170 _cost[f] = ART_COST; 1171 _source[e] = _root; 1172 _target[e] = u; 1173 _cap[e] = INF; 1174 _flow[e] = 0; 1175 _cost[e] = 0; 1176 _state[e] = STATE_LOWER; 1177 ++f; 1178 } 1179 } 1180 _all_arc_num = f; 1118 1181 } 1119 1182 … … 1375 1438 1376 1439 // Check feasibility 1377 if (_sum_supply < 0) { 1378 for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) { 1379 if (_supply[u] >= 0 && _flow[e] != 0) return INFEASIBLE; 1380 } 1381 } 1382 else if (_sum_supply > 0) { 1383 for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) { 1384 if (_supply[u] <= 0 && _flow[e] != 0) return INFEASIBLE; 1385 } 1386 } 1387 else { 1388 for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) { 1389 if (_flow[e] != 0) return INFEASIBLE; 1390 } 1440 for (int e = _search_arc_num; e != _all_arc_num; ++e) { 1441 if (_flow[e] != 0) return INFEASIBLE; 1391 1442 } 1392 1443 … … 1402 1453 } 1403 1454 } 1455 1456 // Shift potentials to meet the requirements of the GEQ/LEQ type 1457 // optimality conditions 1458 if (_sum_supply == 0) { 1459 if (_stype == GEQ) { 1460 Cost max_pot = std::numeric_limits<Cost>::min(); 1461 for (int i = 0; i != _node_num; ++i) { 1462 if (_pi[i] > max_pot) max_pot = _pi[i]; 1463 } 1464 if (max_pot > 0) { 1465 for (int i = 0; i != _node_num; ++i) 1466 _pi[i] = max_pot; 1467 } 1468 } else { 1469 Cost min_pot = std::numeric_limits<Cost>::max(); 1470 for (int i = 0; i != _node_num; ++i) { 1471 if (_pi[i] < min_pot) min_pot = _pi[i]; 1472 } 1473 if (min_pot < 0) { 1474 for (int i = 0; i != _node_num; ++i) 1475 _pi[i] = min_pot; 1476 } 1477 } 1478 } 1404 1479 1405 1480 return OPTIMAL;
Note: See TracChangeset
for help on using the changeset viewer.