gravatar
kpeter (Peter Kovacs)
kpeter@inf.elte.hu
Fix the GEQ/LEQ handling in NetworkSimplex + improve doc (#291) - Fix the optimality conditions for the GEQ/LEQ form. - Fix the initialization of the algortihm. It ensures correct solutions and it is much faster for the inequality forms. - Fix the pivot rules to search all the arcs that have to be allowed to get in the basis. - Better block size for the Block Search pivot rule. - Improve documentation of the problem and move it to a separate page.
0 3 1
default
4 files changed with 338 insertions and 188 deletions:
↑ Collapse diff ↑
Ignore white space 6 line context
1
/* -*- mode: C++; indent-tabs-mode: nil; -*-
2
 *
3
 * This file is a part of LEMON, a generic C++ optimization library.
4
 *
5
 * Copyright (C) 2003-2009
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
namespace lemon {
20

	
21
/**
22
\page min_cost_flow Minimum Cost Flow Problem
23

	
24
\section mcf_def Definition (GEQ form)
25

	
26
The \e minimum \e cost \e flow \e problem is to find a feasible flow of
27
minimum total cost from a set of supply nodes to a set of demand nodes
28
in a network with capacity constraints (lower and upper bounds)
29
and arc costs.
30

	
31
Formally, let \f$G=(V,A)\f$ be a digraph, \f$lower: A\rightarrow\mathbf{R}\f$,
32
\f$upper: A\rightarrow\mathbf{R}\cup\{+\infty\}\f$ denote the lower and
33
upper bounds for the flow values on the arcs, for which
34
\f$lower(uv) \leq upper(uv)\f$ must hold for all \f$uv\in A\f$,
35
\f$cost: A\rightarrow\mathbf{R}\f$ denotes the cost per unit flow
36
on the arcs and \f$sup: V\rightarrow\mathbf{R}\f$ denotes the
37
signed supply values of the nodes.
38
If \f$sup(u)>0\f$, then \f$u\f$ is a supply node with \f$sup(u)\f$
39
supply, if \f$sup(u)<0\f$, then \f$u\f$ is a demand node with
40
\f$-sup(u)\f$ demand.
41
A minimum cost flow is an \f$f: A\rightarrow\mathbf{R}\f$ solution
42
of the following optimization problem.
43

	
44
\f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
45
\f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \geq
46
    sup(u) \quad \forall u\in V \f]
47
\f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
48

	
49
The sum of the supply values, i.e. \f$\sum_{u\in V} sup(u)\f$ must be
50
zero or negative in order to have a feasible solution (since the sum
51
of the expressions on the left-hand side of the inequalities is zero).
52
It means that the total demand must be greater or equal to the total
53
supply and all the supplies have to be carried out from the supply nodes,
54
but there could be demands that are not satisfied.
55
If \f$\sum_{u\in V} sup(u)\f$ is zero, then all the supply/demand
56
constraints have to be satisfied with equality, i.e. all demands
57
have to be satisfied and all supplies have to be used.
58

	
59

	
60
\section mcf_algs Algorithms
61

	
62
LEMON contains several algorithms for solving this problem, for more
63
information see \ref min_cost_flow_algs "Minimum Cost Flow Algorithms".
64

	
65
A feasible solution for this problem can be found using \ref Circulation.
66

	
67

	
68
\section mcf_dual Dual Solution
69

	
70
The dual solution of the minimum cost flow problem is represented by
71
node potentials \f$\pi: V\rightarrow\mathbf{R}\f$.
72
An \f$f: A\rightarrow\mathbf{R}\f$ primal feasible solution is optimal
73
if and only if for some \f$\pi: V\rightarrow\mathbf{R}\f$ node potentials
74
the following \e complementary \e slackness optimality conditions hold.
75

	
76
 - For all \f$uv\in A\f$ arcs:
77
   - if \f$cost^\pi(uv)>0\f$, then \f$f(uv)=lower(uv)\f$;
78
   - if \f$lower(uv)<f(uv)<upper(uv)\f$, then \f$cost^\pi(uv)=0\f$;
79
   - if \f$cost^\pi(uv)<0\f$, then \f$f(uv)=upper(uv)\f$.
80
 - For all \f$u\in V\f$ nodes:
81
   - \f$\pi(u)<=0\f$;
82
   - if \f$\sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \neq sup(u)\f$,
83
     then \f$\pi(u)=0\f$.
84
 
85
Here \f$cost^\pi(uv)\f$ denotes the \e reduced \e cost of the arc
86
\f$uv\in A\f$ with respect to the potential function \f$\pi\f$, i.e.
87
\f[ cost^\pi(uv) = cost(uv) + \pi(u) - \pi(v).\f]
88

	
89
All algorithms provide dual solution (node potentials), as well,
90
if an optimal flow is found.
91

	
92

	
93
\section mcf_eq Equality Form
94

	
95
The above \ref mcf_def "definition" is actually more general than the
96
usual formulation of the minimum cost flow problem, in which strict
97
equalities are required in the supply/demand contraints.
98

	
99
\f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
100
\f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) =
101
    sup(u) \quad \forall u\in V \f]
102
\f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
103

	
104
However if the sum of the supply values is zero, then these two problems
105
are equivalent.
106
The \ref min_cost_flow_algs "algorithms" in LEMON support the general
107
form, so if you need the equality form, you have to ensure this additional
108
contraint manually.
109

	
110

	
111
\section mcf_leq Opposite Inequalites (LEQ Form)
112

	
113
Another possible definition of the minimum cost flow problem is
114
when there are <em>"less or equal"</em> (LEQ) supply/demand constraints,
115
instead of the <em>"greater or equal"</em> (GEQ) constraints.
116

	
117
\f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
118
\f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \leq
119
    sup(u) \quad \forall u\in V \f]
120
\f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
121

	
122
It means that the total demand must be less or equal to the 
123
total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
124
positive) and all the demands have to be satisfied, but there
125
could be supplies that are not carried out from the supply
126
nodes.
127
The equality form is also a special case of this form, of course.
128

	
129
You could easily transform this case to the \ref mcf_def "GEQ form"
130
of the problem by reversing the direction of the arcs and taking the
131
negative of the supply values (e.g. using \ref ReverseDigraph and
132
\ref NegMap adaptors).
133
However \ref NetworkSimplex algorithm also supports this form directly
134
for the sake of convenience.
135

	
136
Note that the optimality conditions for this supply constraint type are
137
slightly differ from the conditions that are discussed for the GEQ form,
138
namely the potentials have to be non-negative instead of non-positive.
139
An \f$f: A\rightarrow\mathbf{R}\f$ feasible solution of this problem
140
is optimal if and only if for some \f$\pi: V\rightarrow\mathbf{R}\f$
141
node potentials the following conditions hold.
142

	
143
 - For all \f$uv\in A\f$ arcs:
144
   - if \f$cost^\pi(uv)>0\f$, then \f$f(uv)=lower(uv)\f$;
145
   - if \f$lower(uv)<f(uv)<upper(uv)\f$, then \f$cost^\pi(uv)=0\f$;
146
   - if \f$cost^\pi(uv)<0\f$, then \f$f(uv)=upper(uv)\f$.
147
 - For all \f$u\in V\f$ nodes:
148
   - \f$\pi(u)>=0\f$;
149
   - if \f$\sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \neq sup(u)\f$,
150
     then \f$\pi(u)=0\f$.
151

	
152
*/
153
}
Ignore white space 6 line context
... ...
@@ -10,2 +10,3 @@
10 10
	doc/migration.dox \
