-- | Maximum Flow algorithm
-- We are given a flow network G=(V,E) with source s and sink t where each
-- edge (u,v) in E has a nonnegative capacity c(u,v)>=0, and we wish to
-- find a flow of maximum value from s to t.
-- 
-- A flow in G=(V,E) is a real-valued function f:VxV->R that satisfies:
--
-- @
-- For all u,v in V, f(u,v)\<=c(u,v)
-- For all u,v in V, f(u,v)=-f(v,u)
-- For all u in V-{s,t}, Sum{f(u,v):v in V } = 0
-- @
--
-- The value of a flow f is defined as |f|=Sum {f(s,v)|v in V}, i.e.,
-- the total net flow out of the source.
--
-- In this module we implement the Edmonds-Karp algorithm, which is the
-- Ford-Fulkerson method but using the shortest path from s to t as the
-- augmenting path along which the flow is incremented.

module Data.Graph.Inductive.Query.MaxFlow(
    getRevEdges, augmentGraph, updAdjList, updateFlow, mfmg, mf, maxFlowgraph,
    maxFlow
) where


import Data.List

import Data.Graph.Inductive.Basic
import Data.Graph.Inductive.Graph
--import Data.Graph.Inductive.Tree
import Data.Graph.Inductive.Query.BFS

-- |
-- @
--                 i                                 0
-- For each edge a--->b this function returns edge b--->a .
--          i
-- Edges a\<--->b are ignored
--          j
-- @
getRevEdges :: (Num b,Ord b) => [(Node,Node)] -> [(Node,Node,b)]
getRevEdges [] = []
getRevEdges ((u,v):es) | notElem (v,u) es = (v,u,0):getRevEdges es
                       | otherwise        = getRevEdges (delete (v,u) es)

-- |
-- @
--                 i                                  0
-- For each edge a--->b insert into graph the edge a\<---b . Then change the
--                            i         (i,0,i)
-- label of every edge from a---->b to a------->b
-- @
-- 
-- where label (x,y,z)=(Max Capacity, Current flow, Residual capacity)
augmentGraph :: (DynGraph gr,Num b,Ord b) => gr a b -> gr a (b,b,b)
augmentGraph g = emap (\i->(i,0,i)) (insEdges (getRevEdges (edges g)) g)
                                                
-- | Given a successor or predecessor list for node u and given node v, find
-- the label corresponding to edge (u,v) and update the flow and residual
-- capacity of that edge's label. Then return the updated list.
updAdjList::(Num b,Ord b) => [((b,b,b),Node)]->Node->b->Bool->[((b,b,b),Node)]
updAdjList s v cf fwd | fwd == True = ((x,y+cf,z-cf),w):rs
                      | otherwise   = ((x,y-cf,z+cf),w):rs
                        where ((x,y,z),w) = head (filter (\(_,w')->v==w') s)
                              rs          = filter (\(_,w')->v/=w') s

-- | Update flow and residual capacity along augmenting path from s to t in
-- graph G. For a path [u,v,w,...] find the node u in G and its successor and
-- predecessor list, then update the corresponding edges (u,v) and (v,u) on
-- those lists by using the minimum residual capacity of the path.
updateFlow :: (DynGraph gr,Num b,Ord b) => Path -> b -> gr a (b,b,b) -> gr a (b,b,b)
updateFlow []	     _ g = g
updateFlow [_]       _ g = g
updateFlow (u:v:vs) cf g = case match u g of
                             (Nothing,g')        -> g'
                             (Just (p,u',l,s),g') -> (p',u',l,s') & g2
                                where g2 = updateFlow (v:vs) cf g'
                                      s' = updAdjList s v cf True
                                      p' = updAdjList p v cf False

-- | Compute the flow from s to t on a graph whose edges are labeled with
-- (x,y,z)=(max capacity,current flow,residual capacity) and all edges
-- are of the form a\<---->b. First compute the residual graph, that is,
-- delete those edges whose residual capacity is zero. Then compute the
-- shortest augmenting path from s to t, and finally update the flow and
-- residual capacity along that path by using the minimum capacity of
-- that path. Repeat this process until no shortest path from s to t exist.
mfmg :: (DynGraph gr,Num b,Ord b) => gr a (b,b,b) -> Node -> Node -> gr a (b,b,b)
mfmg g s t | augPath == [] = g
           | otherwise     = mfmg (updateFlow augPath minC g) s t
             where minC        = minimum (map ((\(_,_,z)->z).snd)(tail augLPath))
                   augPath     = map fst augLPath
                   LP augLPath = lesp s t gf
                   gf          = elfilter (\(_,_,z)->z/=0) g

-- | Compute the flow from s to t on a graph whose edges are labeled with
-- x, which is the max capacity and where not all edges need to be of the
-- form a\<---->b. Return the flow as a grap whose edges are labeled with
-- (x,y,z)=(max capacity,current flow,residual capacity) and all edges
-- are of the form a\<---->b
mf :: (DynGraph gr,Num b,Ord b) => gr a b -> Node -> Node -> gr a (b,b,b)
mf g s t = mfmg (augmentGraph g) s t

-- | Compute the maximum flow from s to t on a graph whose edges are labeled
-- with x, which is the max capacity and where not all edges need to be of
-- the form a\<---->b. Return the flow as a grap whose edges are labeled with
-- (y,x) = (current flow, max capacity).
maxFlowgraph :: (DynGraph gr,Num b,Ord b) => gr a b -> Node -> Node -> gr a (b,b)
maxFlowgraph g s t = emap (\(u,v,_)->(v,u)) g2
                           where g2 = elfilter (\(x,_,_)->x/=0) g1
                                 g1 = mf g s t

-- | Compute the value of a maximumflow
maxFlow :: (DynGraph gr,Num b,Ord b) => gr a b -> Node -> Node -> b
maxFlow g s t = foldr (+) 0 (map (\(_,_,(x,_))->x)(out (maxFlowgraph g s t) s))

------------------------------------------------------------------------------
-- Some test cases: clr595 is from the CLR textbook, page 595. The value of
-- the maximum flow for s=1 and t=6 (23) coincides with the example but the 
-- flow itself is slightly different since the textbook does not compute the 
-- shortest augmenting path from s to t, but just any path. However remember 
-- that for a given flow graph the maximum flow is not unique.
-- (gr595 is defined in GraphData.hs)
------------------------------------------------------------------------------