RhodeCode

  /// \tparam V The value type used for flow amounts, capacity bounds

  /// and supply values in the algorithm. By default it is \c int.

  /// \tparam C The value type used for costs and potentials in the

  /// algorithm. By default it is the same as \c V.

  ///

  /// \warning Both value types must be signed and all input data must

  /// be integer.

  ///

  /// \note %NetworkSimplex provides five different pivot rule

  /// implementations, from which the most efficient one is used

  /// by default. For more information see \ref PivotRule.

  template <typename GR, typename V = int, typename C = V>

  class NetworkSimplex

  {

  public:

    typedef V Value;

    typedef C Cost;

  public:

    /// \brief Problem type constants for the \c run() function.

    ///

    /// Enum type containing the problem type constants that can be

    /// returned by the \ref run() function of the algorithm.

    enum ProblemType {

      /// The problem has no feasible solution (flow).

      INFEASIBLE,

      /// The problem has optimal solution (i.e. it is feasible and

      /// bounded), and the algorithm has found optimal flow and node

      /// potentials (primal and dual solutions).

      OPTIMAL,

      /// The objective function of the problem is unbounded, i.e.

      /// there is a directed cycle having negative total cost and

      /// In a major iteration a candidate list is built from eligible arcs

      /// in a wraparound fashion and in the following minor iterations

      /// the best eligible arc is selected from this list.

      CANDIDATE_LIST,

      /// The Altering Candidate List pivot rule.

      /// It is a modified version of the Candidate List method.

      /// It keeps only the several best eligible arcs from the former

      /// candidate list and extends this list in every iteration.

      ALTERING_LIST

    };

  private:

    TEMPLATE_DIGRAPH_TYPEDEFS(GR);

    typedef std::vector<Arc> ArcVector;

    typedef std::vector<Node> NodeVector;

    typedef std::vector<int> IntVector;

    typedef std::vector<bool> BoolVector;

    typedef std::vector<Cost> CostVector;

    // State constants for arcs

    enum ArcStateEnum {

      STATE_UPPER = -1,

      STATE_TREE  =  0,

      STATE_LOWER =  1

    };

  private:

    // Data related to the underlying digraph

    const GR &_graph;

    int _node_num;

    int _arc_num;

    // Parameters of the problem

    SupplyType _stype;

    Value _sum_supply;

    // Data structures for storing the digraph

    IntNodeMap _node_id;

    IntVector _source;

    IntVector _target;

    // Node and arc data

    CostVector _cost;

    CostVector _pi;

    // Data for storing the spanning tree structure

    IntVector _parent;

    IntVector _pred;

    IntVector _thread;

    IntVector _rev_thread;

    IntVector _succ_num;

    IntVector _last_succ;

    IntVector _dirty_revs;

    BoolVector _forward;

    IntVector _state;

    int _root;

    // Temporary data used in the current pivot iteration

    int in_arc, join, u_in, v_in, u_out, v_out;

        pop_heap( _candidates.begin(), _candidates.begin() + _curr_length,

                  _sort_func );

        _curr_length = std::min(_head_length, _curr_length - 1);

        return true;

      }

    }; //class AlteringListPivotRule

  public:

    /// \brief Constructor.

    ///

    /// The constructor of the class.

    ///

    /// \param graph The digraph the algorithm runs on.

    NetworkSimplex(const GR& graph) :

      INF(std::numeric_limits<Value>::has_infinity ?

          std::numeric_limits<Value>::infinity() :

          std::numeric_limits<Value>::max())

    {

      // Check the value types

      LEMON_ASSERT(std::numeric_limits<Value>::is_signed,

        "The flow type of NetworkSimplex must be signed");

      LEMON_ASSERT(std::numeric_limits<Cost>::is_signed,

        "The cost type of NetworkSimplex must be signed");

    }

    }

    /// \name Parameters

    /// The parameters of the algorithm can be specified using these

    /// functions.

    /// @{

    /// \brief Set the lower bounds on the arcs.

    ///

    /// This function sets the lower bounds on the arcs.

    /// If it is not used before calling \ref run(), the lower bounds

    /// will be set to zero on all arcs.

    ///

    /// \param map An arc map storing the lower bounds.

    /// Its \c Value type must be convertible to the \c Value type

    /// of the algorithm.

    ///

    /// \return <tt>(*this)</tt>

    template <typename LowerMap>

    NetworkSimplex& lowerMap(const LowerMap& map) {

      for (ArcIt a(_graph); a != INVALID; ++a) {

      }

      return *this;

    }

    /// \brief Set the upper bounds (capacities) on the arcs.

    ///

    /// This function sets the upper bounds (capacities) on the arcs.

    /// If it is not used before calling \ref run(), the upper bounds

    /// will be set to \ref INF on all arcs (i.e. the flow value will be

    /// unbounded from above on each arc).

    ///

    /// \param map An arc map storing the upper bounds.

    /// Its \c Value type must be convertible to the \c Value type

    /// of the algorithm.

    ///

    /// \return <tt>(*this)</tt>

    template<typename UpperMap>

    NetworkSimplex& upperMap(const UpperMap& map) {

      for (ArcIt a(_graph); a != INVALID; ++a) {

      }

      return *this;

    }

    /// \brief Set the costs of the arcs.

    ///

    /// This function sets the costs of the arcs.

    /// If it is not used before calling \ref run(), the costs

    /// will be set to \c 1 on all arcs.

    ///

    /// \param map An arc map storing the costs.

    /// Its \c Value type must be convertible to the \c Cost type

    /// of the algorithm.

    ///

    /// \return <tt>(*this)</tt>

    template<typename CostMap>

    NetworkSimplex& costMap(const CostMap& map) {

      for (ArcIt a(_graph); a != INVALID; ++a) {

      }

      return *this;

    }

    /// \brief Set the supply values of the nodes.

    ///

    /// This function sets the supply values of the nodes.

    /// If neither this function nor \ref stSupply() is used before

    /// calling \ref run(), the supply of each node will be set to zero.

    /// (It makes sense only if non-zero lower bounds are given.)

    ///

    /// \param map A node map storing the supply values.

    /// Its \c Value type must be convertible to the \c Value type

    /// of the algorithm.

    ///

    /// \return <tt>(*this)</tt>

    template<typename SupplyMap>

    NetworkSimplex& supplyMap(const SupplyMap& map) {

      for (NodeIt n(_graph); n != INVALID; ++n) {

      }

      return *this;

    }

    /// \brief Set single source and target nodes and a supply value.

    ///

    /// This function sets a single source node and a single target node

    /// and the required flow value.

    /// If neither this function nor \ref supplyMap() is used before

    /// calling \ref run(), the supply of each node will be set to zero.

    /// (It makes sense only if non-zero lower bounds are given.)

    ///

    /// Using this function has the same effect as using \ref supplyMap()

    /// with such a map in which \c k is assigned to \c s, \c -k is

    /// assigned to \c t and all other nodes have zero supply value.

    ///

    /// \param s The source node.

    /// \param t The target node.

    /// \param k The required amount of flow from node \c s to node \c t

    /// (i.e. the supply of \c s and the demand of \c t).

    ///

    /// \return <tt>(*this)</tt>

    NetworkSimplex& stSupply(const Node& s, const Node& t, Value k) {

      return *this;

    }

    /// \brief Set the type of the supply constraints.

    ///

    /// This function sets the type of the supply/demand constraints.

    /// If it is not used before calling \ref run(), the \ref GEQ supply

    /// type will be used.

    ///

    /// For more information see \ref SupplyType.

    ///

    /// \return <tt>(*this)</tt>

    NetworkSimplex& supplyType(SupplyType supply_type) {

      _stype = supply_type;

      return *this;

    }

    /// @}

    /// \name Execution Control

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

    /// @{

    /// \brief Run the algorithm.

    ///

    /// This function runs the algorithm.

    /// The paramters can be specified using functions \ref lowerMap(),

    /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(), 

    /// For example,

    /// \code

    ///   NetworkSimplex<ListDigraph> ns(graph);

    ///   ns.lowerMap(lower).upperMap(upper).costMap(cost)

    ///     .supplyMap(sup).run();

    /// \endcode

    ///

    /// This function can be called more than once. All the parameters

    /// that have been given are kept for the next call, unless

    /// \ref reset() is called, thus only the modified parameters

    /// have to be set again. See \ref reset() for examples.

    ///

    /// \param pivot_rule The pivot rule that will be used during the

    /// algorithm. For more information see \ref PivotRule.

    ///

    /// \return \c INFEASIBLE if no feasible flow exists,

    /// \n \c OPTIMAL if the problem has optimal solution

    /// (i.e. it is feasible and bounded), and the algorithm has found

    /// optimal flow and node potentials (primal and dual solutions),

    /// \n \c UNBOUNDED if the objective function of the problem is

    /// unbounded, i.e. there is a directed cycle having negative total

    /// cost and infinite upper bound.

    ///

    /// \see ProblemType, PivotRule

    ProblemType run(PivotRule pivot_rule = BLOCK_SEARCH) {

      if (!init()) return INFEASIBLE;

      return start(pivot_rule);

    }

    /// \brief Reset all the parameters that have been given before.

    ///

    /// This function resets all the paramaters that have been given

    /// before using functions \ref lowerMap(), \ref upperMap(),

    ///

    /// It is useful for multiple run() calls. If this function is not

    /// used, all the parameters given before are kept for the next

    /// \ref run() call.

    ///

    /// For example,

    /// \code

    ///   NetworkSimplex<ListDigraph> ns(graph);

    ///

    ///   // First run

    ///   ns.lowerMap(lower).upperMap(upper).costMap(cost)

    ///     .supplyMap(sup).run();

    ///

    ///   // Run again with modified cost map (reset() is not called,

    ///   // so only the cost map have to be set again)

    ///   cost[e] += 100;

    ///   ns.costMap(cost).run();

    ///

    ///   // Run again from scratch using reset()

    ///   // (the lower bounds will be set to zero on all arcs)

    ///   ns.reset();

    ///   ns.upperMap(capacity).costMap(cost)

    ///     .supplyMap(sup).run();

    /// \endcode

    ///

    /// \return <tt>(*this)</tt>

    NetworkSimplex& reset() {

      _stype = GEQ;

      return *this;

    }

    /// @}

    /// \name Query Functions

    /// The results of the algorithm can be obtained using these

    /// functions.\n

    /// The \ref run() function must be called before using them.

    /// @{

    /// \brief Return the total cost of the found flow.

    ///

    /// This function returns the total cost of the found flow.

    /// Its complexity is O(e).

    ///

    /// \note The return type of the function can be specified as a

    /// template parameter. For example,

    /// \code

    ///   ns.totalCost<double>();

    /// \endcode

    /// It is useful if the total cost cannot be stored in the \c Cost

    /// type of the algorithm, which is the default return type of the

    /// function.

    ///

    /// \pre \ref run() must be called before using this function.

      }

      return c;

    }

