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 6 line context
... ...
@@ -54,7 +54,8 @@
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.
... ...
@@ -1044,8 +1045,10 @@
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];
... ...
@@ -1069,7 +1072,7 @@
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)
... ...
@@ -1079,6 +1082,18 @@
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) {
... ...
@@ -1100,17 +1115,17 @@
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

	
... ...
@@ -1327,22 +1342,10 @@
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

	
0 comments (0 inline)