COIN-OR::LEMON - Graph Library

Ticket #291: 291-ns-fixes-e239cde7ae57.patch

File 291-ns-fixes-e239cde7ae57.patch, 29.0 KB (added by Peter Kovacs, 12 years ago)
  • doc/Makefile.am

    # HG changeset patch
    # User Peter Kovacs <kpeter@inf.elte.hu>
    # Date 1242122800 -7200
    # Node ID e239cde7ae572839447f31364ff9e67ae86a059a
    # Parent  ca92c2f936b03d0966865203a090abbd2e703cbc
    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.
    
    diff --git a/doc/Makefile.am b/doc/Makefile.am
    a b  
    88        doc/license.dox \
    99        doc/mainpage.dox \
    1010        doc/migration.dox \
     11        doc/min_cost_flow.dox \
    1112        doc/named-param.dox \
    1213        doc/namespaces.dox \
    1314        doc/html \
  • doc/groups.dox

    diff --git a/doc/groups.dox b/doc/groups.dox
    a b  
    345345*/
    346346
    347347/**
    348 @defgroup min_cost_flow Minimum Cost Flow Algorithms
     348@defgroup min_cost_flow_algs Minimum Cost Flow Algorithms
    349349@ingroup algs
    350350
    351351\brief Algorithms for finding minimum cost flows and circulations.
    352352
    353353This group contains the algorithms for finding minimum cost flows and
    354 circulations.
     354circulations. For more information about this problem and its dual
     355solution see \ref min_cost_flow "Minimum Cost Flow Problem".
    355356
    356 The \e minimum \e cost \e flow \e problem is to find a feasible flow of
    357 minimum total cost from a set of supply nodes to a set of demand nodes
    358 in a network with capacity constraints (lower and upper bounds)
    359 and arc costs.
    360 Formally, let \f$G=(V,A)\f$ be a digraph, \f$lower: A\rightarrow\mathbf{Z}\f$,
    361 \f$upper: A\rightarrow\mathbf{Z}\cup\{+\infty\}\f$ denote the lower and
    362 upper bounds for the flow values on the arcs, for which
    363 \f$lower(uv) \leq upper(uv)\f$ must hold for all \f$uv\in A\f$,
    364 \f$cost: A\rightarrow\mathbf{Z}\f$ denotes the cost per unit flow
    365 on the arcs and \f$sup: V\rightarrow\mathbf{Z}\f$ denotes the
    366 signed supply values of the nodes.
    367 If \f$sup(u)>0\f$, then \f$u\f$ is a supply node with \f$sup(u)\f$
    368 supply, if \f$sup(u)<0\f$, then \f$u\f$ is a demand node with
    369 \f$-sup(u)\f$ demand.
    370 A minimum cost flow is an \f$f: A\rightarrow\mathbf{Z}\f$ solution
    371 of the following optimization problem.
    372 
    373 \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
    374 \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \geq
    375     sup(u) \quad \forall u\in V \f]
    376 \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
    377 
    378 The sum of the supply values, i.e. \f$\sum_{u\in V} sup(u)\f$ must be
    379 zero or negative in order to have a feasible solution (since the sum
    380 of the expressions on the left-hand side of the inequalities is zero).
    381 It means that the total demand must be greater or equal to the total
    382 supply and all the supplies have to be carried out from the supply nodes,
    383 but there could be demands that are not satisfied.
    384 If \f$\sum_{u\in V} sup(u)\f$ is zero, then all the supply/demand
    385 constraints have to be satisfied with equality, i.e. all demands
    386 have to be satisfied and all supplies have to be used.
    387 
    388 If you need the opposite inequalities in the supply/demand constraints
    389 (i.e. the total demand is less than the total supply and all the demands
    390 have to be satisfied while there could be supplies that are not used),
    391 then you could easily transform the problem to the above form by reversing
    392 the direction of the arcs and taking the negative of the supply values
    393 (e.g. using \ref ReverseDigraph and \ref NegMap adaptors).
    394 However \ref NetworkSimplex algorithm also supports this form directly
    395 for the sake of convenience.
    396 
    397 A feasible solution for this problem can be found using \ref Circulation.
    398 
    399 Note that the above formulation is actually more general than the usual
    400 definition of the minimum cost flow problem, in which strict equalities
    401 are required in the supply/demand contraints, i.e.
    402 
    403 \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) =
    404     sup(u) \quad \forall u\in V. \f]
    405 
    406 However if the sum of the supply values is zero, then these two problems
    407 are equivalent. So if you need the equality form, you have to ensure this
    408 additional contraint for the algorithms.
    409 
    410 The dual solution of the minimum cost flow problem is represented by node
    411 potentials \f$\pi: V\rightarrow\mathbf{Z}\f$.
    412 An \f$f: A\rightarrow\mathbf{Z}\f$ feasible solution of the problem
    413 is optimal if and only if for some \f$\pi: V\rightarrow\mathbf{Z}\f$
    414 node potentials the following \e complementary \e slackness optimality
    415 conditions hold.
    416 
    417  - For all \f$uv\in A\f$ arcs:
    418    - if \f$cost^\pi(uv)>0\f$, then \f$f(uv)=lower(uv)\f$;
    419    - if \f$lower(uv)<f(uv)<upper(uv)\f$, then \f$cost^\pi(uv)=0\f$;
    420    - if \f$cost^\pi(uv)<0\f$, then \f$f(uv)=upper(uv)\f$.
    421  - For all \f$u\in V\f$ nodes:
    422    - if \f$\sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \neq sup(u)\f$,
    423      then \f$\pi(u)=0\f$.
    424  
    425 Here \f$cost^\pi(uv)\f$ denotes the \e reduced \e cost of the arc
    426 \f$uv\in A\f$ with respect to the potential function \f$\pi\f$, i.e.
    427 \f[ cost^\pi(uv) = cost(uv) + \pi(u) - \pi(v).\f]
    428 
    429 All algorithms provide dual solution (node potentials) as well,
    430 if an optimal flow is found.
    431 
    432 LEMON contains several algorithms for solving minimum cost flow problems.
     357LEMON contains several algorithms for this problem.
    433358 - \ref NetworkSimplex Primal Network Simplex algorithm with various
    434359   pivot strategies.
    435360 - \ref CostScaling Push-Relabel and Augment-Relabel algorithms based on
     
    439364 - \ref CancelAndTighten The Cancel and Tighten algorithm.
    440365 - \ref CycleCanceling Cycle-Canceling algorithms.
    441366
    442 Most of these implementations support the general inequality form of the
    443 minimum cost flow problem, but CancelAndTighten and CycleCanceling
    444 only support the equality form due to the primal method they use.
    445 
    446367In general NetworkSimplex is the most efficient implementation,
    447368but in special cases other algorithms could be faster.
    448369For example, if the total supply and/or capacities are rather small,
  • new file doc/min_cost_flow.dox

    diff --git a/doc/min_cost_flow.dox b/doc/min_cost_flow.dox
    new file mode 100644
    - +  
     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
     19namespace lemon {
     20
     21/**
     22\page min_cost_flow Minimum Cost Flow Problem
     23
     24\section mcf_def Definition
     25
     26The \e minimum \e cost \e flow \e problem is to find a feasible flow of
     27minimum total cost from a set of supply nodes to a set of demand nodes
     28in a network with capacity constraints (lower and upper bounds)
     29and arc costs.
     30
     31Formally, 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
     33upper 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
     36on the arcs and \f$sup: V\rightarrow\mathbf{R}\f$ denotes the
     37signed supply values of the nodes.
     38If \f$sup(u)>0\f$, then \f$u\f$ is a supply node with \f$sup(u)\f$
     39supply, if \f$sup(u)<0\f$, then \f$u\f$ is a demand node with
     40\f$-sup(u)\f$ demand.
     41A minimum cost flow is an \f$f: A\rightarrow\mathbf{R}\f$ solution
     42of 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
     49The sum of the supply values, i.e. \f$\sum_{u\in V} sup(u)\f$ must be
     50zero or negative in order to have a feasible solution (since the sum
     51of the expressions on the left-hand side of the inequalities is zero).
     52It means that the total demand must be greater or equal to the total
     53supply and all the supplies have to be carried out from the supply nodes,
     54but there could be demands that are not satisfied.
     55If \f$\sum_{u\in V} sup(u)\f$ is zero, then all the supply/demand
     56constraints have to be satisfied with equality, i.e. all demands
     57have to be satisfied and all supplies have to be used.
     58
     59
     60\section mcf_algs Algorithms
     61
     62LEMON contains several algorithms for solving this problem, for more
     63information see \ref min_cost_flow_algs "Minimum Cost Flow Algorithms".
     64
     65A feasible solution for this problem can be found using \ref Circulation.
     66
     67
     68\section mcf_dual Dual Solution
     69
     70The dual solution of the minimum cost flow problem is represented by
     71node potentials \f$\pi: V\rightarrow\mathbf{R}\f$.
     72An \f$f: A\rightarrow\mathbf{R}\f$ primal feasible solution is optimal
     73if and only if for some \f$\pi: V\rightarrow\mathbf{R}\f$ node potentials
     74the 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 
     85Here \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
     89All algorithms provide dual solution (node potentials), as well,
     90if an optimal flow is found.
     91
     92
     93\section mcf_eq Equality Form
     94
     95The above \ref mcf_def "definition" is actually more general than the
     96usual formulation of the minimum cost flow problem, in which strict
     97equalities 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
     104However if the sum of the supply values is zero, then these two problems
     105are equivalent.
     106The \ref min_cost_flow_algs "algorithms" in LEMON support the general
     107form, so if you need the equality form, you have to ensure this additional
     108contraint manually.
     109
     110
     111\section mcf_leq Opposite Inequality Form
     112
     113Another possible definition of the minimum cost flow problem is
     114when there are <em>"less or equal"</em> (LEQ) supply/demand constraints,
     115instead 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
     122It means that the total demand must be less or equal to the
     123total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
     124positive) and all the demands have to be satisfied, but there
     125could be supplies that are not carried out from the supply
     126nodes.
     127The equality form is also a special case of this form, of course.
     128
     129You could easily transform this case to the \ref mcf_def "GEQ form"
     130of the problem by reversing the direction of the arcs and taking the
     131negative of the supply values (e.g. using \ref ReverseDigraph and
     132\ref NegMap adaptors).
     133However \ref NetworkSimplex algorithm also supports this form directly
     134for the sake of convenience.
     135
     136Note that the optimality conditions for this supply constraint type are
     137slightly differ from the conditions that are discussed for the GEQ form,
     138namely the potentials have to be non-negative instead of non-positive.
     139An \f$f: A\rightarrow\mathbf{R}\f$ feasible solution of this problem
     140is optimal if and only if for some \f$\pi: V\rightarrow\mathbf{R}\f$
     141node 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}
  • lemon/network_simplex.h

    diff --git a/lemon/network_simplex.h b/lemon/network_simplex.h
    a b  
    1919#ifndef LEMON_NETWORK_SIMPLEX_H
    2020#define LEMON_NETWORK_SIMPLEX_H
    2121
    22 /// \ingroup min_cost_flow
     22/// \ingroup min_cost_flow_algs
    2323///
    2424/// \file
    2525/// \brief Network Simplex algorithm for finding a minimum cost flow.
     
    3333
    3434namespace lemon {
    3535
    36   /// \addtogroup min_cost_flow
     36  /// \addtogroup min_cost_flow_algs
    3737  /// @{
    3838
    3939  /// \brief Implementation of the primal Network Simplex algorithm
     
    102102    /// i.e. the direction of the inequalities in the supply/demand
    103103    /// constraints of the \ref min_cost_flow "minimum cost flow problem".
    104104    ///
    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.
    112108    enum SupplyType {
    113 
    114109      /// 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.
    128111      GEQ,
    129       /// It is just an alias for the \c GEQ option.
    130       CARRY_SUPPLIES = GEQ,
    131 
    132112      /// 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
    149115    };
    150116   
    151117    /// \brief Constants for selecting the pivot rule.
     
    215181    const GR &_graph;
    216182    int _node_num;
    217183    int _arc_num;
     184    int _all_arc_num;
     185    int _search_arc_num;
    218186
    219187    // Parameters of the problem
    220188    bool _have_lower;
     
    277245      const IntVector  &_state;
    278246      const CostVector &_pi;
    279247      int &_in_arc;
    280       int _arc_num;
     248      int _search_arc_num;
    281249
    282250      // Pivot rule data
    283251      int _next_arc;
     
    288256      FirstEligiblePivotRule(NetworkSimplex &ns) :
    289257        _source(ns._source), _target(ns._target),
    290258        _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)
    292261      {}
    293262
    294263      // Find next entering arc
    295264      bool findEnteringArc() {
    296265        Cost c;
    297         for (int e = _next_arc; e < _arc_num; ++e) {
     266        for (int e = _next_arc; e < _search_arc_num; ++e) {
    298267          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
    299268          if (c < 0) {
    300269            _in_arc = e;
     
    328297      const IntVector  &_state;
    329298      const CostVector &_pi;
    330299      int &_in_arc;
    331       int _arc_num;
     300      int _search_arc_num;
    332301
    333302    public:
    334303
     
    336305      BestEligiblePivotRule(NetworkSimplex &ns) :
    337306        _source(ns._source), _target(ns._target),
    338307        _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)
    340309      {}
    341310
    342311      // Find next entering arc
    343312      bool findEnteringArc() {
    344313        Cost c, min = 0;
    345         for (int e = 0; e < _arc_num; ++e) {
     314        for (int e = 0; e < _search_arc_num; ++e) {
    346315          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
    347316          if (c < min) {
    348317            min = c;
     
    367336      const IntVector  &_state;
    368337      const CostVector &_pi;
    369338      int &_in_arc;
    370       int _arc_num;
     339      int _search_arc_num;
    371340
    372341      // Pivot rule data
    373342      int _block_size;
     
    379348      BlockSearchPivotRule(NetworkSimplex &ns) :
    380349        _source(ns._source), _target(ns._target),
    381350        _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)
    383353      {
    384354        // The main parameters of the pivot rule
    385         const double BLOCK_SIZE_FACTOR = 2.0;
     355        const double BLOCK_SIZE_FACTOR = 0.5;
    386356        const int MIN_BLOCK_SIZE = 10;
    387357
    388358        _block_size = std::max( int(BLOCK_SIZE_FACTOR *
    389                                     std::sqrt(double(_arc_num))),
     359                                    std::sqrt(double(_search_arc_num))),
    390360                                MIN_BLOCK_SIZE );
    391361      }
    392362
     
    395365        Cost c, min = 0;
    396366        int cnt = _block_size;
    397367        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) {
    399369          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
    400370          if (c < min) {
    401371            min = c;
     
    440410      const IntVector  &_state;
    441411      const CostVector &_pi;
    442412      int &_in_arc;
    443       int _arc_num;
     413      int _search_arc_num;
    444414
    445415      // Pivot rule data
    446416      IntVector _candidates;
     
    454424      CandidateListPivotRule(NetworkSimplex &ns) :
    455425        _source(ns._source), _target(ns._target),
    456426        _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)
    458429      {
    459430        // The main parameters of the pivot rule
    460431        const double LIST_LENGTH_FACTOR = 1.0;
     
    463434        const int MIN_MINOR_LIMIT = 3;
    464435
    465436        _list_length = std::max( int(LIST_LENGTH_FACTOR *
    466                                      std::sqrt(double(_arc_num))),
     437                                     std::sqrt(double(_search_arc_num))),
    467438                                 MIN_LIST_LENGTH );
    468439        _minor_limit = std::max( int(MINOR_LIMIT_FACTOR * _list_length),
    469440                                 MIN_MINOR_LIMIT );
     
    500471        // Major iteration: build a new candidate list
    501472        min = 0;
    502473        _curr_length = 0;
    503         for (e = _next_arc; e < _arc_num; ++e) {
     474        for (e = _next_arc; e < _search_arc_num; ++e) {
    504475          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
    505476          if (c < 0) {
    506477            _candidates[_curr_length++] = e;
     
    546517      const IntVector  &_state;
    547518      const CostVector &_pi;
    548519      int &_in_arc;
    549       int _arc_num;
     520      int _search_arc_num;
    550521
    551522      // Pivot rule data
    552523      int _block_size, _head_length, _curr_length;
     
    574545      AlteringListPivotRule(NetworkSimplex &ns) :
    575546        _source(ns._source), _target(ns._target),
    576547        _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)
    579550      {
    580551        // The main parameters of the pivot rule
    581552        const double BLOCK_SIZE_FACTOR = 1.5;
     
    584555        const int MIN_HEAD_LENGTH = 3;
    585556
    586557        _block_size = std::max( int(BLOCK_SIZE_FACTOR *
    587                                     std::sqrt(double(_arc_num))),
     558                                    std::sqrt(double(_search_arc_num))),
    588559                                MIN_BLOCK_SIZE );
    589560        _head_length = std::max( int(HEAD_LENGTH_FACTOR * _block_size),
    590561                                 MIN_HEAD_LENGTH );
     
    610581        int last_arc = 0;
    611582        int limit = _head_length;
    612583
    613         for (int e = _next_arc; e < _arc_num; ++e) {
     584        for (int e = _next_arc; e < _search_arc_num; ++e) {
    614585          _cand_cost[e] = _state[e] *
    615586            (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
    616587          if (_cand_cost[e] < 0) {
     
    678649      _node_num = countNodes(_graph);
    679650      _arc_num = countArcs(_graph);
    680651      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;
    682653
    683       _source.resize(all_arc_num);
    684       _target.resize(all_arc_num);
     654      _source.resize(max_arc_num);
     655      _target.resize(max_arc_num);
    685656
    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);
    690661      _supply.resize(all_node_num);
    691       _flow.resize(all_arc_num);
     662      _flow.resize(max_arc_num);
    692663      _pi.resize(all_node_num);
    693664
    694665      _parent.resize(all_node_num);
     
    698669      _rev_thread.resize(all_node_num);
    699670      _succ_num.resize(all_node_num);
    700671      _last_succ.resize(all_node_num);
    701       _state.resize(all_arc_num);
     672      _state.resize(max_arc_num);
    702673
    703674      // Copy the graph (store the arcs in a mixed order)
    704675      int i = 0;
     
    10691040      // Initialize artifical cost
    10701041      Cost ART_COST;
    10711042      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;
    10731044      } else {
    10741045        ART_COST = std::numeric_limits<Cost>::min();
    10751046        for (int i = 0; i != _arc_num; ++i) {
     
    10931064      _succ_num[_root] = _node_num + 1;
    10941065      _last_succ[_root] = _root - 1;
    10951066      _supply[_root] = -_sum_supply;
    1096       _pi[_root] = _sum_supply < 0 ? -ART_COST : ART_COST;
     1067      _pi[_root] = 0;
    10971068
    10981069      // 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          }
    11171098        }
    11181099      }
     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      }
    11191182
    11201183      return true;
    11211184    }
     
    13741437      }
    13751438     
    13761439      // 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;
    13911442      }
    13921443
    13931444      // Transform the solution and the supply map to the original form
     
    14011452          }
    14021453        }
    14031454      }
     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      }
    14041479
    14051480      return OPTIMAL;
    14061481    }