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
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)