Main

March 26, 2012

My talk about Constraint Programming at Google (Paris)

The presentation can be downloaded here: google_talk_20120323.ppt.

In late January this year, I was invited by Laurent Perron - head of the or-tools group - to talk about my view on Constraint Programming and my experience with the or-tools system (I have done quite a few models using the Python, Java, and C# interfaces).

The talk was this Friday (March 23) at the Google's Paris office. It was a lovely day but unfortunately I got a common cold the day before so it was little hard to enjoy all the things Paris can offer.

Friday started with Laurent and me just talking about CP in general, and or-tools in particular and it was really fun and interesting. Later on we where joined by two other guests: Charles Proud'Homme and Nicolas Beldiceanu, both from Ecole des Mines de Nantes and it was great talking with them as well and, above all, listen when they discussed various CP things.

The Google office in Paris was very impressive, very high ceilings and seemed to be build to get lost easily (though neither of us quests got completely lost).

At 1400 I started the talk in front of an audience of about 20 engineers at the Google office (+ the two guests from Ecole des Mines de Nantes) and I think it went quite well considering the cold and all. It was recorded for internal use at Google. I don't know how public it will be but I will blog about this when it has been edited etc. After the 50 minutes talk there was a little Q/A session.

Thanks Laurent for the invitation and a memorable day.

Little more about the talk

The talk was aimed for programmers that don't know very much about Constraint Programming and I especially wanted to convey my own fascination about CP by using this agenda:
  • Introducing the principles of CP (very simplified)
  • Showing the declarativeness of CP by some code in the high level G12 MiniZinc and then in or-tools Python, C#, and sometimes Java.
  • The basic principle of propagation of constraints and domains is shown via a very simple 4x4 Sudoku problem.
  • After that, some of - IMHO - the most fascinating concepts in modeling CP where presented:
    • Global constraints
    • Element constraint
    • Reification
    • Bi-directedness
      Note: After the talk Nicolas Beldiceanu commented that this is more known as "reversibility" in the Prolog world.
    • 1/N/All solutions
    • Symmetry breaking
Here is the talk: google_talk_20120323.ppt.

I would like to thank the following for various degrees of comments, suggestions, and encouragement regarding the presentation:
  • Magnus Bodin
  • Carl Mäsak
  • Mikael Lagerkvist
  • Christian Schulte
  • Laurent Perron
  • Alastair Andrew
And a special thanks to Nikolaj van Omme for his very detailed comments.

March 04, 2012

Talk at Google, Paris (France): My view on Constraint Programming

I'm very honored to announce that I have been invited by Laurent Perron (Google or-tools) to talk about my view on Constraint Programming at Google, Paris later this month.

The talk will be about Constraint Programming and what is so fascinating about it, including - of course - examples in or-tools (which I have tested and modeled in for the Python, Java - and now lately - C# interfaces). There will also be more general views about CP and perhaps comparisons between different CP systems. The audience will be the engineers at Google (Paris) which makes it an extra challenge (and fun). It will also be really great to meet the or-tools group.

After the event I will blog about the details, including publishing the slides.

Also see: My or-tools page.

March 03, 2012

Some newer models (most in MiniZinc, some in or-tools/C#)

For a time I have not blogged about all models I've created, instead just tweeted about them (I'm hakankj at Twitter). And sometimes not even than, just published them directly on my <CP system> pages without any noise.

Well, here I have collect some of these newer and unblogged models in - roughly - time order. Most are in MiniZinc (since I often use MiniZinc for prototyping) but there are also some in other systems. (See Constraint Programming for a list of these CP systems pages.)

  • xkcd_among_diff_0.mzn: Xkcd problem using among_diff_0
    This is another approach of the Xkcd "knapsack problem" where the object is to order dishes for an amount of exactly 15.05.
    This version where inspired by a comment in Helmut Simonis presentation Acquiring Global Constraint Models (page 3), where he use the global constraint among_diff_0 ("Count how many variables are different from 0")
    However, my implementation differs in some ways:
    • it use integers instead of floats.
    • it implements a slightly more general approach
    For more about the global constraint among_diff_0: see my MiniZinc model among_diff_0.mzn. Also see xkcd.mzn, my first approach of the problem.

  • monorail.mzn: Monorail puzzle
    From Aaron Iba's Blog Users of my iOS Game Teach Me a Lesson MIT Didn't:
    The object of Monorail puzzles is to complete a closed-circuit loop through all the stations (dots) by drawing rails. The loop must pass through each station exactly once and close back on itself, like an actual monorail system might in a city.
    Problem:
    
        .   .   .   .       1  2  3  4
    
        .   .___.   .       5  6  7  8
            |
        .   .   .   .       9 10 11 12
    
        .   .___.   .      13 14 15 16
     
    Solution:
        .__.___.___.
        |          |
        .  .___.___.
        |  | 
        .  .___.___.
        |          |
        .__.___.___.
    
    
    Also see
    Problem instances:

  • dennys_menu.mzn
    From Mind Your Decisions (about game theory and personal finance) Denny's math commercial:
    So there’s the question: how many different price combinations will total $10 when menu items priced at $2, $4, $6, and $8?

  • Some Rosetta Code implementations of various knapsack problems


  • Newspaper problem (job-shop)
    Implementations:
    Problem statement from Snehal Patel's CS course
    There are four students: Algy, Bertie, Charlie and Digby, who share a flat. Four newspapers are delivered to the house: the Financial Times, the Guardian, the Daily Express and the Sun. Each of the students reads all of the newspapers, in particular order and for a specified amount of time (see below). Question: Given that Algy gets up at 8:30, Bertie and Charlie at 8:45 and Digby at 9:30, what is the earliest that they can all set off for college?
         Algy           Bertie        Charlie      Digby
    1st  FT       60    Guardian 75   Express  5   Sun      90
    2nd  Guardian 30    Express   3   Guardian 15  FT        1
    3rd  Express   2    FT       25   FT       10  Guardian  1
    4th  Sun       3    Sun      10   Sun      30  Express   1
    
    Extra requirements: All reads the newspaper in a specific order:
    - Algy order   : - FT, Guardian, Express, Sun
    - Bertie order : - Guardian, Express, FT, Sun
    - Charlie order: - Express, Guardian, FT, Sun
    - Digby order  : - Sun, FT, Guardian, Express
    
    This origin of this problem is from S. French: "Sequencing and Scheduling : an introduction to the mathematics of the job-shop", Ellis Horwood Limited, 1982.

    Tim Duncan wrote about it in his paper "Scheduling Problems and Constraint Logic Programming: A Simple Example and its Solution", AIAI-TR-120, 1990, page 5. The paper also includes a program in CHIP solving the problem.)

    Two versions differs by that the first (newpaper0.mzn) is not loaded with so much output stuff as the latter (newpaper.mzn).

  • schedule2.mzn
    Problem from Dennis E. Sasha's book "Puzzles for Programmers and Pros", page 131f:
    In which order do you schedule the tasks starting from current
    day 0?:
    Task  T1 takes 4 days with deadline on day 45
    Task  T2 takes 4 days with deadline on day 48
    Task  T3 takes 5 days with deadline on day 25
    Task  T4 takes 2 days with deadline on day 49
    Task  T5 takes 5 days with deadline on day 36
    Task  T6 takes 2 days with deadline on day 31
    Task  T7 takes 7 days with deadline on day 9
    Task  T8 takes 5 days with deadline on day 39
    Task  T9 takes 4 days with deadline on day 13
    Task T10 takes 6 days with deadline on day 17
    Task T11 takes 4 days with deadline on day 29
    Task T12 takes 1 days with deadline on day 19
    


  • Some Project Euler problems
    Whenever I learn a new programming language, I tend to solve at least the first - say - 20 Project Euler problems. Unfortunately many of the problems requires arbitrary precision or recursive approaches, and neither has good support in MiniZinc. Here are some of the problems:

  • hitting_set.mzn
    From MathWorld: VertexCover
    Let S be a collection of subsets of a finite set X. A subset Y of X that meets every member of S is called the vertex cover, or hitting set. The smallest possible such subset for a given graph G is known as a minimum vertex cover (Skiena 1990, p. 218), and its size is called the vertex cover number, denoted tau(G).
    This model some different problem instances, for example those from the cited article above but also from other sources

    By the way, this model was implemented after I read the paper There is no 16-Clue Sudoku: Solving the Sudoku Minimum Number of Clues Problem by Gary McGuire, Bastian Tugemann, and Gilles Civario. The abstract states: We apply our new hitting set enumeration algorithm to solve the sudoku minimum number of clues problem, which is the following question: What is the smallest number of clues (givens) that a sudoku puzzle may have? [...]
  • maximal_independent_sets.mzn
    From Wikipedia: Maximal independent set:
    In graph theory, a maximal independent set or maximal stable set is an independent set that is not a subset of any other independent set. That is, it is a set S such that every edge of the graph has at least one endpoint not in S and every vertex not in S has at least one neighbor in S. A maximal independent set is also a dominating set in the graph, and every dominating set that is independent must be maximal independent, so maximal independent sets are also called independent dominating sets.

    A graph may have many maximal independent sets of widely varying sizes; a largest maximal independent set is called a maximum independent set.
    The model contains a few problem instances, e.g. those from the above cited Wikipedia article.
    Also see Wikipedia Independent set (graph theory), and compare with the MiniZinc model misp.mzn which use another representation and approach (inspired from GLPK:s model).

  • magic_square_frenicle_form.mzn
    From Wikipedia Frénicle standard form
    A magic square is in Frénicle standard form, named for Bernard Frénicle de Bessy, if the following two conditions apply:
    - the element at position [1,1] (top left corner) is the smallest of the four corner elements; and
    - the element at position [1,2] (top edge, second from left) is smaller than the element in [2,1].
    Activating all these constraints we get the "standard" way of counting the number of solutions:
    (1), 0, 1, 880, 275305224
    
    which is the sequence A006052 from the excellent Online Encyclopedia of Integer Sequence

    Without these symmetry constraints the number of solutions are:
    N  #solutions
    -------------
    1     1
    2     0
    3     8
    4  7040
    5  many...
    
    (Counting the number of the solutions of a CP model is a very good way of ensuring that the model is correct, or rather: if the number of the solutions are not the expected, then the model is wrong.)

  • scheduling_with_assignments.mzn
    This was done directly after I've done the newspaper.mzn (see above). It then struck me that most examples I've seen of scheduling in Constraint Programming (especially when demonstrated the works of cumulative), just showed the times of the jobs not the assignments of the workers.
    In this model, both the job assignments and the worker assignments are shown in different way. For example the solution of one of my own standard problem furniture_moving.mzn (from Marriott and Stuckey: "Programming with constraints", page 112f) is shown as
    earliest_end_time: 60
    num_jobs   : 4
    num_workers: 4
    start_time : [0, 0, 30, 45]
    duration   : [30, 10, 15, 15]
    end_time   : [30, 10, 45, 60]
    resource   : [3, 1, 3, 2]
    allow_idle : true
    collect_workers : false
    do_precendences: false
    
    Assignment matrix (jobs/workers):
    Job1: 1 0 1 1
    Job2: 0 1 0 0
    Job3: 1 0 1 1
    Job4: 1 1 0 0
    
    Assignment matrix (workers/jobs):
    Worker1: 1 0 1 1
    Worker2: 0 1 0 1
    Worker3: 1 0 1 0
    Worker4: 1 0 1 0
    
    Time range for the jobs and the assigned worker(s):
    Job1(0..30): 1 3 4 
    Job2(0..10): 2 
    Job3(30..45): 1 3 4 
    Job4(45..60): 1 2 
    
    ....
    
    Schedule: worker(job), timeline: (earliest_end_time: 60)
    Worker: 1:    1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4 
    Worker: 2:    2  2  2  2  2  2  2  2  2  2  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  -  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4 
    Worker: 3:    1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  -  -  -  -  -  -  -  -  -  -  -  -  -  -  - 
    Worker: 4:    1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  -  -  -  -  -  -  -  -  -  -  -  -  -  -  - 
    
    Time:         0  1  2  3  4  5  6  7  8  9  10 1  2  3  4  5  6  7  8  9  20 1  2  3  4  5  6  7  8  9  30 1  2  3  4  5  6  7  8  9  40 1  2  3  4  5  6  7  8  9  50 1  2  3  4  5  6  7  8  9  
    
    The model presents:
    • start time, duration, and end time for all jobs
    • assignment matrix jobs/workers
    • assignment matrix workers/jobs
    • jobs with time range, and the assigned workers
    • schedule (time line in time units) for the jobs, showing the assigned worker
    • schedule (time line) for the workers, showing the time the workers work
    • and last, a Gantt like schedule showing which job each worker is scheduled to to in each time.

    The model also have some other "bells & whistles" such as
    • handling precedences
    • modeling as a bin pack problem
    • "collecting workers", which may be useful for certain problem, such as perfect square placements

    The two last "features" shows that the job scheduling problem has a family resemblance with bin pack and perfect square placement problem. Unfortunately these are not very fast in this model.

    Well, since I plan to blog about this more, including a benchmark, I leave it for now.

    Here are the problems instances. They has been taken from various sources:

  • equal_sized_groups.mzn
    This is a problem from or-exchange (where many from the area of Operations Research hang out): dividing into roughly equal sized groups, with a sorted list
    I have a problem, and it seems like it should be something that someone has studied before. I have a sorted list of N elements, and I want to divide them into K groups, by choosing K-1 split points between them. There may be elements with the same value, and we want to have items with same value in the same group. Find K groups as close in size to round(N/K) as possible.

    For example, divide these 32 elements in to 4 groups of size 8:
     1 1 1 1 2 2 3 3 3 3 3 3 3 3 4 4 4 4 5 5 5 5 5 6 6 6 6 7 8 9 10 10
    
    One solution would be these 3 break points:
     
     1 1 1 1 2 2 | 3 3 3 3 3 3 3 3 | 4 4 4 4 5 5 5 5 5 | 6 6 6 6 7 8 9 10 10
    [            6                 14                 23                    ]
    [    6 elts      8 elts               9 elts              9 elts          ]
     
     1 1 1 1 2 2         = 6 elements,  error = abs(8-6)=2
     3 3 3 3 3 3 3 3     = 8 elements,  error = abs(8-8)=0
     4 4 4 4 5 5 5 5 5     = 9 elements,  error = abs(8-9)=1
     6 6 6 6 7 8 9 10 10 = 9 elements,  error = abs(8-9)=1
     
     total error = 4
    
    Does this look familiar to anyone? I'd like an approximation algorithm if possible.

    Thanks, Craig Schmidt
    The model contains some examples and results.

  • houses.mzn
    Problem from a Kanren example: houses.scm
    Taken from _Algebra 1_, Glencoe/McGraw-Hill, New York, New York, 1998 pg. 411, Problem 56

    There are 8 houses on McArthur St, all in a row. These houses are numbered from 1 to 8.

    Allison, whose house number is greater than 2, lives next door to her best friend, Adrienne. Belinda, whose house number is greater than 5, lives 2 doors away from her boyfriend, Benito. Cheri, whose house number is greater than Benito's, lives three doors away from her piano teacher, Mr. Crawford. Daryl, whose house number is less than 4, lives 4 doors from his teammate, Don. Who lives in each house?
    One thing to note is the use of the global constraint inverse for channeling the each person to the houses and vice versa.

  • Finding an optimal wedding seating chart
    This problem has been implemented in both MiniZinc and or-tools/C#:
    The problem is from Meghan L. Bellows and J. D. Luc Peterson Finding an optimal seating chart for a wedding (PDF), via Improbable Research Finding an optimal seating chart for a wedding:
    Every year, millions of brides (not to mention their mothers, future mothers-in-law, and occasionally grooms) struggle with one of the most daunting tasks during the wedding-planning process: the seating chart. The guest responses are in, banquet hall is booked, menu choices have been made. You think the hard parts are over, but you have yet to embark upon the biggest headache of them all. In order to make this process easier, we present a mathematical formulation that models the seating chart problem. This model can be solved to find the optimal arrangement of guests at tables. At the very least, it can provide a starting point and hopefully minimize stress and arguments…
    As mentioned before (e.g. in A matching problem, a related seating problem, and some other seating problems) I'm quite fascinated by this type of seating problems.

    And I'm not the only one. After I tweeted about my MiniZinc implementation (wedding_optimal_chart.mzn), Erwin Kalvehagen (with the excellent blog Yet Another Math Programming Consultant showed a GAMS (MIP) model in Weddings and optimal seating . (He also found a bug in my model which was fixed quite easily. Thanks Erwin.)
  • grime_puzzle.mzn
    This problem was taken from the blog Travels in a mathematical world A puzzle from James Grime about abcdef:
    Today James Grime tweeted this question/puzzle:

    Is there a six digit number abcdef such that the following all hold?

    a+b+c+d+e+f = y
    ab+cd+ef=10y
    abc+def=100y


    If not, show why not.

    A little tweeting back and forth verified that "ab" means 10a+b not a×b.

  • balanced_brackets.mzn
    This model generates balanced brackets of size m*2. The number of generated solutions for m:
     m        #
     ----------
     1       1
     2       2
     3       5
     4      14
     5      42
     6     132
     7     429
     8    1430
     9    4862
    10   16796
    11   58786
    12  208012
    13  742900
    
    Which - of course - is the Catalan numbers. See OEIS: 1,2,5,14,42,132,429,1430,4862,16796,58786,208012, and the entry for Catalan numbers: A000108.

  • The "8809 = 6" puzzle
    This problem has been encoded in both MiniZinc and or-tools/C# (and was created yesterday):
    The problem seems to have been around for a couple of years, but it wasn't until I read about it in another excellent blog (Michael Lugo's God Plays Dice): A puzzle that I really tried it (and was lucky to solve it direct without any computational aid). The problem was stated thus:
     
    Here’s a puzzle.
    
    8809 = 6
    7662 = 2
    9312 = 1
    8193 = 3
    8096 = 5
    7756 = 1
    6855 = 3
    9881 = 5
    
    2581 = ?
    
    After a few days a mathematical solution came in the blog post: An answer to a puzzle, but then the problem had been expanded a little:
    8809 = 6
    7111 = 0
    2172 = 0
    6666 = 4
    1111 = 0
    3213 = 0
    7662 = 2
    9312 = 1
    0000 = 4
    2222 = 0
    3333 = 0
    5555 = 0
    8193 = 3
    8096 = 5
    7777 = 0
    9999 = 4
    7756 = 1
    6855 = 3
    9881 = 5
    5531 = 0
    
    2581 = ?
    
    The first version of the problem actually has two different solutions (of the unknown "?", represented as x in the models), since some of the variables are under defined. The second version has a unique solution of x, but there are 10 slightly different solutions since one of the variables is (still) under defined, i.e. the x is the same in all these solutions.

    The approach in my models was inspired by Michael Lugo's mathematical solution (as an equation system), though the or-tools/C# model implements another version using a matrix to represent the problem.

