☰

Constraint Logic Programming

Course Notes (lecture version)


These are the notes accompanying the slides Constraint Logic Programming (7_clp_nosem.pdf) of the Computational Logic course.

This is the lecture version: it follows the slides in order, but keeps much more of what was actually said in class — the asides, the reasons, the second explanations, and the things that were only ever said out loud. Where the slide gives a bullet, the lecture usually gave a paragraph and an example.

How to use this interactive document:

The code boxes, marked by a question mark (?) in the right top corner, allow interaction. Feel free to edit any program and make all the queries you want!

  • Code needs to be loaded before running, by pressing on the question mark — code box should be checked ✔.
  • Under each program there are editable boxes with a right-pointing triangle (▶) that allow interaction with the code (once loaded!) by issuing queries.
  • Expected answers are printed in plain text underneath each query box, so you can read the notes straight through without running anything.
  • When available, pressing the northeast pointing arrow (↗) will load the code in a separate Prolog playground window.

Part I — What constraints are, and why we want them

Slide 2 — Constraints

You may have seen constraints already in other areas. What we are trying to do here is to see the combination of constraints with logic and with programming — and that combination, which is what CLP is, is quite interesting and unique. It has characteristics that none of the components have by themselves.

Constraints for representing problems. The idea of constraint solving was born within artificial intelligence, for doing things like house design, or intelligent design in general. Constraints have been used ever since — in the same way as logic, and as a generalization of logic — to represent problems.

Take a puzzle, say a murder case, with rules like "the man in yellow does not have green eyes", or "the murderer knows that no detective will ever wear dark clothes" — and, as we shall see, quantitative conditions too.

Those statements are a set of constraints over the possible solutions: a description of the state space. And a solution is an assignment of values to variables. You may not see the variables in that sentence, but "the man" is a variable — we do not know his name, so we can assign one — and in particular "the murderer" is a variable.

And now the idea to hold on to. Sometimes we do not have enough information to answer with a particular value — to say "the murderer is López". But we can often still answer with a set of constraints: we can refine the space. The answer might be "the murderer is one of those who had met the entertainer" — more precise than what we started with, a consequence of it, but not solved.

The answer to a constraint problem can itself be a constraint.

And of course we must also cater for there being no solution, because a system can be inconsistent.

A CLP system is characterized by four things:

  • the domain of computation — reals, rationals, integers, booleans, structures, ...
  • the expressions that can be built over it — +, *, ∧, ∨
  • the constraints allowed — equations, inequations, disequations: =, \=, =<, >=, <, >
  • the constraint solving algorithms — simplex, Gauss, and others

Where this came from

There are a great many ancestors for constraint solving: very early algorithms like Sketchpad, which was for drawing; Waltz's algorithm; ThingLab; Macsyma.

Constraint programming — being able to use constraints inside a programming language — was born within logic programming, first as extensions of Prolog: natural extensions to fix the arithmetic, and to fix other issues such as infinite trees and the problems with unification and the occurs check.

Then a lot of these separate efforts were unified by a really interesting theory, which gave the field its name, developed by Jaffar and Lassez in 1987. From there standalone languages such as CLP(R) and CHIP were born, and were very popular for a while.

At first it seemed they had to be different languages from Prolog. Time has shown that the most natural thing is to embed constraints inside Prolog. What happens now is that many Prologs — really all Prologs — have some kind of constraint solving embedded. In Ciao we have clpr, clpq and clpfd: packages that extend the basic language.

This has also affected imperative languages, where there are commercial libraries doing the same thing, so you can call a constraint library from C++ or Python and have constraint solving happen while the imperative program runs. The famous ones are ILOG and Gecode. And in functional languages there have been extensions adding constraints, to evaluate expressions containing free variables and give them values — basically doing logic programming inside functional programming.

Slides 3–5 — A comparison with classic LP

Let us learn it by example, with very simple examples.

Herbrand: you have been doing constraint programming all along

:- module(_, _, [clpq]).

% A) Equality constraints on Herbrand terms (unifications)
hp(X,Y,Z) :- Z = f(X,Y).

% B) Equality constraints on arithmetic expressions, classic Prolog
ap(X,Y,Z) :- Z is X + Y.
This clause does one thing: it takes Z and imposes that Z must be f(X,Y).

?- hp(3, 4, Z).
Expected answer: 

Z = f(3,4) ?
A concrete value. Ask it the other way and it also runs:

?- hp(X, Y, f(3,4)).
Expected answer: 

X = 3,
Y = 4 ?
Now — this is constraint programming. Traditional Prolog programming is constraint programming. Why? Because we imposed a constraint, that Z must equal f(X,Y), and then we can ask questions whose answers, first, have no directionality — constraints are symmetrical, they are like logic, they are logic — and second, can be not only values but constraints.

?- hp(X, Y, Z).
Expected answer: 

Z = f(X,Y) ?
That answer is a constraint: whatever we do, Z must be f(X,Y). Ask hp(M,N,f(X,Y)) and the answer is again a constraint — this is fine provided what is in X is the same as what is in M, and Y the same as N.

So this matches the idea of constraint solving perfectly: we impose constraints, and we answer with constraints, possibly more concrete ones.

Prolog arithmetic: where it breaks

Now use is instead of =. Think about how similar they look: Z = f(X,Y) and Z is X+Y differ in that X+Y is evaluated, whereas f(X,Y) is just a tree that gets constructed.

?- ap(3, 4, Z).
Expected answer: 

Z = 7 ?
A concrete value, just as f(3,4) was. But ask anything else:

?- ap(X, 4, 7).
Expected answer: 

{ERROR: arithmetic:is/2, arg 2 - argument is not sufficiently instantiated}
An arithmetic error — because it does not know how to answer with a constraint. It can only answer with concrete values, and here there is no concrete value to give; the answer would have to be a constraint.

So we are just remembering the limitations of arithmetic in Prolog. What we want is something that works like unification, but for arithmetic. And that is exactly what constraints do.

Constraint arithmetic

Why can we write .=. and have it work? Because of :- use_package(clpr)., which loads an extension of Prolog with new operators and constraints — ones that work backwards.

:- module(_, _, [clpr]).

ac(X,Y,Z) :- Z .=. X + Y.
Forwards it behaves the same:

?- ac(3, 4, Z).
Expected answer: 

Z.=.7.0 ?
— with the funny equals, because it is not the Herbrand equals, it is the constraint equals. But now:

?- ac(X, 4, 7).
Expected answer: 

X.=.3.0 ?
It has done the subtraction. And much more interestingly:

?- ac(X, Y, 7).
Expected answer: 

X.=.7.0-Y ?
An equation. We know more than we did before about X and Y: they have to be related. Not all pairs are possible. That is constraint solving in arithmetic, and it is what we said at the start — we answer with a constraint.

The values come back as 7.0 and 3.0 because clpr computes over the reals and everything in it is a float. Load clpq and the same queries answer Z = 7, X = 3, X.=.7-Y — exact rationals, and no equation printed for a variable the solver has pinned down completely. Slides 11–13 come back to this. (The deck shows Z .=. 7 and X .=. 3, which is clpr's form with clpq's numbers.)

Why bother — and what it costs

Advantages. Programs become more expressive and more flexible. We get an arithmetic that is correct from the point of view of logic. It saves a great deal of coding. And in some cases it is much more efficient when we are searching — for two distinct reasons:

  • because all these solvers are implemented inside, off the shelf and very efficiently: just as Prolog gives you the unification algorithm, this extension gives you Gaussian elimination, simplex and the rest;
  • because of search space reduction. In logic programming you must do generate and test; with CLP you can do constrain and generate. Slides 6–10 are that sentence worked out.
Disadvantages. Solvers are complex to implement — not a problem for the user, of course. And constraint solving can affect performance: writing everything in terms of constraints is not necessarily as fast as plain arithmetic. But you gain a great deal of power, and the performance problem is worked on with better solving algorithms, compile-time optimization, global analysis and parallelism. In modern Prolog systems constraints are optional, so you only pay for them if you use them.

Slide 11 — Constraint systems: CLP(𝒳)

If we look at the answers above, ac/3 answered 7.0 — so it is obviously working on floating point numbers, that is, reals in the mathematical sense. Which is why it is called CLP(ℛ): it works over the reals.

So what are the elements of a constraint domain?

  • First, the data: the domain of computation. Reals, integers, booleans, any domain — and these domains have to meet certain conditions.
  • Then the constraints allowed over that domain: for arithmetic, perhaps plus, times, equals, greater than, less or equal.
  • And then constraint solving algorithms: something that, when you write these things that look like equations, gives you the answer — which may itself be an equation. Simplex, or Gaussian elimination, or propagation and consistency.
Which, if you think about it, is exactly what the unification algorithm does for equations over trees — what we do all the time in ordinary Prolog. So what we need here is to generalize that: algorithms that work with these constraints instead of =, and these domains instead of Herbrand trees.

Classical logic programming is a particular case of CLP, where the domain of computation is the Herbrand terms — the trees, the data structures — and there is a single constraint symbol, =.

Remember the examples: X = f(a,3), f(X,Y) = f(a,B). These are equations. And the solution to such an equation — how do I make these two trees equal — is X = a, and making Y be the same variable as B while leaving it free. We have an algorithm that always gives the right answer, the unification algorithm. And that answer is itself a constraint.

