Dijkstra's algorithm: Difference between revisions

Content added Content deleted
Line 758: Line 758:
[https://en.wikipedia.org/w/index.php?title=Dijkstra%27s_algorithm&oldid=1117533081#Using_a_priority_queue a Wikipedia article about Dijkstra's algorithm], although I use a different method to determine whether a queue entry is obsolete.
[https://en.wikipedia.org/w/index.php?title=Dijkstra%27s_algorithm&oldid=1117533081#Using_a_priority_queue a Wikipedia article about Dijkstra's algorithm], although I use a different method to determine whether a queue entry is obsolete.


I prove that the algorithm terminates.
I prove that the algorithm terminates and that the priority queue has at least enough storage. The priority queue is a binary heap implemented as an array.


<syntaxhighlight lang="ats">
<syntaxhighlight lang="ats">
Line 783: Line 783:
macdef infinity = $extval (flt, "INFINITY")
macdef infinity = $extval (flt, "INFINITY")


prfn
(*------------------------------------------------------------------*)
mul_compare_lte
(* A very inefficient implementation of a priority queue, for the sake
{i, j, n : nat | i <= j}
of demonstrating Dijkstra's algorithm. The queue is represented as
an association list. *)
()
:<prf> [i * n <= j * n] void =

mul_gte_gte_gte {j - i, n} ()
typedef pqueue (priority_t : t@ype+,
value_t : t@ype+,
n : int) =
list (@(priority_t, value_t), n)


prfn
prfn
mul_compare_lt
lemma_pqueue_param
{n : int}
{i, j, n : int | 0 <= i; i < j; 1 <= n}
{priority_t, value_t : t@ype}
(pq : pqueue (priority_t, value_t, n))
:<prf> [0 <= n]
void =
lemma_list_param pq

extern fn {priority_t : t@ype}
pqueue$cmp :
(priority_t, priority_t) -<> int

implement pqueue$cmp<double> (x, y) = compare (x, y)

fn
pqueue_make_empty
{priority_t : t@ype}
{value_t : t@ype}
()
()
:<!wrt> pqueue (priority_t, value_t, 0) =
:<prf> [i * n < j * n] void =
mul_compare_lte {i + 1, j, n} ()
NIL

fn {}
pqueue_is_empty
{n : int}
{priority_t : t@ype}
{value_t : t@ype}
(pq : pqueue (priority_t, value_t, n))
:<> bool (n == 0) =
list_is_nil pq

fn {priority_t : t@ype}
{value_t : t@ype}
pqueue_insert
{n : int}
(pq : &pqueue (priority_t, value_t, n)
>> pqueue (priority_t, value_t, n + 1),
priority : priority_t,
value : value_t)
:<!wrt> void =
let
prval () = lemma_pqueue_param pq
in
pq := @(priority, value) :: pq
end

fn {priority_t : t@ype}
{value_t : t@ype}
pqueue_extract_min
{n : pos}
(pq : &pqueue (priority_t, value_t, n)
>> pqueue (priority_t, value_t, n - 1))
:<!wrt> @(priority_t, value_t) =
let
fun
loop {m1, m2 : nat | m1 + m2 + 1 == n}
.<m1>.
(current_min : @(priority_t, value_t),
to_be_scanned : list (@(priority_t, value_t), m1),
already_scanned : list (@(priority_t, value_t), m2))
:<> @(@(priority_t, value_t),
list (@(priority_t, value_t), n - 1)) =
case+ to_be_scanned of
| NIL => @(current_min, already_scanned)
| head :: tail =>
if pqueue$cmp<priority_t> (head.0, current_min.0) < 0 then
loop (head, tail, current_min :: already_scanned)
else
loop (current_min, tail, head :: already_scanned)

val+ head :: tail = pq
val @(retval, pq1) = loop (head, tail, NIL)
in
pq := pq1;
retval
end

(* A little unit testing of the priority queue. *)
var pq = pqueue_make_empty ()
val- true = pqueue_is_empty pq
val () = pqueue_insert (pq, 3.0, 3)
val () = pqueue_insert (pq, 5.0, 5)
val () = pqueue_insert (pq, 1.0, 1)
val () = pqueue_insert (pq, 2.0, 2)
val () = pqueue_insert (pq, 4.0, 4)
val- false = pqueue_is_empty pq
val- @(1.0, 1) = pqueue_extract_min<double> pq
val- @(2.0, 2) = pqueue_extract_min<double> pq
val- @(3.0, 3) = pqueue_extract_min<double> pq
val- @(4.0, 4) = pqueue_extract_min<double> pq
val- @(5.0, 5) = pqueue_extract_min<double> pq
val- true = pqueue_is_empty pq


(*------------------------------------------------------------------*)
(*------------------------------------------------------------------*)
Line 1,047: Line 957:
end
end
end
end

(*------------------------------------------------------------------*)
(* A binary-heap priority queue, similar to the Pascal in Robert
Sedgewick, "Algorithms", 2nd ed. (reprinted with corrections),
1989. Note that Sedgewick does an extract-max, whereas we do an
extract-min.
Niklaus Wirth, within the heapsort implementation of "Algorithms +
Data Structures = Programs", has, I will note, some Pascal code
that is practically the same as Sedgewick's. Can we trace that code
back farther to Algol?
We do not have "goto" for Sedgewick's "downheap" (or Wirth's
"sift"), but do have mutual tail call as an obvious alternative to
the "goto". Nevertheless, because the code "jumped to" is small, I
simply use a macro to duplicate it. *)

typedef pqueue (priority_t : t@ype+,
value_t : t@ype+,
n : int,
n_max : int) =
[n <= n_max]
@{
arr = arrayref (@(priority_t, value_t), n_max + 1),
n = size_t n,
n_max = size_t n_max
}

prfn
lemma_pqueue_param
{n_max : int}
{n : int}
{priority_t, value_t : t@ype}
(pq : pqueue (priority_t, value_t, n, n_max))
:<prf> [0 <= n; n <= n_max] void =
lemma_g1uint_param (pq.n)

extern praxi
lemma_pqueue_size
{n_max : int}
{n : int}
{priority_t, value_t : t@ype}
(pq : pqueue (priority_t, value_t, n, n_max))
:<prf> [n1 : int | n1 == n] void

extern fn {priority_t : t@ype}
pqueue$cmp :
(priority_t, priority_t) -<> int

extern fn {priority_t : t@ype}
pqueue$priority_min :
() -<> priority_t

implement pqueue$cmp<double> (x, y) = compare (x, y)
implement pqueue$priority_min<double> () = neg infinity

fn {priority_t : t@ype}
{value_t : t@ype}
pqueue_make_empty
{n_max : int}
(n_max : size_t n_max,
arbitrary_entry : @(priority_t, value_t))
:<!wrt> pqueue (priority_t, value_t, 0, n_max) =
let
(* Currently an array is allocated whose size is the proven
bound. It might be better to use a smaller array and allow
reallocation up to this maximum size, or to break the array
into pieces. *)

prval () = lemma_g1uint_param n_max
val arr =
arrayref_make_elt<@(priority_t, value_t)>
(succ n_max, arbitrary_entry)
in
@{arr = arr,
n = i2sz 0,
n_max = n_max}
end

fn {}
pqueue_clear
{n_max : int}
{n : int}
{priority_t : t@ype}
{value_t : t@ype}
(pq : &pqueue (priority_t, value_t, n, n_max)
>> pqueue (priority_t, value_t, 0, n_max))
:<!wrt> void =
let
prval () = lemma_g1uint_param (pq.n_max)
in
pq := @{arr = pq.arr,
n = i2sz 0,
n_max = pq.n_max}
end

fn {}
pqueue_is_empty
{n_max : int}
{n : int}
{priority_t : t@ype}
{value_t : t@ype}
(pq : pqueue (priority_t, value_t, n, n_max))
:<> bool (n == 0) =
(pq.n) = i2sz 0

fn {}
pqueue_size
{n_max : int}
{n : int}
{priority_t : t@ype}
{value_t : t@ype}
(pq : pqueue (priority_t, value_t, n, n_max))
:<> size_t n =
pq.n

fn {priority_t : t@ype}
{value_t : t@ype}
_upheap {n_max : pos}
{n : int | n <= n_max}
{k0 : nat | k0 <= n}
(arr : arrayref (@(priority_t, value_t), n_max + 1),
k0 : size_t k0)
:<!refwrt> void =
let
macdef lt (x, y) = (pqueue$cmp<priority_t> (,(x), ,(y)) < 0)
macdef prio x = ,(x).0

val entry = arr[k0]

fun
loop {k : nat | k <= n}
.<k>.
(k : size_t k)
:<!refwrt> void =
if k = i2sz 0 then
arr[k] := entry
else
let
val kh = half k
in
if (prio entry) \lt (prio arr[kh]) then
begin
arr[k] := arr[kh];
loop kh
end
else
arr[k] := entry
end
in
arr[0] := @(pqueue$priority_min<priority_t> (), arr[0].1);
loop k0
end

fn {priority_t : t@ype}
{value_t : t@ype}
pqueue_insert
{n_max : int}
{n : int | n < n_max}
(pq : &pqueue (priority_t, value_t, n, n_max)
>> pqueue (priority_t, value_t, n + 1, n_max),
entry : @(priority_t, value_t))
:<!refwrt> void =
let
prval () = lemma_g1uint_param (pq.n)
val arr = pq.arr
and n1 = succ (pq.n)
in
arr[n1] := entry;
_upheap {n_max} {n + 1} (arr, n1);
pq := @{arr = arr,
n = n1,
n_max = pq.n_max}
end

fn {priority_t : t@ype}
{value_t : t@ype}
_downheap {n_max : pos}
{n : pos | n <= n_max}
(arr : arrayref (@(priority_t, value_t), n_max + 1),
n : size_t n)
:<!refwrt> void =
let
macdef lt (x, y) = (pqueue$cmp<priority_t> (,(x), ,(y)) < 0)
macdef prio x = ,(x).0

val entry = arr[1]
and nh = half n

fun
loop {k : pos | k <= n}
.<n - k>.
(k : size_t k)
:<!refwrt> void =
let
macdef move_data i =
if (prio entry) \lt (prio arr[,(i)]) then
arr[k] := entry
else
begin
arr[k] := arr[,(i)];
loop ,(i)
end
in
if nh < k then
arr[k] := entry
else
let
stadef j = 2 * k
prval () = prop_verify {j <= n} ()
val j : size_t j = k + k
in
if j < n then
let
stadef j1 = j + 1
prval () = prop_verify {j1 <= n} ()
val j1 : size_t j1 = succ j
in
if ~((prio arr[j]) \lt (prio arr[j1])) then
move_data j1
else
move_data j
end
else
move_data j
end
end
in
loop (i2sz 1)
end

fn {priority_t : t@ype}
{value_t : t@ype}
pqueue_peek
{n_max : int}
{n : pos | n <= n_max}
(pq : pqueue (priority_t, value_t, n, n_max))
:<!ref> @(priority_t, value_t) =
let
val arr = pq.arr
in
arr[1]
end

fn {priority_t : t@ype}
{value_t : t@ype}
pqueue_delete
{n_max : int}
{n : pos | n <= n_max}
(pq : &pqueue (priority_t, value_t, n, n_max)
>> pqueue (priority_t, value_t, n - 1, n_max))
:<!refwrt> void =
let
val @{arr = arr,
n = n,
n_max = n_max} = pq
in
if i2sz 0 < pred n then
begin
arr[1] := arr[n];
_downheap {n_max} {n - 1} (arr, pred n)
end;
pq := @{arr = arr,
n = pred n,
n_max = n_max}
end

fn {priority_t : t@ype}
{value_t : t@ype}
pqueue_extract
{n_max : int}
{n : pos | n <= n_max}
(pq : &pqueue (priority_t, value_t, n, n_max)
>> pqueue (priority_t, value_t, n - 1, n_max))
:<!refwrt> @(priority_t, value_t) =
let
val retval = pqueue_peek<priority_t><value_t> {n_max} {n} pq
in
pqueue_delete<priority_t><value_t> {n_max} {n} pq;
retval
end

local (* A little unit testing of the priority queue
implementation. *)
#define NMAX 10
in
var pq = pqueue_make_empty<double><int> (i2sz NMAX, @(0.0, 0))
val- true = pqueue_is_empty pq
val- true = (pqueue_size pq = i2sz 0)
val () = pqueue_insert (pq, @(3.0, 3))
val () = pqueue_insert (pq, @(5.0, 5))
val () = pqueue_insert (pq, @(1.0, 1))
val () = pqueue_insert (pq, @(2.0, 2))
val () = pqueue_insert (pq, @(4.0, 4))
val- false = pqueue_is_empty pq
val- true = (pqueue_size pq = i2sz 5)
val- @(1.0, 1) = pqueue_extract<double> pq
val- @(2.0, 2) = pqueue_extract<double> pq
val- @(3.0, 3) = pqueue_extract<double> pq
val- @(4.0, 4) = pqueue_extract<double> pq
val- @(5.0, 5) = pqueue_extract<double> pq
val- true = pqueue_is_empty pq
val- true = (pqueue_size pq = i2sz 0)
end


(*------------------------------------------------------------------*)
(*------------------------------------------------------------------*)
Line 1,060: Line 1,274:
(* Returns total-costs and previous-hops arrays. *)
(* Returns total-costs and previous-hops arrays. *)
:<!refwrt> @(arrayref (flt, n),
:<!refwrt> @(arrayref (flt, n),
arrayref ([i : nat | i <= n] size_t i, n)) =
arrayref ([i : nat | i <= n] size_t i, n)) =
let
let
prval () = lemma_matrixref_param adj_matrix
prval () = lemma_matrixref_param adj_matrix
Line 1,067: Line 1,281:
typedef defined_index_t = [i : nat | i < n] size_t i
typedef defined_index_t = [i : nat | i < n] size_t i
val index_t_undefined : size_t n = n
val index_t_undefined : size_t n = n

val arbitrary_pq_entry : @(flt, defined_index_t) =
@(0.0, i2sz 0)


val prev = arrayref_make_elt<index_t> (n, index_t_undefined)
val prev = arrayref_make_elt<index_t> (n, index_t_undefined)
and cost = arrayref_make_elt<flt> (n, infinity)
and cost = arrayref_make_elt<flt> (n, infinity)
val () = cost[source] := zero
val () = cost[source] := zero

(* The priority queue never gets larger than m_max. There is code
below that proves this; thus there is no risk of overrunning
the queue's storage (unless the queue implementation itself is
made unsafe). FIXME: Is it possible to prove a tighter bound on
the size of the priority queue? *)
stadef m_max = (n * n) + n + n
prval () = mul_pos_pos_pos (mul_make {n, n} ())
prval () = prop_verify {n + n < m_max} ()
val m_max : size_t m_max = (n * n) + n + n


