/* -*- C++ -*-

 *

 * This file is a part of LEMON, a generic C++ optimization library

 *

 * Copyright (C) 2003-2006

 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport

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

 *

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

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

 * precise terms see the accompanying LICENSE file.

 *

 * This software is provided "AS IS" with no warranty of any kind,

 * express or implied, and with no claim as to its suitability for any

 * purpose.

 *

 */

#ifndef LEMON_MIN_COST_FLOW_H

#define LEMON_MIN_COST_FLOW_H

///\ingroup flowalgs 

///

///\file \brief An algorithm for finding a flow of value \c k (for

///small values of \c k) having minimal total cost

#include <lemon/dijkstra.h>

#include <lemon/graph_adaptor.h>

#include <lemon/maps.h>

#include <vector>

namespace lemon {

  /// \addtogroup flowalgs

  /// @{

  /// \brief Implementation of an algorithm for finding a flow of

  /// value \c k (for small values of \c k) having minimal total cost

  /// between two nodes

  /// 

  ///

  /// The \ref lemon::SspMinCostFlow "Successive Shortest Path Minimum

  /// Cost Flow" implements an algorithm for finding a flow of value

  /// \c k having minimal total cost from a given source node to a

  /// given target node in a directed graph with a cost function on

  /// the edges. To this end, the edge-capacities and edge-costs have

  /// to be nonnegative.  The edge-capacities should be integers, but

  /// the edge-costs can be integers, reals or of other comparable

  /// numeric type.  This algorithm is intended to be used only for

  /// small values of \c k, since it is only polynomial in k, not in

  /// the length of k (which is log k): in order to find the minimum

  /// cost flow of value \c k it finds the minimum cost flow of value

  /// \c i for every \c i between 0 and \c k.

  ///

  ///\param Graph The directed graph type the algorithm runs on.

  ///\param LengthMap The type of the length map.

  ///\param CapacityMap The capacity map type.

  ///

  ///\author Attila Bernath

  template <typename Graph, typename LengthMap, typename CapacityMap>

  class SspMinCostFlow {

    typedef typename LengthMap::Value Length;

    //Warning: this should be integer type

    typedef typename CapacityMap::Value Capacity;

    typedef typename Graph::Node Node;

    typedef typename Graph::NodeIt NodeIt;

    typedef typename Graph::Edge Edge;

    typedef typename Graph::OutEdgeIt OutEdgeIt;

    typedef typename Graph::template EdgeMap<int> EdgeIntMap;

    typedef ResGraphAdaptor<const Graph,int,CapacityMap,EdgeIntMap> ResGW;

    typedef typename ResGW::Edge ResGraphEdge;

  protected:

    const Graph& g;

    const LengthMap& length;

    const CapacityMap& capacity;

    EdgeIntMap flow; 

    typedef typename Graph::template NodeMap<Length> PotentialMap;

    PotentialMap potential;

    Node s;

    Node t;

    Length total_length;

    class ModLengthMap {   

      typedef typename Graph::template NodeMap<Length> NodeMap;

      const ResGW& g;

      const LengthMap &length;

      const NodeMap &pot;

    public :

      typedef typename LengthMap::Key Key;

      typedef typename LengthMap::Value Value;

      ModLengthMap(const ResGW& _g, 

		   const LengthMap &_length, const NodeMap &_pot) : 

	g(_g), /*rev(_rev),*/ length(_length), pot(_pot) { }

      Value operator[](typename ResGW::Edge e) const {     

	if (g.forward(e))

	  return  length[e]-(pot[g.target(e)]-pot[g.source(e)]);   

	else

	  return -length[e]-(pot[g.target(e)]-pot[g.source(e)]);   

      }     

    }; //ModLengthMap

    ResGW res_graph;

    ModLengthMap mod_length;

    Dijkstra<ResGW, ModLengthMap> dijkstra;

  public :

    /// \brief The constructor of the class.

    ///

    /// \param _g The directed graph the algorithm runs on. 

    /// \param _length The length (cost) of the edges. 

    /// \param _cap The capacity of the edges. 

    /// \param _s Source node.

    /// \param _t Target node.

    SspMinCostFlow(Graph& _g, LengthMap& _length, CapacityMap& _cap, 

                   Node _s, Node _t) : 