So constraints are a complete generalization of logic programming, and logic programming is just this particular case.

The formal statement

CLP is not a language — it is a family of languages, a family of extensions to logic programming. We write CLP(𝒳), and the semantics is parameterized by the constraint domain 𝒳. Fixing 𝒳 gives a particular language: CLP(R), CLP(Q), CLP(FD), CLP(B).

A constraint domain is a tuple 𝒳 = (Σ, D, ℒ, 𝒯):

  • Σ, the signature: the set of predicate and function symbols with their arities. Both your own ps and fs and the constraint symbols themselves — >, <, +.
  • ℒ, a set of Σ-formulae: the constraints you may write, built from those symbols. It has to be a first-order language.
  • D, the set of actual elements — the data you work over: the reals, the integers, the booleans, the finite trees.
  • 𝒟, the Σ-structure giving the meaning of the predicate and function symbols, and hence of the constraints. This is what makes 1 > 2 false.
  • 𝒯, a first-order theory axiomatizing some properties of 𝒟 — think of it as the part that says how to operate with these symbols, in effect the algorithm that solves the constraints.
The pair (D, ℒ) is what we call a constraint domain. So: data; constraints over that data; and algorithms for solving them — a more formal way of saying what slide 2 said informally.

The conditions. There are a number of conditions for all this to work properly, and the full theory — solution compactness, semi-decidability, and the rest of what Jaffar and Lassez imposed so that you get the same results as with traditional first-order logic — is more than we can cover. But the flavor is:

  • ℒ is built on a first-order language. No higher order inside the constraints.
  • = is in Σ and means identity in 𝒟. Every constraint domain has an equality, meaning "the same element of D".
  • There are identically-true and identically-false constraints in ℒ. We need these because every answer has to be a constraint — so we need one that just says "yes" and one that says "no".
  • ℒ is closed under renaming, conjunction and existential quantification, so that we can build up expressions.

Slides 12–14 — Constraint domains

Now the examples. Each is a different choice of Σ and D.

ℛ — arithmetic over the reals. Take the signature 0, 1, +, *, =, <, =<, let D be the reals, and let 𝒟 interpret the symbols as usual. Then you can write things like

x² + 2xy < y/x ∧ x > 0

We do not need symbols for squaring or for the constant 2: x² is x*x, and 2 is 1+1. Keeping Σ minimal is a mathematical simplification to make the presentation canonical — hence the slide's question about whether 0 is needed and how it might be represented.

ℛ_Lin — linear arithmetic. Drop *, so the signature becomes 0, 1, +, =, <, =<. Same data, the reals, but now you may not multiply two variables. Multiplying by a constant is still fine, since 3x is x+x+x. So you can write 3x − y < 3 but not x·y < 3.

ℛ_LinEq — linear equations. Drop the inequalities too, leaving 0, 1, +, =. Now you can write only conjunctions of linear equations, 3x + y = 5 ∧ y = 2x.

Why bother with the restrictions? Because the solving algorithm gets cheaper each time. The general case needs the simplex algorithm, which is relatively complicated. Linear equations alone can be solved by Gaussian elimination, which is much faster. The hierarchy is not mathematical tidiness — it is about what you have to pay at run time.

𝒬 — the rationals. The same set of domains can be defined over the rationals instead: numbers that are a ratio of two integers. And this is not a footnote. Compare:

:- module(_, _, [clpr]).

third(X) :- 3 * X .=. 1.
?- third(X).
Expected answer: 

X.=.0.3333333333333333 ?
against the same program over the rationals:

:- module(_, _, [clpq]).

third(X) :- 3 * X .=. 1.
frac(X, Y) :- X .=. 22/7, Y .=. X * 3.
?- third(X).
Expected answer: 

X = 1/3 ?
?- frac(X, Y).
Expected answer: 

X = 22/7,
Y = 66/7 ?
CLP(Q) gives you 1/3 exactly. This is far more precise than floating point: some numbers are simply not representable as floats, whereas a rational always is — and because Ciao has unbounded integers, numerator and denominator can grow arbitrarily large. In other words, CLP(Q) computes with infinite precision. Looking at the two answers side by side, the rationals are incomparably nicer.

This is a real choice you make in the module header, and it changes your answers. :- module(_, _, [clpr]). gives floats; :- module(_, _, [clpq]). gives exact rationals. 0-Examples/7examples.pl currently loads clpq, so running the deck's examples from that file shows rationals where the older slides show decimals.

There is also a practical trap: because clpq answers are exact rationals, a result that prints as 27/2 is the term /(27,2), and unifying an answer variable with 27/2 will not succeed the way you expect. Compare such answers with .=. rather than =.

Slide 15 — CLP(𝒳) programs

So far we have only described the constraint domain. We have not said what the programs look like.

Let Π ⊆ Σ be the predicate symbols a program is allowed to define. Then:

  • an atom is p(t1, ..., tn) where p is in Π and the ts are terms — exactly what we had in logic programming;
  • a primitive constraint is p(t1, ..., tn) where p is in Σ, a predicate symbol — the same shape, but now p is something like > or <;
  • a constraint is a first-order formula built from primitive constraints — conjunctions, disjunctions.
Which formulae count as constraints varies from domain to domain; generally only a subset is allowed.

A CLP program is a collection of rules a :- b1, ..., bn where a is an atom and each bi is an atom or a constraint. That single change is the whole difference from logic programming: before, a body could only hold calls to predicates; now it can also hold > and <.

A constraint in a body is not a call to a predicate — it is a call to the constraint system. Think of it exactly as you think of a unification: when you write X = f(Y) in a body it is handled by the unification algorithm, not by looking up a predicate. Here, X .>. Y + 1 is handled by the solver. A clause body is a mixture of calls to your predicates and calls to the solver.

A fact is a rule a :- c where c is a constraint. That looks odd at first, so it is worth unpacking. In Herbrand you would write

p(f(X), Y).
which is plainly a fact. But you can equally write it as

p(Z, Y) :- Z = f(X).
All that has happened is that the unifications which were implicit in the head have been peeled out — normalized out — into the body. There is still no call to any other predicate, so it is still a fact: a fact is something that is simply done, it does not go anywhere else. In the same way, p(Z, X, Y) :- Z .>. X + Y. is a fact — it posts constraints and calls nothing.

A goal is a conjunction of constraints and atoms. So you can call p(X) and get values; you can post X .>. 1 and get values; and you can combine the two — give me the ps whose X is greater than 1. And of course the cool thing is that you may put the constraint in front and it works just the same. The Prolog version — the one without the dots — is happy if you put the comparison behind, but not in front: it does not like variables in a comparison, because it does not know what to do with them.

Part II — CLP(ℛ) in practice

Slide 16 — A case study: CLP(ℛ)

We are going to look at two constraint logic programming languages: CLP(ℛ), and CLP over finite domains, which is quite different — and that is exactly why it is worth seeing both.

CLP(ℛ) is a language based on Prolog plus constraint solving over the reals, strictly ℛ_Lin. It keeps Prolog's execution strategy unchanged: depth-first, left to right. It allows linear equations and inequations over the reals. Linear constraints are solved; non-linear ones are passive — delayed until they become linear, or until they become a simple check:

  • X*Y = 7 becomes linear once X is given a definite value;
  • X*X + 2*X + 1 = 0 becomes a check once X is given a definite value.
And the feature that matters most for us: Prolog arithmetic is subsumed by constraint solving.

A little history. At some point CLP(ℛ) was a separate language, meant to replace Prolog. The people who did it concentrated a great deal on the constraint solving — and built an excellent solver — but the Prolog underneath was not very fast, so for ordinary programming it was not so good. Then people put two and two together and realized that what you want is a really good Prolog matched with really good constraint solving, all in the same system. And that is what you get nowadays.

It is really CLP(ℛ, 𝓕𝓣). Because this language was an extension of Prolog, it had both the real constraints and unification — so you could still write your list programs. So in reality it was never CLP(ℛ); it was CLP over two domains at once. Normally nobody mentions that, because you assume you have unification and terms, so you only name the additional domain. In fact in modern systems you can load clpq and clpr, and you can load one constraint domain in one module and a completely different one in another. So it is really CLP of several domains.

Why .=. and not =. Because you have the Prolog parts and the new constraint parts at the same time, and you want to keep the good old ISO Prolog built-ins — otherwise you are not an ISO Prolog — the clpr package uses new primitives, .=. and .>. and so on, to distinguish the CLP(ℛ) constraints from the ISO Prolog arithmetic primitives.

X .=. Y + 5, Y .>. 1        % constraints: solved, and reversible
X is Y + 5,  Y > 1          % ISO arithmetic: evaluated, and one-directional

Exercises

Some constraint programming exercises — opens in the Ciao playground, the same link as the button on the slide.

Slide 17 — Linear equations

The dot product of two vectors, ·: ℛⁿ × ℛⁿ → ℛ, with the vectors represented as lists:

:- module(_, _, [clpr]).

prod([], [], Result) :- Result .=. 0.
prod([X|Xs], [Y|Ys], Result) :-
    Result .=. X * Y + Rest,
    prod(Xs, Ys, Rest).
The dot product of two empty vectors is 0; otherwise multiply X by Y and add it to whatever the dot product of the two remaining lists turns out to be. For fun we represent the vectors as lists — which is worth noticing in itself: the numbers are handled by the constraint solver, the lists by unification, and the two are working together in one clause.