typedef pqueue_t (m : int) =
typedef pqueue_t (m : int) =
[0 <= m] pqueue (flt, defined_index_t, m)
[0 <= m; m <= m_max]
pqueue (flt, defined_index_t, m, m_max)
typedef pqueue_t =
typedef pqueue_t =
[m : int] pqueue_t m
[m : int] pqueue_t m
var pq : pqueue_t = pqueue_make_empty ()


fn
pq_make_empty ()
:<!wrt> pqueue_t 0 =
(* Create a priority queue, whose maximum size is our proven
upper bound on the queue size. *)
pqueue_make_empty<flt><defined_index_t>
(m_max, arbitrary_pq_entry)

var pq = pq_make_empty ()
val active = arrayref_make_elt<bool> (n, true)
val active = arrayref_make_elt<bool> (n, true)
var num_active : [i : nat | i <= n] size_t i = n
var num_active : [i : nat | i <= n] size_t i = n
Line 1,084: Line 1,320:
fill_pq {i : nat | i <= n}
fill_pq {i : nat | i <= n}
.<n - i>.
.<n - i>.
(pq : &pqueue_t >> _,
(pq : &pqueue_t i >> pqueue_t n,
i : size_t i)
i : size_t i)
:<!refwrt> void =
:<!refwrt> void =
if i <> n then
if i <> n then
begin
begin
pqueue_insert (pq, cost[i], i);
pqueue_insert {m_max} {i} (pq, @(cost[i], i));
fill_pq (pq, succ i)
fill_pq {i + 1} (pq, succ i)
end
end


