3 people like it.

Max Flow - push-relabel algorithm

Push-relabel algorithm for Max Flow problems in flow graphs. (Does not check if graph has improper back edges)

  1: 
  2: 
  3: 
  4: 
  5: 
  6: 
  7: 
  8: 
  9: 
 10: 
 11: 
 12: 
 13: 
 14: 
 15: 
 16: 
 17: 
 18: 
 19: 
 20: 
 21: 
 22: 
 23: 
 24: 
 25: 
 26: 
 27: 
 28: 
 29: 
 30: 
 31: 
 32: 
 33: 
 34: 
 35: 
 36: 
 37: 
 38: 
 39: 
 40: 
 41: 
 42: 
 43: 
 44: 
 45: 
 46: 
 47: 
 48: 
 49: 
 50: 
 51: 
 52: 
 53: 
 54: 
 55: 
 56: 
 57: 
 58: 
 59: 
 60: 
 61: 
 62: 
 63: 
 64: 
 65: 
 66: 
 67: 
 68: 
 69: 
 70: 
 71: 
 72: 
 73: 
 74: 
 75: 
 76: 
 77: 
 78: 
 79: 
 80: 
 81: 
 82: 
 83: 
 84: 
 85: 
 86: 
 87: 
 88: 
 89: 
 90: 
 91: 
 92: 
 93: 
 94: 
 95: 
 96: 
 97: 
 98: 
 99: 
100: 
101: 
102: 
103: 
104: 
105: 
106: 
107: 
108: 
109: 
110: 
111: 
112: 
113: 
114: 
115: 
116: 
117: 
118: 
119: 
120: 
121: 
122: 
123: 
124: 
125: 
126: 
127: 
128: 
129: 
130: 
131: 
132: 
133: 
134: 
135: 
136: 
137: 
138: 
139: 
140: 
141: 
142: 
143: 
144: 
145: 
146: 
147: 
148: 
149: 
150: 
151: 
152: 
153: 
154: 
155: 
156: 
157: 
158: 
159: 
160: 
161: 
162: 
(*
Max Flow - push-relabel algorithm
*)
type VType = S | T | V                                             // vertex type S=source; T=sink; V=internal vertex
type V = {Id:string; T : VType}                                    // vertex in input graph
type E = {Target:string; C:float}                                  // edge in input graph C = capacity
type GInp = (V * E list) list                                      // input graph - adjacency list
type Ef = {Src:string; Tgt:string; C:float; Fwd:float; Bck:float}  // flow in residual graph (forward and backward flows are together)
type Gf = Map<string*string,Ef>                                    // residual graph
type Vst = {H:int; E:float}                                        // node Height and overflow (E)
type HF = Map<string,Vst>                                          // height and overflow state kept by node id
type PQ = Set<int*string>                                          // use Set as poor man's priority queue (-height * node id)
type G = Map<string, (V * E list * string list)>                   // cached internal graph structure with back node pointers (static over the alg)
type S = G * (Gf * HF * PQ)                                        // captures current state of the process
type Dir = Fwd | Bck                                               // direction of flow forward (src --> tgt) or backward

let order (a:string,b:string) = if a < b then (a,b) else (b,a)     // res graph key is lexical order of src tgt

//add a flow to flow graph
let setFlow (ef:Ef) (gf:Gf) =
    gf |> Map.add (order(ef.Src,ef.Tgt)) ef
    
let getFlow id tgt (gf:Gf) = gf |> Map.tryFind (order(id,tgt))
    