This program would not run in Prolog. Given the two vectors and asked for the result, the second clause posts Result = X*Y + Rest at a point where Rest is not yet known — with is/2 that is an instantiation error. You would have to move the arithmetic after the recursive call, or carry an accumulator. With constraints you can simply state the relation first and let the solver hold it.

?- prod([2, 3], [4, 5], K).
Expected answer: 

K.=.23.0 ?
Backwards, with a hole in one of the vectors:

?- prod([2, 3], [4, Y2], 23).
Expected answer: 

Y2.=.5.0 ?
And now the interesting one. Given a vector, which vectors have dot product zero with it? There is no single answer — there is a plane of them — so what comes back is a constraint:

?- prod([2, 7, 3], [Vx, Vy, Vz], 0).
Expected answer: 

Vz.=. -0.6666666666666666*Vx-2.333333333333333*Vy ?
Any computed answer is, in general, a constraint over the variables in the query.

The slide gives this answer as Vx .=. -1.5*Vz - 3.5*Vy. That is the same plane, solved for a different variable: multiply Ciao's answer by 3 to get 3Vz = -2Vx - 7Vy, rearrange for Vx, and you have the slide's form. Which variable the solver eliminates is its own business.

Load clpq instead of clpr and the same query answers in exact rationals — -2/3 and -7/3 instead of 0.6666666666666666 and 2.333333333333333. As noted on slide 14, the rationals are incomparably nicer to read.

Slide 18 — Systems of linear equations

Can we solve systems? Take

3x + y = 5
x + 8y = 3
Look at what two dot products actually say: [3,1]·[X,Y] = 5 and [1,8]·[X,Y] = 3 are those two equations. So write them both at the prompt and the solver does the rest:

?- prod([3, 1], [X, Y], 5), prod([1, 8], [X, Y], 3).
Expected answer: 

Y.=.0.17391304347826075,
X.=.1.608695652173913 ?
Better still, mimic the mathematical notation A·x = b directly — a matrix of coefficients, a vector of variables, a vector of independent terms:

:- module(_, _, [clpr]).

prod([], [], Result) :- Result .=. 0.
prod([X|Xs], [Y|Ys], Result) :-
    Result .=. X * Y + Rest,
    prod(Xs, Ys, Rest).

system([], _Vars, []).
system([Co|Coefs], Vars, [Ind|Indeps]) :-
    prod(Co, Vars, Ind),
    system(Coefs, Vars, Indeps).
?- system([[3, 1], [1, 8]], [X, Y], [5, 3]).
Expected answer: 

Y.=.0.17391304347826075,
X.=.1.608695652173913 ?
That is, in a dozen lines, a solver for systems of linear equations — and system/3 never mentions solving at all. It only says what the system is.

Exercises

Exercises: a little analytic geometry — opens in the Ciao playground, the same link as the button on the slide.

Slide 19 — Non-linear equations

Non-linear equations are in principle delayed — meaning the system does not solve them, it just remembers them.

:- module(_, _, [clpr]).

trig(X) :- sin(X) .=. cos(X).
poly(X) :- X*X + 2*X + 1 .=. 0.
?- trig(X).
Expected answer: 

_A.=.sin(X),
_A.=.cos(X) ?
Nothing has been solved. The two constraints are simply pending, tied together through an auxiliary variable — it has not solved it, but it does remember. And then if we do something afterwards — say X .=. 1 — it may solve it, or tell us that it does not hold.

The same happens even where a solving procedure does exist:

?- poly(X).
Expected answer: 

_A.=. -1.0-2.0*X,
_A.=.X*X ?
The reason is not that this particular equation is hard — it is that no general technique is known, so CLP(ℛ) solves only linear (dis)equations and defers the rest.

Once enough information arrives to make them linear, they are handled properly:

:- module(_, _, [clpr]).

cossin(X, Y) :- X .=. cos(sin(Y)), Y .=. 2 + Y*3.
?- cossin(X, Y).
Expected answer: 

Y.=. -1.0,
X.=.0.6663667453928805 ?
Y .=. 2 + Y*3 is linear, so it is solved immediately, giving Y = -1; that makes cos(sin(Y)) a constant, and X follows.

And inequations are solved with a modified, incremental simplex:

:- module(_, _, [clpr]).

box(X, Y) :- X + Y .=<. 4, Y .>=. 4, X .>=. 0.
?- box(X, Y).
Expected answer: 

Y.=.4.0,
X.=. -0.0 ?
Three inequations, and the feasible region turns out to be a single point.

Slide 20 — Fibonacci revisited (standard Prolog)

One thing that is fun is to revisit the standard Prolog programs that use arithmetic and recode them with constraints, because then they become reversible.

The Fibonacci numbers: F₀ = 0, F₁ = 1, Fₙ₊₂ = Fₙ₊₁ + Fₙ. The good old Prolog version, using the standard > and is:

:- module(_, _, [clpq]).

fib_pl(0, 0).
fib_pl(1, 1).
fib_pl(N, F) :-
    N > 1,
    N1 is N - 1,
    N2 is N - 2,
    fib_pl(N1, F1),
    fib_pl(N2, F2),
    F is F1 + F2.
It works, and it can only be used with the first argument instantiated to a number:

?- fib_pl(N, 89).
Expected answer: 

{ERROR: arithmetic:>/2, arg 1 - argument is not sufficiently instantiated}
N > 1 cannot be evaluated when N is unbound, so the question "which Fibonacci number is 89?" cannot even be asked.

Slide 21 — Fibonacci revisited (CLP(ℛ))

Syntactically almost the same program:

:- module(_, _, [clpq]).

fib_clp(N, N) :- N .=. 0.
fib_clp(N, N) :- N .=. 1.
fib_clp(N, R) :-
    N .>. 1, F1 .>=. 0, F2 .>=. 0,
    N1 .=. N - 1, N2 .=. N - 2,
    fib_clp(N1, F1), fib_clp(N2, F2),
    R .=. F1 + F2.
Look at what changed. fib_clp(N,N) :- N .=. 0. does not unify N with 0 in the head — it constrains both arguments to be zero. N .>. 1 replaces the >; the two .=. replace the two is; the last line replaces the final is.

And two constraints are here that were not in the Prolog version at all: F1 .>=. 0 and F2 .>=. 0. They are not needed to compute forwards. They are put in because they are true — we are always talking about numbers greater than zero, since the sequence goes up from zero and has no negative members — and because stating what is true bounds the search and is what lets the program run backwards. Stating everything you know is good practice here in a way it is not in Prolog.

Note also that this program uses only numbers and equations: no lists, no other data structures. This is "pure CLP(ℛ)".

And here is the magic of constraints: run it forwards and it works exactly the same, giving 89 for input 11 — and it also runs backwards. Ask the question the Prolog version could not even be asked:

?- fib_clp(N, 89).
Expected answer: 

N = 11 ?
And with both arguments free it simply enumerates the sequence:

?- fib_clp(N, F).
Expected answer: 

F = 0, N = 0 ? ;
F = 1, N = 1 ? ;
F = 1, N = 2 ? ;
F = 2, N = 3 ? ;
F = 3, N = 4 ? ;
0-Examples/7/040_fibonacci.pl calls these fib_pl/2 and fib_clp/2, so that both can live in one file; the slide calls them both fib. The file loads clpq, which is why the answers above are integers rather than the .=. 0.0 the slide shows for clpr.

Exercises

Exercises: a family of Fibonacci sequences — opens in the Ciao playground, the same link as the button on the slide.

Slides 22–23 — Mortgage calculation

The classic CLP(ℛ) example, and the one that shows most directly what reversibility is worth. Five quantities:

  • P — principal, the balance at the beginning
  • T — term, the number of interest periods
  • I — interest rate per period, where 0.1 means 10%
  • B — balance at the end
  • MP — the payment made each period 
:- module(_, _, [clpr]).

mg(P, T, I, B, MP) :-
    T .=. 1,
    B + MP .=. P * (1 + I).
mg(P, T, I, B, MP) :-
    T .>. 1,
    P1 .=. P * (1 + I) - MP,
    T1 .=. T - 1,
    mg(P1, T1, I, B, MP).
Read it as a statement of what a mortgage is: after one period, the balance plus the payment equals the principal grown by the interest; for a longer term, take one period off and recurse on the reduced principal. Nothing in it says which quantity is the input.

So one program answers four different questions.

What will the payments be on a $1000 loan over 30 periods at 3%, paid off completely?

?- mg(1000, 30, 0.03, 0, MP).
Expected answer: 

MP.=.51.0192593202526 ?
What can we afford if we can pay $20 a period?

?- mg(P, 30, 0.03, 0, 20).
Expected answer: 

P.=.392.00882698939535 ?
Seen as an investment instead, what does $1000 grow to over 30 periods if we take nothing out?

?- mg(1000, 30, 0.03, B, 0).
Expected answer: 

B.=.2427.262471189662 ?
And if we add $50 a period — a negative payment:

?- mg(1000, 30, 0.03, B, -50).
Expected answer: 

B.=.4806.033256505762 ?
Four questions a bank has different software for, from one six-line relation. The program was never told which arguments are inputs, because in a relation there is no such thing.

Exercises

Exercises: simulating a population — opens in the Ciao playground, the same link as the button on the slide.

Slides 24–28 — Analog RLC circuits

