Changeset - r655:6ac5d9ae1d3d
« 654:9ad8d2
656:e6927f »
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:
lemon/network_simplex.h
27
24
↑ Collapse diff ↑
lemon/network_simplex.h
Show white space 6 line context
... ... 
@@ -56,3 +56,4 @@
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 
  ///
... ... 
@@ -1046,4 +1047,6 @@
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) {
... ... 
@@ -1071,3 +1074,3 @@
1071 1074 
          for (int i = 0; i != _arc_num; ++i)
1072 
            _cap[i] = max_cap;
1075 
            _cap[i] = inf_cap;
1073 1076 
        }
... ... 
@@ -1081,2 +1084,14 @@
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 

			




	...
	...
	
		
@@ -1102,4 +1117,4 @@


			




	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;


			




	...
	...
	
		
@@ -1108,3 +1123,3 @@


			




	1108
	1123
	
		
          _forward[u] = true;


			




	1109
	
	
		
          _pi[u] = -max_cost;


			




	
	1124
	
		
          _pi[u] = -art_cost;


			




	1110
	1125
	
		
        } else {


			




	...
	...
	
		
@@ -1112,3 +1127,3 @@


			




	1112
	1127
	
		
          _forward[u] = false;


			




	1113
	
	
		
          _pi[u] = max_cost;


			




	
	1128
	
		
          _pi[u] = art_cost;


			




	1114
	1129
	
		
        }


			




	...
	...
	
		
@@ -1329,18 +1344,6 @@


			




	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
	
		
      }


			




        

    






    
    



    


    

     
    
0 comments (0 inline)