//create initial state from input graph 
let initialize (g:GInp) : S =
    let incoming = 
        (Map.empty,g) 
        ||> List.fold (fun acc (v,es) -> 
            (acc,es) 
            ||> List.fold (fun acc e -> 
                let ls =
                    acc 
                    |> Map.tryFind e.Target
                    |> Option.defaultValue []
                acc |> Map.add e.Target (v.Id::ls)))
    let g = (Map.empty,g) ||> List.fold (fun acc (v,es) -> acc |> Map.add v.Id (v,es, incoming |> Map.tryFind v.Id |> Option.defaultValue []))
    let vertices = Map.count g
    let st =
        ((Map.empty,Map.empty,Set.empty),g)
        ||> Map.fold (fun (gf,hf,pq) _ (n,es,_) ->
            match n.T with            
            | V | T ->
                let gf = 
                    (gf,es) 
                    ||> List.fold (fun gf e -> gf |> setFlow {Src=n.Id; Tgt=e.Target; C=e.C; Fwd=0.0; Bck=0.0})
                let hf = hf |> Map.add n.Id {H=0; E=0.0}
                gf,hf,pq
            | S -> 
                let gf =
                    (gf,es) 
                    ||> List.fold (fun gf e -> gf |> setFlow {Src=n.Id; Tgt=e.Target; C=e.C; Fwd=e.C; Bck=0.0})
                let hf = hf |> Map.add n.Id {H=vertices; E=0.0} 
                let hf = (hf,es) ||> List.fold (fun hf e -> hf |> Map.add e.Target {H=0; E=e.C})
                let pq = (pq,es) ||> List.fold (fun pq e -> pq |> Set.add (-0,e.Target))
                gf,hf,pq)
    g,st
    
//determine if flow - in the given direction - has spare capacity
let hasCap = function  Fwd,f -> f.Fwd < f.C | Bck,f -> f.Bck < f.Fwd

//increase height of node to the given height
let increaseHeight id h2 (s:S) =
    let g,(gf,hf,pq) = s
    let (v,_,_) = g.[id]
    match v.T with S  | T -> failwith "cannot change height of source and sink" | _ -> ()
    let vHf = hf.[id]
    let pq = pq |> Set.remove (-vHf.H,id)             //update height in priority q
    let pq = pq |> Set.add (-h2,id)
    let hf = hf |> Map.add id ({hf.[id] with H=h2})
    g,(gf,hf,pq)

//finds unsaturated edges out from the given node (includes forward and backward edges)
let unsatEs id (s:S) = 
    let g,(gf,_,_) = s
    let (v,es,incming) = g.[id]         
    let unsat_Fwd = 
        es 
        |> List.choose (fun e -> 
            getFlow id e.Target gf
            |> Option.map (fun f -> Fwd,f)
            |> Option.filter hasCap)
    let unsat_Bck =
        incming 
            |> List.choose (fun inc -> 
                getFlow id inc gf
                |> Option.map (fun f -> Bck,f)
                |> Option.filter hasCap)
    unsat_Fwd @ unsat_Bck

//push deltaFlow from src to tgt (does not check constraints)
let pushFlow src tgt deltaFlow (s:S) =
    let g,(gf,hf,pq) = s
    let (vTgt,_,_) = g.[tgt]
    let srcHf = hf.[src]
    let tgtHf = hf.[tgt]
    let f = getFlow src tgt gf |> Option.get 
    let f =
        if f.Src = src then 
            {f with Fwd = f.Fwd + deltaFlow}
        else
            {f with Bck = f.Bck + deltaFlow}
    let gf = setFlow f gf
    let srcHf = {srcHf with E = srcHf.E - deltaFlow}
    let tgtHf = {tgtHf with E = tgtHf.E + deltaFlow}
    let hf = hf |> Map.add src srcHf |> Map.add tgt tgtHf
    let pq = if srcHf.E <= 0.0 then pq |> Set.remove (-srcHf.H,src) else pq
    let pq = match vTgt.T with V when tgtHf.E > 0.0 -> pq |> Set.add (-tgtHf.H,tgt) | _ -> pq
    g,(gf,hf,pq)

//push of push-relable algorithim
let push id (s:S) =
    let g,(gf,hf,pq) = s
    let srcHf = hf.[id]
    let e_unsat = unsatEs id s |> List.tryFind (function Fwd,f -> srcHf.H > hf.[f.Tgt].H | Bck,f -> srcHf.H > hf.[f.Src].H)
    e_unsat
    |> Option.map (fun (dir,f) -> 
        let src,tgt,df =
            match dir with
            | Fwd -> f.Src, f.Tgt, min srcHf.E (f.C - f.Fwd)
            | Bck -> f.Tgt, f.Src, min srcHf.E (f.Fwd - f.Bck)
        pushFlow src tgt df s)
                