The last CLP(ℛ) example, and the most striking, because it does not only analyze circuits — it synthesizes them.

Recall that in pure logic programming we already did circuits — with named transistors and gates. Those programs were topological: they described the shape of the circuit, so we could ask "is there a resistor here?" and even run them backwards to find where the ALU sits in a chip. What we could not do was talk about the currents and the voltages, because that needs arithmetic, and making arithmetic reversible does not work in Prolog — it would work with Peano numbers, very inefficiently. Now we have constraints, so we can go back to this kind of problem and work on the currents and voltages too.

The simplest kind of analog circuit with alternating current is the RLC network in steady state: circuits of resistors, capacitors and inductors, where you feed in an alternating current at some frequency and see how it settles. Each circuit is either a simple component or a connection of simpler circuits; we allow only series and parallel connections, which keeps us to Ohm's law. You could do more, but then you have to write Kirchhoff's laws. The entry point is

circuit(C, V, I, W) — across network C, the voltage is V, the current is I, and the frequency is W.

V and I must be complex numbers: the imaginary part is what carries the angular frequency. And note where the two domains meet — the network C is a Herbrand term, built by unification, while the electrical quantities are real constraints.

Complex numbers. X + Yi is the term c(X, Y):

c_add(c(Re1,Im1), c(Re2,Im2), c(Re12,Im12)) :-
    Re12 .=. Re1 + Re2,
    Im12 .=. Im1 + Im2.

c_mult(c(Re1,Im1), c(Re2,Im2), c(Re3,Im3)) :-
    Re3 .=. Re1 * Re2 - Im1 * Im2,
    Im3 .=. Re1 * Im2 + Re2 * Im1.
Adding two complex numbers adds the real parts and adds the imaginary parts; multiplying is the usual rule. And equality is very simple, because unification does it: if two complex numbers are equal, both are c(_,_) and the two fields have to be equal — which is exactly what c_equal(c(R,I), c(R,I)) says, so it needs no definition at all.

Connections. What does it mean to put two circuits in series? The current through the composed circuit is the same as the current through each part — if things are in series the current has to go through everything. And the voltage across the whole is the drop across the first plus the drop across the second, which is a complex addition.

In parallel it is the other way round: the voltage across each part is the same, but the currents split, and the two currents sum to the current through the whole. The two clauses say exactly that and nothing else:

circuit(series(N1, N2), V, I, W) :-
    c_add(V1, V2, V),
    circuit(N1, V1, I, W),
    circuit(N2, V2, I, W).

circuit(parallel(N1, N2), V, I, W) :-
    c_add(I1, I2, I),
    circuit(N1, V, I1, W),
    circuit(N2, V, I2, W).
Components. Each is one line of physics — and the three are more different than they look.

The voltage across a resistor is the current times c(R,0): a resistor has only a real part. This is the classical Ohm's law V = IR. The resistor is transparent to the alternating current — it could not care less whether the current is DC or AC — so the imaginary part is zero and the frequency plays no role at all.

The inductor and the capacitor are different: they behave like resistors whose value varies with the frequency. For the inductor the imaginary part is W*L, the frequency times the inductance; for the capacitor it is -1/(W*C). So the voltage drop across an inductor gets larger with frequency, and across a capacitor smaller.

circuit(resistor(R), V, I, _W) :-        % V = I * (R + 0i)
    c_mult(I, c(R, 0), V).

circuit(inductor(L), V, I, W) :-         % V = I * (0 + WL i)
    Im .=. W * L,
    c_mult(I, c(0, Im), V).

circuit(capacitor(C), V, I, W) :-        % V = I * (0 - 1/(WC) i)
    Im .=. -1 / (W * C),
    c_mult(I, c(0, Im), V).
That is the entire program. It is a description of what a circuit is. Here it is in one piece, in the clause order of 0-Examples/7/060_rlc_circuits.pl — the order matters, because it is the order in which the synthesis below enumerates topologies:

:- module(_, _, [clpr]).

c_add(c(Re1,Im1), c(Re2,Im2), c(Re12,Im12)) :-
    Re12 .=. Re1 + Re2,
    Im12 .=. Im1 + Im2.

c_mult(c(Re1,Im1), c(Re2,Im2), c(Re3,Im3)) :-
    Re3 .=. Re1 * Re2 - Im1 * Im2,
    Im3 .=. Re1 * Im2 + Re2 * Im1.

circuit(resistor(R), V, I, _W) :-
    c_mult(I, c(R, 0), V).
circuit(inductor(L), V, I, W) :-
    Im .=. W * L,
    c_mult(I, c(0, Im), V).
circuit(capacitor(C), V, I, W) :-
    Im .=. -1 / (W * C),
    c_mult(I, c(0, Im), V).
circuit(parallel(N1, N2), V, I, W) :-
    c_add(I1, I2, I),
    circuit(N1, V, I1, W),
    circuit(N2, V, I2, W).
circuit(series(N1, N2), V, I, W) :-
    c_add(V1, V2, V),
    circuit(N1, V1, I, W),
    circuit(N2, V2, I, W).
Analysis. Now take a circuit, describe it, and get results. This one is a parallel circuit: one branch is an inductor of value 0.073, the other branch is a series circuit of a resistor R and a capacitor C, both unknown. Across the whole thing we have 4.5 volts, 0.65 amps, at 2400 Hz. So the question is: what values must R and C have for this to work? Which is typically what you do when designing a filter.

?- circuit(parallel(inductor(0.073),
                    series(capacitor(C), resistor(R))),
           c(4.5, 0), c(0.65, 0), 2400).
Expected answer: 

R.=.6.912283687298522,
C.=.0.0015254646541517786 ?
Synthesis. Now leave the topology itself unbound and ask which circuits behave this way:

?- circuit(C, c(4.5, 0), c(0.65, 0), 2400).
Expected answer: 

C = resistor(_A),
_A.=.6.9230769230769225 ? ;

C = parallel(resistor(_F),resistor(_A)),
_B.=._A*_C,
_B.=.0.0,
_D.=._E*_A,
_D.=.4.5,
...
The first answer is a single resistor of 6.92 Ω. The second is two resistors in parallel, and the answer is not a pair of numbers but a system of constraints relating them — every pair of resistances satisfying it is a valid circuit. Ask for more and it goes on: a resistor in parallel with an inductor, and so on.

The program enumerates circuit topologies, and for each one hands back the equations its components must satisfy. Nothing in those seven clauses was written with synthesis in mind.

Part III — Puzzles, search, and finite domains

Slide 29 — The N queens problem

Place N chess queens on an N × N board so that none attacks another.

Representing the state space. The first thing to do — as with all problem representations — is to represent the state space: to represent the problem with a data structure. And we start from an observation which is really a heuristic, or a symmetry, or a property of the problem space, whatever you want to call it, that lets us prune a great many possibilities.

Queens attack along the diagonals, but they also attack along the rows and the columns. So we cannot have more than one queen in each row, and we cannot have more than one in each column. That prunes the search space quite a bit already — before we have written a line.

So the only thing we need in order to describe a position is, for each row, the column its queen is in: a list of N numbers. And every solution is therefore a permutation of [1,2,...,N].

The board

. Q . .
. . . Q
Q . . .
. . Q .
is the list [2,4,1,3] — first row column 2, second row column 4, third row column 1, fourth row column 3. Checking it by eye: no two queens share a row, a column, or a diagonal.

The general idea. Start from a partial solution — imagine we are at some intermediate step with some queens already placed. Look at where a new queen could go, choose one of those places non-deterministically, and then check whether that queen attacks the ones already placed. If it does not attack, go forward and try to place another. If it does, backtrack and try somewhere else.

So at each point we try to place just one more queen, and at each point we check whether it attacks the others.

The lecture slows down here deliberately — "I'm going pretty slowly because it is not a simple program. It is a complicated problem, and we have to go slowly so that everybody understands and nobody gets lost." It is worth reading the next two sections at that pace.

Slide 30 — The N queens problem in Prolog

:- module(_,_).

queens_pl(N, Qs) :-
    queens_list_pl(N, Ns),          % e.g. Ns = [4,3,2,1]
    queens_pl_(Ns, [], Qs).

queens_pl_([], Qs, Qs).             % all queens placed
queens_pl_(Unplaced, Placed, Qs) :-
    selectq(Unplaced, Q, NewUnplaced),   % e.g. Q=4, NewUnplaced=[3,2,1]
    no_attack_pl(Placed, Q, 1),          % fail if it attacks
    queens_pl_(NewUnplaced, [Q|Placed], Qs).

no_attack_pl([], _Queen, _Nb).
no_attack_pl([Y|Ys], Queen, Nb) :-
    Queen =\= Y + Nb,  Queen =\= Y - Nb,  Nb1 is Nb + 1,
    no_attack_pl(Ys, Queen, Nb1).

selectq([X|Ys], X, Ys).
selectq([Y|Ys], X, [Y|Zs]) :- selectq(Ys, X, Zs).

queens_list_pl(0, []).
queens_list_pl(N, [N|Ns]) :- N > 0, N1 is N - 1, queens_list_pl(N1, Ns).
Four pieces, none of them clever:

queens_list_pl/2 counts down from N and builds [N, ..., 2, 1] — the columns still available. It is plain, deterministic Prolog.

selectq/3 is member/2 with a third argument: it picks one element non-deterministically and returns the rest. On backtracking it yields 4, then 3, then 2, then 1.

