Dijkstra's algorithm

This article assumes familiarity with Dijkstra's shortest path algorithm. For a refresher, see [1]. The code assumes open Core is in effect and is online here.

The first part of the program organizes our thoughts about what we are setting out to compute. The signature summarizes the notion (for our purposes) of a graph definition in modular form. A module implementing this signature defines a type vertex_t for vertices, a type t for graphs and type extern_t : a representation of a t for interaction between an implemening module and its "outside world".

module type Graph_sig = sig
  type vertex_t [@@deriving sexp]
  type t [@@deriving sexp]
  type extern_t

  type load_error = [ `Duplicate_vertex of vertex_t ] [@@deriving sexp]
  exception Load_error of load_error [@@deriving sexp]

  val of_adjacency : extern_t -> [ `Ok of t | `Load_error of load_error ]
  val to_adjacency : t -> extern_t

  module Dijkstra : sig
    type state

    type error = [
      | `Relax of vertex_t
    ] [@@deriving sexp]
    exception Error of error [@@deriving sexp]

    val dijkstra : vertex_t -> t -> [ `Ok of state | `Error of error ]
    val d : state -> (vertex_t * float) list
    val shortest_paths : state -> (vertex_t * vertex_t list) list
  end

end

A realization of Graph_sig provides "conversion" functions of_adjacency/to_adjacency between the types extern_t and t and nests a module Dijkstra. The signature of the sub-module Dijkstra requires concrete modules provide a type state and an implementation of Dijkstra's algorithm in terms of the function signature val dijkstra : vertex_t -> t -> [ `Ok of state | `Error of error ].

For reusability, the strategy for implementing graphs will be generic programming via functors over modules implementing s vertex type.

An implementation of the module type GRAPH defines a module type VERT which is required to provide a comparable type t. It further defines a module type S that is exactly module type Graph_sig above. Lastly, modules of type GRAPH provide a functor Make that maps any module of type VERT to new module of type S fixing extern_t to an adjacency list representation in terms of the native OCaml type 'a list and float to represent weights on edges.

module type GRAPH = sig
  module type VERT = sig
    type t[@@deriving sexp]
    include Comparable.S with type t := t
  end

  module type S = sig
    include Graph_sig
  end

  module Make : functor (V : VERT) ->
    S with type vertex_t = V.t
       and type extern_t = (V.t * (V.t * float) list) list
end

The two module types Graph_sig and GRAPH together provide the specification for the program. module Graph in the next section implements this specification.

Implementation of module Graph is in outline this.

module Graph : GRAPH = struct
  module type VERT = sig
    type t[@@deriving sexp]
    include Comparable.S with type t := t
  end

  module type S = sig
    include Graph_sig
  end

  module Make : functor (V : VERT) ->
    S with type vertex_t = V.t
       and type extern_t = (V.t * (V.t * float) list) list
    =

    functor (V : VERT) -> struct
       ...
    end
end

As per the requirements of GRAPH the module types VERT and S are provided as is the functor Make. It is the code that is ellided by the ... above in the definition of Make that is now the focus.

Modules produced by applications of Make satisfy S. This requires suitable definitions of types vertext_t, t and extern_t. The modules Map and Set are available due to modules of type VERT being comparable in their type t.

      module Map = V.Map
      module Set = V.Set

      type vertex_t = V.t [@@deriving sexp]
      type t = (vertex_t * float) list Map.t [@@deriving sexp]
      type extern_t = (vertex_t * (vertex_t * float) list) list
      type load_error = [ `Duplicate_vertex of vertex_t ] [@@deriving sexp]
      exception Load_error of load_error [@@deriving sexp]

While the external representation extern_t of graphs is chosen to be an adjacency list representation in terms of association lists, the internal representation t is a vertex map of adjacency lists providing logarithmic loookup complexity. The conversion functions between the two representations "come for free" via module Map.

      let to_adjacency g = Map.to_alist g

      let of_adjacency_exn l =  match Map.of_alist l with
        | `Ok t -> t
        | `Duplicate_key c -> raise (Load_error (`Duplicate_vertex c))

      let of_adjacency l =
        try
          `Ok (of_adjacency_exn l)
        with
        | Load_error err -> `Load_error err

At this point the "scaffolding" for Dijkstra's algorithm, that part of GRAPH dealing with the representation of graphs is implemented.

The interpretation of Dijkstra's algorithm we adopt is functional : the idea is we loop over vertices relaxing their edges until all shortest paths are known. What we know on any recursive iteration of the loop is a current "state" (of the computation) and each iteration produces a new state. This next definition is the formal definition of type state.

      module Dijkstra = struct

        type state = {
          src    :                  vertex_t
        ; g      :                         t
        ; d      :               float Map.t
        ; pred   :            vertex_t Map.t
        ; s      :                     Set.t
        ; v_s    : (vertex_t * float) Heap.t
        }

The fields of this record are:

