COIN-OR::LEMON - Graph Library

Changeset 840:2914b6f0fde0 in lemon-1.2


Ignore:
Timestamp:
02/26/10 14:00:20 (9 years ago)
Author:
Alpar Juttner <alpar@…>
Branch:
default
Parents:
838:2c35bef44dd1 (diff), 839:f3bc4e9b5f3a (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
Phase:
public
Message:

Merge #340

Location:
lemon
Files:
8 edited

Legend:

Unmodified
Added
Removed
  • 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/capacity_scaling.h

    r839 r840  
    7878  /// \tparam GR The digraph type the algorithm runs on.
    7979  /// \tparam V The number type used for flow amounts, capacity bounds
    80   /// and supply values in the algorithm. By default it is \c int.
     80  /// and supply values in the algorithm. By default, it is \c int.
    8181  /// \tparam C The number type used for costs and potentials in the
    82   /// algorithm. By default it is the same as \c V.
     82  /// algorithm. By default, it is the same as \c V.
     83  /// \tparam TR The traits class that defines various types used by the
     84  /// algorithm. By default, it is \ref CapacityScalingDefaultTraits
     85  /// "CapacityScalingDefaultTraits<GR, V, C>".
     86  /// In most cases, this parameter should not be set directly,
     87  /// consider to use the named template parameters instead.
    8388  ///
    8489  /// \warning Both number types must be signed and all input data must
     
    316321        "The cost type of CapacityScaling must be signed");
    317322
     323      // Reset data structures
     324      reset();
     325    }
     326
     327    /// \name Parameters
     328    /// The parameters of the algorithm can be specified using these
     329    /// functions.
     330
     331    /// @{
     332
     333    /// \brief Set the lower bounds on the arcs.
     334    ///
     335    /// This function sets the lower bounds on the arcs.
     336    /// If it is not used before calling \ref run(), the lower bounds
     337    /// will be set to zero on all arcs.
     338    ///
     339    /// \param map An arc map storing the lower bounds.
     340    /// Its \c Value type must be convertible to the \c Value type
     341    /// of the algorithm.
     342    ///
     343    /// \return <tt>(*this)</tt>
     344    template <typename LowerMap>
     345    CapacityScaling& lowerMap(const LowerMap& map) {
     346      _have_lower = true;
     347      for (ArcIt a(_graph); a != INVALID; ++a) {
     348        _lower[_arc_idf[a]] = map[a];
     349        _lower[_arc_idb[a]] = map[a];
     350      }
     351      return *this;
     352    }
     353
     354    /// \brief Set the upper bounds (capacities) on the arcs.
     355    ///
     356    /// This function sets the upper bounds (capacities) on the arcs.
     357    /// If it is not used before calling \ref run(), the upper bounds
     358    /// will be set to \ref INF on all arcs (i.e. the flow value will be
     359    /// unbounded from above).
     360    ///
     361    /// \param map An arc map storing the upper bounds.
     362    /// Its \c Value type must be convertible to the \c Value type
     363    /// of the algorithm.
     364    ///
     365    /// \return <tt>(*this)</tt>
     366    template<typename UpperMap>
     367    CapacityScaling& upperMap(const UpperMap& map) {
     368      for (ArcIt a(_graph); a != INVALID; ++a) {
     369        _upper[_arc_idf[a]] = map[a];
     370      }
     371      return *this;
     372    }
     373
     374    /// \brief Set the costs of the arcs.
     375    ///
     376    /// This function sets the costs of the arcs.
     377    /// If it is not used before calling \ref run(), the costs
     378    /// will be set to \c 1 on all arcs.
     379    ///
     380    /// \param map An arc map storing the costs.
     381    /// Its \c Value type must be convertible to the \c Cost type
     382    /// of the algorithm.
     383    ///
     384    /// \return <tt>(*this)</tt>
     385    template<typename CostMap>
     386    CapacityScaling& costMap(const CostMap& map) {
     387      for (ArcIt a(_graph); a != INVALID; ++a) {
     388        _cost[_arc_idf[a]] =  map[a];
     389        _cost[_arc_idb[a]] = -map[a];
     390      }
     391      return *this;
     392    }
     393
     394    /// \brief Set the supply values of the nodes.
     395    ///
     396    /// This function sets the supply values of the nodes.
     397    /// If neither this function nor \ref stSupply() is used before
     398    /// calling \ref run(), the supply of each node will be set to zero.
     399    ///
     400    /// \param map A node map storing the supply values.
     401    /// Its \c Value type must be convertible to the \c Value type
     402    /// of the algorithm.
     403    ///
     404    /// \return <tt>(*this)</tt>
     405    template<typename SupplyMap>
     406    CapacityScaling& supplyMap(const SupplyMap& map) {
     407      for (NodeIt n(_graph); n != INVALID; ++n) {
     408        _supply[_node_id[n]] = map[n];
     409      }
     410      return *this;
     411    }
     412
     413    /// \brief Set single source and target nodes and a supply value.
     414    ///
     415    /// This function sets a single source node and a single target node
     416    /// and the required flow value.
     417    /// If neither this function nor \ref supplyMap() is used before
     418    /// calling \ref run(), the supply of each node will be set to zero.
     419    ///
     420    /// Using this function has the same effect as using \ref supplyMap()
     421    /// with such a map in which \c k is assigned to \c s, \c -k is
     422    /// assigned to \c t and all other nodes have zero supply value.
     423    ///
     424    /// \param s The source node.
     425    /// \param t The target node.
     426    /// \param k The required amount of flow from node \c s to node \c t
     427    /// (i.e. the supply of \c s and the demand of \c t).
     428    ///
     429    /// \return <tt>(*this)</tt>
     430    CapacityScaling& stSupply(const Node& s, const Node& t, Value k) {
     431      for (int i = 0; i != _node_num; ++i) {
     432        _supply[i] = 0;
     433      }
     434      _supply[_node_id[s]] =  k;
     435      _supply[_node_id[t]] = -k;
     436      return *this;
     437    }
     438   
     439    /// @}
     440
     441    /// \name Execution control
     442    /// The algorithm can be executed using \ref run().
     443
     444    /// @{
     445
     446    /// \brief Run the algorithm.
     447    ///
     448    /// This function runs the algorithm.
     449    /// The paramters can be specified using functions \ref lowerMap(),
     450    /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply().
     451    /// For example,
     452    /// \code
     453    ///   CapacityScaling<ListDigraph> cs(graph);
     454    ///   cs.lowerMap(lower).upperMap(upper).costMap(cost)
     455    ///     .supplyMap(sup).run();
     456    /// \endcode
     457    ///
     458    /// This function can be called more than once. All the given parameters
     459    /// are kept for the next call, unless \ref resetParams() or \ref reset()
     460    /// is used, thus only the modified parameters have to be set again.
     461    /// If the underlying digraph was also modified after the construction
     462    /// of the class (or the last \ref reset() call), then the \ref reset()
     463    /// function must be called.
     464    ///
     465    /// \param factor The capacity scaling factor. It must be larger than
     466    /// one to use scaling. If it is less or equal to one, then scaling
     467    /// will be disabled.
     468    ///
     469    /// \return \c INFEASIBLE if no feasible flow exists,
     470    /// \n \c OPTIMAL if the problem has optimal solution
     471    /// (i.e. it is feasible and bounded), and the algorithm has found
     472    /// optimal flow and node potentials (primal and dual solutions),
     473    /// \n \c UNBOUNDED if the digraph contains an arc of negative cost
     474    /// and infinite upper bound. It means that the objective function
     475    /// is unbounded on that arc, however, note that it could actually be
     476    /// bounded over the feasible flows, but this algroithm cannot handle
     477    /// these cases.
     478    ///
     479    /// \see ProblemType
     480    /// \see resetParams(), reset()
     481    ProblemType run(int factor = 4) {
     482      _factor = factor;
     483      ProblemType pt = init();
     484      if (pt != OPTIMAL) return pt;
     485      return start();
     486    }
     487
     488    /// \brief Reset all the parameters that have been given before.
     489    ///
     490    /// This function resets all the paramaters that have been given
     491    /// before using functions \ref lowerMap(), \ref upperMap(),
     492    /// \ref costMap(), \ref supplyMap(), \ref stSupply().
     493    ///
     494    /// It is useful for multiple \ref run() calls. Basically, all the given
     495    /// parameters are kept for the next \ref run() call, unless
     496    /// \ref resetParams() or \ref reset() is used.
     497    /// If the underlying digraph was also modified after the construction
     498    /// of the class or the last \ref reset() call, then the \ref reset()
     499    /// function must be used, otherwise \ref resetParams() is sufficient.
     500    ///
     501    /// For example,
     502    /// \code
     503    ///   CapacityScaling<ListDigraph> cs(graph);
     504    ///
     505    ///   // First run
     506    ///   cs.lowerMap(lower).upperMap(upper).costMap(cost)
     507    ///     .supplyMap(sup).run();
     508    ///
     509    ///   // Run again with modified cost map (resetParams() is not called,
     510    ///   // so only the cost map have to be set again)
     511    ///   cost[e] += 100;
     512    ///   cs.costMap(cost).run();
     513    ///
     514    ///   // Run again from scratch using resetParams()
     515    ///   // (the lower bounds will be set to zero on all arcs)
     516    ///   cs.resetParams();
     517    ///   cs.upperMap(capacity).costMap(cost)
     518    ///     .supplyMap(sup).run();
     519    /// \endcode
     520    ///
     521    /// \return <tt>(*this)</tt>
     522    ///
     523    /// \see reset(), run()
     524    CapacityScaling& resetParams() {
     525      for (int i = 0; i != _node_num; ++i) {
     526        _supply[i] = 0;
     527      }
     528      for (int j = 0; j != _res_arc_num; ++j) {
     529        _lower[j] = 0;
     530        _upper[j] = INF;
     531        _cost[j] = _forward[j] ? 1 : -1;
     532      }
     533      _have_lower = false;
     534      return *this;
     535    }
     536
     537    /// \brief Reset the internal data structures and all the parameters
     538    /// that have been given before.
     539    ///
     540    /// This function resets the internal data structures and all the
     541    /// paramaters that have been given before using functions \ref lowerMap(),
     542    /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply().
     543    ///
     544    /// It is useful for multiple \ref run() calls. Basically, all the given
     545    /// parameters are kept for the next \ref run() call, unless
     546    /// \ref resetParams() or \ref reset() is used.
     547    /// If the underlying digraph was also modified after the construction
     548    /// of the class or the last \ref reset() call, then the \ref reset()
     549    /// function must be used, otherwise \ref resetParams() is sufficient.
     550    ///
     551    /// See \ref resetParams() for examples.
     552    ///
     553    /// \return <tt>(*this)</tt>
     554    ///
     555    /// \see resetParams(), run()
     556    CapacityScaling& reset() {
    318557      // Resize vectors
    319558      _node_num = countNodes(_graph);
     
    379618     
    380619      // Reset parameters
    381       reset();
    382     }
    383 
    384     /// \name Parameters
    385     /// The parameters of the algorithm can be specified using these
    386     /// functions.
    387 
    388     /// @{
    389 
    390     /// \brief Set the lower bounds on the arcs.
    391     ///
    392     /// This function sets the lower bounds on the arcs.
    393     /// If it is not used before calling \ref run(), the lower bounds
    394     /// will be set to zero on all arcs.
    395     ///
    396     /// \param map An arc map storing the lower bounds.
    397     /// Its \c Value type must be convertible to the \c Value type
    398     /// of the algorithm.
    399     ///
    400     /// \return <tt>(*this)</tt>
    401     template <typename LowerMap>
    402     CapacityScaling& lowerMap(const LowerMap& map) {
    403       _have_lower = true;
    404       for (ArcIt a(_graph); a != INVALID; ++a) {
    405         _lower[_arc_idf[a]] = map[a];
    406         _lower[_arc_idb[a]] = map[a];
    407       }
    408       return *this;
    409     }
    410 
    411     /// \brief Set the upper bounds (capacities) on the arcs.
    412     ///
    413     /// This function sets the upper bounds (capacities) on the arcs.
    414     /// If it is not used before calling \ref run(), the upper bounds
    415     /// will be set to \ref INF on all arcs (i.e. the flow value will be
    416     /// unbounded from above).
    417     ///
    418     /// \param map An arc map storing the upper bounds.
    419     /// Its \c Value type must be convertible to the \c Value type
    420     /// of the algorithm.
    421     ///
    422     /// \return <tt>(*this)</tt>
    423     template<typename UpperMap>
    424     CapacityScaling& upperMap(const UpperMap& map) {
    425       for (ArcIt a(_graph); a != INVALID; ++a) {
    426         _upper[_arc_idf[a]] = map[a];
    427       }
    428       return *this;
    429     }
    430 
    431     /// \brief Set the costs of the arcs.
    432     ///
    433     /// This function sets the costs of the arcs.
    434     /// If it is not used before calling \ref run(), the costs
    435     /// will be set to \c 1 on all arcs.
    436     ///
    437     /// \param map An arc map storing the costs.
    438     /// Its \c Value type must be convertible to the \c Cost type
    439     /// of the algorithm.
    440     ///
    441     /// \return <tt>(*this)</tt>
    442     template<typename CostMap>
    443     CapacityScaling& costMap(const CostMap& map) {
    444       for (ArcIt a(_graph); a != INVALID; ++a) {
    445         _cost[_arc_idf[a]] =  map[a];
    446         _cost[_arc_idb[a]] = -map[a];
    447       }
    448       return *this;
    449     }
    450 
    451     /// \brief Set the supply values of the nodes.
    452     ///
    453     /// This function sets the supply values of the nodes.
    454     /// If neither this function nor \ref stSupply() is used before
    455     /// calling \ref run(), the supply of each node will be set to zero.
    456     ///
    457     /// \param map A node map storing the supply values.
    458     /// Its \c Value type must be convertible to the \c Value type
    459     /// of the algorithm.
    460     ///
    461     /// \return <tt>(*this)</tt>
    462     template<typename SupplyMap>
    463     CapacityScaling& supplyMap(const SupplyMap& map) {
    464       for (NodeIt n(_graph); n != INVALID; ++n) {
    465         _supply[_node_id[n]] = map[n];
    466       }
    467       return *this;
    468     }
    469 
    470     /// \brief Set single source and target nodes and a supply value.
    471     ///
    472     /// This function sets a single source node and a single target node
    473     /// and the required flow value.
    474     /// If neither this function nor \ref supplyMap() is used before
    475     /// calling \ref run(), the supply of each node will be set to zero.
    476     ///
    477     /// Using this function has the same effect as using \ref supplyMap()
    478     /// with such a map in which \c k is assigned to \c s, \c -k is
    479     /// assigned to \c t and all other nodes have zero supply value.
    480     ///
    481     /// \param s The source node.
    482     /// \param t The target node.
    483     /// \param k The required amount of flow from node \c s to node \c t
    484     /// (i.e. the supply of \c s and the demand of \c t).
    485     ///
    486     /// \return <tt>(*this)</tt>
    487     CapacityScaling& stSupply(const Node& s, const Node& t, Value k) {
    488       for (int i = 0; i != _node_num; ++i) {
    489         _supply[i] = 0;
    490       }
    491       _supply[_node_id[s]] =  k;
    492       _supply[_node_id[t]] = -k;
    493       return *this;
    494     }
    495    
    496     /// @}
    497 
    498     /// \name Execution control
    499     /// The algorithm can be executed using \ref run().
    500 
    501     /// @{
    502 
    503     /// \brief Run the algorithm.
    504     ///
    505     /// This function runs the algorithm.
    506     /// The paramters can be specified using functions \ref lowerMap(),
    507     /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply().
    508     /// For example,
    509     /// \code
    510     ///   CapacityScaling<ListDigraph> cs(graph);
    511     ///   cs.lowerMap(lower).upperMap(upper).costMap(cost)
    512     ///     .supplyMap(sup).run();
    513     /// \endcode
    514     ///
    515     /// This function can be called more than once. All the parameters
    516     /// that have been given are kept for the next call, unless
    517     /// \ref reset() is called, thus only the modified parameters
    518     /// have to be set again. See \ref reset() for examples.
    519     /// However, the underlying digraph must not be modified after this
    520     /// class have been constructed, since it copies and extends the graph.
    521     ///
    522     /// \param factor The capacity scaling factor. It must be larger than
    523     /// one to use scaling. If it is less or equal to one, then scaling
    524     /// will be disabled.
    525     ///
    526     /// \return \c INFEASIBLE if no feasible flow exists,
    527     /// \n \c OPTIMAL if the problem has optimal solution
    528     /// (i.e. it is feasible and bounded), and the algorithm has found
    529     /// optimal flow and node potentials (primal and dual solutions),
    530     /// \n \c UNBOUNDED if the digraph contains an arc of negative cost
    531     /// and infinite upper bound. It means that the objective function
    532     /// is unbounded on that arc, however, note that it could actually be
    533     /// bounded over the feasible flows, but this algroithm cannot handle
    534     /// these cases.
    535     ///
    536     /// \see ProblemType
    537     ProblemType run(int factor = 4) {
    538       _factor = factor;
    539       ProblemType pt = init();
    540       if (pt != OPTIMAL) return pt;
    541       return start();
    542     }
    543 
    544     /// \brief Reset all the parameters that have been given before.
    545     ///
    546     /// This function resets all the paramaters that have been given
    547     /// before using functions \ref lowerMap(), \ref upperMap(),
    548     /// \ref costMap(), \ref supplyMap(), \ref stSupply().
    549     ///
    550     /// It is useful for multiple run() calls. If this function is not
    551     /// used, all the parameters given before are kept for the next
    552     /// \ref run() call.
    553     /// However, the underlying digraph must not be modified after this
    554     /// class have been constructed, since it copies and extends the graph.
    555     ///
    556     /// For example,
    557     /// \code
    558     ///   CapacityScaling<ListDigraph> cs(graph);
    559     ///
    560     ///   // First run
    561     ///   cs.lowerMap(lower).upperMap(upper).costMap(cost)
    562     ///     .supplyMap(sup).run();
    563     ///
    564     ///   // Run again with modified cost map (reset() is not called,
    565     ///   // so only the cost map have to be set again)
    566     ///   cost[e] += 100;
    567     ///   cs.costMap(cost).run();
    568     ///
    569     ///   // Run again from scratch using reset()
    570     ///   // (the lower bounds will be set to zero on all arcs)
    571     ///   cs.reset();
    572     ///   cs.upperMap(capacity).costMap(cost)
    573     ///     .supplyMap(sup).run();
    574     /// \endcode
    575     ///
    576     /// \return <tt>(*this)</tt>
    577     CapacityScaling& reset() {
    578       for (int i = 0; i != _node_num; ++i) {
    579         _supply[i] = 0;
    580       }
    581       for (int j = 0; j != _res_arc_num; ++j) {
    582         _lower[j] = 0;
    583         _upper[j] = INF;
    584         _cost[j] = _forward[j] ? 1 : -1;
    585       }
    586       _have_lower = false;
     620      resetParams();
    587621      return *this;
    588622    }
  • 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/cost_scaling.h

    r839 r840  
    105105  /// \tparam GR The digraph type the algorithm runs on.
    106106  /// \tparam V The number type used for flow amounts, capacity bounds
    107   /// and supply values in the algorithm. By default it is \c int.
     107  /// and supply values in the algorithm. By default, it is \c int.
    108108  /// \tparam C The number type used for costs and potentials in the
    109   /// algorithm. By default it is the same as \c V.
     109  /// algorithm. By default, it is the same as \c V.
     110  /// \tparam TR The traits class that defines various types used by the
     111  /// algorithm. By default, it is \ref CostScalingDefaultTraits
     112  /// "CostScalingDefaultTraits<GR, V, C>".
     113  /// In most cases, this parameter should not be set directly,
     114  /// consider to use the named template parameters instead.
    110115  ///
    111116  /// \warning Both number types must be signed and all input data must
     
    137142    ///
    138143    /// The large cost type used for internal computations.
    139     /// Using the \ref CostScalingDefaultTraits "default traits class",
    140     /// it is \c long \c long if the \c Cost type is integer,
     144    /// By default, it is \c long \c long if the \c Cost type is integer,
    141145    /// otherwise it is \c double.
    142146    typedef typename TR::LargeCost LargeCost;
     
    341345      LEMON_ASSERT(std::numeric_limits<Cost>::is_signed,
    342346        "The cost type of CostScaling must be signed");
    343 
     347     
     348      // Reset data structures
     349      reset();
     350    }
     351
     352    /// \name Parameters
     353    /// The parameters of the algorithm can be specified using these
     354    /// functions.
     355
     356    /// @{
     357
     358    /// \brief Set the lower bounds on the arcs.
     359    ///
     360    /// This function sets the lower bounds on the arcs.
     361    /// If it is not used before calling \ref run(), the lower bounds
     362    /// will be set to zero on all arcs.
     363    ///
     364    /// \param map An arc map storing the lower bounds.
     365    /// Its \c Value type must be convertible to the \c Value type
     366    /// of the algorithm.
     367    ///
     368    /// \return <tt>(*this)</tt>
     369    template <typename LowerMap>
     370    CostScaling& lowerMap(const LowerMap& map) {
     371      _have_lower = true;
     372      for (ArcIt a(_graph); a != INVALID; ++a) {
     373        _lower[_arc_idf[a]] = map[a];
     374        _lower[_arc_idb[a]] = map[a];
     375      }
     376      return *this;
     377    }
     378
     379    /// \brief Set the upper bounds (capacities) on the arcs.
     380    ///
     381    /// This function sets the upper bounds (capacities) on the arcs.
     382    /// If it is not used before calling \ref run(), the upper bounds
     383    /// will be set to \ref INF on all arcs (i.e. the flow value will be
     384    /// unbounded from above).
     385    ///
     386    /// \param map An arc map storing the upper bounds.
     387    /// Its \c Value type must be convertible to the \c Value type
     388    /// of the algorithm.
     389    ///
     390    /// \return <tt>(*this)</tt>
     391    template<typename UpperMap>
     392    CostScaling& upperMap(const UpperMap& map) {
     393      for (ArcIt a(_graph); a != INVALID; ++a) {
     394        _upper[_arc_idf[a]] = map[a];
     395      }
     396      return *this;
     397    }
     398
     399    /// \brief Set the costs of the arcs.
     400    ///
     401    /// This function sets the costs of the arcs.
     402    /// If it is not used before calling \ref run(), the costs
     403    /// will be set to \c 1 on all arcs.
     404    ///
     405    /// \param map An arc map storing the costs.
     406    /// Its \c Value type must be convertible to the \c Cost type
     407    /// of the algorithm.
     408    ///
     409    /// \return <tt>(*this)</tt>
     410    template<typename CostMap>
     411    CostScaling& costMap(const CostMap& map) {
     412      for (ArcIt a(_graph); a != INVALID; ++a) {
     413        _scost[_arc_idf[a]] =  map[a];
     414        _scost[_arc_idb[a]] = -map[a];
     415      }
     416      return *this;
     417    }
     418
     419    /// \brief Set the supply values of the nodes.
     420    ///
     421    /// This function sets the supply values of the nodes.
     422    /// If neither this function nor \ref stSupply() is used before
     423    /// calling \ref run(), the supply of each node will be set to zero.
     424    ///
     425    /// \param map A node map storing the supply values.
     426    /// Its \c Value type must be convertible to the \c Value type
     427    /// of the algorithm.
     428    ///
     429    /// \return <tt>(*this)</tt>
     430    template<typename SupplyMap>
     431    CostScaling& supplyMap(const SupplyMap& map) {
     432      for (NodeIt n(_graph); n != INVALID; ++n) {
     433        _supply[_node_id[n]] = map[n];
     434      }
     435      return *this;
     436    }
     437
     438    /// \brief Set single source and target nodes and a supply value.
     439    ///
     440    /// This function sets a single source node and a single target node
     441    /// and the required flow value.
     442    /// If neither this function nor \ref supplyMap() is used before
     443    /// calling \ref run(), the supply of each node will be set to zero.
     444    ///
     445    /// Using this function has the same effect as using \ref supplyMap()
     446    /// with such a map in which \c k is assigned to \c s, \c -k is
     447    /// assigned to \c t and all other nodes have zero supply value.
     448    ///
     449    /// \param s The source node.
     450    /// \param t The target node.
     451    /// \param k The required amount of flow from node \c s to node \c t
     452    /// (i.e. the supply of \c s and the demand of \c t).
     453    ///
     454    /// \return <tt>(*this)</tt>
     455    CostScaling& stSupply(const Node& s, const Node& t, Value k) {
     456      for (int i = 0; i != _res_node_num; ++i) {
     457        _supply[i] = 0;
     458      }
     459      _supply[_node_id[s]] =  k;
     460      _supply[_node_id[t]] = -k;
     461      return *this;
     462    }
     463   
     464    /// @}
     465
     466    /// \name Execution control
     467    /// The algorithm can be executed using \ref run().
     468
     469    /// @{
     470
     471    /// \brief Run the algorithm.
     472    ///
     473    /// This function runs the algorithm.
     474    /// The paramters can be specified using functions \ref lowerMap(),
     475    /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply().
     476    /// For example,
     477    /// \code
     478    ///   CostScaling<ListDigraph> cs(graph);
     479    ///   cs.lowerMap(lower).upperMap(upper).costMap(cost)
     480    ///     .supplyMap(sup).run();
     481    /// \endcode
     482    ///
     483    /// This function can be called more than once. All the given parameters
     484    /// are kept for the next call, unless \ref resetParams() or \ref reset()
     485    /// is used, thus only the modified parameters have to be set again.
     486    /// If the underlying digraph was also modified after the construction
     487    /// of the class (or the last \ref reset() call), then the \ref reset()
     488    /// function must be called.
     489    ///
     490    /// \param method The internal method that will be used in the
     491    /// algorithm. For more information, see \ref Method.
     492    /// \param factor The cost scaling factor. It must be larger than one.
     493    ///
     494    /// \return \c INFEASIBLE if no feasible flow exists,
     495    /// \n \c OPTIMAL if the problem has optimal solution
     496    /// (i.e. it is feasible and bounded), and the algorithm has found
     497    /// optimal flow and node potentials (primal and dual solutions),
     498    /// \n \c UNBOUNDED if the digraph contains an arc of negative cost
     499    /// and infinite upper bound. It means that the objective function
     500    /// is unbounded on that arc, however, note that it could actually be
     501    /// bounded over the feasible flows, but this algroithm cannot handle
     502    /// these cases.
     503    ///
     504    /// \see ProblemType, Method
     505    /// \see resetParams(), reset()
     506    ProblemType run(Method method = PARTIAL_AUGMENT, int factor = 8) {
     507      _alpha = factor;
     508      ProblemType pt = init();
     509      if (pt != OPTIMAL) return pt;
     510      start(method);
     511      return OPTIMAL;
     512    }
     513
     514    /// \brief Reset all the parameters that have been given before.
     515    ///
     516    /// This function resets all the paramaters that have been given
     517    /// before using functions \ref lowerMap(), \ref upperMap(),
     518    /// \ref costMap(), \ref supplyMap(), \ref stSupply().
     519    ///
     520    /// It is useful for multiple \ref run() calls. Basically, all the given
     521    /// parameters are kept for the next \ref run() call, unless
     522    /// \ref resetParams() or \ref reset() is used.
     523    /// If the underlying digraph was also modified after the construction
     524    /// of the class or the last \ref reset() call, then the \ref reset()
     525    /// function must be used, otherwise \ref resetParams() is sufficient.
     526    ///
     527    /// For example,
     528    /// \code
     529    ///   CostScaling<ListDigraph> cs(graph);
     530    ///
     531    ///   // First run
     532    ///   cs.lowerMap(lower).upperMap(upper).costMap(cost)
     533    ///     .supplyMap(sup).run();
     534    ///
     535    ///   // Run again with modified cost map (resetParams() is not called,
     536    ///   // so only the cost map have to be set again)
     537    ///   cost[e] += 100;
     538    ///   cs.costMap(cost).run();
     539    ///
     540    ///   // Run again from scratch using resetParams()
     541    ///   // (the lower bounds will be set to zero on all arcs)
     542    ///   cs.resetParams();
     543    ///   cs.upperMap(capacity).costMap(cost)
     544    ///     .supplyMap(sup).run();
     545    /// \endcode
     546    ///
     547    /// \return <tt>(*this)</tt>
     548    ///
     549    /// \see reset(), run()
     550    CostScaling& resetParams() {
     551      for (int i = 0; i != _res_node_num; ++i) {
     552        _supply[i] = 0;
     553      }
     554      int limit = _first_out[_root];
     555      for (int j = 0; j != limit; ++j) {
     556        _lower[j] = 0;
     557        _upper[j] = INF;
     558        _scost[j] = _forward[j] ? 1 : -1;
     559      }
     560      for (int j = limit; j != _res_arc_num; ++j) {
     561        _lower[j] = 0;
     562        _upper[j] = INF;
     563        _scost[j] = 0;
     564        _scost[_reverse[j]] = 0;
     565      }     
     566      _have_lower = false;
     567      return *this;
     568    }
     569
     570    /// \brief Reset all the parameters that have been given before.
     571    ///
     572    /// This function resets all the paramaters that have been given
     573    /// before using functions \ref lowerMap(), \ref upperMap(),
     574    /// \ref costMap(), \ref supplyMap(), \ref stSupply().
     575    ///
     576    /// It is useful for multiple run() calls. If this function is not
     577    /// used, all the parameters given before are kept for the next
     578    /// \ref run() call.
     579    /// However, the underlying digraph must not be modified after this
     580    /// class have been constructed, since it copies and extends the graph.
     581    /// \return <tt>(*this)</tt>
     582    CostScaling& reset() {
    344583      // Resize vectors
    345584      _node_num = countNodes(_graph);
     
    409648     
    410649      // Reset parameters
    411       reset();
    412     }
    413 
    414     /// \name Parameters
    415     /// The parameters of the algorithm can be specified using these
    416     /// functions.
    417 
    418     /// @{
    419 
    420     /// \brief Set the lower bounds on the arcs.
    421     ///
    422     /// This function sets the lower bounds on the arcs.
    423     /// If it is not used before calling \ref run(), the lower bounds
    424     /// will be set to zero on all arcs.
    425     ///
    426     /// \param map An arc map storing the lower bounds.
    427     /// Its \c Value type must be convertible to the \c Value type
    428     /// of the algorithm.
    429     ///
    430     /// \return <tt>(*this)</tt>
    431     template <typename LowerMap>
    432     CostScaling& lowerMap(const LowerMap& map) {
    433       _have_lower = true;
    434       for (ArcIt a(_graph); a != INVALID; ++a) {
    435         _lower[_arc_idf[a]] = map[a];
    436         _lower[_arc_idb[a]] = map[a];
    437       }
    438       return *this;
    439     }
    440 
    441     /// \brief Set the upper bounds (capacities) on the arcs.
    442     ///
    443     /// This function sets the upper bounds (capacities) on the arcs.
    444     /// If it is not used before calling \ref run(), the upper bounds
    445     /// will be set to \ref INF on all arcs (i.e. the flow value will be
    446     /// unbounded from above).
    447     ///
    448     /// \param map An arc map storing the upper bounds.
    449     /// Its \c Value type must be convertible to the \c Value type
    450     /// of the algorithm.
    451     ///
    452     /// \return <tt>(*this)</tt>
    453     template<typename UpperMap>
    454     CostScaling& upperMap(const UpperMap& map) {
    455       for (ArcIt a(_graph); a != INVALID; ++a) {
    456         _upper[_arc_idf[a]] = map[a];
    457       }
    458       return *this;
    459     }
    460 
    461     /// \brief Set the costs of the arcs.
    462     ///
    463     /// This function sets the costs of the arcs.
    464     /// If it is not used before calling \ref run(), the costs
    465     /// will be set to \c 1 on all arcs.
    466     ///
    467     /// \param map An arc map storing the costs.
    468     /// Its \c Value type must be convertible to the \c Cost type
    469     /// of the algorithm.
    470     ///
    471     /// \return <tt>(*this)</tt>
    472     template<typename CostMap>
    473     CostScaling& costMap(const CostMap& map) {
    474       for (ArcIt a(_graph); a != INVALID; ++a) {
    475         _scost[_arc_idf[a]] =  map[a];
    476         _scost[_arc_idb[a]] = -map[a];
    477       }
    478       return *this;
    479     }
    480 
    481     /// \brief Set the supply values of the nodes.
    482     ///
    483     /// This function sets the supply values of the nodes.
    484     /// If neither this function nor \ref stSupply() is used before
    485     /// calling \ref run(), the supply of each node will be set to zero.
    486     ///
    487     /// \param map A node map storing the supply values.
    488     /// Its \c Value type must be convertible to the \c Value type
    489     /// of the algorithm.
    490     ///
    491     /// \return <tt>(*this)</tt>
    492     template<typename SupplyMap>
    493     CostScaling& supplyMap(const SupplyMap& map) {
    494       for (NodeIt n(_graph); n != INVALID; ++n) {
    495         _supply[_node_id[n]] = map[n];
    496       }
    497       return *this;
    498     }
    499 
    500     /// \brief Set single source and target nodes and a supply value.
    501     ///
    502     /// This function sets a single source node and a single target node
    503     /// and the required flow value.
    504     /// If neither this function nor \ref supplyMap() is used before
    505     /// calling \ref run(), the supply of each node will be set to zero.
    506     ///
    507     /// Using this function has the same effect as using \ref supplyMap()
    508     /// with such a map in which \c k is assigned to \c s, \c -k is
    509     /// assigned to \c t and all other nodes have zero supply value.
    510     ///
    511     /// \param s The source node.
    512     /// \param t The target node.
    513     /// \param k The required amount of flow from node \c s to node \c t
    514     /// (i.e. the supply of \c s and the demand of \c t).
    515     ///
    516     /// \return <tt>(*this)</tt>
    517     CostScaling& stSupply(const Node& s, const Node& t, Value k) {
    518       for (int i = 0; i != _res_node_num; ++i) {
    519         _supply[i] = 0;
    520       }
    521       _supply[_node_id[s]] =  k;
    522       _supply[_node_id[t]] = -k;
    523       return *this;
    524     }
    525    
    526     /// @}
    527 
    528     /// \name Execution control
    529     /// The algorithm can be executed using \ref run().
    530 
    531     /// @{
    532 
    533     /// \brief Run the algorithm.
    534     ///
    535     /// This function runs the algorithm.
    536     /// The paramters can be specified using functions \ref lowerMap(),
    537     /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply().
    538     /// For example,
    539     /// \code
    540     ///   CostScaling<ListDigraph> cs(graph);
    541     ///   cs.lowerMap(lower).upperMap(upper).costMap(cost)
    542     ///     .supplyMap(sup).run();
    543     /// \endcode
    544     ///
    545     /// This function can be called more than once. All the parameters
    546     /// that have been given are kept for the next call, unless
    547     /// \ref reset() is called, thus only the modified parameters
    548     /// have to be set again. See \ref reset() for examples.
    549     /// However, the underlying digraph must not be modified after this
    550     /// class have been constructed, since it copies and extends the graph.
    551     ///
    552     /// \param method The internal method that will be used in the
    553     /// algorithm. For more information, see \ref Method.
    554     /// \param factor The cost scaling factor. It must be larger than one.
    555     ///
    556     /// \return \c INFEASIBLE if no feasible flow exists,
    557     /// \n \c OPTIMAL if the problem has optimal solution
    558     /// (i.e. it is feasible and bounded), and the algorithm has found
    559     /// optimal flow and node potentials (primal and dual solutions),
    560     /// \n \c UNBOUNDED if the digraph contains an arc of negative cost
    561     /// and infinite upper bound. It means that the objective function
    562     /// is unbounded on that arc, however, note that it could actually be
    563     /// bounded over the feasible flows, but this algroithm cannot handle
    564     /// these cases.
    565     ///
    566     /// \see ProblemType, Method
    567     ProblemType run(Method method = PARTIAL_AUGMENT, int factor = 8) {
    568       _alpha = factor;
    569       ProblemType pt = init();
    570       if (pt != OPTIMAL) return pt;
    571       start(method);
    572       return OPTIMAL;
    573     }
    574 
    575     /// \brief Reset all the parameters that have been given before.
    576     ///
    577     /// This function resets all the paramaters that have been given
    578     /// before using functions \ref lowerMap(), \ref upperMap(),
    579     /// \ref costMap(), \ref supplyMap(), \ref stSupply().
    580     ///
    581     /// It is useful for multiple run() calls. If this function is not
    582     /// used, all the parameters given before are kept for the next
    583     /// \ref run() call.
    584     /// However, the underlying digraph must not be modified after this
    585     /// class have been constructed, since it copies and extends the graph.
    586     ///
    587     /// For example,
    588     /// \code
    589     ///   CostScaling<ListDigraph> cs(graph);
    590     ///
    591     ///   // First run
    592     ///   cs.lowerMap(lower).upperMap(upper).costMap(cost)
    593     ///     .supplyMap(sup).run();
    594     ///
    595     ///   // Run again with modified cost map (reset() is not called,
    596     ///   // so only the cost map have to be set again)
    597     ///   cost[e] += 100;
    598     ///   cs.costMap(cost).run();
    599     ///
    600     ///   // Run again from scratch using reset()
    601     ///   // (the lower bounds will be set to zero on all arcs)
    602     ///   cs.reset();
    603     ///   cs.upperMap(capacity).costMap(cost)
    604     ///     .supplyMap(sup).run();
    605     /// \endcode
    606     ///
    607     /// \return <tt>(*this)</tt>
    608     CostScaling& reset() {
    609       for (int i = 0; i != _res_node_num; ++i) {
    610         _supply[i] = 0;
    611       }
    612       int limit = _first_out[_root];
    613       for (int j = 0; j != limit; ++j) {
    614         _lower[j] = 0;
    615         _upper[j] = INF;
    616         _scost[j] = _forward[j] ? 1 : -1;
    617       }
    618       for (int j = limit; j != _res_arc_num; ++j) {
    619         _lower[j] = 0;
    620         _upper[j] = INF;
    621         _scost[j] = 0;
    622         _scost[_reverse[j]] = 0;
    623       }     
    624       _have_lower = false;
     650      resetParams();
    625651      return *this;
    626652    }
  • 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/cycle_canceling.h

    r839 r840  
    252252        "The cost type of CycleCanceling must be signed");
    253253
     254      // Reset data structures
     255      reset();
     256    }
     257
     258    /// \name Parameters
     259    /// The parameters of the algorithm can be specified using these
     260    /// functions.
     261
     262    /// @{
     263
     264    /// \brief Set the lower bounds on the arcs.
     265    ///
     266    /// This function sets the lower bounds on the arcs.
     267    /// If it is not used before calling \ref run(), the lower bounds
     268    /// will be set to zero on all arcs.
     269    ///
     270    /// \param map An arc map storing the lower bounds.
     271    /// Its \c Value type must be convertible to the \c Value type
     272    /// of the algorithm.
     273    ///
     274    /// \return <tt>(*this)</tt>
     275    template <typename LowerMap>
     276    CycleCanceling& lowerMap(const LowerMap& map) {
     277      _have_lower = true;
     278      for (ArcIt a(_graph); a != INVALID; ++a) {
     279        _lower[_arc_idf[a]] = map[a];
     280        _lower[_arc_idb[a]] = map[a];
     281      }
     282      return *this;
     283    }
     284
     285    /// \brief Set the upper bounds (capacities) on the arcs.
     286    ///
     287    /// This function sets the upper bounds (capacities) on the arcs.
     288    /// If it is not used before calling \ref run(), the upper bounds
     289    /// will be set to \ref INF on all arcs (i.e. the flow value will be
     290    /// unbounded from above).
     291    ///
     292    /// \param map An arc map storing the upper bounds.
     293    /// Its \c Value type must be convertible to the \c Value type
     294    /// of the algorithm.
     295    ///
     296    /// \return <tt>(*this)</tt>
     297    template<typename UpperMap>
     298    CycleCanceling& upperMap(const UpperMap& map) {
     299      for (ArcIt a(_graph); a != INVALID; ++a) {
     300        _upper[_arc_idf[a]] = map[a];
     301      }
     302      return *this;
     303    }
     304
     305    /// \brief Set the costs of the arcs.
     306    ///
     307    /// This function sets the costs of the arcs.
     308    /// If it is not used before calling \ref run(), the costs
     309    /// will be set to \c 1 on all arcs.
     310    ///
     311    /// \param map An arc map storing the costs.
     312    /// Its \c Value type must be convertible to the \c Cost type
     313    /// of the algorithm.
     314    ///
     315    /// \return <tt>(*this)</tt>
     316    template<typename CostMap>
     317    CycleCanceling& costMap(const CostMap& map) {
     318      for (ArcIt a(_graph); a != INVALID; ++a) {
     319        _cost[_arc_idf[a]] =  map[a];
     320        _cost[_arc_idb[a]] = -map[a];
     321      }
     322      return *this;
     323    }
     324
     325    /// \brief Set the supply values of the nodes.
     326    ///
     327    /// This function sets the supply values of the nodes.
     328    /// If neither this function nor \ref stSupply() is used before
     329    /// calling \ref run(), the supply of each node will be set to zero.
     330    ///
     331    /// \param map A node map storing the supply values.
     332    /// Its \c Value type must be convertible to the \c Value type
     333    /// of the algorithm.
     334    ///
     335    /// \return <tt>(*this)</tt>
     336    template<typename SupplyMap>
     337    CycleCanceling& supplyMap(const SupplyMap& map) {
     338      for (NodeIt n(_graph); n != INVALID; ++n) {
     339        _supply[_node_id[n]] = map[n];
     340      }
     341      return *this;
     342    }
     343
     344    /// \brief Set single source and target nodes and a supply value.
     345    ///
     346    /// This function sets a single source node and a single target node
     347    /// and the required flow value.
     348    /// If neither this function nor \ref supplyMap() is used before
     349    /// calling \ref run(), the supply of each node will be set to zero.
     350    ///
     351    /// Using this function has the same effect as using \ref supplyMap()
     352    /// with such a map in which \c k is assigned to \c s, \c -k is
     353    /// assigned to \c t and all other nodes have zero supply value.
     354    ///
     355    /// \param s The source node.
     356    /// \param t The target node.
     357    /// \param k The required amount of flow from node \c s to node \c t
     358    /// (i.e. the supply of \c s and the demand of \c t).
     359    ///
     360    /// \return <tt>(*this)</tt>
     361    CycleCanceling& stSupply(const Node& s, const Node& t, Value k) {
     362      for (int i = 0; i != _res_node_num; ++i) {
     363        _supply[i] = 0;
     364      }
     365      _supply[_node_id[s]] =  k;
     366      _supply[_node_id[t]] = -k;
     367      return *this;
     368    }
     369   
     370    /// @}
     371
     372    /// \name Execution control
     373    /// The algorithm can be executed using \ref run().
     374
     375    /// @{
     376
     377    /// \brief Run the algorithm.
     378    ///
     379    /// This function runs the algorithm.
     380    /// The paramters can be specified using functions \ref lowerMap(),
     381    /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply().
     382    /// For example,
     383    /// \code
     384    ///   CycleCanceling<ListDigraph> cc(graph);
     385    ///   cc.lowerMap(lower).upperMap(upper).costMap(cost)
     386    ///     .supplyMap(sup).run();
     387    /// \endcode
     388    ///
     389    /// This function can be called more than once. All the given parameters
     390    /// are kept for the next call, unless \ref resetParams() or \ref reset()
     391    /// is used, thus only the modified parameters have to be set again.
     392    /// If the underlying digraph was also modified after the construction
     393    /// of the class (or the last \ref reset() call), then the \ref reset()
     394    /// function must be called.
     395    ///
     396    /// \param method The cycle-canceling method that will be used.
     397    /// For more information, see \ref Method.
     398    ///
     399    /// \return \c INFEASIBLE if no feasible flow exists,
     400    /// \n \c OPTIMAL if the problem has optimal solution
     401    /// (i.e. it is feasible and bounded), and the algorithm has found
     402    /// optimal flow and node potentials (primal and dual solutions),
     403    /// \n \c UNBOUNDED if the digraph contains an arc of negative cost
     404    /// and infinite upper bound. It means that the objective function
     405    /// is unbounded on that arc, however, note that it could actually be
     406    /// bounded over the feasible flows, but this algroithm cannot handle
     407    /// these cases.
     408    ///
     409    /// \see ProblemType, Method
     410    /// \see resetParams(), reset()
     411    ProblemType run(Method method = CANCEL_AND_TIGHTEN) {
     412      ProblemType pt = init();
     413      if (pt != OPTIMAL) return pt;
     414      start(method);
     415      return OPTIMAL;
     416    }
     417
     418    /// \brief Reset all the parameters that have been given before.
     419    ///
     420    /// This function resets all the paramaters that have been given
     421    /// before using functions \ref lowerMap(), \ref upperMap(),
     422    /// \ref costMap(), \ref supplyMap(), \ref stSupply().
     423    ///
     424    /// It is useful for multiple \ref run() calls. Basically, all the given
     425    /// parameters are kept for the next \ref run() call, unless
     426    /// \ref resetParams() or \ref reset() is used.
     427    /// If the underlying digraph was also modified after the construction
     428    /// of the class or the last \ref reset() call, then the \ref reset()
     429    /// function must be used, otherwise \ref resetParams() is sufficient.
     430    ///
     431    /// For example,
     432    /// \code
     433    ///   CycleCanceling<ListDigraph> cs(graph);
     434    ///
     435    ///   // First run
     436    ///   cc.lowerMap(lower).upperMap(upper).costMap(cost)
     437    ///     .supplyMap(sup).run();
     438    ///
     439    ///   // Run again with modified cost map (resetParams() is not called,
     440    ///   // so only the cost map have to be set again)
     441    ///   cost[e] += 100;
     442    ///   cc.costMap(cost).run();
     443    ///
     444    ///   // Run again from scratch using resetParams()
     445    ///   // (the lower bounds will be set to zero on all arcs)
     446    ///   cc.resetParams();
     447    ///   cc.upperMap(capacity).costMap(cost)
     448    ///     .supplyMap(sup).run();
     449    /// \endcode
     450    ///
     451    /// \return <tt>(*this)</tt>
     452    ///
     453    /// \see reset(), run()
     454    CycleCanceling& resetParams() {
     455      for (int i = 0; i != _res_node_num; ++i) {
     456        _supply[i] = 0;
     457      }
     458      int limit = _first_out[_root];
     459      for (int j = 0; j != limit; ++j) {
     460        _lower[j] = 0;
     461        _upper[j] = INF;
     462        _cost[j] = _forward[j] ? 1 : -1;
     463      }
     464      for (int j = limit; j != _res_arc_num; ++j) {
     465        _lower[j] = 0;
     466        _upper[j] = INF;
     467        _cost[j] = 0;
     468        _cost[_reverse[j]] = 0;
     469      }     
     470      _have_lower = false;
     471      return *this;
     472    }
     473
     474    /// \brief Reset the internal data structures and all the parameters
     475    /// that have been given before.
     476    ///
     477    /// This function resets the internal data structures and all the
     478    /// paramaters that have been given before using functions \ref lowerMap(),
     479    /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply().
     480    ///
     481    /// It is useful for multiple \ref run() calls. Basically, all the given
     482    /// parameters are kept for the next \ref run() call, unless
     483    /// \ref resetParams() or \ref reset() is used.
     484    /// If the underlying digraph was also modified after the construction
     485    /// of the class or the last \ref reset() call, then the \ref reset()
     486    /// function must be used, otherwise \ref resetParams() is sufficient.
     487    ///
     488    /// See \ref resetParams() for examples.
     489    ///
     490    /// \return <tt>(*this)</tt>
     491    ///
     492    /// \see resetParams(), run()
     493    CycleCanceling& reset() {
    254494      // Resize vectors
    255495      _node_num = countNodes(_graph);
     
    317557     
    318558      // Reset parameters
    319       reset();
    320     }
    321 
    322     /// \name Parameters
    323     /// The parameters of the algorithm can be specified using these
    324     /// functions.
    325 
    326     /// @{
    327 
    328     /// \brief Set the lower bounds on the arcs.
    329     ///
    330     /// This function sets the lower bounds on the arcs.
    331     /// If it is not used before calling \ref run(), the lower bounds
    332     /// will be set to zero on all arcs.
    333     ///
    334     /// \param map An arc map storing the lower bounds.
    335     /// Its \c Value type must be convertible to the \c Value type
    336     /// of the algorithm.
    337     ///
    338     /// \return <tt>(*this)</tt>
    339     template <typename LowerMap>
    340     CycleCanceling& lowerMap(const LowerMap& map) {
    341       _have_lower = true;
    342       for (ArcIt a(_graph); a != INVALID; ++a) {
    343         _lower[_arc_idf[a]] = map[a];
    344         _lower[_arc_idb[a]] = map[a];
    345       }
    346       return *this;
    347     }
    348 
    349     /// \brief Set the upper bounds (capacities) on the arcs.
    350     ///
    351     /// This function sets the upper bounds (capacities) on the arcs.
    352     /// If it is not used before calling \ref run(), the upper bounds
    353     /// will be set to \ref INF on all arcs (i.e. the flow value will be
    354     /// unbounded from above).
    355     ///
    356     /// \param map An arc map storing the upper bounds.
    357     /// Its \c Value type must be convertible to the \c Value type
    358     /// of the algorithm.
    359     ///
    360     /// \return <tt>(*this)</tt>
    361     template<typename UpperMap>
    362     CycleCanceling& upperMap(const UpperMap& map) {
    363       for (ArcIt a(_graph); a != INVALID; ++a) {
    364         _upper[_arc_idf[a]] = map[a];
    365       }
    366       return *this;
    367     }
    368 
    369     /// \brief Set the costs of the arcs.
    370     ///
    371     /// This function sets the costs of the arcs.
    372     /// If it is not used before calling \ref run(), the costs
    373     /// will be set to \c 1 on all arcs.
    374     ///
    375     /// \param map An arc map storing the costs.
    376     /// Its \c Value type must be convertible to the \c Cost type
    377     /// of the algorithm.
    378     ///
    379     /// \return <tt>(*this)</tt>
    380     template<typename CostMap>
    381     CycleCanceling& costMap(const CostMap& map) {
    382       for (ArcIt a(_graph); a != INVALID; ++a) {
    383         _cost[_arc_idf[a]] =  map[a];
    384         _cost[_arc_idb[a]] = -map[a];
    385       }
    386       return *this;
    387     }
    388 
    389     /// \brief Set the supply values of the nodes.
    390     ///
    391     /// This function sets the supply values of the nodes.
    392     /// If neither this function nor \ref stSupply() is used before
    393     /// calling \ref run(), the supply of each node will be set to zero.
    394     ///
    395     /// \param map A node map storing the supply values.
    396     /// Its \c Value type must be convertible to the \c Value type
    397     /// of the algorithm.
    398     ///
    399     /// \return <tt>(*this)</tt>
    400     template<typename SupplyMap>
    401     CycleCanceling& supplyMap(const SupplyMap& map) {
    402       for (NodeIt n(_graph); n != INVALID; ++n) {
    403         _supply[_node_id[n]] = map[n];
    404       }
    405       return *this;
    406     }
    407 
    408     /// \brief Set single source and target nodes and a supply value.
    409     ///
    410     /// This function sets a single source node and a single target node
    411     /// and the required flow value.
    412     /// If neither this function nor \ref supplyMap() is used before
    413     /// calling \ref run(), the supply of each node will be set to zero.
    414     ///
    415     /// Using this function has the same effect as using \ref supplyMap()
    416     /// with such a map in which \c k is assigned to \c s, \c -k is
    417     /// assigned to \c t and all other nodes have zero supply value.
    418     ///
    419     /// \param s The source node.
    420     /// \param t The target node.
    421     /// \param k The required amount of flow from node \c s to node \c t
    422     /// (i.e. the supply of \c s and the demand of \c t).
    423     ///
    424     /// \return <tt>(*this)</tt>
    425     CycleCanceling& stSupply(const Node& s, const Node& t, Value k) {
    426       for (int i = 0; i != _res_node_num; ++i) {
    427         _supply[i] = 0;
    428       }
    429       _supply[_node_id[s]] =  k;
    430       _supply[_node_id[t]] = -k;
    431       return *this;
    432     }
    433    
    434     /// @}
    435 
    436     /// \name Execution control
    437     /// The algorithm can be executed using \ref run().
    438 
    439     /// @{
    440 
    441     /// \brief Run the algorithm.
    442     ///
    443     /// This function runs the algorithm.
    444     /// The paramters can be specified using functions \ref lowerMap(),
    445     /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply().
    446     /// For example,
    447     /// \code
    448     ///   CycleCanceling<ListDigraph> cc(graph);
    449     ///   cc.lowerMap(lower).upperMap(upper).costMap(cost)
    450     ///     .supplyMap(sup).run();
    451     /// \endcode
    452     ///
    453     /// This function can be called more than once. All the parameters
    454     /// that have been given are kept for the next call, unless
    455     /// \ref reset() is called, thus only the modified parameters
    456     /// have to be set again. See \ref reset() for examples.
    457     /// However, the underlying digraph must not be modified after this
    458     /// class have been constructed, since it copies and extends the graph.
    459     ///
    460     /// \param method The cycle-canceling method that will be used.
    461     /// For more information, see \ref Method.
    462     ///
    463     /// \return \c INFEASIBLE if no feasible flow exists,
    464     /// \n \c OPTIMAL if the problem has optimal solution
    465     /// (i.e. it is feasible and bounded), and the algorithm has found
    466     /// optimal flow and node potentials (primal and dual solutions),
    467     /// \n \c UNBOUNDED if the digraph contains an arc of negative cost
    468     /// and infinite upper bound. It means that the objective function
    469     /// is unbounded on that arc, however, note that it could actually be
    470     /// bounded over the feasible flows, but this algroithm cannot handle
    471     /// these cases.
    472     ///
    473     /// \see ProblemType, Method
    474     ProblemType run(Method method = CANCEL_AND_TIGHTEN) {
    475       ProblemType pt = init();
    476       if (pt != OPTIMAL) return pt;
    477       start(method);
    478       return OPTIMAL;
    479     }
    480 
    481     /// \brief Reset all the parameters that have been given before.
    482     ///
    483     /// This function resets all the paramaters that have been given
    484     /// before using functions \ref lowerMap(), \ref upperMap(),
    485     /// \ref costMap(), \ref supplyMap(), \ref stSupply().
    486     ///
    487     /// It is useful for multiple run() calls. If this function is not
    488     /// used, all the parameters given before are kept for the next
    489     /// \ref run() call.
    490     /// However, the underlying digraph must not be modified after this
    491     /// class have been constructed, since it copies and extends the graph.
    492     ///
    493     /// For example,
    494     /// \code
    495     ///   CycleCanceling<ListDigraph> cs(graph);
    496     ///
    497     ///   // First run
    498     ///   cc.lowerMap(lower).upperMap(upper).costMap(cost)
    499     ///     .supplyMap(sup).run();
    500     ///
    501     ///   // Run again with modified cost map (reset() is not called,
    502     ///   // so only the cost map have to be set again)
    503     ///   cost[e] += 100;
    504     ///   cc.costMap(cost).run();
    505     ///
    506     ///   // Run again from scratch using reset()
    507     ///   // (the lower bounds will be set to zero on all arcs)
    508     ///   cc.reset();
    509     ///   cc.upperMap(capacity).costMap(cost)
    510     ///     .supplyMap(sup).run();
    511     /// \endcode
    512     ///
    513     /// \return <tt>(*this)</tt>
    514     CycleCanceling& reset() {
    515       for (int i = 0; i != _res_node_num; ++i) {
    516         _supply[i] = 0;
    517       }
    518       int limit = _first_out[_root];
    519       for (int j = 0; j != limit; ++j) {
    520         _lower[j] = 0;
    521         _upper[j] = INF;
    522         _cost[j] = _forward[j] ? 1 : -1;
    523       }
    524       for (int j = limit; j != _res_arc_num; ++j) {
    525         _lower[j] = 0;
    526         _upper[j] = INF;
    527         _cost[j] = 0;
    528         _cost[_reverse[j]] = 0;
    529       }     
    530       _have_lower = false;
     559      resetParams();
    531560      return *this;
    532561    }
  • 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()) {
  • lemon/network_simplex.h

    r839 r840  
    196196    IntVector _source;
    197197    IntVector _target;
     198    bool _arc_mixing;
    198199
    199200    // Node and arc data
     
    635636    NetworkSimplex(const GR& graph, bool arc_mixing = false) :
    636637      _graph(graph), _node_id(graph), _arc_id(graph),
     638      _arc_mixing(arc_mixing),
    637639      MAX(std::numeric_limits<Value>::max()),
    638640      INF(std::numeric_limits<Value>::has_infinity ?
     
    645647        "The cost type of NetworkSimplex must be signed");
    646648       
    647       // Resize vectors
    648       _node_num = countNodes(_graph);
    649       _arc_num = countArcs(_graph);
    650       int all_node_num = _node_num + 1;
    651       int max_arc_num = _arc_num + 2 * _node_num;
    652 
    653       _source.resize(max_arc_num);
    654       _target.resize(max_arc_num);
    655 
    656       _lower.resize(_arc_num);
    657       _upper.resize(_arc_num);
    658       _cap.resize(max_arc_num);
    659       _cost.resize(max_arc_num);
    660       _supply.resize(all_node_num);
    661       _flow.resize(max_arc_num);
    662       _pi.resize(all_node_num);
    663 
    664       _parent.resize(all_node_num);
    665       _pred.resize(all_node_num);
    666       _forward.resize(all_node_num);
    667       _thread.resize(all_node_num);
    668       _rev_thread.resize(all_node_num);
    669       _succ_num.resize(all_node_num);
    670       _last_succ.resize(all_node_num);
    671       _state.resize(max_arc_num);
    672 
    673       // Copy the graph
    674       int i = 0;
    675       for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
    676         _node_id[n] = i;
    677       }
    678       if (arc_mixing) {
    679         // Store the arcs in a mixed order
    680         int k = std::max(int(std::sqrt(double(_arc_num))), 10);
    681         int i = 0, j = 0;
    682         for (ArcIt a(_graph); a != INVALID; ++a) {
    683           _arc_id[a] = i;
    684           _source[i] = _node_id[_graph.source(a)];
    685           _target[i] = _node_id[_graph.target(a)];
    686           if ((i += k) >= _arc_num) i = ++j;
    687         }
    688       } else {
    689         // Store the arcs in the original order
    690         int i = 0;
    691         for (ArcIt a(_graph); a != INVALID; ++a, ++i) {
    692           _arc_id[a] = i;
    693           _source[i] = _node_id[_graph.source(a)];
    694           _target[i] = _node_id[_graph.target(a)];
    695         }
    696       }
    697      
    698       // Reset parameters
     649      // Reset data structures
    699650      reset();
    700651    }
     
    844795    /// \endcode
    845796    ///
    846     /// This function can be called more than once. All the parameters
    847     /// that have been given are kept for the next call, unless
    848     /// \ref reset() is called, thus only the modified parameters
    849     /// have to be set again. See \ref reset() for examples.
    850     /// However, the underlying digraph must not be modified after this
    851     /// class have been constructed, since it copies and extends the graph.
     797    /// This function can be called more than once. All the given parameters
     798    /// are kept for the next call, unless \ref resetParams() or \ref reset()
     799    /// is used, thus only the modified parameters have to be set again.
     800    /// If the underlying digraph was also modified after the construction
     801    /// of the class (or the last \ref reset() call), then the \ref reset()
     802    /// function must be called.
    852803    ///
    853804    /// \param pivot_rule The pivot rule that will be used during the
     
    863814    ///
    864815    /// \see ProblemType, PivotRule
     816    /// \see resetParams(), reset()
    865817    ProblemType run(PivotRule pivot_rule = BLOCK_SEARCH) {
    866818      if (!init()) return INFEASIBLE;
     
    874826    /// \ref costMap(), \ref supplyMap(), \ref stSupply(), \ref supplyType().
    875827    ///
    876     /// It is useful for multiple run() calls. If this function is not
    877     /// used, all the parameters given before are kept for the next
    878     /// \ref run() call.
    879     /// However, the underlying digraph must not be modified after this
    880     /// class have been constructed, since it copies and extends the graph.
     828    /// It is useful for multiple \ref run() calls. Basically, all the given
     829    /// parameters are kept for the next \ref run() call, unless
     830    /// \ref resetParams() or \ref reset() is used.
     831    /// If the underlying digraph was also modified after the construction
     832    /// of the class or the last \ref reset() call, then the \ref reset()
     833    /// function must be used, otherwise \ref resetParams() is sufficient.
    881834    ///
    882835    /// For example,
     
    888841    ///     .supplyMap(sup).run();
    889842    ///
    890     ///   // Run again with modified cost map (reset() is not called,
     843    ///   // Run again with modified cost map (resetParams() is not called,
    891844    ///   // so only the cost map have to be set again)
    892845    ///   cost[e] += 100;
    893846    ///   ns.costMap(cost).run();
    894847    ///
    895     ///   // Run again from scratch using reset()
     848    ///   // Run again from scratch using resetParams()
    896849    ///   // (the lower bounds will be set to zero on all arcs)
    897     ///   ns.reset();
     850    ///   ns.resetParams();
    898851    ///   ns.upperMap(capacity).costMap(cost)
    899852    ///     .supplyMap(sup).run();
     
    901854    ///
    902855    /// \return <tt>(*this)</tt>
    903     NetworkSimplex& reset() {
     856    ///
     857    /// \see reset(), run()
     858    NetworkSimplex& resetParams() {
    904859      for (int i = 0; i != _node_num; ++i) {
    905860        _supply[i] = 0;
     
    915870    }
    916871
     872    /// \brief Reset the internal data structures and all the parameters
     873    /// that have been given before.
     874    ///
     875    /// This function resets the internal data structures and all the
     876    /// paramaters that have been given before using functions \ref lowerMap(),
     877    /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(),
     878    /// \ref supplyType().
     879    ///
     880    /// It is useful for multiple \ref run() calls. Basically, all the given
     881    /// parameters are kept for the next \ref run() call, unless
     882    /// \ref resetParams() or \ref reset() is used.
     883    /// If the underlying digraph was also modified after the construction
     884    /// of the class or the last \ref reset() call, then the \ref reset()
     885    /// function must be used, otherwise \ref resetParams() is sufficient.
     886    ///
     887    /// See \ref resetParams() for examples.
     888    ///
     889    /// \return <tt>(*this)</tt>
     890    ///
     891    /// \see resetParams(), run()
     892    NetworkSimplex& reset() {
     893      // Resize vectors
     894      _node_num = countNodes(_graph);
     895      _arc_num = countArcs(_graph);
     896      int all_node_num = _node_num + 1;
     897      int max_arc_num = _arc_num + 2 * _node_num;
     898
     899      _source.resize(max_arc_num);
     900      _target.resize(max_arc_num);
     901
     902      _lower.resize(_arc_num);
     903      _upper.resize(_arc_num);
     904      _cap.resize(max_arc_num);
     905      _cost.resize(max_arc_num);
     906      _supply.resize(all_node_num);
     907      _flow.resize(max_arc_num);
     908      _pi.resize(all_node_num);
     909
     910      _parent.resize(all_node_num);
     911      _pred.resize(all_node_num);
     912      _forward.resize(all_node_num);
     913      _thread.resize(all_node_num);
     914      _rev_thread.resize(all_node_num);
     915      _succ_num.resize(all_node_num);
     916      _last_succ.resize(all_node_num);
     917      _state.resize(max_arc_num);
     918
     919      // Copy the graph
     920      int i = 0;
     921      for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
     922        _node_id[n] = i;
     923      }
     924      if (_arc_mixing) {
     925        // Store the arcs in a mixed order
     926        int k = std::max(int(std::sqrt(double(_arc_num))), 10);
     927        int i = 0, j = 0;
     928        for (ArcIt a(_graph); a != INVALID; ++a) {
     929          _arc_id[a] = i;
     930          _source[i] = _node_id[_graph.source(a)];
     931          _target[i] = _node_id[_graph.target(a)];
     932          if ((i += k) >= _arc_num) i = ++j;
     933        }
     934      } else {
     935        // Store the arcs in the original order
     936        int i = 0;
     937        for (ArcIt a(_graph); a != INVALID; ++a, ++i) {
     938          _arc_id[a] = i;
     939          _source[i] = _node_id[_graph.source(a)];
     940          _target[i] = _node_id[_graph.target(a)];
     941        }
     942      }
     943     
     944      // Reset parameters
     945      resetParams();
     946      return *this;
     947    }
     948   
    917949    /// @}
    918950
Note: See TracChangeset for help on using the changeset viewer.