// -*- c++ -*-
#ifndef HUGO_MINCOSTFLOW_H
#define HUGO_MINCOSTFLOW_H

///\ingroup galgs
///\file
///\brief An algorithm for finding the minimum cost flow of given value in an uncapacitated network

#include <hugo/dijkstra.h>
#include <hugo/graph_wrapper.h>
#include <hugo/maps.h>
#include <vector>
#include <list>
#include <hugo/for_each_macros.h>
#include <hugo/unionfind.h>

namespace hugo {

/// \addtogroup galgs
/// @{

  ///\brief Implementation of an algorithm for solving the minimum cost general
  /// flow problem in an uncapacitated network
  /// 
  ///
  /// The class \ref hugo::MinCostFlow "MinCostFlow" implements
  /// an algorithm for solving the following general minimum cost flow problem>
  /// 
  ///
  ///
  /// \warning It is assumed here that the problem has a feasible solution
  ///
  /// The range of the cost (weight) function is nonnegative reals but 
  /// the range of capacity function is the set of nonnegative integers. 
  /// It is not a polinomial time algorithm for counting the minimum cost
  /// maximal flow, since it counts the minimum cost flow for every value 0..M
  /// where \c M is the value of the maximal flow.
  ///
  ///\author Attila Bernath
  template <typename Graph, typename CostMap, typename SupplyDemandMap>
  class MinCostFlow {

    typedef typename CostMap::ValueType Cost;


    typedef typename SupplyDemandMap::ValueType SupplyDemand;
    
    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<SupplyDemand> FlowMap;
    typedef ConstMap<Edge,SupplyDemand> ConstEdgeMap;

    //    typedef ConstMap<Edge,int> ConstMap;

    typedef ResGraphWrapper<const Graph,int,ConstEdgeMap,FlowMap> ResGraph;
    typedef typename ResGraph::Edge ResGraphEdge;

    class ModCostMap {   
      //typedef typename ResGraph::template NodeMap<Cost> NodeMap;
      typedef typename Graph::template NodeMap<Cost> NodeMap;
      const ResGraph& res_graph;
      //      const EdgeIntMap& rev;
      const CostMap &ol;
      const NodeMap &pot;
    public :
      typedef typename CostMap::KeyType KeyType;
      typedef typename CostMap::ValueType ValueType;
	
      ValueType operator[](typename ResGraph::Edge e) const {     
	if (res_graph.forward(e))
	  return  ol[e]-(pot[res_graph.head(e)]-pot[res_graph.tail(e)]);   
	else
	  return -ol[e]-(pot[res_graph.head(e)]-pot[res_graph.tail(e)]);   
      }     
	
      ModCostMap(const ResGraph& _res_graph,
		   const CostMap &o,  const NodeMap &p) : 
	res_graph(_res_graph), /*rev(_rev),*/ ol(o), pot(p){}; 
    };//ModCostMap


  protected:
    
    //Input
    const Graph& graph;
    const CostMap& cost;
    const SupplyDemandMap& supply_demand;//supply or demand of nodes


    //auxiliary variables

    //To store the flow
    FlowMap flow; 
    //To store the potentila (dual variables)
    typename Graph::template NodeMap<Cost> potential;
    //To store excess-deficit values
    SupplyDemandMap excess_deficit;
    

    Cost total_cost;


  public :


    MinCostFlow(Graph& _graph, CostMap& _cost, SupplyDemandMap& _supply_demand) : graph(_graph), 
      cost(_cost), supply_demand(_supply_demand), flow(_graph), potential(_graph){ }

    
    ///Runs the algorithm.

    ///Runs the algorithm.

