« A first look at G12 Zinc: Basic learning models etc | Main | Improvements of some Google CP Solver models »

A first look at Google CP Solver/Python (Google or-tools)

Here is my Google CP Solver page.

Two weeks ago there was a new Google project announced: Google or-tools (Operations Research Tools developed at Google) which consists of a constraint programming solver. It also include support for local search (Tabu search, simulated annealing, etc) but I have not studied these features yet.

The last week I started to implement my standard "learning problems" in CP Solver/Python. These models (with some omissions and additions) are listed below. Note that I have only used at the Python interface and not at all studied the C++ interface.

Here is the announcement on the project group: Welcome, presenting the tool:
The code we have released is a subset of all the OR tools we are developing at Google. This code is used internally. We will maintain this open-source branch in sync with our internal code and will likely contribute more tools, more technology.

The documentation is missing. We are working hard on that. In the meantime, do not hesitate to contact us, and please have a look at our C++ and python examples.

We are currently very limited in our ports. In particular, we are not supporting xcode on mac OS X Leopard. The STL is not standard there. We do not plan to support xcode < 3.2.3 on Snow Leopard. Sorry for the trouble. We are trying to find a work around using custom GCC on mac OS X.
The solver is available via SVN: http://or-tools.googlecode.com/svn/trunk/, and is obtained by the SVN command:
svn checkout http://or-tools.googlecode.com/svn/trunk/ or-tools-read-only
Note: it requires google-gflags and zlib.

Documentation: or-tools Wiki

Getting Started

Forum: Google group or-tools-discuss

Examples: There is not yet any tutorial or complete reference of the Python interface to the constraints. The Wiki page UsingCPSolverCppAPI lists the C++ version of the methods with some comments. The Python methods are (almost) always the C++ version minus the leading Make. E.g. the Python version of the C++ method MakeAllDifferent is AllDifferent (in the models: solver.AllDifferent).

There are classes of constraints that I have not (yet) used in the models below, especially the support for Scheduling, and Local Search.

See also:
Jean-Charles Régin and Pierre Schaus has already blogged about the CP Solver on their blog Constraint Programming: I guess (and hope) that Pierre will write more about the Google solver and give nice examples how to use the Local Search support.

Some general comments on the Python interface

Here are some general comments about the Python interface. Please note that I have only used the system for a week and may have missed a thing or two. Since I'm not an expert Python programmer there may also be some shortcuts missed in the models.

