COIN-OR::LEMON - Graph Library

Ticket #254: ns-numerical-6ac5d9ae1d3d.patch

File ns-numerical-6ac5d9ae1d3d.patch, 3.9 KB (added by Peter Kovacs, 10 years ago)
  • lemon/network_simplex.h

    # HG changeset patch
    # User Peter Kovacs <kpeter@inf.elte.hu>
    # Date 1238777955 -7200
    # Node ID 6ac5d9ae1d3dfdfcc48802e03bc2160d5881a78c
    # Parent  9ad8d2122b50707d35ae3709c49e0c9368a2a6c4
    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.
    
    diff --git a/lemon/network_simplex.h b/lemon/network_simplex.h
    a b  
    5454  /// \tparam C The value type used for costs and potentials in the
    5555  /// algorithm. By default it is the same as \c F.
    5656  ///
    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.
    5859  ///
    5960  /// \note %NetworkSimplex provides five different pivot rule
    6061  /// implementations. For more information see \ref PivotRule.
     
    10441045      }
    10451046
    10461047      // 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();
    10491052      if (_pupper && _pcost) {
    10501053        for (int i = 0; i != _arc_num; ++i) {
    10511054          Arc e = _arc_ref[i];
     
    10691072            _cap[i] = (*_pupper)[_arc_ref[i]];
    10701073        } else {
    10711074          for (int i = 0; i != _arc_num; ++i)
    1072             _cap[i] = max_cap;
     1075            _cap[i] = inf_cap;
    10731076        }
    10741077        if (_pcost) {
    10751078          for (int i = 0; i != _arc_num; ++i)
     
    10791082            _cost[i] = 1;
    10801083        }
    10811084      }
     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      }
    10821097
    10831098      // Remove non-zero lower bounds
    10841099      if (_plower) {
     
    11001115        _last_succ[u] = u;
    11011116        _parent[u] = _root;
    11021117        _pred[u] = e;
    1103         _cost[e] = max_cost;
    1104         _cap[e] = max_cap;
     1118        _cost[e] = art_cost;
     1119        _cap[e] = inf_cap;
    11051120        _state[e] = STATE_TREE;
    11061121        if (_supply[u] >= 0) {
    11071122          _flow[e] = _supply[u];
    11081123          _forward[u] = true;
    1109           _pi[u] = -max_cost;
     1124          _pi[u] = -art_cost;
    11101125        } else {
    11111126          _flow[e] = -_supply[u];
    11121127          _forward[u] = false;
    1113           _pi[u] = max_cost;
     1128          _pi[u] = art_cost;
    11141129        }
    11151130      }
    11161131
     
    13271342      Cost sigma = _forward[u_in] ?
    13281343        _pi[v_in] - _pi[u_in] - _cost[_pred[u_in]] :
    13291344        _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;
    13461349      }
    13471350    }
    13481351