gravatar
kpeter (Peter Kovacs)
kpeter@inf.elte.hu
Support negative costs and bounds in NetworkSimplex (#270) * The interface is reworked to support negative costs and bounds. - ProblemType and problemType() are renamed to SupplyType and supplyType(), see also #234. - ProblemType type is introduced similarly to the LP interface. - 'bool run()' is replaced by 'ProblemType run()' to handle unbounded problem instances, as well. - Add INF public member constant similarly to the LP interface. * Remove capacityMap() and boundMaps(), see also #266. * Update the problem definition in the MCF module. * Remove the usage of Circulation (and adaptors) for checking feasibility. Check feasibility by examining the artifical arcs instead (after solving the problem). * Additional check for unbounded negative cycles found during the algorithm (it is possible now, since negative costs are allowed). * Fix in the constructor (the value types needn't be integer any more), see also #254. * Improve and extend the doc. * Rework the test file and add test cases for negative costs and bounds.
0 4 0
default
4 files changed with 348 insertions and 330 deletions:
↑ Collapse diff ↑
Show white space 6 line context
... ...
@@ -352,17 +352,17 @@
352 352
minimum total cost from a set of supply nodes to a set of demand nodes
353 353
in a network with capacity constraints (lower and upper bounds)
354 354
and arc costs.
355
Formally, let \f$G=(V,A)\f$ be a digraph,
356
\f$lower, upper: A\rightarrow\mathbf{Z}^+_0\f$ denote the lower and
355
Formally, let \f$G=(V,A)\f$ be a digraph, \f$lower: A\rightarrow\mathbf{Z}\f$,
356
\f$upper: A\rightarrow\mathbf{Z}\cup\{+\infty\}\f$ denote the lower and
357 357
upper bounds for the flow values on the arcs, for which
358
\f$0 \leq lower(uv) \leq upper(uv)\f$ holds for all \f$uv\in A\f$.
359
\f$cost: A\rightarrow\mathbf{Z}^+_0\f$ denotes the cost per unit flow
360
on the arcs, and \f$sup: V\rightarrow\mathbf{Z}\f$ denotes the
358
\f$lower(uv) \leq upper(uv)\f$ must hold for all \f$uv\in A\f$,
359
\f$cost: A\rightarrow\mathbf{Z}\f$ denotes the cost per unit flow
360
on the arcs and \f$sup: V\rightarrow\mathbf{Z}\f$ denotes the
361 361
signed supply values of the nodes.
362 362
If \f$sup(u)>0\f$, then \f$u\f$ is a supply node with \f$sup(u)\f$
363 363
supply, if \f$sup(u)<0\f$, then \f$u\f$ is a demand node with
364 364
\f$-sup(u)\f$ demand.
365
A minimum cost flow is an \f$f: A\rightarrow\mathbf{Z}^+_0\f$ solution
365
A minimum cost flow is an \f$f: A\rightarrow\mathbf{Z}\f$ solution
366 366
of the following optimization problem.
367 367

	
368 368
\f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
... ...
@@ -404,7 +404,7 @@
404 404

	
405 405
The dual solution of the minimum cost flow problem is represented by node 
406 406
potentials \f$\pi: V\rightarrow\mathbf{Z}\f$.
407
An \f$f: A\rightarrow\mathbf{Z}^+_0\f$ feasible solution of the problem
407
An \f$f: A\rightarrow\mathbf{Z}\f$ feasible solution of the problem
408 408
is optimal if and only if for some \f$\pi: V\rightarrow\mathbf{Z}\f$
409 409
node potentials the following \e complementary \e slackness optimality
410 410
conditions hold.
... ...
@@ -413,15 +413,15 @@
413 413
   - if \f$cost^\pi(uv)>0\f$, then \f$f(uv)=lower(uv)\f$;
414 414
   - if \f$lower(uv)<f(uv)<upper(uv)\f$, then \f$cost^\pi(uv)=0\f$;
415 415
   - if \f$cost^\pi(uv)<0\f$, then \f$f(uv)=upper(uv)\f$.
416
 - For all \f$u\in V\f$:
416
 - For all \f$u\in V\f$ nodes:
417 417
   - if \f$\sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \neq sup(u)\f$,
418 418
     then \f$\pi(u)=0\f$.
419 419
 
420 420
Here \f$cost^\pi(uv)\f$ denotes the \e reduced \e cost of the arc
421
\f$uv\in A\f$ with respect to the node potentials \f$\pi\f$, i.e.
421
\f$uv\in A\f$ with respect to the potential function \f$\pi\f$, i.e.
422 422
\f[ cost^\pi(uv) = cost(uv) + \pi(u) - \pi(v).\f]
423 423

	
424
All algorithms provide dual solution (node potentials) as well
424
All algorithms provide dual solution (node potentials) as well,
425 425
if an optimal flow is found.
426 426

	
427 427
LEMON contains several algorithms for solving minimum cost flow problems.
Show white space 6 line context
... ...
@@ -30,9 +30,6 @@
30 30

	
31 31
#include <lemon/core.h>
32 32
#include <lemon/math.h>
33
#include <lemon/maps.h>
34
#include <lemon/circulation.h>
35
#include <lemon/adaptors.h>
36 33

	
37 34
namespace lemon {
38 35

	
... ...
@@ -50,8 +47,13 @@
50 47
  ///
51 48
  /// In general this class is the fastest implementation available
52 49
  /// in LEMON for the minimum cost flow problem.
53
  /// Moreover it supports both direction of the supply/demand inequality
54
  /// constraints. For more information see \ref ProblemType.
50
  /// Moreover it supports both directions of the supply/demand inequality
51
  /// constraints. For more information see \ref SupplyType.
52
  ///
53
  /// Most of the parameters of the problem (except for the digraph)
54
  /// can be given using separate functions, and the algorithm can be
55
  /// executed using the \ref run() function. If some parameters are not
56
  /// specified, then default values will be used.
55 57
  ///
56 58
  /// \tparam GR The digraph type the algorithm runs on.
57 59
  /// \tparam F The value type used for flow amounts, capacity bounds
... ...
@@ -88,11 +90,80 @@
88 90

	
89 91
  public:
90 92

	
91
    /// \brief Enum type for selecting the pivot rule.
93
    /// \brief Problem type constants for the \c run() function.
92 94
    ///
93
    /// Enum type for selecting the pivot rule for the \ref run()
95
    /// Enum type containing the problem type constants that can be
96
    /// returned by the \ref run() function of the algorithm.
97
    enum ProblemType {
98
      /// The problem has no feasible solution (flow).
99
      INFEASIBLE,
100
      /// The problem has optimal solution (i.e. it is feasible and
101
      /// bounded), and the algorithm has found optimal flow and node
102
      /// potentials (primal and dual solutions).
103
      OPTIMAL,
104
      /// The objective function of the problem is unbounded, i.e.
105
      /// there is a directed cycle having negative total cost and
106
      /// infinite upper bound.
107
      UNBOUNDED
108
    };
109
    
110
    /// \brief Constants for selecting the type of the supply constraints.
111
    ///
112
    /// Enum type containing constants for selecting the supply type,
113
    /// i.e. the direction of the inequalities in the supply/demand
114
    /// constraints of the \ref min_cost_flow "minimum cost flow problem".
115
    ///
116
    /// The default supply type is \c GEQ, since this form is supported
117
    /// by other minimum cost flow algorithms and the \ref Circulation
118
    /// algorithm, as well.
119
    /// The \c LEQ problem type can be selected using the \ref supplyType()
94 120
    /// function.
95 121
    ///
122
    /// Note that the equality form is a special case of both supply types.
123
    enum SupplyType {
124

	
125
      /// This option means that there are <em>"greater or equal"</em>
126
      /// supply/demand constraints in the definition, i.e. the exact
127
      /// formulation of the problem is the following.
128
      /**
129
          \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
130
          \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \geq
131
              sup(u) \quad \forall u\in V \f]
132
          \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
133
      */
134
      /// It means that the total demand must be greater or equal to the 
135
      /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
136
      /// negative) and all the supplies have to be carried out from 
137
      /// the supply nodes, but there could be demands that are not 
138
      /// satisfied.
139
      GEQ,
140
      /// It is just an alias for the \c GEQ option.
141
      CARRY_SUPPLIES = GEQ,
142

	
143
      /// This option means that there are <em>"less or equal"</em>
144
      /// supply/demand constraints in the definition, i.e. the exact
145
      /// formulation of the problem is the following.
146
      /**
147
          \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
148
          \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \leq
149
              sup(u) \quad \forall u\in V \f]
150
          \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
151
      */
152
      /// It means that the total demand must be less or equal to the 
153
      /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
154
      /// positive) and all the demands have to be satisfied, but there
155
      /// could be supplies that are not carried out from the supply
156
      /// nodes.
157
      LEQ,
158
      /// It is just an alias for the \c LEQ option.
159
      SATISFY_DEMANDS = LEQ
160
    };
161
    
162
    /// \brief Constants for selecting the pivot rule.
163
    ///
164
    /// Enum type containing constants for selecting the pivot rule for
165
    /// the \ref run() function.
166
    ///
96 167
    /// \ref NetworkSimplex provides five different pivot rule
97 168
    /// implementations that significantly affect the running time
98 169
    /// of the algorithm.
... ...
@@ -131,58 +202,6 @@
131 202
      ALTERING_LIST
132 203
    };
133 204
    
134
    /// \brief Enum type for selecting the problem type.
135
    ///
136
    /// Enum type for selecting the problem type, i.e. the direction of
137
    /// the inequalities in the supply/demand constraints of the
138
    /// \ref min_cost_flow "minimum cost flow problem".
139
    ///
140
    /// The default problem type is \c GEQ, since this form is supported
141
    /// by other minimum cost flow algorithms and the \ref Circulation
142
    /// algorithm as well.
143
    /// The \c LEQ problem type can be selected using the \ref problemType()
144
    /// function.
145
    ///
146
    /// Note that the equality form is a special case of both problem type.
147
    enum ProblemType {
148

	
149
      /// This option means that there are "<em>greater or equal</em>"
150
      /// constraints in the defintion, i.e. the exact formulation of the
151
      /// problem is the following.
152
      /**
153
          \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
154
          \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \geq
155
              sup(u) \quad \forall u\in V \f]
156
          \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
157
      */
158
      /// It means that the total demand must be greater or equal to the 
159
      /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
160
      /// negative) and all the supplies have to be carried out from 
161
      /// the supply nodes, but there could be demands that are not 
162
      /// satisfied.
163
      GEQ,
164
      /// It is just an alias for the \c GEQ option.
165
      CARRY_SUPPLIES = GEQ,
166

	
167
      /// This option means that there are "<em>less or equal</em>"
168
      /// constraints in the defintion, i.e. the exact formulation of the
169
      /// problem is the following.
170
      /**
171
          \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
172
          \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \leq
173
              sup(u) \quad \forall u\in V \f]
174
          \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
175
      */
176
      /// It means that the total demand must be less or equal to the 
177
      /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
178
      /// positive) and all the demands have to be satisfied, but there
179
      /// could be supplies that are not carried out from the supply
180
      /// nodes.
181
      LEQ,
182
      /// It is just an alias for the \c LEQ option.
183
      SATISFY_DEMANDS = LEQ
184
    };
185

	
186 205
  private:
187 206

	
188 207
    TEMPLATE_DIGRAPH_TYPEDEFS(GR);
... ...
@@ -220,7 +239,9 @@
220 239
    bool _pstsup;
221 240
    Node _psource, _ptarget;
222 241
    Flow _pstflow;
223
    ProblemType _ptype;
242
    SupplyType _stype;
243
    
244
    Flow _sum_supply;
224 245

	
225 246
    // Result maps
226 247
    FlowMap *_flow_map;
... ...
@@ -259,6 +280,15 @@
259 280
    int stem, par_stem, new_stem;
260 281
    Flow delta;
261 282

	
283
  public:
284
  
285
    /// \brief Constant for infinite upper bounds (capacities).
286
    ///
287
    /// Constant for infinite upper bounds (capacities).
288
    /// It is \c std::numeric_limits<Flow>::infinity() if available,
289
    /// \c std::numeric_limits<Flow>::max() otherwise.
290
    const Flow INF;
291

	
262 292
  private:
263 293

	
264 294
    // Implementation of the First Eligible pivot rule
... ...
@@ -661,17 +691,19 @@
661 691
    NetworkSimplex(const GR& graph) :
662 692
      _graph(graph),
663 693
      _plower(NULL), _pupper(NULL), _pcost(NULL),
664
      _psupply(NULL), _pstsup(false), _ptype(GEQ),
694
      _psupply(NULL), _pstsup(false), _stype(GEQ),
665 695
      _flow_map(NULL), _potential_map(NULL),
666 696
      _local_flow(false), _local_potential(false),
667
      _node_id(graph)
697
      _node_id(graph),
698
      INF(std::numeric_limits<Flow>::has_infinity ?
699
          std::numeric_limits<Flow>::infinity() :
700
          std::numeric_limits<Flow>::max())
668 701
    {
669
      LEMON_ASSERT(std::numeric_limits<Flow>::is_integer &&
670
                   std::numeric_limits<Flow>::is_signed,
671
        "The flow type of NetworkSimplex must be signed integer");
672
      LEMON_ASSERT(std::numeric_limits<Cost>::is_integer &&
673
                   std::numeric_limits<Cost>::is_signed,
674
        "The cost type of NetworkSimplex must be signed integer");
702
      // Check the value types
703
      LEMON_ASSERT(std::numeric_limits<Flow>::is_signed,
704
        "The flow type of NetworkSimplex must be signed");
705
      LEMON_ASSERT(std::numeric_limits<Cost>::is_signed,
706
        "The cost type of NetworkSimplex must be signed");
675 707
    }
676 708

	
677 709
    /// Destructor.
... ...
@@ -689,17 +721,16 @@
689 721
    /// \brief Set the lower bounds on the arcs.
690 722
    ///
691 723
    /// This function sets the lower bounds on the arcs.
692
    /// If neither this function nor \ref boundMaps() is used before
693
    /// calling \ref run(), the lower bounds will be set to zero
694
    /// on all arcs.
724
    /// If it is not used before calling \ref run(), the lower bounds
725
    /// will be set to zero on all arcs.
695 726
    ///
696 727
    /// \param map An arc map storing the lower bounds.
697 728
    /// Its \c Value type must be convertible to the \c Flow type
698 729
    /// of the algorithm.
699 730
    ///
700 731
    /// \return <tt>(*this)</tt>
701
    template <typename LOWER>
702
    NetworkSimplex& lowerMap(const LOWER& map) {
732
    template <typename LowerMap>
733
    NetworkSimplex& lowerMap(const LowerMap& map) {
703 734
      delete _plower;
704 735
      _plower = new FlowArcMap(_graph);
705 736
      for (ArcIt a(_graph); a != INVALID; ++a) {
... ...
@@ -711,18 +742,17 @@
711 742
    /// \brief Set the upper bounds (capacities) on the arcs.
712 743
    ///
713 744
    /// This function sets the upper bounds (capacities) on the arcs.
714
    /// If none of the functions \ref upperMap(), \ref capacityMap()
715
    /// and \ref boundMaps() is used before calling \ref run(),
716
    /// the upper bounds (capacities) will be set to
717
    /// \c std::numeric_limits<Flow>::max() on all arcs.
745
    /// If it is not used before calling \ref run(), the upper bounds
746
    /// will be set to \ref INF on all arcs (i.e. the flow value will be
747
    /// unbounded from above on each arc).
718 748
    ///
719 749
    /// \param map An arc map storing the upper bounds.
720 750
    /// Its \c Value type must be convertible to the \c Flow type
721 751
    /// of the algorithm.
722 752
    ///
723 753
    /// \return <tt>(*this)</tt>
724
    template<typename UPPER>
725
    NetworkSimplex& upperMap(const UPPER& map) {
754
    template<typename UpperMap>
755
    NetworkSimplex& upperMap(const UpperMap& map) {
726 756
      delete _pupper;
727 757
      _pupper = new FlowArcMap(_graph);
728 758
      for (ArcIt a(_graph); a != INVALID; ++a) {
... ...
@@ -731,43 +761,6 @@
731 761
      return *this;
732 762
    }
733 763

	
734
    /// \brief Set the upper bounds (capacities) on the arcs.
735
    ///
736
    /// This function sets the upper bounds (capacities) on the arcs.
737
    /// It is just an alias for \ref upperMap().
738
    ///
739
    /// \return <tt>(*this)</tt>
740
    template<typename CAP>
741
    NetworkSimplex& capacityMap(const CAP& map) {
742
      return upperMap(map);
743
    }
744

	
745
    /// \brief Set the lower and upper bounds on the arcs.
746
    ///
747
    /// This function sets the lower and upper bounds on the arcs.
748
    /// If neither this function nor \ref lowerMap() is used before
749
    /// calling \ref run(), the lower bounds will be set to zero
750
    /// on all arcs.
751
    /// If none of the functions \ref upperMap(), \ref capacityMap()
752
    /// and \ref boundMaps() is used before calling \ref run(),
753
    /// the upper bounds (capacities) will be set to
754
    /// \c std::numeric_limits<Flow>::max() on all arcs.
755
    ///
756
    /// \param lower An arc map storing the lower bounds.
757
    /// \param upper An arc map storing the upper bounds.
758
    ///
759
    /// The \c Value type of the maps must be convertible to the
760
    /// \c Flow type of the algorithm.
761
    ///
762
    /// \note This function is just a shortcut of calling \ref lowerMap()
763
    /// and \ref upperMap() separately.
764
    ///
765
    /// \return <tt>(*this)</tt>
766
    template <typename LOWER, typename UPPER>
767
    NetworkSimplex& boundMaps(const LOWER& lower, const UPPER& upper) {
768
      return lowerMap(lower).upperMap(upper);
769
    }
770

	
771 764
    /// \brief Set the costs of the arcs.
772 765
    ///
773 766
    /// This function sets the costs of the arcs.
... ...
@@ -779,8 +772,8 @@
779 772
    /// of the algorithm.
780 773
    ///
781 774
    /// \return <tt>(*this)</tt>
782
    template<typename COST>
783
    NetworkSimplex& costMap(const COST& map) {
775
    template<typename CostMap>
776
    NetworkSimplex& costMap(const CostMap& map) {
784 777
      delete _pcost;
785 778
      _pcost = new CostArcMap(_graph);
786 779
      for (ArcIt a(_graph); a != INVALID; ++a) {
... ...
@@ -801,8 +794,8 @@
801 794
    /// of the algorithm.
802 795
    ///
803 796
    /// \return <tt>(*this)</tt>
804
    template<typename SUP>
805
    NetworkSimplex& supplyMap(const SUP& map) {
797
    template<typename SupplyMap>
798
    NetworkSimplex& supplyMap(const SupplyMap& map) {
806 799
      delete _psupply;
807 800
      _pstsup = false;
808 801
      _psupply = new FlowNodeMap(_graph);
... ...
@@ -820,6 +813,10 @@
820 813
    /// calling \ref run(), the supply of each node will be set to zero.
821 814
    /// (It makes sense only if non-zero lower bounds are given.)
822 815
    ///
816
    /// Using this function has the same effect as using \ref supplyMap()
817
    /// with such a map in which \c k is assigned to \c s, \c -k is
818
    /// assigned to \c t and all other nodes have zero supply value.
819
    ///
823 820
    /// \param s The source node.
824 821
    /// \param t The target node.
825 822
    /// \param k The required amount of flow from node \c s to node \c t
... ...
@@ -836,17 +833,17 @@
836 833
      return *this;
837 834
    }
838 835
    
839
    /// \brief Set the problem type.
836
    /// \brief Set the type of the supply constraints.
840 837
    ///
841
    /// This function sets the problem type for the algorithm.
842
    /// If it is not used before calling \ref run(), the \ref GEQ problem
838
    /// This function sets the type of the supply/demand constraints.
839
    /// If it is not used before calling \ref run(), the \ref GEQ supply
843 840
    /// type will be used.
844 841
    ///
845
    /// For more information see \ref ProblemType.
842
    /// For more information see \ref SupplyType.
846 843
    ///
847 844
    /// \return <tt>(*this)</tt>
848
    NetworkSimplex& problemType(ProblemType problem_type) {
849
      _ptype = problem_type;
845
    NetworkSimplex& supplyType(SupplyType supply_type) {
846
      _stype = supply_type;
850 847
      return *this;
851 848
    }
852 849

	
... ...
@@ -896,13 +893,12 @@
896 893
    ///
897 894
    /// This function runs the algorithm.
898 895
    /// The paramters can be specified using functions \ref lowerMap(),
899
    /// \ref upperMap(), \ref capacityMap(), \ref boundMaps(),
900
    /// \ref costMap(), \ref supplyMap(), \ref stSupply(), 
901
    /// \ref problemType(), \ref flowMap() and \ref potentialMap().
896
    /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(), 
897
    /// \ref supplyType(), \ref flowMap() and \ref potentialMap().
902 898
    /// For example,
903 899
    /// \code
904 900
    ///   NetworkSimplex<ListDigraph> ns(graph);
905
    ///   ns.boundMaps(lower, upper).costMap(cost)
901
    ///   ns.lowerMap(lower).upperMap(upper).costMap(cost)
906 902
    ///     .supplyMap(sup).run();
907 903
    /// \endcode
908 904
    ///
... ...
@@ -914,17 +910,25 @@
914 910
    /// \param pivot_rule The pivot rule that will be used during the
915 911
    /// algorithm. For more information see \ref PivotRule.
916 912
    ///
917
    /// \return \c true if a feasible flow can be found.
918
    bool run(PivotRule pivot_rule = BLOCK_SEARCH) {
919
      return init() && start(pivot_rule);
913
    /// \return \c INFEASIBLE if no feasible flow exists,
914
    /// \n \c OPTIMAL if the problem has optimal solution
915
    /// (i.e. it is feasible and bounded), and the algorithm has found
916
    /// optimal flow and node potentials (primal and dual solutions),
917
    /// \n \c UNBOUNDED if the objective function of the problem is
918
    /// unbounded, i.e. there is a directed cycle having negative total
919
    /// cost and infinite upper bound.
920
    ///
921
    /// \see ProblemType, PivotRule
922
    ProblemType run(PivotRule pivot_rule = BLOCK_SEARCH) {
923
      if (!init()) return INFEASIBLE;
924
      return start(pivot_rule);
920 925
    }
921 926

	
922 927
    /// \brief Reset all the parameters that have been given before.
923 928
    ///
924 929
    /// This function resets all the paramaters that have been given
925 930
    /// before using functions \ref lowerMap(), \ref upperMap(),
926
    /// \ref capacityMap(), \ref boundMaps(), \ref costMap(),
927
    /// \ref supplyMap(), \ref stSupply(), \ref problemType(), 
931
    /// \ref costMap(), \ref supplyMap(), \ref stSupply(), \ref supplyType(),
928 932
    /// \ref flowMap() and \ref potentialMap().
929 933
    ///
930 934
    /// It is useful for multiple run() calls. If this function is not
... ...
@@ -936,7 +940,7 @@
936 940
    ///   NetworkSimplex<ListDigraph> ns(graph);
937 941
    ///
938 942
    ///   // First run
939
    ///   ns.lowerMap(lower).capacityMap(cap).costMap(cost)
943
    ///   ns.lowerMap(lower).upperMap(upper).costMap(cost)
940 944
    ///     .supplyMap(sup).run();
941 945
    ///
942 946
    ///   // Run again with modified cost map (reset() is not called,
... ...
@@ -947,7 +951,7 @@
947 951
    ///   // Run again from scratch using reset()
948 952
    ///   // (the lower bounds will be set to zero on all arcs)
949 953
    ///   ns.reset();
950
    ///   ns.capacityMap(cap).costMap(cost)
954
    ///   ns.upperMap(capacity).costMap(cost)
951 955
    ///     .supplyMap(sup).run();
952 956
    /// \endcode
953 957
    ///
... ...
@@ -962,7 +966,7 @@
962 966
      _pcost = NULL;
963 967
      _psupply = NULL;
964 968
      _pstsup = false;
965
      _ptype = GEQ;
969
      _stype = GEQ;
966 970
      if (_local_flow) delete _flow_map;
967 971
      if (_local_potential) delete _potential_map;
968 972
      _flow_map = NULL;
... ...
@@ -985,7 +989,7 @@
985 989
    /// \brief Return the total cost of the found flow.
986 990
    ///
987 991
    /// This function returns the total cost of the found flow.
988
    /// The complexity of the function is O(e).
992
    /// Its complexity is O(e).
989 993
    ///
990 994
    /// \note The return type of the function can be specified as a
991 995
    /// template parameter. For example,
... ...
@@ -997,9 +1001,9 @@
997 1001
    /// function.
998 1002
    ///
999 1003
    /// \pre \ref run() must be called before using this function.
1000
    template <typename Num>
1001
    Num totalCost() const {
1002
      Num c = 0;
1004
    template <typename Value>
1005
    Value totalCost() const {
1006
      Value c = 0;
1003 1007
      if (_pcost) {
1004 1008
        for (ArcIt e(_graph); e != INVALID; ++e)
1005 1009
          c += (*_flow_map)[e] * (*_pcost)[e];
... ...
@@ -1050,7 +1054,7 @@
1050 1054
    ///
1051 1055
    /// This function returns a const reference to a node map storing
1052 1056
    /// the found potentials, which form the dual solution of the
1053
    /// \ref min_cost_flow "minimum cost flow" problem.
1057
    /// \ref min_cost_flow "minimum cost flow problem".
1054 1058
    ///
1055 1059
    /// \pre \ref run() must be called before using this function.
1056 1060
    const PotentialMap& potentialMap() const {
... ...
@@ -1101,7 +1105,7 @@
1101 1105

	
1102 1106
      // Initialize node related data
1103 1107
      bool valid_supply = true;
1104
      Flow sum_supply = 0;
1108
      _sum_supply = 0;
1105 1109
      if (!_pstsup && !_psupply) {
1106 1110
        _pstsup = true;
1107 1111
        _psource = _ptarget = NodeIt(_graph);
... ...
@@ -1112,10 +1116,10 @@
1112 1116
        for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
1113 1117
          _node_id[n] = i;
1114 1118
          _supply[i] = (*_psupply)[n];
1115
          sum_supply += _supply[i];
1119
          _sum_supply += _supply[i];
1116 1120
        }
1117
        valid_supply = (_ptype == GEQ && sum_supply <= 0) ||
1118
                       (_ptype == LEQ && sum_supply >= 0);
1121
        valid_supply = (_stype == GEQ && _sum_supply <= 0) ||
1122
                       (_stype == LEQ && _sum_supply >= 0);
1119 1123
      } else {
1120 1124
        int i = 0;
1121 1125
        for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
... ...
@@ -1127,92 +1131,18 @@
1127 1131
      }
1128 1132
      if (!valid_supply) return false;
1129 1133

	
1130
      // Infinite capacity value
1131
      Flow inf_cap =
1132
        std::numeric_limits<Flow>::has_infinity ?
1133
        std::numeric_limits<Flow>::infinity() :
1134
        std::numeric_limits<Flow>::max();
1135

	
1136 1134
      // Initialize artifical cost
1137
      Cost art_cost;
1135
      Cost ART_COST;
1138 1136
      if (std::numeric_limits<Cost>::is_exact) {
1139
        art_cost = std::numeric_limits<Cost>::max() / 4 + 1;
1137
        ART_COST = std::numeric_limits<Cost>::max() / 4 + 1;
1140 1138
      } else {
1141
        art_cost = std::numeric_limits<Cost>::min();
1139
        ART_COST = std::numeric_limits<Cost>::min();
1142 1140
        for (int i = 0; i != _arc_num; ++i) {
1143
          if (_cost[i] > art_cost) art_cost = _cost[i];
1141
          if (_cost[i] > ART_COST) ART_COST = _cost[i];
1144 1142
        }
1145
        art_cost = (art_cost + 1) * _node_num;
1143
        ART_COST = (ART_COST + 1) * _node_num;
1146 1144
      }
1147 1145

	
1148
      // Run Circulation to check if a feasible solution exists
1149
      typedef ConstMap<Arc, Flow> ConstArcMap;
1150
      ConstArcMap zero_arc_map(0), inf_arc_map(inf_cap);
1151
      FlowNodeMap *csup = NULL;
1152
      bool local_csup = false;
1153
      if (_psupply) {
1154
        csup = _psupply;
1155
      } else {
1156
        csup = new FlowNodeMap(_graph, 0);
1157
        (*csup)[_psource] =  _pstflow;
1158
        (*csup)[_ptarget] = -_pstflow;
1159
        local_csup = true;
1160
      }
1161
      bool circ_result = false;
1162
      if (_ptype == GEQ || (_ptype == LEQ && sum_supply == 0)) {
1163
        // GEQ problem type
1164
        if (_plower) {
1165
          if (_pupper) {
1166
            Circulation<GR, FlowArcMap, FlowArcMap, FlowNodeMap>
1167
              circ(_graph, *_plower, *_pupper, *csup);
1168
            circ_result = circ.run();
1169
          } else {
1170
            Circulation<GR, FlowArcMap, ConstArcMap, FlowNodeMap>
1171
              circ(_graph, *_plower, inf_arc_map, *csup);
1172
            circ_result = circ.run();
1173
          }
1174
        } else {
1175
          if (_pupper) {
1176
            Circulation<GR, ConstArcMap, FlowArcMap, FlowNodeMap>
1177
              circ(_graph, zero_arc_map, *_pupper, *csup);
1178
            circ_result = circ.run();
1179
          } else {
1180
            Circulation<GR, ConstArcMap, ConstArcMap, FlowNodeMap>
1181
              circ(_graph, zero_arc_map, inf_arc_map, *csup);
1182
            circ_result = circ.run();
1183
          }
1184
        }
1185
      } else {
1186
        // LEQ problem type
1187
        typedef ReverseDigraph<const GR> RevGraph;
1188
        typedef NegMap<FlowNodeMap> NegNodeMap;
1189
        RevGraph rgraph(_graph);
1190
        NegNodeMap neg_csup(*csup);
1191
        if (_plower) {
1192
          if (_pupper) {
1193
            Circulation<RevGraph, FlowArcMap, FlowArcMap, NegNodeMap>
1194
              circ(rgraph, *_plower, *_pupper, neg_csup);
1195
            circ_result = circ.run();
1196
          } else {
1197
            Circulation<RevGraph, FlowArcMap, ConstArcMap, NegNodeMap>
1198
              circ(rgraph, *_plower, inf_arc_map, neg_csup);
1199
            circ_result = circ.run();
1200
          }
1201
        } else {
1202
          if (_pupper) {
1203
            Circulation<RevGraph, ConstArcMap, FlowArcMap, NegNodeMap>
1204
              circ(rgraph, zero_arc_map, *_pupper, neg_csup);
1205
            circ_result = circ.run();
1206
          } else {
1207
            Circulation<RevGraph, ConstArcMap, ConstArcMap, NegNodeMap>
1208
              circ(rgraph, zero_arc_map, inf_arc_map, neg_csup);
1209
            circ_result = circ.run();
1210
          }
1211
        }
1212
      }
1213
      if (local_csup) delete csup;
1214
      if (!circ_result) return false;
1215

	
1216 1146
      // Set data for the artificial root node
1217 1147
      _root = _node_num;
1218 1148
      _parent[_root] = -1;
... ...
@@ -1221,11 +1151,11 @@
1221 1151
      _rev_thread[0] = _root;
1222 1152
      _succ_num[_root] = all_node_num;
1223 1153
      _last_succ[_root] = _root - 1;
1224
      _supply[_root] = -sum_supply;
1225
      if (sum_supply < 0) {
1226
        _pi[_root] = -art_cost;
1154
      _supply[_root] = -_sum_supply;
1155
      if (_sum_supply < 0) {
1156
        _pi[_root] = -ART_COST;
1227 1157
      } else {
1228
        _pi[_root] = art_cost;
1158
        _pi[_root] = ART_COST;
1229 1159
      }
1230 1160

	
1231 1161
      // Store the arcs in a mixed order
... ...
@@ -1260,7 +1190,7 @@
1260 1190
            _cap[i] = (*_pupper)[_arc_ref[i]];
1261 1191
        } else {
1262 1192
          for (int i = 0; i != _arc_num; ++i)
1263
            _cap[i] = inf_cap;
1193
            _cap[i] = INF;
1264 1194
        }
1265 1195
        if (_pcost) {
1266 1196
          for (int i = 0; i != _arc_num; ++i)
... ...
@@ -1275,8 +1205,17 @@
1275 1205
      if (_plower) {
1276 1206
        for (int i = 0; i != _arc_num; ++i) {
1277 1207
          Flow c = (*_plower)[_arc_ref[i]];
1278
          if (c != 0) {
1208
          if (c > 0) {
1209
            if (_cap[i] < INF) _cap[i] -= c;
1210
            _supply[_source[i]] -= c;
1211
            _supply[_target[i]] += c;
1212
          }
1213
          else if (c < 0) {
1214
            if (_cap[i] < INF + c) {
1279 1215
            _cap[i] -= c;
1216
            } else {
1217
              _cap[i] = INF;
1218
            }
1280 1219
            _supply[_source[i]] -= c;
1281 1220
            _supply[_target[i]] += c;
1282 1221
          }
... ...
@@ -1291,17 +1230,17 @@
1291 1230
        _last_succ[u] = u;
1292 1231
        _parent[u] = _root;
1293 1232
        _pred[u] = e;
1294
        _cost[e] = art_cost;
1295
        _cap[e] = inf_cap;
1233
        _cost[e] = ART_COST;
1234
        _cap[e] = INF;
1296 1235
        _state[e] = STATE_TREE;
1297
        if (_supply[u] > 0 || (_supply[u] == 0 && sum_supply <= 0)) {
1236
        if (_supply[u] > 0 || (_supply[u] == 0 && _sum_supply <= 0)) {
1298 1237
          _flow[e] = _supply[u];
1299 1238
          _forward[u] = true;
1300
          _pi[u] = -art_cost + _pi[_root];
1239
          _pi[u] = -ART_COST + _pi[_root];
1301 1240
        } else {
1302 1241
          _flow[e] = -_supply[u];
1303 1242
          _forward[u] = false;
1304
          _pi[u] = art_cost + _pi[_root];
1243
          _pi[u] = ART_COST + _pi[_root];
1305 1244
        }
1306 1245
      }
1307 1246

	
... ...
@@ -1342,7 +1281,8 @@
1342 1281
      // Search the cycle along the path form the first node to the root
1343 1282
      for (int u = first; u != join; u = _parent[u]) {
1344 1283
        e = _pred[u];
1345
        d = _forward[u] ? _flow[e] : _cap[e] - _flow[e];
1284
        d = _forward[u] ?
1285
          _flow[e] : (_cap[e] == INF ? INF : _cap[e] - _flow[e]);
1346 1286
        if (d < delta) {
1347 1287
          delta = d;
1348 1288
          u_out = u;
... ...
@@ -1352,7 +1292,8 @@
1352 1292
      // Search the cycle along the path form the second node to the root
1353 1293
      for (int u = second; u != join; u = _parent[u]) {
1354 1294
        e = _pred[u];
1355
        d = _forward[u] ? _cap[e] - _flow[e] : _flow[e];
1295
        d = _forward[u] ? 
1296
          (_cap[e] == INF ? INF : _cap[e] - _flow[e]) : _flow[e];
1356 1297
        if (d <= delta) {
1357 1298
          delta = d;
1358 1299
          u_out = u;
... ...
@@ -1526,7 +1467,7 @@
1526 1467
    }
1527 1468

	
1528 1469
    // Execute the algorithm
1529
    bool start(PivotRule pivot_rule) {
1470
    ProblemType start(PivotRule pivot_rule) {
1530 1471
      // Select the pivot rule implementation
1531 1472
      switch (pivot_rule) {
1532 1473
        case FIRST_ELIGIBLE:
... ...
@@ -1540,17 +1481,18 @@
1540 1481
        case ALTERING_LIST:
1541 1482
          return start<AlteringListPivotRule>();
1542 1483
      }
1543
      return false;
1484
      return INFEASIBLE; // avoid warning
1544 1485
    }
1545 1486

	
1546 1487
    template <typename PivotRuleImpl>
1547
    bool start() {
1488
    ProblemType start() {
1548 1489
      PivotRuleImpl pivot(*this);
1549 1490

	
1550 1491
      // Execute the Network Simplex algorithm
1551 1492
      while (pivot.findEnteringArc()) {
1552 1493
        findJoinNode();
1553 1494
        bool change = findLeavingArc();
1495
        if (delta >= INF) return UNBOUNDED;
1554 1496
        changeFlow(change);
1555 1497
        if (change) {
1556 1498
          updateTreeStructure();
... ...
@@ -1558,6 +1500,23 @@
1558 1500
        }
1559 1501
      }
1560 1502

	
1503
      // Check feasibility
1504
      if (_sum_supply < 0) {
1505
        for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
1506
          if (_supply[u] >= 0 && _flow[e] != 0) return INFEASIBLE;
1507
        }
1508
      }
1509
      else if (_sum_supply > 0) {
1510
        for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
1511
          if (_supply[u] <= 0 && _flow[e] != 0) return INFEASIBLE;
1512
        }
1513
      }
1514
      else {
1515
        for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
1516
          if (_flow[e] != 0) return INFEASIBLE;
1517
        }
1518
      }
1519

	
1561 1520
      // Copy flow values to _flow_map
1562 1521
      if (_plower) {
1563 1522
        for (int i = 0; i != _arc_num; ++i) {
... ...
@@ -1574,7 +1533,7 @@
1574 1533
        _potential_map->set(n, _pi[_node_id[n]]);
1575 1534
      }
1576 1535

	
1577
      return true;
1536
      return OPTIMAL;
1578 1537
    }
1579 1538

	
1580 1539
  }; //class NetworkSimplex
Show white space 6 line context
... ...
@@ -18,6 +18,7 @@
18 18

	
19 19
#include <iostream>
20 20
#include <fstream>
21
#include <limits>
21 22

	
22 23
#include <lemon/list_graph.h>
23 24
#include <lemon/lgf_reader.h>
... ...
@@ -33,50 +34,50 @@
33 34

	
34 35
char test_lgf[] =
35 36
  "@nodes\n"
36
  "label  sup1 sup2 sup3 sup4 sup5\n"
37
  "    1    20   27    0   20   30\n"
38
  "    2    -4    0    0   -8   -3\n"
39
  "    3     0    0    0    0    0\n"
40
  "    4     0    0    0    0    0\n"
41
  "    5     9    0    0    6   11\n"
42
  "    6    -6    0    0   -5   -6\n"
43
  "    7     0    0    0    0    0\n"
44
  "    8     0    0    0    0    3\n"
45
  "    9     3    0    0    0    0\n"
46
  "   10    -2    0    0   -7   -2\n"
47
  "   11     0    0    0  -10    0\n"
48
  "   12   -20  -27    0  -30  -20\n"
37
  "label  sup1 sup2 sup3 sup4 sup5 sup6\n"
38
  "    1    20   27    0   30   20   30\n"
39
  "    2    -4    0    0    0   -8   -3\n"
40
  "    3     0    0    0    0    0    0\n"
41
  "    4     0    0    0    0    0    0\n"
42
  "    5     9    0    0    0    6   11\n"
43
  "    6    -6    0    0    0   -5   -6\n"
44
  "    7     0    0    0    0    0    0\n"
45
  "    8     0    0    0    0    0    3\n"
46
  "    9     3    0    0    0    0    0\n"
47
  "   10    -2    0    0    0   -7   -2\n"
48
  "   11     0    0    0    0  -10    0\n"
49
  "   12   -20  -27    0  -30  -30  -20\n"
49 50
  "\n"
50 51
  "@arcs\n"
51
  "       cost  cap low1 low2\n"
52
  " 1  2    70   11    0    8\n"
53
  " 1  3   150    3    0    1\n"
54
  " 1  4    80   15    0    2\n"
55
  " 2  8    80   12    0    0\n"
56
  " 3  5   140    5    0    3\n"
57
  " 4  6    60   10    0    1\n"
58
  " 4  7    80    2    0    0\n"
59
  " 4  8   110    3    0    0\n"
60
  " 5  7    60   14    0    0\n"
61
  " 5 11   120   12    0    0\n"
62
  " 6  3     0    3    0    0\n"
63
  " 6  9   140    4    0    0\n"
64
  " 6 10    90    8    0    0\n"
65
  " 7  1    30    5    0    0\n"
66
  " 8 12    60   16    0    4\n"
67
  " 9 12    50    6    0    0\n"
68
  "10 12    70   13    0    5\n"
69
  "10  2   100    7    0    0\n"
70
  "10  7    60   10    0    0\n"
71
  "11 10    20   14    0    6\n"
72
  "12 11    30   10    0    0\n"
52
  "       cost  cap low1 low2 low3\n"
53
  " 1  2    70   11    0    8    8\n"
54
  " 1  3   150    3    0    1    0\n"
55
  " 1  4    80   15    0    2    2\n"
56
  " 2  8    80   12    0    0    0\n"
57
  " 3  5   140    5    0    3    1\n"
58
  " 4  6    60   10    0    1    0\n"
59
  " 4  7    80    2    0    0    0\n"
60
  " 4  8   110    3    0    0    0\n"
61
  " 5  7    60   14    0    0    0\n"
62
  " 5 11   120   12    0    0    0\n"
63
  " 6  3     0    3    0    0    0\n"
64
  " 6  9   140    4    0    0    0\n"
65
  " 6 10    90    8    0    0    0\n"
66
  " 7  1    30    5    0    0   -5\n"
67
  " 8 12    60   16    0    4    3\n"
68
  " 9 12    50    6    0    0    0\n"
69
  "10 12    70   13    0    5    2\n"
70
  "10  2   100    7    0    0    0\n"
71
  "10  7    60   10    0    0   -3\n"
72
  "11 10    20   14    0    6  -20\n"
73
  "12 11    30   10    0    0  -10\n"
73 74
  "\n"
74 75
  "@attributes\n"
75 76
  "source 1\n"
76 77
  "target 12\n";
77 78

	
78 79

	
79
enum ProblemType {
80
enum SupplyType {
80 81
  EQ,
81 82
  GEQ,
82 83
  LEQ
... ...
@@ -98,8 +99,6 @@
98 99
      b = mcf.reset()
99 100
             .lowerMap(lower)
100 101
             .upperMap(upper)
101
             .capacityMap(upper)
102
             .boundMaps(lower, upper)
103 102
             .costMap(cost)
104 103
             .supplyMap(sup)
105 104
             .stSupply(n, n, k)
... ...
@@ -112,10 +111,12 @@
112 111
      const typename MCF::FlowMap &fm = const_mcf.flowMap();
113 112
      const typename MCF::PotentialMap &pm = const_mcf.potentialMap();
114 113

	
115
      v = const_mcf.totalCost();
114
      c = const_mcf.totalCost();
116 115
      double x = const_mcf.template totalCost<double>();
117 116
      v = const_mcf.flow(a);
118
      v = const_mcf.potential(n);
117
      c = const_mcf.potential(n);
118
      
119
      v = const_mcf.INF;
119 120

	
120 121
      ignore_unused_variable_warning(fm);
121 122
      ignore_unused_variable_warning(pm);
... ...
@@ -137,6 +138,7 @@
137 138
    const Arc &a;
138 139
    const Flow &k;
139 140
    Flow v;
141
    Cost c;
140 142
    bool b;
141 143

	
142 144
    typename MCF::FlowMap &flow;
... ...
@@ -151,7 +153,7 @@
151 153
           typename SM, typename FM >
152 154
bool checkFlow( const GR& gr, const LM& lower, const UM& upper,
153 155
                const SM& supply, const FM& flow,
154
                ProblemType type = EQ )
156
                SupplyType type = EQ )
155 157
{
156 158
  TEMPLATE_DIGRAPH_TYPEDEFS(GR);
157 159

	
... ...
@@ -208,16 +210,17 @@
208 210
// Run a minimum cost flow algorithm and check the results
209 211
template < typename MCF, typename GR,
210 212
           typename LM, typename UM,
211
           typename CM, typename SM >
212
void checkMcf( const MCF& mcf, bool mcf_result,
213
           typename CM, typename SM,
214
           typename PT >
215
void checkMcf( const MCF& mcf, PT mcf_result,
213 216
               const GR& gr, const LM& lower, const UM& upper,
214 217
               const CM& cost, const SM& supply,
215
               bool result, typename CM::Value total,
218
               PT result, bool optimal, typename CM::Value total,
216 219
               const std::string &test_id = "",
217
               ProblemType type = EQ )
220
               SupplyType type = EQ )
218 221
{
219 222
  check(mcf_result == result, "Wrong result " + test_id);
220
  if (result) {
223
  if (optimal) {
221 224
    check(checkFlow(gr, lower, upper, supply, mcf.flowMap(), type),
222 225
          "The flow is not feasible " + test_id);
223 226
    check(mcf.totalCost() == total, "The flow is not optimal " + test_id);
... ...
@@ -244,8 +247,8 @@
244 247

	
245 248
  // Read the test digraph
246 249
  Digraph gr;
247
  Digraph::ArcMap<int> c(gr), l1(gr), l2(gr), u(gr);
248
  Digraph::NodeMap<int> s1(gr), s2(gr), s3(gr), s4(gr), s5(gr);
250
  Digraph::ArcMap<int> c(gr), l1(gr), l2(gr), l3(gr), u(gr);
251
  Digraph::NodeMap<int> s1(gr), s2(gr), s3(gr), s4(gr), s5(gr), s6(gr);
249 252
  ConstMap<Arc, int> cc(1), cu(std::numeric_limits<int>::max());
250 253
  Node v, w;
251 254

	
... ...
@@ -255,15 +258,57 @@
255 258
    .arcMap("cap", u)
256 259
    .arcMap("low1", l1)
257 260
    .arcMap("low2", l2)
261
    .arcMap("low3", l3)
258 262
    .nodeMap("sup1", s1)
259 263
    .nodeMap("sup2", s2)
260 264
    .nodeMap("sup3", s3)
261 265
    .nodeMap("sup4", s4)
262 266
    .nodeMap("sup5", s5)
267
    .nodeMap("sup6", s6)
263 268
    .node("source", v)
264 269
    .node("target", w)
265 270
    .run();
266 271

	
272
  // Build a test digraph for testing negative costs
273
  Digraph ngr;
274
  Node n1 = ngr.addNode();
275
  Node n2 = ngr.addNode();
276
  Node n3 = ngr.addNode();
277
  Node n4 = ngr.addNode();
278
  Node n5 = ngr.addNode();
279
  Node n6 = ngr.addNode();
280
  Node n7 = ngr.addNode();
281
  
282
  Arc a1 = ngr.addArc(n1, n2);
283
  Arc a2 = ngr.addArc(n1, n3);
284
  Arc a3 = ngr.addArc(n2, n4);
285
  Arc a4 = ngr.addArc(n3, n4);
286
  Arc a5 = ngr.addArc(n3, n2);
287
  Arc a6 = ngr.addArc(n5, n3);
288
  Arc a7 = ngr.addArc(n5, n6);
289
  Arc a8 = ngr.addArc(n6, n7);
290
  Arc a9 = ngr.addArc(n7, n5);
291
  
292
  Digraph::ArcMap<int> nc(ngr), nl1(ngr, 0), nl2(ngr, 0);
293
  ConstMap<Arc, int> nu1(std::numeric_limits<int>::max()), nu2(5000);
294
  Digraph::NodeMap<int> ns(ngr, 0);
295
  
296
  nl2[a7] =  1000;
297
  nl2[a8] = -1000;
298
  
299
  ns[n1] =  100;
300
  ns[n4] = -100;
301
  
302
  nc[a1] =  100;
303
  nc[a2] =   30;
304
  nc[a3] =   20;
305
  nc[a4] =   80;
306
  nc[a5] =   50;
307
  nc[a6] =   10;
308
  nc[a7] =   80;
309
  nc[a8] =   30;
310
  nc[a9] = -120;
311

	
267 312
  // A. Test NetworkSimplex with the default pivot rule
268 313
  {
269 314
    NetworkSimplex<Digraph> mcf(gr);
... ...
@@ -271,63 +316,77 @@
271 316
    // Check the equality form
272 317
    mcf.upperMap(u).costMap(c);
273 318
    checkMcf(mcf, mcf.supplyMap(s1).run(),
274
             gr, l1, u, c, s1, true,  5240, "#A1");
319
             gr, l1, u, c, s1, mcf.OPTIMAL, true,   5240, "#A1");
275 320
    checkMcf(mcf, mcf.stSupply(v, w, 27).run(),
276
             gr, l1, u, c, s2, true,  7620, "#A2");
321
             gr, l1, u, c, s2, mcf.OPTIMAL, true,   7620, "#A2");
277 322
    mcf.lowerMap(l2);
278 323
    checkMcf(mcf, mcf.supplyMap(s1).run(),
279
             gr, l2, u, c, s1, true,  5970, "#A3");
324
             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#A3");
280 325
    checkMcf(mcf, mcf.stSupply(v, w, 27).run(),
281
             gr, l2, u, c, s2, true,  8010, "#A4");
326
             gr, l2, u, c, s2, mcf.OPTIMAL, true,   8010, "#A4");
282 327
    mcf.reset();
283 328
    checkMcf(mcf, mcf.supplyMap(s1).run(),
284
             gr, l1, cu, cc, s1, true,  74, "#A5");
329
             gr, l1, cu, cc, s1, mcf.OPTIMAL, true,   74, "#A5");
285 330
    checkMcf(mcf, mcf.lowerMap(l2).stSupply(v, w, 27).run(),
286
             gr, l2, cu, cc, s2, true,  94, "#A6");
331
             gr, l2, cu, cc, s2, mcf.OPTIMAL, true,   94, "#A6");
287 332
    mcf.reset();
288 333
    checkMcf(mcf, mcf.run(),
289
             gr, l1, cu, cc, s3, true,   0, "#A7");
290
    checkMcf(mcf, mcf.boundMaps(l2, u).run(),
291
             gr, l2, u, cc, s3, false,   0, "#A8");
334
             gr, l1, cu, cc, s3, mcf.OPTIMAL, true,    0, "#A7");
335
    checkMcf(mcf, mcf.lowerMap(l2).upperMap(u).run(),
336
             gr, l2, u, cc, s3, mcf.INFEASIBLE, false, 0, "#A8");
337
    mcf.reset().lowerMap(l3).upperMap(u).costMap(c).supplyMap(s4);
338
    checkMcf(mcf, mcf.run(),
339
             gr, l3, u, c, s4, mcf.OPTIMAL, true,   6360, "#A9");
292 340

	
293 341
    // Check the GEQ form
294
    mcf.reset().upperMap(u).costMap(c).supplyMap(s4);
342
    mcf.reset().upperMap(u).costMap(c).supplyMap(s5);
295 343
    checkMcf(mcf, mcf.run(),
296
             gr, l1, u, c, s4, true,  3530, "#A9", GEQ);
297
    mcf.problemType(mcf.GEQ);
344
             gr, l1, u, c, s5, mcf.OPTIMAL, true,   3530, "#A10", GEQ);
345
    mcf.supplyType(mcf.GEQ);
298 346
    checkMcf(mcf, mcf.lowerMap(l2).run(),
299
             gr, l2, u, c, s4, true,  4540, "#A10", GEQ);
300
    mcf.problemType(mcf.CARRY_SUPPLIES).supplyMap(s5);
347
             gr, l2, u, c, s5, mcf.OPTIMAL, true,   4540, "#A11", GEQ);
348
    mcf.supplyType(mcf.CARRY_SUPPLIES).supplyMap(s6);
301 349
    checkMcf(mcf, mcf.run(),
302
             gr, l2, u, c, s5, false,    0, "#A11", GEQ);
350
             gr, l2, u, c, s6, mcf.INFEASIBLE, false,  0, "#A12", GEQ);
303 351

	
304 352
    // Check the LEQ form
305
    mcf.reset().problemType(mcf.LEQ);
306
    mcf.upperMap(u).costMap(c).supplyMap(s5);
353
    mcf.reset().supplyType(mcf.LEQ);
354
    mcf.upperMap(u).costMap(c).supplyMap(s6);
307 355
    checkMcf(mcf, mcf.run(),
308
             gr, l1, u, c, s5, true,  5080, "#A12", LEQ);
356
             gr, l1, u, c, s6, mcf.OPTIMAL, true,   5080, "#A13", LEQ);
309 357
    checkMcf(mcf, mcf.lowerMap(l2).run(),
310
             gr, l2, u, c, s5, true,  5930, "#A13", LEQ);
311
    mcf.problemType(mcf.SATISFY_DEMANDS).supplyMap(s4);
358
             gr, l2, u, c, s6, mcf.OPTIMAL, true,   5930, "#A14", LEQ);
359
    mcf.supplyType(mcf.SATISFY_DEMANDS).supplyMap(s5);
312 360
    checkMcf(mcf, mcf.run(),
313
             gr, l2, u, c, s4, false,    0, "#A14", LEQ);
361
             gr, l2, u, c, s5, mcf.INFEASIBLE, false,  0, "#A15", LEQ);
362

	
363
    // Check negative costs
364
    NetworkSimplex<Digraph> nmcf(ngr);
365
    nmcf.lowerMap(nl1).costMap(nc).supplyMap(ns);
366
    checkMcf(nmcf, nmcf.run(),
367
      ngr, nl1, nu1, nc, ns, nmcf.UNBOUNDED, false,    0, "#A16");
368
    checkMcf(nmcf, nmcf.upperMap(nu2).run(),
369
      ngr, nl1, nu2, nc, ns, nmcf.OPTIMAL, true,  -40000, "#A17");
370
    nmcf.reset().lowerMap(nl2).costMap(nc).supplyMap(ns);
371
    checkMcf(nmcf, nmcf.run(),
372
      ngr, nl2, nu1, nc, ns, nmcf.UNBOUNDED, false,    0, "#A18");
314 373
  }
315 374

	
316 375
  // B. Test NetworkSimplex with each pivot rule
317 376
  {
318 377
    NetworkSimplex<Digraph> mcf(gr);
319
    mcf.supplyMap(s1).costMap(c).capacityMap(u).lowerMap(l2);
378
    mcf.supplyMap(s1).costMap(c).upperMap(u).lowerMap(l2);
320 379

	
321 380
    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::FIRST_ELIGIBLE),
322
             gr, l2, u, c, s1, true,  5970, "#B1");
381
             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B1");
323 382
    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::BEST_ELIGIBLE),
324
             gr, l2, u, c, s1, true,  5970, "#B2");
383
             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B2");
325 384
    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::BLOCK_SEARCH),
326
             gr, l2, u, c, s1, true,  5970, "#B3");
385
             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B3");
327 386
    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::CANDIDATE_LIST),
328
             gr, l2, u, c, s1, true,  5970, "#B4");
387
             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B4");
329 388
    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::ALTERING_LIST),
330
             gr, l2, u, c, s1, true,  5970, "#B5");
389
             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B5");
331 390
  }
332 391

	
333 392
  return 0;
Show white space 6 line context
... ...
@@ -119,8 +119,8 @@
119 119

	
120 120
  ti.restart();
121 121
  NetworkSimplex<Digraph, Value> ns(g);
122
  ns.lowerMap(lower).capacityMap(cap).costMap(cost).supplyMap(sup);
123
  if (sum_sup > 0) ns.problemType(ns.LEQ);
122
  ns.lowerMap(lower).upperMap(cap).costMap(cost).supplyMap(sup);
123
  if (sum_sup > 0) ns.supplyType(ns.LEQ);
124 124
  if (report) std::cerr << "Setup NetworkSimplex class: " << ti << '\n';
125 125
  ti.restart();
126 126
  bool res = ns.run();
0 comments (0 inline)