//relabel (increase height) of the given node
let relabel id (s:S) =
    let g,(gf,hf,pq) = s
    let srcHf = hf.[id]
    let v,es,incmng = g.[id]
    let unsates = unsatEs id s
    let minHNbr = unsates |> List.map (function Fwd,f -> hf.[f.Tgt].H | Bck,f->hf.[f.Src].H) |> List.min
    increaseHeight id (minHNbr+1) s
        
//push or relabel active node with the highest height
let pushRelabel (s:S) =
    let g,(gf,hf,pq) = s
    let _,nActive = Set.minElement pq //highest 
    match push nActive s with
    | Some s' -> s'
    | None    -> relabel nActive s

//push-relabel main loop
let maxFlow (s:S) =
    let rec loop ((_,(_,_,pq)) as s) =
        if Set.isEmpty pq then
            s
        else
            let s' = pushRelabel s
            loop s'
    loop s

let printFlow (s:S) =
    let g,(gf,_,_) = s
    let (v,es,_) = g |> Map.toSeq |> Seq.map snd |> Seq.find (fun (v,es,_) -> v.T = S) //source
    let maxFlow = es |> List.choose (fun e -> getFlow v.Id e.Target gf) |> List.sumBy (fun f->f.Fwd-f.Bck)
    printfn "** Maxflow: %0.02f **" maxFlow
    gf 
    |> Map.toSeq
    |> Seq.map snd
    |> Seq.filter (fun f -> f.Fwd - f.Bck > 0.0)
    |> Seq.iter (fun f -> 
        printfn "%s --> %s : %0.02f" f.Src f.Tgt (f.Fwd - f.Bck)
    )

Test case *

 1: 
 2: 
 3: 
 4: 
 5: 
 6: 
 7: 
 8: 
 9: 
10: 
11: 
12: 
13: 
14: 
15: 
16: 
17: 
18: 
19: 
20: 
21: 
22: 
23: 
24: 
25: 
26: 
let ginp1 =
    [
        {Id="Vancouver"; T=S}, [{Target="Edmonton"; C=16.0}; {Target="Calgary"; C=13.0}]
        {Id="Edmonton"; T=V}, [{Target="Saskatoon"; C=12.}]
        {Id="Calgary"; T=V}, [{Target="Edmonton";C=4.0}; {Target="Regina"; C=14.0}]
        {Id="Regina"; T=V}, [{Target="Saskatoon"; C=7.0}; {Target="Winnipeg"; C=4.0;}]
        {Id="Saskatoon";T=V}, [{Target="Winnipeg"; C=20.}; {Target="Calgary";C=9.0}]
        {Id="Winnipeg"; T=T}, []
    ]

let s1 = initialize ginp1

let s1g = maxFlow s1

printFlow s1g

(*
** Maxflow: 23.00 **
Calgary --> Regina : 11.00
Vancouver --> Calgary : 11.00
Edmonton --> Saskatoon : 12.00
Vancouver --> Edmonton : 12.00
Regina --> Saskatoon : 7.00
Regina --> Winnipeg : 4.00
Saskatoon --> Winnipeg : 19.00
*)
union case VType.S: VType
union case VType.T: VType
union case VType.V: VType
type V =
  { Id: string
    T: VType }
V.Id: string
Multiple items
val string : value:'T -> string

--------------------
type string = System.String
V.T: VType
type VType =
  | S
  | T
  | V
type E =
  { Target: string
    C: float }
E.Target: string
E.C: float
Multiple items
val float : value:'T -> float (requires member op_Explicit)

--------------------
type float = System.Double

--------------------
type float<'Measure> = float
type GInp = (V * E list) list
type 'T list = List<'T>
type Ef =
  { Src: string
    Tgt: string
    C: float
    Fwd: float
    Bck: float }
