{-# LANGUAGE LambdaCase #-}
----------------------------------------------------------------------------
-- |
-- Module     : Algebra.Graph.Bipartite.AdjacencyMap.Algorithm
-- Copyright  : (c) Andrey Mokhov 2016-2022
-- License    : MIT (see the file LICENSE)
-- Maintainer : andrey.mokhov@gmail.com
-- Stability  : experimental
--
-- __Alga__ is a library for algebraic construction and manipulation of graphs
-- in Haskell. See <https://github.com/snowleopard/alga-paper this paper> for
-- the motivation behind the library, the underlying theory, and
-- implementation details.
--
-- This module provides several basic algorithms on undirected bipartite graphs.
----------------------------------------------------------------------------
module Algebra.Graph.Bipartite.AdjacencyMap.Algorithm (
    -- * Bipartiteness test
    OddCycle, detectParts,

    -- * Matchings
    Matching, pairOfLeft, pairOfRight, matching, isMatchingOf, matchingSize,
    maxMatching,

    -- * Vertex covers
    VertexCover, isVertexCoverOf, vertexCoverSize, minVertexCover,

    -- * Independent sets
    IndependentSet, isIndependentSetOf, independentSetSize, maxIndependentSet,

    -- * Miscellaneous
    augmentingPath, consistentMatching
    ) where

import Algebra.Graph.Bipartite.AdjacencyMap

import Control.Monad             (guard, when)
import Control.Monad.Trans.Class (lift)
import Control.Monad.Trans.Maybe (MaybeT(..))
import Control.Monad.Trans.State (State, runState, get, put, modify)
import Control.Monad.ST          (ST, runST)
import Data.Either               (fromLeft)
import Data.Foldable             (asum, foldl')
import Data.Functor              (($>))
import Data.List                 (sort)
import Data.Maybe                (fromJust)
import Data.STRef                (STRef, newSTRef, readSTRef, writeSTRef, modifySTRef)
import GHC.Generics

import qualified Algebra.Graph.AdjacencyMap as AM

import qualified Data.Map.Strict as Map
import qualified Data.Set        as Set
import qualified Data.Sequence   as Seq

import Data.Map.Strict (Map)
import Data.Set        (Set)
import Data.Sequence   (Seq, ViewL (..), (|>))

-- TODO: Make this representation type-safe
-- | A cycle of odd length. For example, @[1,2,3]@ represents the cycle
-- @1@ @->@ @2@ @->@ @3@ @->@ @1@.
type OddCycle a = [a]

data Part = LeftPart | RightPart deriving (Show, Eq)

otherPart :: Part -> Part
otherPart LeftPart  = RightPart
otherPart RightPart = LeftPart

-- | Test the bipartiteness of a given "Algebra.Graph.AdjacencyMap". In case of
-- success, return an 'AdjacencyMap' with the same set of edges and each vertex
-- marked with the part it belongs to. In case of failure, return any cycle of
-- odd length in the graph.
--
-- The returned partition is lexicographically smallest, assuming that vertices
-- of the left part precede all the vertices of the right part.
--
-- The returned cycle is optimal in the following sense: there exists a path
-- that is either empty or ends in a vertex adjacent to the first vertex in the
-- cycle, such that all vertices in @path@ @++@ @cycle@ are distinct and
-- @path@ @++@ @cycle@ is lexicographically smallest among all such pairs of
-- paths and cycles.
--
-- /Note/: since 'AdjacencyMap' represents /undirected/ bipartite graphs, all
-- edges in the input graph are treated as undirected. See the examples and the
-- correctness property for a clarification.
--
-- Complexity: /O((n + m) * log(n))/ time and /O(n + m)/ memory.
--
-- @
-- detectParts 'Algebra.Graph.AdjacencyMap.empty'                                       == Right 'empty'
-- detectParts ('Algebra.Graph.AdjacencyMap.vertex' x)                                  == Right ('leftVertex' x)
-- detectParts ('Algebra.Graph.AdjacencyMap.edge' x x)                                  == Left [x]
-- detectParts ('Algebra.Graph.AdjacencyMap.edge' 1 2)                                  == Right ('edge' 1 2)
-- detectParts (1 * (2 + 3))                               == Right ('edges' [(1,2), (1,3)])
-- detectParts (1 * 2 * 3)                                 == Left [1, 2, 3]
-- detectParts ((1 + 3) * (2 + 4) + 6 * 5)                 == Right ('swap' (1 + 3) * (2 + 4) + 'swap' 5 * 6)
-- detectParts ((1 * 3 * 4) + 2 * (1 + 2))                 == Left [2]
-- detectParts ('Algebra.Graph.AdjacencyMap.clique' [1..10])                            == Left [1, 2, 3]
-- detectParts ('Algebra.Graph.AdjacencyMap.circuit' [1..10])                           == Right ('circuit' [(x, x + 1) | x <- [1,3,5,7,9]])
-- detectParts ('Algebra.Graph.AdjacencyMap.circuit' [1..11])                           == Left [1..11]
-- detectParts ('Algebra.Graph.AdjacencyMap.biclique' [] xs)                            == Right ('vertices' xs [])
-- detectParts ('Algebra.Graph.AdjacencyMap.biclique' ('map' Left (x:xs)) ('map' Right ys)) == Right ('biclique' ('map' Left (x:xs)) ('map' Right ys))
-- 'isRight' (detectParts ('Algebra.Graph.AdjacencyMap.star' x ys))                       == 'notElem' x ys
-- 'isRight' (detectParts ('fromBipartite' ('toBipartite' x)))   == True
-- @
--
-- The correctness of 'detectParts' can be expressed by the following property:
--
-- @
-- let undirected = 'Algebra.Graph.AdjacencyMap.symmetricClosure' input in
-- case detectParts input of
--     Left cycle -> 'mod' (length cycle) 2 == 1 && 'Algebra.Graph.AdjacencyMap.isSubgraphOf' ('Algebra.Graph.AdjacencyMap.circuit' cycle) undirected
--     Right result -> 'Algebra.Graph.AdjacencyMap.gmap' 'Data.Either.Extra.fromEither' ('fromBipartite' result) == undirected
-- @
detectParts :: Ord a => AM.AdjacencyMap a -> Either (OddCycle a) (AdjacencyMap a a)
detectParts x = case runState (runMaybeT dfs) Map.empty of
    (Nothing, partMap) -> Right $ toBipartiteWith (toEither partMap) g
    (Just c , _      ) -> Left  $ oddCycle c
  where
    -- g :: AM.AdjacencyMap a
    g = AM.symmetricClosure x

    -- type PartMap a = Map a Part
    -- type PartMonad a = MaybeT (State (PartMap a)) [a]
    -- dfs :: PartMonad a
    dfs = asum [ processVertex v | v <- AM.vertexList g ]

    -- processVertex :: a -> PartMonad a
    processVertex v = do partMap <- lift get
                         guard (Map.notMember v partMap)
                         inVertex LeftPart v

    -- inVertex :: Part -> a -> PartMonad a
    inVertex vertexPart v = (v :) <$> do
        lift $ modify (Map.insert v vertexPart)
        let otherVertexPart = otherPart vertexPart
        asum [ onEdge otherVertexPart u | u <- Set.toAscList (AM.postSet v g) ]

    {-# INLINE onEdge #-}
    -- onEdge :: Part -> a -> PartMonad a
    onEdge vertexPart v = do partMap <- lift get
                             case Map.lookup v partMap of
                                 Nothing   -> inVertex vertexPart v
                                 Just part -> do guard (vertexPart /= part)
                                                 return [v] -- found a cycle!

    -- toEither :: PartMap a -> a -> Either a a
    toEither partMap v = case fromJust (Map.lookup v partMap) of
                             LeftPart  -> Left  v
                             RightPart -> Right v

    -- oddCycle :: [a] -> [a]
    oddCycle pathToCycle = init $ dropWhile (/= lastVertex) pathToCycle
      where
        lastVertex = last pathToCycle

-- | A /matching/ is a set of pairwise non-adjacent edges between the two parts
-- of a bipartite graph.
--
-- The 'Show' instance is defined using the 'matching' function, with the edges
-- listed in the ascending order of left vertices.
--
-- @
-- show ('matching' [])                 == "matching []"
-- show ('matching' [(2,\'a\'), (1,\'b\')]) == "matching [(1,\'b\'),(2,\'a\')]"
-- @
data Matching a b = Matching {
    -- | The map of vertices covered by the matching in the left part to their
    -- neighbours in the right part.
    -- Complexity: /O(1)/ time.
    --
    -- @
    -- pairOfLeft ('matching' [])                 == Map.'Data.Map.Strict.empty'
    -- pairOfLeft ('matching' [(2,\'a\'), (1,\'b\')]) == Map.'Data.Map.Strict.fromList' [(1,\'b\'), (2,\'a\')]
    -- Map.'Map.size' . pairOfLeft                    == Map.'Map.size' . pairOfRight
    -- @
    pairOfLeft  :: Map a b,

    -- | The map of vertices covered by the matching in the right part to their
    -- neighbours in the left part.
    -- Complexity: /O(1)/.
    --
    -- @
    -- pairOfRight ('matching' [])                 == Map.'Data.Map.Strict.empty'
    -- pairOfRight ('matching' [(2,\'a\'), (1,\'b\')]) == Map.'Data.Map.Strict.fromList' [(\'a\',2), (\'b\',1)]
    -- Map.'Map.size' . pairOfRight                    == Map.'Map.size' . pairOfLeft
    -- @
    pairOfRight :: Map b a
} deriving Generic

instance (Show a, Show b) => Show (Matching a b) where
    showsPrec _ m = showString "matching " . showList (Map.toAscList $ pairOfLeft m)

instance (Eq a, Eq b) => Eq (Matching a b) where
    x == y = pairOfLeft x == pairOfLeft y

instance (Ord a, Ord b) => Ord (Matching a b) where
    compare x y = compare (pairOfLeft x) (pairOfLeft y)

addEdgeUnsafe :: (Ord a, Ord b) => a -> b -> Matching a b -> Matching a b
addEdgeUnsafe a b (Matching ab ba) = Matching (Map.insert a b ab) (Map.insert b a ba)

addEdge :: (Ord a, Ord b) => a -> b -> Matching a b -> Matching a b
addEdge a b (Matching ab ba) = addEdgeUnsafe a b (Matching ab' ba')
    where
        ab' = case b `Map.lookup` ba of
                  Nothing -> Map.delete a ab
                  Just a' -> Map.delete a (Map.delete a' ab)
        ba' = case a `Map.lookup` ab of
                  Nothing -> Map.delete b ba
                  Just b' -> Map.delete b (Map.delete b' ba)

leftCovered :: Ord a => a -> Matching a b -> Bool
leftCovered a = Map.member a . pairOfLeft

-- | Construct a 'Matching' from a list of edges.
-- Complexity: /O(L * log(L))/ time, where /L/ is the length of the given list.
--
-- Edges that appear closer to the end of the list supersede all previous edges.
-- That is, if two edges from the list share a vertex, the one that appears
-- closer to the beginning is ignored.
--
-- @
-- 'pairOfLeft'  (matching [])                     == Map.'Data.Map.Strict.empty'
-- 'pairOfRight' (matching [])                     == Map.'Data.Map.Strict.empty'
-- 'pairOfLeft'  (matching [(2,\'a\'), (1,\'b\')])     == Map.'Data.Map.Strict.fromList' [(2,\'a\'), (1,\'b\')]
-- 'pairOfLeft'  (matching [(1,\'a\'), (1,\'b\')])     == Map.'Data.Map.Strict.singleton' 1 \'b\'
-- matching [(1,\'a\'), (1,\'b\'), (2,\'b\'), (2,\'a\')] == matching [(2,\'a\')]
-- @
matching :: (Ord a, Ord b) => [(a, b)] -> Matching a b
matching = foldl' (flip (uncurry addEdge)) (Matching Map.empty Map.empty)

-- | Check if a given 'Matching' is a valid /matching/ of a bipartite graph.
-- Complexity: /O(S * log(n))/, where /S/ is the size of the matching.
--
-- @
-- isMatchingOf ('matching' []) x               == True
-- isMatchingOf ('matching' xs) 'empty'           == 'null' xs
-- isMatchingOf ('matching' [(x,y)]) ('edge' x y) == True
-- isMatchingOf ('matching' [(1,2)]) ('edge' 2 1) == False
-- @
isMatchingOf :: (Ord a, Ord b) => Matching a b -> AdjacencyMap a b -> Bool
isMatchingOf m@(Matching ab _) g = consistentMatching m
    && and [ hasEdge a b g | (a, b) <- Map.toList ab ]

-- | The number of edges in a matching.
-- Complexity: /O(1)/ time.
--
-- @
-- matchingSize ('matching' [])                 == 0
-- matchingSize ('matching' [(2,\'a\'), (1,\'b\')]) == 2
-- matchingSize ('matching' [(1,\'a\'), (1,\'b\')]) == 1
-- matchingSize ('matching' xs)                 <= 'length' xs
-- matchingSize                               == Map.'Data.Map.Strict.size' . 'pairOfLeft'
-- @
matchingSize :: Matching a b -> Int
matchingSize = Map.size . pairOfLeft

-- | Find a /maximum matching/ in a bipartite graph. A matching is maximum if it
-- has the largest possible size.
-- Complexity: /O(m * sqrt(n) * log(n))/ time.
--
-- @
-- maxMatching 'empty'                                          == 'matching' []
-- maxMatching ('vertices' xs ys)                               == 'matching' []
-- maxMatching ('path' [1,2,3,4])                               == 'matching' [(1,2), (3,4)]
-- 'matchingSize' (maxMatching ('circuit' [(1,2), (3,4), (5,6)])) == 3
-- 'matchingSize' (maxMatching ('star' x (y:ys)))                 == 1
-- 'matchingSize' (maxMatching ('biclique' xs ys))                == 'min' ('length' ('Data.List.nub' xs)) ('length' ('Data.List.nub' ys))
-- 'isMatchingOf' (maxMatching x) x                             == True
-- @
maxMatching :: (Ord a, Ord b) => AdjacencyMap a b -> Matching a b
maxMatching graph = runST (maxMatchingHK graph)

-- TODO: Should we use a more efficient data structure for the queue?
-- TODO: We could try speeding this up by representing vertices with 'Int's.
-- The state maintained by the Hopcroft-Karp algorithm implemented below
data HKState s a b = HKState
    { distance    :: STRef s (Map a Int)
    , curMatching :: STRef s (Matching a b)
    , queue       :: STRef s (Seq a)
    , visited     :: STRef s (Set a) }

-- See https://en.wikipedia.org/wiki/Hopcroft-Karp_algorithm
maxMatchingHK :: forall a b s. (Ord a, Ord b) => AdjacencyMap a b -> ST s (Matching a b)
maxMatchingHK g = do
    distance    <- newSTRef Map.empty
    curMatching <- newSTRef (Matching Map.empty Map.empty)
    queue       <- newSTRef Seq.empty
    visited     <- newSTRef Set.empty
    runHK (HKState distance curMatching queue visited)
    readSTRef curMatching
  where
    runHK :: HKState s a b -> ST s ()
    runHK state = do writeSTRef (distance state) Map.empty
                     foundAugmentingPath <- bfs state
                     when foundAugmentingPath $ do
                         writeSTRef (visited state) Set.empty
                         dfs state
                         runHK state

    currentlyUncovered :: HKState s a b -> ST s [a]
    currentlyUncovered state = do
        m <- readSTRef (curMatching state)
        return [ v | v <- leftVertexList g, not (leftCovered v m) ]


    bfs :: HKState s a b -> ST s Bool
    bfs state = do
        uncovered <- currentlyUncovered state
        mapM_ (enqueue state 1) uncovered
        bfsLoop state

    enqueue :: HKState s a b -> Int -> a -> ST s ()
    enqueue state d v = do modifySTRef (distance state) (Map.insert v d)
                           modifySTRef (queue    state) (|> v)

    dequeue :: HKState s a b -> ST s (Maybe a)
    dequeue state = do q <- readSTRef (queue state)
                       case Seq.viewl q of
                           a :< q -> writeSTRef (queue state) q $> Just a
                           EmptyL -> return Nothing

    bfsLoop :: HKState s a b -> ST s Bool
    bfsLoop state = dequeue state >>= \case
                        Just v  -> do p <- bfsVertex state v
                                      q <- bfsLoop state
                                      return (p || q)
                        Nothing -> return False

    bfsVertex :: HKState s a b -> a -> ST s Bool
    bfsVertex state v = do dist <- readSTRef (distance state)
                           let d = fromJust (v `Map.lookup` dist) + 1
                           or <$> mapM (bfsEdge state d) (neighbours v)

    checkEnqueue :: HKState s a b -> Int -> a -> ST s ()
    checkEnqueue state d v = do dist <- readSTRef (distance state)
                                when (v `Map.notMember` dist) (enqueue state d v)

    bfsEdge :: HKState s a b -> Int -> b -> ST s Bool
    bfsEdge state d u = do m <- readSTRef (curMatching state)
                           case u `Map.lookup` pairOfRight m of
                               Just v  -> checkEnqueue state d v $> False
                               Nothing -> return True

    dfs :: HKState s a b -> ST s ()
    dfs state = currentlyUncovered state >>= mapM_ (dfsVertex state 0)

    dfsVertex :: HKState s a b -> Int -> a -> ST s Bool
    dfsVertex state d v = do dist <- readSTRef (distance state)
                             vis  <- readSTRef (visited state)
                             let dv = fromJust (v `Map.lookup` dist)
                             case (d + 1 == dv) && (v `Set.notMember` vis) of
                                 False -> return False
                                 True  -> do modifySTRef (visited state) (Set.insert v)
                                             dfsEdges state dv v (neighbours v)

    dfsEdges :: HKState s a b -> Int -> a -> [b] -> ST s Bool
    dfsEdges _     _ _ []     = return False
    dfsEdges state d a (b:bs) = do m <- readSTRef (curMatching state)
                                   case b `Map.lookup` pairOfRight m of
                                       Nothing -> addEdge state a b $> True
                                       Just w  -> dfsVertex state d w >>= \case
                                            True  -> addEdge state a b $> True
                                            False -> dfsEdges state d a bs

    addEdge :: HKState s a b -> a -> b -> ST s ()
    addEdge state a b = modifySTRef (curMatching state) (addEdgeUnsafe a b)

    neighbours :: a -> [b]
    neighbours a = Set.toAscList $ fromJust $ Map.lookup a $ leftAdjacencyMap g

-- | A /vertex cover/ of a bipartite graph.
--
-- A /vertex cover/ is a subset of vertices such that every edge is incident to
-- some vertex in the subset. We represent vertex covers by storing two sets of
-- vertices, one for each part. An equivalent representation, which is slightly
-- less memory efficient, is @Set@ @(Either@ @a@ @b)@.
type VertexCover a b = (Set a, Set b)

-- | Check if a given pair of sets is a /vertex cover/ of a bipartite graph.
-- Complexity: /O(m * log(n))/.
--
-- @
-- isVertexCoverOf (xs             , ys             ) 'empty'          == Set.'Set.null' xs && Set.'Set.null' ys
-- isVertexCoverOf (xs             , ys             ) ('leftVertex' x) == Set.'Set.isSubsetOf' xs (Set.'Set.singleton' x) && Set.'Set.null' ys
-- isVertexCoverOf (Set.'Set.empty'      , Set.'Set.empty'      ) ('edge' x y)     == False
-- isVertexCoverOf (Set.'Set.singleton' x, ys             ) ('edge' x y)     == Set.'Set.isSubsetOf' ys (Set.'Set.singleton' y)
-- isVertexCoverOf (xs             , Set.'Set.singleton' y) ('edge' x y)     == Set.'Set.isSubsetOf' xs (Set.'Set.singleton' x)
-- @
isVertexCoverOf :: (Ord a, Ord b) => (Set a, Set b) -> AdjacencyMap a b -> Bool
isVertexCoverOf (as, bs) g = as `Set.isSubsetOf` leftVertexSet g
    && bs `Set.isSubsetOf` rightVertexSet g
    && and [ a `Set.member` as || b `Set.member` bs | (a, b) <- edgeList g ]

-- | The number of vertices in a vertex cover.
-- Complexity: /O(1)/ time.
vertexCoverSize :: VertexCover a b -> Int
vertexCoverSize (as, bs) = Set.size as + Set.size bs

-- | Find a /minimum vertex cover/ in a bipartite graph. A vertex cover is
-- minimum if it has the smallest possible size.
-- Complexity: /O(m * sqrt(n) * log(n))/.
--
-- @
-- minVertexCover 'empty'                              == (Set.'Set.empty', Set.'Set.empty')
-- minVertexCover ('vertices' xs ys)                   == (Set.'Set.empty', Set.'Set.empty')
-- minVertexCover ('path' [1,2,3])                     == (Set.'Set.empty', Set.'Set.singleton' 2)
-- minVertexCover ('star' x (1:2:ys))                  == (Set.'Set.singleton' x, Set.'Set.empty')
-- 'vertexCoverSize' (minVertexCover ('biclique' xs ys)) == 'min' ('length' ('Data.List.nub' xs)) ('length' ('Data.List.nub' ys))
-- 'vertexCoverSize' . minVertexCover                  == 'matchingSize' . 'maxMatching'
-- 'isVertexCoverOf' (minVertexCover x) x              == True
-- @
minVertexCover :: (Ord a, Ord b) => AdjacencyMap a b -> VertexCover a b
minVertexCover g = fromLeft panic $ augmentingPath (maxMatching g) g
  where
    panic = error "minVertexCover: internal error (found augmenting path)"

-- | An /independent set/ of a bipartite graph.
--
-- An /independent set/ is a subset of vertices such that no two of them are
-- adjacent. We represent independent sets by storing two sets of vertices, one
-- for each part. An equivalent representation, which is slightly less memory
-- efficient, is @Set@ @(Either@ @a@ @b)@.
type IndependentSet a b = (Set a, Set b)

-- | Check if a given pair of sets is an /independent set/ of a bipartite graph.
-- Complexity: /O(m * log(n))/.
--
-- @
-- isIndependentSetOf (xs             , ys             ) 'empty'          == Set.'Set.null' xs && Set.'Set.null' ys
-- isIndependentSetOf (xs             , ys             ) ('leftVertex' x) == Set.'Set.isSubsetOf' xs (Set.'Set.singleton' x) && Set.'Set.null' ys
-- isIndependentSetOf (Set.'Set.empty'      , Set.'Set.empty'      ) ('edge' x y)     == True
-- isIndependentSetOf (Set.'Set.singleton' x, ys             ) ('edge' x y)     == Set.'Set.null' ys
-- isIndependentSetOf (xs             , Set.'Set.singleton' y) ('edge' x y)     == Set.'Set.null' xs
-- @
isIndependentSetOf :: (Ord a, Ord b) => (Set a, Set b) -> AdjacencyMap a b -> Bool
isIndependentSetOf (as, bs) g = as `Set.isSubsetOf` leftVertexSet g
    && bs `Set.isSubsetOf` rightVertexSet g
    && and [ not (a `Set.member` as && b `Set.member` bs) | (a, b) <- edgeList g ]

-- | The number of vertices in an independent set.
-- Complexity: /O(1)/ time.
independentSetSize :: IndependentSet a b -> Int
independentSetSize (as, bs) = Set.size as + Set.size bs

-- | Find a /maximum independent set/ in a bipartite graph. An independent set
-- is maximum if it has the largest possible size.
-- Complexity: /O(m * sqrt(n) * log(n))/.
--
-- @
-- maxIndependentSet 'empty'                                 == (Set.'Set.empty', Set.'Set.empty')
-- maxIndependentSet ('vertices' xs ys)                      == (Set.'Set.fromList' xs, Set.'Set.fromList' ys)
-- maxIndependentSet ('path' [1,2,3])                        == (Set.'Set.fromList' [1,3], Set.'Set.empty')
-- maxIndependentSet ('star' x (1:2:ys))                     == (Set.'Set.empty', Set.'Set.fromList' (1:2:ys))
-- 'independentSetSize' (maxIndependentSet ('biclique' xs ys)) == 'max' ('length' ('Data.List.nub' xs)) ('length' ('Data.List.nub' ys))
-- 'independentSetSize' (maxIndependentSet x)                == 'vertexCount' x - 'vertexCoverSize' ('minVertexCover' x)
-- 'isIndependentSetOf' (maxIndependentSet x) x              == True
-- @
maxIndependentSet :: (Ord a, Ord b) => AdjacencyMap a b -> IndependentSet a b
maxIndependentSet g =
    (leftVertexSet g `Set.difference` as, rightVertexSet g `Set.difference` bs)
  where
    (as, bs) = minVertexCover g

-- | Given a matching in a bipartite graph, find either a /vertex cover/ of the
-- same size or an /augmenting path/ with respect to the matching, thereby
-- demonstrating that the matching is not maximum.
-- Complexity: /O((m + n) * log(n))/.
--
-- An /alternating path/ is a path whose edges belong alternately to the
-- matching and not to the matching. An /augmenting path/ is an alternating path
-- that starts from and ends on the vertices that are not covered by the
-- matching. A matching is maximum if and only if there is no augmenting path
-- with respect to it.
--
-- @
-- augmentingPath ('matching' [])      'empty'            == Left (Set.'Set.empty', Set.'Set.empty')
-- augmentingPath ('matching' [])      ('edge' 1 2)       == Right [1,2]
-- augmentingPath ('matching' [(1,2)]) ('path' [1,2,3])   == Left (Set.'Set.empty', Set.'Set.singleton' 2)
-- augmentingPath ('matching' [(3,2)]) ('path' [1,2,3,4]) == Right [1,2,3,4]
-- isLeft (augmentingPath ('maxMatching' x) x)          == True
-- @
augmentingPath :: (Ord a, Ord b) => Matching a b -> AdjacencyMap a b -> Either (VertexCover a b) (List a b)
augmentingPath = augmentingPathImpl

type AugPathMonad a b = MaybeT (State (VertexCover a b)) (List a b)

-- The implementation is in a separate function to avoid the "forall" in docs.
augmentingPathImpl :: forall a b. (Ord a, Ord b) => Matching a b -> AdjacencyMap a b -> Either (VertexCover a b) (List a b)
augmentingPathImpl m g = case runState (runMaybeT dfs) (leftVertexSet g, Set.empty) of
    (Nothing  , cover) -> Left cover
    (Just path, _    ) -> Right path
  where
    dfs :: AugPathMonad a b
    dfs = asum [ inVertex v | v <- leftVertexList g, not (leftCovered v m) ]

    inVertex :: a -> AugPathMonad a b
    inVertex a = do (as, bs) <- lift get
                    guard (a `Set.member` as)
                    lift $ put (Set.delete a as, bs)
                    asum [ onEdge a b | b <- neighbours a ]

    onEdge :: a -> b -> AugPathMonad a b
    onEdge a b = addEdge a b <$> do (as, bs) <- lift get
                                    lift $ put (as, Set.insert b bs)
                                    case b `Map.lookup` pairOfRight m of
                                        Just a  -> inVertex a
                                        Nothing -> return Nil

    addEdge :: a -> b -> List a b -> List a b
    addEdge a b = Cons a . Cons b

    neighbours :: a -> [b]
    neighbours a = Set.toAscList $ fromJust $ Map.lookup a $ leftAdjacencyMap g

-- | Check if the internal representation of a matching is consistent, i.e. that
-- every edge that is present in 'pairOfLeft' is also present in 'pairOfRight'.
-- Complexity: /O(S * log(S))/, where /S/ is the size of the matching.
--
-- @
-- consistentMatching ('matching' xs)   == True
-- consistentMatching ('maxMatching' x) == True
-- @
consistentMatching :: (Ord a, Ord b) => Matching a b -> Bool
consistentMatching (Matching ab ba) =
    Map.toAscList ab == sort [ (a, b) | (b, a) <- Map.toAscList ba ]