deba@2276: /* -*- C++ -*-
deba@2276:  *
deba@2276:  * This file is a part of LEMON, a generic C++ optimization library
deba@2276:  *
alpar@2391:  * Copyright (C) 2003-2007
deba@2276:  * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
deba@2276:  * (Egervary Research Group on Combinatorial Optimization, EGRES).
deba@2276:  *
deba@2276:  * Permission to use, modify and distribute this software is granted
deba@2276:  * provided that this copyright notice appears in all copies. For
deba@2276:  * precise terms see the accompanying LICENSE file.
deba@2276:  *
deba@2276:  * This software is provided "AS IS" with no warranty of any kind,
deba@2276:  * express or implied, and with no claim as to its suitability for any
deba@2276:  * purpose.
deba@2276:  *
deba@2276:  */
deba@2276: 
deba@2276: #ifndef LEMON_MIN_COST_FLOW_H
deba@2276: #define LEMON_MIN_COST_FLOW_H
deba@2276: 
deba@2376: ///\ingroup min_cost_flow 
deba@2276: ///
alpar@2350: ///\file
alpar@2350: ///\brief An algorithm for finding a flow of value \c k (for
deba@2276: ///small values of \c k) having minimal total cost
deba@2276: 
deba@2276: 
deba@2276: #include <lemon/dijkstra.h>
deba@2276: #include <lemon/graph_adaptor.h>
deba@2276: #include <lemon/maps.h>
deba@2276: #include <vector>
deba@2276: 
deba@2276: namespace lemon {
deba@2276: 
deba@2376:   /// \addtogroup min_cost_flow
deba@2276:   /// @{
deba@2276: 
deba@2276:   /// \brief Implementation of an algorithm for finding a flow of
deba@2276:   /// value \c k (for small values of \c k) having minimal total cost
deba@2276:   /// between two nodes
deba@2276:   /// 
deba@2276:   ///
deba@2276:   /// The \ref lemon::SspMinCostFlow "Successive Shortest Path Minimum
deba@2276:   /// Cost Flow" implements an algorithm for finding a flow of value
deba@2276:   /// \c k having minimal total cost from a given source node to a
deba@2276:   /// given target node in a directed graph with a cost function on
deba@2276:   /// the edges. To this end, the edge-capacities and edge-costs have
deba@2276:   /// to be nonnegative.  The edge-capacities should be integers, but
deba@2276:   /// the edge-costs can be integers, reals or of other comparable
deba@2276:   /// numeric type.  This algorithm is intended to be used only for
deba@2276:   /// small values of \c k, since it is only polynomial in k, not in
deba@2276:   /// the length of k (which is log k): in order to find the minimum
deba@2276:   /// cost flow of value \c k it finds the minimum cost flow of value
deba@2276:   /// \c i for every \c i between 0 and \c k.
deba@2276:   ///
deba@2276:   ///\param Graph The directed graph type the algorithm runs on.
deba@2276:   ///\param LengthMap The type of the length map.
deba@2276:   ///\param CapacityMap The capacity map type.
deba@2276:   ///
deba@2276:   ///\author Attila Bernath
deba@2276:   template <typename Graph, typename LengthMap, typename CapacityMap>
deba@2276:   class SspMinCostFlow {
deba@2276: 
deba@2276:     typedef typename LengthMap::Value Length;
deba@2276: 
deba@2276:     //Warning: this should be integer type
deba@2276:     typedef typename CapacityMap::Value Capacity;
deba@2276:     
deba@2276:     typedef typename Graph::Node Node;
deba@2276:     typedef typename Graph::NodeIt NodeIt;
deba@2276:     typedef typename Graph::Edge Edge;
deba@2276:     typedef typename Graph::OutEdgeIt OutEdgeIt;
deba@2276:     typedef typename Graph::template EdgeMap<int> EdgeIntMap;
deba@2276: 
deba@2276:     typedef ResGraphAdaptor<const Graph,int,CapacityMap,EdgeIntMap> ResGW;
deba@2276:     typedef typename ResGW::Edge ResGraphEdge;
deba@2276: 
deba@2276:   protected:
deba@2276: 
deba@2276:     const Graph& g;
deba@2276:     const LengthMap& length;
deba@2276:     const CapacityMap& capacity;
deba@2276: 
deba@2276:     EdgeIntMap flow; 
deba@2276:     typedef typename Graph::template NodeMap<Length> PotentialMap;
deba@2276:     PotentialMap potential;
deba@2276: 
deba@2276:     Node s;
deba@2276:     Node t;
deba@2276:     
deba@2276:     Length total_length;
deba@2276: 
deba@2276:     class ModLengthMap {   
deba@2276:       typedef typename Graph::template NodeMap<Length> NodeMap;
deba@2276:       const ResGW& g;
deba@2276:       const LengthMap &length;
deba@2276:       const NodeMap &pot;
deba@2276:     public :
deba@2276:       typedef typename LengthMap::Key Key;
deba@2276:       typedef typename LengthMap::Value Value;
deba@2276: 
deba@2276:       ModLengthMap(const ResGW& _g, 
deba@2276: 		   const LengthMap &_length, const NodeMap &_pot) : 
deba@2276: 	g(_g), /*rev(_rev),*/ length(_length), pot(_pot) { }
deba@2276: 	
deba@2276:       Value operator[](typename ResGW::Edge e) const {     
deba@2276: 	if (g.forward(e))
deba@2276: 	  return  length[e]-(pot[g.target(e)]-pot[g.source(e)]);   
deba@2276: 	else
deba@2276: 	  return -length[e]-(pot[g.target(e)]-pot[g.source(e)]);   
deba@2276:       }     
deba@2276: 	
deba@2276:     }; //ModLengthMap
deba@2276: 
deba@2276:     ResGW res_graph;
deba@2276:     ModLengthMap mod_length;
deba@2276:     Dijkstra<ResGW, ModLengthMap> dijkstra;
deba@2276: 
deba@2276:   public :
deba@2276: 
deba@2276:     /// \brief The constructor of the class.
deba@2276:     ///
deba@2276:     /// \param _g The directed graph the algorithm runs on. 
deba@2276:     /// \param _length The length (cost) of the edges. 
deba@2276:     /// \param _cap The capacity of the edges. 
deba@2276:     /// \param _s Source node.
deba@2276:     /// \param _t Target node.
deba@2276:     SspMinCostFlow(Graph& _g, LengthMap& _length, CapacityMap& _cap, 
deba@2276:                    Node _s, Node _t) : 
deba@2276:       g(_g), length(_length), capacity(_cap), flow(_g), potential(_g), 
deba@2276:       s(_s), t(_t), 
deba@2276:       res_graph(g, capacity, flow), 
deba@2276:       mod_length(res_graph, length, potential),
deba@2276:       dijkstra(res_graph, mod_length) { 
deba@2276:       reset();
deba@2276:     }
deba@2276:     
deba@2276:     /// \brief Tries to augment the flow between s and t by 1. The
deba@2276:     /// return value shows if the augmentation is successful.
deba@2276:     bool augment() {
deba@2276:       dijkstra.run(s);
deba@2276:       if (!dijkstra.reached(t)) {
deba@2276: 
deba@2276: 	//Unsuccessful augmentation.
deba@2276: 	return false;
deba@2276:       } else {
deba@2276: 
deba@2276: 	//We have to change the potential
deba@2276: 	for(typename ResGW::NodeIt n(res_graph); n!=INVALID; ++n)
deba@2276: 	  potential.set(n, potential[n]+dijkstra.distMap()[n]);
deba@2276: 	
deba@2276: 	//Augmenting on the shortest path
deba@2276: 	Node n=t;
deba@2276: 	ResGraphEdge e;
deba@2276: 	while (n!=s){
deba@2276: 	  e = dijkstra.predEdge(n);
deba@2276: 	  n = dijkstra.predNode(n);
deba@2276: 	  res_graph.augment(e,1);
deba@2276: 	  //Let's update the total length
deba@2276: 	  if (res_graph.forward(e))
deba@2276: 	    total_length += length[e];
deba@2276: 	  else 
deba@2276: 	    total_length -= length[e];	    
deba@2276: 	}
deba@2276: 
deba@2276: 	return true;
deba@2276:       }
deba@2276:     }
deba@2276:     
deba@2276:     /// \brief Runs the algorithm.
deba@2276:     ///
deba@2276:     /// Runs the algorithm.
deba@2276:     /// Returns k if there is a flow of value at least k from s to t.
deba@2276:     /// Otherwise it returns the maximum value of a flow from s to t.
deba@2276:     ///
deba@2276:     /// \param k The value of the flow we are looking for.
deba@2276:     ///
deba@2276:     /// \todo May be it does make sense to be able to start with a
deba@2276:     /// nonzero feasible primal-dual solution pair as well.
deba@2276:     ///
deba@2276:     /// \todo If the actual flow value is bigger than k, then
deba@2276:     /// everything is cleared and the algorithm starts from zero
deba@2276:     /// flow. Is it a good approach?
deba@2276:     int run(int k) {
deba@2276:       if (flowValue()>k) reset();
deba@2276:       while (flowValue()<k && augment()) { }
deba@2276:       return flowValue();
deba@2276:     }
deba@2276: 
deba@2276:     /// \brief The class is reset to zero flow and potential. The
deba@2276:     /// class is reset to zero flow and potential.
deba@2276:     void reset() {
deba@2276:       total_length=0;
deba@2276:       for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) flow.set(e, 0);
deba@2276:       for (typename Graph::NodeIt n(g); n!=INVALID; ++n) potential.set(n, 0);  
deba@2276:     }
deba@2276: 
deba@2276:     /// \brief Returns the value of the actual flow. 
deba@2276:     int flowValue() const {
deba@2276:       int i=0;
deba@2276:       for (typename Graph::OutEdgeIt e(g, s); e!=INVALID; ++e) i+=flow[e];
deba@2276:       for (typename Graph::InEdgeIt e(g, s); e!=INVALID; ++e) i-=flow[e];
deba@2276:       return i;
deba@2276:     }
deba@2276: 
deba@2276:     /// \brief Total cost of the found flow.
deba@2276:     ///
deba@2276:     /// This function gives back the total cost of the found flow.
deba@2276:     Length totalLength(){
deba@2276:       return total_length;
deba@2276:     }
deba@2276: 
deba@2276:     /// \brief Returns a const reference to the EdgeMap \c flow. 
deba@2276:     ///
deba@2276:     /// Returns a const reference to the EdgeMap \c flow. 
deba@2276:     const EdgeIntMap &getFlow() const { return flow;}
deba@2276: 
deba@2276:     /// \brief Returns a const reference to the NodeMap \c potential
deba@2276:     /// (the dual solution).
deba@2276:     ///
deba@2276:     /// Returns a const reference to the NodeMap \c potential (the
deba@2276:     /// dual solution).
deba@2276:     const PotentialMap &getPotential() const { return potential;}
deba@2276: 
deba@2276:     /// \brief Checking the complementary slackness optimality criteria.
deba@2276:     ///
deba@2276:     /// This function checks, whether the given flow and potential 
deba@2276:     /// satisfy the complementary slackness conditions (i.e. these are optimal).
deba@2276:     /// This function only checks optimality, doesn't bother with feasibility.
deba@2276:     /// For testing purpose.
deba@2276:     bool checkComplementarySlackness(){
deba@2276:       Length mod_pot;
deba@2276:       Length fl_e;
deba@2276:         for(typename Graph::EdgeIt e(g); e!=INVALID; ++e) {
deba@2276: 	//C^{\Pi}_{i,j}
deba@2276: 	mod_pot = length[e]-potential[g.target(e)]+potential[g.source(e)];
deba@2276: 	fl_e = flow[e];
deba@2276: 	if (0<fl_e && fl_e<capacity[e]) {
deba@2276: 	  /// \todo better comparison is needed for real types, moreover, 
deba@2276: 	  /// this comparison here is superfluous.
deba@2276: 	  if (mod_pot != 0)
deba@2276: 	    return false;
deba@2276: 	} 
deba@2276: 	else {
deba@2276: 	  if (mod_pot > 0 && fl_e != 0)
deba@2276: 	    return false;
deba@2276: 	  if (mod_pot < 0 && fl_e != capacity[e])
deba@2276: 	    return false;
deba@2276: 	}
deba@2276:       }
deba@2276:       return true;
deba@2276:     }
deba@2276:     
deba@2276:   }; //class SspMinCostFlow
deba@2276: 
deba@2276:   ///@}
deba@2276: 
deba@2276: } //namespace lemon
deba@2276: 
deba@2276: #endif //LEMON_MIN_COST_FLOW_H