11
	doc/min_cost_flow.dox \
11 12
	doc/named-param.dox \
Ignore white space 2 line context
... ...
@@ -337,3 +337,3 @@
337 337
/**
338
@defgroup min_cost_flow Minimum Cost Flow Algorithms
338
@defgroup min_cost_flow_algs Minimum Cost Flow Algorithms
339 339
@ingroup algs
... ...
@@ -343,81 +343,6 @@
343 343
This group contains the algorithms for finding minimum cost flows and
344
circulations.
344
circulations. For more information about this problem and its dual
345
solution see \ref min_cost_flow "Minimum Cost Flow Problem".
345 346

	
346
The \e minimum \e cost \e flow \e problem is to find a feasible flow of
347
minimum total cost from a set of supply nodes to a set of demand nodes
348
in a network with capacity constraints (lower and upper bounds)
349
and arc costs.
350
Formally, let \f$G=(V,A)\f$ be a digraph, \f$lower: A\rightarrow\mathbf{Z}\f$,
351
\f$upper: A\rightarrow\mathbf{Z}\cup\{+\infty\}\f$ denote the lower and
352
upper bounds for the flow values on the arcs, for which
353
\f$lower(uv) \leq upper(uv)\f$ must hold for all \f$uv\in A\f$,
354
\f$cost: A\rightarrow\mathbf{Z}\f$ denotes the cost per unit flow
355
on the arcs and \f$sup: V\rightarrow\mathbf{Z}\f$ denotes the
356
signed supply values of the nodes.
357
If \f$sup(u)>0\f$, then \f$u\f$ is a supply node with \f$sup(u)\f$
358
supply, if \f$sup(u)<0\f$, then \f$u\f$ is a demand node with
359
\f$-sup(u)\f$ demand.
360
A minimum cost flow is an \f$f: A\rightarrow\mathbf{Z}\f$ solution
361
of the following optimization problem.
362

	
363
\f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
364
\f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \geq
365
    sup(u) \quad \forall u\in V \f]
366
\f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
367

	
368
The sum of the supply values, i.e. \f$\sum_{u\in V} sup(u)\f$ must be
369
zero or negative in order to have a feasible solution (since the sum
370
of the expressions on the left-hand side of the inequalities is zero).
371
It means that the total demand must be greater or equal to the total
372
supply and all the supplies have to be carried out from the supply nodes,
373
but there could be demands that are not satisfied.
374
If \f$\sum_{u\in V} sup(u)\f$ is zero, then all the supply/demand
375
constraints have to be satisfied with equality, i.e. all demands
376
have to be satisfied and all supplies have to be used.
377

	
378
If you need the opposite inequalities in the supply/demand constraints
379
(i.e. the total demand is less than the total supply and all the demands
380
have to be satisfied while there could be supplies that are not used),
381
then you could easily transform the problem to the above form by reversing
382
the direction of the arcs and taking the negative of the supply values
383
(e.g. using \ref ReverseDigraph and \ref NegMap adaptors).
384
However \ref NetworkSimplex algorithm also supports this form directly
385
for the sake of convenience.
386

	
387
A feasible solution for this problem can be found using \ref Circulation.
388

	
389
Note that the above formulation is actually more general than the usual
390
definition of the minimum cost flow problem, in which strict equalities
391
are required in the supply/demand contraints, i.e.
392

	
393
\f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) =
394
    sup(u) \quad \forall u\in V. \f]
395

	
396
However if the sum of the supply values is zero, then these two problems
397
are equivalent. So if you need the equality form, you have to ensure this
398
additional contraint for the algorithms.
399

	
400
The dual solution of the minimum cost flow problem is represented by node 
401
potentials \f$\pi: V\rightarrow\mathbf{Z}\f$.
402
An \f$f: A\rightarrow\mathbf{Z}\f$ feasible solution of the problem
403
is optimal if and only if for some \f$\pi: V\rightarrow\mathbf{Z}\f$
404
node potentials the following \e complementary \e slackness optimality
405
conditions hold.
406

	
407
 - For all \f$uv\in A\f$ arcs:
408
   - if \f$cost^\pi(uv)>0\f$, then \f$f(uv)=lower(uv)\f$;
409
   - if \f$lower(uv)<f(uv)<upper(uv)\f$, then \f$cost^\pi(uv)=0\f$;
410
   - if \f$cost^\pi(uv)<0\f$, then \f$f(uv)=upper(uv)\f$.
411
 - For all \f$u\in V\f$ nodes:
412
   - if \f$\sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \neq sup(u)\f$,
413
     then \f$\pi(u)=0\f$.
414
 
415
Here \f$cost^\pi(uv)\f$ denotes the \e reduced \e cost of the arc
416
\f$uv\in A\f$ with respect to the potential function \f$\pi\f$, i.e.
417
\f[ cost^\pi(uv) = cost(uv) + \pi(u) - \pi(v).\f]
418

	
419
All algorithms provide dual solution (node potentials) as well,
420
if an optimal flow is found.
421

	
422
LEMON contains several algorithms for solving minimum cost flow problems.
347
LEMON contains several algorithms for this problem.
423 348
 - \ref NetworkSimplex Primal Network Simplex algorithm with various
... ...
@@ -431,6 +356,2 @@
431 356

	
432
Most of these implementations support the general inequality form of the
433
minimum cost flow problem, but CancelAndTighten and CycleCanceling
434
only support the equality form due to the primal method they use.
435

	
436 357
In general NetworkSimplex is the most efficient implementation,
Ignore white space 6 line context
... ...
@@ -21,3 +21,3 @@
21 21

	
22
/// \ingroup min_cost_flow
22
/// \ingroup min_cost_flow_algs
23 23
///
... ...
@@ -35,3 +35,3 @@
35 35

	
36
  /// \addtogroup min_cost_flow
36
  /// \addtogroup min_cost_flow_algs
37 37
  /// @{
... ...
@@ -104,46 +104,12 @@
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
    };
... ...
@@ -217,2 +183,4 @@
217 183
    int _arc_num;
184
    int _all_arc_num;
185
    int _search_arc_num;
218 186

	
... ...
@@ -279,3 +247,3 @@
279 247
      int &_in_arc;
280
      int _arc_num;
248
      int _search_arc_num;
281 249

	
... ...
@@ -290,3 +258,4 @@
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
      {}
... ...
@@ -296,3 +265,3 @@
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]]);
... ...
@@ -330,3 +299,3 @@
330 299
      int &_in_arc;
331
      int _arc_num;
300
      int _search_arc_num;
332 301

	
... ...
@@ -338,3 +307,3 @@
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
      {}
... ...
@@ -344,3 +313,3 @@
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]]);
... ...
@@ -369,3 +338,3 @@
369 338
      int &_in_arc;
370
      int _arc_num;
339
      int _search_arc_num;
371 340

	
... ...
@@ -381,6 +350,7 @@
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;
... ...
@@ -388,3 +358,3 @@
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 );
... ...
@@ -397,3 +367,3 @@
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]]);
... ...
@@ -442,3 +412,3 @@
442 412
      int &_in_arc;
443
      int _arc_num;
413
      int _search_arc_num;
444 414

	
... ...
@@ -456,3 +426,4 @@
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
      {
... ...
@@ -465,3 +436,3 @@
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 );
... ...
@@ -502,3 +473,3 @@
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]]);
... ...
@@ -548,3 +519,3 @@
548 519
      int &_in_arc;
549
      int _arc_num;
520
      int _search_arc_num;
550 521

	
... ...
@@ -576,4 +547,4 @@
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
      {
... ...
@@ -586,3 +557,3 @@
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 );
... ...
@@ -612,3 +583,3 @@
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] *
... ...
@@ -680,13 +651,13 @@
680 651
      int all_node_num = _node_num + 1;
681
      int all_arc_num = _arc_num + _node_num;
652
      int max_arc_num = _arc_num + 2 * _node_num;
682 653

	
683
      _source.resize(all_arc_num);
684
      _target.resize(all_arc_num);
654
      _source.resize(max_arc_num);
655
      _target.resize(max_arc_num);
685 656

	
686
      _lower.resize(all_arc_num);
687
      _upper.resize(all_arc_num);
688
      _cap.resize(all_arc_num);
689
      _cost.resize(all_arc_num);
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);
... ...
@@ -700,3 +671,3 @@
700 671
      _last_succ.resize(all_node_num);
701
      _state.resize(all_arc_num);
672
      _state.resize(max_arc_num);
702 673

	
... ...
@@ -1071,3 +1042,3 @@
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 {
... ...
@@ -1095,25 +1066,117 @@
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];
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
          }
1117 1098
        }
1118 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;
1181
      }
1119 1182

	
... ...
@@ -1376,16 +1439,4 @@
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
      }
... ...
@@ -1403,2 +1454,26 @@
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

	
0 comments (0 inline)