« SweConsNet Workshop 2013 (Call for speakers) | Main | Gecode version 4.0.0 released »

A first look at B-Prolog

First, here's my B-Prolog page with a lot of CLP models.

B-Prolog is a Prolog implementation that, besides extensions such as CLP(FD), CLP(Set), and CLP(Bool), also shines with the following extensions: B-Prolog also has built in support for memoization (table) which make dynamic programming easier (see some examples here). It also has an almost seamless support for changing between CLP(FD), LP/MIP, and SAT solver.

However, here I have almost exclusively concentrated on the CLP(FD) extensions and with heavily use of foreach loops, list comprehensions, and arrays/matrices (since this is how I tend to think when modeling these kind of combinatorial problems). Most problems has been implemented in other C(L)P systems before, see Common constraint programming problems.

More info about B-Prolog

First model: SEND+MORE=MONEY

Model: send_more_money.pl

Let's start with one of the standard problems just to get a feeling for the syntax of B-Prolog 
sendmore(Digits) :-
    Digits = [S,E,N,D,M,O,R,Y],
    Digits :: 0..9,
    alldifferent(Digits),
    S #> 0, M #> 0,
                 1000*S + 100*E + 10*N + D
    +            1000*M + 100*O + 10*R + E
    #= 10000*M + 1000*O + 100*N + 10*E + Y,       
    labeling(Digits).