February 05, 2012

A first look at Google or-tools C# interface

In early January Laurent Perron announced that or-tools now supports .NET (C#). A few minutes after I read that, I (and probably many others) wrote to him and asked when there will be a support of Mono (.NET support on Linux, Mac etc). After some week Laurent mailed me to announce that he and Mark Farkas had fixed Mono support as well. Great work!

So, this is what I have done the last weeks: Learning more about C# and modeled some models in or-tools/C#. My models are at my or-tools/C# page, and most of them has been pushed to the project's SVN: code.google.com/p/or-tools/source/browse/trunk/csharp/.

This blog post can be seen as the third part in the "First looks" at Google or-tools tools, where the previous two where A first look at Google CP Solver/Python (Google or-tools) and A first look at Google or-tools' Java interface. Since much of the things in or-tools/C# works the same as for Python and Java I will here mention more about of the differences from these versions.

or-tools/C#

Google or-tools is Google's Operations Research Tools developed at Google and is an open source (though committing and contribution require special handling). The kernel of the solvers is written in C++, and it now supports Constraint Programming (which is - of course - my main interest here) and Mathematical Programming with support for the external solvers Gurubi, GLPK, CBC, and SCIP/Zib Optimization Suite, and modeling interface for C++, Python, Java, and now .NET/C#.

Syntactic sugar - a first example

The first version of or-tools/C# (.NET) was released about January 19, but didn't have much syntactic sugar; instead the models looked very much like the Java models which was kind of lucky for me since it made the porting of these models quite easy. However, after a while Laurent added quite a few shortcuts (i.e. syntactic sugar) which made the modeling much easier. During the following weeks, Laurent added more sugaring, some probably after my complaint, suggestions, or - in some case - confusion.

For example: Here is the first version of least_diff.cs (Least Diff is a simple alphametic optimization problem where the object is to minimize the difference ABCDEF-FGHIJ where A..J is distinct digits). You can see the SVN history of this model here.
using System;
using Google.OrTools.ConstraintSolver;

public class LeastDiff
{
  private static void Solve()
  {
    Solver solver = new Solver("LeastDiff");

    IntVar A = solver.MakeIntVar(0, 9, "A");
    IntVar B = solver.MakeIntVar(0, 9, "B");
    IntVar C = solver.MakeIntVar(0, 9, "C");
    IntVar D = solver.MakeIntVar(0, 9, "D");
    IntVar E = solver.MakeIntVar(0, 9, "E");
    IntVar F = solver.MakeIntVar(0, 9, "F");
    IntVar G = solver.MakeIntVar(0, 9, "G");
    IntVar H = solver.MakeIntVar(0, 9, "H");
    IntVar I = solver.MakeIntVar(0, 9, "I");
    IntVar J = solver.MakeIntVar(0, 9, "J");

    IntVar[] all = new IntVar[] {A,B,C,D,E,F,G,H,I,J};
    int[] coeffs = {10000,1000,100,10,1};
    IntVar x =
        solver.MakeScalProd(new IntVar[]{A,B,C,D,E}, coeffs).Var();
    IntVar y =
        solver.MakeScalProd(new IntVar[]{F,G,H,I,J}, coeffs).Var();
    IntVar diff =
        solver.MakeDifference(x,y).VarWithName("diff");

    solver.Add(solver.MakeAllDifferent(all));
    solver.Add(A > 0);
    solver.Add(F > 0);
    solver.Add(diff > 0);

    OptimizeVar obj = solver.MakeMinimize(diff, 1);

    DecisionBuilder db = solver.MakePhase(all,
                                          Solver.CHOOSE_PATH,
                                          Solver.ASSIGN_MIN_VALUE);

    solver.NewSearch(db, obj);
    while (solver.NextSolution()) {
      Console.WriteLine("{0} - {1} = {2}  ({3}",
           x.Value(), y.Value(), diff.Value(), diff.ToString());
    }

    Console.WriteLine("\nSolutions: {0}", solver.Solutions());
    Console.WriteLine("WallTime: {0}ms", solver.WallTime());
    Console.WriteLine("Failures: {0}", solver.Failures());
    Console.WriteLine("Branches: {0} ", solver.Branches());

    solver.EndSearch();

  }

  public static void Main(String[] args)
  {
    Solve();
  }
}
After Laurent added some syntactic sugar, the constraint and solve part looks like this which in my book is better and easier to work with: all the Make* constructs are removed. Also, equality (==) was sugared after a while but it's not shown here.
// ....
IntVar x = new IntVar[]{A,B,C,D,E}.ScalProd(coeffs).Var();
IntVar y = new IntVar[]{F,G,H,I,J}.ScalProd(coeffs).Var();
IntVar diff = (x - y).VarWithName("diff");

solver.Add(all.AllDifferent());
solver.Add(A > 0);
solver.Add(F > 0);
solver.Add(diff > 0);

OptimizeVar obj = diff.Minimize(1);

DecisionBuilder db = solver.MakePhase(all,
                                      Solver.CHOOSE_PATH,
                                      Solver.ASSIGN_MIN_VALUE);
solver.NewSearch(db, obj);
Here's how to compile and run this model with Mono on a Linux box.
$ export MONO_PATH=$MONO_PATH:/path-to-or-tools/
$ mono-csc /r:Google.OrTools.ConstraintSolver.dll least_diff.cs
$ ./least_diff.exe
....
50123 - 49876 = 247  (diff(247)

Solutions: 30
WallTime: 27ms
Failures: 20
Branches: 83 

Performance comparing to Python and Java versions

Since the Python, Java, and C# models use the same engine there is very little difference between of the performance of models written in these languages, at least when using the same approach and the same heuristics.

As an example, we can continue with the comparison in A first look at Google or-tools' Java interface where the Python and Java versions of a N-queens problem where benchmarked. The results for (the first solution of) N=8, 100, 1000, and 10000 using the C# model nqueens.cs is the same as for the other two implementations:
     N    Failures       Branches      Wall time Required
                                                 memory
  -------------------------------------------------------
     8      3               9            5ms
   100     12             120           10ms 
  1000      0             996          470ms
 10000      0            9999        78923ms     2.1Gb
This is exactly the same result as for Python and Java. Again, this is not surpising since they all use the same underlying solver engine. Comparing Java and C# (using Mono) it seems that the Mono version is slightly faster both in compiling and the time for upstart, though this has not been tested systematic.

The (interesting part of the) model nqueens.cs is:
IntVar[] q = solver.MakeIntVarArray(n, 0, n-1, "q");

solver.Add(q.AllDifferent());

IntVar[] q1 = new IntVar[n];
IntVar[] q2 = new IntVar[n];
for(int i = 0; i < n; i++) {
  q1[i] = (q[i] + i).Var();
  q2[i] = (q[i] - i).Var();
}
solver.Add(q1.AllDifferent());
solver.Add(q2.AllDifferent());

DecisionBuilder db = solver.MakePhase(q,
           Solver.CHOOSE_MIN_SIZE_LOWEST_MAX,
           Solver.ASSIGN_CENTER_VALUE);

LINQ

The nqueens model shown above is a very straightforward implementation, but the interesting part of C# is that we can use LINQ to do array handling in another way using "query comprehensions" (array comprehensions). Although it was one of my objects to learn more about LINQ, I'm not sure if this is better (= more clear) in this particular case. Anyway, here is one way of stating the constraints using LINQ:
solver.Add(q.AllDifferent());
solver.Add((from i in Enumerable.Range(0, n)
            select (q[i] + i).Var()).ToArray().AllDifferent());
solver.Add((from i in Enumerable.Range(0, n)
            select (q[i] - i).Var()).ToArray().AllDifferent());
I think that LINQ shines more in hidato_table.cs (a port of Laurent Perron's Python version using a Table/AllowedAssignments constraint):
  public static int[,] BuildPairs(int rows, int cols)
  {
    int[] ix = {-1, 0, 1};
    var result_tmp = (from x in Enumerable.Range(0, rows)
                      from y in Enumerable.Range(0, cols)
                      from dx in ix
                      from dy in ix
                      where
                      x + dx >= 0 &&
                      x + dx < rows &&
                      y + dy >= 0 &&
                      y + dy < cols &&
                      (dx != 0 || dy != 0)
                   select new int[] {x * cols + y, (x + dx) * cols + (y + dy)}
                      ).ToArray();

    // Convert to len x 2 matrix
    // ....


    return result;
  }
It's not as nice as the Python version, since we must also convert to an int[,] matrix because AllowedAssignments requires this type, but still. (And - of course - there might be another and even simpler version of this that don't require this conversion....)

A side note: one can see how old my C# models are in the use of LINQ, since in later models I tend to use LINQ (and foreach) much more than in earlier models. Also, the models use a quite restricted part of LINQ: almost only from, select, and sometimes where, but not much of the more fancy stuff that LINQ supports. This is because I use it mostly as a layers of handling IntVar's where grouping etc is not appropriate (as I have understood it with some experiments).

MakeIntVarMatrix

One thing that is quite nice - and a sign of Laurent's willingness to help lazy programmers such as myself - is MakeIntVarMatrix. Here is the way one had to do in earlier C#-models to declare an IntVars matrix (i.e. matrix of decision variables) as well as creating an 1D array (x_flat) to be used for the search.
IntVar[,] x = new IntVar[r,c];
IntVar[] x_flat = new IntVar[r*c];
for(int i = 0; i < r; i++) {
  for(int j = 0; j < c; j++) {
    x[i,j] = solver.MakeIntVar(0, 1);
    x_flat[i*c+j] = x[i,j];
  }
} 	
(Note: this is still the way one has to make matrices of decision variables in Python and Java.) However, now one can use these much simpler constructs in C#:
IntVar[,] x = solver.MakeIntVarMatrix(r, c, 0, 1, "x");
IntVar[] x_flat = x.Flatten();
Very handy!

Element

It would not be a complete First look-report if there wasn't any comment about the support of Element in the tested CP system (see for example Learning constraint programming - part II: Modeling with the Element constraint). In one way, how Element is supported can be seen as a proxy of the "expressiveness" of a CP system and its host language. (And this expressiveness is quite important, at least when learning a CP system or perhaps learning CP from the start.)

As with other constraints, Laurent simplified Element after a while. The Who killed Agatha problem is often a good test how easy it is to use Element in a system. The C# model who_killed_agatha.cs has a description of the problem statements.

One of the constraint in this model can be expressed as simple as this in - say - MiniZinc or Comet where hates is a matrix if IntVar, and both the_killer and the_victim are IntVar.
 % MiniZinc/Comet code
 hates[the_killer, the_victim] == 1
This constraint was at first expressed in C# as the code below, which is a rather faithful port from the Java version (WhoKilledAgatha.java):
solver.Add(
   solver.MakeEquality(
      solver.MakeElement(
            hates_flat,
            solver.MakeSum(
            solver.MakeProd(the_killer, n).Var(),
   the_victim).Var()).Var(), 1));
After Laurent's sugaring it can now be stated as the much simpler:
solver.Add(hates_flat.Element(the_killer * n + the_victim) == 1);
This is perhaps as good as one can get when it's not possible to state the constraints in the matrix direct by overloading index access ([,]) of a matrix in the host language. Almost all CP systems embedded in a host programming language requires this kind of handling. (Compare with Gecode that supports constraints directly with the matrix, see Gecode: Modeling with Element for matrices -- revisited; Christian Schulte fixed this support after my complaints during a CP lecture...)

Decompositions

Another thing that it's nice to have in a CP system it to be able to write decompositions that can be used similar as the built-in constraints where - again - simplified syntax helps much. Here is my standard decomposition of the channeling between an IntVar array representing digits and and IntVar variable representing the numbers with these digits (see below for a LINQed version)
private static Constraint ToNum(IntVar[] a, IntVar num, int bbase) {
  int len = a.Length;
  IntVar[] tmp = new IntVar[len];
  for(int i = 0; i < len; i++) {
    tmp[i] = (a[i]*(int)Math.Pow(bbase,(len-i-1))).Var();
  }
  return tmp.Sum() == num;
}
It can now be used much as any constraint, as in pandigital_numbers.cs (yet another alphametic problem). The object is to find all possible multiplications of the form term1 * term2 = product which use all digits. Here we see how the first term (called num1) is calculated:
IntVar[] x = solver.MakeIntVarArray(x_len, start, max_d, "x");
// ...
// First term
solver.Add(ToNum((from i in Enumerable.Range(0, len1)
                  select x[i].Var()).ToArray(),
                  num1, 
                  bbase));
This is also another example how LINQ can simplify the model.

In to_num.cs its use is more direct:
IntVar[] x = solver.MakeIntVarArray(n, 0, bbase - 1, "x");
IntVar num = solver.MakeIntVar(0, (int)Math.Pow(bbase, n) - 1, "num");
solver.Add(x.AllDifferent());
solver.Add(ToNum(x, num, bbase));


As an afterthought: The constraint ToNum could be written like this instead using LINQ:
// ...
int len = a.Length;
return (from i in Enumerable.Range(0, len)
        select (a[i]*(int)Math.Pow(bbase,len-i-1)).Var()
        ).ToArray().Sum() == num;

Mathematical programming

There is also support for Mathematical programming in or-tools (there is also much support for linear programming that I haven't tested much in my earlier or-tools projects). I have just tested this with three simple models which solves almost the same problem using different approaches: then Volsay problem (from OPL).
  • volsay.cs: Volsay LP problem (using "Natural API")
  • volsay2.cs: Volsay LP problem, using arrays for representing the data
  • volsay3.cs: Volsay LP problem, arrays/matrices with components added.
The first model (volsay.cs) looks like this using the "Natural API" for stating the constraints:
// ...
Solver solver = new Solver("Volsay",
                               Solver.CLP_LINEAR_PROGRAMMING);

Variable Gas = solver.MakeNumVar(0, 100000, "Gas");
Variable Chloride = solver.MakeNumVar(0, 100000, "Cloride");
Constraint c1 = solver.Add(Gas + Chloride <= 50);
Constraint c2 = solver.Add(3 * Gas + 4 * Chloride <= 180);

solver.Maximize(40 * Gas + 50 * Chloride);
int resultStatus = solver.Solve();
if (resultStatus != Solver.OPTIMAL) {
  Console.WriteLine("The problem don't have an optimal solution.");
  return;
}
Console.WriteLine("Objective: {0}", solver.ObjectiveValue());
Console.WriteLine("Gas      : {0} ReducedCost: {1}", 
                  Gas.SolutionValue(), 
                  Gas.ReducedCost());
// ...

Search, labeling

Searching and labeling are very much the same in C# as in Python and Java. As can be seen in the Least Diff model above, the C# interface use a simplified way of stating the optimization variable:
OptimizeVar obj = solver.MakeMinimize(diff, 1);
But otherwise it was no surprises compared to the Python and Java models.

Not tested (yet)

There are some things I have yet to test for C#, especially Nonogram using DefaultSearch, where the Python version was quite fast after some tweaking: Google CP Solver: A much faster Nonogram solver using DefaultSearch. Perhaps I will port this to C# as well. Also, there are certain MIP models in or-tools/Python that perhaps will be converted to C#.

Summary

I like this .NET/C# interface to or-tools, much because there is more syntactic sugar than in Java version making it easier to write models, especially when prototyping models; though with Python, been the VHLL language it, is can be even easier. It was quite easy to port from my earlier Java and Python models, especially from the Java models since Java and C# are very similar on the host language level and they use almost exactly the same way of stating the constraint API. Though it took some more time to LINQ'ing these models.

There where very few models where I have any real problems porting from Python/Java, though I had to search for certain things, for example the best way of converting matrices written in Java/Python to C#'s matrices; often it required a change to jagged arrays especially when using Element on a row/column.

When writing models in C# for or-tools I found extremely few bugs, which is a sign of the stability of the support for external programming languages in or-tools.

Also, this was my first serious project to use Mono/C# and I'm impressed how smooth Mono is. This has given me a much better ground for developing more C# models/programs.

Also see

See earlier writings about or-tools, in the category google_cp_solver (which was the name I earlier used for or-tools), and my Google or-tools page which has quite a few Python, Java, and now C# models. My page Common constraint programming problems shows listings of the models that is implemented in other CP systems.

For more information about or-tools:

My or-tools/C# models

Here are all the (about 100) published or-tools/C# models so far. Most of them are ports from the Java or Python models, and most are also at the SVN: code.google.com/p/or-tools/source/browse/trunk/csharp/. Note that Laurent Perron has simplified and/or beautified some of these models. Thanks for this Laurent!

January 03, 2012

Google or-tools supports .NET

Laurent Perron announced yesterday that Google or-tools now has .NET support: .NET support :
Hello,

Thanks to the invaluable help of Mark Farkas, we have .NET support in or-tools.
This has been tested with visual studio 10.

We have 4 new targets:
make csharp
make csharpexe
make cleancsharp
make test_csharp


There are examples in or-tools/csharp and a complete solution in or-tools/csharp/solutions/.


CP Library with arithmetic operators and writing custom decision builder is supported.
LP is supported, but I have not yet written the arithmetic layer as in the python libraries
. Knapsack and graph algorithms are supported.


Please tell me if anything goes wrong.

--Laurent
Some example models are here.

September 09, 2011

Summary of new developments in Google or-tools

Laurent Perron wrote the following summary of the developments in Google or-tools:
Here is what we have been doing in Q3.

Graph algorithms:
- a lot of work went into making the algorithms (min cost flow, max flow) incremental after small modifications.

Linear solver:
- Added utility of constraints, primal and dual tolerances.
- Basis status of variables.
- Added solving model from protobuf, outputting solutions as protobuf. Useful for having MP solvers in another process or machine for instance.

Constraint solver:
- Added complete model visitor
- Applications to pretty print the model or display special stats
- Main application is to implement model read/write to file. I will send another message to detail all this.
- Added model_util binary to query/test/modify exported models.
- Change the internal way of storing constraints and expressions in caches to avoid duplicating equivalent classes.
- Lots of cleanup for the reference manual. Size has been divided by nearly 4 by hiding all implementation classes: check http://or-tools.googlecode.com/svn/doc/full_reference_manual/or-tools/index.html

And of course the usual batch of bug fixes.

--Laurent (on behalf of the complete Operations Research team at Google).

June 22, 2011

News from Google or-tools (CP solver)

Laurent Perron mailed two announcement about Google or-tools to the or-tools-discuss (discussion group): New developments
Hello,

It is time for a status update on our code!

Since the last announcements, we have added:
  • linear assignment (including dimacs challenge support): A. V. Goldberg and R. Kennedy, "An Efficient Cost Scaling Algorithm for the Assignment Problem." Mathematical Programming, Vol. 71, pages 153-178, December 1995.
  • Min cost flow: R.K. Ahuja, A.V. Goldberg, J.B. Orlin, and R.E. Tarjan, "Finding minimum-cost flows by double scaling," Mathematical Programming, (1992) 53:243-266. http://www.springerlink.com/index/gu7404218u6kt166.pdf
  • Max flow (many references, see graph/max_flow.h for details).
  • SCIP support (see scip.zib.de). We will make a separate announcement on how to compile scip to be used in or-tools. [See below.]
  • Deviation constraint in the Constraint Programming solver : Pierre Schaus et. al., "Bound Consistent Deviation Constraint", CP07.
  • Initial support for no good management in the CP search tree.
  • 30-50% speedup on large sums in constraint programming.
--Laurent (on behalf of the complete Operations Research team at Google).
Regarding SCIP support, Laurent also posted Using SCIP with or-tools.

Related:
Some days ago I posted a comment on or-tools-discuss about how to run Google or-tools on Linux Ubuntu 11.04: can I use or-tools with Ubuntu 11.04? which may be of some help.

March 31, 2011

Google or-tools has released support for LP/IP

Today Google or-tools released support for linear programming (LP) and integer programming (IP). Compiled versions of GLPK and CBC can be downloaded here.

Some small examples in
  • C++, integer_solver_example.cpp, linear_solver_example.cpp
  • Python, integer_solver_example.py, linear_solver_example.py
  • Java, IntegerSolverExample.cpp, LinearSolverExample.cpp

Coins grid LP version

Just to test this I converted one of my CP models coins_grid.py to a IP model: coins_grid_mip.py. This was quite straightforward, except that some constrains are not exactly like the CP interface. The model solves the problem in no time.

Here is the full model (edited)
from linear_solver import pywraplp
def main(unused_argv):
    # using GLPK
  solver = pywraplp.Solver('CoinsGridGLPK',
                          pywraplp.Solver.GLPK_LINEAR_PROGRAMMING)

  # Using CLP
  # solver = pywraplp.Solver('CoinsGridCLP',
  #                          pywraplp.Solver.CLP_LINEAR_PROGRAMMING)

  n =  31  # the grid size
  c =  14  # number of coins per row/column

  x = {}
  for i in range(n):
    for j in range(n):
      x[(i,j)] = solver.IntVar(0, 1, 'x[%i,%i]' % (i, j))

  # sum rows/columns == c
  for i in range(n):
    solver.Add(solver.Sum(
        [x[(i, j)] for j in range(n)]) == c)      # sum rows
    solver.Add(solver.Sum(
        [x[(j, i)] for j in range(n)]) == c) # sum cols

  # quadratic horizonal distance var
  objective_var = solver.Sum(
        [x[(i, j)] * (i - j) * (i - j)
          for i in range(n) for j in range(n)])

  objective = solver.Minimize(objective_var)

  solver.Solve()

  for i in range(n):
    for j in range(n):
      # int representation
      print int(x[(i, j)].solution_value()),
    print
  print

if __name__ == '__main__':
    main("coin grids")
Note: the solver outputs floats (0.0 and 1.0) but I strip this to 0 and 1 in the output.

Installation

On my Linux Ubuntu box I have to recompile the source code (from SVN version) to get this to work. Before that some things had to be changed.

Add the following two directories to LD_LIBRARY_PATH so the solver can find the .so files
  • (some path)/glpk-4.45/lib
  • (some path)/cbc-2.6.2/lib/coin
some_path is the directory where the GLPK and CBC packages are installed. Also, one have to set the proper path to GLPK and CBC in the Makefile.

After that everything worked fine.

For testing the Java models, type the following in the or-tools directory (where the Makefile is):
make run_LinearSolverExample
make run_IntegerSolverExample

March 29, 2011

A first look at Google or-tools' Java interface

The Java interface for Google or-tools was released some week ago, and I have now ported some of my Python models to Java (about 35 models). See my Google CP Solver page / Java models for a list of them. Some of them has also been committed to the SVN repository of the or-tools project.

In A first look at Google CP Solver/Python (Google or-tools) I described the Python interface for Google or-tools. Since much of what is written there is applicable with the Java interface as well, I here concentrate on the differences between these two versions.

The porting was surprisingly easy to do for most models, since all the work of looking up the proper API call to use was already done for the Python models. See below for some of the differences.

I found very few unexpected things in the Java version. The Python interface (being slightly older) had some bugs when I first tested it, but they where fixed very fast by Laurent Perron and his colleagues.

Example: SEND + MORE = MONEY

One difference between the Java interface and the Python interface is that the Python interface has some more syntactic sugar, e.g. it recognizes the arithmetic operators (+, -, /, *), and (dis)equality: ==, !=. The Java interface has none of this which caused me some pain to translate.

For example: To get the result of two numbers (a + b == c) in Python:
  solver.Add(a + b == c);
In the Java interface one has to do use two methods (addSum and makeEquality):
  solver.addConstraint(solver.makeEquality(solver.makeSum(a, b), c));
Here is a table of corresponding operations:
Python       Java
----------------------------
a + b        makeSum(a, b)
a * b        makeProd(a, b)
a - b        makeDifference(a, b)
- a          makeOpposite(a)
a / b        makeDiv(a, b)
a == b       makeEquality(a, b)
a != b       makeDifference(a, b)
Example: the usual core of the SEND + MORE == MONEY problem can be written in a Python model (and similar in other constraint solver systems):
# define variables etc
# ...
#
solver.Add(      s*1000 + e*100 + n*10 + d +
                 m*1000 + o*100 + r*10 + e ==
           m*10000 + o*1000 + n*100 + e*10 + y
However, using this approach is not really doable with the Java interface, since it would require a lot of makeSum and makeProd, etc. A better way (and this is probably a general advice for the Python interface as well) is to instead use a scalar product: makeScalProd or makeScalProdEquality. Here is one way, where we also see how the decision variables is declared.

Solver solver = new Solver("SendMoreMoney");

int base = 10;
IntVar s = solver.makeIntVar(0, base - 1, "s");
IntVar e = solver.makeIntVar(0, base - 1, "e");
IntVar n = solver.makeIntVar(0, base - 1, "n");
IntVar d = solver.makeIntVar(0, base - 1, "d");
IntVar m = solver.makeIntVar(0, base - 1, "m");
IntVar o = solver.makeIntVar(0, base - 1, "o");
IntVar r = solver.makeIntVar(0, base - 1, "r");
IntVar y = solver.makeIntVar(0, base - 1, "y");

IntVar[] x = {s,e,n,d,m,o,r,y};

IntVar[] eq = {s,e,n,d,  m,o,r,e, m,o,n,e,y};
int[] coeffs = {1000, 100, 10, 1,            //    S E N D +
                1000, 100, 10, 1,            //    M O R E
                -10000, -1000, -100, -10, -1 // == M O N E Y
};
solver.addConstraint(
      solver.makeScalProdEquality(eq, coeffs, 0));
// s > 0
solver.addConstraint(
      solver.makeGreater(s, 0));
// m > 0
solver.addConstraint(
      solver.makeGreater(m, 0));
// all different
solver.addConstraint(
      solver.makeAllDifferent(x, true));
This approach, using scalar products, is implemented in SendMoreMoney.java. In SendMoreMoney2.java I experimented with some different approaches of coding this problem, some quite hairy.

Search, labeling, output

There is not really much difference how to do search and labeling in the Java version compared to the Python version. Here is the search section of SendMoreMoney.java:
DecisionBuilder db = solver.makePhase(x,
                                       solver.INT_VAR_DEFAULT,
                                       solver.INT_VALUE_DEFAULT);
solver.newSearch(db);
while (solver.nextSolution()) {
   for(int i = 0; i < 8; i++) {
     System.out.print(x[i].toString() + " ");
   }
   System.out.println();
 }
solver.endSearch();

// Statistics
System.out.println();
System.out.println("Solutions: " + solver.solutions());
System.out.println("Failures: " + solver.failures());
System.out.println("Branches: " + solver.branches());
System.out.println("Wall time: " + solver.wall_time() + "ms");
One thing that is different, is how to define the variables for labeling. In Python it is quite easy to concatenate two arrays of decision variables like this:
solver.makePhase(x + y, ...);
But this don't work as easy in Java (being Java). I tend to call such a labeling array for all and fill in all the IntVar from x and y; there are other way of doing this of course. On the other hand, the Python interface cannot handle labeling on matrices (multidimensional tuples), and in these models I also had to rely on a flattened "view" of this matrix (often called x_flat or something like that).

Python also has niftier shortcuts for creating temporary arrays using array comprehensions. For example in the Python version of set_covering.py one can write the following:
# ensure that all cities are covered
for i in range(num_cities):
   b = [x[j] for j in range(num_cities) if distance[i][j] <= min_distance]
   solver.Add(solver.SumGreaterOrEqual(b, 1))
This I have translated to Java code using an ArrayList instead (SetCovering.java):
// ensure that all cities are covered
for(int i = 0; i < num_cities; i++) {
  ArrayList b = new ArrayList();
  for(int j = 0; j < num_cities; j++) {
    if (distance[i][j] <= min_distance) {
      b.add(x[j]);
    }
  }
  solver.addConstraint(
      solver.makeSumGreaterOrEqual(b.toArray(new IntVar[1]), 1));
}

Decompositions, and somewhat under the hood

As far as I know, there is not (yet) support for modulo operations with decision variables (IntVar's). Since modulo is a favorite operator of mine, I implemented a version in DivisibleBy9Through1.java which corresponds quite closely to the Python version divisible_by_9_through_1.py. Here is the Python version:
def my_mod(solver, x, y, r):

    if not isinstance(y, int):
        solver.Add(y != 0)

    lbx = x.Min()
    ubx = x.Max()
    ubx_neg = -ubx
    lbx_neg = -lbx
    min_x = min(lbx, ubx_neg)
    max_x = max(ubx, lbx_neg)

    d = solver.IntVar(min_x, max_x, 'd')    

    if not isinstance(r, int):
        solver.Add(r >= 0)
        solver.Add(x*r >= 0)

    if not isinstance(r, int) and not isinstance(r, int):
        solver.Add(-abs(y) < r)
        solver.Add(r < abs(y))
        
    solver.Add(min_x <= d)
    solver.Add(d <= max_x)
    solver.Add(x == y*d+r)
The Java version is not very different, however I simplified it somewhat and requires that all arguments are IntVar's:
public static void my_mod(Solver solver, IntVar x, IntVar y, IntVar r) {
        
  long lbx = x.min();
  long ubx = x.max();
  long ubx_neg = -ubx;
  long lbx_neg = -lbx;
  int min_x = (int)Math.min(lbx, ubx_neg);
  int max_x = (int)Math.max(ubx, lbx_neg);

  IntVar d = solver.makeIntVar(min_x, max_x, "d");
    
  // r >= 0
  solver.addConstraint(
        solver.makeGreaterOrEqual(r,0));

  // x*r >= 0
  solver.addConstraint(
      solver.makeGreaterOrEqual(
          solver.makeProd(x,r).Var(), 0));

  // -abs(y) < r
  solver.addConstraint(
      solver.makeLess(
          solver.makeOpposite(
              solver.makeAbs(y).Var()).Var(), r));

  // r < abs(y)
  solver.addConstraint(
      solver.makeLess(r,
          solver.makeAbs(y).Var().Var()));

  // min_x <= d, i.e. d > min_x
  solver.addConstraint(
        solver.makeGreater(d, min_x));


  // d <= max_x
  solver.addConstraint(
        solver.makeLessOrEqual(d, max_x));

  // x == y*d+r
  solver.addConstraint(
        solver.makeEquality(x,
            solver.makeSum(
               solver.makeProd(y,d).Var(),r).Var()));
}
Note: This implementation is based on the ECLiPSe version mentioned in A Modulo propagator for ECLiPSE. Here is the ECLiPSe Prolog source code: modulo_propagator.ecl.

Another example is the decomposition of the global constraint circuit (from Circuit model) based on some observations about the "orbit" of x.
public static void circuit(Solver solver, IntVar[] x) {

  int n = x.length;
  IntVar[] z = solver.makeIntVarArray(n, 0, n - 1, "z");

  solver.addConstraint(
      solver.makeAllDifferent(x, true));
  solver.addConstraint(
      solver.makeAllDifferent(z, true));

  // put the orbit of x[0] in z[0..n-1]
  solver.addConstraint(
        solver.makeEquality(z[0], x[0]));
  for(int i = 1; i < n-1; i++) {
    solver.addConstraint(
        solver.makeEquality(z[i], 
            solver.makeElement(x, z[i-1]).Var()));
  }

  // z may not be 0 for i < n-1
  for(int i = 1; i < n - 1; i++) {
    solver.addConstraint(
          solver.makeNonEquality(z[i], 0));
  }
        
  // when i = n-1 it must be 0
  solver.addConstraint(
        solver.makeEquality(z[n - 1], 0));

}

Some benchmark: n-queens

Below is some statistics for the the n-queens problem, implemented in NQueens2 model. The constraint section is:
IntVar[] q = solver.makeIntVarArray(n, 0, n-1, "q");

solver.addConstraint(solver.makeAllDifferent(q, true));

IntVar[] q1 = new IntVar[n];
IntVar[] q2 = new IntVar[n];
for(int i = 0; i < n; i++) {
  q1[i] = solver.makeSum(q[i], i).Var();
  q2[i] = solver.makeSum(q[i], -i).Var();
}
solver.addConstraint(
    solver.makeAllDifferent(q1, true));
solver.addConstraint(
    solver.makeAllDifferent(q2, true));
Here is the timings (for one solution):
     N    Failures       Branches      Wall time Required
                                                 memory
  -------------------------------------------------------
     8      3             9              2ms
    10      5            16              3ms
    50     15            70              4ms
   100     12           120             10ms
   200      1           198             22ms
   500      2           496            126ms
   501      6           503            128ms
  1000      0           996            504ms 
  2000      0          1995           2164ms
  3000      2          3001           5119ms
  4000      0          3999           9493ms
  5000      1          4994          15803ms
  6000      1          5996          24511ms
  7000      0          6997          34354ms   1.1Gb
  8000      0          7996          46985ms   1.5Gb
  9000      0          8996          62905ms   1.8Gb
 10000      9          9999          81312ms   2.3Gb
I ran these tests on a Linux Ubuntu 10.4 Lynx, with 12Gb RAM, 8 cores, (Intel Core i7 930 @ 2.80GHz). Note: the "Wall time" is the time reported by the solver, and the "Required memory" is manually via the "top" command for the running process.

I wondered if the Python version nqueens3.py that use the same constraints and labeling as the Java version should have any different result. Well, it has exactly the same failures and branches as the Java version (no surprise since it use the same solver engine). Also the wall times and memory usage are almost identical.

Related

See also these blog posts, all about the Python interface:

Some kind of conclusion

As mentioned above, it was quite easy to port the Python models to Java since they share a lot of common API (with some difference in details). If there is no requirement about the language to use, I personally would probably prefer the Python interface since it's more high level and is less verbose than the Java interface (I like brevity).

Maybe (repeat maybe) I will also have a take on Google or-tools' C++ interface.

My Google or-tools Java models

Here are my Java models. All of them where ported from existing Python models. They are also available at Google CP Solver page / Java models (the Python models is available from the same page). For implementation of these problems in other constraint programming systems, see Common constraint programming problems.

November 11, 2010

Google CP Solver: A much faster Nonogram solver using DefaultSearch

This is an extension of the comparison done in Comparison of some Nonogram solvers: Google CP Solver vs Gecode and two MiniZinc solvers with includes a new and much faster variant in Google CP Solver.

The day after the above blog post was written, Laurent Perron and I discussed alternative ways of solving these Nonogram problems. Laurent suggested another solving strategy: DefaultSearch (solver.DefaultPhase) which gave excellent overall result. Then I did some more experiment and tweaking which gave the results shown in the table below.

The new Nonogram model

Model: nonogram_default_search.py

The changes compared to the previous version nonogram_table2.py was very small in the Python program (there where also some changes in the C++ code). Here are the changes:
parameters = pywrapcp.DefaultPhaseParameters()
parameters.heuristic_period = 200000
db = solver.DefaultPhase(board_label, parameters)  

DefaultSearch

To be honest, I am not exactly sure why DefaultSearch is so fast for the Nonogram problems. It is a variant of Impact-Based Search (see Philippe Refalo: Impact-Based Search Strategies for Constraint Programming, 2004) combined with periodic heuristics. The source code is here: constraint_solver/default_search.cc.

heuristic_period

The value of 200000 for parameters.heuristic_period was found by experimenting. For quite long I used the value of 11000 and later found that even though increasing this value didn't do much difference for the simpler problems (same time and same number of failures), for the harder problems - such as Lion, Marley, and Nature - the number of failures decreased and the time was improved a bit. Larger values above that, such as 2000000, decreased the number of failures somewhat for the Nature problem, but the time seemed to be slightly worse. So I settled with 200000.

This may or may not be a good value for general Nonogram problems. Doing this kind of testing on a limited - albeit mixed - problem set may have biased the search in some unfortunate way.

After deciding on the above heuristics, I finally tested (as a verification set) with some other problem instances and found that they still was easy solved:
  • P200: 0.25s / 1637 failures
  • P199: 0.06s / 242 failures
  • Dragonfly: 0.03s / 0 failures
P200 was slightly worse in this version: the previous version solved it in 0.14s and with 520 failures. The other two problem was better or about the same with the current version.

As an experiment, I changed the order of variables to board_flat (instead of board_label), then Lion was solved in 2 seconds (compared to 4.45 seconds) and Marley in 11 seconds (compared to 7:18 minutes). On the down side, 9-Dom then took 9 seconds (compared to 3.5 seconds) and Nature timed out (compared to 20:25 minutes), and other problems took slightly longer time. So this is probably not a good way to go.

I also experimented with restarts (which is reported to work well with Impact-Based Search), and found some initial promising results. I may come back to this later.

I should also note that the current model is quite different from my first model nonogram_regular.py presented in Google CP Solver: Regular constraint, Nonogram solver, nurse rostering, etc. Laurent has done a great job of optimizing all sort of things, both in the Python code and in the C++ code. For example, he added the "regular like" constraint solver.TransitionConstraint which obsoleted my decomposition of regular for this problem.

A comment on the results

None of the Nonogram solvers that Jan Wolter had tested in his excellent Survey of Paint-by-Number Puzzle Solvers has solved all the problem instances under the 30 minutes time limit. The table below shows that the DefaultSearch variant solves all the problems. At least on my computer, and hopefully also on Wolter's benchmark computer if he want to test this solver. All problems - except Marley and Nature - where solved well under 10 seconds, including the infamous Lion problem (4.46s).

In the previous version of the Nonogram solver a couple of the problems timed out, marked + in the table. Some of these, such as Hot and Light, is now solved under 2 seconds. All these "remarkable" improvements in bold.

Test machine

The machine used for the test is the same as in the former test: Linux Ubuntu 10.4 LTS with 12Gb RAM, 8 core (Intel i7/930, 2.80GHz). All programs ran on a single processor. Python version 2.6. Google CP Solver used was revision 276.

Comparison - table

Here is the full table with the new version (DefaultSearch) in the last column. The number in parenthesis is Wolter's times. The second row in each entry is my own timing. The number of failures where applicable; LazyFD always returned 0 choice points so these are skipped. A + indicates time out (30 minutes).

The links for the puzzle is to Wolter's puzzle pages.
Puzzle Size Notes Gecode MiniZinc
Gecode/fz
MiniZinc
LazyFD
Google
CP Solver
Google
CP Solver
DefaultSearch
#1: Dancer* 5 x 10 Line (0.01s)
0.01s/0 failures
(0,01s)
0.08s/0 failures
(0.04s)
0.1s
0.02s
0 failures
0.02s
0 failures
#6: Cat* 20 x 20 Line (0.01s)
0.1s/0 failures
(0.02s)
0.1s/0 failures
(0.44s)
0.95s
0.12s
0 failures
0.03s
0 failures
#21: Skid* 14 x 25 Line, Blanks (0.01s)
0.01s/0 failures
(0,02s)
0.09s/13 failures
(0.64s)
0.9s
0.03s
0 failures
0.03s
0 failures
#27: Bucks* 27 x 23 Blanks (0.01s)
0.2s/2 failures
(0.02s)
0.1s/3 failures
(1.1s)
1.6s
0.05s
3 failures
0.04s
6 failures
#23: Edge* 10 x 11   (0.01s)
0.01s/15 failures
(0.01s)
0.09s/25 failures
(0.08s)
0.16s
0.03s
25 failures
0.03s
58 failures
#2413: Smoke 20 x 20   (0.01s)
N/A
(0.02s)
0.11s/5 failures
(0.60s)
1.18s
0.05s
8 failures
0.04s
20 failures
#16: Knot* 34 x 34 Line (0.01s)
0.01s/0 failures
(0.02s)
0.14s/0 failures
(5.5s)
5.6s
0.12s
0 failures
0.06s
0 failures
#529: Swing* 45 x 45 Line (0.02s)
0.02s/0 failures
(0.04s)
0.24s/0 failures
(13.2s)
22.3s
0.3s
0 failures
0.11s
0 failures
#65: Mum* 34 x 40   (0.02s)
0.02s/17
(0.04s)
0.18s/20 failures
(11.0s)
12.4s
0.15s
22 failures
0.10s
46 failures
#7604: DiCap 52 x 63   (0.02s)
N/A
(0.06s)
0.29s/0 failures
(+)
1:48m
0.45s
0 failures
0.18s
0 failures
#1694: Tragic 45 x 50   (0.14s)
N/A
(3.6m)
2:14m/394841 failures
(32.1s)
34.57s
2:33m
394841 failures
0.35s
1862 failures
#1611: Merka* 55 x 60 Blanks (0.03s)
N/A
(+)
+
(1.1m)
1:10m
0.47s
27 failures
0.20s
96 failures
#436: Petro* 40 x 35   (1.4m[s?])
0.05s/48 failures
(1.6s)
1.15s/1738 failures
(6.2s)
9.3s
1.19s
1738 failures
0.245s
770 failures
#4645: M&M 50 x 70 Blanks (0.93s)
N/A
(0.28s)
0.41s/89 failures
(48.1s)
1:14m
0.48s
82 failures
0.46s
376 failures
#3541: Signed 60 x 50   (0.57s)
N/A
(0.56s)
0.61s/929 failures
(1.0m)
1:02m
11.7s
6484 failures
0.92s
3821 failures
#803: Light* 50 x 45 Blanks (+)
N/A
(+)
+
(4.7s)
18.32s
+ 1.5s
4052 failures
#6574: Forever* 25 x 25   (4.7s)
N/A
(2.3s)
1.5s/17143 failures
(1,7s)
1.99s
1.6s
17143 failures
0.26s
2421 failures
#2040: Hot 55 x 60   (+)
N/A
(+)
+
(2.6m)
1:05m
+ 0.62s
2394 failures
#6739: Karate 40 x 40 Multiple (56.0s)
N/A
(57.8s)
38.0s/215541 failures
(13.6s)
22.52s
56.1s
215541 failures
0.2s
1361 failures
#8098: 9-Dom* 19 x 19   (12.6m)
N/A
(3.8s)
2.59s/45226 failures
(1.8s)
3.11s
3.9s
45226 failures
3.54s
27650 failures
#2556: Flag 65 x 45 Multiple, Blanks (3.4s)
N/A
(+)
+
(4.2s)
16.18s
1.6s
14859 failures
0.15s
442 failures
#2712: Lion 47 x 47 Multiple (+)
N/
(+)
+
(+)
8:54m
+ 4.34s
15429 failures
#10088: Marley 52 x 63 Multiple (+)
N/A
(+)
+
(2.9m)
+
+ 7:18m
309293 failures
#9892: Nature 50 x 40 Multiple (+)
N/A
(+)
+
(2.4m)
25.29s
+ 20:25m
7743891 failures

November 07, 2010

Google CP Solver: Regular constraint, Nonogram solver, nurse rostering, etc

Models mentioned: After some sweat and confusions, I have now implemented a Nonogram solver in Google CP Solver using a (decomposition of) regular constraint. The model is here: nonogram_regular.py. Also, a test model for the regular constraint is in regular.py

The hardest part by far was to translate the regular constraint from Comet to Google CP solver. The Comet version (nonogram_regular.co) was described in a sequence of blog posts: As we will see below, other problems where implemented by the regular constraint.

Update: Laurent Perron suggested an improvement of the regular constraint, which makes some of the models much faster. I have commented this below.

Translation of regular

As mentioned above, the hardest part was the translation of the regular constraint from Comet (from nonogram_regular.co) to Google CP solver (the Comet version was a translated from the MiniZinc version). I started with the simpler example regular.py) which basically consists of two parts also included in the nonogram model:

make_transition_matrix

The function translates an array of numbers to an finite state automaton which representing and regular expression of 0's and 1's. E.g. the array [3,2,1] represents the regular expression 0*1{3}0+1{2}0+1{1}0* and is translated into the automaton (the valid states) which must start in state 1 and end in state 9 (0 is a failing state):
1 2
0 3
0 4
5 0
5 6
0 7
8 0
8 9
9 0
The hardest part here was that the Comet constraint is 1-based, but Python is 0-based. There are many things that can - and will - go wrong when translating from a 1-based system to a 0-base system. And it don't help that state 0 is used as a fail state. regular constraint

Then, given the automaton (the valid states), the regular constraint transforms the states into an array of decision variables. The domains are 1..2 (representing the values of 0 and 1).

Update Laurent Perron suggested an improvement using the table constraint (called solver.AllowedAssignments in Google CP Solver). I have changed to this version here.

This was his first improvement (last in the method):
#
# Global constraint regular
#
# This is a translation of MiniZinc's regular constraint (defined in
# lib/zinc/globals.mzn), via the Comet code refered above.
# All comments are from the MiniZinc code.
# '''
# The sequence of values in array 'x' (which must all be in the range 1..S)
# is accepted by the DFA of 'Q' states with input 1..S and transition
# function 'd' (which maps (1..Q, 1..S) -> 0..Q)) and initial state 'q0'
# (which must be in 1..Q) and accepting states 'F' (which all must be in
# 1..Q).  We reserve state 0 to be an always failing state.
# '''
#
# x : IntVar array
# Q : number of states
# S : input_max
# d : transition matrix
# q0: initial state
# F : accepting states
def regular(x, Q, S, d, q0, F):

  solver = x[0].solver()

  assert Q > 0, 'regular: "Q" must be greater than zero'
  assert S > 0, 'regular: "S" must be greater than zero'

  # d2 is the same as d, except we add one extra transition for
  # each possible input;  each extra transition is from state zero
  # to state zero.  This allows us to continue even if we hit a
  # non-accepted input.

  d2 = []
  for i in range(Q + 1):
    for j in range(S):
      if i == 0:
        d2.append((0, j, 0))
      else:
        d2.append((i, j, d[i - 1][j]))

  # If x has index set m..n, then a[m-1] holds the initial state
  # (q0), and a[i+1] holds the state we're in after processing
  # x[i].  If a[n] is in F, then we succeed (ie. accept the
  # string).
  x_range = range(0, len(x))
  m = 0
  n = len(x)

  a = [solver.IntVar(0, Q + 1, 'a[%i]' % i) for i in range(m, n + 1)]

    # Check that the final state is in F
  solver.Add(solver.MemberCt(a[-1], F))
    # First state is q0
  solver.Add(a[m] == q0)
  for i in x_range:
    solver.Add(x[i] >= 1)
    solver.Add(x[i] <= S)
    # Determine a[i+1]: a[i+1] == d2[a[i], x[i]]
    # Laurent Perron suggested using this: 
    solver.Add(solver.AllowedAssignments((a[i], x[i] - 1, a[i + 1]), d2))
Update 2: Some hour later, Laurent Perron added C++ support for simplifying this by adding solver.TransitionConstraint, so the definition of regular is now much neater:
def regular(x, Q, S, d, q0, F):

  solver = x[0].solver()

  assert Q > 0, 'regular: "Q" must be greater than zero'
  assert S > 0, 'regular: "S" must be greater than zero'

  # d2 is the same as d, except we add one extra transition for
  # each possible input;  each extra transition is from state zero
  # to state zero.  This allows us to continue even if we hit a
  # non-accepted input.

  d2 = []
  for i in range(Q + 1):
    for j in range(1, S + 1):
      if i == 0:
        d2.append((0, j, 0))
      else:
        d2.append((i, j, d[i - 1][j - 1]))

  solver.Add(solver.TransitionConstraint(x, d2, q0, F))
This improvement is in the model: regular_table2.py and the Nonogram model: nonogram_table2.py

Nonogram

After some more problems was fixed it was really a simple matter of translating the rest of the Nonogram model (nonogram_regular.py).

Labeling: I has the same labeling, or as approximate as possible, as in the Comet model:
# This labeling was inspired by a suggestion from
# Pascal Van Hentenryck about my Comet nonogram model.
board_label = []
if rows * row_rule_len < cols * col_rule_len:
    for i in range(rows):
        for j in range(cols):
            board_label.append(board[i,j])
else:
    for j in range(cols):        
        for i in range(rows):
            board_label.append(board[i,j])
Update: Here are Laurent Perron's versions with the improved regular: nonogram_table.py and nonogram_table2.py

Testing and problem instances for Nonogram

I have compared two problem instances, p200 (nonogram_p200.py) and p199 (nonogram_p199.py), with some of my other Nonogram solver: the Comet model nonogram_regular.co (with my decomposition of regular), the MiniZinc model nonogram_create_automaton2.mzn with Gecode/fz, using Gecode's built in regular.

Update: I have added Laurent Perron's version using table (solver.AllowedAssignments) as an separate entry for comparison:

Problem Google CP Solver Google CP Solver with table Comet Gecode/fz
P200 1.8s/11615 failures 0.2s/520 failures 0.4s/794 failures 0.2s/520 failures
P199 17.607s/104239 failures 0.3s/772 failures 0.6s/1175 failures 0.1s/772 failures

The result for P200 is not too bad for Google CP Solver, but I am a little surprised by the large times/failures for the P199, almost 18 seconds. However, one could note that using a plain ordering of the variables for this problem instance, i.e.
board_label = [board[i,j] for i in range(rows) for j in range(cols)]
then P199 is solved in 0.4 seconds with 2005 failures. However, with this ordering, P200 takes very long time.

Update: It's interesting to see that the improved table version has exactly the same number of failures as the version using Gecode/fz.
Update 2:The second improvement of Laurent Perron didn't change the number of failures, but the time was slightly improved with a couple of ms's: 169ms vs 209ms for P200, and 264ms vs 296ms for P199.

Thanks for the improvements, Laurent!

Here are some other problem instances:

Nurse rostering

Model: nurse_rostering.py

A standard application of regular constraint is in different types of scheduling, e.g. nurse rostering. The model nurse_rostering.py is a small example based on the model in the MiniZinc Tutorial (PDF), page 35f (the graph of the DFA is shown at page 37).

The example is for 7 nurses and 14 days using three states: day shift, night shift, and day off. There are two kind of constraints:

* Rules using the regular constraint:
  • - one day off every 4 days
  • - no 3 nights in a row.
These are translated as the following transition function:
transition_fn =  [
   # d,n,o
    [2,3,1], # state 1
    [4,4,1], # state 2
    [4,5,1], # state 3 
    [6,6,1], # state 4
    [6,0,1], # state 5
    [0,0,1]  # state 6
    ]
Note: This transition function can be used for any number of nurses and days. * Other rules
Here I went wild and required some extra rules for this specific number of nurses and days:
  • Each nurse must work between 7 and 10 days during the period (14 days)
  • On weekends there must be 2 nurses on day shift, 1 on night shift, and 4 has day off
  • On working days: 3 nurses on day shift, 2 on night shift, and 2 has day off
Below is one solution, including statistics for each nurse and day, respectively:
Nurse0:  d d d o d d d o d d d o d o  day_stat: [('d', 10), ('o', 4)] total: 10 workdays
Nurse1:  d d d o d d d o d d d o d o  day_stat: [('d', 10), ('o', 4)] total: 10 workdays
Nurse2:  d d o d d n o d d d o d n o  day_stat: [('d', 8), ('o', 4), ('n', 2)] total: 10 workdays
Nurse3:  n n o d n o n d n o d d o d  day_stat: [('d', 5), ('o', 4), ('n', 5)] total: 10 workdays
Nurse4:  n o d n n o o d n n o n o o  day_stat: [('d', 2), ('o', 6), ('n', 6)] total: 8 workdays
Nurse5:  o n n d o o o n o o n n o d  day_stat: [('d', 2), ('o', 7), ('n', 5)] total: 7 workdays
Nurse6:  o o n n o o o n o n n d o n  day_stat: [('d', 1), ('o', 7), ('n', 6)] total: 7 workdays

Statistics per day:
        d n o
Day 0:  3 2 2
Day 1:  3 2 2
Day 2:  3 2 2
Day 3:  3 2 2
Day 4:  3 2 2
Day 5:  2 1 4
Day 6:  2 1 4
Day 7:  3 2 2
Day 8:  3 2 2
Day 9:  3 2 2
Day10:  3 2 2
Day11:  3 2 2
Day12:  2 1 4
Day13:  2 1 4

Solving 3 jugs problem using regular

Model: 3_jugs_regular.py
Model: 3_jugs_regular2.py

This is a model solving the classic 3 jugs problem (a.k.a. water jugs problem), based on Taha's network flow model (Introduction to Operations Research, page 245f) but using regular instead. Compare, for example, with the Comet model 3_jugs.co using a network flow approach.

The different possible states in the problem is represented as nodes, and we start at state 1 (node 1). The name of the nodes represent the amount of water in each water jug; e.g. 8,0,0 means that 8 liter/gallon is in the first jug, and 0 in both jug 2 and 3. The goal is then to go from state 8,0,0 to 4,4,0. Here are all the nodes:
nodes = [
    '8,0,0', # 1 start
    '5,0,3', # 2
    '5,3,0', # 3 
    '2,3,3', # 4 
    '2,5,1', # 5
    '7,0,1', # 6
    '7,1,0', # 7
    '4,1,3', # 8
    '3,5,0', # 9
    '3,2,3', # 10
    '6,2,0', # 11
    '6,0,2', # 12
    '1,5,2', # 13
    '1,4,3', # 14
    '4,4,0'  # 15 goal
    ]
The valid transitions between these states are as follows.
states = [
    [2,9],  # state 1
    [3],    # state 2
    [4, 9], # state 3
    [5],    # state 4
    [6,9],  # state 5
    [7],    # state 6
    [8,9],  # state 7
    [15],   # state 8
    [10],   # state 9
    [11],   # state 10
    [12],   # state 11
    [13],   # state 12
    [14],   # state 13
    [15]    # state 14
    ]
And from this adjacency lists we can build the DFA:
transition_fn =  [
   # 1  2  3  4  5  6  7  8  9  0  1  2  3  4  5 
    [0, 2, 0, 0, 0, 0, 0, 0, 9, 0, 0, 0, 0, 0, 0],  # 1
    [0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],  # 2
    [0, 0, 0, 4, 0, 0, 0, 0, 9, 0, 0, 0, 0, 0, 0],  # 3
    [0, 0, 0, 0, 5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],  # 4
    [0, 0, 0, 0, 0, 6, 0, 0, 9, 0, 0, 0, 0, 0, 0],  # 5
    [0, 0, 0, 0, 0, 0, 7, 0, 0, 0, 0, 0, 0, 0, 0],  # 6
    [0, 0, 0, 0, 0, 0, 0, 8, 9, 0, 0, 0, 0, 0, 0],  # 7
    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 15], # 8
    [0, 0, 0, 0, 0, 0, 0, 0, 0, 10, 0, 0, 0, 0, 0], # 9
    [0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 11, 0, 0, 0, 0], # 10
    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 12, 0, 0, 0], # 11
    [0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 13, 0, 0], # 12
    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 14, 0], # 13
    [0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 15], # 14
                                                    # 15
    ]
As you may see this is very much like the representation in the network flow model, but the connections is a state number instead of just 0/1 (connected/not connected).

Running this program gives the following solution where x contains the chosen states, and then the nodes of the states are shown.
x: [1, 9, 10, 11, 12, 13, 14, 15]
8,0,0 -> 3,5,0
3,5,0 -> 3,2,3
3,2,3 -> 6,2,0
6,2,0 -> 6,0,2
6,0,2 -> 1,5,2
1,5,2 -> 1,4,3
1,4,3 -> 4,4,0

num_solutions: 1
failures: 1
branches: 2
wall_time: 0 ms

Found a solution of length 8: [1, 9, 10, 11, 12, 13, 14, 15]
Note: This model don't have any explicit Minimization objective to minimize the number of states (which is implied in the problem). Instead it calls the main function with the length of the array x, ranging from length 1 to 15, and thus find the first, minimal, solution of length 8.

3 jugs problem, alternative solution

Model: 3_jugs_regular2.py

However, I was not satisfied with this calling of main repeatedly. Instead an alternative solution was created which use an extra (last) state which may be looped into (15->15->15...):
# ....
[15],   # state 14
[15]    # state 15: loop ("fix-point")
With this state we can fix the length of the array of decision variables (x) to 15, the number of states, and add an objective variable cost which is then minimized. This variable counts the number of elements in x that is not in the final state (state 15):
cost = solver.IntVar(0, input_max, 'cost')
# ...
sum_b = [solver.IsLessCstVar(x[i], input_max-1) for i in range(input_max)]
solver.Add(cost == solver.Sum(sum_b))

objective = solver.Minimize(cost, 1)

Global contiguity

Model: contiguity_regular.py

The definition of the global constraint global_contiguity from Global Constraint Catalog is:
Enforce all variables of the VARIABLES collection to be assigned to 0 or 1. In addition, all variables assigned to value 1 appear contiguously.

Example:
(<0, 1, 1, 0>)

The global_contiguity constraint holds since the sequence 0 1 1 0 contains no more than one group of contiguous 1.
The regular expression for this is ^0*1*0*$ (or maybe ^0*1+0*$ if there is an assumption that there must be at least one 1). The transition function is thus:
transition_fn =  [
    [1,2], # state 1: 0*
    [3,2], # state 2: 1* 
    [3,0], # state 3: 0* 
    ]
For an array of length 4 there are 11 solutions:
[0, 0, 0, 0]
[0, 0, 0, 1]
[0, 0, 1, 0]
[0, 0, 1, 1]
[0, 1, 0, 0]
[0, 1, 1, 0]
[0, 1, 1, 1]
[1, 0, 0, 0]
[1, 1, 0, 0]
[1, 1, 1, 0]
[1, 1, 1, 1]

num_solutions: 11
failures: 0
branches: 20
wall_time: 0 ms
For an array of length 100 there are 5051 solutions, still with 0 failures (and 10100 branches). It takes about 4 seconds to generate and print all solutions; without printing it takes 0.2 seconds.

Regular expressions

Model: regexp.py

The history behind this is that my last name (Kjellerstrand) is quite often misspelled, in ways according to this regular expression:
k(je|ä)ll(er|ar)?(st|b)r?an?d
Note: In Swedish the two constructs "kje(ll)" and "kä(ll)" sounds alike.

This model generates all the words that can be construed by this regular expression. As with the other models, we create a DFA. Note that the states are groups of characters (maybe it should be just single characters):
transition_fn =  [
   # 1 2 3 4 5 6 7 8 9 0 1 2     # 
    [0,2,3,0,0,0,0,0,0,0,0,0],   #  1 k 
    [0,0,0,4,0,0,0,0,0,0,0,0],   #  2 je
    [0,0,0,4,0,0,0,0,0,0,0,0],   #  3 ä
    [0,0,0,0,5,6,7,8,0,0,0,0],   #  4 ll
    [0,0,0,0,0,0,7,8,0,0,0,0],   #  5 er
    [0,0,0,0,0,0,7,8,0,0,0,0],   #  6 ar
    [0,0,0,0,0,0,0,0,9,10,0,0],  #  7 st 
    [0,0,0,0,0,0,0,0,9,10,0,0],  #  8 b
    [0,0,0,0,0,0,0,0,0,10,0,0],  #  9 r
    [0,0,0,0,0,0,0,0,0,0,11,12], # 10 a
    [0,0,0,0,0,0,0,0,0,0,0,12],  # 11 n
                                 # 12 d 
    ]

s = ['k','je','ä','ll','er','ar','st','b','r','a','n','d']
Running the program generates all the 48 variants of my last name. Though, I have to confess that I have never been called "Kjellbad", "Källbad" or "Kjellbrad". Yet.
kjellstad
kjellbad
källstad
källbad
kjellerstad
kjellerbad
kjellarstad
kjellarbad
kjellstrad
kjellstand
kjellbrad
kjellband
källerstad
källerbad
källarstad
källarbad
källstrad
källstand
källbrad
källband
kjellerstrad
kjellerstand
kjellerbrad
kjellerband
kjellarstrad
kjellarstand
kjellarbrad
kjellarband
kjellstrand
kjellbrand
källerstrad
källerstand
källerbrad
källerband
källarstrad
källarstand
källarbrad
källarband
källstrand
källbrand
kjellerstrand
kjellerbrand
kjellarstrand
kjellarbrand
källerstrand
källerbrand
källarstrand
källarbrand
Compare with the Gecode model all_regexp.cpp, and was described in Regular expressions in Gecode. Here I also describe my fascination with regular expressions.

October 27, 2010

Some new Google CP Solver/Python models

In A first look at Google CP Solver/Python (Google or-tools) and Improvements of some Google CP Solver models I presented my first Google CP Solver models. However, some of my traditional learning problems was not implemented then since I didn't know how to implement certain things with the solver.

Here are implementations of some of these omitted models as well as some other new models. They are is listed below, in all about 40.

Furniture moving (cumulative)

Model: furniture_moving.py

The reason I didn't implement this problem before was that there is no built-in cumulative constraint. This model includes an experimental version of this constraint based on the MiniZinc implementation:
# Parameters:
#
# s: start_times    assumption: array of IntVar
# d: durations      assumption: array of int
# r: resources      assumption: array of int
# b: resource limit assumption: IntVar or int
#
def my_cumulative(solver, s, d, r, b):
    
    # tasks = [i for i in range(len(s))]
    tasks = [i for i in range(len(s)) if r[i] > 0 and d[i] > 0]
    times_min = min([s[i].Min() for i in tasks])
    times_max = max([s[i].Max() + max(d) for i in tasks])
    for t in xrange(times_min, times_max+1):
        bb = []
        for i in tasks:
            c1 = solver.IsLessOrEqualCstVar(s[i], t)  # s[i] <= t
            c2 = solver.IsGreaterCstVar(s[i]+d[i], t) # t < s[i] + d[i]              
            bb.append(c1*c2*r[i])
        solver.Add(solver.Sum(bb) <= b)

    # Somewhat experimental:
    # This constraint is needed to constrain the upper limit of b.
    if not isinstance(b, int):
        solver.Add(b <= sum(r))       
As noted above, this is experimental but seems to solve the problem at hand.

Who killed Agatha: Element and reifications

Model: who_killed_agatha.py

This was another problem I postponed of the standard set of problems, much because I thought it should be hard (or tedious) to convert the many element constraints and reifications to Google CP Solver. However, this was not really true: even if it's not as easy as in MiniZinc or Comet it was quite straightforward.

The problem formulation is from The h1 Tool Suite (originally from F. J. Pelletier: Seventy-five problems for testing automatic theorem provers, Journal of Automated Reasoning, 2: 216, 1986):
Someone in Dreadsbury Mansion killed Aunt Agatha. Agatha, the butler, and Charles live in Dreadsbury Mansion, and are the only ones to live there. A killer always hates, and is no richer than his victim. Charles hates noone that Agatha hates. Agatha hates everybody except the butler. The butler hates everyone not richer than Aunt Agatha. The butler hates everyone whom Agatha hates. Noone hates everyone. Who killed Agatha?
Some examples:

The constraint "A killer always hates, and is no richer than his victim." is translated to this constraint in MiniZinc and Comet: hates[the_killer, the_victim] == 1. In Google CP Solver we state this with an Element constraint:
solver.Add(solver.Element(hates_flat,the_killer*n+the_victim) == 1)
There are also quite a few reifications and logical operators in this model. As I wrote in Improvements of some Google CP Solver models, there are translated straight forward into calls like this (for the constraint Charles hates noone that Agatha hates. :
#  (hates[agatha, i] = 1) => (hates[charles, i] = 0),
for i in range(n):
    b1a = solver.IsEqualCstVar(hates[agatha][i], 1)
    b1b = solver.IsEqualCstVar(hates[charles][i], 0)
    solver.Add(b1a-b1b <= 0)
(We don't have to use Element here since the index i is of type int.)

The actual constraint to the solver store (b1a-b1b<=0) is standard OR translations of the condition (see below).

Another model with a lot of Element is a_round_of_golf.py.

Mr Smith problem (logical operators)

Model: mr_smith.py

This problem is a logic puzzle taken from an IF Prolog example:
The Smith family and their three children want to pay a visit but they do not all have the time to do so. Following are few hints who will go and who will not: - If Mr Smith comes, his wife will come too. - At least one of their two sons Matt and John will come. - Either Mrs Smith or Tim will come, but not both. - Either Tim and John will come, or neither will come. - If Matt comes, then John and his father will also come.
As mentioned earlier, Google CP solver don't have a syntactic sugar for boolean (logic) operations on decision variables, such as AND, OR, IF, EQ. Instead one have to rely on traditional translations ("tricks") used in IP modeling. For example the constraint If Mr Smith comes, his wife will come too (which is stated as Mr_Smith -&gr; Mrs_Smith in other systems), is translated here into:
solver.Add(Mr_Smith - Mrs_Smith <= 0)
Here is a list of some other similar translations. a and b are two BoolVar/IntVar variables (with domain 0..1).
a AND b: a * b == 1
a OR b: a + b >= 1
a XOR b: a + b == 1
IF a THEN b: a - b <= 0
IF AND ONLY IF a THEN b: a == b (i.e. A EQ B)
There are other ways of formulating these with may be more natural for the problem at hand. Also, for more complex constraints one may have to tweak a little more. For more about this, see the excellent book Model Building in Mathematical Programming by H. Paul Williams (ISBN: 9780471997887), especially chapter 9.

The new models

Here are all the new models, also available via my Google CP Solver page. They are also at the Google Code SVN repository.

October 05, 2010

Improvements of some Google CP Solver models

In A first look at Google CP Solver/Python (Google or-tools) I showed my first models in Google CP Solver (or-tools).

Laurent Perron (the designer of the Google CP Solver) was kind enough to comment and improve some of these models. Thanks, Laurent!

Decomposition of alldifferent_except_0

Model: alldifferent_except_0.py

Laurent mentioned that reifications can be stated much simpler than my approach. My first version was:
def alldifferent_except_0(solver, a):
    n = len(a)
    for i in range(n):
        for j in range(i):
            # This syntactic sugar would be nice to have:
            # solver.Add( ((a[i] != 0) and (a[j] != 0)) -> (a[i] != a[j] ) )            
            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.IsDifferentCt(a[i],a[j], bij))
            solver.Add(bi*bj <= bij)
However, these temporary variables don't have to be declared and added to the solver. Instead one can write
bi = solver.IsDifferentCstVar(a[i],0)
Note the slightly different method name. The Var at the end means that the method returns a variable.

So, the decomposition is now much nicer and easier to understand.
def alldifferent_except_0(solver, a):
   n = len(a)
   for i in range(n):
      for j in range(i):
          bi = solver.IsDifferentCstVar(a[i],0)
          bj = solver.IsDifferentCstVar(a[j],0)
          bij = solver.IsDifferentVar(a[i],a[j])
          solver.Add(bi*bj <= bij)
This following one liner version works but is less readable. (Yes, I'm trying to push the limits. :-)
# ...
solver.Add(solver.IsDifferentCstVar(a[i],0)*
           solver.IsDifferentCstVar(a[j],0) <=
           solver.IsDifferentVar(a[i],a[j]))

Young tableaux and partition

Model: young_tableaux.py

The improvement here was about the same as for alldifferent_except_0. The following constraints:
        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))
is replaced with
         b = [solver.IsLessOrEqualCstVar(x[(i, j)], n) for j in range(n)]
         solver.Add(p[i] == solver.Sum(b))
Again, there are not so many variables sent to the constraint store. Alternative it could be written in a more compact form:
solver.Add(p[i] == solver.Sum([solver.IsLessOrEqualCstVar(x[(i, j)], n)
                               for j in range(n)]))

Coins grid

Model: coins_grid.py

Laurent suggested several improvements in this model:
  • First, I had forgot to add objective to the monitor list.
  • The objective (z) is not a (declared) decision variable any more. It was defined as
    z = solver.IntVar(0,1000000, 'z')
    # ...
    solver.Add(solver.Sum([x[(i,j)]*abs(i-j)*abs(i-j) 
                          for i in range(n) for j in range(n)]) == z)
    # ...
    
    Now the return value of solver.Sum is used. The sum is also slightly redefined (skipping two abs calls):
    objective_var = solver.Sum([x[(i, j)] * (i - j) * (i - j)
                               for i in range(n) for j in range(n)])
    
  • Added solution.AddObjective(objective_var)
  • Using the value strategy solver.ASSIGN_MAX_VALUE instead of solver.ASSIGN_MIN_VALUE.
  • The rows/column summation use the faster solver.SumEquality instead of solver.Sum

Hidato

My own Hidato model was quite ugly and slow.

Laurent wrote a much more elegant version (hidato_laurent.py) using, amongst other things, the global constraint solver.AllowedAssignments. This solves the hardest problem in the file in a blink with just a few failures.

A humbling lesson.

New model: A "general" alphametic solver

Model: alphametic.py

Yesterday I wrote a "general" alphametic solver which solves used defined (from the command-line) problems of the form "SEND+MORE=MONEY". It only solves addition problems of the form NUMBER_1+NUMBER_2...+NUMBER_N_1 = NUMBER_N where the last number is the sum (so it's not very general).

The program has the following usage:
python alphametic.py
         ->  solves SEND+MORE=MONEY in base 10
    
python alphametic.py  'SEND+MOST=MONEY' 11
         -> solver SEND+MOST=MONEY in base 11

python alphametic.py TEST 
         -> solve some test problems in base 
                   (defined in test_problems())
Compare with the Zinc model alphametic.zinc which use the same basic ideas. However Zinc don't have Python's excellent support for strings and sets.

October 03, 2010

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

Here is my Google CP Solver page.

Update: With help from Laurent Perron, I have improved some of these models. See Improvements of some Google CP Solver models

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()
    solution.Add(x_flat)
    solution.Add(s)
    db = solver.Phase(x_flat,
                      solver.CHOOSE_FIRST_UNBOUND,                     
                      solver.ASSIGN_CENTER_VALUE
                      )    
    solver.NewSearch(db)

    # 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(),
            print
        print
        num_solutions += 1        
    solver.EndSearch()

    print
    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],
                            solver.CHOOSE_PATH,
                            solver.ASSIGN_MIN_VALUE),
                            [objective, search_log, collector])