no_attack_pl/3 is the only part that needs thought — and it is short because the representation has already done most of the work. Rows and columns need no checking at all: one queen per row is built into the list, one per column into the permutation. Only the diagonals remain. The third argument starts at 1 and grows by one per recursive step, so at depth k we are comparing against a queen k rows away and the two forbidden columns are Y + k and Y - k.

queens_pl_/3 ties it together, carrying the placed queens as an accumulating parameter and pushing the finished list out through the third argument when the unplaced list runs out.

?- queens_pl(4, Qs).
Expected answer: 

Qs = [2,4,1,3] ? ;
Qs = [3,1,4,2] ? ;
no
Exactly two solutions on the 4 × 4 board.

This is about as simple as the problem gets, and that matters for what comes next. You could do very sophisticated things, but look at the program: you really cannot make it much simpler, because it is almost the definition. no_attack_pl/3 is just the definition of what it is to attack; selectq/3 is just the definition of what the possibilities are. It is a few logical rules that give the formulation of the problem — and they serve to solve it.

So when the timings below turn bad, the fault is not that the program was written carelessly.

Slides 32–33 — The N queens problem in CLP(ℛ)

Now the same problem written to constrain first and generate afterwards. no_attack/3 moves in front of the placing, and because it now runs on a list of fresh variables it is folded into the list construction — so constrain_values/3 builds the list and posts the constraints, and place_queens/2 does the generating at the end.

:- module(_,_).
:- use_package(clpr).

queens(N, Qs) :-
    constrain_values(N, N, Qs),
    place_queens(N, Qs).

constrain_values(0, _N, []).
constrain_values(I, N, [X|Xs]) :-
    I .>. 0, X .>. 0, X .=<. N, I1 .=. I - 1,
    constrain_values(I1, N, Xs),
    no_attack(Xs, X, 1).

no_attack([], _Queen, _Nb).            % identical to the Prolog version,
no_attack([Y|Ys], Queen, Nb) :-        % but with constraints
    Queen .<>. Y + Nb,  Queen .<>. Y - Nb,  Nb1 .=. Nb + 1,
    no_attack(Ys, Queen, Nb1).

place_queens(0, _).
place_queens(N, Q) :-
    N .>. 0,
    memberq(N, Q),
    N1 .=. N - 1,
    place_queens(N1, Q).

memberq(X, [X|_]).
memberq(X, [_|Xs]) :- memberq(X, Xs).
no_attack/3 is character for character the same relation as before, with =\= replaced by .<>. and is by .=..

?- queens(4, Qs).
Expected answer: 

Qs = [_A,_B,_C,_D],
_A.=.2.0, _B.=.4.0, _C.=.1.0, _D.=.3.0 ? ;

Qs = [_A,_B,_C,_D],
_A.=.3.0, _B.=.1.0, _C.=.4.0, _D.=.2.0 ? ;

no
The same two solutions.

And now the payoff of writing it as a relation. Leave the board size unbound too:

?- queens(N, L).
Expected answer: 

L = [], N = 0 ? ;
L = [_A], _A.=.1.0, N.=.1.0 ? ;
L = [_A,_B,_C,_D], _A.=.2.0, _B.=.4.0, _C.=.1.0, _D.=.3.0, N.=.4.0 ? ;
L = [_A,_B,_C,_D], _A.=.3.0, _B.=.1.0, _C.=.4.0, N.=.4.0, _D.=.2.0 ? ;
...
It enumerates boards, skipping N = 2 and N = 3 because they have no solutions. Nothing in the program was written with that query in mind.

Note also what is happening underneath: the board is a Herbrand term (a list, built by unification) and the positions are real constraints. As on slide 16, we are using ℛ and 𝓕𝓣 together — and the length of the partially built list Xs in no_attack(Xs, X, 1) is itself doing work.

Slides 34–36 — What the constraints actually look like

It is worth stopping to look at what constrain_values/3 has built before any queen is placed. Run it on its own for a 4 × 4 board:

?- constrain_values(4, 4, Qs).
Expected answer: 

Qs = [_A,_B,_C,_D],
nonzero(_E), nonzero(_F), nonzero(_G), nonzero(_H),
nonzero(_I), nonzero(_J), nonzero(_K), nonzero(_L),
nonzero(_M), nonzero(_N), nonzero(_O), nonzero(_P),
_D.=<.4.0, _D.>.0,
_C.=<.4.0, _C.>.0,
_B.=<.4.0, _B.>.0,
_A.=<.4.0, _A.>.0,
_E.=.3.0+_A-_D,   _F.=. -3.0+_A-_D,
_G.=.2.0+_A-_C,   _H.=. -2.0+_A-_C,
_I.=.1+_A-_B,     _J.=. -1+_A-_B,
_K.=.2.0+_B-_D,   _L.=. -2.0+_B-_D,
_M.=.1+_B-_C,     _N.=. -1+_B-_C,
_O.=.1+_C-_D,     _P.=. -1+_C-_D ?
That is the whole 4 × 4 board expressed as arithmetic: the answer is a list of four variables together with every diagonal relation that must hold between them. Each nonzero(_X) paired with an equation says "these two queens must not differ by exactly this amount".

This is worth doing deliberately with constraint programs: run them with variables, the way you run ordinary Prolog programs backwards, and look at what comes out. With constraints you see the equations coming out, and here they are the whole 4 × 4 board expressed as arithmetic — every diagonal, for every placement. The equations are the board.

And now the point. With all that in the store, putting a number into any position makes every relevant equation fire at once — and the answer comes back immediately, without doing a lot of search. That is the trick.

It is worth pausing on how little code this took. Imagine having to program it.

?- constrain_values(4, 4, Qs), Qs = [3,1|OQs].
Expected answer: 

OQs = [_A,_B],
Qs = [_C,_D,_A,_B],
nonzero(_E), ..., nonzero(_N),
_B.=<.4.0, _B.>.0, _A.=<.4.0, _A.>.0,
_E.=.6.0-_B,  _C.=.3.0,  _D.=.1.0,  _F.=. -_B,
_G.=.5.0-_A,  _H.=.1.0-_A,  _I.=.3.0-_B,  _J.=. -1.0-_B,
_K.=.2.0-_A,  _L.=. -_A,  _M.=.1+_A-_B,  _N.=. -1+_A-_B ?
Notice what changed. _C and _D are now definite values, not equations; twelve nonzeros have become ten; and the surviving equations are simpler. The placement is still feasible — this is the beginning of a possible solution.

Whereas:

?- constrain_values(4, 4, Qs), Qs = [3,2|OQs].
Expected answer: 

no
Rejected outright. Not searched and then rejected — the solver simply found the system inconsistent. Continuing from the feasible one, [3,1,4|_] works and leaves A with only one possible value, and [3,1,4,2] closes it out with every variable determined.

Compare with the tree on slide 31. The constraints sit at the root, so a bad choice is contradicted the moment it is made, at depth one, instead of after a full descent. Most of the tree is never built.

Does it help? Measured in Ciao 1.25, honestly: not much, on this problem.

all solutions            N=8      N=9      N=10     N=11
queens_pl  (Prolog)      0.0006   0.0026   0.012    0.063 s
queens     (CLP(R))      0.030    0.146    0.757    4.17  s
The constrain-and-generate version explores a smaller tree but pays simplex prices at every node, and here the second effect wins. That is not a failure of the technique — it is the wrong solver. CLP(ℛ) is built for linear arithmetic over the reals; this is a problem about placing a finite number of things in a finite number of slots. Which is exactly what the next slides are for.

Slide 37 — Finite domains (I)

Finite domains are a particular type of constraint domain. Remember that we talked about many: arithmetic over the reals, arithmetic over the rationals, strings, booleans. A finite domain constraint solver associates each variable with a finite subset of ℤ — a finite set of integers, positive or negative. So one might write

E ∈ {−123, −10..4, 10}

which in ECLiPSe notation is E :: [-123, -10..4, 10] and in Ciao is

E in -123 \/ (-10..4) \/ 10
where the \/ is a union.

With such variables you can perform arithmetic — +, -, *, / — and establish linear relationships between arithmetic expressions, written #=, #<, #=< and so on. The distinct symbols matter for the same reason .=. did: you may have real or rational constraints in the same program, and the two solvers must not be confused with each other.

We use the distinct symbols so that in the same program you can have some constraints over the reals or rationals and some over finite domains without the two getting confused.

The purpose of all these operations is to narrow the domains of the variables. Post E #> 0 on the variable above and −123 disappears while −10..4 shrinks to 0..4, leaving 0..4 ∪ {10}. That is what it means to say the operations are intended to narrow.

In Ciao this is loaded with :- use_package(clpfd).

Slide 38 — Finite domains (II)

:- module(_,_).
:- use_package(clpfd).

e1(X, A, B) :- X #= A + B, A in 1..3, B in 3..7.
e2(X, A, B) :- X #= A - B, A in 1..3, B in 3..7.
e3(X, A, B) :- X #= A - B, A in 1..3, B in 3..7, X #>= 0.
Addition. 

?- e1(X, A, B).
Expected answer: 

X in 4..10,
A in 1..3,
B in 3..7 ?
The respective minimums are added and the maximums are added. There is no unique solution, and in finite domains this happens all the time: you get back these sets of values.

