-- | 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) => [Edge] -> [LEdge b]
getRevEdges [] = []
getRevEdges ((u,v):es) | (v,u) `notElem` 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) => 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) => Adj (b,b,b) -> Node -> b -> Bool -> Adj (b,b,b)
updAdjList s v cf fwd = rs ++ ((x,y+cf',z-cf'),w) : rs'
  where
    (rs, ((x,y,z),w):rs') = break ((v==) . snd) s

    cf' = if fwd
             then cf
             else negate cf

-- | 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) => 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
  | null 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 = mfmg (augmentGraph g)

-- | 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 graph 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))
                     . elfilter (\(x,_,_) -> x/=0 )
                     $ 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 = sum (map (fst . edgeLabel) (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)
------------------------------------------------------------------------------