/*******************************************************************************
+
+  LEDA 3.5
+
+  _min_cost_flow.c
+
+  This file is part of the LEDA research version (LEDA-R) that can be 
+  used free of charge in academic research and teaching. Any commercial
+  use of this software requires a license which is distributed by the
+  LEDA Software GmbH, Postfach 151101, 66041 Saarbruecken, FRG
+  (fax +49 681 31104).
+
+  Copyright (c) 1991-1997  by  Max-Planck-Institut fuer Informatik
+  Im Stadtwald, 66123 Saarbruecken, Germany     
+  All rights reserved.
+ 
*******************************************************************************/
//------------------------------------------------------------------------------
// bool MIN_COST_FLOW(G, lcap, ucap, cost, supply, flow)
//------------------------------------------------------------------------------
//
// computes a flow of minimal total cost
// w.r.t. lower capacities      lcap
//        upper capacities      ucap
//        edge  costs           cost
//        node supplies/demands supply
// 
// returns 
//         true  if there exists a feasible flow
//         false otherwise
//
// Algorithm: successive shortest path augmentation + capacity scaling
//
// based on   J. Edmonds and R. M. Karp: 
//            "Theoretical Improvements in Algorithmic Efficiency for Network 
//            Flow Problems", Journal of the  ACM, Vol. 19, No. 2, 1972
//
//            R.K. Ahuja, T.L. Magnanti, J.B. Orlin: 
//            "Network Flows", Section 10.2, Prentice Hall Publ. Company, 1993
//
//
// implementation based on MIN_COST_MAX_FLOW by S. Naeher  (December 1993)
//
// G. Klau, S. Naeher  (1996)
//------------------------------------------------------------------------------


#include <LEDA/graph_alg.h>
#include <LEDA/node_pq.h>

typedef int num_type;


static node DIJKSTRA_IN_RESIDUAL(const graph& G, node s, int delta, 
                                 const edge_array<num_type>& cost,
                                 const edge_array<num_type>& cap_minus_flow,
                                 const edge_array<num_type>& flow,
                                 const node_array<num_type>& excess,
                                 node_array<num_type>& pi,
                                 node_array<edge>& pred)
{
 /* Computes a shortest path from node s until a node t with excess <= -delta
    is reached in the delta-residual graph of network G. Edges are traversed
    in both directions using the forall_in/out_edges iteration macros.The 
    residual capacity of an edge e is cap_minus_flow[e] if e is used from 
    source(e) to target(e) and flow[e] otherwise. The reduced cost of an
    edge e = (u,v) is equal to cost[e] + pi[u] - pi[v] if it is used from 
    u and -cost[e] + pi[v] - pi[u] if it is used from v.

    precondition: every edge in the delta-residual network has non-negative 
                  reduced cost.
    
    returns t (nil,if there is no path from s to a node with excess<-delta)

    and

    pred[v] = predecessor edge of v on the shortest path from s to t

    pi[v]   = pi'[v] + dist[v] - dist[t], if v is permanently labeled
              pi'[v],                     if v is temporarily labeled
              (here pi' is the old potential vector)

 */

  node_array<num_type> dist(G,MAXINT);
  node_pq<num_type> PQ(G);
  node v;
  edge e;

  node t = nil;

  node_list perm_labeled; // list of permanently labeled nodes

  dist[s] = 0;
  PQ.insert(s,0);

  while (!PQ.empty())
  { node u = PQ.del_min();

    perm_labeled.append(u);

    if (excess[u] <= -delta) 
    { t = u;
      break;
     }

    num_type du = dist[u] + pi[u]; 

    forall_out_edges(e,u) 
      if (cap_minus_flow[e] >= delta)  // e in delta-residual graph
      { v = target(e);
        num_type c = du - pi[v] + cost[e];
        if (c < dist[v])
        { if (dist[v] == MAXINT) // v has not been reached before
            PQ.insert(v,c);
          else
            PQ.decrease_p(v,c);
          dist[v] = c;
          pred[v] = e;
         }
       }

    forall_in_edges(e,u) 
      if (flow[e] >= delta)  // e in delta-residual graph
      { v = source(e);
        num_type c = du - pi[v] - cost[e];
        if (c < dist[v])
        { if (dist[v] == MAXINT) // v has not been reached before
            PQ.insert(v,c);
          else
            PQ.decrease_p(v,c);
          dist[v] = c;
          pred[v] = e;
         }
       }

   }

   if (t != nil) // update potentials of all permanently labeled nodes
   { num_type dt = dist[t];
     forall(v,perm_labeled) pi[v] += dist[v] - dt;
    }

   return t;
}