fun
fun
extract_min
extract_min
{m0 : nat}
{m0 : pos | m0 + n <= m_max}
.<m0>.
.<m0>.
(pq : &pqueue_t m0 >> pqueue_t m1)
(pq : &pqueue_t m0 >> pqueue_t m1)
:<!refwrt> #[m1 : nat | m1 <= m0]
:<!refwrt> #[m1 : nat | m1 < m0]
Option @(flt, defined_index_t) =
@(flt, defined_index_t) =
let
if pqueue_is_empty pq then
None ()
val @(priority, vertex) =
pqueue_extract<flt><defined_index_t> {m_max} {m0} pq
else
let
in
val @(priority, vertex) =
if active[vertex] then
@(priority, vertex)
pqueue_extract_min<flt><defined_index_t> pq
else if pqueue_is_empty {m_max} pq then
in
if ~active[vertex] then
arbitrary_pq_entry
extract_min pq
else
else
extract_min pq
end
Some @(priority, vertex)
end


fun
fun
main_loop {num_active0 : nat | num_active0 <= n}
main_loop {num_active0 : nat | num_active0 <= n}
{qsize0 : nat}
{qlimit0 : int | 0 <= qlimit0;
qsize0 <= qlimit0 + n}
.<num_active0>.
.<num_active0>.
(pq : &pqueue_t >> pqueue_t,
(* The pf_qlimit0 prop helps us prove a bound on the
size of the priority queue. We need it because the
proof capabilities built into ATS have very limited
ability to handle multiplication. *)
(pf_qlimit0 : MUL (n - num_active0, n, qlimit0) |
pq : &pqueue_t qsize0 >> pqueue_t 0,
num_active : &size_t num_active0 >> size_t num_active1)
num_active : &size_t num_active0 >> size_t num_active1)
:<!refwrt> #[num_active1 : nat | num_active1 <= num_active0]
:<!refwrt> #[num_active1 : nat | num_active1 <= num_active0]
void =
void =
if num_active <> i2sz 0 then
if num_active = i2sz 0 then
case+ extract_min pq of
pqueue_clear pq
| None () =>
else if pqueue_is_empty {m_max} {qsize0} pq then
(* This branch should never be reached, but the
let (* This should not happen. *)
val- false = true
corresponding sanity check is done further below. Thus
in
the proof of termination of this loop does not depend on
an exception. *)
end
()
else
| Some @(priority, u) =>
let
let
prval () = mul_elim pf_qlimit0
val () = active[u] := false
prval () =
and () = num_active := pred num_active
prop_verify {qsize0 <= ((n - num_active0) * n) + n} ()
prval () = mul_compare_lt {n - num_active0, n, n} ()
prval () = prop_verify {qsize0 < m_max} ()


