Floyd-Warshall algorithm in Prolog

[emiruz]

Floyd-Warshall is a well known and simple algorithm for finding all the shortest paths in a graph through dynamic programming. Here is the pseudo-code:

let dist be a |V| × |V| array, initialized to ∞ (infinity)
for each edge (u, v) do
    dist[u][v] ← w(u, v)  // The weight of the edge (u, v)
for each vertex v do
    dist[v][v] ← 0
for k from 1 to |V|
    for i from 1 to |V|
        for j from 1 to |V|
            if dist[i][j] > dist[i][k] + dist[k][j] 
                dist[i][j] ← dist[i][k] + dist[k][j]
            end if

I’d like an implementation in a Prolog sympathetic way: immutable, perhaps using tabling, whilst hopefully conserving the conciseness of the original.

Does anyone have any such code, or might you fancy a stab at it?

[kwon-young]

Hello,
Here is my version using a list of list for a matrix with fully ground elements, meaning that you have to copy the matrix every time you do an update.
I think this is kind of the worst implementation anyone can think of, therefore kind of a good baseline ?

It is immutable as requested but the conciseness is not quite as good as the original because of for loops, especially nested ones…

floyd_warshall(V, Graph, M4) :-
   length(M1, V),
   maplist(length_(V), M1),
   maplist(maplist(=(inf)), M1),
   foldl(edges, Graph, M1, M2),
   numlist(1, V, Vertices),
   foldl(vertices, Vertices, M2, M3),
   transitive(Vertices, M3, M4).

edges((U1-V1)-W, M1, M2) :-
   replace(U1, V1, M1, _, W, M2).

vertices(I, M1, M2) :-
   replace(I, I, M1, _, 0, M2).

transitive(Vertices, M1, M2) :-
   foldl(transitive(Vertices), Vertices, M1, M2).
transitive(Vertices, K, M1, M2) :-
   foldl(transitive(Vertices, K), Vertices, M1, M2).
transitive(Vertices, K, I, M1, M2) :-
   foldl(transitive_(K, I), Vertices, M1, M2).

transitive_(K, I, J, M1, M2) :-
   replace(I, J, M1, IJ1, IJ2, M2),
   nth2d(I, K, M1, IK),
   nth2d(K, J, M1, KJ),
   IKJ is IK + KJ,
   (  IJ1 > IKJ
   -> IJ2 = IKJ
   ;  IJ2 = IJ1
   ).

Here is a example query, (the graph comes from wikipedia):

?- graph(G), gtrace, floyd_warshall(4, G, M).
G = [2-1-4, 2-3-3, 1-3- -2, 3-4-2, 4-2- -1],
M = [[0, -1, -2, 0], [4, 0, 2, 4], [5, 1, 0, 2], [3, -1, 1, 0]].

prolog code with graph fact and helper predicates

graph([
   (2-1)-4,
   (2-3)-3,
   (1-3)-(-2),
   (3-4)-2,
   (4-2)-(-1)
]).

length_(Len, List) :-
   length(List, Len).

nth2d(I, J, M, V) :-
   nth1(I, M, Row),
   nth1(J, Row, V).

replace(I, L1, Old, New, L3) :-
   nth1(I, L1, Old, L2),
   nth1(I, L3, New, L2).

replace(I1, I2, M1, Old, New, M2) :-
   replace(I1, M1, Row1, Row2, M2),
   replace(I2, Row1, Old, New, Row2).

PS: Seriously, why do you have to do this to me, you can’t expect me not to spend my holiday not optimizing this like I did with dijkstra

[kwon-young]

By the way:

@emiruz

I’d like an implementation in a Prolog sympathetic way: immutable

Does using clpfd variables for tracking distances count as immutable ?

@emiruz

perhaps using tabling

I don’t think tabling will help here with the three nested for loops.

[peter.ludemann]

