Max flow: Edmonds-Karp

README.md

@@ -1,14 +1,16 @@
-# Loom
+![Loom logo](https://raw.github.com/aysylu/loom/master/doc/loom_logo.png "Loom")
 
 **Caveat coder**: this lib is alpha-stage. The API may change in future versions.
 
-That said, feedback welcome. Email: [jkkramer@gmail.com](mailto:jkkramer@gmail.com), Twitter: [jkkramer](http://twitter.com/jkkramer)
+This library is based off of Justin Kramer's Loom library.
+
+Contact: email [robertlachlan@gmail.com](mailto:robertlachlan@gmail.com), [aysylu7 [at] gmail [dot] com](mailto:aysylu7@gmail.com), twitter: [aysylu22](http://twitter.com/aysylu22)
 
 ## Usage
 
 ### Leiningen/Clojars [group-id/name version]
 
-    [jkkramer/loom "0.2.0"]
+    [haeffalump/loom "0.3.0-SNAPSHOT"]
 
 ### Namespaces
 
@@ -139,4 +141,4 @@ Nothing but Clojure. There is optional support for visualization via [GrapViz](h
 
 Copyright (C) 2010 Justin Kramer jkkramer@gmail.com
 
-Distributed under the Eclipse Public License, the same as Clojure.
+Distributed under the [Eclipse Public License](http://opensource.org/licenses/eclipse-1.0.php), the same as Clojure.

project.clj

@@ -1,4 +1,7 @@
-(defproject jkkramer/loom "0.2.0"
+(defproject heffalump/loom "0.3.0-SNAPSHOT"
   :description "Graph library for Clojure"
-  :author "Justin Kramer"
-  :dependencies [[org.clojure/clojure "1.3.0"]])
+  :license {:name "Eclipse Public License"
+            :url "http://www.eclipse.org/legal/epl-v10.html"}
+  :dependencies [[org.clojure/clojure "1.4.0"]]
+  :profiles {:dev 
+             {:dependencies [[org.clojure/clojure "1.5.1"]]}})

src/loom/alg.clj

@@ -3,7 +3,8 @@ Graph, Digraph, or WeightedGraph protocols (as appropriate per algorithm)
 can use these functions."
       :author "Justin Kramer"}
   loom.alg
-  (:require [loom.alg-generic :as gen])
+  (:require [loom.alg-generic :as gen]
+            [loom.flow :as flow])
   (:use [loom.graph
          :only [add-edges nodes edges neighbors weight incoming degree
                 in-degree weighted? directed? graph transpose]
@@ -158,11 +159,90 @@ can use these functions."
   [g start end]
   (first (dijkstra-path-dist g start end)))
 
+(defn- can-relax-edge?
+  "Test for whether we can improve the shortest path to v found so far
+   by going through u."
+  [[u v :as edge] weight costs]
+  (let [vd (get costs v)
+        ud (get costs u)
+        sum (+ ud weight)]
+    (if (> vd sum)
+      true
+      false)))
+
+(defn- relax-edge
+  "If there's a shorter path from s to v via u,
+    update our map of estimated path costs and
+   map of paths from source to vertex v"
+  [[u v :as edge] weight [costs paths :as estimates]]
+  (let [ud (get costs u)
+        sum (+ ud weight)]
+    (if (can-relax-edge? edge weight costs)
+      [(assoc costs v sum) (assoc paths v u)]
+      estimates
+      )))
+
+(defn- relax-edges
+  "Performs edge relaxation on all edges in weighted directed graph"
+  [g start estimates]
+  (->> (edges g)
+       (reduce (fn [estimates [u v :as edge]]
+                 (relax-edge edge (wt g u v) estimates))
+               estimates)))
+
+(defn- init-estimates
+  "Initializes path cost estimates and paths from source to all vertices,
+   for Bellman-Ford algorithm"
+  [graph start]
+  (let [nodes (disj (nodes graph) start)
+        path-costs {start 0}
+        paths {start nil}
+        infinities (repeat Double/POSITIVE_INFINITY)
+        nils (repeat nil)
+        init-costs (interleave nodes infinities)
+        init-paths (interleave nodes nils)]
+    [(apply assoc path-costs init-costs)
+     (apply assoc paths init-paths)]))
+
 
 ;;;
 ;;; Graph algorithms
 ;;;
 
+(defn bellman-ford
+  "Given a weighted, directed graph G = (V, E) with source start,
+   the Bellman-Ford algorithm produces map of single source shortest paths and their costs
+   if no negative-weight cycle that is reachable from the source exits,
+   and false otherwise, indicating that no solution exists."
+  [g start]
+  (let [initial-estimates (init-estimates g start)
+        ;relax-edges is calculated for all edges V-1 times
+        [costs paths] (reduce (fn [estimates _]
+                                (relax-edges g start estimates))
+                              initial-estimates
+                              (-> g nodes count dec range))
+        edges (edges g)]
+      (if (some
+            (fn [[u v :as edge]]
+              (can-relax-edge? edge (wt g u v) costs))
+            edges)
+        false
+        [costs 
+         (->> (keys paths)
+              ;remove vertices that are unreachable from source
+              (remove #(= Double/POSITIVE_INFINITY (get costs %))) 
+              (reduce
+                (fn [final-paths v]
+                  (assoc final-paths v
+                         ; follows the parent pointers
+                         ; to construct path from source to node v
+                         (loop [node v
+                                path ()]
+                           (if node
+                             (recur (get paths node) (cons node path))
+                             path))))
+                {}))])))
+
 (defn dag?
   "Return true if g is a directed acyclic graph"
   [g]
@@ -324,4 +404,23 @@ can use these functions."
      [#{} #{}]
      coloring)))
 
-;; TODO: MST, coloring, matching, etc etc
+(defn max-flow
+  "Returns [flow-map flow-value], where flow-map is a weighted adjacency map
+   representing the maximum flow.  The argument should be a weighted digraph,
+   where the edge weights are flow capacities.  Source and sink are the vertices
+   representing the flow source and sink vertices.  Optionally, pass in
+     :method :algorithm to use.  Currently, the only option is :edmonds-karp ."
+  [g source sink & {:keys [method] :or {method :edmonds-karp}}]
+  (let [method-set #{:edmonds-karp}
+        n (nb g), i (incoming g), c (wt g), s source, t sink
+        [flow-map flow-value] (case method
+                                    :edmonds-karp (flow/edmonds-karp n i c s t)
+                                    (throw
+                                     (java.lang.RuntimeException.
+                                      (str "Method not found.  Choose from: "
+                                                 method-set))))]
+    [flow-map flow-value]))
+
+
+
+;; TODO: MST, coloring, matching, etc etc

src/loom/flow.clj

@@ -0,0 +1,116 @@
+(ns ^{:doc "Algorithms for solving network flow"
+      :author "Robert Lachlan"}
+  loom.flow
+  (:require [loom.alg-generic :as gen :only [bf-path]]))
+
+
+(defn residual-capacity [capacity flow v1 v2]
+  "Computes the residual capacity between nodes v1 and v2.  Capacity is a function
+   that takes two nodes, and returns the capacity on the edge between them, if
+   any.  Flow is the adjacency map which represents the current flow in the
+   network."
+  (+
+   (or (get-in flow [v2 v1]) 0)
+   (- (or (capacity v1 v2) 0)
+      (or (get-in flow [v1 v2]) 0))))
+
+(defn flow-balance [flow]
+  "Given a flow, returns a map of {node (sum(in weight) - sum(out weight))}"
+  (loop [out {}, in {}, adj-list (seq flow)]
+    (if-let [[node neighbours] (first adj-list)]
+      (recur (assoc out node (- (reduce + (vals neighbours))))
+             (merge-with + in neighbours)
+             (next adj-list))
+      (merge-with + out in))))
+
+(defn satisfies-mass-balance? [flow source sink]
+  "Given a flow, verifies whether at each node the sum of in edge weights is equal
+   to the sum of out edge weights, except at the source and sink. The source
+   should have positive net outflow, the sink negative, and together they
+   should balance."
+  (let [balance (flow-balance flow)]
+    (and (<= (or (get balance source) 0) 0)
+         (zero? (+ (or (get balance source) 0)
+                   (or (get balance sink) 0)))
+         (every? zero? (vals (dissoc balance source sink))))))
+
+(defn satisfies-capacity-constraints? [flow capacity]
+  "Given a flow map, and a capacity function, verifies that the flow on each edge
+   is <= capacity of that edge."
+  (empty?
+   (remove (fn [[node flow-to-neighbors]]
+             (every?
+              (fn [[neighbor flow-value]] (<= flow-value (capacity node neighbor)))
+              (seq flow-to-neighbors)))
+           (seq flow))))
+
+(defn is-admissible-flow? [flow capacity source sink]
+  "Verifies that a flow satisfies capacity and mass balance constraints.  Does
+   verify that a flow is maximum."
+  (and (satisfies-mass-balance? flow source sink)
+       (satisfies-capacity-constraints? flow capacity)))
+
+(defn min-weight-along-path [path weight-fn]
+  "Given a path, represented by a sequence of nodes, and weight-function, computes
+   the minimum of the edge weights along the path.  If an edge on the path is
+   missing, returns 0."
+  (reduce min (map #(or (apply weight-fn %) 0)  (partition 2 1 path))))
+
+(defn bf-find-augmenting-path [neighbors incoming capacity flow s t]
+  "Finds a shortest path in the flow network along which there remains residual
+   capacity.  Neighbours is a function which, given a vertex, returns the
+   vertices connected by outgoing edges.  Incoming, similarly is a function to get
+   vertices connected by incoming edges.  Capacity is a function which takes two
+   vertices and returns the capacity between them.  Flow is an adjacency map which
+   contains the current value of network flow.  s is the source node, t the sink."
+  (gen/bf-path
+   (fn [vertex] (distinct (filter #(> (residual-capacity capacity flow vertex %) 0)
+                                  (concat (neighbors vertex) (incoming vertex)))))
+   s t))
+
+(defn augment-along-path [flow capacity path increase]
+  "Given a flow represented as an adjacency map, returns an updated flow.
+   Capacity is a function of two vertices, path is a sequence of nodes, and
+   increase is the amount by which the flow should be augmented on this path.  An
+   exception is thrown if the augmentation is impossible given capacity
+   constraints."
+  (let [vn0 (first path)
+        vn1 (second path)
+        forward-flow (or (get-in flow [vn0 vn1]) 0)
+        forward-capacity (- (or (capacity vn0 vn1) 0) forward-flow)
+        reverse-flow (or (get-in flow [vn1 vn0]) 0)
+        forward-increase (min forward-capacity increase)
+        pushback (- increase forward-increase)
+        flow_1 (if (pos? forward-increase)
+                 (assoc-in flow [vn0 vn1] (+ forward-flow forward-increase)) 
+                 flow)
+        flow_2 (if (pos? pushback)
+                 (assoc-in flow_1 [vn1 vn0] (- reverse-flow pushback))
+                 flow_1)]
+    (cond (> pushback reverse-flow) (throw (java.lang.RuntimeException.
+                                            (str "Path augmentation failure: " vn0 " " vn1)))
+          (> (count path) 2) (recur flow_2 capacity (next path) increase)
+          :else flow_2)))
+
+(defn edmonds-karp
+  "Computes the maximum flow on a network, using the edmonds-karp algorithm.
+   Neighbors is a function that returns the outgoing neighbor vertices of a
+   vertex.  Incoming is a function that returns the incoming neighbor vertices
+   for a vertex.  Capacity is a function of two vertices that returns the
+   capacity on the edge between them.  Source and sink are the unique vertices
+   which supply and consume flow respectively.
+
+   Returns a vector [flow value], where flow is an adjacency map that represents
+   flows between vertices, and value is the quantity of flow passing from source
+   to sink."
+  ([neighbors incoming capacity source sink]
+     (edmonds-karp neighbors incoming capacity source sink {}))
+  ([neighbors incoming capacity source sink flow]
+     (if-let [path (bf-find-augmenting-path neighbors
+                                            incoming capacity flow source sink)]
+       (recur neighbors incoming capacity source sink
+              (augment-along-path flow capacity path
+                                  (min-weight-along-path
+                                   path (partial residual-capacity capacity flow))))
+       (let [value (reduce + (vals (get flow source)))]
+       [flow value]))))

src/loom/graph.clj

@@ -7,8 +7,6 @@ on adjacency lists."
   loom.graph
   (:use [loom.alg-generic :only [bf-traverse]]))
 
-;(set! *warn-on-reflection* true)
-
 ;;;
 ;;; Protocols
 ;;;
@@ -34,6 +32,8 @@ on adjacency lists."
 (defprotocol WeightedGraph
   (weight [g] [g n1 n2] "Return weight of edge [n1 n2] or (partial weight g)"))
 
+(defprotocol VertexLabelled)
+
 ;; Variadic wrappers
 
 (defn add-nodes
@@ -338,6 +338,9 @@ on adjacency lists."
 (defrecord WeightedFlyGraph [fnodes fedges fneighbors fweight start])
 (defrecord WeightedFlyDigraph [fnodes fedges fneighbors fincoming fweight start])
 
+;; Deprecate the flygraphs?  Instead provide interfaces on algorithms to 
+;; run the algorithm on 
+
 (extend FlyGraph
   Graph default-flygraph-graph-impl)
 

test/loom/test/alg.clj

@@ -67,6 +67,25 @@
             :g [:f]
             :h [:g :d]}))
 
+;; Weighted directed graph with a negative-weight cycle
+;; which is reachable from sources :a, :b, :d, and :e.
+;; http://www.seas.gwu.edu/~simhaweb/alg/lectures/module9/module9.html
+(def g11
+  (weighted-digraph [:a :b 3]
+                    [:b :c 4]
+                    [:b :d 5]
+                    [:d :e 2]
+                    [:e :b -8]))
+
+;; Weighted directed graph with a non-negative-weight cycle,
+;; similar to g11, but with the edge [:e :b] reweighed.
+(def g12
+  (weighted-digraph [:a :b 3]
+                    [:b :c 4]
+                    [:b :d 5]
+                    [:d :e 2]
+                    [:e :b -7]))
+
 (deftest depth-first-test
   (are [expected got] (= expected got)
     #{1 2 3 5 6 7} (set (pre-traverse g7))
@@ -139,6 +158,63 @@
     ;; TODO: the rest
     ))
 
+(deftest bellman-ford-test
+  (are [expected graph start]
+       (= expected (bellman-ford graph start))
+
+       false g11 :a
+       false g11 :b
+       [{:e Double/POSITIVE_INFINITY,
+         :d Double/POSITIVE_INFINITY,
+         :b Double/POSITIVE_INFINITY,
+         :a Double/POSITIVE_INFINITY,
+         :c 0}{:c [:c]}] g11 :c 
+       false g11 :d
+       false g11 :e
+       [{:e 10,
+         :d 8,
+         :b 3,
+         :c 7,
+         :a 0}
+        {:a [:a],
+         :c [:a :b :c],
+         :b [:a :b],
+         :d [:a :b :d],
+         :e [:a :b :d :e]}] g12 :a
+       [{:e 7,
+         :d 5,
+         :c 4,
+         :a Double/POSITIVE_INFINITY,
+         :b 0}
+        {:b [:b],
+         :c [:b :c],
+         :d [:b :d],
+         :e [:b :d :e]}] g12 :b
+       [{:e Double/POSITIVE_INFINITY,
+         :d Double/POSITIVE_INFINITY,
+         :b Double/POSITIVE_INFINITY,
+         :a Double/POSITIVE_INFINITY,
+         :c 0}
+        {:c [:c]}] g12 :c 
+       [{:e 2,
+         :b -5,
+         :c -1,
+         :a Double/POSITIVE_INFINITY,
+         :d 0}
+        {:d [:d],
+         :c [:d :e :b :c],
+         :b [:d :e :b],
+         :e [:d :e]}] g12 :d
+       [{:d -2,
+         :b -7,
+         :c -3,
+         :a Double/POSITIVE_INFINITY,
+         :e 0}
+        {:e [:e],
+         :c [:e :b :c],
+         :b [:e :b],
+         :d [:e :b :d]}] g12 :e))
+
 (deftest bipartite-test
   (are [expected got] (= expected got)
     {0 1, 1 0, 5 0, 2 1, 3 1, 4 0} (bipartite-color g6)