From 5a07db44406bad03321a90b0814cc4496c6b7d63 Mon Sep 17 00:00:00 2001 From: Matt McCutchen Date: Thu, 10 Jul 2008 23:46:07 -0400 Subject: [PATCH] Rewrite Bellman-Ford and min-cost flow, especially to stop the latter from crashing. - Rewrite Bellman-Ford and min-cost flow solvers using the state-transformer feature, adding infrastructure as necessary. - Make Bellman-Ford detect the presence of a negative cycle. - Min-cost flow now uses edges of arbitrary capacity, not just 1, although it's the same naive algorithm that augments along the lowest-cost path no matter how little flow that adds. Adjust the reduction accordingly; now it uses fewer edges. Additional benefit: the change brings my min-cost flower's interface closer to that of an eventual CS2 adaptor. - Min-cost flow was crashing on larger instances because round-off error prevented it from identifying the edge to be flipped by subtracting the costs to its endpoints. To avoid this problem, redesign min-cost flow so edges are identified by unique indices. Now doMatching "pulls" flow values for the edgesD it generated instead of being "pushed" by edges in the residual graph. - Make the Show implementation for Instances prettier. --- program/ArrayQueue.hs | 53 ++++++++++++ program/BellmanFord.hs | 164 +++++++++++++++++++++++------------- program/GraphStuff.hs.old | 75 +++++++++++++++++ program/Instance.hs | 10 ++- program/MonadStuff.hs | 5 ++ program/NaiveMinCostFlow.hs | 100 ++++++++++++++++++++++ program/ProposalMatcher.hs | 94 ++++++++++++++------- program/Test.hs | 26 ++++-- program/UnitMinCostFlow.hs | 60 ------------- 9 files changed, 425 insertions(+), 162 deletions(-) create mode 100644 program/ArrayQueue.hs create mode 100644 program/GraphStuff.hs.old create mode 100644 program/MonadStuff.hs create mode 100644 program/NaiveMinCostFlow.hs delete mode 100644 program/UnitMinCostFlow.hs diff --git a/program/ArrayQueue.hs b/program/ArrayQueue.hs new file mode 100644 index 0000000..252b5fb --- /dev/null +++ b/program/ArrayQueue.hs @@ -0,0 +1,53 @@ +module ArrayQueue ( + ArrayQueue, newArrayQueue, + aqEnqueue, aqIsEmpty, aqDequeue +) where +import Control.Monad.ST +import Data.Array.ST +import MonadStuff + +data ArrayQueue s = ArrayQueue { + aqArr :: STUArray s Int Int, + aqHeadI :: Int, -- Also used as an "end of queue" sentinel + aqTailI :: Int -- Element value can be aqHeadI; also used as a "not queued" sentinel +} + +newArrayQueue :: (Int, Int) -> ST s (ArrayQueue s) +newArrayQueue (lo, hi) = do + let headI = lo - 1 + let tailI = lo - 2 + arr <- newArray (tailI, hi) tailI + writeArray arr headI headI -- queue is empty + writeArray arr tailI headI -- tail is head + return $ ArrayQueue arr headI tailI + +aqEnqueue :: ArrayQueue s -> Int -> ST s Bool -- Was it added? +aqEnqueue (ArrayQueue arr headI tailI) newI = do + newIval <- readArray arr newI + if newIval == tailI + then do + lst <- readArray arr tailI + writeArray arr lst newI -- Append newI. + writeArray arr newI headI -- newI is now the tail. + writeArray arr tailI newI -- The tail is now at newI. + return True + else return False -- Already on queue. + +aqIsEmpty :: ArrayQueue s -> ST s Bool +aqIsEmpty (ArrayQueue arr headI tailI) = do + first <- readArray arr headI + return (first == headI) + +aqDequeue :: ArrayQueue s -> ST s (Maybe Int) +aqDequeue (ArrayQueue arr headI tailI) = do + first <- readArray arr headI + if first == headI + then return Nothing + else do + next <- readArray arr first + writeArray arr headI next + if next == headI + then writeArray arr tailI headI -- Emptied queue. + else nop + writeArray arr first tailI -- No longer on queue. + return $ Just first diff --git a/program/BellmanFord.hs b/program/BellmanFord.hs index 506acd9..da0fdf8 100644 --- a/program/BellmanFord.hs +++ b/program/BellmanFord.hs @@ -1,69 +1,113 @@ -module BellmanFord {-(spTree)-} where +module BellmanFord (bellmanFord, BFPath(BFPath)) where import Data.Graph.Inductive.Graph import Data.Graph.Inductive.Internal.Queue -import Data.Graph.Inductive.Internal.RootPath -import Data.Array.Diff +import Data.Array.IArray +import Data.Array.ST +import Data.STRef +import Control.Monad.ST +import ArrayQueue +import MonadStuff +import Data.List -data NodeInfo b = NodeInfo { - path :: Maybe [LNode b], - changed :: Bool -} -data Graph gr => BFState gr a b = BFState { - theGraph :: gr a b, - nis :: DiffArray Node (NodeInfo b), - changedQ :: Queue Node +-- Path to a node +data Real w => BFPath b w = BFPath { + pLen :: w, -- Total distance at the end of the path + pDest :: Node, -- Destination of the path + pFrom :: Maybe (b, BFPath b w) -- Last edge and remaining path (Nothing for source) +} deriving Show + +data (Graph gr, Real w) => BFState s gr a b w = BFState { + bfsGraph :: gr a b, + bfsEdgeWt :: b -> w, + bfsMPaths :: STArray s Node (Maybe (BFPath b w)), + bfsChanged :: ArrayQueue s, + -- Sentinel in the queue that ends a pass. + bfsPassEnder :: Node, + -- The length of the paths we are currently considering. + -- If this reaches the number of nodes in the graph, we must have a negative cycle. + bfsPass :: STRef s Int } -nisToLRTree nis = do - ni <- elems nis - case path ni of - Just lnl -> return (LP lnl) - Nothing -> fail "Node is unreachable" +{-- +negativeCycleCheck :: (Graph gr, Real w) => BFState s gr a b w -> + BFPath b w -> ST s () +negativeCycleCheck state path = + let getNodes (BFPath _ dst from) = dst : case from of + Nothing -> [] + Just (_, p1) -> getNodes p1 in + let nodes = getNodes path in + if length (nub nodes) < length nodes + then error ("Negative cycle detected: " ++ show path) + else nop +--} + +offerPath :: (Graph gr, Real w) => BFState s gr a b w -> + BFPath b w -> ST s () +offerPath state newPath@(BFPath newLen dest _) = do + oldMPath <- readArray (bfsMPaths state) dest + -- Is newPath the first path to dest or shorter than the existing one? + let adoptPath = case oldMPath of + Nothing -> True + Just (BFPath oldLen _ _) -> newLen < oldLen + if adoptPath + then do + -- Save the new path. + writeArray (bfsMPaths state) dest (Just newPath) + -- Mark dest as changed. + aqEnqueue (bfsChanged state) dest + nop + else nop -- Don't update anything. -offerPath :: (Graph gr, Real b) => BFState gr a b -> [LNode b] -> BFState gr a b -offerPath bfs newPath@((dest, newDist): _) = - -- Is newPath the first path to dest, or better than the previous one? - let adoptPath = - case path (nis bfs ! dest) of - Nothing -> True - Just ((_, oldDist) : _) -> newDist < oldDist - in - if adoptPath then - bfs{ - -- Update NodeInfo with the new path. - nis = nis bfs // [(dest, NodeInfo (Just newPath) True)], - changedQ = if changed (nis bfs ! dest) - then changedQ bfs -- Already on the queue; leave as is. - else queuePut dest (changedQ bfs) -- Add to queue. - } - else bfs -- Don't update anything. +processEdge :: (Graph gr, Real w) => BFState s gr a b w -> + BFPath b w -> LEdge b -> ST s () +processEdge state path1@(BFPath len1 _ _) (_, dest, label) = do + let edgeLen = (bfsEdgeWt state) label + let newPath = BFPath (len1 + edgeLen) dest (Just (label, path1)) + offerPath state newPath -processEdge :: (Graph gr, Real b) => [LNode b] -> LEdge b -> BFState gr a b -> BFState gr a b -processEdge srcPath@((_, srcDist) : _) (_, dest, edgeLen) bfs = - let newPath = (dest, srcDist + edgeLen) : srcPath in - offerPath bfs newPath +endPass :: (Graph gr, Real w) => BFState s gr a b w -> ST s () +endPass state = do + qEmpty <- aqIsEmpty (bfsChanged state) + if qEmpty + then nop -- No nodes to visit on the next pass. We're done. + else do + -- Increment the pass number. + modifySTRef (bfsPass state) (+ 1) + p <- readSTRef (bfsPass state) + if p == noNodes (bfsGraph state) + then error "BellmanFord: Negative cycle detected!" + else do + -- Re-enqueue the pass ender. + aqEnqueue (bfsChanged state) (bfsPassEnder state) + search state -- Continue with the next pass. -search :: (Graph gr, Real b) => BFState gr a b -> LRTree b -search bfs = - if queueEmpty (changedQ bfs) then - -- Finished. - nisToLRTree (nis bfs) - else - -- Process a changed node from the queue. - let (src, moreQ) = queueGet (changedQ bfs) in - let srcNI = nis bfs ! src in - -- Clear src's changed flag. - let bfs1 = bfs{nis = nis bfs // [(src, srcNI{changed = False})], changedQ = moreQ} in - let Just srcPath = path srcNI in - let outEdges = out (theGraph bfs) src in - let newBFS = foldr (processEdge srcPath) bfs1 outEdges in - search newBFS +search :: forall s gr a b w. (Graph gr, Real w) => BFState s gr a b w -> ST s () +search state = do + Just node <- aqDequeue (bfsChanged state) + if node == bfsPassEnder state + then endPass state + else do + mPath1 <- readArray (bfsMPaths state) node + let Just path1 = mPath1 + sequence $ map (processEdge state path1) $ + out (bfsGraph state) node + search state -- Keep going. -spTree :: (Graph gr, Real b) => Node -> gr a b -> LRTree b -spTree source theGraph = let theNodes = range (nodeRange theGraph) in - let emptyBFS = BFState theGraph - (array (nodeRange theGraph) (map (\n -> (n, NodeInfo Nothing False)) theNodes)) - mkQueue in - -- Start with a zero-length path to the source. - let initBFS = offerPath emptyBFS [(source, 0)] in - search initBFS +bellmanFord :: (Graph gr, Real w) => + (b -> w) -> -- Edge label -> weight + Node -> -- Source node + gr a b -> -- Graph + Array Node (Maybe (BFPath b w)) -- ! node -> maybe a path +bellmanFord edgeWt source theGraph = runSTArray (do + mPaths <- newArray (nodeRange theGraph) Nothing + let (nlo, nhi) = nodeRange theGraph + let passEnder = nlo - 1 + changed <- newArrayQueue (passEnder, nhi) -- Queue can contain the pass ender. + pass <- newSTRef 0 + let state = BFState theGraph edgeWt mPaths changed passEnder pass + -- Pass is 0, and the queue contains a node that was offered a path with 0 edges. + offerPath state (BFPath 0 source Nothing) + aqEnqueue (bfsChanged state) (bfsPassEnder state) + search state + return (bfsMPaths state) + ) diff --git a/program/GraphStuff.hs.old b/program/GraphStuff.hs.old new file mode 100644 index 0000000..341090b --- /dev/null +++ b/program/GraphStuff.hs.old @@ -0,0 +1,75 @@ +-- I wrote this to help with NaiveMinCostFlow but abandoned it. Wrapping is +-- more trouble than it's worth, and NaiveMinCostFlow needs access to its ST +-- monad in the edge-filtering test, which would have required more mess here. + +module GraphStuff where +import Data.Graph.Inductive.Graph + +data EdgeMapper b1 b2 = EdgeMapper { + -- Map the predecessor and successor edges of the given node for a context. + mapAdj :: Node -> (Adj b1, Adj b1) -> (Adj b2, Adj b2), + -- For labEdges. Should be consistent with the above two!!! + mapLEdge :: LEdge b1 -> [LEdge b2] +} + +data EdgeMappedGraph a b = forall gr b1. Graph gr => EdgeMappedGraph { + emgMapper :: EdgeMapper b1 b, + emgOrig :: gr a b1 +} + +edgeMapContext :: EdgeMapper b1 b2 -> Context a b1 -> Context a b2 +edgeMapContext mapper (p, n, nl, s) = + let (p2, s2) = (mapAdj mapper) n (p, s) in (p2, n, nl, s2) + +instance Graph EdgeMappedGraph where + isEmpty g = null (labEdges g) + labNodes (EdgeMappedGraph mapper orig) = labNodes orig + noNodes (EdgeMappedGraph mapper orig) = noNodes orig + nodeRange (EdgeMappedGraph mapper orig) = nodeRange orig + labEdges (EdgeMappedGraph mapper orig) = concatMap (mapLEdge mapper) (labEdges orig) + + match n (EdgeMappedGraph mapper orig) = + let + (mCtx, g1) = match n orig + mCtx2 = do + ctx <- mCtx + return (edgeMapContext mapper ctx) + in (mCtx2, EdgeMappedGraph mapper g1) + + matchAny (EdgeMappedGraph mapper orig) = + let (ctx, g1) = matchAny orig in + (edgeMapContext mapper ctx, EdgeMappedGraph mapper g1) + + -- Graph construction: not supported. + empty = undefined + mkGraph nodes edges = undefined + +buildEdgeMappedGraph :: Graph gr => EdgeMapper b1 b2 -> gr a b1 -> gr a b2 +buildEdgeMappedGraph mapper g = + mkGraph (labNodes g) (concatMap (mapLEdge mapper) $ labEdges g) + +edgeFilterMapper :: (LEdge b -> Bool) -> EdgeMapper b b +edgeFilterMapper ff = EdgeMapper + (\n (p, s) -> (filter (\(el, pn) -> ff (pn, n, el)) p, + filter (\(el, sn) -> ff ( n, sn, el)) s)) + (\edge -> if ff edge then [edge] else []) + +-- Returns a wrapper (whose type does not support graph construction) instead +-- of building a new graph of the original type. +filterEdgesLite :: Graph gr => (LEdge b -> Bool) -> gr a b -> EdgeMappedGraph a b +filterEdgesLite filterF graph = EdgeMappedGraph (edgeFilterMapper filterF) graph + +edgeBidirectorMapper :: (Bool -> b1 -> b2) -> EdgeMapper b1 b2 +edgeBidirectorMapper flipF = EdgeMapper + (\n (p, s) -> ( + map (\(el, pn) -> (flipF False el, pn)) p ++ + map (\(el, sn) -> (flipF True el, sn)) s, + map (\(el, sn) -> (flipF False el, sn)) s ++ + map (\(el, pn) -> (flipF True el, pn)) p)) + (\(pn, sn, el) -> [(pn, sn, flipF False el), (sn, pn, flipF True el)]) + +bidirectEdgesLite :: Graph gr => (Bool -> b1 -> b2) -> gr a b1 -> EdgeMappedGraph a b2 +bidirectEdgesLite flipF graph = EdgeMappedGraph (edgeBidirectorMapper flipF) graph + +buildBidirectedEdgeGraph :: Graph gr => (Bool -> b1 -> b2) -> gr a b1 -> gr a b2 +buildBidirectedEdgeGraph flipF graph = buildEdgeMappedGraph (edgeBidirectorMapper flipF) graph diff --git a/program/Instance.hs b/program/Instance.hs index 95e00ee..07841e3 100644 --- a/program/Instance.hs +++ b/program/Instance.hs @@ -15,7 +15,9 @@ data Instance = Instance instance Show Instance where show (Instance numRvrs numProps loadA prefA) = - "Instance: " ++ show numRvrs ++ " reviewers, " ++ show numProps ++ " proposals\n" - ++ "Reviewer relative load: " ++ show loadA ++ "\n" - ++ "Preferences:\n" - ++ formatTable (array2DtoListOfLists (amap2 show prefA :: Array (Int, Int) String)) + let theRvrs = [0..numRvrs-1]; theProps = [0..numProps-1] in + "Instance with " ++ show numRvrs ++ " reviewers and " ++ show numProps ++ " proposals:\n" ++ formatTable ( + ( "" : map (\i -> "R#" ++ show i ) theRvrs) : + ( "RLoad" : map (\i -> show (loadA ! i) ) theRvrs) : + map (\j -> ("P#" ++ show j) : map (\i -> show (prefA ! (i, j))) theRvrs) theProps + ) diff --git a/program/MonadStuff.hs b/program/MonadStuff.hs new file mode 100644 index 0000000..904fa83 --- /dev/null +++ b/program/MonadStuff.hs @@ -0,0 +1,5 @@ +module MonadStuff (nop) where + +-- Useful for "if this, then mutate that, else do nothing" +nop :: Monad m => m () +nop = return () diff --git a/program/NaiveMinCostFlow.hs b/program/NaiveMinCostFlow.hs new file mode 100644 index 0000000..1347837 --- /dev/null +++ b/program/NaiveMinCostFlow.hs @@ -0,0 +1,100 @@ +module NaiveMinCostFlow (minCostFlow) where +import BellmanFord +import MonadStuff +import Data.Array.IArray +import Data.Array.ST +import Control.Monad +import Control.Monad.ST +import Data.Graph.Inductive.Graph +import Data.Graph.Inductive.Internal.RootPath +import Data.List + +data MCFEdge i f c = MCFEdge { + edgeIdx :: i, + edgeCap :: f, + edgeCost :: c, + edgeIsRev :: Bool +} +data MCFState s gr a b i f c = MCFState { + mcfGraph :: gr a (MCFEdge i f c), + mcfSource :: Node, + mcfSink :: Node, + mcfFlow :: STArray s i f +} + +edgeCapLeft :: (Graph gr, Ix i, Real f, Real c) => MCFState s gr a b i f c -> + MCFEdge i f c -> ST s f +edgeCapLeft state (MCFEdge i cap _ isRev) = do + fwdFlow <- readArray (mcfFlow state) i + return (if isRev then fwdFlow else cap - fwdFlow) + +edgePush :: (Graph gr, Ix i, Real f, Real c) => MCFState s gr a b i f c -> + MCFEdge i f c -> f -> ST s () +edgePush state (MCFEdge i _ _ isRev) nf = do + oldFwdFlow <- readArray (mcfFlow state) i + let newFwdFlow = if isRev then oldFwdFlow - nf else oldFwdFlow + nf + writeArray (mcfFlow state) i newFwdFlow + +pathCapLeft :: (Graph gr, Ix i, Real f, Real c) => MCFState s gr a b i f c -> + (MCFEdge i f c, BFPath (MCFEdge i f c) c) -> ST s f +pathCapLeft state (lastEdge, BFPath _ _ mFrom) = do + lastCL <- edgeCapLeft state lastEdge + case mFrom of + Nothing -> return lastCL + Just from -> do + fromCL <- pathCapLeft state from + return (min lastCL fromCL) + +augment :: (Graph gr, Ix i, Real f, Real c) => MCFState s gr a b i f c -> + f -> BFPath (MCFEdge i f c) c -> ST s () +augment state augAmt (BFPath _ _ mFrom) = case mFrom of + Nothing -> nop + Just (lastEdge, path1) -> do + edgePush state lastEdge augAmt + augment state augAmt path1 + +doFlow :: forall s gr a b i f c. (Graph gr, Ix i, Real f, Real c) => MCFState s gr a b i f c -> + ST s () +doFlow state = do + filteredEdges <- filterM (\(_, _, l) -> do + ecl <- edgeCapLeft state l + return (ecl /= 0) + ) (labEdges (mcfGraph state)) + let filteredGraph = mkGraph (labNodes (mcfGraph state)) filteredEdges :: gr a (MCFEdge i f c) + -- Why won't we get a negative cycle? The original graph from the + -- proposal matcher is acyclic (so no negative cycle), and if we + -- created a negative cycle, that would contradict the fact that we + -- always augment along the lowest-cost path. + let mAugPath = bellmanFord edgeCost (mcfSource state) filteredGraph + ! (mcfSink state) + case mAugPath of + Nothing -> nop -- Done. + -- source /= sink, so augPasth is not empty. + Just augPath@(BFPath _ _ (Just from)) -> do + augAmt <- pathCapLeft state from + augment state augAmt augPath + doFlow state + +minCostFlow :: forall s gr a b i f c. (Graph gr, Ix i, Real f, Real c) => + (i, i) -> -- Range of edge indices + (b -> i) -> -- Edge label -> unique edge index + (b -> f) -> -- Edge label -> flow capacity + (b -> c) -> -- Edge label -> cost per unit of flow + gr a b -> -- Graph + (Node, Node) -> -- (source, sink) + Array i f -- ! edge index -> flow value +minCostFlow idxBounds edgeIdx edgeCap edgeCost theGraph (source, sink) = runSTArray (do + let ourFlipF isRev l = + MCFEdge (edgeIdx l) (edgeCap l) + (if isRev then -(edgeCost l) else edgeCost l) + isRev + let graph2 = mkGraph (labNodes theGraph) (concatMap + (\(n1, n2, l) -> [ -- Capacity of reverse edge is never used. + (n1, n2, MCFEdge (edgeIdx l) (edgeCap l) ( edgeCost l ) False), + (n2, n1, MCFEdge (edgeIdx l) undefined (-(edgeCost l)) True ) + ]) $ labEdges theGraph) :: gr a (MCFEdge i f c) + flow <- newArray idxBounds 0 + let state = MCFState graph2 source sink flow + doFlow state + return (mcfFlow state) + ) diff --git a/program/ProposalMatcher.hs b/program/ProposalMatcher.hs index c0bd7a0..f523c27 100644 --- a/program/ProposalMatcher.hs +++ b/program/ProposalMatcher.hs @@ -1,5 +1,5 @@ module ProposalMatcher where -import UnitMinCostFlow +import NaiveMinCostFlow import Data.Array.IArray import Data.Graph.Inductive.Graph import Data.Graph.Inductive.Tree @@ -13,7 +13,36 @@ prefBoringness p = if prefIsVeryBoring p then 2 prefExpertness p = if prefIsExpert p then 2 else if prefIsKnowledgeable p then 1 else 0 -doReduction :: Instance -> Gr () Wt +data REdge = REdge { + reIdx :: Int, + reCap :: Int, + reCost :: Wt +} + +instance Show REdge where + show (REdge idx cap cost) = "#" ++ (show idx) ++ ": " + ++ (show cap) ++ " @ " ++ (show cost) + +data ReductionResult = ReductionResult { + rrGraph :: Gr () REdge, + rrSource :: Node, + rrSink :: Node, + rrEIdxBounds :: (Int, Int), + rrEDIdx :: (Int, Int) -> Int +} + +-- Hack: show as much of the reduction result as we easily can +data RR1 = RR1 (Gr () REdge) Node Node (Int, Int) deriving Show +instance Show ReductionResult where + show (ReductionResult g so si eib _) = show (RR1 g so si eib) + +indexEdges :: Int -> [(Int, Int, REdge)] -> (Int, [(Int, Int, REdge)]) +indexEdges i [] = (i, []) +indexEdges i ((v1, v2, re):es) = + let (imax, ies) = indexEdges (i+1) es in + (imax, (v1, v2, re{ reIdx = i }) : ies) + +doReduction :: Instance -> ReductionResult doReduction (Instance numRvrs numProps rloadA prefA) = let source = 0 @@ -21,6 +50,7 @@ doReduction (Instance numRvrs numProps rloadA prefA) = rvrNode i boringness = 2 + 3*i + boringness propNode j expertness = 2 + 3*numRvrs + 3*j + expertness numNodes = 2 + 3*numRvrs + 3*numProps + edIdx (i, j) = i*numProps + j in let totalReviews = reviewsEachProposal * numProps @@ -30,16 +60,19 @@ doReduction (Instance numRvrs numProps rloadA prefA) = edgesABC = do i <- [0 .. numRvrs - 1] let tl = targetLoad i - l <- [0 .. tl + loadTolerance - 1] - let costA = if l < tl - then 0 - else marginalLoadCost ((numAsWt (l - tl) + 1/2) / numAsWt loadTolerance) - let edgeA = (source, rvrNode i 0, costA) - let costB = marginalBoringCost ((numAsWt l + 1/2) / numAsWt tl) - let edgeB = (rvrNode i 0, rvrNode i 1, costB) - let costC = marginalVeryBoringCost ((numAsWt l + 1/2) / numAsWt tl) - let edgeC = (rvrNode i 1, rvrNode i 2, costC) - [edgeA, edgeB, edgeC] + let freeEdgeA = (source, rvrNode i 0, REdge undefined tl 0) + let nonfreeEdgesA = do + l <- [tl .. tl + loadTolerance - 1] + let costA = marginalLoadCost ((numAsWt (l - tl) + 1/2) / numAsWt loadTolerance) + [(source, rvrNode i 0, REdge undefined 1 costA)] + let edgesBC = do + l <- [0 .. tl + loadTolerance - 1] + let costB = marginalBoringCost ((numAsWt l + 1/2) / numAsWt tl) + let edgeB = (rvrNode i 0, rvrNode i 1, REdge undefined 1 costB) + let costC = marginalVeryBoringCost ((numAsWt l + 1/2) / numAsWt tl) + let edgeC = (rvrNode i 1, rvrNode i 2, REdge undefined 1 costC) + [edgeB, edgeC] + [freeEdgeA] ++ nonfreeEdgesA ++ edgesBC edgesD = do i <- [0 .. numRvrs - 1] j <- [0 .. numProps - 1] @@ -48,22 +81,22 @@ doReduction (Instance numRvrs numProps rloadA prefA) = then [] else [(rvrNode i (prefBoringness pref), propNode j (prefExpertness pref), - assignmentCost pref)] - edgesE = do - j <- [0 .. numProps - 1] - [(propNode j 2, propNode j 0, -expertBonus)] - edgesFGH = do + REdge (edIdx (i, j)) 1 (assignmentCost pref))] + edgesEFGH = do j <- [0 .. numProps - 1] - l <- [0 .. reviewsEachProposal - 1] - let edgeF = (propNode j 2, propNode j 1, 0) - let edgeG = (propNode j 1, propNode j 0, - if l == 0 then -knowledgeableBonus else 0) - let edgeH = (propNode j 0, sink, 0) - [edgeF, edgeG, edgeH] + let edgeE = (propNode j 2, propNode j 0, REdge undefined 1 (-expertBonus)) + let edgeF = (propNode j 2, propNode j 1, REdge undefined reviewsEachProposal 0) + let edgeGFirst = (propNode j 1, propNode j 0, REdge undefined 1 (-knowledgeableBonus)) + let edgeGRest = (propNode j 1, propNode j 0, REdge undefined (reviewsEachProposal-1) 0) + let edgeH = (propNode j 0, sink, REdge undefined reviewsEachProposal 0) + [edgeE, edgeF, edgeGFirst, edgeGRest, edgeH] theNodes = [(x, ()) | x <- [0 .. numNodes - 1]] - theEdges = edgesABC ++ edgesD ++ edgesE ++ edgesFGH + -- Index the non-D edges + unindexedEdges = edgesABC ++ edgesEFGH + (imax, reindexedEdges) = indexEdges (numRvrs*numProps) unindexedEdges + theEdges = edgesD ++ reindexedEdges in - mkGraph theNodes theEdges + ReductionResult (mkGraph theNodes theEdges) source sink (0, imax-1) edIdx todo = undefined -- Returns a list of reviews as ordered pairs (reviewer#, proposal#). @@ -79,14 +112,13 @@ doMatching inst@(Instance numRvrs numProps _ _) = idPropNode n = (n - (2 + 3*numRvrs)) `div` 3 numNodes = 2 + 3*numRvrs + 3*numProps in - let graph1 = doReduction inst in - let flow1 = flowDiff graph1 (snd (umcf source sink graph1)) in + let ReductionResult graph source sink idxBounds edIdx = doReduction inst in + let flowArray = minCostFlow idxBounds reIdx reCap reCost graph (source, sink) in let pairs = do i <- [0 .. numRvrs - 1] - boringness <- [0, 1, 2] - n <- suc flow1 (rvrNode i boringness) - if n >= firstPropNode - then [(i, idPropNode n)] + j <- [0 .. numProps - 1] + if flowArray ! edIdx (i, j) == 1 + then [(i, j)] else [] in sort pairs -- for prettiness diff --git a/program/Test.hs b/program/Test.hs index f8cb4ef..260ce0b 100644 --- a/program/Test.hs +++ b/program/Test.hs @@ -29,7 +29,7 @@ import Data.Graph.Inductive.Graphviz -- Other imports we need import BellmanFord -import UnitMinCostFlow +import NaiveMinCostFlow import Data.Array.IArray import Data.Array.Unboxed import Data.Graph.Inductive.Graph @@ -37,11 +37,22 @@ import Data.Graph.Inductive.Tree import ArrayStuff myGraph = mkGraph [(0, ()), (1, ()), (2, ())] - [(0, 1, 2), (0, 2, 3), (2, 1, -2)] :: Gr () Double + [(0, 1, (0, 2)), (0, 2, (1, 3)), (2, 1, (2, -2))] :: Gr () (Int, Int) -spTree1 = spTree 0 myGraph +bfResult = bellmanFord snd 0 myGraph -(flowVal, flowResid) = umcf 0 1 myGraph +flowArray = minCostFlow (0, 2) fst (const 1) snd myGraph (0, 1) + +myNCGraph = mkGraph [(0, ())] [(0, 0, -1)] :: Gr () Int +bfNCResult = bellmanFord id 0 myNCGraph + +data REdgeF = REdgeF Int Int Int Wt +instance Show REdgeF where + show (REdgeF idx cap flow cost) = "#" ++ (show idx) ++ ": " + ++ (show flow) ++ " of " ++ (show cap) ++ " @ " ++ (show cost) +flowAnnotate g fa = + mkGraph (labNodes g) (map (\(n1, n2, REdge i ca co) -> + (n1, n2, REdgeF i ca (fa ! i) co)) $ labEdges g) :: Gr () REdgeF -- Example from idea book p. 425 {- @@ -64,7 +75,8 @@ myPrefs = transposeArray $ listArray ((0,0), (myNumProps-1,myNumRvrs-1)) [ myInst = Instance myNumRvrs myNumProps (funcArray (0, myNumRvrs-1) $ const 1) myPrefs -rdnGraph = doReduction myInst -(rdnFlowVal, rdnFlowResid) = umcf 0 1 rdnGraph -rdnFlow = flowDiff rdnGraph rdnFlowResid +rdnResult = doReduction myInst +ReductionResult rrg rrso rrsi rreib rredi = rdnResult +rdnFlowArray = minCostFlow rreib reIdx reCap reCost rrg (rrso, rrsi) +rrg2 = flowAnnotate rrg rdnFlowArray myMatching = doMatching myInst diff --git a/program/UnitMinCostFlow.hs b/program/UnitMinCostFlow.hs deleted file mode 100644 index 24e44bc..0000000 --- a/program/UnitMinCostFlow.hs +++ /dev/null @@ -1,60 +0,0 @@ -module UnitMinCostFlow (umcf, flowDiff) where -import BellmanFord -import Data.Graph.Inductive.Graph -import Data.Graph.Inductive.Internal.RootPath -import Data.List - -maybeDelete :: Eq a => a -> [a] -> Maybe [a] -maybeDelete _ [] = Nothing -maybeDelete e (h:t) = if e == h - then return t - else do t1 <- maybeDelete e t; return (h:t1) - --- If the edge occurs in the graph, return Just the graph with one occurrence --- deleted; otherwise return Nothing. (delLEdge deletes all occurrences.) -maybeDelOneLEdge :: (DynGraph gr, Eq b) => LEdge b -> gr a b -> Maybe (gr a b) -maybeDelOneLEdge (src, dest, lbl) theGraph = - let (mc, moreGraph) = match src theGraph in do - (p, v, l, s) <- mc - s2 <- maybeDelete (lbl, dest) s - return ((p, v, l, s2) & moreGraph) - -flipEdge (src, dest, lbl) = (dest, src, -lbl) - -flipEdgeIn :: (DynGraph gr, Real b) => LEdge b -> gr a b -> gr a b -flipEdgeIn edge theGraph = - let Just graph1 = maybeDelOneLEdge edge theGraph in - insEdge (flipEdge edge) graph1 - -augment :: (DynGraph gr, Real b) => [LNode b] -> gr a b -> gr a b -augment augPath@((v1, d1) : t1) theGraph = - case t1 of - [] -> theGraph - (v2, d2) : t2 -> augment (tail augPath) - (flipEdgeIn (v1, v2, d2 - d1) theGraph) - --- Find a min-cost flow from s to t in theGraph. --- Each edge of the input graph has unit capacity and cost given by its label. --- Returns: flow value, residual graph. -umcf :: (DynGraph gr, Real b) => Node -> Node -> gr a b -> (b, gr a b) -umcf s t theGraph = - -- Use Bellman-Ford to find an augmenting path from s to t, if one exists. - -- NOTE: getLPath reverses it into s-to-t order! - let LP augPath = getLPath t (spTree s theGraph) in - if null augPath then - -- Finished. - (0, theGraph) - else - -- Augment, and continue flowing in the resulting graph. - let graph2 = augment augPath theGraph in - let (fval1, resid) = umcf s t graph2 in (fval1 + 1, resid) - --- Diffs an original graph and a residual graph, producing the flow graph. -flowDiff :: (DynGraph gr, Real b) => gr a b -> gr a b -> gr a b -flowDiff ograph rgraph = case labEdges ograph of - [] -> mkGraph (labNodes ograph) [] - oedge:_ -> let Just ograph2 = maybeDelOneLEdge oedge ograph in - case maybeDelOneLEdge oedge rgraph of - Just rgraph2 -> flowDiff ograph2 rgraph2 - Nothing -> let Just rgraph2 = maybeDelOneLEdge (flipEdge oedge) rgraph in - insEdge oedge (flowDiff ograph2 rgraph2) \ No newline at end of file -- 2.34.1