The Craft of Prolog by Richard O’Keefe, page 171 (except it says “Warshall’s algorithm” rather than “Floyd-Warshall algorithm”.
And there’s a discussion of the problem starting at page 167.

[j4n_bur53]

@emiruz

hopefully conserving the conciseness of the original

Does “in Prolog” allow using Prolog to develop a little DSL?

?- warshall(4,D).
D = ''(''(0, -1, -2, 0), ''(4, 0, 2, 4), ''(5, 1, 0, 2), ''(3, -1, 1, 0)).

Here is the DSL and the Warshall algorithm:

:- op(100, yf, []).
:- op(800, fx, new).
:- op(800, fx, let).
:- op(800, fx, if).

warshall(N, D) :-
   new D[N,N],
   (between(1,N,U), between(1,N,V), let D[U,V] = 999, fail; true),
   (edge(U,V,W), let D[U,V] = W, fail; true),
   (vertex(V), let D[V,V] = 0, fail; true),
   (between(1,N,K),
       between(1,N,I),
          between(1,N,J),
              let H = D[I,K] + D[K,J],
              (if D[I,J] > H -> let D[I,J] = H; true),
              fail; true).

And here are the helpers of the DSL:

new D[N,M] :-
   functor(D, '', N),
   D =.. [_|L],
   new2(L, M).

new2([], _).
new2([X|L], N) :-
   functor(X, '', N),
   new2(L, N).

let V = E :- var(V), !,
   let2(E,V).
let D[R,C] = E :-
   let2(E,V),
   arg(R, D, H),
   nb_setarg(C, H, V).

let2(D[R,C], V) :- !,
   arg(R, D, H),
   arg(C, H, V).
let2(E+F, R) :- !,
   let2(E, V),
   let2(F, W),
   R is V+W.
let2(V, V).

if E > F :-
   let2(E, V),
   let2(F, W),
   V > W.

[emiruz]

Thanks @kwon-young, I suspect the algorithm may need some re-thinking to be Prolog friendly, if it is at all possible. Yours is a functional version of the imperative code, and roughly O(V^5) I think.

I like the DSL @j4n_bur53 but the actual code looks imperative, mutable and non-backtrackable.

@peter.ludemann its the transitive closure algorithm. It is given on page 5 here. I think that the main difference is the min condition from above is replaced with a logical OR. O’Keefe accumulates results using a ordered set since the values are binary, otherwise it is similar to @kwon-young approach. I don’t think this works when distance is required because we cannot make the value of a member in a set implicit (e.g. if its in the set then true); there needs to be some kind of updating structure, or a reenvisioning of how the algorithm calculates.

[emiruz]

To minimise the pain of Aoc 2024 Day 16 in Prolog, I suspect I’ll end up just biting the bullet and using nb_setarg/3. I’m really loathed to implement Djisktra yet again.

[emiruz]

So I’m thinking of keeping distance state a bit like this. Say I have N vertices:

NN is N**2,
length(Dist0, NN),
maplist(=(10**5), Dist0),
Dist =.. [dist|Dist0].

I’ll pass Dist around. It can then be accessed using arg/3 or changed using setarg/3 or nb_setarg/3 (which I expect is faster). Pretty horrible, but needs must.

[jan]

@emiruz

maplist(=(10**5), Dist0),

Note that 10**5 is not evaluated. I’d consider using a variable as initial value.

@emiruz

Pretty horrible, but needs must.

Yes. It is probably as good as it can be in Prolog. My only other hope would be using mode directed tabling. Possibly you an tweak that to do what you want, but it will always be a lot more expensive.

[j4n_bur53]

@emiruz

I like the DSL @j4n_bur53 but the actual code looks imperative, mutable and non-backtrackable.

Well I made the example to challenge the relevance of the notion
“declarativity” and “logical pureness”. Thats an implementation details.
A client of the predicate warshall/2 doesn’t notice anything.

The contract is as follows:

Input: N, argument in warshall/2
Input: edge/3
Input: vertex/1
Output: D, argument in warshall/2

The contract has no immutable or backtracking criterias that
would be violated, my implementation doesn’t mutate anything,
the array is returned fresh, its also not the case it would

prevent backtracking over warshall/2, since backtracking over a
deterministic predicate anyway does nothing. Last but not least,
its is much faster than a list implementation, did you test it?

[j4n_bur53]

Ok you wrote:

@emiruz

To minimise the pain of Aoc 2024 Day 16 in Prolog, I suspect I’ll end up just biting the bullet and using nb_setarg/3. I’m really loathed to implement Djisktra yet again.

Is this based on performance measure or elegancy of code?
The idea to have predicates implemented with nb_setarg/3
and still have an innocent contract is not new.

There is/was a wave of interested researches
concering “persistent” data structures:

In computing, a persistent data structure or not ephemeral data
structure is a data structure that always preserves the previous version
of itself when it is modified. Such data structures are effectively immutable,
as their operations do not (visibly) update the structure in-place, but
instead always yield a new updated structure. The term was introduced
in Driscoll, Sarnak, Sleator, and Tarjan’s 1986 article.
https://en.wikipedia.org/wiki/Persistent_data_structure

It applies to a whole range of programming languages
such as functional Haskell etc…, and relational Prolog etc…
The idea is that using something destructive inside an implementation

could be a “green” case just like there are “green” cuts.

[emiruz]

@j4n_bur53

Is this based on performance measure or elegancy of code?
The idea to have predicates implemented with nb_setarg/3
and still have an innocent contract is not new.

Each to their own I suppose, but under those guises you might as well implement it in C.

[j4n_bur53]

@emiruz

Each to their own I suppose, but under those guises you might as well implement it in C.

Not really. What is the garbage collection of C? I mean I do have this here:

@j4n_bur53

new D[N,M] :-
   functor(D, '', N),
   D =.. [_|L],
   new2(L, M).

In C one would do malloc(). But where is a free() in my Prolog code?
The closest to what SWI-Prolog does under the hood is probably
stack allocate. But you cannot use stack allocate in C by the

warshall/2 realization, because you want to return a value. You would
most likely loose the memory when you return the routine. So the
overall C code usually becomes more clumsy, delegating the

allocation to the client.

[emiruz]

That rather depends on the SWI Prolog C API, but presumably it caters for that? Further, if memory is copied to Prolog at the end of the call then presumably all allocation can happen within the function. It largely depends on how the API works rather than whether its possible to free memory in an interop

[j4n_bur53]

Here is a performance comparison:

?- time((between(1,1000,_), graph(G), floyd_warshall(4, G, M), fail; true)).
% 3,803,998 inferences, 0.156 CPU in 0.183 seconds (85% CPU, 24345587 Lips)
true.

?- time((between(1,1000,_), warshall(4,D), fail; true)).
% 1,046,998 inferences, 0.062 CPU in 0.062 seconds (100% CPU, 16751968 Lips)
true.

The nb_setarg/3 solution is ca 3 times faster than the list solution.

The nb_setarg/3 solution although it is already faster, it could be
made even an inch faster, if the DSL were not interpreted but
compiled. Currently the DSL is interpreted like here:

@j4n_bur53

let2(D[R,C], V) :- !,
   arg(R, D, H),
   arg(C, H, V).
let2(E+F, R) :- !,
   let2(E, V),
   let2(F, W),
   R is V+W.
let2(V, V).

But one could expand such things during compile time,
making it execute again faster. A programming language
that is basically a compiled DSL is Picat, which is

compiled on top of B-Prolog.

[j4n_bur53]

Ok, DSL compilation gives again a factor 3 speed:

/* DSL interpreted */
?- time((between(1,1000,_), warshall(4,D), fail; true)).
% 1,046,998 inferences, 0.062 CPU in 0.062 seconds (100% CPU, 16751968 Lips)
true.

/* DSL compiled */
?- time((between(1,1000,_), warshall(4,D), fail; true)).
% 336,998 inferences, 0.000 CPU in 0.020 seconds (0% CPU, Infinite Lips)
true.

The compiler was coded by only changing the DSL interpretation
rules into some goal expansion rules. Make sure to read the
expansion rules before the predicate warshall/2,

otherwise it doesn’t get DSL compiled:

goal_expansion(let V = E,
   let2(E,V)) :- var(V), !.
goal_expansion(let D[R,C] = E,
  (let2(E,V),
   arg(R, D, H),
   nb_setarg(C, H, V))).

goal_expansion(let2(V, V), true) :- var(V), !.
goal_expansion(let2(D[R,C], V),
  (arg(R, D, H),
   arg(C, H, V))) :- !.
goal_expansion(let2(E+F, R),
  (let2(E, V),
   let2(F, W),
   R is V+W)) :- !.
goal_expansion(let2(V, V), true).

goal_expansion(if E > F,
  (let2(E, V),
   let2(F, W),
   V > W)).