CLPFD and recursive formulations

Question

Given a field, find the "cheapest" path from (0, 0) to (3, 3). Only the moves up and to the right are allowed.

field([
   [1, 1, 1, 1],
   [3, 3, 4, 3],
   [3, 3, 1, 3],
   [3, 3, 1, 1] 
]).

nth([A|_], 0, A).
nth([_|T], I, A) :- I #> 0, J #= I - 1, nth(T, J, A).

at(I, J, E) :- field(M), nth(M, I, R), nth(R, J, E).

Here's my recursive solution:

:- use_module(library(clpfd)).

price(0, 0, P) :- at(0, 0, P).
price(I, 0, P) :-
  I #=< 3,
  I #> 0,
  I0 #= I - 1,
  price(I0, 0, P0),
  at(I, 0, P1),
  P #= P0 + P1.
price(0, J, P) :-
  J #=< 3,
  J #> 0,
  J0 #= J - 1,
  price(0, J0, P0),
  at(0, J, P1),
  P #= P0 + P1.
price(I, J, P) :-
  I #=< 3,
  J #=< 3,
  I #> 0,
  J #> 0,
  I0 #= I - 1,
  J0 #= J - 1,
  price(I0, J, P0),
  price(I, J0, P1),
  at(I, J, P2),
  P #= min(P0, P1) + P2.

I can run it in different directions:

?- price(3, 3, P).
?- price(I, J, P), I #= 3, J #= 3.
?- price(I, J, P), P #> 10.

1. Is it possible to memorize the intermediate steps, and is it strictly necessary?
2. Is there a way to use memorization only when I and J are grounded?
3. How does CLPFD deal with the fact that there are 4 rules for the predicate price? Does the library merge the different rules into one big if?
4. Is it better to write a single rule price with if_'s inside?

Answer 1

Ad 3, current implementations of clpfd/clpz are all a library that has very few connections to the actual compilation process ("preparation for execution"). So there is no such predicate-wise optimization that you are hoping for. The clauses remain and are tried out one after the other — maybe with a little bit indexing due to the 0es in the head. If you want to get efficient code, you need to combine the common cases directly. E.g., you might first consider in a single rule I in 1..3, J in 1..3 and then avoid all those other comparisons. This means you would need some intermediary auxiliary predicate.

Also note the fine difference between I #=< 3, I #> 0 and I in 1..3. The former succeeds for I = 1+0 while the latter produces a type error. And for 0+0 OTOH, your program fails. You probably want at least consistent behaviour wherefore you may also use (#)/1 in places where general expressions are expected.

Ad 4, if_/3 or rather reif might be of interest, because you have just the special case of 0. Also, please note that often in stead of nesting if_/3 you may use a more complex reified condition as here on page 8 in treememberd_t/3.

Ad 1&2 there is some support for this in some systems, but it still is very bleeding edge. In particular if you want to use it together with constraints.

---

One final remark about your code in general. What you did was to take one-by-one an existing algorithm and you recoded it directly. Think of the last clause where you explore two paths independently and then join their results. You only can do this, because you know a lot about the actual problem, in particular both paths will have a solution. The only case were actual non-determinism happens is the instantiation of the coordinates.

Answer 2

As I show in my other answer here your recursive solution behaves exponentially. The usual approach to counter this is to add memoization. But memoized recursion's results are straightforwardly achieved, by flipping the time arrow, with a dynamic programming approach which just creates the field of resulting optimal prices from the given field right away; that field of prices can then be freely queried using nth or its built-in equivalents.

In the following I will assume field is a fully present NxN list structure. The elements of these lists might or might not be instantiated.

% field(Field), prices(Field,Prices).
prices([],[]).
prices([R1|Rs],[Q1|Qs]):-
   a(R1,0,Q1),  
   bs(Rs,Q1,Qs). 

a([],_,[]).
a([A|AT],P0,[P|PT]):- P #= P0+A, a(AT,P,PT).

bs([],_,[]).
bs([R|Rs],[P0|Q0],[Q|Qs]):- b(R,P0,[P0|Q0],Q), bs(Rs,Q,Qs).

b([],_,[],[]).
b([A|AT],P01,[P0|PT0],[P|PT]):- P #= min(P0,P01)+A, b(AT,P,PT0,PT).

Testing:

69 ?- time(( field(_F),maplist(writeln,_F),writeln('----'),
             prices(_F,_X),maplist(writeln,_X) )).
[1,1,1,1,1,1]
[3,3,4,3,2,5]
[3,3,1,3,1,3]
[3,3,1,1,4,3]
[3,3,1,1,4,2]
[3,3,7,1,4,1]
----
[1,  2, 3, 4, 5, 6]
[4,  5, 7, 7, 7,11]
[7,  8, 8,10, 8,11]   % (hand formatted)
[10,11, 9,10,12,14]
[13,14,10,11,15,16]
[16,17,17,12,16,17]
% 289 inferences, 0.000 CPU in 0.027 seconds (0% CPU, Infinite Lips)
true.

This doesn't literally answer any of your four questions, but it shows how the memoized recursion can be coded with a dynamic programming approach.

Answer 3

To the last half of your question 1, of course memoization is worth it here, your code clocks as exponential instead of quadratic (which obviates your 2, I think):

30 ?- time( (I=5, J=5, price(I, J, P)) ).
% 19,120 inferences, 0.000 CPU in 0.007 seconds (0% CPU, Infinite Lips)
I = J, J = 5,
P = 17 ;
% 10,163 inferences, 0.000 CPU in 0.002 seconds (0% CPU, Infinite Lips)
false.

31 ?- time( (I=4, J=4, price(I, J, P)) ).
% 5,120 inferences, 0.000 CPU in 0.003 seconds (0% CPU, Infinite Lips)
I = J, J = 4,
P = 15 ;
% 2,771 inferences, 0.000 CPU in 0.001 seconds (0% CPU, Infinite Lips)
false.

32 ?- time( (I=3, J=3, price(I, J, P)) ).
% 1,366 inferences, 0.000 CPU in 0.000 seconds (?% CPU, Infinite Lips)
I = J, J = 3,
P = 10 ;
% 769 inferences, 0.000 CPU in 0.000 seconds (?% CPU, Infinite Lips)
false.

33 ?- _X = [ 19.12/5.12, 5.12/1.37 ], maplist(is,A,_X).
A = [3.734375, 3.7372262773722627].

34 ?- _X = [ 29.3/7.9, 7.9/2.15 ], maplist(is,A,_X).
A = [3.7088607594936707, 3.674418604651163].

I ran the above with your code which I modified to enlarge the field, as follows:

field([
   [1, 1, 1, 1, 1, 1],
   [3, 3, 4, 3, 2, 5],
   [3, 3, 1, 3, 1, 3],
   [3, 3, 1, 1, 4, 3],
   [3, 3, 1, 1, 4, 2],
   [3, 3, 7, 1, 4, 1]
]).

price(0, 0, P) :- at(0, 0, P).
price(I, 0, P) :-
  I #=< 5,
  .... .
price(0, J, P) :-
  J #=< 5,
  .... .
price(I, J, P) :-
  I #=< 5,
  J #=< 5,
  .... .