This is no different from most CLP(FD) systems. Note that there is no need to load a specific clpfd module since everything is loaded already (and B-Prolog currently don't support modules).

Here's a short explanation:
  • Digits :: 0..9: This means that all the elements in the list Digits must be in the domain between 0..9.
  • The specific CLP(FD) operators starts with #, e.g. #= and #>.
  • labeling(Digits): starts the search of the solution.

N-Queens: introducing list comprehensions

Model: queens.pl

Tthis N-Queens model show more of the differences between B-Prolog and standard Prolog with CLP(FD) support, namely it's use of list comprehensions. 
queens2(N, Q) :-
        length(Q, N),
        Q :: 1..N,
        Q2 @= [Q[I]+I : I in 1..N],
        Q3 @= [Q[I]-I : I in 1..N],
        alldifferent(Q),
        alldifferent(Q2),
        alldifferent(Q3),       
        labeling([ff],Q).
The extraction the two diagonals is done via list comprehensions (Q2 and Q3), i.e. 
  Q2 @= [Q[I]+I : I in 1..N],
  alldifferent(Q2),
This works since the special operator V @= [....] extracts the elements Q[I] (note the use of array index here).

Here is a small benchmark using the above encoding (predicate queens2 in the file). queens/2 is the same as the above but use alldistinct instead of alldifferent:
Both take very long time for N=500, but is much faster for N=499 and N=501. As we see, using alldifferent is faster than using alldistinct. The latter use a stronger consistency checking, but that don't help in this encoding: 
  For N=400:
    queens2/2: 0.62s 10 backtracks
    queens5/2: 0.71s 10 backtracks
    queens/2:  2.16s 1 backtrack
  For N=499:
    queens2/2: 0.76 1 backtracks
    queens5/2: 1.1s  1 backtracks
    queens/2:  4.168 0 backtracks
  For N=500: too slow 
  For N=501:
    queens5/2: 1.1s 1 backtracks
    queens2/2: 3.14s 1 backtrack
    queens/2:  3.524s 2 backtracks
  For N=1000:
    queens5/2: 8.9s 2 backtracks
    queens2/2: 24.2s 2 backtracks
    queens/2:  34.146s 0 backtracks
[I don't understand why N=500 is so hard when 499 and 501 is solved quite fast. No other CP solver/system I've tested show this behaviour.]

B-Prolog is not the only Prolog based system using foreach loops. See below for some discussion of the differences between B-Prolog and ECliPSe CLP's do-loop construct (also used in SICStus Prolog).

alphametic.pl: a fairly general alphametic solver

Model: alphametic.pl

Using list comprehensions, foreach loops and arrays makes it quite easy to implement a fairly general solver for this kind of alphametic problems. It's not completely general since it (still) use Prolog's feature of handling variables. 
go :-
    L = [[_S,E,N,_D],[M,O,_R,E],[M,O,N,E,_Y]],
    alphametic(L, Base, Res),
    writeln(Res).

alphametic(L,Base, Vars) :- 
    reverse(L,Rev),
    Rev = [Last|Sums],

    term_variables(L, Vars),
    Vars :: 0..Base-1,

    alldifferent(Vars),
    Vals #= sum([Val : S in Sums,[Val],calc(S,Base,Val)]),
    calc(Last,Base,Vals),
    foreach(S in Sums,S[1]#>0),

    labeling([ff,split], Vars).

calc(X,Base,Y) :-
    length(X,Len),
    Y #= sum([X[I]*Base**(Len-I) : I in 1..Len]).
The main work here is done with term_variables/2 which extracts the variables used, the sum/1 which sums all the terms, and calc that is really a special case of scalar product.

Sudoku solver: more lists/arrays

Model: sudoku.pl

Here is a Sudoku solver in B-Prolog: 
go :-
   time(solve(1)).

solve(ProblemName) :-
   problem(ProblemName, Board),
   print_board(Board),
   sudoku(3, Board),
   print_board(Board).

print_board(Board) :-
   N @= Board^length,
   foreach(I in 1..N, 
      (foreach(J in 1..N, [X],
         (X @= Board[I,J],
         (var(X) -> write('  _')
           ;
          format('  ~q', [X]))
         )),
      nl)),
    nl.

problem(1, [](
    [](_, _, 2, _, _, 5, _, 7, 9),
    [](1, _, 5, _, _, 3, _, _, _),
    [](_, _, _, _, _, _, 6, _, _),
    [](_, 1, _, 4, _, _, 9, _, _),
    [](_, 9, _, _, _, _, _, 8, _),
    [](_, _, 4, _, _, 9, _, 1, _),
    [](_, _, 9, _, _, _, _, _, _),
    [](_, _, _, 1, _, _, 3, _, 6),
    [](6, 8, _, 3, _, _, 4, _, _))).

alldifferent_matrix(Board) :-
   Rows @= Board^rows,
   foreach(Row in Rows, alldifferent(Row)),
   Columns @= Board^columns,
   foreach(Column in Columns, alldifferent(Column)).

sudoku(N, Board) :-
   N2 is N*N,
   array_to_list(Board,BoardVar),
   BoardVar :: 1..N2,

   alldifferent_matrix(Board),
   foreach(I in 1..N..N2, J in 1..N..N2, 
          [SubSquare],
          (SubSquare @= [Board[I+K,J+L] : 
                        K in 0..N-1, L in 0..N-1],
          alldifferent(SubSquare))
   ),
   labeling([ff,down], BoardVar).
One of the most interesting things here is the code for the sub squares (SubSquare): (SubSquare @= [Board[I+K,J+L] :
K in 0..N-1, L in 0..N-1]
Unlike most Prolog implementations, it works to access the entries in the matrix with I+K and J+L. However this works only when the indices (I, J, K, L) are ground integers. If any of the indices would be decision variables, then one have to use element/3 (or perhaps nth1/3) to obtain the same thing. See below for more about this.

Running this models yield the following result: 
  _  _  2  _  _  5  _  7  9
  1  _  5  _  _  3  _  _  _
  _  _  _  _  _  _  6  _  _
  _  1  _  4  _  _  9  _  _
  _  9  _  _  _  _  _  8  _
  _  _  4  _  _  9  _  1  _
  _  _  9  _  _  _  _  _  _
  _  _  _  1  _  _  3  _  6
  6  8  _  3  _  _  4  _  _

  3  6  2  8  4  5  1  7  9
  1  7  5  9  6  3  2  4  8
  9  4  8  2  1  7  6  3  5
  7  1  3  4  5  8  9  6  2
  2  9  6  7  3  1  5  8  4
  8  5  4  6  2  9  7  1  3
  4  3  9  5  7  6  8  2  1
  5  2  7  1  8  4  3  9  6
  6  8  1  3  9  2  4  5  7
Time and backtracks: Since there is no built-in predicate in B-Prolog for getting both time and number of backtracks (failures etc), let's add this: 
time2(Goal):-
   cputime(Start),
   statistics(backtracks, Backtracks1),
   call(Goal),
   statistics(backtracks, Backtracks2),
   cputime(End),
   T is (End-Start)/1000,
   Backtracks is Backtracks2 - Backtracks1,
   format('CPU time ~w seconds. Backtracks: ~d\n', [T, Backtracks]).
Running this on the command line: 
 $ time bp -g "[sudoku], time2(solve(1)),halt"
writes the following after the solution, i.e. that it's a 0 backtrack problem. 
CPU time 0.0 seconds. Backtracks: 0
In the model sudoku.pl there is a predicate go3/0 which runs all the 13 problem instances and shows the statistics (time and backtracks). The len:1 in the output is for ensuring that the solution is unique (a requirement on traditional Sudoku problems). This check caught two problem instances that was either typed in erroneously by me or by the source I got them from (and thus I removed these two problems from this model). 
go3 :-
   foreach(P in 1..13,
      [L,Len],
      (writeln(problem:P),
      time2(findall(P,solve2(P),L)),
      length(L,Len),
      writeln(len:Len),
      (Len > 1 ->
         writeln('This has more than one solution!') 
        ; 
      true
      ), nl)).
The result is: 
problem:1
CPU time 0.004 seconds. Backtracks: 0
len:1

....

problem:13
CPU time 0.04 seconds. Backtracks: 0
len:1

What you can - and cannot - do with list comprehensions

When using @= it works to get a list comprehension (provided that there's no problem with array access etc). If one want to use it directly in a (global) constraint, then the requirement is that it must be in an arithmetic context (I first missed that in the documentation).

The following works since sum/1 is a special function in an arithmetic context: 
  Sum #= sum([X[I] : I in 1..N])
However, using comprehensions in global constraints such as alldifferent don't work since it's not in an arithmetic context: 
  % this don't work
  alldifferent([X[I]+I : I in 1..N])
The following works (as we saw above in the N-Queens model) since we extract the elements to a list using @= before using it. 
  % this works
  Q1 @= [X[I]+I : I in 1..N],
  alldifferent(Q1),
  % ...
As mentioned above, array indexing don't work when the indices are decision variables. Then one must use element/3 (or nth/3). See more about element below.

foreach loops

The foreach loop extension i B-Prolog is quite intuitive. They are a often easier to use than the do-loops in ECLiPSe CLP (and SICStus Prolog) since in B-Prolog one declare the local variables in the loop (i.e. "global by default"), whereas in ECLiPSe do-loops one declares the global variables to be used ("local by default") which perhaps is not as intuitive. Both these systems have helpful warnings where a local/global variable is used out of context.

Note that B-Prolog's foreach construct don't have all features found in ECLiPSe's do-loop so certain things might be easier to state in ECLiPSe CLP and SICStus Prolog than in B-Prolog.

In B-Prolog there is also support for accumulators via ac(Var, StartValue) (or ac1/2 which I have not used much), which collects the local variables in the loop to a global variable (here Var).

An example of a foreach loop which use ac is in photo_problem.pl 
preferences(1, 7, [[1,5],
                   [1,6],
                   [2,1],
                   [2,5],
                   [4,6],
                   [4,3],
                   [7,4],
                   [7,3]]).

foreach([Pref1,Pref2] in Preferences,
        ac(Diffs,[]), [P1,P2,Diff],
        (
            element(Pref1,Positions,P1),
            element(Pref2,Positions,P2),
            Diff #= (abs(P1-P2) #= 1),
            Diffs^1 = [Diff|Diffs^0]
        )),
Here the Diff list collects the differences of the positions (from the local Diff variable). The local declarations [P1,P2,Diff] is needed here, otherwise they would be considered global variables in the model (with disastrous result).

Also note that the foreach loop handles the list of lists in Preferences, so the foreach is not limited to integer ranges that is shown in a couple of the other examples here.

More about foreach loops is found in the manual, section Declarative Loops and List Comprehensions.

new_array and lists

Creating an array (of "any" dimension) in B-Prolog is done with new_array(Var,Dimensions): 
   N = 3,
   new_array(X,[N,N]),
where Dimensions is a list of the dimensions (of "any" size; I don't know if there is any limitations of the number of dimensions).

The result of new_array is an array structure such as 
  | ?- new_array(X,[3,3])
  X = []([](_7a8,_7b0,_7b8),[](_7c8,_7d0,_7d8),[](_7e8,_7f0,_7f8))
Converting an array structure to a list is done with array_to_list(Array,List): 
| ?- new_array(X,[3,3]), array_to_list(X,List)
X = []([](_8e0,_8e8,_8f0),[](_900,_908,_910),[](_920,_928,_930))
List = [_8e0,_8e8,_8f0,_900,_908,_910,_920,_928,_930]
Setting the domains for the variables in an array must be done via the converted list (i.e. not directly on the array structure): 
   N = 3,
   new_array(X, [N,N]),
   array_to_list(X, Vars),
   Vars :: 1..N*N,
   % ...
After that conversion, one can freely use constraints on either the array structure (e.g. as a matrix) or on the list representation, such as alldifferent(Vars). This "dual representation" sometimes simplifies modeling much. Note that labeling must be done via the list representation, or by using term_variables/2.

Reifications: alldifferent_except_0

Model: alldifferent_except_0.pl

Using reifications (reasoning about satisfing constraints using boolean decision variables) is quite direct, using #= (equality), #=> (implication) and #<=> (equivalence). Here is an decomposition of the global constraint alldifferent_except_0.pl: 
alldifferent_except_0(Xs) :-
   Len @= Xs^length,
   foreach(I in 1..Len, J in 1..Len,
      (I #\=J #/\ Xs[I] #\= 0 #/\ Xs[J] #\= 0)
       #=> 
      (Xs[I] #\= Xs[J])
     ).
One thing I noticed when playing with reifications was that foreach don't seems be in boolean context. Example: this don't work as I expected (or hoped), i.e. that if not all X[I] are larger than 0 then B is 0 (and vice versa): 
  % don't work
  foreach(I in 1..N, X[I] #> 0) #<=> (B #= 1)

Table constraint

Model: hidato_table.pl

The table constraint is used in hidato_table.pl and is quite easy to use.

First, define the valid connections that can be used. This is - of course - done via a list comprehension: 
Connections
   @= [(I1,J1,I2,J2) :
       I1 in 1..N, J1 in 1..N, 
       I2 in 1..N, J2 in 1..N,
       (abs(I1-I2) =< 1,
       abs(J1-J2) =< 1,
       (I1 \=I2; J1 \= J2)
       )],
Then place all the numbers 1..N*N in the grid where the places are restriced by the connections in Connections. Note that since we are using decision variables for the indices (I1,J1,I2,J2) of the grid (X), we cannot use array extraction. My solution of this is to use the XVar list (which we must have for the labeling anyway), and calculate the corresponding position in this list for each (I,J) postition in the grid. 
XVar @= [X[I,J] : I in 1..N, J in 1..N],
XVar :: 1..N*N,

% ....

% place all integers from 1..N*N
alldifferent(XVar),
foreach(K in 1..(N*N)-1,
   [I1,J1,I2,J2,K2,Offset1,Offset2],
   [ac(AllConn,[]),ac(Offsets,[])],
   (
     K2 is K+1,
     [I1,J1,I2,J2] :: 1..N,
     % [(I1,J1,I2,J2)] in Connections,
     [Offset1,Offset2] :: 1..N*N,
     Offset1 #= (I1-1)*N+J1,
     element(Offset1,XVar,K),
     Offset2 #= (I2-1)*N+J2,
     element(Offset2,XVar,K2),
     AllConn^1 = [(I1,J1,I2,J2)|AllConn^0],
     Offsets^1 = [[Offset1,Offset2]|Offsets^0]
    )
   ),

% the table constraint
AllConn in Connections,
...
The table constraint is used by simply putting all the variables in a list of parenthesis structure (I1,J1,I2,J2) and then checking all these using AllConn in Connections.

Also, note that both AllConn (the collection of selected connections) as well as the offsets (collected in the Offsets list) is included in the labeling to get faster results. 
% ...
term_variables([Offsets,XVar,AllConn],Vars),
labeling([ff],Vars),
% ...
In my experiments with this model, it seems to be faster to put Offsets before XVar and AllConn in the labeling.

This model is not blazingly fast, though it's much faster than the model without the table constraint (the model hidato.pl). For example, problem instance #6 takes 6.5 seconds without the table constraint, and 0.17 seconds with (both versions have 0 backtracks).

Note: B-Prolog also have a powerful built-in memoisation feature called table. Please don't confuse these two usages of "table".

Element

As usual, when porting models from my MiniZinc models, I got the most problem with the element constraint. (Porting to B-Prolog has not been the worst case in this matter, though.)

B-Prolog's support for arrays with decision variables is not as good as compared with the ECLiPSe CLP system, and both have quite a way to go before it's as nice as MiniZinc's and Comet's support (where element is encoded very natural). It's especially when using element in a matrix where I had to stop in my "porting flow" (often the port was quite "flowy") and then use element/nth1 constraints and often a few temporary variables.

When there are much matrices with decision variable indices, I tend to use these two predicates instead (and skipping B-Prolog's arrays): 
matrix_element(X, I, J, Val) :-
   element(I, X, Row),
   element(J, Row, Val).
 
matrix(_, []) :- !.
matrix(L, [Dim|Dims]) :-
   length(L, Dim),
   foreach(X in L, matrix(X, Dims)).
(The latter predicate was suggested by Mats Carlsson when I implemented a lot of SICStus Prolog models some years ago.)

Here is an example in the stable marriage model (stable_marriage.pl).

In MiniZinc one can code element within decision variables like this (where both husband and wife are lists of decision variables). Here's a snippet from the MiniZinc model stable_marriage.mzn: 
forall(m in 1..num_men) ( 
  husband[wife[m]] = m 
)
In B-Prolog, I use the following (using the same overall approach): 
foreach(M in 1..NumMen,[WifeM,HusbandWifeM],
   (element(M,Wife,WifeM),
     element(WifeM,Husband,HusbandWifeM),
     HusbandWifeM #= M )),
Another example is the Five element problem (from Charles W. Trigg): five_elements.pl

The problem statement:
Charles W. Trigg, PI MU EPSILON JOURNAL, Valume 6, Fall 1977, Number 5
"""
From the following square array of the first 25 positive integers, choose five, no two of the same row or column, so that the maximum of the five elements is as small as possible.

2 13 16 11 23
15 1 9 7 10
14 12 21 24 8
3 25 22 18 4
20 19 6 5 17

Ensuring the unicity of rows and columns is simple: 
foreach(I in 1..N, 
       (sum([X[I,J] : J in 1..N]) #= 1,
        sum([X[J,I] : J in 1..N]) #= 1
        ))
However, we also want to find the specific numbers in the matrix for which X[I,J] #= 1. The following way do not work: 
% This don't work
foreach(J in 1..N, [I],ac(Is,[]),
   (I :: 1..N,
    1 #= X[I,J],
    Y[J] #= Matrix[I,J]
   ))
If we - still - want to use this matrix approach, then freeze/2 can be used. However, we then have to collect the I indices in an accumulator list Is to be used in the labeling: 
% Find the specific row I for which X[I,J] is 1.
foreach(J in 1..N, [I],ac(Is,[]),
    (I :: 1..N,
    freeze(I, 1 #= X[I,J]),
    freeze(I, Y[J] #= Matrix[I,J]),
    Is^1 = [I|Is^0])
    ),

MaxY #= max(Y),

term_variables([XVar,Y,MaxY,Is], Vars),
minof(labeling([ff], Vars),MaxY),
labeling(Vars),
% ....
Using freeze/2 is not ideal, but it makes the code more intuitive than with a lot of element/3 or nth1/3, if that even work.

MIP/SAT

Model: coins_grid.pl
As mentioned in the introduction above, B-Prolog also have support for LP/MIP problems (using GLPK as the solver) and an interface to a SAT solver.

Here I just show the MIP solver by one of my standard problems for which CP solvers tend to be worse than MIP solvers, namely coins_grid.pl. The model also implements a CLP(FD) and a cp_solve.

Some differences between the other CLP(FD) models shown here:
  • The CP, MIP, and SAT solvers has a little different syntax than the CLP(FD) solver:
    • The numeric constraints use $= instead of #= (i.e. "$" instead of "#") to mark that it's a CLP constraint.
    • The labeling is different: cp_solve, lp_solve, ip_solve, and sat_solve.
  • There are just a few global constraint that is supported using this: $alldifferent and $element, and the usual arithmetic sum/1, min/1, max/1, min/2, max/2. The CLP(FD) solver has support for many more global constraints. 
coins(N, C) :-
   % N = 10, % 31 the grid size
   % C = 6,  % 14, number of coins per row/column
   
   new_array(X, [N,N]),
   array_to_list(X, Vars),
   Vars :: 0..1,

   Sum $>= 0,
   foreach(I in 1..N, 
       (C $= sum([T : J in 1..N, [T], T @= X[I,J]]), % rows
        C $= sum([T : J in 1..N, [T], T @= X[J,I]]) % columns
       )),

   % quadratic horizontal distance
   Sum $= sum([
               T : I in 1..N, J in 1..N, [T],
               T @= (X[I,J] * abs(I-J)*abs(I-J))
              ]),

   ip_solve([min(Sum),ff,reverse_split], Vars),
   % cp_solve([min(Sum),ff,reverse_split], Vars),
   % sat_solve([min(Sum),ff,reverse_split], Vars),
   writeln(sum:Sum),
   pretty_print(X).
By changing ip_solve to cp_solve (or sat_solve) we use the CP solver (or SAT solver). This is quite nifty.

And as usual the MIP solver solves this problem immediately, whereas CP solvers tend to be much slower.

Things not tested or mentioned

Here are some other things that I have not mentioned above:
  • Action rules and events. These are the used to implement propagators and more effective constraints. However, I have not played with these much. Also, see Programming Constraint Propagators.
  • B-Prolog is not an open source project, i.e. the free distribution just contains the executable, documentation and examples. (This was a little unfortunately since then I had not opportunity to study certain features, e.g. how the global constraints where implemented in Action rules.) For individual, academic and non-commercial use, B-Prolog is free to use. There is commercial licenses where source code is included. See more here: Licenses.
  • And then there are other features that I have not mentioned. See the Manual for detailed descriptions about these.

Summary

I like B-Prolog, especially the nice support for list comprehensions, array access and foreach loops. They are a nodge neater than the one used in ECLiPSe CLP and SICStus Prolog, though not as good as in MiniZinc (or Comet).

It was quite easy to port my ECLiPSe and SICStus Prolog models to B-Prolog, mostly because these already used the do-loop constructs, and it was mostly easy to port directly from my MiniZinc models, except when using matrices with decision variables as indices. Right now I think that B-Prolog's approach is often neater to use than the other Prolog's, especially list comprehensions.

Picat

Finally, I must mention Picat (Predicates, Imperative, Constraints, Actors, and Tabling). Picat is a C(L)P language also created by B-Prolog's creator Neng-Fa Zhou (together with Jonathan Fruhman), which is inspired by Prolog (especially B-Prolog). From the Picat site:
Picat is a general-purpose hybrid language that incorporates many declarative language features for better productivity of software development, including explicit non-determinism, explicit unification, functions, constraints, and tabling. Picat also provides imperative language constructs for programming everyday things. The Picat system will be used for not only symbolic computations, such as knowledge engineering, NLP, and search problems, but also for scripting tasks for the Web, games, and mobile applications.

....

The Picat implementation will be based on the B-Prolog engine and the first version that has the basic functionality is expected to be released by May, 2013. The project is open to anybody and you are welcome to join, as a developer, a sponsor, a user, or a reviewer. Please contact picat@picat-lang.org .
This sounds very interesting, and I will definitely give it a try when it's released (and perhaps blog about it as well).

My B-Prolog models/programs

Here are my public B-Prolog encodings. As of writing it's over 200, most using CLP(FD), some use CLP(Set), and just a few are non-CP Prolog + foreach/list comprehension.

Posted by hakank on March 10, 2013 08:55 AM |