1.1 --- a/doc/Makefile.am Mon May 11 16:38:21 2009 +0100
1.2 +++ b/doc/Makefile.am Tue May 12 12:06:40 2009 +0200
1.3 @@ -8,6 +8,7 @@
1.4 doc/license.dox \
1.5 doc/mainpage.dox \
1.6 doc/migration.dox \
1.7 + doc/min_cost_flow.dox \
1.8 doc/named-param.dox \
1.9 doc/namespaces.dox \
1.10 doc/html \
2.1 --- a/doc/groups.dox Mon May 11 16:38:21 2009 +0100
2.2 +++ b/doc/groups.dox Tue May 12 12:06:40 2009 +0200
2.3 @@ -335,91 +335,16 @@
2.4 */
2.5
2.6 /**
2.7 -@defgroup min_cost_flow Minimum Cost Flow Algorithms
2.8 +@defgroup min_cost_flow_algs Minimum Cost Flow Algorithms
2.9 @ingroup algs
2.10
2.11 \brief Algorithms for finding minimum cost flows and circulations.
2.12
2.13 This group contains the algorithms for finding minimum cost flows and
2.14 -circulations.
2.15 +circulations. For more information about this problem and its dual
2.16 +solution see \ref min_cost_flow "Minimum Cost Flow Problem".
2.17
2.18 -The \e minimum \e cost \e flow \e problem is to find a feasible flow of
2.19 -minimum total cost from a set of supply nodes to a set of demand nodes
2.20 -in a network with capacity constraints (lower and upper bounds)
2.21 -and arc costs.
2.22 -Formally, let \f$G=(V,A)\f$ be a digraph, \f$lower: A\rightarrow\mathbf{Z}\f$,
2.23 -\f$upper: A\rightarrow\mathbf{Z}\cup\{+\infty\}\f$ denote the lower and
2.24 -upper bounds for the flow values on the arcs, for which
2.25 -\f$lower(uv) \leq upper(uv)\f$ must hold for all \f$uv\in A\f$,
2.26 -\f$cost: A\rightarrow\mathbf{Z}\f$ denotes the cost per unit flow
2.27 -on the arcs and \f$sup: V\rightarrow\mathbf{Z}\f$ denotes the
2.28 -signed supply values of the nodes.
2.29 -If \f$sup(u)>0\f$, then \f$u\f$ is a supply node with \f$sup(u)\f$
2.30 -supply, if \f$sup(u)<0\f$, then \f$u\f$ is a demand node with
2.31 -\f$-sup(u)\f$ demand.
2.32 -A minimum cost flow is an \f$f: A\rightarrow\mathbf{Z}\f$ solution
2.33 -of the following optimization problem.
2.34 -
2.35 -\f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
2.36 -\f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \geq
2.37 - sup(u) \quad \forall u\in V \f]
2.38 -\f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
2.39 -
2.40 -The sum of the supply values, i.e. \f$\sum_{u\in V} sup(u)\f$ must be
2.41 -zero or negative in order to have a feasible solution (since the sum
2.42 -of the expressions on the left-hand side of the inequalities is zero).
2.43 -It means that the total demand must be greater or equal to the total
2.44 -supply and all the supplies have to be carried out from the supply nodes,
2.45 -but there could be demands that are not satisfied.
2.46 -If \f$\sum_{u\in V} sup(u)\f$ is zero, then all the supply/demand
2.47 -constraints have to be satisfied with equality, i.e. all demands
2.48 -have to be satisfied and all supplies have to be used.
2.49 -
2.50 -If you need the opposite inequalities in the supply/demand constraints
2.51 -(i.e. the total demand is less than the total supply and all the demands
2.52 -have to be satisfied while there could be supplies that are not used),
2.53 -then you could easily transform the problem to the above form by reversing
2.54 -the direction of the arcs and taking the negative of the supply values
2.55 -(e.g. using \ref ReverseDigraph and \ref NegMap adaptors).
2.56 -However \ref NetworkSimplex algorithm also supports this form directly
2.57 -for the sake of convenience.
2.58 -
2.59 -A feasible solution for this problem can be found using \ref Circulation.
2.60 -
2.61 -Note that the above formulation is actually more general than the usual
2.62 -definition of the minimum cost flow problem, in which strict equalities
2.63 -are required in the supply/demand contraints, i.e.
2.64 -
2.65 -\f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) =
2.66 - sup(u) \quad \forall u\in V. \f]
2.67 -
2.68 -However if the sum of the supply values is zero, then these two problems
2.69 -are equivalent. So if you need the equality form, you have to ensure this
2.70 -additional contraint for the algorithms.
2.71 -
2.72 -The dual solution of the minimum cost flow problem is represented by node
2.73 -potentials \f$\pi: V\rightarrow\mathbf{Z}\f$.
2.74 -An \f$f: A\rightarrow\mathbf{Z}\f$ feasible solution of the problem
2.75 -is optimal if and only if for some \f$\pi: V\rightarrow\mathbf{Z}\f$
2.76 -node potentials the following \e complementary \e slackness optimality
2.77 -conditions hold.
2.78 -
2.79 - - For all \f$uv\in A\f$ arcs:
2.80 - - if \f$cost^\pi(uv)>0\f$, then \f$f(uv)=lower(uv)\f$;
2.81 - - if \f$lower(uv)<f(uv)<upper(uv)\f$, then \f$cost^\pi(uv)=0\f$;
2.82 - - if \f$cost^\pi(uv)<0\f$, then \f$f(uv)=upper(uv)\f$.
2.83 - - For all \f$u\in V\f$ nodes:
2.84 - - if \f$\sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \neq sup(u)\f$,
2.85 - then \f$\pi(u)=0\f$.
2.86 -
2.87 -Here \f$cost^\pi(uv)\f$ denotes the \e reduced \e cost of the arc
2.88 -\f$uv\in A\f$ with respect to the potential function \f$\pi\f$, i.e.
2.89 -\f[ cost^\pi(uv) = cost(uv) + \pi(u) - \pi(v).\f]
2.90 -
2.91 -All algorithms provide dual solution (node potentials) as well,
2.92 -if an optimal flow is found.
2.93 -
2.94 -LEMON contains several algorithms for solving minimum cost flow problems.
2.95 +LEMON contains several algorithms for this problem.
2.96 - \ref NetworkSimplex Primal Network Simplex algorithm with various
2.97 pivot strategies.
2.98 - \ref CostScaling Push-Relabel and Augment-Relabel algorithms based on
2.99 @@ -429,10 +354,6 @@
2.100 - \ref CancelAndTighten The Cancel and Tighten algorithm.
2.101 - \ref CycleCanceling Cycle-Canceling algorithms.
2.102
2.103 -Most of these implementations support the general inequality form of the
2.104 -minimum cost flow problem, but CancelAndTighten and CycleCanceling
2.105 -only support the equality form due to the primal method they use.
2.106 -
2.107 In general NetworkSimplex is the most efficient implementation,
2.108 but in special cases other algorithms could be faster.
2.109 For example, if the total supply and/or capacities are rather small,
3.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
3.2 +++ b/doc/min_cost_flow.dox Tue May 12 12:06:40 2009 +0200
3.3 @@ -0,0 +1,153 @@
3.4 +/* -*- mode: C++; indent-tabs-mode: nil; -*-
3.5 + *
3.6 + * This file is a part of LEMON, a generic C++ optimization library.
3.7 + *
3.8 + * Copyright (C) 2003-2009
3.9 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
3.10 + * (Egervary Research Group on Combinatorial Optimization, EGRES).
3.11 + *
3.12 + * Permission to use, modify and distribute this software is granted
3.13 + * provided that this copyright notice appears in all copies. For
3.14 + * precise terms see the accompanying LICENSE file.
3.15 + *
3.16 + * This software is provided "AS IS" with no warranty of any kind,
3.17 + * express or implied, and with no claim as to its suitability for any
3.18 + * purpose.
3.19 + *
3.20 + */
3.21 +
3.22 +namespace lemon {
3.23 +
3.24 +/**
3.25 +\page min_cost_flow Minimum Cost Flow Problem
3.26 +
3.27 +\section mcf_def Definition (GEQ form)
3.28 +
3.29 +The \e minimum \e cost \e flow \e problem is to find a feasible flow of
3.30 +minimum total cost from a set of supply nodes to a set of demand nodes
3.31 +in a network with capacity constraints (lower and upper bounds)
3.32 +and arc costs.
3.33 +
3.34 +Formally, let \f$G=(V,A)\f$ be a digraph, \f$lower: A\rightarrow\mathbf{R}\f$,
3.35 +\f$upper: A\rightarrow\mathbf{R}\cup\{+\infty\}\f$ denote the lower and
3.36 +upper bounds for the flow values on the arcs, for which
3.37 +\f$lower(uv) \leq upper(uv)\f$ must hold for all \f$uv\in A\f$,
3.38 +\f$cost: A\rightarrow\mathbf{R}\f$ denotes the cost per unit flow
3.39 +on the arcs and \f$sup: V\rightarrow\mathbf{R}\f$ denotes the
3.40 +signed supply values of the nodes.
3.41 +If \f$sup(u)>0\f$, then \f$u\f$ is a supply node with \f$sup(u)\f$
3.42 +supply, if \f$sup(u)<0\f$, then \f$u\f$ is a demand node with
3.43 +\f$-sup(u)\f$ demand.
3.44 +A minimum cost flow is an \f$f: A\rightarrow\mathbf{R}\f$ solution
3.45 +of the following optimization problem.
3.46 +
3.47 +\f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
3.48 +\f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \geq
3.49 + sup(u) \quad \forall u\in V \f]
3.50 +\f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
3.51 +
3.52 +The sum of the supply values, i.e. \f$\sum_{u\in V} sup(u)\f$ must be
3.53 +zero or negative in order to have a feasible solution (since the sum
3.54 +of the expressions on the left-hand side of the inequalities is zero).
3.55 +It means that the total demand must be greater or equal to the total
3.56 +supply and all the supplies have to be carried out from the supply nodes,
3.57 +but there could be demands that are not satisfied.
3.58 +If \f$\sum_{u\in V} sup(u)\f$ is zero, then all the supply/demand
3.59 +constraints have to be satisfied with equality, i.e. all demands
3.60 +have to be satisfied and all supplies have to be used.
3.61 +
3.62 +
3.63 +\section mcf_algs Algorithms
3.64 +
3.65 +LEMON contains several algorithms for solving this problem, for more
3.66 +information see \ref min_cost_flow_algs "Minimum Cost Flow Algorithms".
3.67 +
3.68 +A feasible solution for this problem can be found using \ref Circulation.
3.69 +
3.70 +
3.71 +\section mcf_dual Dual Solution
3.72 +
3.73 +The dual solution of the minimum cost flow problem is represented by
3.74 +node potentials \f$\pi: V\rightarrow\mathbf{R}\f$.
3.75 +An \f$f: A\rightarrow\mathbf{R}\f$ primal feasible solution is optimal
3.76 +if and only if for some \f$\pi: V\rightarrow\mathbf{R}\f$ node potentials
3.77 +the following \e complementary \e slackness optimality conditions hold.
3.78 +
3.79 + - For all \f$uv\in A\f$ arcs:
3.80 + - if \f$cost^\pi(uv)>0\f$, then \f$f(uv)=lower(uv)\f$;
3.81 + - if \f$lower(uv)<f(uv)<upper(uv)\f$, then \f$cost^\pi(uv)=0\f$;
3.82 + - if \f$cost^\pi(uv)<0\f$, then \f$f(uv)=upper(uv)\f$.
3.83 + - For all \f$u\in V\f$ nodes:
3.84 + - \f$\pi(u)<=0\f$;
3.85 + - if \f$\sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \neq sup(u)\f$,
3.86 + then \f$\pi(u)=0\f$.
3.87 +
3.88 +Here \f$cost^\pi(uv)\f$ denotes the \e reduced \e cost of the arc
3.89 +\f$uv\in A\f$ with respect to the potential function \f$\pi\f$, i.e.
3.90 +\f[ cost^\pi(uv) = cost(uv) + \pi(u) - \pi(v).\f]
3.91 +
3.92 +All algorithms provide dual solution (node potentials), as well,
3.93 +if an optimal flow is found.
3.94 +
3.95 +
3.96 +\section mcf_eq Equality Form
3.97 +
3.98 +The above \ref mcf_def "definition" is actually more general than the
3.99 +usual formulation of the minimum cost flow problem, in which strict
3.100 +equalities are required in the supply/demand contraints.
3.101 +
3.102 +\f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
3.103 +\f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) =
3.104 + sup(u) \quad \forall u\in V \f]
3.105 +\f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
3.106 +
3.107 +However if the sum of the supply values is zero, then these two problems
3.108 +are equivalent.
3.109 +The \ref min_cost_flow_algs "algorithms" in LEMON support the general
3.110 +form, so if you need the equality form, you have to ensure this additional
3.111 +contraint manually.
3.112 +
3.113 +
3.114 +\section mcf_leq Opposite Inequalites (LEQ Form)
3.115 +
3.116 +Another possible definition of the minimum cost flow problem is
3.117 +when there are <em>"less or equal"</em> (LEQ) supply/demand constraints,
3.118 +instead of the <em>"greater or equal"</em> (GEQ) constraints.
3.119 +
3.120 +\f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
3.121 +\f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \leq
3.122 + sup(u) \quad \forall u\in V \f]
3.123 +\f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
3.124 +
3.125 +It means that the total demand must be less or equal to the
3.126 +total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
3.127 +positive) and all the demands have to be satisfied, but there
3.128 +could be supplies that are not carried out from the supply
3.129 +nodes.
3.130 +The equality form is also a special case of this form, of course.
3.131 +
3.132 +You could easily transform this case to the \ref mcf_def "GEQ form"
3.133 +of the problem by reversing the direction of the arcs and taking the
3.134 +negative of the supply values (e.g. using \ref ReverseDigraph and
3.135 +\ref NegMap adaptors).
3.136 +However \ref NetworkSimplex algorithm also supports this form directly
3.137 +for the sake of convenience.
3.138 +
3.139 +Note that the optimality conditions for this supply constraint type are
3.140 +slightly differ from the conditions that are discussed for the GEQ form,
3.141 +namely the potentials have to be non-negative instead of non-positive.
3.142 +An \f$f: A\rightarrow\mathbf{R}\f$ feasible solution of this problem
3.143 +is optimal if and only if for some \f$\pi: V\rightarrow\mathbf{R}\f$
3.144 +node potentials the following conditions hold.
3.145 +
3.146 + - For all \f$uv\in A\f$ arcs:
3.147 + - if \f$cost^\pi(uv)>0\f$, then \f$f(uv)=lower(uv)\f$;
3.148 + - if \f$lower(uv)<f(uv)<upper(uv)\f$, then \f$cost^\pi(uv)=0\f$;
3.149 + - if \f$cost^\pi(uv)<0\f$, then \f$f(uv)=upper(uv)\f$.
3.150 + - For all \f$u\in V\f$ nodes:
3.151 + - \f$\pi(u)>=0\f$;
3.152 + - if \f$\sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \neq sup(u)\f$,
3.153 + then \f$\pi(u)=0\f$.
3.154 +
3.155 +*/
3.156 +}
4.1 --- a/lemon/network_simplex.h Mon May 11 16:38:21 2009 +0100
4.2 +++ b/lemon/network_simplex.h Tue May 12 12:06:40 2009 +0200
4.3 @@ -19,7 +19,7 @@
4.4 #ifndef LEMON_NETWORK_SIMPLEX_H
4.5 #define LEMON_NETWORK_SIMPLEX_H
4.6
4.7 -/// \ingroup min_cost_flow
4.8 +/// \ingroup min_cost_flow_algs
4.9 ///
4.10 /// \file
4.11 /// \brief Network Simplex algorithm for finding a minimum cost flow.
4.12 @@ -33,7 +33,7 @@
4.13
4.14 namespace lemon {
4.15
4.16 - /// \addtogroup min_cost_flow
4.17 + /// \addtogroup min_cost_flow_algs
4.18 /// @{
4.19
4.20 /// \brief Implementation of the primal Network Simplex algorithm
4.21 @@ -102,50 +102,16 @@
4.22 /// i.e. the direction of the inequalities in the supply/demand
4.23 /// constraints of the \ref min_cost_flow "minimum cost flow problem".
4.24 ///
4.25 - /// The default supply type is \c GEQ, since this form is supported
4.26 - /// by other minimum cost flow algorithms and the \ref Circulation
4.27 - /// algorithm, as well.
4.28 - /// The \c LEQ problem type can be selected using the \ref supplyType()
4.29 - /// function.
4.30 - ///
4.31 - /// Note that the equality form is a special case of both supply types.
4.32 + /// The default supply type is \c GEQ, the \c LEQ type can be
4.33 + /// selected using \ref supplyType().
4.34 + /// The equality form is a special case of both supply types.
4.35 enum SupplyType {
4.36 -
4.37 /// This option means that there are <em>"greater or equal"</em>
4.38 - /// supply/demand constraints in the definition, i.e. the exact
4.39 - /// formulation of the problem is the following.
4.40 - /**
4.41 - \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
4.42 - \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \geq
4.43 - sup(u) \quad \forall u\in V \f]
4.44 - \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
4.45 - */
4.46 - /// It means that the total demand must be greater or equal to the
4.47 - /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
4.48 - /// negative) and all the supplies have to be carried out from
4.49 - /// the supply nodes, but there could be demands that are not
4.50 - /// satisfied.
4.51 + /// supply/demand constraints in the definition of the problem.
4.52 GEQ,
4.53 - /// It is just an alias for the \c GEQ option.
4.54 - CARRY_SUPPLIES = GEQ,
4.55 -
4.56 /// This option means that there are <em>"less or equal"</em>
4.57 - /// supply/demand constraints in the definition, i.e. the exact
4.58 - /// formulation of the problem is the following.
4.59 - /**
4.60 - \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
4.61 - \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \leq
4.62 - sup(u) \quad \forall u\in V \f]
4.63 - \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
4.64 - */
4.65 - /// It means that the total demand must be less or equal to the
4.66 - /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
4.67 - /// positive) and all the demands have to be satisfied, but there
4.68 - /// could be supplies that are not carried out from the supply
4.69 - /// nodes.
4.70 - LEQ,
4.71 - /// It is just an alias for the \c LEQ option.
4.72 - SATISFY_DEMANDS = LEQ
4.73 + /// supply/demand constraints in the definition of the problem.
4.74 + LEQ
4.75 };
4.76
4.77 /// \brief Constants for selecting the pivot rule.
4.78 @@ -215,6 +181,8 @@
4.79 const GR &_graph;
4.80 int _node_num;
4.81 int _arc_num;
4.82 + int _all_arc_num;
4.83 + int _search_arc_num;
4.84
4.85 // Parameters of the problem
4.86 bool _have_lower;
4.87 @@ -277,7 +245,7 @@
4.88 const IntVector &_state;
4.89 const CostVector &_pi;
4.90 int &_in_arc;
4.91 - int _arc_num;
4.92 + int _search_arc_num;
4.93
4.94 // Pivot rule data
4.95 int _next_arc;
4.96 @@ -288,13 +256,14 @@
4.97 FirstEligiblePivotRule(NetworkSimplex &ns) :
4.98 _source(ns._source), _target(ns._target),
4.99 _cost(ns._cost), _state(ns._state), _pi(ns._pi),
4.100 - _in_arc(ns.in_arc), _arc_num(ns._arc_num), _next_arc(0)
4.101 + _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num),
4.102 + _next_arc(0)
4.103 {}
4.104
4.105 // Find next entering arc
4.106 bool findEnteringArc() {
4.107 Cost c;
4.108 - for (int e = _next_arc; e < _arc_num; ++e) {
4.109 + for (int e = _next_arc; e < _search_arc_num; ++e) {
4.110 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
4.111 if (c < 0) {
4.112 _in_arc = e;
4.113 @@ -328,7 +297,7 @@
4.114 const IntVector &_state;
4.115 const CostVector &_pi;
4.116 int &_in_arc;
4.117 - int _arc_num;
4.118 + int _search_arc_num;
4.119
4.120 public:
4.121
4.122 @@ -336,13 +305,13 @@
4.123 BestEligiblePivotRule(NetworkSimplex &ns) :
4.124 _source(ns._source), _target(ns._target),
4.125 _cost(ns._cost), _state(ns._state), _pi(ns._pi),
4.126 - _in_arc(ns.in_arc), _arc_num(ns._arc_num)
4.127 + _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num)
4.128 {}
4.129
4.130 // Find next entering arc
4.131 bool findEnteringArc() {
4.132 Cost c, min = 0;
4.133 - for (int e = 0; e < _arc_num; ++e) {
4.134 + for (int e = 0; e < _search_arc_num; ++e) {
4.135 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
4.136 if (c < min) {
4.137 min = c;
4.138 @@ -367,7 +336,7 @@
4.139 const IntVector &_state;
4.140 const CostVector &_pi;
4.141 int &_in_arc;
4.142 - int _arc_num;
4.143 + int _search_arc_num;
4.144
4.145 // Pivot rule data
4.146 int _block_size;
4.147 @@ -379,14 +348,15 @@
4.148 BlockSearchPivotRule(NetworkSimplex &ns) :
4.149 _source(ns._source), _target(ns._target),
4.150 _cost(ns._cost), _state(ns._state), _pi(ns._pi),
4.151 - _in_arc(ns.in_arc), _arc_num(ns._arc_num), _next_arc(0)
4.152 + _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num),
4.153 + _next_arc(0)
4.154 {
4.155 // The main parameters of the pivot rule
4.156 - const double BLOCK_SIZE_FACTOR = 2.0;
4.157 + const double BLOCK_SIZE_FACTOR = 0.5;
4.158 const int MIN_BLOCK_SIZE = 10;
4.159
4.160 _block_size = std::max( int(BLOCK_SIZE_FACTOR *
4.161 - std::sqrt(double(_arc_num))),
4.162 + std::sqrt(double(_search_arc_num))),
4.163 MIN_BLOCK_SIZE );
4.164 }
4.165
4.166 @@ -395,7 +365,7 @@
4.167 Cost c, min = 0;
4.168 int cnt = _block_size;
4.169 int e, min_arc = _next_arc;
4.170 - for (e = _next_arc; e < _arc_num; ++e) {
4.171 + for (e = _next_arc; e < _search_arc_num; ++e) {
4.172 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
4.173 if (c < min) {
4.174 min = c;
4.175 @@ -440,7 +410,7 @@
4.176 const IntVector &_state;
4.177 const CostVector &_pi;
4.178 int &_in_arc;
4.179 - int _arc_num;
4.180 + int _search_arc_num;
4.181
4.182 // Pivot rule data
4.183 IntVector _candidates;
4.184 @@ -454,7 +424,8 @@
4.185 CandidateListPivotRule(NetworkSimplex &ns) :
4.186 _source(ns._source), _target(ns._target),
4.187 _cost(ns._cost), _state(ns._state), _pi(ns._pi),
4.188 - _in_arc(ns.in_arc), _arc_num(ns._arc_num), _next_arc(0)
4.189 + _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num),
4.190 + _next_arc(0)
4.191 {
4.192 // The main parameters of the pivot rule
4.193 const double LIST_LENGTH_FACTOR = 1.0;
4.194 @@ -463,7 +434,7 @@
4.195 const int MIN_MINOR_LIMIT = 3;
4.196
4.197 _list_length = std::max( int(LIST_LENGTH_FACTOR *
4.198 - std::sqrt(double(_arc_num))),
4.199 + std::sqrt(double(_search_arc_num))),
4.200 MIN_LIST_LENGTH );
4.201 _minor_limit = std::max( int(MINOR_LIMIT_FACTOR * _list_length),
4.202 MIN_MINOR_LIMIT );
4.203 @@ -500,7 +471,7 @@
4.204 // Major iteration: build a new candidate list
4.205 min = 0;
4.206 _curr_length = 0;
4.207 - for (e = _next_arc; e < _arc_num; ++e) {
4.208 + for (e = _next_arc; e < _search_arc_num; ++e) {
4.209 c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
4.210 if (c < 0) {
4.211 _candidates[_curr_length++] = e;
4.212 @@ -546,7 +517,7 @@
4.213 const IntVector &_state;
4.214 const CostVector &_pi;
4.215 int &_in_arc;
4.216 - int _arc_num;
4.217 + int _search_arc_num;
4.218
4.219 // Pivot rule data
4.220 int _block_size, _head_length, _curr_length;
4.221 @@ -574,8 +545,8 @@
4.222 AlteringListPivotRule(NetworkSimplex &ns) :
4.223 _source(ns._source), _target(ns._target),
4.224 _cost(ns._cost), _state(ns._state), _pi(ns._pi),
4.225 - _in_arc(ns.in_arc), _arc_num(ns._arc_num),
4.226 - _next_arc(0), _cand_cost(ns._arc_num), _sort_func(_cand_cost)
4.227 + _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num),
4.228 + _next_arc(0), _cand_cost(ns._search_arc_num), _sort_func(_cand_cost)
4.229 {
4.230 // The main parameters of the pivot rule
4.231 const double BLOCK_SIZE_FACTOR = 1.5;
4.232 @@ -584,7 +555,7 @@
4.233 const int MIN_HEAD_LENGTH = 3;
4.234
4.235 _block_size = std::max( int(BLOCK_SIZE_FACTOR *
4.236 - std::sqrt(double(_arc_num))),
4.237 + std::sqrt(double(_search_arc_num))),
4.238 MIN_BLOCK_SIZE );
4.239 _head_length = std::max( int(HEAD_LENGTH_FACTOR * _block_size),
4.240 MIN_HEAD_LENGTH );
4.241 @@ -610,7 +581,7 @@
4.242 int last_arc = 0;
4.243 int limit = _head_length;
4.244
4.245 - for (int e = _next_arc; e < _arc_num; ++e) {
4.246 + for (int e = _next_arc; e < _search_arc_num; ++e) {
4.247 _cand_cost[e] = _state[e] *
4.248 (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
4.249 if (_cand_cost[e] < 0) {
4.250 @@ -678,17 +649,17 @@
4.251 _node_num = countNodes(_graph);
4.252 _arc_num = countArcs(_graph);
4.253 int all_node_num = _node_num + 1;
4.254 - int all_arc_num = _arc_num + _node_num;
4.255 + int max_arc_num = _arc_num + 2 * _node_num;
4.256
4.257 - _source.resize(all_arc_num);
4.258 - _target.resize(all_arc_num);
4.259 + _source.resize(max_arc_num);
4.260 + _target.resize(max_arc_num);
4.261
4.262 - _lower.resize(all_arc_num);
4.263 - _upper.resize(all_arc_num);
4.264 - _cap.resize(all_arc_num);
4.265 - _cost.resize(all_arc_num);
4.266 + _lower.resize(_arc_num);
4.267 + _upper.resize(_arc_num);
4.268 + _cap.resize(max_arc_num);
4.269 + _cost.resize(max_arc_num);
4.270 _supply.resize(all_node_num);
4.271 - _flow.resize(all_arc_num);
4.272 + _flow.resize(max_arc_num);
4.273 _pi.resize(all_node_num);
4.274
4.275 _parent.resize(all_node_num);
4.276 @@ -698,7 +669,7 @@
4.277 _rev_thread.resize(all_node_num);
4.278 _succ_num.resize(all_node_num);
4.279 _last_succ.resize(all_node_num);
4.280 - _state.resize(all_arc_num);
4.281 + _state.resize(max_arc_num);
4.282
4.283 // Copy the graph (store the arcs in a mixed order)
4.284 int i = 0;
4.285 @@ -1069,7 +1040,7 @@
4.286 // Initialize artifical cost
4.287 Cost ART_COST;
4.288 if (std::numeric_limits<Cost>::is_exact) {
4.289 - ART_COST = std::numeric_limits<Cost>::max() / 4 + 1;
4.290 + ART_COST = std::numeric_limits<Cost>::max() / 2 + 1;
4.291 } else {
4.292 ART_COST = std::numeric_limits<Cost>::min();
4.293 for (int i = 0; i != _arc_num; ++i) {
4.294 @@ -1093,29 +1064,121 @@
4.295 _succ_num[_root] = _node_num + 1;
4.296 _last_succ[_root] = _root - 1;
4.297 _supply[_root] = -_sum_supply;
4.298 - _pi[_root] = _sum_supply < 0 ? -ART_COST : ART_COST;
4.299 + _pi[_root] = 0;
4.300
4.301 // Add artificial arcs and initialize the spanning tree data structure
4.302 - for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
4.303 - _parent[u] = _root;
4.304 - _pred[u] = e;
4.305 - _thread[u] = u + 1;
4.306 - _rev_thread[u + 1] = u;
4.307 - _succ_num[u] = 1;
4.308 - _last_succ[u] = u;
4.309 - _cost[e] = ART_COST;
4.310 - _cap[e] = INF;
4.311 - _state[e] = STATE_TREE;
4.312 - if (_supply[u] > 0 || (_supply[u] == 0 && _sum_supply <= 0)) {
4.313 - _flow[e] = _supply[u];
4.314 - _forward[u] = true;
4.315 - _pi[u] = -ART_COST + _pi[_root];
4.316 - } else {
4.317 - _flow[e] = -_supply[u];
4.318 - _forward[u] = false;
4.319 - _pi[u] = ART_COST + _pi[_root];
4.320 + if (_sum_supply == 0) {
4.321 + // EQ supply constraints
4.322 + _search_arc_num = _arc_num;
4.323 + _all_arc_num = _arc_num + _node_num;
4.324 + for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
4.325 + _parent[u] = _root;
4.326 + _pred[u] = e;
4.327 + _thread[u] = u + 1;
4.328 + _rev_thread[u + 1] = u;
4.329 + _succ_num[u] = 1;
4.330 + _last_succ[u] = u;
4.331 + _cap[e] = INF;
4.332 + _state[e] = STATE_TREE;
4.333 + if (_supply[u] >= 0) {
4.334 + _forward[u] = true;
4.335 + _pi[u] = 0;
4.336 + _source[e] = u;
4.337 + _target[e] = _root;
4.338 + _flow[e] = _supply[u];
4.339 + _cost[e] = 0;
4.340 + } else {
4.341 + _forward[u] = false;
4.342 + _pi[u] = ART_COST;
4.343 + _source[e] = _root;
4.344 + _target[e] = u;
4.345 + _flow[e] = -_supply[u];
4.346 + _cost[e] = ART_COST;
4.347 + }
4.348 }
4.349 }
4.350 + else if (_sum_supply > 0) {
4.351 + // LEQ supply constraints
4.352 + _search_arc_num = _arc_num + _node_num;
4.353 + int f = _arc_num + _node_num;
4.354 + for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
4.355 + _parent[u] = _root;
4.356 + _thread[u] = u + 1;
4.357 + _rev_thread[u + 1] = u;
4.358 + _succ_num[u] = 1;
4.359 + _last_succ[u] = u;
4.360 + if (_supply[u] >= 0) {
4.361 + _forward[u] = true;
4.362 + _pi[u] = 0;
4.363 + _pred[u] = e;
4.364 + _source[e] = u;
4.365 + _target[e] = _root;
4.366 + _cap[e] = INF;
4.367 + _flow[e] = _supply[u];
4.368 + _cost[e] = 0;
4.369 + _state[e] = STATE_TREE;
4.370 + } else {
4.371 + _forward[u] = false;
4.372 + _pi[u] = ART_COST;
4.373 + _pred[u] = f;
4.374 + _source[f] = _root;
4.375 + _target[f] = u;
4.376 + _cap[f] = INF;
4.377 + _flow[f] = -_supply[u];
4.378 + _cost[f] = ART_COST;
4.379 + _state[f] = STATE_TREE;
4.380 + _source[e] = u;
4.381 + _target[e] = _root;
4.382 + _cap[e] = INF;
4.383 + _flow[e] = 0;
4.384 + _cost[e] = 0;
4.385 + _state[e] = STATE_LOWER;
4.386 + ++f;
4.387 + }
4.388 + }
4.389 + _all_arc_num = f;
4.390 + }
4.391 + else {
4.392 + // GEQ supply constraints
4.393 + _search_arc_num = _arc_num + _node_num;
4.394 + int f = _arc_num + _node_num;
4.395 + for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
4.396 + _parent[u] = _root;
4.397 + _thread[u] = u + 1;
4.398 + _rev_thread[u + 1] = u;
4.399 + _succ_num[u] = 1;
4.400 + _last_succ[u] = u;
4.401 + if (_supply[u] <= 0) {
4.402 + _forward[u] = false;
4.403 + _pi[u] = 0;
4.404 + _pred[u] = e;
4.405 + _source[e] = _root;
4.406 + _target[e] = u;
4.407 + _cap[e] = INF;
4.408 + _flow[e] = -_supply[u];
4.409 + _cost[e] = 0;
4.410 + _state[e] = STATE_TREE;
4.411 + } else {
4.412 + _forward[u] = true;
4.413 + _pi[u] = -ART_COST;
4.414 + _pred[u] = f;
4.415 + _source[f] = u;
4.416 + _target[f] = _root;
4.417 + _cap[f] = INF;
4.418 + _flow[f] = _supply[u];
4.419 + _state[f] = STATE_TREE;
4.420 + _cost[f] = ART_COST;
4.421 + _source[e] = _root;
4.422 + _target[e] = u;
4.423 + _cap[e] = INF;
4.424 + _flow[e] = 0;
4.425 + _cost[e] = 0;
4.426 + _state[e] = STATE_LOWER;
4.427 + ++f;
4.428 + }
4.429 + }
4.430 + _all_arc_num = f;
4.431 + }
4.432
4.433 return true;
4.434 }
4.435 @@ -1374,20 +1437,8 @@
4.436 }
4.437
4.438 // Check feasibility
4.439 - if (_sum_supply < 0) {
4.440 - for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
4.441 - if (_supply[u] >= 0 && _flow[e] != 0) return INFEASIBLE;
4.442 - }
4.443 - }
4.444 - else if (_sum_supply > 0) {
4.445 - for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
4.446 - if (_supply[u] <= 0 && _flow[e] != 0) return INFEASIBLE;
4.447 - }
4.448 - }
4.449 - else {
4.450 - for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
4.451 - if (_flow[e] != 0) return INFEASIBLE;
4.452 - }
4.453 + for (int e = _search_arc_num; e != _all_arc_num; ++e) {
4.454 + if (_flow[e] != 0) return INFEASIBLE;
4.455 }
4.456
4.457 // Transform the solution and the supply map to the original form
4.458 @@ -1401,6 +1452,30 @@
4.459 }
4.460 }
4.461 }
4.462 +
4.463 + // Shift potentials to meet the requirements of the GEQ/LEQ type
4.464 + // optimality conditions
4.465 + if (_sum_supply == 0) {
4.466 + if (_stype == GEQ) {
4.467 + Cost max_pot = std::numeric_limits<Cost>::min();
4.468 + for (int i = 0; i != _node_num; ++i) {
4.469 + if (_pi[i] > max_pot) max_pot = _pi[i];
4.470 + }
4.471 + if (max_pot > 0) {
4.472 + for (int i = 0; i != _node_num; ++i)
4.473 + _pi[i] -= max_pot;
4.474 + }
4.475 + } else {
4.476 + Cost min_pot = std::numeric_limits<Cost>::max();
4.477 + for (int i = 0; i != _node_num; ++i) {
4.478 + if (_pi[i] < min_pot) min_pot = _pi[i];
4.479 + }
4.480 + if (min_pot < 0) {
4.481 + for (int i = 0; i != _node_num; ++i)
4.482 + _pi[i] -= min_pot;
4.483 + }
4.484 + }
4.485 + }
4.486
4.487 return OPTIMAL;
4.488 }