But note carefully what the answer does not say. It does not mean that every value in the first range works with every value in the second and every value in the third. Some combinations may be incompatible. It says only that the ranges have been narrowed — it has not really solved the problem, because there are not enough constraints. Getting from here to actual values needs enumeration, which is the next slide.

Subtraction. 

?- e2(X, A, B).
Expected answer: 

X in-6..0,
A in 1..3,
B in 3..7 ?
The minimum of X is the minimum of A minus the maximum of B.

Enough constraints for a unique answer. 

?- e3(X, A, B).
Expected answer: 

A = 3,
B = 3,
X = 0 ?
Adding X #>= 0 leaves exactly one point, and propagation alone finds it.

Slide 39 — Finite domains (III): domain, labeling and minimize

domain(Variables, Min, Max) is shorthand for a series of in constraints, so that a list of variables can be given a range in one go.

labeling(Options, VarList). When you get an answer like the ones above, it does not mean every value there is usable; without labeling we do not know what the values actually are. What labeling does is try them: it says, I know X is between 1 and 3, so let me instantiate X to 1 and see what is then valid for the others. It is a backtracking search over the current domains of a set of variables, and it is the step that turns narrowed ranges into answers. Options controls the search order.

:- module(_,_).
:- use_package(clpfd).

pyth(N, X, Y, Z) :-
    domain([X,Y,Z], 1, N),
    X*X + Y*Y #= Z*Z,
    X #>= Y,
    labeling([], [X,Y,Z]).

nolab(X, Y, Z) :-
    domain([X,Y,Z], 1, 1000),
    X*X + Y*Y #= Z*Z,
    X #>= Y.
Pythagorean triples. First without labeling, to see how little propagation achieves on its own:

?- nolab(X, Y, Z).
Expected answer: 

X in 1..1000,
Y in 1..1000,
Z in 1..1000 ?
Nothing at all has been narrowed — a non-linear constraint over three variables gives propagation nothing to bite on. Add the labeling and the answers appear:

?- pyth(1000, X, Y, Z).
Expected answer: 

X = 4,  Y = 3,  Z = 5  ? ;
X = 8,  Y = 6,  Z = 10 ? ;
X = 12, Y = 5,  Z = 13 ? ;
Try pyth(10, ...), then pyth(100, ...), then pyth(1000, ...) and watch the cost: 2 solutions almost instantly, 52 in about a second, and at 1000 it does not finish in several minutes. With nothing for propagation to prune, labeling/2 is enumerating a cube.

minimize(Goal, X) solves Goal while minimizing the variable X — a particular kind of labeling. It is the subject of the project-management slides below.

The lecture describes the labeling options as "depth-first, breadth-first, random", by analogy with Prolog's search rule. That analogy does not survive contact with the real option set. Ciao's labeling/2 (core/library/clpfd/fd_labeling.pl) takes three independent groups:

  • which variable next: leftmost (default), ff, ffc, min, max, random_variable(Seed)
  • how to split it: step (default), enum, bisect
  • which value first: up (default), down, random_value(Seed)
There is no breadth-first option. ff — "first fail", pick the variable with the smallest remaining domain — is a variable selection heuristic, and it is the one that matters on slide 46.

Slide 40 — A classic: SEND + MORE = MONEY

Each letter stands for a digit; find the assignment that makes the sum work.

  S E N D
+ M O R E
---------
M O N E Y
:- module(_,_).
:- use_package(clpfd).

smm([S,E,N,D,M,O,R,Y]) :-
    domain([S,E,N,D,M,O,R,Y], 0, 9),   % all digits 0..9
    0 #< S, 0 #< M,                    % no leading zeros
    all_different([S,E,N,D,M,O,R,Y]),  % all digits different
              S*1000 + E*100 + N*10 + D +
              M*1000 + O*100 + R*10 + E #=
    M*10000 + O*1000 + N*100 + E*10 + Y,
    labeling([], [S,E,N,D,M,O,R,Y]).   % instantiate

smm_no_labeling([S,E,N,D,M,O,R,Y]) :-
    domain([S,E,N,D,M,O,R,Y], 0, 9),
    0 #< S, 0 #< M,
    all_different([S,E,N,D,M,O,R,Y]),
              S*1000 + E*100 + N*10 + D +
              M*1000 + O*100 + R*10 + E #=
    M*10000 + O*1000 + N*100 + E*10 + Y.
Read the program and notice how little of it is doing anything. domain/3 says they are digits; 0 #< S, 0 #< M rules out leading zeros; all_different/1 is shorthand for every disequality between every pair; the long arithmetic constraint is the sum itself, written exactly as it appears on paper. Only the last line — labeling/2 — actually goes looking.

?- smm(V).
Expected answer: 

V = [9,5,6,7,1,0,8,2] ?
That is 9567 + 1085 = 10652.

What propagation alone achieves. Drop the labeling/2 and look at the store:

?- smm_no_labeling(V).
Expected answer: 

V = [_A,_B,_C,_D,1,0,_E,_F],
_A in 8..9,
_B in 2..9,
_C in 2..9,
_D in 2..9,
_E in 2..9,
_F in 2..9 ?
M = 1 and O = 0 fall out with no search at all, S is down to two candidates, and every other letter has lost the digits 0 and 1. Labeling only has to search what is left. This is the whole technique in one query.

The example file's smm/1 calls pp_smm/1 as its last goal, which is what produces the pretty layout:

  9 5 6 7
+ 1 0 8 5
---------
1 0 6 5 2
It is written with format/2, and behaves much as the C function of the same name.

Exercises

Exercise: the same puzzle, now with finite domains — opens in the Ciao playground, the same link as the button on the slide.

Slides 41–42 — A project management problem

A job made of tasks A to G, with dependencies and durations given by this graph:

B, C and D cannot start until A is finished; E and F wait on those; G is the end of the job. Each task takes some time: B takes 1 unit, C takes 2, D takes 3, and so on. The whole thing should finish in 10 time units or less.

Modeling it is almost mechanical. Give an upper bound for the whole project — if in doubt, add up every duration, since the worst case is doing them all in sequence — and let every variable range over 0..10. The value of each variable is the time at which that task starts. Then:

  • A #>= 0 — nothing starts before the beginning; G #=< 10 — nothing ends after the end.
  • B #>= A, C #>= A, D #>= A — the three tasks that wait on A.
  • E #>= B + 1 — E cannot start until B has been running for its full 1 unit. Every arrow in the graph becomes one line of exactly this shape. 
:- module(_,_).
:- use_package(clpfd).

pn1(A,B,C,D,E,F,G) :-
    domain([A,B,C,D,E,F,G], 0, 10),
    A #>= 0, G #=< 10,
    B #>= A, C #>= A, D #>= A,
    E #>= B + 1, E #>= C + 2,
    F #>= C + 2, F #>= D + 3,
    G #>= E + 4, G #>= F + 1.

shortest(A,B,C,D,E,F,G) :-
    minimize(pn1(A,B,C,D,E,F,G), G).

shortest_labelled(A,B,C,D,E,F,G) :-
    minimize(pn1(A,B,C,D,E,F,G), G),
    labeling([], [B,D,F]).
That is the scheduling problem. Note that it never mentions scheduling — it is a transcription of the arrows.

?- pn1(A,B,C,D,E,F,G).
Expected answer: 

A in 0..4,
B in 0..5,
C in 0..4,
D in 0..6,
E in 2..6,
F in 3..9,
G in 6..10 ?
A can start any time from 0 to 4; G can finish any time from 6 to 10. What we are looking at is the slack of every task. (Some further constraints are being respected but are not shown by default.)

But this is not yet an answer: it does not say that any combination of these works. A = 4 together with G = 6 is plainly impossible. To get a concrete schedule you must label.

Minimizing. Labeling gives a schedule; usually what we want is the shortest one — minimize G, the time the project ends:

?- shortest(A,B,C,D,E,F,G).
Expected answer: 

A = 0,
C = 0,
E = 2,
G = 6,
B in 0..1,
D in 0..2,
F in 3..5 ?
The project takes 6 units, and the answer tells you much more than that. A, C, E and G have no slack left: they are the critical tasks, the ones where any delay delays the whole project. B, D and F still have room. That is a critical-path analysis, and we never wrote one.

If you want concrete times for the non-critical tasks too, label what is left:

?- shortest_labelled(A,B,C,D,E,F,G).
Expected answer: 

A = 0, B = 0, C = 0, D = 0, E = 2, F = 3, G = 6 ?
One thing to know before you try this yourself. minimize/2 comes from the clpfd package, and in Ciao packages are loaded per module — and the top level is just another module. So loading clpfd inside 100_project_management_1.pl does not make minimize/2 available at the prompt; importing that module gives you pn1/7, not the solver's control predicates. Load the package at the top level first and everything works:

?- use_package(clpfd).
?- use_module('0-Examples/7/100_project_management_1.pl').
?- minimize(pn1(A,B,C,D,E,F,G), G).

A = 0, C = 0, E = 2, G = 6,
B in 0..1, D in 0..2, F in 3..5 ?
This locality is deliberate: it is what makes the module system compositional, and it is frequently not done by other Prolog systems. The same applies to how answers are printed — the in/2 and .. operators come from the package too, so without it at the top level the same answer prints as in(B,..(0,1)) rather than B in 0..1.

