algebraic-graphs-0.6.1: src/Algebra/Graph/Bipartite/AdjacencyMap/Algorithm.hs
{-# 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 ]