#ifndef DOXYGEN

    Cost totalCost() const {

      return totalCost<Cost>();

    }

#endif

    /// \brief Return the flow on the given arc.

    ///

    /// This function returns the flow on the given arc.

    ///

    /// \pre \ref run() must be called before using this function.

    Value flow(const Arc& a) const {

    }

    ///

    ///

    /// \pre \ref run() must be called before using this function.

    }

    /// \brief Return the potential (dual value) of the given node.

    ///

    /// This function returns the potential (dual value) of the

    /// given node.

    ///

    /// \pre \ref run() must be called before using this function.

    Cost potential(const Node& n) const {

    }

    ///

    ///

    /// \pre \ref run() must be called before using this function.

    }

    /// @}

  private:

    // Initialize internal data structures

    bool init() {

      if (_node_num == 0) return false;

      _sum_supply = 0;

          _sum_supply += _supply[i];

        }

      } else {

        }

      }

      // Initialize artifical cost

      Cost ART_COST;

      if (std::numeric_limits<Cost>::is_exact) {

        ART_COST = std::numeric_limits<Cost>::max() / 4 + 1;

      } else {

        ART_COST = std::numeric_limits<Cost>::min();

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

          if (_cost[i] > ART_COST) ART_COST = _cost[i];

        }

        ART_COST = (ART_COST + 1) * _node_num;

      }

      // Set data for the artificial root node

      _root = _node_num;

      _parent[_root] = -1;

      _pred[_root] = -1;

      _thread[_root] = 0;

      _rev_thread[0] = _root;

      _last_succ[_root] = _root - 1;

      _supply[_root] = -_sum_supply;

      // Add artificial arcs and initialize the spanning tree data structure

      for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {

        _thread[u] = u + 1;

        _rev_thread[u + 1] = u;

        _succ_num[u] = 1;

        _last_succ[u] = u;

        _cost[e] = ART_COST;

        _cap[e] = INF;

        _state[e] = STATE_TREE;

        if (_supply[u] > 0 || (_supply[u] == 0 && _sum_supply <= 0)) {

          _flow[e] = _supply[u];

          _forward[u] = true;

          _pi[u] = -ART_COST + _pi[_root];

        } else {

          _flow[e] = -_supply[u];

          _forward[u] = false;

          _pi[u] = ART_COST + _pi[_root];

        }

      }

      return true;

    }

      if (_sum_supply < 0) {

        for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {

          if (_supply[u] >= 0 && _flow[e] != 0) return INFEASIBLE;

        }

      }

      else if (_sum_supply > 0) {

        for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {

          if (_supply[u] <= 0 && _flow[e] != 0) return INFEASIBLE;

        }

      }

      else {

        for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {

          if (_flow[e] != 0) return INFEASIBLE;

        }

      }

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

        }

        }

      }

      return OPTIMAL;

    }

  }; //class NetworkSimplex

  ///@}

} //namespace lemon