The wrapper predicates above (shortest/7, shortest_labelled/7) are therefore a convenience, not a workaround: they let you run the optimization without loading anything at the prompt, which is what the runnable cells here do.

Slides 43–45 — Two variants

Buying acceleration (slide 43).

Suppose task F can be made faster, at a price. Its duration becomes a variable X, so G #>= F + X. We want the project as short as possible — but we do not want to pay for more speed than we need, so we also want X as large as it can be:

:- module(_,_).
:- use_package(clpfd).

pn2(A,B,C,D,E,F,G,X) :-
    domain([A,B,C,D,E,F,G,X], 0, 10),
    A #>= 0, G #=< 10,
    B #>= A, C #>= A, D #>= A,
    E #>= B + 1, E #>= C + 2,
    F #>= C + 2, F #>= D + 3,
    G #>= E + 4, G #>= F + X.

cheapest(A,B,C,D,E,F,G,X) :-
    minimize(pn2(A,B,C,D,E,F,G,X), G),
    maximize(true, X).
?- cheapest(A,B,C,D,E,F,G,X).
Expected answer: 

A = 0, C = 0, D = 0, E = 2, F = 3, G = 6, X = 3,
B in 0..1 ?
The order of the two optimizations carries the meaning. Done this way the project still takes 6 units and F needs no more acceleration than X = 3. Reverse them — maximize X first — and you get X = 7 with G = 10: the longest possible project with the slowest possible F, which is a correct answer to a question nobody asked.

Both optimizations run from the top level too, once clpfd is loaded there: ?- use_package(clpfd). then ?- minimize(pn2(A,B,C,D,E,F,G,X), G), maximize(true,X). gives exactly the answer above.

A shared resource (slides 44–45).

B and D can each be done in 2 units at best, but they compete for something, so they cannot both be quick: together they take 6. That is three more constraints — X #>= 2, Y #>= 2, X + Y #= 6 — and the durations of B and D become X and Y:

:- module(_,_).
:- use_package(clpfd).

pn3(A,B,C,D,E,F,G,X,Y) :-
    domain([A,B,C,D,E,F,G,X,Y], 0, 10),
    A #>= 0, G #=< 10,
    X #>= 2, Y #>= 2, X + Y #= 6,
    B #>= A, C #>= A, D #>= A,
    E #>= B + X, E #>= C + 2,
    F #>= C + 2, F #>= D + Y,
    G #>= E + 4, G #>= F + 1.

shortest3(A,B,C,D,E,F,G,X,Y) :-
    minimize(pn3(A,B,C,D,E,F,G,X,Y), G).

shortest3_labelled(A,B,C,D,E,F,G,X,Y) :-
    minimize(pn3(A,B,C,D,E,F,G,X,Y), G),
    labeling([], [D,F]).
?- shortest3(A,B,C,D,E,F,G,X,Y).
Expected answer: 

A = 0, B = 0, C = 0, E = 2, G = 6, X = 2, Y = 4,
D in 0..1,
F in 4..5 ?
The solver decides to spend the resource on B (X = 2, the fast option) and let D take 4. Every task except D and F is now critical. And as the slide notes, minimize/2 alone sometimes leaves constraints pending, so a labeling may still be wanted:

?- shortest3_labelled(A,B,C,D,E,F,G,X,Y).
Expected answer: 

A = 0, B = 0, C = 0, D = 0, E = 2, F = 4, G = 6, X = 2, Y = 4 ?

Exercises

Exercises: summary of scheduling with finite domains — opens in the Ciao playground, the same link as the button on the slide.

Slide 46 — N queens with finite domains

Look back at the queens problem: place N numbers, each between 1 and N, subject to disequalities. Small finite sets, a permutation to find — that is a finite-domain problem if anything is.

:- module(_,_).
:- use_package(clpfd).

queens_fd(N, Qs, Type) :-           % Type is the labeling strategy
    constrain_values(N, N, Qs),     % constrain before placing
    all_different(Qs),              % built-in constraint
    labeling(Type, Qs).             % labeling places the queens

constrain_values(0, _N, []).
constrain_values(N, Range, [X|Xs]) :-
    N > 0, N1 is N - 1, X in 1 .. Range,
    constrain_values(N1, Range, Xs),
    no_attack(Xs, X, 1).

no_attack([], _Queen, _Nb).         % same as the CLP(R) version,
no_attack([Y|Ys], Queen, Nb) :-     % with clpfd primitives
    Queen #\= Y + Nb, Queen #\= Y - Nb, Nb1 is Nb + 1,
    no_attack(Ys, Queen, Nb1).
It is essentially the same program again, in finite-domain style: constrain, then generate. Three things changed:

  • The range is stated with in rather than by building a list of candidates — finite domains do that for us. And it is done on the fly: each new variable gets 1..Range, not 1..N, so as queens are placed the later variables are created with smaller ranges already.
  • all_different/1 replaces the hand-written column check.
  • labeling/2 does the generating that place_queens/2 used to do.
no_attack/3 is unchanged except for the operators.

The extra argument. Type is the labeling strategy, there so we can try different ones.

?- queens_fd(20, Q, [ff]).
Expected answer: 

Q = [1,3,5,14,17,4,16,7,12,18,15,19,6,10,20,11,8,2,13,9] ?
Instantaneous. So try a hundred:

?- queens_fd(100, Q, [ff]).
Expected answer: 

Q = [1,3,5,57,59,4,64,7,58,71,81,60,6,91,82,90,8,83,77,65,73,26,9,45,37,63,66,
     62,44,10,48,54,43,69,42,47,18,11,72,68,50,56,61,36,33,17,12,51,100,93,97,
     88,35,84,78,19,13,99,67,76,92,75,87,96,94,85,20,14,95,32,98,55,40,80,49,
     52,46,53,21,15,41,2,27,34,22,70,74,29,25,30,38,86,16,79,24,39,28,23,31,89]
In under a second — a problem plain Prolog could not touch beyond about 22.

Now you see why finite domains are popular. They are a little uglier, but they are amazing: they can do things we cannot do with other techniques.

It is the labeling strategy doing the work, not just the domain. Measured in Ciao 1.25:

N        [ff]        default []
8      0.0019 s      0.0018 s
16     0.0025 s      0.288  s
20     0.0091 s      7.40   s
30     0.016  s      —
50     0.373  s      —
100    0.169  s      —
With the default leftmost selection the same program is already in trouble at N = 20. ff — always branch on the variable with the fewest remaining values — is what makes 100 queens possible.

And for all solutions it does not help at all: 

all solutions       N=8      N=9      N=10     N=11
queens_pl (Prolog)  0.0006   0.0026   0.012    0.063 s
queens_fd (clpfd)   0.033    0.155    0.659    3.33  s
If you have to go through the whole search space, then heuristics and smart things are not going to help you much — the tree is the tree. Constraint propagation buys you the first solution.

And normally that is all anyone wants. Who needs 2,680 different ways to place eleven queens? You only need one.

Exercises

Exercise: the heptagon puzzle — opens in the Ciao playground, the same link as the button on the slide.

Part IV — Other domains, and the wider picture

The lecture skips from the finite-domain queens straight to the survey of systems, for want of time. What follows is the material it passed over: the other constraint domains, what the implementations have to do, and where the ideas came from.

Slide 47 — CLP(𝓕𝓣), otherwise known as logic programming

Plain logic programming is a CLP language: the domain is finite trees, and the solver is the unification algorithm. Nothing needs adding to make the point — here is a program that would be awkward in most languages and is three lines here.

Two trees are isomorphic if they have the same elements at each level, allowing subtrees to be swapped:

:- module(_,_).

iso(Tree, Tree).
iso(t(R, I1, D1), t(R, I2, D2)) :-
    iso(I1, D2),
    iso(D1, I2).
The second clause says: same root, and the left of one matches the right of the other, and vice versa. Now ask it to complete two partial trees so that they become isomorphic:

?- iso(t(a, b, t(X, Y, Z)), t(a, t(u, v, W), L)).
Expected answer: 

L = b, X = u, Y = v, Z = W ? ;
L = b, X = u, Y = W, Z = v ? ;
L = b, W = t(_A,_C,_B), X = u, Y = t(_A,_B,_C), Z = v ? ;
L = b, W = t(_A,t(_C,_E,_D),_B), X = u, Y = t(_A,_B,t(_C,_D,_E)), Z = v ? ;
The answers are constraints on the trees, and there are infinitely many of them — each one a larger pair of mirror-image shapes. That is the same behavior we saw from .=. on slide 17, in the domain we have had all along.

Slides 48–49 — CLP(𝓦𝓔) and CLP(𝓦𝓔, 𝒬)

Word equations. Take the domain to be finite strings, with concatenation (.) and string length (::) as the primitive constraints. Then you can ask for strings with a given property:

?- "123".Z = Z."231", Z::0.      ?- "123".Z = Z."231", Z::3.
no                               no

?- "123".Z = Z."231", Z::1.      ?- "123".Z = Z."231", Z::4.
Z = "1"                          Z = "1231"

?- "123".Z = Z."231", Z::2.
no
Which strings Z satisfy "123"Z = Z"231"? Of length 1, "1"; of length 4, "1231"; of lengths 0, 2 and 3, none. These solvers are very complex, and the algorithms used are often incomplete.

