gravatar
kpeter (Peter Kovacs)
kpeter@inf.elte.hu
Support real types + numerical stability fix in NS (#254) - Real types are supported by appropriate inicialization. - A feature of the XTI spanning tree structure is removed to ensure numerical stability (could cause problems using integer types). The node potentials are updated always on the lower subtree, in order to prevent overflow problems. The former method isn't notably faster during to our tests.
0 1 0
default
1 file changed with 27 insertions and 24 deletions:
↑ Collapse diff ↑
Ignore white space 12 line context
... ...
@@ -51,13 +51,14 @@
51 51
  /// \tparam GR The digraph type the algorithm runs on.
52 52
  /// \tparam F The value type used for flow amounts, capacity bounds
53 53
  /// and supply values in the algorithm. By default it is \c int.
54 54
  /// \tparam C The value type used for costs and potentials in the
55 55
  /// algorithm. By default it is the same as \c F.
56 56
  ///
57
  /// \warning Both value types must be signed integer types.
57
  /// \warning Both value types must be signed and all input data must
58
  /// be integer.
58 59
  ///
59 60
  /// \note %NetworkSimplex provides five different pivot rule
60 61
  /// implementations. For more information see \ref PivotRule.
61 62
  template <typename GR, typename F = int, typename C = F>
62 63
  class NetworkSimplex
63 64
  {
... ...
@@ -1041,14 +1042,16 @@
1041 1042
      for (ArcIt e(_graph); e != INVALID; ++e) {
1042 1043
        _arc_ref[i] = e;
1043 1044
        if ((i += k) >= _arc_num) i = (i % k) + 1;
1044 1045
      }
1045 1046

	
1046 1047
      // Initialize arc maps
1047
      Flow max_cap = std::numeric_limits<Flow>::max();
1048
      Cost max_cost = std::numeric_limits<Cost>::max() / 4;
1048
      Flow inf_cap =
1049
        std::numeric_limits<Flow>::has_infinity ?
1050
        std::numeric_limits<Flow>::infinity() :
1051
        std::numeric_limits<Flow>::max();
1049 1052
      if (_pupper && _pcost) {
1050 1053
        for (int i = 0; i != _arc_num; ++i) {
1051 1054
          Arc e = _arc_ref[i];
1052 1055
          _source[i] = _node_id[_graph.source(e)];
1053 1056
          _target[i] = _node_id[_graph.target(e)];
1054 1057
          _cap[i] = (*_pupper)[e];
... ...
@@ -1066,22 +1069,34 @@
1066 1069
        }
1067 1070
        if (_pupper) {
1068 1071
          for (int i = 0; i != _arc_num; ++i)
1069 1072
            _cap[i] = (*_pupper)[_arc_ref[i]];
1070 1073
        } else {
1071 1074
          for (int i = 0; i != _arc_num; ++i)
1072
            _cap[i] = max_cap;
1075
            _cap[i] = inf_cap;
1073 1076
        }
1074 1077
        if (_pcost) {
1075 1078
          for (int i = 0; i != _arc_num; ++i)
1076 1079
            _cost[i] = (*_pcost)[_arc_ref[i]];
1077 1080
        } else {
1078 1081
          for (int i = 0; i != _arc_num; ++i)
1079 1082
            _cost[i] = 1;
1080 1083
        }
1081 1084
      }
1085
      
1086
      // Initialize artifical cost
1087
      Cost art_cost;
1088
      if (std::numeric_limits<Cost>::is_exact) {
1089
        art_cost = std::numeric_limits<Cost>::max() / 4 + 1;
1090
      } else {
1091
        art_cost = std::numeric_limits<Cost>::min();
1092
        for (int i = 0; i != _arc_num; ++i) {
1093
          if (_cost[i] > art_cost) art_cost = _cost[i];
1094
        }
1095
        art_cost = (art_cost + 1) * _node_num;
1096
      }
1082 1097

	
1083 1098
      // Remove non-zero lower bounds
1084 1099
      if (_plower) {
1085 1100
        for (int i = 0; i != _arc_num; ++i) {
1086 1101
          Flow c = (*_plower)[_arc_ref[i]];
1087 1102
          if (c != 0) {
... ...
@@ -1097,23 +1112,23 @@
1097 1112
        _thread[u] = u + 1;
1098 1113
        _rev_thread[u + 1] = u;
1099 1114
        _succ_num[u] = 1;
1100 1115
        _last_succ[u] = u;
1101 1116
        _parent[u] = _root;
1102 1117
        _pred[u] = e;
1103
        _cost[e] = max_cost;
1104
        _cap[e] = max_cap;
1118
        _cost[e] = art_cost;
1119
        _cap[e] = inf_cap;
1105 1120
        _state[e] = STATE_TREE;
1106 1121
        if (_supply[u] >= 0) {
1107 1122
          _flow[e] = _supply[u];
1108 1123
          _forward[u] = true;
1109
          _pi[u] = -max_cost;
1124
          _pi[u] = -art_cost;
1110 1125
        } else {
1111 1126
          _flow[e] = -_supply[u];
1112 1127
          _forward[u] = false;
1113
          _pi[u] = max_cost;
1128
          _pi[u] = art_cost;
1114 1129
        }
1115 1130
      }
1116 1131

	
1117 1132
      return true;
1118 1133
    }
1119 1134

	
... ...
@@ -1324,28 +1339,16 @@
1324 1339

	
1325 1340
    // Update potentials
1326 1341
    void updatePotential() {
1327 1342
      Cost sigma = _forward[u_in] ?
1328 1343
        _pi[v_in] - _pi[u_in] - _cost[_pred[u_in]] :
1329 1344
        _pi[v_in] - _pi[u_in] + _cost[_pred[u_in]];
1330
      if (_succ_num[u_in] > _node_num / 2) {
1331
        // Update in the upper subtree (which contains the root)
1332
        int before = _rev_thread[u_in];
1333
        int after = _thread[_last_succ[u_in]];
1334
        _thread[before] = after;
1335
        _pi[_root] -= sigma;
1336
        for (int u = _thread[_root]; u != _root; u = _thread[u]) {
1337
          _pi[u] -= sigma;
1338
        }
1339
        _thread[before] = u_in;
1340
      } else {
1341
        // Update in the lower subtree (which has been moved)
1342
        int end = _thread[_last_succ[u_in]];
1343
        for (int u = u_in; u != end; u = _thread[u]) {
1344
          _pi[u] += sigma;
1345
        }
1345
      // Update potentials in the subtree, which has been moved
1346
      int end = _thread[_last_succ[u_in]];
1347
      for (int u = u_in; u != end; u = _thread[u]) {
1348
        _pi[u] += sigma;
1346 1349
      }
1347 1350
    }
1348 1351

	
1349 1352
    // Execute the algorithm
1350 1353
    bool start(PivotRule pivot_rule) {
1351 1354
      // Select the pivot rule implementation
0 comments (0 inline)