bool MIN_COST_FLOW(graph& G, const edge_array<num_type>& lcap,
                             const edge_array<num_type>& ucap,
                             const edge_array<num_type>& cost0,
                             const node_array<num_type>& supply,
                             edge_array<num_type>&       flow0)
{
  int n = G.number_of_nodes();
  int m = G.number_of_edges();

  int count = 0;		// number of augmentations
  node v;
  edge e;

  num_type C = 0;		// maximal cost of any path in G
  num_type U = 0;		// maximal supply/demand or capacity

  forall_edges(e, G) { 
    if (cost0[e] > C)  C = cost0[e];
    if (cost0[e] < -C) C = -cost0[e];
    if (ucap[e] > U)   U = ucap[e];
  }

  C = C * G.number_of_nodes();

  
  list<edge> art_edges;		// list of artificial edges (see below)
  list<edge> neg_edges;		// list of negative cost edges

  // add artificial edges to guarantee existence of a path of infinite
  // capacity (MAXINT) and huge cost (C) between any pair of nodes  (v,w)
  // with b(v) > 0 and b(w) < 0 (Ahuja/Magnanti/Orlin, section 9.1, p. 297)
 
  node first_node = G.first_node(); 
  forall_nodes(v, G) if (v != first_node) 
  { art_edges.append(G.new_edge(v, first_node));
    art_edges.append(G.new_edge(first_node, v));
   }

  node_array<num_type>  pi(G, 0);	 // node potential
  node_array<num_type>  excess(G, 0);    // excess
  node_array<edge>      pred(G, nil);    // predecessor edge on shortest path
  edge_array<num_type>  flow(G, 0);
  edge_array<num_type>  cap(G, 0);
  edge_array<num_type>  cost(G, 0);
  edge_array<num_type>  cap_minus_flow(G, 0);

  forall_nodes(v, G) 
  { num_type s = supply[v];
    if ( s > U) U =  s;
    if (-s > U) U = -s;
    excess[v] = s;
   }
  
  forall(e, art_edges) 
  { cap[e] = MAXINT;
    cost[e] = C;
    cap_minus_flow[e] = MAXINT;
   }


  forall_edges(e, G) 
    if (cap[e] != MAXINT) 
    { // original edge
      cap[e] = ucap[e];
      cost[e] = cost0[e];
      if (cost[e] < 0) neg_edges.append(e);

      // removing nonzero lower bounds (Ahuja/Magnanti/Orlin, section 2.4,
      // p. 39 [fig. 2.19 is not correct!]
      
      num_type bound = lcap[e];

      if (bound != 0) 
      {	// nonzero lower bound
	node i = source(e);
	node j = target(e);
	excess[i] -= bound;
	excess[j] += bound;
	cap[e] -= bound;
       }
      cap_minus_flow[e] = cap[e];
     }


  // eliminate negative cost edges by edge reversals 
  // (Ahuja/Magnanti/Orlin, section 2.4, p. 40)

  forall(e, neg_edges) 
  { node u = source(e);
    node v = target(e);
    excess[u] -= cap[e];
    excess[v] += cap[e];
    cost[e] = -cost[e];
    G.rev_edge(e);
   }


  int delta = 1;

  U = (U*m)/(n*n);

  while (U >= delta) delta *= 2;

  while (delta > 0) 
  { // delta scaling phase

    forall_edges(e, G) 
    { // Saturate all edges entering the delta residual network which have
      // a negative reduced edge cost. Then only the reverse edge (with positive
      // reduced edge cost) will stay in the residual network. 

      node u = source(e);
      node v = target(e);
      num_type c = cost[e] + pi[u] - pi[v];

      if (cap_minus_flow[e] >= delta && c < 0)
        { excess[v] += cap_minus_flow[e];
          excess[u] -= cap_minus_flow[e];
          flow[e] = cap[e];
          cap_minus_flow[e] = 0;
         }
      else
        if (flow[e] >= delta && c > 0)
        { excess[u] += flow[e];
          excess[v] -= flow[e];
          flow[e] = 0;
          cap_minus_flow[e] = cap[e];
         }

     }


    // As long as there is a node s with excess >= delta compute a minimum
    // cost augmenting path from s to a sink node t with excess <= -delta in 
    // the delta-residual network and augment the flow by at least delta units 
    // along P (if there are nodes with excess > delta and excess < -delta
    // respectively, the existence of P is guaranteed by the artificial edges 
    // inserted at the begining of the algorithm)

    node s;
    node t;
  
    forall_nodes(s, G)
    while (excess[s] >= delta)
    { 
      count++;
  
      t = DIJKSTRA_IN_RESIDUAL(G,s,delta,cost,cap_minus_flow,flow,excess,pi,pred);
      if (t == nil) break; // there is no node with excess < -delta
  
      // compute maximal amount f of flow that can be pushed through P
  
      num_type f = excess[s];
  
      v = t; 
      while (v != s)
      { e = pred[v]; 
        num_type rcap;
        if (v == target(e))
           { rcap = cap_minus_flow[e];
             v = source(e);
            }
        else
           { rcap = flow[e];
             v = target(e);
            }
        if (rcap < f) f = rcap;
       }
  
      if (f > -excess[t]) f = -excess[t];
  
  
      // add f units of flow along the augmenting path 
  
      v = t; 
      while (v != s)
      { e = pred[v]; 
        if (v==target(e))
           { flow[e] += f;
             cap_minus_flow[e] -= f;
             v = source(e);
            }
        else
           { flow[e] -= f;
             cap_minus_flow[e] += f;
             v = target(e);
            }
       }
  
       excess[s] -= f;
       excess[t] += f;
  
    } // end of delta-phase
  
   delta /= 2;
  
  } // end of scaling
  


  edge x;

  // restore negative edges
  forall(x, neg_edges) 
  { flow[x] = cap[x] - flow[x];
    G.rev_edge(x);
   } 


  // delete artificial edges and check feasibility of flow

  bool feasible = true;

  forall(x, art_edges) 
  { if (flow[x] != 0) feasible = false; 
    G.del_edge(x); 
   }

  // adjust flow

  forall_edges(x, G) flow0[x] = flow[x] + lcap[x];

  return feasible;
}



//------------------------------------------------------------------------------
// Variants
//------------------------------------------------------------------------------


//------------------------------------------------------------------------------
// without lower capacity bounds
//------------------------------------------------------------------------------

bool MIN_COST_FLOW(graph& G, const edge_array<int>& cap,
                             const edge_array<int>& cost,
                             const node_array<int>& supply,
                                   edge_array<int>& flow)
{ edge_array<int> lcap(G,0);
  return MIN_COST_FLOW(G,lcap,cap,cost,supply,flow); 
 }
   



//------------------------------------------------------------------------------
//  min cost maxflow
//------------------------------------------------------------------------------

int MIN_COST_MAX_FLOW(graph& G, node s, node t, const edge_array<int>& cap, 
                                                const edge_array<int>& cost, 
                                                      edge_array<int>& flow)
{ edge_array<int> lcap(G,0);
  node_array<int> supply(G,0);
  int f = MAX_FLOW(G,s,t,cap,flow);
  supply[s] =  f;
  supply[t] = -f;
  MIN_COST_FLOW(G,lcap,cap,cost,supply,flow); 
  return f;
 }