Ef.Src: string
Ef.Tgt: string
Ef.C: float
Ef.Fwd: float
Ef.Bck: float
type Gf = Map<(string * string),Ef>
Multiple items
module Map

from Microsoft.FSharp.Collections

--------------------
type Map<'Key,'Value (requires comparison)> =
  interface IReadOnlyDictionary<'Key,'Value>
  interface IReadOnlyCollection<KeyValuePair<'Key,'Value>>
  interface IEnumerable
  interface IComparable
  interface IEnumerable<KeyValuePair<'Key,'Value>>
  interface ICollection<KeyValuePair<'Key,'Value>>
  interface IDictionary<'Key,'Value>
  new : elements:seq<'Key * 'Value> -> Map<'Key,'Value>
  member Add : key:'Key * value:'Value -> Map<'Key,'Value>
  member ContainsKey : key:'Key -> bool
  ...

--------------------
new : elements:seq<'Key * 'Value> -> Map<'Key,'Value>
type Vst =
  { H: int
    E: float }
Vst.H: int
Multiple items
val int : value:'T -> int (requires member op_Explicit)

--------------------
type int = int32

--------------------
type int<'Measure> = int
Vst.E: float
type HF = Map<string,Vst>
type PQ = Set<int * string>
Multiple items
module Set

from Microsoft.FSharp.Collections

--------------------
type Set<'T (requires comparison)> =
  interface IReadOnlyCollection<'T>
  interface IComparable
  interface IEnumerable
  interface IEnumerable<'T>
  interface ICollection<'T>
  new : elements:seq<'T> -> Set<'T>
  member Add : value:'T -> Set<'T>
  member Contains : value:'T -> bool
  override Equals : obj -> bool
  member IsProperSubsetOf : otherSet:Set<'T> -> bool
  ...

--------------------
new : elements:seq<'T> -> Set<'T>
type G = Map<string,(V * E list * string list)>
type S = G * (Gf * HF * PQ)
type Dir =
  | Fwd
  | Bck
union case Dir.Fwd: Dir
union case Dir.Bck: Dir
val order : a:string * b:string -> string * string
val a : string
val b : string
val setFlow : ef:Ef -> gf:Gf -> Map<(string * string),Ef>
val ef : Ef
val gf : Gf
val add : key:'Key -> value:'T -> table:Map<'Key,'T> -> Map<'Key,'T> (requires comparison)
val getFlow : id:string -> tgt:string -> gf:Gf -> Ef option
val id : string
val tgt : string
val tryFind : key:'Key -> table:Map<'Key,'T> -> 'T option (requires comparison)
val initialize : g:GInp -> S
val g : GInp
val incoming : Map<string,string list>
val empty<'Key,'T (requires comparison)> : Map<'Key,'T> (requires comparison)
Multiple items
module List

from Microsoft.FSharp.Collections

--------------------
type List<'T> =
  | ( [] )
  | ( :: ) of Head: 'T * Tail: 'T list
    interface IReadOnlyList<'T>
    interface IReadOnlyCollection<'T>
    interface IEnumerable
    interface IEnumerable<'T>
    member GetReverseIndex : rank:int * offset:int -> int
    member GetSlice : startIndex:int option * endIndex:int option -> 'T list
    member Head : 'T
    member IsEmpty : bool
    member Item : index:int -> 'T with get
    member Length : int
    ...