fun
val @(priority, u) = extract_min pq
prval [qsize : int] () = lemma_pqueue_size {m_max} pq
loop_over_neighbors
{v : nat | v <= n}
prval () = lemma_pqueue_param {m_max} {qsize} pq
.<n - v>.
prval () = prop_verify {qsize < qsize0} ()
(pq : &pqueue_t >> pqueue_t,
prval () = prop_verify {qsize < m_max} ()
v : size_t v)
:<!refwrt> void =
if v = n then
()
else if ~active[v] then
loop_over_neighbors (pq, succ v)
else
let
val alternative = cost[u] + adj_matrix[u, n, v]
in
if alternative < cost[v] then
begin
cost[v] := alternative;
prev[v] := u;


(* Rather than lower the priority of v, this
val () = active[u] := false
and () = num_active := pred num_active
implementation inserts a new entry for v and
ignores obsolete queue entries. Queue entries
are obsolete if the vertex's entry in the
"active" array is false. *)


fun
pqueue_insert<flt><defined_index_t>
(pq, alternative, v)
loop_over_vertices
end;
{v : nat | v <= n}
loop_over_neighbors (pq, succ v)
{m0 : nat | qsize <= m0; m0 <= qsize + v}
end
.<n - v>.
(pq : &pqueue_t m0 >> pqueue_t m1,
in
loop_over_neighbors (pq, i2sz 0);
v : size_t v)
:<!refwrt> #[m1 : int | qsize <= m1; m1 <= qsize + n]
main_loop {num_active0 - 1} (pq, num_active)
end
void =
if v = n then
()
else if ~active[v] then
loop_over_vertices {v + 1} {m0} (pq, succ v)
else
let
val alternative = cost[u] + adj_matrix[u, n, v]
in
if alternative < cost[v] then
let
prval () = prop_verify {m0 < m_max} ()
in
cost[v] := alternative;
prev[v] := u;