    ///\todo May be it does make sense to be able to start with a nonzero 
    /// feasible primal-dual solution pair as well.
    void run() {

      //Resetting variables from previous runs
      //total_cost = 0;

      typedef typename Graph::template NodeMap<int> HeapMap;
      typedef Heap< Node, SupplyDemand, typename Graph::template NodeMap<int>,
	std::greater<SupplyDemand> > 	HeapType;

      //A heap for the excess nodes
      HeapMap excess_nodes_map(graph,-1);
      HeapType excess_nodes(excess_nodes_map);

      //A heap for the deficit nodes
      HeapMap deficit_nodes_map(graph,-1);
      HeapType deficit_nodes(deficit_nodes_map);

      //A container to store nonabundant arcs
      list<Edge> nonabundant_arcs;

	
      FOR_EACH_LOC(typename Graph::EdgeIt, e, graph){
	flow.set(e,0);
	nonabundant_arcs.push_back(e);
      }

      //Initial value for delta
      SupplyDemand delta = 0;

      typedef UnionFindEnum<Node, Graph::template NodeMap> UFE;

      //A union-find structure to store the abundant components
      UFE::MapType abund_comp_map(graph);
      UFE abundant_components(abund_comp_map);



      FOR_EACH_LOC(typename Graph::NodeIt, n, graph){
       	excess_deficit.set(n,supply_demand[n]);
	//A supply node
	if (excess_deficit[n] > 0){
	  excess_nodes.push(n,excess_deficit[n]);
	}
	//A demand node
	if (excess_deficit[n] < 0){
	  deficit_nodes.push(n, - excess_deficit[n]);
	}
	//Finding out starting value of delta
	if (delta < abs(excess_deficit[n])){
	  delta = abs(excess_deficit[n]);
	}
	//Initialize the copy of the Dijkstra potential to zero
	potential.set(n,0);
	//Every single point is an abundant component initially 
	abundant_components.insert(n);
      }

      //It'll be allright as an initial value, though this value 
      //can be the maximum deficit here
      SupplyDemand max_excess = delta;
      
      ///\bug This is a serious cheat here, before we have an uncapacitated ResGraph
      ConstEdgeMap const_inf_map(MAX_INT);
      
      //We need a residual graph which is uncapacitated
      ResGraph res_graph(graph, const_inf_map, flow);
      
      //An EdgeMap to tell which arcs are abundant
      template typename Graph::EdgeMap<bool> abundant_arcs(graph);

      //Let's construct the sugraph consisting only of the abundant edges
      typedef ConstMap< typename Graph::Node, bool > ConstNodeMap;
      ConstNodeMap const_true_map(true);
      typedef SubGraphWrapper< Graph, ConstNodeMap, 
	 template typename Graph::EdgeMap<bool> > 
	AbundantGraph;
      AbundantGraph abundant_graph(graph, const_true_map, abundant_arcs );
      
      //Let's construct the residual graph for the abundant graph
      typedef ResGraphWrapper<const AbundantGraph,int,CapacityMap,EdgeIntMap> 
	ResAbGraph;
      //Again uncapacitated
      ResAbGraph res_ab_graph(abundant_graph, const_inf_map, flow);
      
      //We need things for the bfs
      typename ResAbGraph::NodeMap<bool> bfs_reached(res_ab_graph);
      typename ResAbGraph::NodeMap<typename ResAbGraph::Edge> 
	bfs_pred(res_ab_graph); 
      NullMap<typename ResAbGraph::Node, int> bfs_dist_dummy(res_ab_graph);
      //We want to run bfs-es (more) on this graph 'res_ab_graph'
      Bfs < ResAbGraph , 
	typename ResAbGraph::NodeMap<bool>, 
	typename ResAbGraph::NodeMap<typename ResAbGraph::Edge>,
	NullMap<typename ResAbGraph::Node, int> > 
	bfs(res_ab_graph, bfs_reached, bfs_pred, bfs_dist_dummy);
      
      ModCostMap mod_cost(res_graph, cost, potential);

      Dijkstra<ResGraph, ModCostMap> dijkstra(res_graph, mod_cost);


      while (max_excess > 0){

	//Reset delta if still too big
	if (8*number_of_nodes*max_excess <= delta){
	  delta = max_excess;
	  
	}

	/*
	 * Beginning of the delta scaling phase 
	*/
	//Merge and stuff
	{
	  SupplyDemand buf=8*number_of_nodes*delta;
	  list<Edge>::iterator i = nonabundant_arcs.begin();
	  while ( i != nonabundant_arcs.end() ){
	    if (flow[i]>=buf){
	      Node a = abundant_components.find(res_graph.head(i));
	      Node b = abundant_components.find(res_graph.tail(i));
	      //Merge
	      if (a != b){
		abundant_components.join(a,b);
		//We want to push the smaller
		//Which has greater absolut value excess/deficit
		Node root=(abs(excess_deficit[a])>abs(excess_deficit[b]))?a:b;
		//Which is the other
		Node non_root = ( a == root ) ? b : a ;
		abundant_components.makeRep(root);
		SupplyDemand qty_to_augment = abs(excess_deficit[non_root]); 
		//Push the positive value
		if (excess_deficit[non_root] < 0)
		  swap(root, non_root);
		//If the non_root node has excess/deficit at all
		if (qty_to_augment>0){
		  //Find path and augment
		  bfs.run(non_root);
		  //root should be reached
		  
		  //Augmenting on the found path
		  Node n=root;
		  ResGraphEdge e;
		  while (n!=non_root){
		    e = bfs_pred(n);
		    n = res_graph.tail(e);
		    res_graph.augment(e,qty_to_augment);
		  }
	  
		  //We know that non_root had positive excess
		  excess_nodes[non_root] -= qty_to_augment;
		  //But what about root node
		  //It might have been positive and so became larger
		  if (excess_deficit[root]>0){
		    excess_nodes[root] += qty_to_augment;
		  }
		  else{
		    //Or negative but not turned into positive
		    deficit_nodes[root] -= qty_to_augment;
		  }

		  //Update the excess_deficit map
		  excess_deficit[non_root] -= qty_to_augment;
		  excess_deficit[root] += qty_to_augment;

		  
		}
	      }
	      //What happens to i?
	      //Marci and Zsolt says I shouldn't do such things
	      nonabundant_arcs.erase(i++);
	      abundant_arcs[i] = true;
	    }
	    else
	      ++i;
	  }
	}


	Node s = excess_nodes.top(); 
	SupplyDemand max_excess = excess_nodes[s];
	Node t = deficit_nodes.top(); 
	if (max_excess < deficit_nodes[t]){
	  max_excess = deficit_nodes[t];
	}


	while(max_excess > (n-1)*delta/n){
	  
	  
	  //s es t valasztasa
	  
	  //Dijkstra part	
	  dijkstra.run(s);
	  
	  /*We know from theory that t can be reached
	  if (!dijkstra.reached(t)){
	    //There are no k paths from s to t
	    break;
	  };
	  */
	  
	  //We have to change the potential
	  FOR_EACH_LOC(typename ResGraph::NodeIt, n, res_graph){
	    potential[n] += dijkstra.distMap()[n];
	  }


	  //Augmenting on the sortest path
	  Node n=t;
	  ResGraphEdge e;
	  while (n!=s){
	    e = dijkstra.pred(n);
	    n = dijkstra.predNode(n);
	    res_graph.augment(e,delta);
	    /*
	    //Let's update the total cost
	    if (res_graph.forward(e))
	      total_cost += cost[e];
	    else 
	      total_cost -= cost[e];	    
	    */
	  }
	  
	  //Update the excess_deficit map
	  excess_deficit[s] -= delta;
	  excess_deficit[t] += delta;
	  

	  //Update the excess_nodes heap
	  if (delta >= excess_nodes[s]){
	    if (delta > excess_nodes[s])
	      deficit_nodes.push(s,delta - excess_nodes[s]);
	    excess_nodes.pop();
	    
	  } 
	  else{
	    excess_nodes[s] -= delta;
	  }
	  //Update the deficit_nodes heap
	  if (delta >= deficit_nodes[t]){
	    if (delta > deficit_nodes[t])
	      excess_nodes.push(t,delta - deficit_nodes[t]);
	    deficit_nodes.pop();
	    
	  } 
	  else{
	    deficit_nodes[t] -= delta;
	  }
	  //Dijkstra part ends here
	  
	  //Choose s and t again
	  s = excess_nodes.top(); 
	  max_excess = excess_nodes[s];
	  t = deficit_nodes.top(); 
	  if (max_excess < deficit_nodes[t]){
	    max_excess = deficit_nodes[t];
	  }

	}

	/*
	 * End of the delta scaling phase 
	*/

	//Whatever this means
	delta = delta / 2;

	/*This is not necessary here
	//Update the max_excess
	max_excess = 0;
	FOR_EACH_LOC(typename Graph::NodeIt, n, graph){
	  if (max_excess < excess_deficit[n]){
	    max_excess = excess_deficit[n];
	  }
	}
	*/

	  
      }//while(max_excess > 0)
      

      return i;
    }




