gravatar
kpeter (Peter Kovacs)
kpeter@inf.elte.hu
Support multiple run() calls in NetworkSimplex (#234)
0 2 0
default
2 files changed with 90 insertions and 31 deletions:
↑ Collapse diff ↑
Ignore white space 24 line context
... ...
@@ -32,30 +32,36 @@
32 32
#include <lemon/math.h>
33 33

	
34 34
namespace lemon {
35 35

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

	
39 39
  /// \brief Implementation of the primal Network Simplex algorithm
40 40
  /// for finding a \ref min_cost_flow "minimum cost flow".
41 41
  ///
42 42
  /// \ref NetworkSimplex implements the primal Network Simplex algorithm
43 43
  /// for finding a \ref min_cost_flow "minimum cost flow".
44
  /// This algorithm is a specialized version of the linear programming
45
  /// simplex method directly for the minimum cost flow problem.
46
  /// It is one of the most efficient solution methods.
47
  ///
48
  /// In general this class is the fastest implementation available
49
  /// in LEMON for the minimum cost flow problem.
44 50
  ///
45 51
  /// \tparam GR The digraph type the algorithm runs on.
46 52
  /// \tparam V The value type used in the algorithm.
47 53
  /// By default it is \c int.
48 54
  ///
49
  /// \warning \c V must be a signed integer type.
55
  /// \warning The value type must be a signed integer type.
50 56
  ///
51 57
  /// \note %NetworkSimplex provides five different pivot rule
52 58
  /// implementations. For more information see \ref PivotRule.
53 59
  template <typename GR, typename V = int>
54 60
  class NetworkSimplex
55 61
  {
56 62
  public:
57 63

	
58 64
    /// The value type of the algorithm
59 65
    typedef V Value;
60 66
    /// The type of the flow map
61 67
    typedef typename GR::template ArcMap<Value> FlowMap;
... ...
@@ -780,41 +786,91 @@
780 786
      return *this;
781 787
    }
782 788

	
783 789
    /// \name Execution Control
784 790
    /// The algorithm can be executed using \ref run().
785 791

	
786 792
    /// @{
787 793

	
788 794
    /// \brief Run the algorithm.
789 795
    ///
790 796
    /// This function runs the algorithm.
791 797
    /// The paramters can be specified using \ref lowerMap(),
792
    /// \ref upperMap(), \ref capacityMap(), \ref boundMaps(), 
798
    /// \ref upperMap(), \ref capacityMap(), \ref boundMaps(),
793 799
    /// \ref costMap(), \ref supplyMap() and \ref stSupply()
794 800
    /// functions. For example,
795 801
    /// \code
796 802
    ///   NetworkSimplex<ListDigraph> ns(graph);
797 803
    ///   ns.boundMaps(lower, upper).costMap(cost)
798 804
    ///     .supplyMap(sup).run();
799 805
    /// \endcode
800 806
    ///
807
    /// This function can be called more than once. All the parameters
808
    /// that have been given are kept for the next call, unless
809
    /// \ref reset() is called, thus only the modified parameters
810
    /// have to be set again. See \ref reset() for examples.
811
    ///
801 812
    /// \param pivot_rule The pivot rule that will be used during the
802 813
    /// algorithm. For more information see \ref PivotRule.
803 814
    ///
804 815
    /// \return \c true if a feasible flow can be found.
805 816
    bool run(PivotRule pivot_rule = BLOCK_SEARCH) {
806 817
      return init() && start(pivot_rule);
807 818
    }
808 819

	
820
    /// \brief Reset all the parameters that have been given before.
821
    ///
822
    /// This function resets all the paramaters that have been given
823
    /// using \ref lowerMap(), \ref upperMap(), \ref capacityMap(),
824
    /// \ref boundMaps(), \ref costMap(), \ref supplyMap() and
825
    /// \ref stSupply() functions before.
826
    ///
827
    /// It is useful for multiple run() calls. If this function is not
828
    /// used, all the parameters given before are kept for the next
829
    /// \ref run() call.
830
    ///
831
    /// For example,
832
    /// \code
833
    ///   NetworkSimplex<ListDigraph> ns(graph);
834
    ///
835
    ///   // First run
836
    ///   ns.lowerMap(lower).capacityMap(cap).costMap(cost)
837
    ///     .supplyMap(sup).run();
838
    ///
839
    ///   // Run again with modified cost map (reset() is not called,
840
    ///   // so only the cost map have to be set again)
841
    ///   cost[e] += 100;
842
    ///   ns.costMap(cost).run();
843
    ///
844
    ///   // Run again from scratch using reset()
845
    ///   // (the lower bounds will be set to zero on all arcs)
846
    ///   ns.reset();
847
    ///   ns.capacityMap(cap).costMap(cost)
848
    ///     .supplyMap(sup).run();
849
    /// \endcode
850
    ///
851
    /// \return <tt>(*this)</tt>
852
    NetworkSimplex& reset() {
853
      delete _plower;
854
      delete _pupper;
855
      delete _pcost;
856
      delete _psupply;
857
      _plower = NULL;
858
      _pupper = NULL;
859
      _pcost = NULL;
860
      _psupply = NULL;
861
      _pstsup = false;
862
      return *this;
863
    }
864

	
809 865
    /// @}
810 866

	
811 867
    /// \name Query Functions
812 868
    /// The results of the algorithm can be obtained using these
813 869
    /// functions.\n
814 870
    /// The \ref run() function must be called before using them.
815 871

	
816 872
    /// @{
817 873

	
818 874
    /// \brief Return the total cost of the found flow.
819 875
    ///
820 876
    /// This function returns the total cost of the found flow.
... ...
@@ -911,35 +967,35 @@
911 967
      _arc_num = countArcs(_graph);
912 968
      int all_node_num = _node_num + 1;
913 969
      int all_arc_num = _arc_num + _node_num;
914 970
      if (_node_num == 0) return false;
915 971

	
916 972
      _arc_ref.resize(_arc_num);
917 973
      _source.resize(all_arc_num);
918 974
      _target.resize(all_arc_num);
919 975

	
920 976
      _cap.resize(all_arc_num);
921 977
      _cost.resize(all_arc_num);
922 978
      _supply.resize(all_node_num);
923
      _flow.resize(all_arc_num, 0);
924
      _pi.resize(all_node_num, 0);
979
      _flow.resize(all_arc_num);
980
      _pi.resize(all_node_num);
925 981

	
926 982
      _parent.resize(all_node_num);
927 983
      _pred.resize(all_node_num);
928 984
      _forward.resize(all_node_num);
929 985
      _thread.resize(all_node_num);
930 986
      _rev_thread.resize(all_node_num);
931 987
      _succ_num.resize(all_node_num);
932 988
      _last_succ.resize(all_node_num);
933
      _state.resize(all_arc_num, STATE_LOWER);
989
      _state.resize(all_arc_num);
934 990

	
935 991
      // Initialize node related data
936 992
      bool valid_supply = true;
937 993
      if (!_pstsup && !_psupply) {
938 994
        _pstsup = true;
939 995
        _psource = _ptarget = NodeIt(_graph);
940 996
        _pstflow = 0;
941 997
      }
942 998
      if (_psupply) {
943 999
        Value sum = 0;
944 1000
        int i = 0;
945 1001
        for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
... ...
@@ -977,30 +1033,34 @@
977 1033
        _arc_ref[i] = e;
978 1034
        if ((i += k) >= _arc_num) i = (i % k) + 1;
979 1035
      }
980 1036

	
981 1037
      // Initialize arc maps
982 1038
      if (_pupper && _pcost) {
983 1039
        for (int i = 0; i != _arc_num; ++i) {
984 1040
          Arc e = _arc_ref[i];
985 1041
          _source[i] = _node_id[_graph.source(e)];
986 1042
          _target[i] = _node_id[_graph.target(e)];
987 1043
          _cap[i] = (*_pupper)[e];
988 1044
          _cost[i] = (*_pcost)[e];
1045
          _flow[i] = 0;
1046
          _state[i] = STATE_LOWER;
989 1047
        }
990 1048
      } else {
991 1049
        for (int i = 0; i != _arc_num; ++i) {
992 1050
          Arc e = _arc_ref[i];
993 1051
          _source[i] = _node_id[_graph.source(e)];
994 1052
          _target[i] = _node_id[_graph.target(e)];
1053
          _flow[i] = 0;
1054
          _state[i] = STATE_LOWER;
995 1055
        }
996 1056
        if (_pupper) {
997 1057
          for (int i = 0; i != _arc_num; ++i)
998 1058
            _cap[i] = (*_pupper)[_arc_ref[i]];
999 1059
        } else {
1000 1060
          Value val = std::numeric_limits<Value>::max();
1001 1061
          for (int i = 0; i != _arc_num; ++i)
1002 1062
            _cap[i] = val;
1003 1063
        }
1004 1064
        if (_pcost) {
1005 1065
          for (int i = 0; i != _arc_num; ++i)
1006 1066
            _cost[i] = (*_pcost)[_arc_ref[i]];
... ...
@@ -1023,36 +1083,36 @@
1023 1083
      }
1024 1084

	
1025 1085
      // Add artificial arcs and initialize the spanning tree data structure
1026 1086
      Value max_cap = std::numeric_limits<Value>::max();
1027 1087
      Value max_cost = std::numeric_limits<Value>::max() / 4;
1028 1088
      for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
1029 1089
        _thread[u] = u + 1;
1030 1090
        _rev_thread[u + 1] = u;
1031 1091
        _succ_num[u] = 1;
1032 1092
        _last_succ[u] = u;
1033 1093
        _parent[u] = _root;
1034 1094
        _pred[u] = e;
1095
        _cost[e] = max_cost;
1096
        _cap[e] = max_cap;
1097
        _state[e] = STATE_TREE;
1035 1098
        if (_supply[u] >= 0) {
1036 1099
          _flow[e] = _supply[u];
1037 1100
          _forward[u] = true;
1038 1101
          _pi[u] = -max_cost;
1039 1102
        } else {
1040 1103
          _flow[e] = -_supply[u];
1041 1104
          _forward[u] = false;
1042 1105
          _pi[u] = max_cost;
1043 1106
        }
1044
        _cost[e] = max_cost;
1045
        _cap[e] = max_cap;
1046
        _state[e] = STATE_TREE;
1047 1107
      }
1048 1108

	
1049 1109
      return true;
1050 1110
    }
1051 1111

	
1052 1112
    // Find the join node
1053 1113
    void findJoinNode() {
1054 1114
      int u = _source[in_arc];
1055 1115
      int v = _target[in_arc];
1056 1116
      while (u != v) {
1057 1117
        if (_succ_num[u] < _succ_num[v]) {
1058 1118
          u = _parent[u];
Ignore white space 6 line context
... ...
@@ -80,25 +80,26 @@
80 80
template <typename GR, typename Value>
81 81
class McfClassConcept
82 82
{
83 83
public:
84 84

	
85 85
  template <typename MCF>
86 86
  struct Constraints {
87 87
    void constraints() {
88 88
      checkConcept<concepts::Digraph, GR>();
89 89

	
90 90
      MCF mcf(g);
91 91

	
92
      b = mcf.lowerMap(lower)
92
      b = mcf.reset()
93
             .lowerMap(lower)
93 94
             .upperMap(upper)
94 95
             .capacityMap(upper)
95 96
             .boundMaps(lower, upper)
96 97
             .costMap(cost)
97 98
             .supplyMap(sup)
98 99
             .stSupply(n, n, k)
99 100
             .run();
100 101

	
101 102
      const typename MCF::FlowMap &fm = mcf.flowMap();
102 103
      const typename MCF::PotentialMap &pm = mcf.potentialMap();
103 104

	
104 105
      v = mcf.totalCost();
... ...
@@ -233,57 +234,55 @@
233 234
    .arcMap("cap", u)
234 235
    .arcMap("low1", l1)
235 236
    .arcMap("low2", l2)
236 237
    .nodeMap("sup1", s1)
237 238
    .nodeMap("sup2", s2)
238 239
    .nodeMap("sup3", s3)
239 240
    .node("source", v)
240 241
    .node("target", w)
241 242
    .run();
242 243

	
243 244
  // A. Test NetworkSimplex with the default pivot rule
244 245
  {
245
    NetworkSimplex<Digraph> mcf1(gr), mcf2(gr), mcf3(gr), mcf4(gr),
246
                            mcf5(gr), mcf6(gr), mcf7(gr), mcf8(gr);
246
    NetworkSimplex<Digraph> mcf(gr);
247 247

	
248
    checkMcf(mcf1, mcf1.upperMap(u).costMap(c).supplyMap(s1).run(),
248
    mcf.upperMap(u).costMap(c);
249
    checkMcf(mcf, mcf.supplyMap(s1).run(),
249 250
             gr, l1, u, c, s1, true,  5240, "#A1");
250
    checkMcf(mcf2, mcf2.upperMap(u).costMap(c).stSupply(v, w, 27).run(),
251
    checkMcf(mcf, mcf.stSupply(v, w, 27).run(),
251 252
             gr, l1, u, c, s2, true,  7620, "#A2");
252
    checkMcf(mcf3, mcf3.boundMaps(l2, u).costMap(c).supplyMap(s1).run(),
253
    mcf.lowerMap(l2);
254
    checkMcf(mcf, mcf.supplyMap(s1).run(),
253 255
             gr, l2, u, c, s1, true,  5970, "#A3");
254
    checkMcf(mcf4, mcf4.boundMaps(l2, u).costMap(c).stSupply(v, w, 27).run(),
256
    checkMcf(mcf, mcf.stSupply(v, w, 27).run(),
255 257
             gr, l2, u, c, s2, true,  8010, "#A4");
256
    checkMcf(mcf5, mcf5.supplyMap(s1).run(),
258
    mcf.reset();
259
    checkMcf(mcf, mcf.supplyMap(s1).run(),
257 260
             gr, l1, cu, cc, s1, true,  74, "#A5");
258
    checkMcf(mcf6, mcf6.stSupply(v, w, 27).lowerMap(l2).run(),
261
    checkMcf(mcf, mcf.lowerMap(l2).stSupply(v, w, 27).run(),
259 262
             gr, l2, cu, cc, s2, true,  94, "#A6");
260
    checkMcf(mcf7, mcf7.run(),
263
    mcf.reset();
264
    checkMcf(mcf, mcf.run(),
261 265
             gr, l1, cu, cc, s3, true,   0, "#A7");
262
    checkMcf(mcf8, mcf8.boundMaps(l2, u).run(),
266
    checkMcf(mcf, mcf.boundMaps(l2, u).run(),
263 267
             gr, l2, u, cc, s3, false,   0, "#A8");
264 268
  }
265 269

	
266 270
  // B. Test NetworkSimplex with each pivot rule
267 271
  {
268
    NetworkSimplex<Digraph> mcf1(gr), mcf2(gr), mcf3(gr), mcf4(gr), mcf5(gr);
269
    NetworkSimplex<Digraph>::PivotRule pr;
272
    NetworkSimplex<Digraph> mcf(gr);
273
    mcf.supplyMap(s1).costMap(c).capacityMap(u).lowerMap(l2);
270 274

	
271
    pr = NetworkSimplex<Digraph>::FIRST_ELIGIBLE;
272
    checkMcf(mcf1, mcf1.boundMaps(l2, u).costMap(c).supplyMap(s1).run(pr),
275
    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::FIRST_ELIGIBLE),
273 276
             gr, l2, u, c, s1, true,  5970, "#B1");
274
    pr = NetworkSimplex<Digraph>::BEST_ELIGIBLE;
275
    checkMcf(mcf2, mcf2.boundMaps(l2, u).costMap(c).supplyMap(s1).run(pr),
277
    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::BEST_ELIGIBLE),
276 278
             gr, l2, u, c, s1, true,  5970, "#B2");
277
    pr = NetworkSimplex<Digraph>::BLOCK_SEARCH;
278
    checkMcf(mcf3, mcf3.boundMaps(l2, u).costMap(c).supplyMap(s1).run(pr),
279
    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::BLOCK_SEARCH),
279 280
             gr, l2, u, c, s1, true,  5970, "#B3");
280
    pr = NetworkSimplex<Digraph>::CANDIDATE_LIST;
281
    checkMcf(mcf4, mcf4.boundMaps(l2, u).costMap(c).supplyMap(s1).run(pr),
281
    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::CANDIDATE_LIST),
282 282
             gr, l2, u, c, s1, true,  5970, "#B4");
283
    pr = NetworkSimplex<Digraph>::ALTERING_LIST;
284
    checkMcf(mcf5, mcf5.boundMaps(l2, u).costMap(c).supplyMap(s1).run(pr),
283
    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::ALTERING_LIST),
285 284
             gr, l2, u, c, s1, true,  5970, "#B5");
286 285
  }
287 286

	
288 287
  return 0;
289 288
}
0 comments (0 inline)