      g(_g), length(_length), capacity(_cap), flow(_g), potential(_g), 

      s(_s), t(_t), 

      res_graph(g, capacity, flow), 

      mod_length(res_graph, length, potential),

      dijkstra(res_graph, mod_length) { 

      reset();

    }

    /// \brief Tries to augment the flow between s and t by 1. The

    /// return value shows if the augmentation is successful.

    bool augment() {

      dijkstra.run(s);

      if (!dijkstra.reached(t)) {

	//Unsuccessful augmentation.

	return false;

      } else {

	//We have to change the potential

	for(typename ResGW::NodeIt n(res_graph); n!=INVALID; ++n)

	  potential.set(n, potential[n]+dijkstra.distMap()[n]);

	//Augmenting on the shortest path

	Node n=t;

	ResGraphEdge e;

	while (n!=s){

	  e = dijkstra.predEdge(n);

	  n = dijkstra.predNode(n);

	  res_graph.augment(e,1);

	  //Let's update the total length

	  if (res_graph.forward(e))

	    total_length += length[e];

	  else 

	    total_length -= length[e];	    

	}

	return true;

      }

    }

    /// \brief Runs the algorithm.

    ///

    /// Runs the algorithm.

    /// Returns k if there is a flow of value at least k from s to t.

    /// Otherwise it returns the maximum value of a flow from s to t.

    ///

    /// \param k The value of the flow we are looking for.

    ///

    /// \todo May be it does make sense to be able to start with a

    /// nonzero feasible primal-dual solution pair as well.

    ///

    /// \todo If the actual flow value is bigger than k, then

    /// everything is cleared and the algorithm starts from zero

    /// flow. Is it a good approach?

    int run(int k) {

      if (flowValue()>k) reset();

      while (flowValue()<k && augment()) { }

      return flowValue();

    }

    /// \brief The class is reset to zero flow and potential. The

    /// class is reset to zero flow and potential.

    void reset() {

      total_length=0;

      for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) flow.set(e, 0);

      for (typename Graph::NodeIt n(g); n!=INVALID; ++n) potential.set(n, 0);  

    }

    /// \brief Returns the value of the actual flow. 

    int flowValue() const {

      int i=0;

      for (typename Graph::OutEdgeIt e(g, s); e!=INVALID; ++e) i+=flow[e];

      for (typename Graph::InEdgeIt e(g, s); e!=INVALID; ++e) i-=flow[e];

      return i;

    }

    /// \brief Total cost of the found flow.

    ///

    /// This function gives back the total cost of the found flow.

    Length totalLength(){

      return total_length;

    }

    /// \brief Returns a const reference to the EdgeMap \c flow. 

    ///

    /// Returns a const reference to the EdgeMap \c flow. 

    const EdgeIntMap &getFlow() const { return flow;}

    /// \brief Returns a const reference to the NodeMap \c potential

    /// (the dual solution).

    ///

    /// Returns a const reference to the NodeMap \c potential (the

    /// dual solution).

    const PotentialMap &getPotential() const { return potential;}

    /// \brief Checking the complementary slackness optimality criteria.

    ///

    /// This function checks, whether the given flow and potential 

    /// satisfy the complementary slackness conditions (i.e. these are optimal).

    /// This function only checks optimality, doesn't bother with feasibility.

    /// For testing purpose.

    bool checkComplementarySlackness(){

      Length mod_pot;

      Length fl_e;

        for(typename Graph::EdgeIt e(g); e!=INVALID; ++e) {

	//C^{\Pi}_{i,j}

	mod_pot = length[e]-potential[g.target(e)]+potential[g.source(e)];

	fl_e = flow[e];

	if (0<fl_e && fl_e<capacity[e]) {

	  /// \todo better comparison is needed for real types, moreover, 

	  /// this comparison here is superfluous.

	  if (mod_pot != 0)

	    return false;

	} 

	else {

	  if (mod_pot > 0 && fl_e != 0)

	    return false;

	  if (mod_pot < 0 && fl_e != capacity[e])

	    return false;

	}

      }

      return true;

    }

  }; //class SspMinCostFlow

  ///@}

} //namespace lemon

#endif //LEMON_MIN_COST_FLOW_H