Rewrite Bellman-Ford and min-cost flow, especially to stop the latter from crashing.
authorMatt McCutchen <matt@mattmccutchen.net>
Fri, 11 Jul 2008 03:46:07 +0000 (23:46 -0400)
committerMatt McCutchen <matt@mattmccutchen.net>
Fri, 11 Jul 2008 03:46:07 +0000 (23:46 -0400)
- 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 [new file with mode: 0644]
program/BellmanFord.hs
program/GraphStuff.hs.old [new file with mode: 0644]
program/Instance.hs
program/MonadStuff.hs [new file with mode: 0644]
program/NaiveMinCostFlow.hs [new file with mode: 0644]
program/ProposalMatcher.hs
program/Test.hs
program/UnitMinCostFlow.hs [deleted file]

diff --git a/program/ArrayQueue.hs b/program/ArrayQueue.hs
new file mode 100644 (file)
index 0000000..252b5fb
--- /dev/null
@@ -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
index 506acd9..da0fdf8 100644 (file)
-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 (file)
index 0000000..341090b
--- /dev/null
@@ -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
index 95e00ee..07841e3 100644 (file)
@@ -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 (file)
index 0000000..904fa83
--- /dev/null
@@ -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 (file)
index 0000000..1347837
--- /dev/null
@@ -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)
+       )
index c0bd7a0..f523c27 100644 (file)
@@ -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
index f8cb4ef..260ce0b 100644 (file)
@@ -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 (file)
index 24e44bc..0000000
+++ /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