const ResGW& G;

       const ResGW& g;

       const LengthMap &ol;

       const LengthMap &length;

 	if (G.forward(e))

 	if (g.forward(e))

 	  return  ol[e]-(pot[G.head(e)]-pot[G.tail(e)]);   

 	  return  length[e]-(pot[g.head(e)]-pot[g.tail(e)]);   

 	  return -ol[e]-(pot[G.head(e)]-pot[G.tail(e)]);   

 	  return -length[e]-(pot[g.head(e)]-pot[g.tail(e)]);   

       ModLengthMap(const ResGW& _G,

     }; //ModLengthMap

 		   const LengthMap &o,  const NodeMap &p) : 

 	G(_G), /*rev(_rev),*/ ol(o), pot(p){}; 

     ResGW res_graph;

     };//ModLengthMap

     ModLengthMap mod_length;

     Dijkstra<ResGW, ModLengthMap> dijkstra;

   protected:

     //Input

     const Graph& G;

     const LengthMap& length;

     const CapacityMap& capacity;

     //auxiliary variables

     //To store the flow

     EdgeIntMap flow; 

     //To store the potential (dual variables)

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

     PotentialMap potential;

     Length total_length;

     /// The constructor of the class.

     /*! \brief The constructor of the class.

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

     \param _g The directed graph the algorithm runs on. 

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

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

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

     \param _cap The capacity of the edges. 

     MinCostFlow(Graph& _G, LengthMap& _length, CapacityMap& _cap) : G(_G), 

     \param _s Source node.

       length(_length), capacity(_cap), flow(_G), potential(_G){ }

     \param _t Target node.

     */

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

     ///Runs the algorithm.

 		Node _s, Node _t) : 

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

     ///Runs the algorithm.

       s(_s), t(_t), 

     ///Returns k if there is a flow of value at least k edge-disjoint 

       res_graph(g, capacity, flow), 

     ///from s to t.

       mod_length(res_graph, length, potential),

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

       dijkstra(res_graph, mod_length) { 

     ///

       reset();

     ///\param s The source node.

       }

     ///\param t The target node.

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

     /*! Tries to augment the flow between s and t by 1.

     ///

       The return value shows if the augmentation is successful.

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

      */

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

     bool augment() {

     int run(Node s, Node t, int k) {

       dijkstra.run(s);

       if (!dijkstra.reached(t)) {

       //Resetting variables from previous runs

       total_length = 0;

 	//Unsuccessful augmentation.

 	return false;

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

       } else {

       //Initialize the potential to zero

 	//We have to change the potential

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

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

 	  potential[n] += dijkstra.distMap()[n];

       //We need a residual graph

       ResGW res_graph(G, capacity, flow);

       ModLengthMap mod_length(res_graph, length, potential);

       Dijkstra<ResGW, ModLengthMap> dijkstra(res_graph, mod_length);

       int i;

       for (i=0; i<k; ++i){

 	dijkstra.run(s);

 	if (!dijkstra.reached(t)){

 	  //There are no flow of value k from s to t

 	  break;

 	};

 	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);  

     }

     /*! 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];

     /// Total weight of the found flow.

     /// Gives back the total weight of the found flow.

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

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

     ///Assumes that \c run() has been run and nothing changed since then.

     ///Returns a const reference to the NodeMap \c potential (the dual solution).

     /*! \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).

     Returns a const reference to the NodeMap \c potential (the dual solution).

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

     */

     /// Checking the complementary slackness optimality criteria

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

     ///This function checks, whether the given solution is optimal

     This function checks, whether the given flow and potential 

     ///If executed after the call of \c run() then it should return with true.

     satisfiy the complementary slackness cnditions (i.e. these are optimal).

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

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

     ///It is meant for testing purposes.

     For testing purpose.

     ///

     */

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

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

 	mod_pot = length[e]-potential[G.head(e)]+potential[G.tail(e)];

 	mod_pot = length[e]-potential[g.head(e)]+potential[g.tail(e)];