*

  *

  *

  *

  * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport

  * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport

  * (Egervary Research Group on Combinatorial Optimization, EGRES).

  * (Egervary Research Group on Combinatorial Optimization, EGRES).

  *

  *

  * Permission to use, modify and distribute this software is granted

  * Permission to use, modify and distribute this software is granted

  * provided that this copyright notice appears in all copies. For

  * provided that this copyright notice appears in all copies. For

     };

     };

   private:

   private:

     TEMPLATE_DIGRAPH_TYPEDEFS(GR);

     TEMPLATE_DIGRAPH_TYPEDEFS(GR);

     typedef std::vector<int> IntVector;

     typedef std::vector<int> IntVector;

     typedef std::vector<double> DoubleVector;

     typedef std::vector<double> DoubleVector;

     typedef std::vector<Value> ValueVector;

     typedef std::vector<Value> ValueVector;

     typedef std::vector<Cost> CostVector;

     typedef std::vector<Cost> CostVector;

     typedef std::vector<char> BoolVector;

     typedef std::vector<char> BoolVector;

     // Note: vector<char> is used instead of vector<bool> for efficiency reasons

     // Note: vector<char> is used instead of vector<bool> for efficiency reasons

   private:

   private:

     template <typename KT, typename VT>

     template <typename KT, typename VT>

     class StaticVectorMap {

     class StaticVectorMap {

     public:

     public:

       typedef KT Key;

       typedef KT Key;

       typedef VT Value;

       typedef VT Value;

       StaticVectorMap(std::vector<Value>& v) : _v(v) {}

       StaticVectorMap(std::vector<Value>& v) : _v(v) {}

       const Value& operator[](const Key& key) const {

       const Value& operator[](const Key& key) const {

         return _v[StaticDigraph::id(key)];

         return _v[StaticDigraph::id(key)];

       }

       }

       Value& operator[](const Key& key) {

       Value& operator[](const Key& key) {

         return _v[StaticDigraph::id(key)];

         return _v[StaticDigraph::id(key)];

       }

       }

       void set(const Key& key, const Value& val) {

       void set(const Key& key, const Value& val) {

         _v[StaticDigraph::id(key)] = val;

         _v[StaticDigraph::id(key)] = val;

       }

       }

     private:

     private:

     std::vector<IntPair> _arc_vec;

     std::vector<IntPair> _arc_vec;

     std::vector<Cost> _cost_vec;

     std::vector<Cost> _cost_vec;

     IntVector _id_vec;

     IntVector _id_vec;

     CostArcMap _cost_map;

     CostArcMap _cost_map;

     CostNodeMap _pi_map;

     CostNodeMap _pi_map;

   public:

   public:

     /// \brief Constant for infinite upper bounds (capacities).

     /// \brief Constant for infinite upper bounds (capacities).

     ///

     ///

     /// Constant for infinite upper bounds (capacities).

     /// Constant for infinite upper bounds (capacities).

     /// It is \c std::numeric_limits<Value>::infinity() if available,

     /// It is \c std::numeric_limits<Value>::infinity() if available,

     /// \c std::numeric_limits<Value>::max() otherwise.

     /// \c std::numeric_limits<Value>::max() otherwise.

       }

       }

       _supply[_node_id[s]] =  k;

       _supply[_node_id[s]] =  k;

       _supply[_node_id[t]] = -k;

       _supply[_node_id[t]] = -k;

       return *this;

       return *this;

     }

     }

     /// @}

     /// @}

     /// \name Execution control

     /// \name Execution control

     /// The algorithm can be executed using \ref run().

     /// The algorithm can be executed using \ref run().

       for (int j = limit; j != _res_arc_num; ++j) {

       for (int j = limit; j != _res_arc_num; ++j) {

         _lower[j] = 0;

         _lower[j] = 0;

         _upper[j] = INF;

         _upper[j] = INF;

         _cost[j] = 0;

         _cost[j] = 0;

         _cost[_reverse[j]] = 0;

         _cost[_reverse[j]] = 0;

       _have_lower = false;

       _have_lower = false;

       return *this;

       return *this;

     }

     }

     /// \brief Reset the internal data structures and all the parameters

     /// \brief Reset the internal data structures and all the parameters

       _lower.resize(_res_arc_num);

       _lower.resize(_res_arc_num);

       _upper.resize(_res_arc_num);

       _upper.resize(_res_arc_num);

       _cost.resize(_res_arc_num);

       _cost.resize(_res_arc_num);

       _supply.resize(_res_node_num);

       _supply.resize(_res_node_num);

       _res_cap.resize(_res_arc_num);

       _res_cap.resize(_res_arc_num);

       _pi.resize(_res_node_num);

       _pi.resize(_res_node_num);

       _arc_vec.reserve(_res_arc_num);

       _arc_vec.reserve(_res_arc_num);

       _cost_vec.reserve(_res_arc_num);

       _cost_vec.reserve(_res_arc_num);

         int fi = _arc_idf[a];

         int fi = _arc_idf[a];

         int bi = _arc_idb[a];

         int bi = _arc_idb[a];

         _reverse[fi] = bi;

         _reverse[fi] = bi;

         _reverse[bi] = fi;

         _reverse[bi] = fi;

       }

       }

       // Reset parameters

       // Reset parameters

       resetParams();

       resetParams();

       return *this;

       return *this;

     }

     }

       _sum_supply = 0;

       _sum_supply = 0;

       for (int i = 0; i != _root; ++i) {

       for (int i = 0; i != _root; ++i) {

         _sum_supply += _supply[i];

         _sum_supply += _supply[i];

       }

       }

       if (_sum_supply > 0) return INFEASIBLE;

       if (_sum_supply > 0) return INFEASIBLE;

       // Initialize vectors

       // Initialize vectors

       for (int i = 0; i != _res_node_num; ++i) {

       for (int i = 0; i != _res_node_num; ++i) {

         _pi[i] = 0;

         _pi[i] = 0;

       }

       }

       ValueVector excess(_supply);

       ValueVector excess(_supply);

       // Remove infinite upper bounds and check negative arcs

       // Remove infinite upper bounds and check negative arcs

       const Value MAX = std::numeric_limits<Value>::max();

       const Value MAX = std::numeric_limits<Value>::max();

       int last_out;

       int last_out;

       if (_have_lower) {

       if (_have_lower) {

         for (int i = 0; i != _root; ++i) {

         for (int i = 0; i != _root; ++i) {

           _res_cap[ra] = 0;

           _res_cap[ra] = 0;

           _cost[a] = 0;

           _cost[a] = 0;

           _cost[ra] = 0;

           _cost[ra] = 0;

         }

         }

       }

       }

       return OPTIMAL;

       return OPTIMAL;

     }

     }

     // Build a StaticDigraph structure containing the current

     // Build a StaticDigraph structure containing the current

     // residual network

     // residual network

     void buildResidualNetwork() {

     void buildResidualNetwork() {

       _arc_vec.clear();

       _arc_vec.clear();

       _cost_vec.clear();

       _cost_vec.clear();

     // Execute the "Simple Cycle Canceling" method

     // Execute the "Simple Cycle Canceling" method

     void startSimpleCycleCanceling() {

     void startSimpleCycleCanceling() {

       // Constants for computing the iteration limits

       // Constants for computing the iteration limits

       const int BF_FIRST_LIMIT  = 2;

       const int BF_FIRST_LIMIT  = 2;

       const double BF_LIMIT_FACTOR = 1.5;

       const double BF_LIMIT_FACTOR = 1.5;

       typedef StaticVectorMap<StaticDigraph::Arc, Value> FilterMap;

       typedef StaticVectorMap<StaticDigraph::Arc, Value> FilterMap;

       typedef FilterArcs<StaticDigraph, FilterMap> ResDigraph;

       typedef FilterArcs<StaticDigraph, FilterMap> ResDigraph;

       typedef StaticVectorMap<StaticDigraph::Node, StaticDigraph::Arc> PredMap;

       typedef StaticVectorMap<StaticDigraph::Node, StaticDigraph::Arc> PredMap;

       typedef typename BellmanFord<ResDigraph, CostArcMap>

       typedef typename BellmanFord<ResDigraph, CostArcMap>

         ::template SetDistMap<CostNodeMap>

         ::template SetDistMap<CostNodeMap>

         ::template SetPredMap<PredMap>::Create BF;

         ::template SetPredMap<PredMap>::Create BF;

       // Build the residual network

       // Build the residual network

       _arc_vec.clear();

       _arc_vec.clear();

       _cost_vec.clear();

       _cost_vec.clear();

       for (int j = 0; j != _res_arc_num; ++j) {

       for (int j = 0; j != _res_arc_num; ++j) {

         _arc_vec.push_back(IntPair(_source[j], _target[j]));

         _arc_vec.push_back(IntPair(_source[j], _target[j]));

     void startMinMeanCycleCanceling() {

     void startMinMeanCycleCanceling() {

       typedef SimplePath<StaticDigraph> SPath;

       typedef SimplePath<StaticDigraph> SPath;

       typedef typename SPath::ArcIt SPathArcIt;

       typedef typename SPath::ArcIt SPathArcIt;

       typedef typename HowardMmc<StaticDigraph, CostArcMap>

       typedef typename HowardMmc<StaticDigraph, CostArcMap>

         ::template SetPath<SPath>::Create MMC;

         ::template SetPath<SPath>::Create MMC;

       SPath cycle;

       SPath cycle;

       MMC mmc(_sgr, _cost_map);

       MMC mmc(_sgr, _cost_map);

       mmc.cycle(cycle);

       mmc.cycle(cycle);

       buildResidualNetwork();

       buildResidualNetwork();

       while (mmc.findCycleMean() && mmc.cycleCost() < 0) {

       while (mmc.findCycleMean() && mmc.cycleCost() < 0) {

           int j = _id_vec[_sgr.id(a)];

           int j = _id_vec[_sgr.id(a)];

           _res_cap[j] -= delta;

           _res_cap[j] -= delta;

           _res_cap[_reverse[j]] += delta;

           _res_cap[_reverse[j]] += delta;

         }

         }

         buildResidualNetwork();

         buildResidualNetwork();

       }

       }

     }

     }

     // Execute the "Cancel And Tighten" method

     // Execute the "Cancel And Tighten" method

           MMC mmc(_sgr, _cost_map);

           MMC mmc(_sgr, _cost_map);

           mmc.findCycleMean();

           mmc.findCycleMean();

           epsilon = -mmc.cycleMean();

           epsilon = -mmc.cycleMean();

           Cost cycle_cost = mmc.cycleCost();

           Cost cycle_cost = mmc.cycleCost();

           int cycle_size = mmc.cycleSize();

           int cycle_size = mmc.cycleSize();

           // Compute feasible potentials for the current epsilon

           // Compute feasible potentials for the current epsilon

           for (int i = 0; i != int(_cost_vec.size()); ++i) {

           for (int i = 0; i != int(_cost_vec.size()); ++i) {

             _cost_vec[i] = cycle_size * _cost_vec[i] - cycle_cost;

             _cost_vec[i] = cycle_size * _cost_vec[i] - cycle_cost;

           }

           }

           BF bf(_sgr, _cost_map);

           BF bf(_sgr, _cost_map);

           bf.init(0);

           bf.init(0);

           bf.start();

           bf.start();

           for (int u = 0; u != _res_node_num; ++u) {

           for (int u = 0; u != _res_node_num; ++u) {

             pi[u] = static_cast<double>(_pi[u]) / cycle_size;

             pi[u] = static_cast<double>(_pi[u]) / cycle_size;

           }

           }

           iter = limit;

           iter = limit;

         }

         }

       }

       }

     }

     }