#endif //LEMON_NETWORK_SIMPLEX_H

  "10  7    60   10    0    0   -3\n"

  "11 10    20   14    0    6  -20\n"

  "12 11    30   10    0    0  -10\n"

  "\n"

  "@attributes\n"

  "source 1\n"

  "target 12\n";

enum SupplyType {

  EQ,

  GEQ,

  LEQ

};

// Check the interface of an MCF algorithm

class McfClassConcept

{

public:

  template <typename MCF>

  struct Constraints {

    void constraints() {

      checkConcept<concepts::Digraph, GR>();

      MCF mcf(g);

      b = mcf.reset()

             .lowerMap(lower)

             .upperMap(upper)

             .costMap(cost)

             .supplyMap(sup)

             .stSupply(n, n, k)

             .run();

      c = const_mcf.totalCost();

      v = const_mcf.flow(a);

      c = const_mcf.potential(n);

    }

    typedef typename GR::Node Node;

    typedef typename GR::Arc Arc;

    typedef concepts::ReadMap<Arc, Cost> CAM;

    const GR &g;

    const CAM &cost;

    const NM ⊃

    const Node &n;

    const Arc &a;

    bool b;

  };

};

// Check the feasibility of the given flow (primal soluiton)

template < typename GR, typename LM, typename UM,

           typename SM, typename FM >