• src : vertex_t, the source vertex;
• g : t, G the graph;
• d : float Map.t, d the shortest path weight estimates;
• pre : vertex_t Map.t, π the predecessor relation;
• s : Set.t, the set S of nodes for which the lower bound shortest path weight is known;
• v_s : (vertex_t * float) Heap.t, V - {S}, , the set of nodes of g for which the lower bound of the shortest path weight is not yet known ordered on their estimates.

Function invocation init src g compuates an initial state for the graph g containing the source node src. In the initial state, d is everywhere ∞ except for src which is 0. Set S (i.e. s) and the predecessor relation π (i.e. pred) are empty and the set V - {S} (i.e. v_s) contains all nodes.

        let init src g =
          let init x = match V.equal src x with
            | true -> 0.0 | false -> Float.infinity in
          let d = List.fold (Map.keys g) ~init:Map.empty
              ~f:(fun acc x -> Map.set acc ~key:x ~data:(init x)) in
          {
            src
          ; g
          ; s = Set.empty
          ; d
          ; pred = Map.empty
          ; v_s = Heap.of_list (Map.to_alist d)
                ~cmp:(fun (_, e1) (_, e2) -> Float.compare e1 e2)
          }

Relaxing an edge (u, v) with weight w (u, v) tests whether the shortest path to v so far can be improved by going through u and if so, updating d (v) and π (v) accordingly.

        type error = [
          | `Relax of vertex_t
        ] [@@deriving sexp]
        exception Error of error [@@deriving sexp]

        let relax state (u, v, w) =
          let {d; pred; v_s; _} = state in
          let dv = match Map.find d v with
            | Some dv -> dv
            | None -> raise (Error (`Relax v)) in
          let du = match Map.find d u with
            | Some du -> du
            | None -> raise (Error (`Relax u)) in
          if dv > du +. w then
            let dv = du +. w in
            (match Heap.find_elt v_s ~f:(fun (n, _) -> V.equal n v) with
            | Some tok -> ignore (Heap.update v_s tok (v, dv))
            | None -> raise (Error (`Relax v))
            );
            { state with
              d = Map.change d v
                  ~f:(function
                      | Some _ -> Some dv
                      | None -> raise (Error (`Relax v))
                    )
            ; pred = Map.set (Map.remove pred v) ~key:v ~data:u
            }
          else state

Here, relaxation can result in a linear heap update operation. A better implementation might seek to avoid that.

One iteration of the body of the loop of Dijkstra's algorithm consists of the node in V - {S} with the least shortest path weight estimate being moved to S and its edges relaxed.

        let dijkstra_exn src g =
          let rec loop ({s; v_s; _} as state) =
            match Heap.is_empty v_s with
            | true -> state
            | false ->
              let u = fst (Heap.pop_exn v_s) in
              loop (
                List.fold (Map.find_exn g u)
                  ~init:{ state with s = Set.add s u }
                  ~f:(fun state (v, w) -> relax state (u, v, w))
              )
          in loop (init src g)

        let dijkstra src g =
          try
            `Ok (dijkstra_exn src g)
          with
          | Error err -> `Error err

The shortest path estimates contained by a value of state is given by the projection d.

        let d state = Map.to_alist (state.d)

The shortest paths themselves are easily computed as,

   let path state n =
          let rec loop acc x =
            (match V.equal x state.src with
            | true -> x :: acc
            | false -> loop (x :: acc) (Map.find_exn state.pred x)
            ) in
          loop [] n

        let shortest_paths state =
          List.map (Map.keys state.g) ~f:(fun n -> (n, path state n))
      end
    end

which completes the implementation of Make.

The following program produces a concrete instance of the shortest path problem (with some evaluation output from the top-level).

module G : Graph.S with
  type vertex_t = char and type extern_t = (char * (char * float) list) list
  =
  Graph.Make (Char)

let g : G.t =
  match G.of_adjacency
          [ 's', ['u',  3.0; 'x', 5.0]
          ; 'u', ['x',  2.0; 'v', 6.0]
          ; 'x', ['v',  4.0; 'y', 6.0; 'u', 1.0]
          ; 'v', ['y',  2.0]
          ; 'y', ['v',  7.0]
          ]
  with
  | `Ok g -> g
  | `Load_error e -> failwiths "Graph load error : %s" e G.sexp_of_load_error
;;
let s = match (G.Dijkstra.dijkstra 's' g) with
  | `Ok s -> s
  | `Error e -> failwiths "Error : %s" e G.Dijkstra.sexp_of_error

;; G.Dijkstra.d s
- : (char * float) list =
[('s', 0.); ('u', 3.); ('v', 9.); ('x', 5.); ('y', 11.)]

;; G.Dijkstra.shortest_paths s
- : (char * char list) list =
[('s', ['s']); ('u', ['s'; 'u']); ('v', ['s'; 'u'; 'v']); ('x', ['s'; 'x']);
 ('y', ['s'; 'x'; 'y'])]

---

References:
[1] "Introduction to Algorithms" Section 24.3:Dijkstra's algorithm -- Cormen et. al. (Second ed.) 2001.