val fold : folder:('State -> 'T -> 'State) -> state:'State -> list:'T list -> 'State
val acc : Map<string,string list>
val v : V
val es : E list
val e : E
val ls : string list
module Option

from Microsoft.FSharp.Core
val defaultValue : value:'T -> option:'T option -> 'T
val g : Map<string,(V * E list * string list)>
val acc : Map<string,(V * E list * string list)>
val vertices : int
val count : table:Map<'Key,'T> -> int (requires comparison)
val st : Map<(string * string),Ef> * Map<string,Vst> * Set<int * string>
val empty<'T (requires comparison)> : Set<'T> (requires comparison)
val fold : folder:('State -> 'Key -> 'T -> 'State) -> state:'State -> table:Map<'Key,'T> -> 'State (requires comparison)
val gf : Map<(string * string),Ef>
val hf : Map<string,Vst>
val pq : Set<int * string>
val n : V
val add : value:'T -> set:Set<'T> -> Set<'T> (requires comparison)
val hasCap : Dir * Ef -> bool
val f : Ef
val increaseHeight : id:string -> h2:int -> G * (Gf * HF * PQ) -> G * (Gf * Map<string,Vst> * Set<int * string>)
val h2 : int
val s : S
val g : G
val hf : HF
val pq : PQ
val failwith : message:string -> 'T
val vHf : Vst
val remove : value:'T -> set:Set<'T> -> Set<'T> (requires comparison)
val unsatEs : id:string -> G * (Gf * HF * PQ) -> (Dir * Ef) list
val incming : string list
val unsat_Fwd : (Dir * Ef) list
val choose : chooser:('T -> 'U option) -> list:'T list -> 'U list
val map : mapping:('T -> 'U) -> option:'T option -> 'U option
val filter : predicate:('T -> bool) -> option:'T option -> 'T option
val unsat_Bck : (Dir * Ef) list
val inc : string
val pushFlow : src:string -> tgt:string -> deltaFlow:float -> G * (Gf * HF * PQ) -> G * (Map<(string * string),Ef> * Map<string,Vst> * Set<int * string>)
val src : string
val deltaFlow : float
val vTgt : V
val srcHf : Vst
val tgtHf : Vst
val get : option:'T option -> 'T
val push : id:string -> G * (Gf * HF * PQ) -> (G * (Map<(string * string),Ef> * Map<string,Vst> * Set<int * string>)) option
val e_unsat : (Dir * Ef) option
val tryFind : predicate:('T -> bool) -> list:'T list -> 'T option
val dir : Dir
val df : float
val min : e1:'T -> e2:'T -> 'T (requires comparison)
val relabel : id:string -> G * (Gf * HF * PQ) -> G * (Gf * Map<string,Vst> * Set<int * string>)
val incmng : string list
val unsates : (Dir * Ef) list
val minHNbr : int
val map : mapping:('T -> 'U) -> list:'T list -> 'U list
val min : list:'T list -> 'T (requires comparison)
val pushRelabel : G * (Gf * HF * PQ) -> G * (Map<(string * string),Ef> * Map<string,Vst> * Set<int * string>)
val nActive : string
val minElement : set:Set<'T> -> 'T (requires comparison)
union case Option.Some: Value: 'T -> Option<'T>
val s' : G * (Map<(string * string),Ef> * Map<string,Vst> * Set<int * string>)
union case Option.None: Option<'T>
val maxFlow : G * (Gf * HF * PQ) -> G * (Gf * HF * Set<int * string>)
val loop : (G * (Gf * HF * Set<int * string>) -> G * (Gf * HF * Set<int * string>))
val s : G * (Gf * HF * Set<int * string>)
val isEmpty : set:Set<'T> -> bool (requires comparison)
val printFlow : G * (Gf * HF * PQ) -> unit
val toSeq : table:Map<'Key,'T> -> seq<'Key * 'T> (requires comparison)
module Seq

from Microsoft.FSharp.Collections
val map : mapping:('T -> 'U) -> source:seq<'T> -> seq<'U>
val snd : tuple:('T1 * 'T2) -> 'T2
val find : predicate:('T -> bool) -> source:seq<'T> -> 'T
val maxFlow : float
val sumBy : projection:('T -> 'U) -> list:'T list -> 'U (requires member ( + ) and member get_Zero)
val printfn : format:Printf.TextWriterFormat<'T> -> 'T
val filter : predicate:('T -> bool) -> source:seq<'T> -> seq<'T>
val iter : action:('T -> unit) -> source:seq<'T> -> unit
val ginp1 : (V * E list) list
val s1 : S
val s1g : G * (Gf * HF * Set<int * string>)
Raw view Test code New version

More information

Link:http://fssnip.net/85K
Posted:28 days ago
Author:Faisal Waris
Tags: graphs , graph theory; , max flow