bool checkFlow( const GR& gr, const LM& lower, const UM& upper,

                const SM& supply, const FM& flow,

                SupplyType type = EQ )

{

  TEMPLATE_DIGRAPH_TYPEDEFS(GR);

  for (ArcIt e(gr); e != INVALID; ++e) {

    if (flow[e] < lower[e] || flow[e] > upper[e]) return false;

}

// Run a minimum cost flow algorithm and check the results

template < typename MCF, typename GR,

           typename LM, typename UM,

           typename CM, typename SM,

           typename PT >

void checkMcf( const MCF& mcf, PT mcf_result,

               const GR& gr, const LM& lower, const UM& upper,

               const CM& cost, const SM& supply,

               PT result, bool optimal, typename CM::Value total,

               const std::string &test_id = "",

               SupplyType type = EQ )

{

  check(mcf_result == result, "Wrong result " + test_id);

  if (optimal) {

          "The flow is not feasible " + test_id);

    check(mcf.totalCost() == total, "The flow is not optimal " + test_id);

          "Wrong potentials " + test_id);

  }

}

int main()

{

  // Check the interfaces

  {

    typedef concepts::Digraph GR;

  }

  // Run various MCF tests

  typedef ListDigraph Digraph;

  DIGRAPH_TYPEDEFS(ListDigraph);

  // Read the test digraph

  Digraph gr;

  Digraph::ArcMap<int> c(gr), l1(gr), l2(gr), l3(gr), u(gr);

  Digraph::NodeMap<int> s1(gr), s2(gr), s3(gr), s4(gr), s5(gr), s6(gr);

  ConstMap<Arc, int> cc(1), cu(std::numeric_limits<int>::max());

  Node v, w;

  std::istringstream input(test_lgf);

  DigraphReader<Digraph>(gr, input)

    .arcMap("cost", c)