    ///This function gives back the total cost of the found paths.
    ///Assumes that \c run() has been run and nothing changed since then.
    Cost totalCost(){
      return total_cost;
    }

    ///Returns a const reference to the EdgeMap \c flow. \pre \ref run() must
    ///be called before using this function.
    const EdgeIntMap &getFlow() const { return flow;}

  ///Returns a const reference to the NodeMap \c potential (the dual solution).
    /// \pre \ref run() must be called before using this function.
    const EdgeIntMap &getPotential() const { return potential;}

    ///This function checks, whether the given solution is optimal
    ///Running after a \c run() should return with true
    ///In this "state of the art" this only check optimality, doesn't bother with feasibility
    ///
    ///\todo Is this OK here?
    bool checkComplementarySlackness(){
      Cost mod_pot;
      Cost fl_e;
      FOR_EACH_LOC(typename Graph::EdgeIt, e, graph){
	//C^{\Pi}_{i,j}
	mod_pot = cost[e]-potential[graph.head(e)]+potential[graph.tail(e)];
	fl_e = flow[e];
	//	std::cout << fl_e << std::endl;
	if (0<fl_e && fl_e<capacity[e]){
	  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 MinCostFlow

  ///@}

} //namespace hugo

#endif //HUGO_MINCOSTFLOW_H