(* Rather than lower the priority of v, this
val () = fill_pq (pq, i2sz 0)
implementation inserts a new entry for v and
val () = main_loop {n} (pq, num_active)
ignores obsolete queue entries. Queue entries
are obsolete if the vertex's entry in the
"active" array is false. *)


pqueue_insert<flt><defined_index_t>
(* Sanity check. *)
{m_max} {m0}
val- true = iseqz num_active
(pq, @(alternative, v));

loop_over_vertices {v + 1} {m0 + 1}
(pq, succ v)
end
else
loop_over_vertices {v + 1} {m0} (pq, succ v)
end

val () = loop_over_vertices {0} {qsize} (pq, i2sz 0)
in
main_loop {num_active0 - 1}
(MULind pf_qlimit0 | pq, num_active)
end
in
in
fill_pq {0} (pq, i2sz 0);
main_loop {n} {n} (MULbas () | pq, num_active);
@(cost, prev)
@(cost, prev)
end
end
Line 1,276: Line 1,537:


{{out}}
{{out}}
<pre>$ patscc -DATS_MEMALLOC_GCBDW dijkstra-algorithm.dats -lgc && ./a.out
<pre>$ patscc -O3 -DATS_MEMALLOC_GCBDW dijkstra-algorithm.dats -lgc && ./a.out
The requested paths:
The requested paths:
a -> c -> d -> e (cost = 26.000000)
a -> c -> d -> e (cost = 26.000000)
Retrieved from "https://rosettacode.org/wiki/Dijkstra%27s_algorithm"