With arithmetic added. CLP(𝓦𝓔, 𝒬) combines word equations with arithmetic over the rationals, and the classic demonstration is a small piece of real mathematics. The sequence xᵢ₊₂ = |xᵢ₊₁| − xᵢ has period 9, whatever x₀ and x₁ you start from. Prove it as follows: describe the sequence; then look for a subsequence that would violate the period condition.

seq(<Y, X>).                          abs(Y, Y) :- Y >= 0.
seq(<Y1 - X, Y, X>.U) :-              abs(Y,-Y) :- Y < 0.
    seq(<Y, X>.U),
    abs(Y, Y1).
(Prolog III syntax, slightly modified.) And then the question — is there any 11-element sequence whose initial 2-tuple differs from the final one? — is

?- seq(U.V.W), U::2, V::7, W::2, U#W.
fail
The failure is the proof. U::2 and W::2 say the first and last chunks are two elements long, V::7 fixes seven in between, and U#W demands that the seed and the trailing pair differ. No such sequence exists, so the period is 9.

Slide 50 — Summarizing

In general. Data structures (Herbrand terms) come for free, and every logical variable can carry constraints relating it to others.

Problem modeling. Rules represent the problem at a high level — they give you program structure and modularity, and recursion is what sets the constraints up. Constraints encode the conditions. And solutions are themselves expressed as constraints.

Combinatorial search. CLP languages give you backtracking, so enumeration is easy to write; the constraints are what keep the search space manageable. Slides 29–46 are that sentence worked out at length.

Tackling a problem. Keep an open mind: new approaches are often possible. The circuit program on slide 28 is the example to remember — nobody wrote it in order to synthesize circuits.

Slide 51 — Complex constraints

Systems also offer complex constraints, which express many simpler constraints compactly. Operationally they are often treated as passive.

The cardinality operator #(L, [c₁, ..., cₙ], U) says that the number of constraints in the list that hold lies between L and U — and L and U may themselves be variables. It subsumes a surprising amount:

  • L = U = n — all of them must hold (conjunction)
  • L = U = 1 — exactly one is true
  • U = 0 — none is true, i.e. the conjunction of the negations
  • L > 0 — at least one is true (disjunction)
A disjunctive constructive constraint c₁ ∨ c₂, if properly handled, avoids search and backtracking altogether. Compare writing "X is non-zero" as two clauses —

nz(X) :- X .>. 0.
nz(X) :- X .<. 0.
— which creates a choice point, with posting it as one disjunctive constraint, which does not.

Slide 52 — Other primitives

Beyond constraints proper, CLP systems provide:

  • enum(X) — enumerate X within its current domain.
  • maximize(X), minimize(X) — find the maximum or minimum of X under the active constraints. We used minimize/2 on slides 41–45.
  • delay Goal until Condition — run Goal only once the variables are instantiated enough for it to do something useful. Using it well needs a deep understanding of the constraint system; it is widely available in Prolog systems too; and strictly it is not a constraint at all but a control primitive.

Slide 53 — Implementation: satisfiability

Deciding satisfiability has to be incremental to be practical — constraints arrive one at a time as the program runs, and re-solving from scratch at each step is hopeless. Good average-case behavior matters more than good worst-case behavior.

The common technique is to keep satisfiable constraints in a solved form. This is not possible in every domain, but 𝓕𝓣 is the familiar example: constraints are held as

x₁ = t₁(ỹ), ..., xₙ = tₙ(ỹ)

where each tᵢ(ỹ) is a term over the variables ỹ, and no xᵢ occurs in ỹ. Those are exactly the idempotent substitutions that the unification algorithm produces — so Prolog has been keeping its constraint store in solved form all along.

Slide 54 — Implementation: backtracking

Backtracking is harder than in Prolog, because changes to constraints must be undone as well as bindings. Constraints are typically stored as an association from a variable to an expression, and trailing those expressions is costly — you cannot afford to record every change.

The saving is that a change only needs trailing if there is a choice point between it and the previous one; otherwise both will be undone together anyway. The standard technique is time stamping: compare the age of the choice point with the age of the variable when it was last trailed, and skip the trailing when no choice point intervenes.

Slides 55–56 — Implementation: extensibility, and attributed variables

Modern Prolog systems implement constraint domains in Prolog itself, using two mechanisms.

Attributed variables provide a hook into unification: attach an attribute to a variable, and when that variable is unified, user code is called. They are how constraint solvers get written — and they have other uses too, such as distributed execution.

Constraint handling rules (CHR) are the higher-level abstraction: they let you define propagation algorithms declaratively, and are usually translated down to attributed-variable code.

The primitives are attach_attribute/2, get_attribute/2, detach_attribute/1, update_attribute/2, verify_attribute/2 and combine_attributes/2. The classic small example is freeze — delay a goal until a variable is bound:

freeze(X, Goal) :-
    attach_attribute(V, frozen(V, Goal)),
    X = V.

verify_attribute(frozen(Var, Goal), Value) :-
    detach_attribute(Var),
    Var = Value,
    call(Goal).

combine_attributes(frozen(V1, G1), frozen(V2, G2)) :-
    detach_attribute(V1),
    detach_attribute(V2),
    V1 = V2,
    attach_attribute(V1, frozen(V1, (G1, G2))).
freeze/2 attaches the goal to a fresh variable and unifies it with X. verify_attribute/2 is what unification calls when the variable finally gets a value: detach, bind, run the goal. combine_attributes/2 handles the case where two frozen variables are unified with each other — the two goals are conjoined onto the merged variable. Thirteen lines, and coroutining exists.

Slide 57 — Programming tips

Over-constrain deliberately. This looks like a violation of "don't do unnecessary work", but stating extra constraints — even redundant ones — can cut much more of the search space than they cost to check. It matters especially when propagation is weak, or when some of your conditions live outside the constraint domain. The F1 .>=. 0, F2 .>=. 0 on slide 21 was exactly this.

Use cut very sparingly, and carefully. Determinacy is subtler than in Prolog, partly because constraints may sit in non-solved form: choosing a clause does not rule out the other clauses the way plain unification does. Compare —

:- module(_,_).
:- use_package(clpr).

max(X, Y, X) :- X .>. Y.
max(X, Y, Y) :- X .=<. Y.

maxcut(X, Y, X) :- X .>. Y, !.
maxcut(X, Y, Y) :- X .=<. Y.
?- max(X, Y, Z).
Expected answer: 

Z = X, Y.<.X ? ;
Z = Y, X.=<.Y ? ;
no
?- maxcut(X, Y, Z).
Expected answer: 

Z = X, Y.<.X ? ;
no
Both answers of max/3 are correct and they are mutually exclusive — but nothing in the machinery knows that, so the choice point stays. The cut removes it, and in doing so throws away a correct answer in every call where the first clause happens to match first. In Prolog you would reach for the cut here without much thought; in CLP that instinct is worth resisting.

Slide 58 — CLP systems

CLP defines a class of languages, obtained by choosing

  • the particular constraint system or systems, and
  • the computation and selection rules.
Most practical systems include the Herbrand domain with = and then add further domains and solver algorithms on top; and most keep Prolog's computation and selection rules.

Slides 59–60 — The original systems, and the current ones

The originals. CLP(ℛ) did linear arithmetic over the reals with incremental Gaussian elimination and incremental simplex. Prolog III had CLP(R), two-valued booleans, infinite (rational) trees, and equations over finite strings — the CLP(𝓦𝓔) of slide 48. CHIP, and its successor the ILOG library, had CLP(FD), CLP(B) and CLP(Q), with user-defined constraints and solver algorithms. BNR-Prolog did arithmetic over closed real intervals, plus FD and booleans. RISC-CLP handled non-linear real arithmetic. clp(FD)/gprolog did finite domains.

Today: most Prolog systems support constraints. SICStus, ECLiPSe, SWI and Ciao all have CLP(R), CLP(Q) and CLP(FD); SWI adds CLP(B). SICStus, SWI and Ciao offer attributed variables and CHR for adding your own domains — ECLiPSe does not.

And one thing unique to Ciao: because it can load different syntax and semantics per module, you can have different constraint domains active in different modules of the same program. Most systems let you work with one domain at a time, or require you to keep them apart.

The lecture adds an aside here: "SICStus, which has CLP(R), CLP(Q) and CLP(FD) — and also CLP(B), if I remember well; I think it's missing here." He is right, and it is recoverable: the deck source has % , CLP(B). commented out on the SICStus line.

Slide 61 — Origins, and other instances

Ancestors. SKETCHPAD (1963), Waltz's algorithm (1965), ThingLab (1981), Macsyma (1983).

Constraints in logic languages — the origin of "constraint programming" as we know it. The general theory was developed by Jaffar and Lassez in 1987; standalone systems came first (clpr, CHIP, ...) and were later folded into mainstream Prolog implementations. It has given rise to an entire research area.

Constraints in imperative languages. Equation-solving libraries such as ILOG and Gecode. And note the trick that makes constraints meaningful there at all: time stamping variables, so that x := x + 1 becomes xᵢ₊₁ := xᵢ + 1 — a relation between two distinct variables rather than a destructive update. It is the same device used by iterative methods in numerical analysis.

Constraints in functional languages, via extensions: evaluating expressions containing free variables, and absolute set abstraction.

The general theory is Jaffar and Lassez, "Constraint Logic Programming", POPL 1987. Both the slide and the lecture originally said '97; the slide has been corrected.

Generated with LPdoc using Ciao