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 12 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
... ...
@@ -5,12 +5,13 @@
5 5
	doc/dirs.dox \
6 6
	doc/groups.dox \
7 7
	doc/lgf.dox \
8 8
	doc/license.dox \
9 9
	doc/mainpage.dox \
10 10
	doc/migration.dox \
11
	doc/min_cost_flow.dox \
11 12
	doc/named-param.dox \
12 13
	doc/namespaces.dox \
13 14
	doc/html \
14 15
	doc/CMakeLists.txt
15 16

	
16 17
DOC_EPS_IMAGES18 = \
Ignore white space 6 line context
... ...
@@ -332,110 +332,31 @@
332 332
for finding feasible circulations, which is a somewhat different problem,
333 333
but it is strongly related to maximum flow.
334 334
For more information, see \ref Circulation.
335 335
*/
336 336

	
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
340 340

	
341 341
\brief Algorithms for finding minimum cost flows and circulations.
342 342

	
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
424 349
   pivot strategies.
425 350
 - \ref CostScaling Push-Relabel and Augment-Relabel algorithms based on
426 351
   cost scaling.
427 352
 - \ref CapacityScaling Successive Shortest %Path algorithm with optional
428 353
   capacity scaling.
429 354
 - \ref CancelAndTighten The Cancel and Tighten algorithm.
430 355
 - \ref CycleCanceling Cycle-Canceling algorithms.
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,
437 358
but in special cases other algorithms could be faster.
438 359
For example, if the total supply and/or capacities are rather small,
439 360
CapacityScaling is usually the fastest algorithm (without effective scaling).
440 361
*/
441 362

	
Ignore white space 6 line context
... ...
@@ -16,13 +16,13 @@
16 16
 *
17 17
 */
18 18

	
19 19
#ifndef LEMON_NETWORK_SIMPLEX_H
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
25 25
/// \brief Network Simplex algorithm for finding a minimum cost flow.
26 26

	
27 27
#include <vector>
28 28
#include <limits>
... ...
@@ -30,13 +30,13 @@
30 30

	
31 31
#include <lemon/core.h>
32 32
#include <lemon/math.h>
33 33

	
34 34
namespace lemon {
35 35

	
36
  /// \addtogroup min_cost_flow
36
  /// \addtogroup min_cost_flow_algs
37 37
  /// @{
38 38

	
39 39
  /// \brief Implementation of the primal Network Simplex algorithm
40 40
  /// for finding a \ref min_cost_flow "minimum cost flow".
41 41
  ///
42 42
  /// \ref NetworkSimplex implements the primal Network Simplex algorithm
... ...
@@ -99,56 +99,22 @@
99 99
    /// \brief Constants for selecting the type of the supply constraints.
100 100
    ///
101 101
    /// Enum type containing constants for selecting the supply type,
102 102
    /// i.e. the direction of the inequalities in the supply/demand
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
    
151 117
    /// \brief Constants for selecting the pivot rule.
152 118
    ///
153 119
    /// Enum type containing constants for selecting the pivot rule for
154 120
    /// the \ref run() function.
... ...
@@ -212,12 +178,14 @@
212 178
  private:
213 179

	
214 180
    // Data related to the underlying digraph
215 181
    const GR &_graph;
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
220 188
    bool _have_lower;
221 189
    SupplyType _stype;
222 190
    Value _sum_supply;
223 191

	
... ...
@@ -274,30 +242,31 @@
274 242
      const IntVector  &_source;
275 243
      const IntVector  &_target;
276 244
      const CostVector &_cost;
277 245
      const IntVector  &_state;
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
283 251
      int _next_arc;
284 252

	
285 253
    public:
286 254

	
287 255
      // Constructor
288 256
      FirstEligiblePivotRule(NetworkSimplex &ns) :
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

	
294 263
      // Find next entering arc
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) {
300 269
            _in_arc = e;
301 270
            _next_arc = e + 1;
302 271
            return true;
303 272
          }
... ...
@@ -325,27 +294,27 @@
325 294
      const IntVector  &_source;
326 295
      const IntVector  &_target;
327 296
      const CostVector &_cost;
328 297
      const IntVector  &_state;
329 298
      const CostVector &_pi;
330 299
      int &_in_arc;
331
      int _arc_num;
300
      int _search_arc_num;
332 301

	
333 302
    public:
334 303

	
335 304
      // Constructor
336 305
      BestEligiblePivotRule(NetworkSimplex &ns) :
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

	
342 311
      // Find next entering arc
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) {
348 317
            min = c;
349 318
            _in_arc = e;
350 319
          }
351 320
        }
... ...
@@ -364,41 +333,42 @@
364 333
      const IntVector  &_source;
365 334
      const IntVector  &_target;
366 335
      const CostVector &_cost;
367 336
      const IntVector  &_state;
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
373 342
      int _block_size;
374 343
      int _next_arc;
375 344

	
376 345
    public:
377 346

	
378 347
      // Constructor
379 348
      BlockSearchPivotRule(NetworkSimplex &ns) :
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
      }
392 362

	
393 363
      // Find next entering arc
