# 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


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  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  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  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 nonzero lower bounds 
1084  1099  if (_plower) { 
… 
… 

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