COIN-OR::LEMON - Graph Library

Changeset 910:f3bc4e9b5f3a in lemon


Ignore:
Timestamp:
02/20/10 18:39:03 (8 years ago)
Author:
Peter Kovacs <kpeter@…>
Branch:
default
Message:

New heuristics for MCF algorithms (#340)
and some implementation improvements.

  • A useful heuristic is added to NetworkSimplex? to make the initial pivots faster.
  • A powerful global update heuristic is added to CostScaling? and the implementation is reworked with various improvements.
  • Better relabeling in CostScaling? to improve numerical stability and make the code faster.
  • A small improvement is made in CapacityScaling? for better delta computation.
  • Add notes to the classes about the usage of vector<char> instead of vector<bool> for efficiency reasons.
Location:
lemon
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • lemon/capacity_scaling.h

    r887 r910  
    135135 
    136136    typedef std::vector<int> IntVector; 
    137     typedef std::vector<char> BoolVector; 
    138137    typedef std::vector<Value> ValueVector; 
    139138    typedef std::vector<Cost> CostVector; 
     139    typedef std::vector<char> BoolVector; 
     140    // Note: vector<char> is used instead of vector<bool> for efficiency reasons 
    140141 
    141142  private: 
     
    765766      if (_factor > 1) { 
    766767        // With scaling 
    767         Value max_sup = 0, max_dem = 0; 
    768         for (int i = 0; i != _node_num; ++i) { 
     768        Value max_sup = 0, max_dem = 0, max_cap = 0; 
     769        for (int i = 0; i != _root; ++i) { 
    769770          Value ex = _excess[i]; 
    770771          if ( ex > max_sup) max_sup =  ex; 
    771772          if (-ex > max_dem) max_dem = -ex; 
    772         } 
    773         Value max_cap = 0; 
    774         for (int j = 0; j != _res_arc_num; ++j) { 
    775           if (_res_cap[j] > max_cap) max_cap = _res_cap[j]; 
     773          int last_out = _first_out[i+1] - 1; 
     774          for (int j = _first_out[i]; j != last_out; ++j) { 
     775            if (_res_cap[j] > max_cap) max_cap = _res_cap[j]; 
     776          } 
    776777        } 
    777778        max_sup = std::min(std::min(max_sup, max_dem), max_cap); 
  • lemon/cost_scaling.h

    r887 r910  
    198198 
    199199    typedef std::vector<int> IntVector; 
    200     typedef std::vector<char> BoolVector; 
    201200    typedef std::vector<Value> ValueVector; 
    202201    typedef std::vector<Cost> CostVector; 
    203202    typedef std::vector<LargeCost> LargeCostVector; 
     203    typedef std::vector<char> BoolVector; 
     204    // Note: vector<char> is used instead of vector<bool> for efficiency reasons 
    204205 
    205206  private: 
     
    245246    bool _have_lower; 
    246247    Value _sum_supply; 
     248    int _sup_node_num; 
    247249 
    248250    // Data structures for storing the digraph 
     
    273275    int _alpha; 
    274276 
     277    IntVector _buckets; 
     278    IntVector _bucket_next; 
     279    IntVector _bucket_prev; 
     280    IntVector _rank; 
     281    int _max_rank; 
     282   
    275283    // Data for a StaticDigraph structure 
    276284    typedef std::pair<int, int> IntPair; 
     
    803811      } 
    804812 
     813      _sup_node_num = 0; 
     814      for (NodeIt n(_graph); n != INVALID; ++n) { 
     815        if (sup[n] > 0) ++_sup_node_num; 
     816      } 
     817 
    805818      // Find a feasible flow using Circulation 
    806819      Circulation<Digraph, ConstMap<Arc, Value>, ValueArcMap, ValueNodeMap> 
     
    837850        for (int a = _first_out[_root]; a != _res_arc_num; ++a) { 
    838851          int ra = _reverse[a]; 
    839           _res_cap[a] = 1; 
     852          _res_cap[a] = 0; 
    840853          _res_cap[ra] = 0; 
    841854          _cost[a] = 0; 
     
    851864      // Maximum path length for partial augment 
    852865      const int MAX_PATH_LENGTH = 4; 
    853        
     866 
     867      // Initialize data structures for buckets       
     868      _max_rank = _alpha * _res_node_num; 
     869      _buckets.resize(_max_rank); 
     870      _bucket_next.resize(_res_node_num + 1); 
     871      _bucket_prev.resize(_res_node_num + 1); 
     872      _rank.resize(_res_node_num + 1); 
     873   
    854874      // Execute the algorithm 
    855875      switch (method) { 
     
    890910      } 
    891911    } 
     912     
     913    // Initialize a cost scaling phase 
     914    void initPhase() { 
     915      // Saturate arcs not satisfying the optimality condition 
     916      for (int u = 0; u != _res_node_num; ++u) { 
     917        int last_out = _first_out[u+1]; 
     918        LargeCost pi_u = _pi[u]; 
     919        for (int a = _first_out[u]; a != last_out; ++a) { 
     920          int v = _target[a]; 
     921          if (_res_cap[a] > 0 && _cost[a] + pi_u - _pi[v] < 0) { 
     922            Value delta = _res_cap[a]; 
     923            _excess[u] -= delta; 
     924            _excess[v] += delta; 
     925            _res_cap[a] = 0; 
     926            _res_cap[_reverse[a]] += delta; 
     927          } 
     928        } 
     929      } 
     930       
     931      // Find active nodes (i.e. nodes with positive excess) 
     932      for (int u = 0; u != _res_node_num; ++u) { 
     933        if (_excess[u] > 0) _active_nodes.push_back(u); 
     934      } 
     935 
     936      // Initialize the next arcs 
     937      for (int u = 0; u != _res_node_num; ++u) { 
     938        _next_out[u] = _first_out[u]; 
     939      } 
     940    } 
     941     
     942    // Early termination heuristic 
     943    bool earlyTermination() { 
     944      const double EARLY_TERM_FACTOR = 3.0; 
     945 
     946      // Build a static residual graph 
     947      _arc_vec.clear(); 
     948      _cost_vec.clear(); 
     949      for (int j = 0; j != _res_arc_num; ++j) { 
     950        if (_res_cap[j] > 0) { 
     951          _arc_vec.push_back(IntPair(_source[j], _target[j])); 
     952          _cost_vec.push_back(_cost[j] + 1); 
     953        } 
     954      } 
     955      _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end()); 
     956 
     957      // Run Bellman-Ford algorithm to check if the current flow is optimal 
     958      BellmanFord<StaticDigraph, LargeCostArcMap> bf(_sgr, _cost_map); 
     959      bf.init(0); 
     960      bool done = false; 
     961      int K = int(EARLY_TERM_FACTOR * std::sqrt(double(_res_node_num))); 
     962      for (int i = 0; i < K && !done; ++i) { 
     963        done = bf.processNextWeakRound(); 
     964      } 
     965      return done; 
     966    } 
     967 
     968    // Global potential update heuristic 
     969    void globalUpdate() { 
     970      int bucket_end = _root + 1; 
     971     
     972      // Initialize buckets 
     973      for (int r = 0; r != _max_rank; ++r) { 
     974        _buckets[r] = bucket_end; 
     975      } 
     976      Value total_excess = 0; 
     977      for (int i = 0; i != _res_node_num; ++i) { 
     978        if (_excess[i] < 0) { 
     979          _rank[i] = 0; 
     980          _bucket_next[i] = _buckets[0]; 
     981          _bucket_prev[_buckets[0]] = i; 
     982          _buckets[0] = i; 
     983        } else { 
     984          total_excess += _excess[i]; 
     985          _rank[i] = _max_rank; 
     986        } 
     987      } 
     988      if (total_excess == 0) return; 
     989 
     990      // Search the buckets 
     991      int r = 0; 
     992      for ( ; r != _max_rank; ++r) { 
     993        while (_buckets[r] != bucket_end) { 
     994          // Remove the first node from the current bucket 
     995          int u = _buckets[r]; 
     996          _buckets[r] = _bucket_next[u]; 
     997           
     998          // Search the incomming arcs of u 
     999          LargeCost pi_u = _pi[u]; 
     1000          int last_out = _first_out[u+1]; 
     1001          for (int a = _first_out[u]; a != last_out; ++a) { 
     1002            int ra = _reverse[a]; 
     1003            if (_res_cap[ra] > 0) { 
     1004              int v = _source[ra]; 
     1005              int old_rank_v = _rank[v]; 
     1006              if (r < old_rank_v) { 
     1007                // Compute the new rank of v 
     1008                LargeCost nrc = (_cost[ra] + _pi[v] - pi_u) / _epsilon; 
     1009                int new_rank_v = old_rank_v; 
     1010                if (nrc < LargeCost(_max_rank)) 
     1011                  new_rank_v = r + 1 + int(nrc); 
     1012                   
     1013                // Change the rank of v 
     1014                if (new_rank_v < old_rank_v) { 
     1015                  _rank[v] = new_rank_v; 
     1016                  _next_out[v] = _first_out[v]; 
     1017                   
     1018                  // Remove v from its old bucket 
     1019                  if (old_rank_v < _max_rank) { 
     1020                    if (_buckets[old_rank_v] == v) { 
     1021                      _buckets[old_rank_v] = _bucket_next[v]; 
     1022                    } else { 
     1023                      _bucket_next[_bucket_prev[v]] = _bucket_next[v]; 
     1024                      _bucket_prev[_bucket_next[v]] = _bucket_prev[v]; 
     1025                    } 
     1026                  } 
     1027                   
     1028                  // Insert v to its new bucket 
     1029                  _bucket_next[v] = _buckets[new_rank_v]; 
     1030                  _bucket_prev[_buckets[new_rank_v]] = v; 
     1031                  _buckets[new_rank_v] = v; 
     1032                } 
     1033              } 
     1034            } 
     1035          } 
     1036 
     1037          // Finish search if there are no more active nodes 
     1038          if (_excess[u] > 0) { 
     1039            total_excess -= _excess[u]; 
     1040            if (total_excess <= 0) break; 
     1041          } 
     1042        } 
     1043        if (total_excess <= 0) break; 
     1044      } 
     1045       
     1046      // Relabel nodes 
     1047      for (int u = 0; u != _res_node_num; ++u) { 
     1048        int k = std::min(_rank[u], r); 
     1049        if (k > 0) { 
     1050          _pi[u] -= _epsilon * k; 
     1051          _next_out[u] = _first_out[u]; 
     1052        } 
     1053      } 
     1054    } 
    8921055 
    8931056    /// Execute the algorithm performing augment and relabel operations 
    8941057    void startAugment(int max_length = std::numeric_limits<int>::max()) { 
    8951058      // Paramters for heuristics 
    896       const int BF_HEURISTIC_EPSILON_BOUND = 1000; 
    897       const int BF_HEURISTIC_BOUND_FACTOR  = 3; 
    898  
     1059      const int EARLY_TERM_EPSILON_LIMIT = 1000; 
     1060      const double GLOBAL_UPDATE_FACTOR = 3.0; 
     1061 
     1062      const int global_update_freq = int(GLOBAL_UPDATE_FACTOR * 
     1063        (_res_node_num + _sup_node_num * _sup_node_num)); 
     1064      int next_update_limit = global_update_freq; 
     1065       
     1066      int relabel_cnt = 0; 
     1067       
    8991068      // Perform cost scaling phases 
    900       IntVector pred_arc(_res_node_num); 
    901       std::vector<int> path_nodes; 
     1069      std::vector<int> path; 
    9021070      for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ? 
    9031071                                        1 : _epsilon / _alpha ) 
    9041072      { 
    905         // "Early Termination" heuristic: use Bellman-Ford algorithm 
    906         // to check if the current flow is optimal 
    907         if (_epsilon <= BF_HEURISTIC_EPSILON_BOUND) { 
    908           _arc_vec.clear(); 
    909           _cost_vec.clear(); 
    910           for (int j = 0; j != _res_arc_num; ++j) { 
    911             if (_res_cap[j] > 0) { 
    912               _arc_vec.push_back(IntPair(_source[j], _target[j])); 
    913               _cost_vec.push_back(_cost[j] + 1); 
    914             } 
    915           } 
    916           _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end()); 
    917  
    918           BellmanFord<StaticDigraph, LargeCostArcMap> bf(_sgr, _cost_map); 
    919           bf.init(0); 
    920           bool done = false; 
    921           int K = int(BF_HEURISTIC_BOUND_FACTOR * sqrt(_res_node_num)); 
    922           for (int i = 0; i < K && !done; ++i) 
    923             done = bf.processNextWeakRound(); 
    924           if (done) break; 
    925         } 
    926  
    927         // Saturate arcs not satisfying the optimality condition 
    928         for (int a = 0; a != _res_arc_num; ++a) { 
    929           if (_res_cap[a] > 0 && 
    930               _cost[a] + _pi[_source[a]] - _pi[_target[a]] < 0) { 
    931             Value delta = _res_cap[a]; 
    932             _excess[_source[a]] -= delta; 
    933             _excess[_target[a]] += delta; 
    934             _res_cap[a] = 0; 
    935             _res_cap[_reverse[a]] += delta; 
    936           } 
     1073        // Early termination heuristic 
     1074        if (_epsilon <= EARLY_TERM_EPSILON_LIMIT) { 
     1075          if (earlyTermination()) break; 
    9371076        } 
    9381077         
    939         // Find active nodes (i.e. nodes with positive excess) 
    940         for (int u = 0; u != _res_node_num; ++u) { 
    941           if (_excess[u] > 0) _active_nodes.push_back(u); 
    942         } 
    943  
    944         // Initialize the next arcs 
    945         for (int u = 0; u != _res_node_num; ++u) { 
    946           _next_out[u] = _first_out[u]; 
    947         } 
    948  
     1078        // Initialize current phase 
     1079        initPhase(); 
     1080         
    9491081        // Perform partial augment and relabel operations 
    9501082        while (true) { 
     
    9561088          if (_active_nodes.size() == 0) break; 
    9571089          int start = _active_nodes.front(); 
    958           path_nodes.clear(); 
    959           path_nodes.push_back(start); 
    9601090 
    9611091          // Find an augmenting path from the start node 
     1092          path.clear(); 
    9621093          int tip = start; 
    963           while (_excess[tip] >= 0 && 
    964                  int(path_nodes.size()) <= max_length) { 
     1094          while (_excess[tip] >= 0 && int(path.size()) < max_length) { 
    9651095            int u; 
    966             LargeCost min_red_cost, rc; 
    967             int last_out = _sum_supply < 0 ? 
    968               _first_out[tip+1] : _first_out[tip+1] - 1; 
     1096            LargeCost min_red_cost, rc, pi_tip = _pi[tip]; 
     1097            int last_out = _first_out[tip+1]; 
    9691098            for (int a = _next_out[tip]; a != last_out; ++a) { 
    970               if (_res_cap[a] > 0 && 
    971                   _cost[a] + _pi[_source[a]] - _pi[_target[a]] < 0) { 
    972                 u = _target[a]; 
    973                 pred_arc[u] = a; 
     1099              u = _target[a]; 
     1100              if (_res_cap[a] > 0 && _cost[a] + pi_tip - _pi[u] < 0) { 
     1101                path.push_back(a); 
    9741102                _next_out[tip] = a; 
    9751103                tip = u; 
    976                 path_nodes.push_back(tip); 
    9771104                goto next_step; 
    9781105              } 
     
    9801107 
    9811108            // Relabel tip node 
    982             min_red_cost = std::numeric_limits<LargeCost>::max() / 2; 
     1109            min_red_cost = std::numeric_limits<LargeCost>::max(); 
     1110            if (tip != start) { 
     1111              int ra = _reverse[path.back()]; 
     1112              min_red_cost = _cost[ra] + pi_tip - _pi[_target[ra]]; 
     1113            } 
    9831114            for (int a = _first_out[tip]; a != last_out; ++a) { 
    984               rc = _cost[a] + _pi[_source[a]] - _pi[_target[a]]; 
     1115              rc = _cost[a] + pi_tip - _pi[_target[a]]; 
    9851116              if (_res_cap[a] > 0 && rc < min_red_cost) { 
    9861117                min_red_cost = rc; 
     
    9881119            } 
    9891120            _pi[tip] -= min_red_cost + _epsilon; 
    990  
    991             // Reset the next arc of tip 
    9921121            _next_out[tip] = _first_out[tip]; 
     1122            ++relabel_cnt; 
    9931123 
    9941124            // Step back 
    9951125            if (tip != start) { 
    996               path_nodes.pop_back(); 
    997               tip = path_nodes.back(); 
     1126              tip = _source[path.back()]; 
     1127              path.pop_back(); 
    9981128            } 
    9991129 
     
    10031133          // Augment along the found path (as much flow as possible) 
    10041134          Value delta; 
    1005           int u, v = path_nodes.front(), pa; 
    1006           for (int i = 1; i < int(path_nodes.size()); ++i) { 
     1135          int pa, u, v = start; 
     1136          for (int i = 0; i != int(path.size()); ++i) { 
     1137            pa = path[i]; 
    10071138            u = v; 
    1008             v = path_nodes[i]; 
    1009             pa = pred_arc[v]; 
     1139            v = _target[pa]; 
    10101140            delta = std::min(_res_cap[pa], _excess[u]); 
    10111141            _res_cap[pa] -= delta; 
     
    10161146              _active_nodes.push_back(v); 
    10171147          } 
     1148 
     1149          // Global update heuristic 
     1150          if (relabel_cnt >= next_update_limit) { 
     1151            globalUpdate(); 
     1152            next_update_limit += global_update_freq; 
     1153          } 
    10181154        } 
    10191155      } 
     
    10231159    void startPush() { 
    10241160      // Paramters for heuristics 
    1025       const int BF_HEURISTIC_EPSILON_BOUND = 1000; 
    1026       const int BF_HEURISTIC_BOUND_FACTOR  = 3; 
    1027  
     1161      const int EARLY_TERM_EPSILON_LIMIT = 1000; 
     1162      const double GLOBAL_UPDATE_FACTOR = 2.0; 
     1163 
     1164      const int global_update_freq = int(GLOBAL_UPDATE_FACTOR * 
     1165        (_res_node_num + _sup_node_num * _sup_node_num)); 
     1166      int next_update_limit = global_update_freq; 
     1167 
     1168      int relabel_cnt = 0; 
     1169       
    10281170      // Perform cost scaling phases 
    10291171      BoolVector hyper(_res_node_num, false); 
     1172      LargeCostVector hyper_cost(_res_node_num); 
    10301173      for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ? 
    10311174                                        1 : _epsilon / _alpha ) 
    10321175      { 
    1033         // "Early Termination" heuristic: use Bellman-Ford algorithm 
    1034         // to check if the current flow is optimal 
    1035         if (_epsilon <= BF_HEURISTIC_EPSILON_BOUND) { 
    1036           _arc_vec.clear(); 
    1037           _cost_vec.clear(); 
    1038           for (int j = 0; j != _res_arc_num; ++j) { 
    1039             if (_res_cap[j] > 0) { 
    1040               _arc_vec.push_back(IntPair(_source[j], _target[j])); 
    1041               _cost_vec.push_back(_cost[j] + 1); 
    1042             } 
    1043           } 
    1044           _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end()); 
    1045  
    1046           BellmanFord<StaticDigraph, LargeCostArcMap> bf(_sgr, _cost_map); 
    1047           bf.init(0); 
    1048           bool done = false; 
    1049           int K = int(BF_HEURISTIC_BOUND_FACTOR * sqrt(_res_node_num)); 
    1050           for (int i = 0; i < K && !done; ++i) 
    1051             done = bf.processNextWeakRound(); 
    1052           if (done) break; 
    1053         } 
    1054  
    1055         // Saturate arcs not satisfying the optimality condition 
    1056         for (int a = 0; a != _res_arc_num; ++a) { 
    1057           if (_res_cap[a] > 0 && 
    1058               _cost[a] + _pi[_source[a]] - _pi[_target[a]] < 0) { 
    1059             Value delta = _res_cap[a]; 
    1060             _excess[_source[a]] -= delta; 
    1061             _excess[_target[a]] += delta; 
    1062             _res_cap[a] = 0; 
    1063             _res_cap[_reverse[a]] += delta; 
    1064           } 
    1065         } 
    1066  
    1067         // Find active nodes (i.e. nodes with positive excess) 
    1068         for (int u = 0; u != _res_node_num; ++u) { 
    1069           if (_excess[u] > 0) _active_nodes.push_back(u); 
    1070         } 
    1071  
    1072         // Initialize the next arcs 
    1073         for (int u = 0; u != _res_node_num; ++u) { 
    1074           _next_out[u] = _first_out[u]; 
    1075         } 
     1176        // Early termination heuristic 
     1177        if (_epsilon <= EARLY_TERM_EPSILON_LIMIT) { 
     1178          if (earlyTermination()) break; 
     1179        } 
     1180         
     1181        // Initialize current phase 
     1182        initPhase(); 
    10761183 
    10771184        // Perform push and relabel operations 
    10781185        while (_active_nodes.size() > 0) { 
    1079           LargeCost min_red_cost, rc; 
     1186          LargeCost min_red_cost, rc, pi_n; 
    10801187          Value delta; 
    10811188          int n, t, a, last_out = _res_arc_num; 
    10821189 
     1190        next_node: 
    10831191          // Select an active node (FIFO selection) 
    1084         next_node: 
    10851192          n = _active_nodes.front(); 
    1086           last_out = _sum_supply < 0 ? 
    1087             _first_out[n+1] : _first_out[n+1] - 1; 
    1088  
     1193          last_out = _first_out[n+1]; 
     1194          pi_n = _pi[n]; 
     1195           
    10891196          // Perform push operations if there are admissible arcs 
    10901197          if (_excess[n] > 0) { 
    10911198            for (a = _next_out[n]; a != last_out; ++a) { 
    10921199              if (_res_cap[a] > 0 && 
    1093                   _cost[a] + _pi[_source[a]] - _pi[_target[a]] < 0) { 
     1200                  _cost[a] + pi_n - _pi[_target[a]] < 0) { 
    10941201                delta = std::min(_res_cap[a], _excess[n]); 
    10951202                t = _target[a]; 
     
    10971204                // Push-look-ahead heuristic 
    10981205                Value ahead = -_excess[t]; 
    1099                 int last_out_t = _sum_supply < 0 ? 
    1100                   _first_out[t+1] : _first_out[t+1] - 1; 
     1206                int last_out_t = _first_out[t+1]; 
     1207                LargeCost pi_t = _pi[t]; 
    11011208                for (int ta = _next_out[t]; ta != last_out_t; ++ta) { 
    11021209                  if (_res_cap[ta] > 0 &&  
    1103                       _cost[ta] + _pi[_source[ta]] - _pi[_target[ta]] < 0) 
     1210                      _cost[ta] + pi_t - _pi[_target[ta]] < 0) 
    11041211                    ahead += _res_cap[ta]; 
    11051212                  if (ahead >= delta) break; 
     
    11081215 
    11091216                // Push flow along the arc 
    1110                 if (ahead < delta) { 
     1217                if (ahead < delta && !hyper[t]) { 
    11111218                  _res_cap[a] -= ahead; 
    11121219                  _res_cap[_reverse[a]] += ahead; 
     
    11151222                  _active_nodes.push_front(t); 
    11161223                  hyper[t] = true; 
     1224                  hyper_cost[t] = _cost[a] + pi_n - pi_t; 
    11171225                  _next_out[n] = a; 
    11181226                  goto next_node; 
     
    11371245          // Relabel the node if it is still active (or hyper) 
    11381246          if (_excess[n] > 0 || hyper[n]) { 
    1139             min_red_cost = std::numeric_limits<LargeCost>::max() / 2; 
     1247             min_red_cost = hyper[n] ? -hyper_cost[n] : 
     1248               std::numeric_limits<LargeCost>::max(); 
    11401249            for (int a = _first_out[n]; a != last_out; ++a) { 
    1141               rc = _cost[a] + _pi[_source[a]] - _pi[_target[a]]; 
     1250              rc = _cost[a] + pi_n - _pi[_target[a]]; 
    11421251              if (_res_cap[a] > 0 && rc < min_red_cost) { 
    11431252                min_red_cost = rc; 
     
    11451254            } 
    11461255            _pi[n] -= min_red_cost + _epsilon; 
     1256            _next_out[n] = _first_out[n]; 
    11471257            hyper[n] = false; 
    1148  
    1149             // Reset the next arc 
    1150             _next_out[n] = _first_out[n]; 
     1258            ++relabel_cnt; 
    11511259          } 
    11521260         
     
    11581266            _active_nodes.pop_front(); 
    11591267          } 
     1268           
     1269          // Global update heuristic 
     1270          if (relabel_cnt >= next_update_limit) { 
     1271            globalUpdate(); 
     1272            for (int u = 0; u != _res_node_num; ++u) 
     1273              hyper[u] = false; 
     1274            next_update_limit += global_update_freq; 
     1275          } 
    11601276        } 
    11611277      } 
  • lemon/cycle_canceling.h

    r886 r910  
    145145     
    146146    typedef std::vector<int> IntVector; 
    147     typedef std::vector<char> CharVector; 
    148147    typedef std::vector<double> DoubleVector; 
    149148    typedef std::vector<Value> ValueVector; 
    150149    typedef std::vector<Cost> CostVector; 
     150    typedef std::vector<char> BoolVector; 
     151    // Note: vector<char> is used instead of vector<bool> for efficiency reasons 
    151152 
    152153  private: 
     
    199200    IntArcMap _arc_idb; 
    200201    IntVector _first_out; 
    201     CharVector _forward; 
     202    BoolVector _forward; 
    202203    IntVector _source; 
    203204    IntVector _target; 
     
    934935      DoubleVector pi(_res_node_num, 0.0); 
    935936      IntVector level(_res_node_num); 
    936       CharVector reached(_res_node_num); 
    937       CharVector processed(_res_node_num); 
     937      BoolVector reached(_res_node_num); 
     938      BoolVector processed(_res_node_num); 
    938939      IntVector pred_node(_res_node_num); 
    939940      IntVector pred_arc(_res_node_num); 
  • lemon/network_simplex.h

    r878 r910  
    165165 
    166166    typedef std::vector<int> IntVector; 
    167     typedef std::vector<char> CharVector; 
    168167    typedef std::vector<Value> ValueVector; 
    169168    typedef std::vector<Cost> CostVector; 
     169    typedef std::vector<char> BoolVector; 
     170    // Note: vector<char> is used instead of vector<bool> for efficiency reasons 
    170171 
    171172    // State constants for arcs 
     
    213214    IntVector _last_succ; 
    214215    IntVector _dirty_revs; 
    215     CharVector _forward; 
    216     CharVector _state; 
     216    BoolVector _forward; 
     217    BoolVector _state; 
    217218    int _root; 
    218219 
     
    245246      const IntVector  &_target; 
    246247      const CostVector &_cost; 
    247       const CharVector &_state; 
     248      const BoolVector &_state; 
    248249      const CostVector &_pi; 
    249250      int &_in_arc; 
     
    266267      bool findEnteringArc() { 
    267268        Cost c; 
    268         for (int e = _next_arc; e < _search_arc_num; ++e) { 
     269        for (int e = _next_arc; e != _search_arc_num; ++e) { 
    269270          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); 
    270271          if (c < 0) { 
     
    274275          } 
    275276        } 
    276         for (int e = 0; e < _next_arc; ++e) { 
     277        for (int e = 0; e != _next_arc; ++e) { 
    277278          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); 
    278279          if (c < 0) { 
     
    297298      const IntVector  &_target; 
    298299      const CostVector &_cost; 
    299       const CharVector &_state; 
     300      const BoolVector &_state; 
    300301      const CostVector &_pi; 
    301302      int &_in_arc; 
     
    314315      bool findEnteringArc() { 
    315316        Cost c, min = 0; 
    316         for (int e = 0; e < _search_arc_num; ++e) { 
     317        for (int e = 0; e != _search_arc_num; ++e) { 
    317318          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); 
    318319          if (c < min) { 
     
    336337      const IntVector  &_target; 
    337338      const CostVector &_cost; 
    338       const CharVector &_state; 
     339      const BoolVector &_state; 
    339340      const CostVector &_pi; 
    340341      int &_in_arc; 
     
    355356      { 
    356357        // The main parameters of the pivot rule 
    357         const double BLOCK_SIZE_FACTOR = 0.5; 
     358        const double BLOCK_SIZE_FACTOR = 1.0; 
    358359        const int MIN_BLOCK_SIZE = 10; 
    359360 
     
    368369        int cnt = _block_size; 
    369370        int e; 
    370         for (e = _next_arc; e < _search_arc_num; ++e) { 
     371        for (e = _next_arc; e != _search_arc_num; ++e) { 
    371372          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); 
    372373          if (c < min) { 
     
    379380          } 
    380381        } 
    381         for (e = 0; e < _next_arc; ++e) { 
     382        for (e = 0; e != _next_arc; ++e) { 
    382383          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); 
    383384          if (c < min) { 
     
    409410      const IntVector  &_target; 
    410411      const CostVector &_cost; 
    411       const CharVector &_state; 
     412      const BoolVector &_state; 
    412413      const CostVector &_pi; 
    413414      int &_in_arc; 
     
    470471        min = 0; 
    471472        _curr_length = 0; 
    472         for (e = _next_arc; e < _search_arc_num; ++e) { 
     473        for (e = _next_arc; e != _search_arc_num; ++e) { 
    473474          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); 
    474475          if (c < 0) { 
     
    481482          } 
    482483        } 
    483         for (e = 0; e < _next_arc; ++e) { 
     484        for (e = 0; e != _next_arc; ++e) { 
    484485          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); 
    485486          if (c < 0) { 
     
    512513      const IntVector  &_target; 
    513514      const CostVector &_cost; 
    514       const CharVector &_state; 
     515      const BoolVector &_state; 
    515516      const CostVector &_pi; 
    516517      int &_in_arc; 
     
    565566        // Check the current candidate list 
    566567        int e; 
    567         for (int i = 0; i < _curr_length; ++i) { 
     568        for (int i = 0; i != _curr_length; ++i) { 
    568569          e = _candidates[i]; 
    569570          _cand_cost[e] = _state[e] * 
     
    578579        int limit = _head_length; 
    579580 
    580         for (e = _next_arc; e < _search_arc_num; ++e) { 
     581        for (e = _next_arc; e != _search_arc_num; ++e) { 
    581582          _cand_cost[e] = _state[e] * 
    582583            (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); 
     
    590591          } 
    591592        } 
    592         for (e = 0; e < _next_arc; ++e) { 
     593        for (e = 0; e != _next_arc; ++e) { 
    593594          _cand_cost[e] = _state[e] * 
    594595            (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); 
     
    13291330 
    13301331      // Update _rev_thread using the new _thread values 
    1331       for (int i = 0; i < int(_dirty_revs.size()); ++i) { 
     1332      for (int i = 0; i != int(_dirty_revs.size()); ++i) { 
    13321333        u = _dirty_revs[i]; 
    13331334        _rev_thread[_thread[u]] = u; 
     
    13991400        _pi[u] += sigma; 
    14001401      } 
     1402    } 
     1403 
     1404    // Heuristic initial pivots 
     1405    bool initialPivots() { 
     1406      Value curr, total = 0; 
     1407      std::vector<Node> supply_nodes, demand_nodes; 
     1408      for (NodeIt u(_graph); u != INVALID; ++u) { 
     1409        curr = _supply[_node_id[u]]; 
     1410        if (curr > 0) { 
     1411          total += curr; 
     1412          supply_nodes.push_back(u); 
     1413        } 
     1414        else if (curr < 0) { 
     1415          demand_nodes.push_back(u); 
     1416        } 
     1417      } 
     1418      if (_sum_supply > 0) total -= _sum_supply; 
     1419      if (total <= 0) return true; 
     1420 
     1421      IntVector arc_vector; 
     1422      if (_sum_supply >= 0) { 
     1423        if (supply_nodes.size() == 1 && demand_nodes.size() == 1) { 
     1424          // Perform a reverse graph search from the sink to the source 
     1425          typename GR::template NodeMap<bool> reached(_graph, false); 
     1426          Node s = supply_nodes[0], t = demand_nodes[0]; 
     1427          std::vector<Node> stack; 
     1428          reached[t] = true; 
     1429          stack.push_back(t); 
     1430          while (!stack.empty()) { 
     1431            Node u, v = stack.back(); 
     1432            stack.pop_back(); 
     1433            if (v == s) break; 
     1434            for (InArcIt a(_graph, v); a != INVALID; ++a) { 
     1435              if (reached[u = _graph.source(a)]) continue; 
     1436              int j = _arc_id[a]; 
     1437              if (_cap[j] >= total) { 
     1438                arc_vector.push_back(j); 
     1439                reached[u] = true; 
     1440                stack.push_back(u); 
     1441              } 
     1442            } 
     1443          } 
     1444        } else { 
     1445          // Find the min. cost incomming arc for each demand node 
     1446          for (int i = 0; i != int(demand_nodes.size()); ++i) { 
     1447            Node v = demand_nodes[i]; 
     1448            Cost c, min_cost = std::numeric_limits<Cost>::max(); 
     1449            Arc min_arc = INVALID; 
     1450            for (InArcIt a(_graph, v); a != INVALID; ++a) { 
     1451              c = _cost[_arc_id[a]]; 
     1452              if (c < min_cost) { 
     1453                min_cost = c; 
     1454                min_arc = a; 
     1455              } 
     1456            } 
     1457            if (min_arc != INVALID) { 
     1458              arc_vector.push_back(_arc_id[min_arc]); 
     1459            } 
     1460          } 
     1461        } 
     1462      } else { 
     1463        // Find the min. cost outgoing arc for each supply node 
     1464        for (int i = 0; i != int(supply_nodes.size()); ++i) { 
     1465          Node u = supply_nodes[i]; 
     1466          Cost c, min_cost = std::numeric_limits<Cost>::max(); 
     1467          Arc min_arc = INVALID; 
     1468          for (OutArcIt a(_graph, u); a != INVALID; ++a) { 
     1469            c = _cost[_arc_id[a]]; 
     1470            if (c < min_cost) { 
     1471              min_cost = c; 
     1472              min_arc = a; 
     1473            } 
     1474          } 
     1475          if (min_arc != INVALID) { 
     1476            arc_vector.push_back(_arc_id[min_arc]); 
     1477          } 
     1478        } 
     1479      } 
     1480 
     1481      // Perform heuristic initial pivots 
     1482      for (int i = 0; i != int(arc_vector.size()); ++i) { 
     1483        in_arc = arc_vector[i]; 
     1484        if (_state[in_arc] * (_cost[in_arc] + _pi[_source[in_arc]] - 
     1485            _pi[_target[in_arc]]) >= 0) continue; 
     1486        findJoinNode(); 
     1487        bool change = findLeavingArc(); 
     1488        if (delta >= MAX) return false; 
     1489        changeFlow(change); 
     1490        if (change) { 
     1491          updateTreeStructure(); 
     1492          updatePotential(); 
     1493        } 
     1494      } 
     1495      return true; 
    14011496    } 
    14021497 
     
    14231518      PivotRuleImpl pivot(*this); 
    14241519 
     1520      // Perform heuristic initial pivots 
     1521      if (!initialPivots()) return UNBOUNDED; 
     1522 
    14251523      // Execute the Network Simplex algorithm 
    14261524      while (pivot.findEnteringArc()) { 
Note: See TracChangeset for help on using the changeset viewer.