394 364
      bool findEnteringArc() {
395 365
        Cost c, min = 0;
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) {
401 371
            min = c;
402 372
            min_arc = e;
403 373
          }
404 374
          if (--cnt == 0) {
... ...
@@ -437,13 +407,13 @@
437 407
      const IntVector  &_source;
438 408
      const IntVector  &_target;
439 409
      const CostVector &_cost;
440 410
      const IntVector  &_state;
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
446 416
      IntVector _candidates;
447 417
      int _list_length, _minor_limit;
448 418
      int _curr_length, _minor_count;
449 419
      int _next_arc;
... ...
@@ -451,22 +421,23 @@
451 421
    public:
452 422

	
453 423
      /// Constructor
454 424
      CandidateListPivotRule(NetworkSimplex &ns) :
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
460 431
        const double LIST_LENGTH_FACTOR = 1.0;
461 432
        const int MIN_LIST_LENGTH = 10;
462 433
        const double MINOR_LIMIT_FACTOR = 0.1;
463 434
        const int MIN_MINOR_LIMIT = 3;
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),
469 440
                                 MIN_MINOR_LIMIT );
470 441
        _curr_length = _minor_count = 0;
471 442
        _candidates.resize(_list_length);
472 443
      }
... ...
@@ -497,13 +468,13 @@
497 468
          }
498 469
        }
499 470

	
500 471
        // Major iteration: build a new candidate list
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) {
506 477
            _candidates[_curr_length++] = e;
507 478
            if (c < min) {
508 479
              min = c;
509 480
              min_arc = e;
... ...
@@ -543,13 +514,13 @@
543 514
      const IntVector  &_source;
544 515
      const IntVector  &_target;
545 516
      const CostVector &_cost;
546 517
      const IntVector  &_state;
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
552 523
      int _block_size, _head_length, _curr_length;
553 524
      int _next_arc;
554 525
      IntVector _candidates;
555 526
      CostVector _cand_cost;
... ...
@@ -571,23 +542,23 @@
571 542
    public:
572 543

	
573 544
      // Constructor
574 545
      AlteringListPivotRule(NetworkSimplex &ns) :
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
581 552
        const double BLOCK_SIZE_FACTOR = 1.5;
582 553
        const int MIN_BLOCK_SIZE = 10;
583 554
        const double HEAD_LENGTH_FACTOR = 0.1;
584 555
        const int MIN_HEAD_LENGTH = 3;
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),
590 561
                                 MIN_HEAD_LENGTH );
591 562
        _candidates.resize(_head_length + _block_size);
592 563
        _curr_length = 0;
593 564
      }
... ...
@@ -607,13 +578,13 @@
607 578

	
608 579
        // Extend the list
609 580
        int cnt = _block_size;
610 581
        int last_arc = 0;
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]]);
616 587
          if (_cand_cost[e] < 0) {
617 588
            _candidates[_curr_length++] = e;
618 589
            last_arc = e;
619 590
          }
... ...
@@ -675,33 +646,33 @@
675 646
        "The cost type of NetworkSimplex must be signed");
676 647
        
677 648
      // Resize vectors
678 649
      _node_num = countNodes(_graph);
679 650
      _arc_num = countArcs(_graph);
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);
693 664

	
694 665
      _parent.resize(all_node_num);
695 666
      _pred.resize(all_node_num);
696 667
      _forward.resize(all_node_num);
697 668
      _thread.resize(all_node_num);
698 669
      _rev_thread.resize(all_node_num);
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)
704 675
      int i = 0;
705 676
      for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
706 677
        _node_id[n] = i;
707 678
      }
... ...
@@ -1066,13 +1037,13 @@
1066 1037
        }
1067 1038
      }
1068 1039

	
1069 1040
      // Initialize artifical cost
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();
1075 1046
        for (int i = 0; i != _arc_num; ++i) {
1076 1047
          if (_cost[i] > ART_COST) ART_COST = _cost[i];
1077 1048
        }
1078 1049
        ART_COST = (ART_COST + 1) * _node_num;
... ...
@@ -1090,35 +1061,127 @@
1090 1061
      _pred[_root] = -1;
1091 1062
      _thread[_root] = 0;
1092 1063
      _rev_thread[0] = _root;
1093 1064
      _succ_num[_root] = _node_num + 1;
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];
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

	
1120 1183
      return true;
1121 1184
    }
1122 1185

	
1123 1186
    // Find the join node
1124 1187
    void findJoinNode() {
... ...
@@ -1371,26 +1434,14 @@
1371 1434
          updateTreeStructure();
1372 1435
          updatePotential();
1373 1436
        }
1374 1437
      }
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

	
1393 1444
      // Transform the solution and the supply map to the original form
1394 1445
      if (_have_lower) {
1395 1446
        for (int i = 0; i != _arc_num; ++i) {
1396 1447
          Value c = _lower[i];
... ...
@@ -1398,12 +1449,36 @@
1398 1449
            _flow[i] += c;
1399 1450
            _supply[_source[i]] += c;
1400 1451
            _supply[_target[i]] -= c;
1401 1452
          }
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;
1406 1481
    }
1407 1482

	
1408 1483
  }; //class NetworkSimplex
1409 1484

	
0 comments (0 inline)