Diet1

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()
  solution.Add(cost)
  solution.Add(x)

  # last solution since it's a minimization problem
  collector = solver.LastSolutionCollector(solution)
  search_log = solver.SearchLog(100, cost)
  solver.Solve(solver.Phase(x + [cost],
                            solver.INT_VAR_SIMPLE,
                            solver.ASSIGN_MIN_VALUE),
                            [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)],
                           solver.CHOOSE_PATH,
                           solver.ASSIGN_MIN_VALUE),
                           [collector])
 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)]), 
         print
     print


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)],
                          solver.CHOOSE_PATH,
                          solver.ASSIGN_MIN_VALUE)
solver.NewSearch(db)

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(),
        print
    print
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)])
solution.Add(total_sum)

solver.Solve(solver.Phase([x[(i,j)] for i in range(n) for j in range(n)],
                          solver.CHOOSE_PATH,
                          solver.ASSIGN_MIN_VALUE),
                          [collector])
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()
solution.Add(x_flat)

# etc ...

coins_grid

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.

toNum

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
20.03016
30.041108
40.070384
50.2311000/td>
60.7362160/td>
72.24116/td>
867168/td>
91611664/td>
104218000/td>

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)

nqueens

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).
solver.Add(solver.AllDifferent(q,True))
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.

Minesweeper

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 = [
    [2,3,X,2,2,X,2,1],
    [X,X,4,X,X,4,X,2],
    [X,X,X,X,X,X,4,X],
    [X,5,X,6,X,X,X,2],
    [2,X,X,X,5,5,X,2],
    [1,3,4,X,X,X,4,X],
    [0,1,X,4,X,X,X,3],
    [0,1,2,X,2,3,X,2]
    ]

# ....
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
           solver.Add(
               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))

xkcd

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)
        cargs.extend([objective])

    # ....

    solver.Solve(solver.Phase(x,
                              solver.CHOOSE_FIRST_UNBOUND,
                              solver.ASSIGN_MIN_VALUE),
                              cargs)
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.

Nonogram

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
    solver.Add(
        solver.Element(A_flat,E[overlapping[I][0]]*word_len+overlapping[I][1])
        == 
        solver.Element(A_flat,E[overlapping[I][2]]*word_len+overlapping[I][3]))

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)))
    x.append(t)
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)]

Discrete_tomography

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:

Sicherman_dice

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.

Hidato

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)