I like the Google CP Solver/Python. It is easy to use the solver, especially if one has used some other CP system before and is comfortable in Python. Here is a simple example of a standard CP problem: magic square (the complete model is magic_square.py):
from constraint_solver import pywrapcp
def main(n=3):
    # Create the solver.
    solver = pywrapcp.Solver('Magic square')
    # declare variables
    x = {}
    for i in range(n):
        for j in range(n):
            x[(i, j)] = solver.IntVar(1, n*n, 'x(%i,%i)' % (i, j))
    x_flat = [x[(i,j)] for i in range(n) for j in range(n)]

    # the sum
    s = solver.IntVar(1, n*n*n,'s')
    # constraints  
    solver.Add(solver.AllDifferent(x_flat, True))
    [solver.Add(solver.Sum([x[(i,j)] for j in range(n)]) == s) for i in range(n)]    
    [solver.Add(solver.Sum([x[(i,j)] for i in range(n)]) == s) for j in range(n)]        
    solver.Add(solver.Sum([ x[(i,i)]     for i in range(n)]) == s) # diag 1
    solver.Add(solver.Sum([ x[(i,n-i-1)] for i in range(n)]) == s) # diag 2       

    # solution and search
    solution = solver.Assignment()
    db = solver.Phase(x_flat,

    # output
    num_solutions = 0
    while solver.NextSolution():
        print "s:", s.Value()
        for i in range(n):
            for j in range(n):
                print "%2i" % x[(i,j)].Value(),
        num_solutions += 1        

    print "num_solutions:", num_solutions
    print "failures:", solver.failures()
    print "branches:", solver.branches()
    print "wall_time:", solver.wall_time()
The output for n=3 is
s: 15
 6  7  2
 1  5  9
 8  3  4


num_solutions: 8
failures: 448
branches: 910
wall_time: 0
I haven't stumbled upon any bugs at all in the system, which is really great. Maybe it's not that surprising, given that Google has used this tools internally for some times, and that Google is ... Google.

A few of my models are not very fast (see below), but I think this is because I either screwed up something or didn't found any good search heuristics.

And here are my (mandatory) complaints about syntactic sugar and support for different things:
  • It would be nice to have built-in support for IntVar matrices, e.g.
    • making an IntVar matrix
    • flattening an IntVar matrix
  • It would also be nice with some syntactic sugar for
    • reifications
    • boolean operations on constraints
    • Element constraints in arrays/matrices
    Note: there are support for reifications and Elements; I just would like to have then "sugared".
  • Unfortunately, Ctrl-C don't propagates to the solver so I have to Ctrl-Z and kill the process (kill %%). Well, this is really not a problem with the CP Solver, but how Python interacts with C++.
In my description of my models below, I'll write some more about these things. For a general discussion (rant) about the support of Element constraints in different CP systems, see Learning constraint programming - part II: Modeling with the Element constraint; please note that some systems, especially Gecode, has much better support for the Element constraints in later versions.

My Google CP Solver Models

Here are my first models in Google CP Solver so far. As mentioned above they are most of my "learning problems", which some omissions (see comments below). The order of the models is about the same as I implemented them so later models may have features that earlier lacks. Also, I include some comments of the problem and/or the objective of the "lesson", i.e. what I test with the problem.

The files includes problem description and references to other implementations. (See also the page Common CP models for a listing of all my models that are implemented in at least two different CP systems.)

Except for a few problems, I have not tried to tweaked the search strategies; instead I settled for a "good enough" principle.

Least Diff

Model: least_diff.py
  • minimize the difference ABCDE - FGHIJ (digits all distinct)
  • this was actually my first "real" (i.e. non-standard) constraint programming problem
  • how to optimize (here minimize) a solution
This model is based on the sendmore.py (SEND+MORE=MONEY) model in the distribution. Added is the objective to minimize the difference:
solver.Add(x == 10000*a +1000*b +100*c +10*d + e)
solver.Add(y == 10000*f +1000*g +100*h +10*i + j)
solver.Add(diff == x - y)
solver.Add(diff > 0)  
solver.Add(solver.AllDifferent(letters, True))

# objective
objective = solver.Minimize(diff, 1)

# ...
solver.Solve(solver.Phase(letters + [x,y,diff],
                            [objective, search_log, collector])


Model: diet1.py
Model: diet1_b.py
  • simple OR problem
  • how to interact with integer arrays and variable arrays
This is a simple diet problem. Here is the complete model:
def main(): 
  n = 4
  price  = [ 50, 20, 30, 80] # in cents
  limits = [500,  6, 10,  8] # requirements for each nutrition type

  # nutritions for each product
  calories  = [400, 200, 150, 500]
  chocolate = [3,2,0,0]
  sugar     = [2,2,4,4]
  fat       = [2,4,1,5]

  # declare variables
  x = [solver.IntVar(0, 100, 'x%d' % i) for i in range(n)]
  cost = solver.IntVar(0,10000, 'cost')

  # constraints
  solver.Add(solver.Sum([x[i]*calories[i]  for i in range(n)])   >= limits[0])
  solver.Add(solver.Sum([x[i]*chocolate[i] for i in range(n)])   >= limits[1])
  solver.Add(solver.Sum([x[i]*sugar[i]     for i in range(n)])   >= limits[2])
  solver.Add(solver.Sum([x[i]*fat[i]       for i in range(n)])   >= limits[3])
  # objective
  objective = solver.Minimize(cost, 1)

  # solution
  solution = solver.Assignment()

  # last solution since it's a minimization problem
  collector = solver.LastSolutionCollector(solution)
  search_log = solver.SearchLog(100, cost)
  solver.Solve(solver.Phase(x + [cost],
                            [objective, search_log, collector])

  # get the first (and only) solution
  current = collector.solution(0)

  print "cost:", current.Value(cost)
  print [("abcdefghij"[i],current.Value(x[i])) for i in range(n)]
The constraints could be stated more simple using ScalProd (scalar product) instead of Sum (and is implemented like this in diet1_b.py):
solver.Add(solver.ScalProd(x,calories)  >= limits[0])
solver.Add(solver.ScalProd(x,chocolate) >= limits[1])
solver.Add(solver.ScalProd(x,sugar)     >= limits[2])
solver.Add(solver.ScalProd(x,fat)       >= limits[3])
An interesting feature is the use of solver.SearchLog(100, cost) (called a search monitor) which prints what happens during search. For this small problem instance it may not be very useful, but for bigger problems it may be very convenient.
[10:56:12] constraint_solver/search.cc:204: Start search, memory used = 0
[10:56:12] constraint_solver/search.cc:204: Root node processed (time = 0 ms, constraints = 4, memory = memory used = 0)
[10:56:12] constraint_solver/search.cc:113: Solution #0 (objective value = 0, time = 0 ms, branches = 5, failures = 0, depth = 5, memory used = 0)
[10:56:12] constraint_solver/search.cc:204: Finished search tree, time = 0 ms, branches = 5, failures = 6, memory used = 0)
[10:56:12] constraint_solver/search.cc:204: End search (time = 0 ms, branches = 5, failures = 6, memory used = 0)
cost: 0
[('a', 0), ('b', 3), ('c', 0), ('d', 1)]

failures: 6
branches: 5
wall_time: 0

Seseman Convent problem

Model: seseman.py
Model: seseman_b.py
  • this is one of my earliest constraint programming (modeled in ECLiPSe)
  • generating one or all solutions
  • handling of matrices
For the first version (seseman.py) it was a problem with even for a quite small problems (n=4): When looking for all solutions (using collector = solver.AllSolutionCollector(solution)) the RAM was eaten up since all solution is first collected before the output. Here is the search part:
 collector = solver.AllSolutionCollector(solution)
 solver.Solve(solver.Phase([x[(i,j)] for i in range(n) for j in range(n)],
 num_solutions = collector.solution_count()
 for s in range(num_solutions):
     current = collector.solution(s)
     print "total_sum:", current.Value(total_sum)
     for i in range(n):
         for j in range(n):
             print current.Value(x[(i,j)]), 

However, in Getting all silutions with google solver, Pierre Schaus showed a better way of showing all solutions. Here is the alternative way (which nowadays is my prefered approach) which is implemented in seseman.py:
db = solver.Phase([x[(i,j)] for i in range(n) for j in range(n)],

num_solutions = 0

while solver.NextSolution():
    num_solutions += 1
    print "total_sum:", total_sum.Value()
    for i in range(n):
        for j in range(n):
            print x[(i,j)].Value(),
Matrices is not directly supported in Google CP Solver. Instead one had to build these manually in some way. Here is the tuple variant:
x = {}
  for i in range(n):
     for j in range(n):
         x[(i,j)] = solver.IntVar(0,n*n, 'x %i %i' % (i, j))
However, when using this "matrix" different context, it need to be flattened.
# ...
solution = solver.Assignment()
solution.Add([x[(i,j)] for i in range(n) for j in range(n)])

solver.Solve(solver.Phase([x[(i,j)] for i in range(n) for j in range(n)],
In later models I tend to define a flattened version directly and use that instead:
x_flat [x[(i,j)] for i in range(n) for j in range(n)]

# ...
solution = solver.Assignment()

# etc ...


Model: coins_grid.py
  • Tony Hubermann's grid puzzle, minimize distances in the grid
  • how does a CP system handle this MIP problem (mostly quite bad compared to MIP)
No surprises here. As with other CP solver, this is not very fast in Google CP Solver either.


Model: toNum.py

This model is just a test case for the method toNum which converts an integer to/from an array of digits (where the domain is 0..base-1):
def toNum(solver, t, s, base):
    tlen = len(t)
    solver.Add(s == solver.Sum([(base**(tlen-i-1))*t[i] for i in range(tlen)]))
This is used in de Bruijn sequence model below.

de bruijn sequence

Model: debruijn_binary.py
  • a personal favorite
  • see how large problem an implementation can handle
  • command line options to the model
This model creates de Bruijn sequences.

Since this is a Python program, handling command line options is a cinch, just use the standard Python way.
base = 2
n    = 3
m    = base**n
if __name__ == '__main__':
    if len(sys.argv) > 1:
        base = string.atoi(sys.argv[1])
    if len(sys.argv) > 2:
        n = string.atoi(sys.argv[2])
    if len(sys.argv) > 3:
        m = string.atoi(sys.argv[3])
    main(base, n, m)

Also, here I use the global constraint Distribute (global cardinality count) for showing and restricting the distribution of the digits.
# ensure that all the numbers in the de Bruijn sequence
# (bin_code) has the same occurrences (if check_same_gcc is True
# and mathematically possible)
gcc = [solver.IntVar(0,m,'gcc%i' %i) for i in range(base)  ]
solver.Add(solver.Distribute(bin_code, range(base), gcc))
if check_same_gcc and m % base == 0:
    for i in range(1,base):
        solver.Add(gcc[i] == gcc[i-1])
Some statistics for n = 4 with different value of base:
Base Seconds Failures

alldifferent_except_0: Decomposition of a constraint

Model: alldifferent_except_0.py
  • how to implement a (decomposition of a) global constraint
  • logical operators on variables: implication (if then), &&, !=
This is a decomposition of the global constraint alldifferent_except_0. See the Global Constraint Catalog alldifferent_except_0 for more about this constraint.

In this implementation, I use reifications and logical operations on variables. However, CP Solver has no syntactic sugar in the Python interface for stating reifications or logical operations on variables. It has (which Laurent Pierron kindly pointed out to me) support for it via C++ methods listed in the Constraints section of the C++ API documentation.

In "sugered" CP systems like MiniZinc, Zinc, Comet, and Gecode the constraint can be stated like this (using MiniZinc/Zinc syntax):
forall(i,j in index_set(x) where i != j) (
        (x[i] > 0 /\ x[j] > 0) -> x[i] != x[j]      
The CP Solver is more verbose using temporary BoolVar for the reification and a standard OR way of stating the condition:
# Decomposition of alldifferent_except_0
def alldifferent_except_0(solver, a):
    n = len(a)
    for i in range(n):
        for j in range(i):
            bi = solver.BoolVar('bi')
            bj = solver.BoolVar('bj')
            bij = solver.BoolVar('nij')
            solver.Add(solver.IsDifferentCstCt(a[i],0, bi))
            solver.Add(solver.IsDifferentCstCt(a[j],0, bj))
            solver.Add(solver.MakeIsDifferentCt(a[i],a[j], bij))
            solver.Add(bi*bj <= bij)


Model: nqueens.py
Model: nqueens2.py

The problem is stated with this AllDifferent constraint (where the second parameter True means that a Bounds version of the constraint should be used; when False a Value based version is used).
together with either of these variants:
for i in range(n):
    for j in range(i):
        solver.Add(q[i] != q[j])
        solver.Add(q[i] + i != q[j] + j)
        solver.Add(q[i] - i != q[j] - j)
for i in range(n):
     for j in range(i):
         solver.Add(abs(q[i]-q[j]) != abs(i-j))
The first version (nqueens.py) use solver.AllSolutionCollector(solution), while the second (nqueens2.py) use while solver.NextSolution() approach. See the Seseman problem above for more discussions about this.

Later add: Sylvain Soliman asked me on Twitter for some benchmark. Using the nqueens2.py model I got the following using the variable strategy CHOOSE_MIN_SIZE_LOWEST_MIN and value strategy ASSIGN_CENTER_VALUE:
  • N=300 almost immediately
  • N=500: some second
  • N=1000 takes about 18 seconds using almost 1Gb
System: I'm on a Linux Ubuntu (10.4), 8 core 64 bit with 12 Gb RAM.

Furniture moving (simple scheduling problem)

Model: N/A

This is a small scheduling problem how to best move some furnitures. The traditional way is to use cumulative. However, the CP Solver has no direct support for cumulative. Instead it has a quite advanced and general support for scheduling.

Being lazy, I have not (yet) done a cumulative wrapper, and skipped this model so far.


Model: minesweeper.py
  • more advanced example (quite simple to implement, though)
  • reading problems from a file or handling different problems in some other way
As always, I'm amazed how easy it is to state this problem in a CP system. The main constraint is:
X = -1
game = [

# ....
for i in range(r):
   for j in range(c):
       if game[i][j] >= 0:
           solver.Add(mines[i,j] == 0)
           # this cell is the sum of all the surrounding cells
               game[i][j] == solver.Sum([mines[i+a,j+b]
                                  for a in S for b in S
                                  if i + a >=0 and
                                     j + b >=0 and
                                     i + a < r  and
                                     j + b < c
       if game[i][j] > X:
           # This cell cannot be a mine
           solver.Add(mines[i,j] == 0)

Problem instances: Reading the files are easy with Python file and string handling.

Quasigroup completion

Model: quasigroup_completion.py
  • yet another grid puzzle/problem
  • alldifferent on rows/columns
  • handling matrices
  • handling of different problems (file or some other way)
Problem instances:

Survo puzzle

Model: survo_puzzle.py
  • handling of different problems (file etc )
As a learning problem, the Survo Puzzle is about the same things as quasigroup_completion, i.e. reading files etc, and may be redundant. I may skip one of these two some day.

Problem instances:

Young tableaux and partitions

Model: young_tableaux.py
  • combinatorial problem, easy to state and sometimes tricky to model
This was a straight forward conversion from other implementations. The thing that got me confused somewhat was to convert the 1-based MiniZinc model (young_tableaux.mzn) to the 0-based Python model (and using range(low, high) where the second (or only) parameter (high) is non-inclusive). Well, this is just confusing when looking at two implementations at the same time...

Another thing is reification. Instead of MiniZinc's succint way of stating reification:
forall(1..n) (
    p[i] == sum(j in 1..n) (bool2int(x[i,j] <= n))
we have to create a temporary array (b) and check the requirement with IsLessOrEqualCstVar:
for i in range(n):
    b = [solver.BoolVar("b%i"%j) for j in range(n)]
    for j in range(n):
        solver.Add(b[j] == solver.IsLessOrEqualCstVar(x[(i,j)],n))
    solver.Add(p[i] == solver.Sum(b))


Model: xkcd.py
  • simple problem (easy to understand)
  • knapsack / set covering
No surprises here.

Who killed Agatha?

Model: N/A
  • logical problem
  • reification
  • Element constraint in matrices etc
Well, this is another model I skip for now since the lack of syntactic sugar for reifications, logical operations on variables and stating more complex ways of Element constraints.

Send Most Money

Model: send_most_money.py
  • in some versions: first optimize, then generating all solutions
Again, being a Python program, it is very easy to use the same method as both minimization problem and an "all-solution" model.

Send More Money in "any" base

Model: send_more_money_any_base.py
This is another model where we test how to get parameters from command line.

Simple map coloring

Model: map.py
  • using graph/matrix
  • optimization

Decomposition of global constraint circuit

Model: circuit.py

This is an decomposition of the global constraint circuit which is based on some observation on permutation orbits.
  def circuit(solver, x):  
    n = len(x)
    z = [solver.IntVar(0, n-1, 'z%i' % i) for i in range(n)]

    solver.Add(solver.AllDifferent(x, True))
    solver.Add(solver.AllDifferent(z, True))    

    # put the orbit of x[1] in in z[1..n]        
    solver.Add(z[0] == x[0])
    for i in range(1,n-1):
        # TypeError: list indices must be integers, not IntVar
        # solver.Add(z[i] == x[z[i-1]])

        # using Element
        solver.Add(z[i] == solver.Element(x,z[i-1]))
    # may not be 0 for i < n-1
    for i in range(1, n-1):
        solver.Add(z[i] != 0)
    # when i = n-1 it must be 0
    solver.Add(z[n-1] == 0)
The biggest issue here was that the following (which is the "intuitive" way of stating the constraint):
solver.Add(z[i] == x[z[i-1]])
gives the error TypeError: list indices must be integers, not IntVar

Solution: use solver.Element instead:
solver.Add(z[i] == solver.Element(x,z[i-1]))

Bus scheduling

Model: bus_schedule.py

Here is another problem where we first checks for the optimal value and then find all solutions which this value. This double handling (which is also in the SEND+MOST=MONEY problem) involve some (little) trickery since the minimization is a parameter to solver.Solve:
def main(num_buses_check=0):

    # ...
    if num_buses_check > 0:
        solver.Add(num_buses == num_buses_check)

    collector = solver.AllSolutionCollector(solution)
    cargs = [collector]

    # objective
    if num_buses_check == 0:
        objective = solver.Minimize(num_buses, 1)

    # ....

Note: an alternative way is the following:
solver.NewSearch(db, [objective])

Place Number Puzzle

Model: place_number_puzzle.py

This is a simple graph problem where two connected node can have consecutive numbers. Note: The graph is 1-based, and I compensate for this later in the model.


Model: N/A

I have skipped this problem since there is no built-in support for the regular constraint. Maybe I write a decomposition some day. (It's not very hard but involves some trickeries, and was done for Comet the last year, see Comet: regular constraint, a much faster Nonogram with the regular constraint, some OPL models, and more)

Crossword problem

Model: crossword2.py

  • very simple crossword example
  • using strings/chars combined with int variables
  • a study in non-trivial Element constraint (array/matrix access)
I was just a little surprised that I had to use Element explicitly (and not just accessing the elements via []) . Note that we operate on the flatted version of the A matrix (A_flat).
A = {}
for I in range(num_words):
    for J in range(word_len):
        A[(I,J)] = solver.IntVar(0,26, 'A(%i,%i)' % (I, J))

A_flat = [A[(I,J)] for I in range(num_words) for J in range(word_len)]

# ...

for I in range(num_overlapping):
    # This is what I would do:
    # solver.Add(A[(E[overlapping[I][0]], overlapping[I][1])] ==  
    #            A[(E[overlapping[I][2]], overlapping[I][3])])

    # But we must use Element explicitly

Word square

Model: word_square.py
  • continue from Crosswords
  • another study in non-trivial Element (which is actually simpler than for Crosswords)
  • how to read a file (the word list)
Element again. A and A_flat is defined the same way as in the Crosswords problem (see above).
for i in range(word_len):
    for j in range(word_len):
        # This is what I would like to do:
        # solver.Add(A[(E[i],j)] == A[(E[j],i)])

        # We must use Element explicitly
        solver.Add(solver.Element(A_flat, E[i]*word_len+j) ==
                   solver.Element(A_flat, E[j]*word_len+i))

Combinatorial auction

Model: combinatorial_auction2.py

The most interesting here is probably the use of the scalar product method solver.ScalProd which multiplies the elements of X (an array of Boolean variables) and the bid_amount:
solver.Add(obj == solver.ScalProd(X,bid_amount))
for item in items_t:
   solver.Add(solver.Sum([X[bid] for bid in items_t[item]]) <= 1)

Assignment problem

Model: assignment.py

Here I changed the representation of the variable matrix (x) from the tuple based version to an array based, since it seemed to be easier using the different matrix manipulations such as scalar products, operations on each row and each column:
x = []
for i in range(rows):
    t = []
    for j in range(cols):
       t.append(solver.IntVar(0,1, 'x(%i,%i)' % (i, j)))
The three constraint are then quite easy to state:
# calculate total_cost
solver.Add(total_cost == solver.Sum([solver.ScalProd(x_row,  cost_row) for (x_row, cost_row) in zip(x, cost)]))
# exacly one assignment per row, all rows must be assigned
[solver.Add(solver.Sum([x[row][j] for j in range(cols)]) == 1) for row in range(rows)]

# zero or one assignments per column
[solver.Add(solver.Sum([x[i][col] for i in range(rows)]) <= 1) for col in range(cols)]


Model: discrete_tomography.py

I used the array based matrix represenation here as well, but that was not really needed. Here are the two principal constraints that the rows should sum to the different values in row_sums and col_sums
[solver.Add(solver.Sum([x[i][j] for j in range(c)]) == row_sums[i]) for i in range(r)]
[solver.Add(solver.Sum([x[i][j] for i in range(r)]) == col_sums[j]) for j in range(c)]        
Problem instances:


Model: sicherman_dice.py

Again, we have little more trickeries due to reifications. Here is the constraint that ensures that additions of the values of the two dice is the standard distribution:
standard_dist = [1,2,3,4,5,6,5,4,3,2,1]
n = 6
m = 10
# the two dice
x1 = [solver.IntVar(0, m,"x1(%i)"%i) for i in range(n)]
x2 = [solver.IntVar(0, m,"x2(%i)"%i) for i in range(n)]

for k in range(len(standard_dist)):
    tmp = [solver.BoolVar() for i in range(n) for j in range(n)]
    for i in range(n):
        for j in range(n):
            solver.Add(tmp[i*n+j] == solver.IsEqualCstVar(x1[i] + x2[j],k+2))
    solver.Add(standard_dist[k] == solver.Sum(tmp))
# ...
This can compared how it is stated in another Python based CP system, Numberjack, where we can write the constraint like this, the reification is marked in bold:
# ...
model = Model (
      [ standard_dist[k] == Sum([x1[i] + x2[j] == k+2 for i in range(n) for j in range(n)])
                                                       for k in range(len(standard_dist))],
      # ...

Knapsack problem

Model: knapsack.py

Simple knapsack problem.

Magic square

Model: magic_square.py

Quite straightforward implementation...

Magic square and cards

Model: magic_square_and_cards.py

From Martin Gardner (July 1971):
Allowing duplicates values, what is the largest constant sum for an order-3 magic square that can be formed with nine cards from the deck.
This is implemented as a magic square problem, but without the all different constraint. Instead we must add the constraint that there only 4 cards of each value, using Distribute again:
counts = [solver.IntVar(0,4,"counts(%i)"%i) for i in range(14)]
solver.Add(solver.Distribute(x_flat, range(14), counts))

Subset sum

Model: subset_sum.py A simple subset sum problem.


Model: hidato.py

Well, this takes long time for non-trivial problems. This may be caused by the model is wrong, or ineffecient, or that I have not found any good search strategy. In this problem we have Element and reification:
for k in range(1,r*c):
    i = solver.IntVar(0,r)
    j = solver.IntVar(0,c)
    a = solver.IntVar(-1,1)
    b = solver.IntVar(-1,1)

    # 1) First: fix "this" k
    # 2) and then find the position of the next value (k+1)
    solver.Add(k == solver.Element(x_flat, i*c+j))        
    solver.Add(k + 1 == solver.Element(x_flat, (i+a)*c+(j+b)))        

    solver.Add(i+a >= 0)
    solver.Add(j+b >= 0)
    solver.Add(i+a < r)
    solver.Add(j+b < c)
    # solver.Add(((a != 0) | (b != 0)))
    a_nz = solver.BoolVar()
    b_nz = solver.BoolVar()
    solver.Add(a_nz == solver.IsDifferentCstVar(a,0))
    solver.Add(b_nz == solver.IsDifferentCstVar(b,0))
    solver.Add(a_nz + b_nz >= 1)