COIN-OR::LEMON - Graph Library

Changes in / [846:9d380bf27194:847:754272f20318] in lemon-main


Ignore:
Files:
13 edited

Legend:

Unmodified
Added
Removed
  • demo/arg_parser_demo.cc

    r440 r842  
    6666    .other("...");
    6767
     68  // Throw an exception when problems occurs. The default behavior is to
     69  // exit(1) on these cases, but this makes Valgrind falsely warn
     70  // about memory leaks.
     71  ap.throwOnProblems();
     72 
    6873  // Perform the parsing process
    6974  // (in case of any error it terminates the program)
    70   ap.parse();
     75  // The try {} construct is necessary only if the ap.trowOnProblems()
     76  // setting is in use.
     77  try {
     78    ap.parse();
     79  } catch (ArgParserException &) { return 1; }
    7180
    7281  // Check each option if it has been given and print its value
  • lemon/arg_parser.cc

    r440 r842  
    2121namespace lemon {
    2222
     23  void ArgParser::_terminate(ArgParserException::Reason reason) const
     24  {
     25    if(_exit_on_problems)
     26      exit(1);
     27    else throw(ArgParserException(reason));
     28  }
     29 
     30 
    2331  void ArgParser::_showHelp(void *p)
    2432  {
    2533    (static_cast<ArgParser*>(p))->showHelp();
    26     exit(1);
     34    (static_cast<ArgParser*>(p))->_terminate(ArgParserException::HELP);
    2735  }
    2836
    2937  ArgParser::ArgParser(int argc, const char * const *argv)
    30     :_argc(argc), _argv(argv), _command_name(argv[0]) {
     38    :_argc(argc), _argv(argv), _command_name(argv[0]),
     39    _exit_on_problems(true) {
    3140    funcOption("-help","Print a short help message",_showHelp,this);
    3241    synonym("help","-help");
     
    343352        i!=_others_help.end();++i) showHelp(i);
    344353    for(Opts::const_iterator i=_opts.begin();i!=_opts.end();++i) showHelp(i);
    345     exit(1);
     354    _terminate(ArgParserException::HELP);
    346355  }
    347356
     
    352361    std::cerr << "\nType '" << _command_name <<
    353362      " --help' to obtain a short summary on the usage.\n\n";
    354     exit(1);
     363    _terminate(ArgParserException::UNKNOWN_OPT);
    355364  }
    356365
     
    415424      std::cerr << "\nType '" << _command_name <<
    416425        " --help' to obtain a short summary on the usage.\n\n";
    417       exit(1);
     426      _terminate(ArgParserException::INVALID_OPT);
    418427    }
    419428  }
  • lemon/arg_parser.h

    r440 r842  
    3535namespace lemon {
    3636
     37  ///Exception used by ArgParser
     38  class ArgParserException : public Exception {
     39  public:
     40    enum Reason {
     41      HELP,         /// <tt>--help</tt> option was given
     42      UNKNOWN_OPT,  /// Unknown option was given
     43      INVALID_OPT   /// Invalid combination of options
     44    };
     45   
     46  private:
     47    Reason _reason;
     48   
     49  public:
     50    ///Constructor
     51    ArgParserException(Reason r) throw() : _reason(r) {}
     52    ///Virtual destructor
     53    virtual ~ArgParserException() throw() {}
     54    ///A short description of the exception
     55    virtual const char* what() const throw() {
     56      switch(_reason)
     57        {
     58        case HELP:
     59          return "lemon::ArgParseException: ask for help";
     60          break;
     61        case UNKNOWN_OPT:
     62          return "lemon::ArgParseException: unknown option";
     63          break;
     64        case INVALID_OPT:
     65          return "lemon::ArgParseException: invalid combination of options";
     66          break;
     67        }
     68      return "";
     69    }
     70    ///Return the reason for the failure
     71    Reason reason() const {return _reason; }
     72  };
     73
     74
    3775  ///Command line arguments parser
    3876
     
    104142    std::string _command_name;
    105143
    106 
     144   
    107145  private:
    108146    //Bind a function to an option.
     
    116154                    const std::string &help,
    117155                    void (*func)(void *),void *data);
     156
     157    bool _exit_on_problems;
     158   
     159    void _terminate(ArgParserException::Reason reason) const;
    118160
    119161  public:
     
    381423    const std::vector<std::string> &files() const { return _file_args; }
    382424
     425    ///Throw instead of exit in case of problems
     426    void throwOnProblems()
     427    {
     428      _exit_on_problems=false;
     429    }
    383430  };
    384431}
  • lemon/bellman_ford.h

    r825 r844  
    2929#include <lemon/error.h>
    3030#include <lemon/maps.h>
     31#include <lemon/tolerance.h>
    3132#include <lemon/path.h>
    3233
     
    3536namespace lemon {
    3637
    37   /// \brief Default OperationTraits for the BellmanFord algorithm class.
     38  /// \brief Default operation traits for the BellmanFord algorithm class.
    3839  /// 
    3940  /// This operation traits class defines all computational operations
     
    4243  /// If the numeric type does not have infinity value, then the maximum
    4344  /// value is used as extremal infinity value.
     45  ///
     46  /// \see BellmanFordToleranceOperationTraits
    4447  template <
    4548    typename V,
    4649    bool has_inf = std::numeric_limits<V>::has_infinity>
    4750  struct BellmanFordDefaultOperationTraits {
    48     /// \e
     51    /// \brief Value type for the algorithm.
    4952    typedef V Value;
    5053    /// \brief Gives back the zero value of the type.
     
    8588  };
    8689 
     90  /// \brief Operation traits for the BellmanFord algorithm class
     91  /// using tolerance.
     92  ///
     93  /// This operation traits class defines all computational operations
     94  /// and constants that are used in the Bellman-Ford algorithm.
     95  /// The only difference between this implementation and
     96  /// \ref BellmanFordDefaultOperationTraits is that this class uses
     97  /// the \ref Tolerance "tolerance technique" in its \ref less()
     98  /// function.
     99  ///
     100  /// \tparam V The value type.
     101  /// \tparam eps The epsilon value for the \ref less() function.
     102  /// By default, it is the epsilon value used by \ref Tolerance
     103  /// "Tolerance<V>".
     104  ///
     105  /// \see BellmanFordDefaultOperationTraits
     106#ifdef DOXYGEN
     107  template <typename V, V eps>
     108#else
     109  template <
     110    typename V,
     111    V eps = Tolerance<V>::def_epsilon>
     112#endif
     113  struct BellmanFordToleranceOperationTraits {
     114    /// \brief Value type for the algorithm.
     115    typedef V Value;
     116    /// \brief Gives back the zero value of the type.
     117    static Value zero() {
     118      return static_cast<Value>(0);
     119    }
     120    /// \brief Gives back the positive infinity value of the type.
     121    static Value infinity() {
     122      return std::numeric_limits<Value>::infinity();
     123    }
     124    /// \brief Gives back the sum of the given two elements.
     125    static Value plus(const Value& left, const Value& right) {
     126      return left + right;
     127    }
     128    /// \brief Gives back \c true only if the first value is less than
     129    /// the second.
     130    static bool less(const Value& left, const Value& right) {
     131      return left + eps < right;
     132    }
     133  };
     134
    87135  /// \brief Default traits class of BellmanFord class.
    88136  ///
     
    108156    /// It defines the used operations and the infinity value for the
    109157    /// given \c Value type.
    110     /// \see BellmanFordDefaultOperationTraits
     158    /// \see BellmanFordDefaultOperationTraits,
     159    /// BellmanFordToleranceOperationTraits
    111160    typedef BellmanFordDefaultOperationTraits<Value> OperationTraits;
    112161 
     
    838887    /// It defines the used operations and the infinity value for the
    839888    /// given \c Value type.
    840     /// \see BellmanFordDefaultOperationTraits
     889    /// \see BellmanFordDefaultOperationTraits,
     890    /// BellmanFordToleranceOperationTraits
    841891    typedef BellmanFordDefaultOperationTraits<Value> OperationTraits;
    842892
  • lemon/capacity_scaling.h

    r831 r840  
    140140
    141141    typedef std::vector<int> IntVector;
    142     typedef std::vector<char> BoolVector;
    143142    typedef std::vector<Value> ValueVector;
    144143    typedef std::vector<Cost> CostVector;
     144    typedef std::vector<char> BoolVector;
     145    // Note: vector<char> is used instead of vector<bool> for efficiency reasons
    145146
    146147  private:
     
    799800      if (_factor > 1) {
    800801        // With scaling
    801         Value max_sup = 0, max_dem = 0;
    802         for (int i = 0; i != _node_num; ++i) {
     802        Value max_sup = 0, max_dem = 0, max_cap = 0;
     803        for (int i = 0; i != _root; ++i) {
    803804          Value ex = _excess[i];
    804805          if ( ex > max_sup) max_sup =  ex;
    805806          if (-ex > max_dem) max_dem = -ex;
    806         }
    807         Value max_cap = 0;
    808         for (int j = 0; j != _res_arc_num; ++j) {
    809           if (_res_cap[j] > max_cap) max_cap = _res_cap[j];
     807          int last_out = _first_out[i+1] - 1;
     808          for (int j = _first_out[i]; j != last_out; ++j) {
     809            if (_res_cap[j] > max_cap) max_cap = _res_cap[j];
     810          }
    810811        }
    811812        max_sup = std::min(std::min(max_sup, max_dem), max_cap);
  • lemon/cost_scaling.h

    r831 r840  
    202202
    203203    typedef std::vector<int> IntVector;
    204     typedef std::vector<char> BoolVector;
    205204    typedef std::vector<Value> ValueVector;
    206205    typedef std::vector<Cost> CostVector;
    207206    typedef std::vector<LargeCost> LargeCostVector;
     207    typedef std::vector<char> BoolVector;
     208    // Note: vector<char> is used instead of vector<bool> for efficiency reasons
    208209
    209210  private:
     
    249250    bool _have_lower;
    250251    Value _sum_supply;
     252    int _sup_node_num;
    251253
    252254    // Data structures for storing the digraph
     
    277279    int _alpha;
    278280
     281    IntVector _buckets;
     282    IntVector _bucket_next;
     283    IntVector _bucket_prev;
     284    IntVector _rank;
     285    int _max_rank;
     286 
    279287    // Data for a StaticDigraph structure
    280288    typedef std::pair<int, int> IntPair;
     
    829837      }
    830838
     839      _sup_node_num = 0;
     840      for (NodeIt n(_graph); n != INVALID; ++n) {
     841        if (sup[n] > 0) ++_sup_node_num;
     842      }
     843
    831844      // Find a feasible flow using Circulation
    832845      Circulation<Digraph, ConstMap<Arc, Value>, ValueArcMap, ValueNodeMap>
     
    863876        for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
    864877          int ra = _reverse[a];
    865           _res_cap[a] = 1;
     878          _res_cap[a] = 0;
    866879          _res_cap[ra] = 0;
    867880          _cost[a] = 0;
     
    877890      // Maximum path length for partial augment
    878891      const int MAX_PATH_LENGTH = 4;
    879      
     892
     893      // Initialize data structures for buckets     
     894      _max_rank = _alpha * _res_node_num;
     895      _buckets.resize(_max_rank);
     896      _bucket_next.resize(_res_node_num + 1);
     897      _bucket_prev.resize(_res_node_num + 1);
     898      _rank.resize(_res_node_num + 1);
     899 
    880900      // Execute the algorithm
    881901      switch (method) {
     
    916936      }
    917937    }
     938   
     939    // Initialize a cost scaling phase
     940    void initPhase() {
     941      // Saturate arcs not satisfying the optimality condition
     942      for (int u = 0; u != _res_node_num; ++u) {
     943        int last_out = _first_out[u+1];
     944        LargeCost pi_u = _pi[u];
     945        for (int a = _first_out[u]; a != last_out; ++a) {
     946          int v = _target[a];
     947          if (_res_cap[a] > 0 && _cost[a] + pi_u - _pi[v] < 0) {
     948            Value delta = _res_cap[a];
     949            _excess[u] -= delta;
     950            _excess[v] += delta;
     951            _res_cap[a] = 0;
     952            _res_cap[_reverse[a]] += delta;
     953          }
     954        }
     955      }
     956     
     957      // Find active nodes (i.e. nodes with positive excess)
     958      for (int u = 0; u != _res_node_num; ++u) {
     959        if (_excess[u] > 0) _active_nodes.push_back(u);
     960      }
     961
     962      // Initialize the next arcs
     963      for (int u = 0; u != _res_node_num; ++u) {
     964        _next_out[u] = _first_out[u];
     965      }
     966    }
     967   
     968    // Early termination heuristic
     969    bool earlyTermination() {
     970      const double EARLY_TERM_FACTOR = 3.0;
     971
     972      // Build a static residual graph
     973      _arc_vec.clear();
     974      _cost_vec.clear();
     975      for (int j = 0; j != _res_arc_num; ++j) {
     976        if (_res_cap[j] > 0) {
     977          _arc_vec.push_back(IntPair(_source[j], _target[j]));
     978          _cost_vec.push_back(_cost[j] + 1);
     979        }
     980      }
     981      _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end());
     982
     983      // Run Bellman-Ford algorithm to check if the current flow is optimal
     984      BellmanFord<StaticDigraph, LargeCostArcMap> bf(_sgr, _cost_map);
     985      bf.init(0);
     986      bool done = false;
     987      int K = int(EARLY_TERM_FACTOR * std::sqrt(double(_res_node_num)));
     988      for (int i = 0; i < K && !done; ++i) {
     989        done = bf.processNextWeakRound();
     990      }
     991      return done;
     992    }
     993
     994    // Global potential update heuristic
     995    void globalUpdate() {
     996      int bucket_end = _root + 1;
     997   
     998      // Initialize buckets
     999      for (int r = 0; r != _max_rank; ++r) {
     1000        _buckets[r] = bucket_end;
     1001      }
     1002      Value total_excess = 0;
     1003      for (int i = 0; i != _res_node_num; ++i) {
     1004        if (_excess[i] < 0) {
     1005          _rank[i] = 0;
     1006          _bucket_next[i] = _buckets[0];
     1007          _bucket_prev[_buckets[0]] = i;
     1008          _buckets[0] = i;
     1009        } else {
     1010          total_excess += _excess[i];
     1011          _rank[i] = _max_rank;
     1012        }
     1013      }
     1014      if (total_excess == 0) return;
     1015
     1016      // Search the buckets
     1017      int r = 0;
     1018      for ( ; r != _max_rank; ++r) {
     1019        while (_buckets[r] != bucket_end) {
     1020          // Remove the first node from the current bucket
     1021          int u = _buckets[r];
     1022          _buckets[r] = _bucket_next[u];
     1023         
     1024          // Search the incomming arcs of u
     1025          LargeCost pi_u = _pi[u];
     1026          int last_out = _first_out[u+1];
     1027          for (int a = _first_out[u]; a != last_out; ++a) {
     1028            int ra = _reverse[a];
     1029            if (_res_cap[ra] > 0) {
     1030              int v = _source[ra];
     1031              int old_rank_v = _rank[v];
     1032              if (r < old_rank_v) {
     1033                // Compute the new rank of v
     1034                LargeCost nrc = (_cost[ra] + _pi[v] - pi_u) / _epsilon;
     1035                int new_rank_v = old_rank_v;
     1036                if (nrc < LargeCost(_max_rank))
     1037                  new_rank_v = r + 1 + int(nrc);
     1038                 
     1039                // Change the rank of v
     1040                if (new_rank_v < old_rank_v) {
     1041                  _rank[v] = new_rank_v;
     1042                  _next_out[v] = _first_out[v];
     1043                 
     1044                  // Remove v from its old bucket
     1045                  if (old_rank_v < _max_rank) {
     1046                    if (_buckets[old_rank_v] == v) {
     1047                      _buckets[old_rank_v] = _bucket_next[v];
     1048                    } else {
     1049                      _bucket_next[_bucket_prev[v]] = _bucket_next[v];
     1050                      _bucket_prev[_bucket_next[v]] = _bucket_prev[v];
     1051                    }
     1052                  }
     1053                 
     1054                  // Insert v to its new bucket
     1055                  _bucket_next[v] = _buckets[new_rank_v];
     1056                  _bucket_prev[_buckets[new_rank_v]] = v;
     1057                  _buckets[new_rank_v] = v;
     1058                }
     1059              }
     1060            }
     1061          }
     1062
     1063          // Finish search if there are no more active nodes
     1064          if (_excess[u] > 0) {
     1065            total_excess -= _excess[u];
     1066            if (total_excess <= 0) break;
     1067          }
     1068        }
     1069        if (total_excess <= 0) break;
     1070      }
     1071     
     1072      // Relabel nodes
     1073      for (int u = 0; u != _res_node_num; ++u) {
     1074        int k = std::min(_rank[u], r);
     1075        if (k > 0) {
     1076          _pi[u] -= _epsilon * k;
     1077          _next_out[u] = _first_out[u];
     1078        }
     1079      }
     1080    }
    9181081
    9191082    /// Execute the algorithm performing augment and relabel operations
    9201083    void startAugment(int max_length = std::numeric_limits<int>::max()) {
    9211084      // Paramters for heuristics
    922       const int BF_HEURISTIC_EPSILON_BOUND = 1000;
    923       const int BF_HEURISTIC_BOUND_FACTOR  = 3;
    924 
     1085      const int EARLY_TERM_EPSILON_LIMIT = 1000;
     1086      const double GLOBAL_UPDATE_FACTOR = 3.0;
     1087
     1088      const int global_update_freq = int(GLOBAL_UPDATE_FACTOR *
     1089        (_res_node_num + _sup_node_num * _sup_node_num));
     1090      int next_update_limit = global_update_freq;
     1091     
     1092      int relabel_cnt = 0;
     1093     
    9251094      // Perform cost scaling phases
    926       IntVector pred_arc(_res_node_num);
    927       std::vector<int> path_nodes;
     1095      std::vector<int> path;
    9281096      for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?
    9291097                                        1 : _epsilon / _alpha )
    9301098      {
    931         // "Early Termination" heuristic: use Bellman-Ford algorithm
    932         // to check if the current flow is optimal
    933         if (_epsilon <= BF_HEURISTIC_EPSILON_BOUND) {
    934           _arc_vec.clear();
    935           _cost_vec.clear();
    936           for (int j = 0; j != _res_arc_num; ++j) {
    937             if (_res_cap[j] > 0) {
    938               _arc_vec.push_back(IntPair(_source[j], _target[j]));
    939               _cost_vec.push_back(_cost[j] + 1);
    940             }
    941           }
    942           _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end());
    943 
    944           BellmanFord<StaticDigraph, LargeCostArcMap> bf(_sgr, _cost_map);
    945           bf.init(0);
    946           bool done = false;
    947           int K = int(BF_HEURISTIC_BOUND_FACTOR * sqrt(_res_node_num));
    948           for (int i = 0; i < K && !done; ++i)
    949             done = bf.processNextWeakRound();
    950           if (done) break;
    951         }
    952 
    953         // Saturate arcs not satisfying the optimality condition
    954         for (int a = 0; a != _res_arc_num; ++a) {
    955           if (_res_cap[a] > 0 &&
    956               _cost[a] + _pi[_source[a]] - _pi[_target[a]] < 0) {
    957             Value delta = _res_cap[a];
    958             _excess[_source[a]] -= delta;
    959             _excess[_target[a]] += delta;
    960             _res_cap[a] = 0;
    961             _res_cap[_reverse[a]] += delta;
    962           }
     1099        // Early termination heuristic
     1100        if (_epsilon <= EARLY_TERM_EPSILON_LIMIT) {
     1101          if (earlyTermination()) break;
    9631102        }
    9641103       
    965         // Find active nodes (i.e. nodes with positive excess)
    966         for (int u = 0; u != _res_node_num; ++u) {
    967           if (_excess[u] > 0) _active_nodes.push_back(u);
    968         }
    969 
    970         // Initialize the next arcs
    971         for (int u = 0; u != _res_node_num; ++u) {
    972           _next_out[u] = _first_out[u];
    973         }
    974 
     1104        // Initialize current phase
     1105        initPhase();
     1106       
    9751107        // Perform partial augment and relabel operations
    9761108        while (true) {
     
    9821114          if (_active_nodes.size() == 0) break;
    9831115          int start = _active_nodes.front();
    984           path_nodes.clear();
    985           path_nodes.push_back(start);
    9861116
    9871117          // Find an augmenting path from the start node
     1118          path.clear();
    9881119          int tip = start;
    989           while (_excess[tip] >= 0 &&
    990                  int(path_nodes.size()) <= max_length) {
     1120          while (_excess[tip] >= 0 && int(path.size()) < max_length) {
    9911121            int u;
    992             LargeCost min_red_cost, rc;
    993             int last_out = _sum_supply < 0 ?
    994               _first_out[tip+1] : _first_out[tip+1] - 1;
     1122            LargeCost min_red_cost, rc, pi_tip = _pi[tip];
     1123            int last_out = _first_out[tip+1];
    9951124            for (int a = _next_out[tip]; a != last_out; ++a) {
    996               if (_res_cap[a] > 0 &&
    997                   _cost[a] + _pi[_source[a]] - _pi[_target[a]] < 0) {
    998                 u = _target[a];
    999                 pred_arc[u] = a;
     1125              u = _target[a];
     1126              if (_res_cap[a] > 0 && _cost[a] + pi_tip - _pi[u] < 0) {
     1127                path.push_back(a);
    10001128                _next_out[tip] = a;
    10011129                tip = u;
    1002                 path_nodes.push_back(tip);
    10031130                goto next_step;
    10041131              }
     
    10061133
    10071134            // Relabel tip node
    1008             min_red_cost = std::numeric_limits<LargeCost>::max() / 2;
     1135            min_red_cost = std::numeric_limits<LargeCost>::max();
     1136            if (tip != start) {
     1137              int ra = _reverse[path.back()];
     1138              min_red_cost = _cost[ra] + pi_tip - _pi[_target[ra]];
     1139            }
    10091140            for (int a = _first_out[tip]; a != last_out; ++a) {
    1010               rc = _cost[a] + _pi[_source[a]] - _pi[_target[a]];
     1141              rc = _cost[a] + pi_tip - _pi[_target[a]];
    10111142              if (_res_cap[a] > 0 && rc < min_red_cost) {
    10121143                min_red_cost = rc;
     
    10141145            }
    10151146            _pi[tip] -= min_red_cost + _epsilon;
    1016 
    1017             // Reset the next arc of tip
    10181147            _next_out[tip] = _first_out[tip];
     1148            ++relabel_cnt;
    10191149
    10201150            // Step back
    10211151            if (tip != start) {
    1022               path_nodes.pop_back();
    1023               tip = path_nodes.back();
     1152              tip = _source[path.back()];
     1153              path.pop_back();
    10241154            }
    10251155
     
    10291159          // Augment along the found path (as much flow as possible)
    10301160          Value delta;
    1031           int u, v = path_nodes.front(), pa;
    1032           for (int i = 1; i < int(path_nodes.size()); ++i) {
     1161          int pa, u, v = start;
     1162          for (int i = 0; i != int(path.size()); ++i) {
     1163            pa = path[i];
    10331164            u = v;
    1034             v = path_nodes[i];
    1035             pa = pred_arc[v];
     1165            v = _target[pa];
    10361166            delta = std::min(_res_cap[pa], _excess[u]);
    10371167            _res_cap[pa] -= delta;
     
    10421172              _active_nodes.push_back(v);
    10431173          }
     1174
     1175          // Global update heuristic
     1176          if (relabel_cnt >= next_update_limit) {
     1177            globalUpdate();
     1178            next_update_limit += global_update_freq;
     1179          }
    10441180        }
    10451181      }
     
    10491185    void startPush() {
    10501186      // Paramters for heuristics
    1051       const int BF_HEURISTIC_EPSILON_BOUND = 1000;
    1052       const int BF_HEURISTIC_BOUND_FACTOR  = 3;
    1053 
     1187      const int EARLY_TERM_EPSILON_LIMIT = 1000;
     1188      const double GLOBAL_UPDATE_FACTOR = 2.0;
     1189
     1190      const int global_update_freq = int(GLOBAL_UPDATE_FACTOR *
     1191        (_res_node_num + _sup_node_num * _sup_node_num));
     1192      int next_update_limit = global_update_freq;
     1193
     1194      int relabel_cnt = 0;
     1195     
    10541196      // Perform cost scaling phases
    10551197      BoolVector hyper(_res_node_num, false);
     1198      LargeCostVector hyper_cost(_res_node_num);
    10561199      for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?
    10571200                                        1 : _epsilon / _alpha )
    10581201      {
    1059         // "Early Termination" heuristic: use Bellman-Ford algorithm
    1060         // to check if the current flow is optimal
    1061         if (_epsilon <= BF_HEURISTIC_EPSILON_BOUND) {
    1062           _arc_vec.clear();
    1063           _cost_vec.clear();
    1064           for (int j = 0; j != _res_arc_num; ++j) {
    1065             if (_res_cap[j] > 0) {
    1066               _arc_vec.push_back(IntPair(_source[j], _target[j]));
    1067               _cost_vec.push_back(_cost[j] + 1);
    1068             }
    1069           }
    1070           _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end());
    1071 
    1072           BellmanFord<StaticDigraph, LargeCostArcMap> bf(_sgr, _cost_map);
    1073           bf.init(0);
    1074           bool done = false;
    1075           int K = int(BF_HEURISTIC_BOUND_FACTOR * sqrt(_res_node_num));
    1076           for (int i = 0; i < K && !done; ++i)
    1077             done = bf.processNextWeakRound();
    1078           if (done) break;
    1079         }
    1080 
    1081         // Saturate arcs not satisfying the optimality condition
    1082         for (int a = 0; a != _res_arc_num; ++a) {
    1083           if (_res_cap[a] > 0 &&
    1084               _cost[a] + _pi[_source[a]] - _pi[_target[a]] < 0) {
    1085             Value delta = _res_cap[a];
    1086             _excess[_source[a]] -= delta;
    1087             _excess[_target[a]] += delta;
    1088             _res_cap[a] = 0;
    1089             _res_cap[_reverse[a]] += delta;
    1090           }
    1091         }
    1092 
    1093         // Find active nodes (i.e. nodes with positive excess)
    1094         for (int u = 0; u != _res_node_num; ++u) {
    1095           if (_excess[u] > 0) _active_nodes.push_back(u);
    1096         }
    1097 
    1098         // Initialize the next arcs
    1099         for (int u = 0; u != _res_node_num; ++u) {
    1100           _next_out[u] = _first_out[u];
    1101         }
     1202        // Early termination heuristic
     1203        if (_epsilon <= EARLY_TERM_EPSILON_LIMIT) {
     1204          if (earlyTermination()) break;
     1205        }
     1206       
     1207        // Initialize current phase
     1208        initPhase();
    11021209
    11031210        // Perform push and relabel operations
    11041211        while (_active_nodes.size() > 0) {
    1105           LargeCost min_red_cost, rc;
     1212          LargeCost min_red_cost, rc, pi_n;
    11061213          Value delta;
    11071214          int n, t, a, last_out = _res_arc_num;
    11081215
     1216        next_node:
    11091217          // Select an active node (FIFO selection)
    1110         next_node:
    11111218          n = _active_nodes.front();
    1112           last_out = _sum_supply < 0 ?
    1113             _first_out[n+1] : _first_out[n+1] - 1;
    1114 
     1219          last_out = _first_out[n+1];
     1220          pi_n = _pi[n];
     1221         
    11151222          // Perform push operations if there are admissible arcs
    11161223          if (_excess[n] > 0) {
    11171224            for (a = _next_out[n]; a != last_out; ++a) {
    11181225              if (_res_cap[a] > 0 &&
    1119                   _cost[a] + _pi[_source[a]] - _pi[_target[a]] < 0) {
     1226                  _cost[a] + pi_n - _pi[_target[a]] < 0) {
    11201227                delta = std::min(_res_cap[a], _excess[n]);
    11211228                t = _target[a];
     
    11231230                // Push-look-ahead heuristic
    11241231                Value ahead = -_excess[t];
    1125                 int last_out_t = _sum_supply < 0 ?
    1126                   _first_out[t+1] : _first_out[t+1] - 1;
     1232                int last_out_t = _first_out[t+1];
     1233                LargeCost pi_t = _pi[t];
    11271234                for (int ta = _next_out[t]; ta != last_out_t; ++ta) {
    11281235                  if (_res_cap[ta] > 0 &&
    1129                       _cost[ta] + _pi[_source[ta]] - _pi[_target[ta]] < 0)
     1236                      _cost[ta] + pi_t - _pi[_target[ta]] < 0)
    11301237                    ahead += _res_cap[ta];
    11311238                  if (ahead >= delta) break;
     
    11341241
    11351242                // Push flow along the arc
    1136                 if (ahead < delta) {
     1243                if (ahead < delta && !hyper[t]) {
    11371244                  _res_cap[a] -= ahead;
    11381245                  _res_cap[_reverse[a]] += ahead;
     
    11411248                  _active_nodes.push_front(t);
    11421249                  hyper[t] = true;
     1250                  hyper_cost[t] = _cost[a] + pi_n - pi_t;
    11431251                  _next_out[n] = a;
    11441252                  goto next_node;
     
    11631271          // Relabel the node if it is still active (or hyper)
    11641272          if (_excess[n] > 0 || hyper[n]) {
    1165             min_red_cost = std::numeric_limits<LargeCost>::max() / 2;
     1273             min_red_cost = hyper[n] ? -hyper_cost[n] :
     1274               std::numeric_limits<LargeCost>::max();
    11661275            for (int a = _first_out[n]; a != last_out; ++a) {
    1167               rc = _cost[a] + _pi[_source[a]] - _pi[_target[a]];
     1276              rc = _cost[a] + pi_n - _pi[_target[a]];
    11681277              if (_res_cap[a] > 0 && rc < min_red_cost) {
    11691278                min_red_cost = rc;
     
    11711280            }
    11721281            _pi[n] -= min_red_cost + _epsilon;
     1282            _next_out[n] = _first_out[n];
    11731283            hyper[n] = false;
    1174 
    1175             // Reset the next arc
    1176             _next_out[n] = _first_out[n];
     1284            ++relabel_cnt;
    11771285          }
    11781286       
     
    11841292            _active_nodes.pop_front();
    11851293          }
     1294         
     1295          // Global update heuristic
     1296          if (relabel_cnt >= next_update_limit) {
     1297            globalUpdate();
     1298            for (int u = 0; u != _res_node_num; ++u)
     1299              hyper[u] = false;
     1300            next_update_limit += global_update_freq;
     1301          }
    11861302        }
    11871303      }
  • lemon/cycle_canceling.h

    r830 r840  
    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;
     
    963964      DoubleVector pi(_res_node_num, 0.0);
    964965      IntVector level(_res_node_num);
    965       CharVector reached(_res_node_num);
    966       CharVector processed(_res_node_num);
     966      BoolVector reached(_res_node_num);
     967      BoolVector processed(_res_node_num);
    967968      IntVector pred_node(_res_node_num);
    968969      IntVector pred_arc(_res_node_num);
  • lemon/graph_to_eps.h

    r786 r838  
    685685#else
    686686      os << bits::getWinFormattedDate();
     687      os << std::endl;
    687688#endif
    688689    }
    689     os << std::endl;
    690690
    691691    if (_autoArcWidthScale) {
  • lemon/hartmann_orlin.h

    r825 r841  
    406406    /// \pre \ref run() or \ref findMinMean() must be called before
    407407    /// using this function.
    408     LargeValue cycleLength() const {
    409       return _best_length;
     408    Value cycleLength() const {
     409      return static_cast<Value>(_best_length);
    410410    }
    411411
  • lemon/howard.h

    r825 r841  
    385385    /// \pre \ref run() or \ref findMinMean() must be called before
    386386    /// using this function.
    387     LargeValue cycleLength() const {
    388       return _best_length;
     387    Value cycleLength() const {
     388      return static_cast<Value>(_best_length);
    389389    }
    390390
  • lemon/karp.h

    r825 r841  
    393393    /// \pre \ref run() or \ref findMinMean() must be called before
    394394    /// using this function.
    395     LargeValue cycleLength() const {
    396       return _cycle_length;
     395    Value cycleLength() const {
     396      return static_cast<Value>(_cycle_length);
    397397    }
    398398
  • lemon/network_simplex.h

    r830 r840  
    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
     
    214215    IntVector _last_succ;
    215216    IntVector _dirty_revs;
    216     CharVector _forward;
    217     CharVector _state;
     217    BoolVector _forward;
     218    BoolVector _state;
    218219    int _root;
    219220
     
    246247      const IntVector  &_target;
    247248      const CostVector &_cost;
    248       const CharVector &_state;
     249      const BoolVector &_state;
    249250      const CostVector &_pi;
    250251      int &_in_arc;
     
    267268      bool findEnteringArc() {
    268269        Cost c;
    269         for (int e = _next_arc; e < _search_arc_num; ++e) {
     270        for (int e = _next_arc; e != _search_arc_num; ++e) {
    270271          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
    271272          if (c < 0) {
     
    275276          }
    276277        }
    277         for (int e = 0; e < _next_arc; ++e) {
     278        for (int e = 0; e != _next_arc; ++e) {
    278279          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
    279280          if (c < 0) {
     
    298299      const IntVector  &_target;
    299300      const CostVector &_cost;
    300       const CharVector &_state;
     301      const BoolVector &_state;
    301302      const CostVector &_pi;
    302303      int &_in_arc;
     
    315316      bool findEnteringArc() {
    316317        Cost c, min = 0;
    317         for (int e = 0; e < _search_arc_num; ++e) {
     318        for (int e = 0; e != _search_arc_num; ++e) {
    318319          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
    319320          if (c < min) {
     
    337338      const IntVector  &_target;
    338339      const CostVector &_cost;
    339       const CharVector &_state;
     340      const BoolVector &_state;
    340341      const CostVector &_pi;
    341342      int &_in_arc;
     
    356357      {
    357358        // The main parameters of the pivot rule
    358         const double BLOCK_SIZE_FACTOR = 0.5;
     359        const double BLOCK_SIZE_FACTOR = 1.0;
    359360        const int MIN_BLOCK_SIZE = 10;
    360361
     
    369370        int cnt = _block_size;
    370371        int e;
    371         for (e = _next_arc; e < _search_arc_num; ++e) {
     372        for (e = _next_arc; e != _search_arc_num; ++e) {
    372373          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
    373374          if (c < min) {
     
    380381          }
    381382        }
    382         for (e = 0; e < _next_arc; ++e) {
     383        for (e = 0; e != _next_arc; ++e) {
    383384          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
    384385          if (c < min) {
     
    410411      const IntVector  &_target;
    411412      const CostVector &_cost;
    412       const CharVector &_state;
     413      const BoolVector &_state;
    413414      const CostVector &_pi;
    414415      int &_in_arc;
     
    471472        min = 0;
    472473        _curr_length = 0;
    473         for (e = _next_arc; e < _search_arc_num; ++e) {
     474        for (e = _next_arc; e != _search_arc_num; ++e) {
    474475          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
    475476          if (c < 0) {
     
    482483          }
    483484        }
    484         for (e = 0; e < _next_arc; ++e) {
     485        for (e = 0; e != _next_arc; ++e) {
    485486          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
    486487          if (c < 0) {
     
    513514      const IntVector  &_target;
    514515      const CostVector &_cost;
    515       const CharVector &_state;
     516      const BoolVector &_state;
    516517      const CostVector &_pi;
    517518      int &_in_arc;
     
    566567        // Check the current candidate list
    567568        int e;
    568         for (int i = 0; i < _curr_length; ++i) {
     569        for (int i = 0; i != _curr_length; ++i) {
    569570          e = _candidates[i];
    570571          _cand_cost[e] = _state[e] *
     
    579580        int limit = _head_length;
    580581
    581         for (e = _next_arc; e < _search_arc_num; ++e) {
     582        for (e = _next_arc; e != _search_arc_num; ++e) {
    582583          _cand_cost[e] = _state[e] *
    583584            (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
     
    591592          }
    592593        }
    593         for (e = 0; e < _next_arc; ++e) {
     594        for (e = 0; e != _next_arc; ++e) {
    594595          _cand_cost[e] = _state[e] *
    595596            (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
     
    13611362
    13621363      // Update _rev_thread using the new _thread values
    1363       for (int i = 0; i < int(_dirty_revs.size()); ++i) {
     1364      for (int i = 0; i != int(_dirty_revs.size()); ++i) {
    13641365        u = _dirty_revs[i];
    13651366        _rev_thread[_thread[u]] = u;
     
    14311432        _pi[u] += sigma;
    14321433      }
     1434    }
     1435
     1436    // Heuristic initial pivots
     1437    bool initialPivots() {
     1438      Value curr, total = 0;
     1439      std::vector<Node> supply_nodes, demand_nodes;
     1440      for (NodeIt u(_graph); u != INVALID; ++u) {
     1441        curr = _supply[_node_id[u]];
     1442        if (curr > 0) {
     1443          total += curr;
     1444          supply_nodes.push_back(u);
     1445        }
     1446        else if (curr < 0) {
     1447          demand_nodes.push_back(u);
     1448        }
     1449      }
     1450      if (_sum_supply > 0) total -= _sum_supply;
     1451      if (total <= 0) return true;
     1452
     1453      IntVector arc_vector;
     1454      if (_sum_supply >= 0) {
     1455        if (supply_nodes.size() == 1 && demand_nodes.size() == 1) {
     1456          // Perform a reverse graph search from the sink to the source
     1457          typename GR::template NodeMap<bool> reached(_graph, false);
     1458          Node s = supply_nodes[0], t = demand_nodes[0];
     1459          std::vector<Node> stack;
     1460          reached[t] = true;
     1461          stack.push_back(t);
     1462          while (!stack.empty()) {
     1463            Node u, v = stack.back();
     1464            stack.pop_back();
     1465            if (v == s) break;
     1466            for (InArcIt a(_graph, v); a != INVALID; ++a) {
     1467              if (reached[u = _graph.source(a)]) continue;
     1468              int j = _arc_id[a];
     1469              if (_cap[j] >= total) {
     1470                arc_vector.push_back(j);
     1471                reached[u] = true;
     1472                stack.push_back(u);
     1473              }
     1474            }
     1475          }
     1476        } else {
     1477          // Find the min. cost incomming arc for each demand node
     1478          for (int i = 0; i != int(demand_nodes.size()); ++i) {
     1479            Node v = demand_nodes[i];
     1480            Cost c, min_cost = std::numeric_limits<Cost>::max();
     1481            Arc min_arc = INVALID;
     1482            for (InArcIt a(_graph, v); a != INVALID; ++a) {
     1483              c = _cost[_arc_id[a]];
     1484              if (c < min_cost) {
     1485                min_cost = c;
     1486                min_arc = a;
     1487              }
     1488            }
     1489            if (min_arc != INVALID) {
     1490              arc_vector.push_back(_arc_id[min_arc]);
     1491            }
     1492          }
     1493        }
     1494      } else {
     1495        // Find the min. cost outgoing arc for each supply node
     1496        for (int i = 0; i != int(supply_nodes.size()); ++i) {
     1497          Node u = supply_nodes[i];
     1498          Cost c, min_cost = std::numeric_limits<Cost>::max();
     1499          Arc min_arc = INVALID;
     1500          for (OutArcIt a(_graph, u); a != INVALID; ++a) {
     1501            c = _cost[_arc_id[a]];
     1502            if (c < min_cost) {
     1503              min_cost = c;
     1504              min_arc = a;
     1505            }
     1506          }
     1507          if (min_arc != INVALID) {
     1508            arc_vector.push_back(_arc_id[min_arc]);
     1509          }
     1510        }
     1511      }
     1512
     1513      // Perform heuristic initial pivots
     1514      for (int i = 0; i != int(arc_vector.size()); ++i) {
     1515        in_arc = arc_vector[i];
     1516        if (_state[in_arc] * (_cost[in_arc] + _pi[_source[in_arc]] -
     1517            _pi[_target[in_arc]]) >= 0) continue;
     1518        findJoinNode();
     1519        bool change = findLeavingArc();
     1520        if (delta >= MAX) return false;
     1521        changeFlow(change);
     1522        if (change) {
     1523          updateTreeStructure();
     1524          updatePotential();
     1525        }
     1526      }
     1527      return true;
    14331528    }
    14341529
     
    14551550      PivotRuleImpl pivot(*this);
    14561551
     1552      // Perform heuristic initial pivots
     1553      if (!initialPivots()) return UNBOUNDED;
     1554
    14571555      // Execute the Network Simplex algorithm
    14581556      while (pivot.findEnteringArc()) {
  • test/bellman_ford_test.cc

    r791 r844  
    105105      ::SetDistMap<concepts::ReadWriteMap<Node,Value> >
    106106      ::SetOperationTraits<BellmanFordDefaultOperationTraits<Value> >
     107      ::SetOperationTraits<BellmanFordToleranceOperationTraits<Value, 0> >
    107108      ::Create bf_test(gr,length);
    108109
Note: See TracChangeset for help on using the changeset viewer.