From 9026145d65fbc87e685dd51cdc2e1bb88d106fd1 Mon Sep 17 00:00:00 2001 From: hakank Date: Sun, 10 Oct 2010 07:17:39 +0000 Subject: [PATCH] Added models crew.py all_interval.py divisible_by_9_through_1.py coins3.py. Added a note in hidato.py to use Laurent's hidato_table.py instead. --- python/all_interval.py | 113 +++++++++++++++ python/coins3.py | 104 ++++++++++++++ python/crew.py | 224 +++++++++++++++++++++++++++++ python/divisible_by_9_through_1.py | 178 +++++++++++++++++++++++ python/hidato.py | 3 +- 5 files changed, 621 insertions(+), 1 deletion(-) create mode 100644 python/all_interval.py create mode 100644 python/coins3.py create mode 100644 python/crew.py create mode 100644 python/divisible_by_9_through_1.py diff --git a/python/all_interval.py b/python/all_interval.py new file mode 100644 index 0000000000..61f64cd574 --- /dev/null +++ b/python/all_interval.py @@ -0,0 +1,113 @@ +# Copyright 2010 Hakan Kjellerstrand hakank@bonetmail.com +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +""" + + All interval problem in Google CP Solver. + + CSPLib problem number 7 + http://www.cs.st-andrews.ac.uk/~ianm/CSPLib/prob/prob007/index.html + ''' + Given the twelve standard pitch-classes (c, c , d, ...), represented by + numbers 0,1,...,11, find a series in which each pitch-class occurs exactly + once and in which the musical intervals between neighbouring notes cover + the full set of intervals from the minor second (1 semitone) to the major + seventh (11 semitones). That is, for each of the intervals, there is a + pair of neigbhouring pitch-classes in the series, between which this + interval appears. The problem of finding such a series can be easily + formulated as an instance of a more general arithmetic problem on Z_n, + the set of integer residues modulo n. Given n in N, find a vector + s = (s_1, ..., s_n), such that (i) s is a permutation of + Z_n = {0,1,...,n-1}; and (ii) the interval vector + v = (|s_2-s_1|, |s_3-s_2|, ... |s_n-s_{n-1}|) is a permutation of + Z_n-{0} = {1,2,...,n-1}. A vector v satisfying these conditions is + called an all-interval series of size n; the problem of finding such + a series is the all-interval series problem of size n. We may also be + interested in finding all possible series of a given size. + ''' + + Compare with the following models: + * MiniZinc: http://www.hakank.org/minizinc/all_interval.mzn + * Comet : http://www.hakank.org/comet/all_interval.co + * Gecode/R: http://www.hakank.org/gecode_r/all_interval.rb + * ECLiPSe : http://www.hakank.org/eclipse/all_interval.ecl + * SICStus : http://www.hakank.org/sicstus/all_interval.pl + + + This model was created by Hakan Kjellerstrand (hakank@bonetmail.com) + Also see my other Google CP Solver models: http://www.hakank.org/google_or_tools/ + +""" +import string, sys + +from constraint_solver import pywrapcp + +def main(n=12): + + # Create the solver. + solver = pywrapcp.Solver('All interval') + + # + # data + # + print "n:", n + + # + # declare variables + # + x = [solver.IntVar(1, n, 'x[%i]' % i) for i in range(n)] + diffs = [solver.IntVar(1, n-1, 'diffs[%i]' % i) for i in range(n-1)] + + # + # constraints + # + solver.Add(solver.AllDifferent(x, True)) + solver.Add(solver.AllDifferent(diffs, True)) + + for k in range(n-1): + solver.Add(diffs[k] == abs(x[k+1]-x[k])) + + # symmetry breaking + solver.Add(x[0] < x[n-1]) + solver.Add(diffs[0] < diffs[1]) + + # + # solution and search + # + solution = solver.Assignment() + solution.Add(x) + solution.Add(diffs) + + db = solver.Phase(x, + solver.CHOOSE_FIRST_UNBOUND, + solver.ASSIGN_MIN_VALUE) + + solver.NewSearch(db) + num_solutions = 0 + while solver.NextSolution(): + print "x:", [x[i].Value() for i in range(n)] + print "diffs:", [diffs[i].Value() for i in range(n-1)] + num_solutions += 1 + print + + print "num_solutions:", num_solutions + print "failures:", solver.failures() + print "branches:", solver.branches() + print "wall_time:", solver.wall_time() + +n=12 +if __name__ == '__main__': + if len(sys.argv) > 1: + n = string.atoi(sys.argv[1]) + main(n) diff --git a/python/coins3.py b/python/coins3.py new file mode 100644 index 0000000000..a726f10023 --- /dev/null +++ b/python/coins3.py @@ -0,0 +1,104 @@ +# Copyright 2010 Hakan Kjellerstrand hakank@bonetmail.com +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +""" + + Coin application in Google CP Solver. + + From 'Constraint Logic Programming using ECLiPSe' + pages 99f and 234 ff. + The solution in ECLiPSe is at page 236. + + ''' + What is the minimum number of coins that allows one to pay _exactly_ + any amount smaller than one Euro? Recall that there are six different + euro cents, of denomination 1, 2, 5, 10, 20, 50 + ''' + + Compare with the following models: + * MiniZinc: http://hakank.org/minizinc/coins3.mzn + * Comet : http://www.hakank.org/comet/coins3.co + * Gecode : http://hakank.org/gecode/coins3.cpp + * SICStus : http://hakank.org/sicstus/coins3.pl + + + This model was created by Hakan Kjellerstrand (hakank@bonetmail.com) + Also see my other Google CP Solver models: http://www.hakank.org/google_or_tools/ +""" +import sys, string +from constraint_solver import pywrapcp + + +def main(): + # Create the solver. + solver = pywrapcp.Solver('Coins') + + # + # data + # + n = 6 # number of different coins + variables = [1, 2, 5, 10, 25, 50] + + + # declare variables + x = [solver.IntVar(0, 99, 'x%i' % i) for i in range(n)] + num_coins = solver.IntVar(0, 99, 'num_coins') + + # + # constraints + # + + # number of used coins, to be minimized + solver.Add(num_coins == solver.Sum(x)) + + # Check that all changes from 1 to 99 can be made. + for j in range(1, 100): + tmp = [solver.IntVar(0, 99, 'b%i'%i) for i in range(n)] + solver.Add(solver.ScalProd(tmp, variables) == j) + [solver.Add(tmp[i] <= x[i]) for i in range(n)] + + + # objective + objective = solver.Minimize(num_coins, 1) + + # + # solution and search + # + solution = solver.Assignment() + solution.Add(x) + solution.Add(num_coins) + solution.AddObjective(num_coins) + + db = solver.Phase(x, + solver.CHOOSE_MIN_SIZE_LOWEST_MAX, + solver.ASSIGN_MIN_VALUE) + + solver.NewSearch(db, [objective]) + num_solutions = 0 + while solver.NextSolution(): + print "x: ", [x[i].Value() for i in range(n)] + print "num_coins:", num_coins.Value() + 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() + + +if __name__ == '__main__': + main() diff --git a/python/crew.py b/python/crew.py new file mode 100644 index 0000000000..e3f7cbfcc3 --- /dev/null +++ b/python/crew.py @@ -0,0 +1,224 @@ +# Copyright 2010 Hakan Kjellerstrand hakank@bonetmail.com +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +""" + + Crew allocation problem in Google CP Solver. + + From Gecode example crew + examples/crew.cc + ''' + * Example: Airline crew allocation + * + * Assign 20 flight attendants to 10 flights. Each flight needs a certain + * number of cabin crew, and they have to speak certain languages. + * Every cabin crew member has two flights off after an attended flight. + * + ''' + + Compare with the following models: + * MiniZinc: http://www.hakank.org/minizinc/crew.mzn + * Comet : http://www.hakank.org/comet/crew.co + * ECLiPSe : http://hakank.org/eclipse/crew.ecl + * SICStus : http://hakank.org/sicstus/crew.pl + + + This model was created by Hakan Kjellerstrand (hakank@bonetmail.com) + Also see my other Google CP Solver models: http://www.hakank.org/google_or_tools/ +""" +import sys, string +from constraint_solver import pywrapcp + +def main(sols=1): + + # Create the solver. + solver = pywrapcp.Solver('Crew') + + # + # data + # + names = ["Tom", + "David", + "Jeremy", + "Ron", + "Joe", + "Bill", + "Fred", + "Bob", + "Mario", + "Ed", + "Carol", + "Janet", + "Tracy", + "Marilyn", + "Carolyn", + "Cathy", + "Inez", + "Jean", + "Heather", + "Juliet" + ] + + num_persons = len(names) # number of persons + + attributes = [ + # steward, hostess, french, spanish, german + [1,0,0,0,1], # Tom = 1 + [1,0,0,0,0], # David = 2 + [1,0,0,0,1], # Jeremy = 3 + [1,0,0,0,0], # Ron = 4 + [1,0,0,1,0], # Joe = 5 + [1,0,1,1,0], # Bill = 6 + [1,0,0,1,0], # Fred = 7 + [1,0,0,0,0], # Bob = 8 + [1,0,0,1,1], # Mario = 9 + [1,0,0,0,0], # Ed = 10 + [0,1,0,0,0], # Carol = 11 + [0,1,0,0,0], # Janet = 12 + [0,1,0,0,0], # Tracy = 13 + [0,1,0,1,1], # Marilyn = 14 + [0,1,0,0,0], # Carolyn = 15 + [0,1,0,0,0], # Cathy = 16 + [0,1,1,1,1], # Inez = 17 + [0,1,1,0,0], # Jean = 18 + [0,1,0,1,1], # Heather = 19 + [0,1,1,0,0] # Juliet = 20 + ] + + # The columns are in the following order: + # staff : Overall number of cabin crew needed + # stewards : How many stewards are required + # hostesses : How many hostesses are required + # french : How many French speaking employees are required + # spanish : How many Spanish speaking employees are required + # german : How many German speaking employees are required + required_crew = [ + [4,1,1,1,1,1], # Flight 1 + [5,1,1,1,1,1], # Flight 2 + [5,1,1,1,1,1], # .. + [6,2,2,1,1,1], + [7,3,3,1,1,1], + [4,1,1,1,1,1], + [5,1,1,1,1,1], + [6,1,1,1,1,1], + [6,2,2,1,1,1], # ... + [7,3,3,1,1,1] # Flight 10 + ] + + num_flights = len(required_crew) # number of flights + + # + # declare variables + # + crew = {} + for i in range(num_flights): + for j in range(num_persons): + crew[(i,j)] = solver.IntVar(0, 1, 'crew[%i,%i]' % (i, j)) + crew_flat = [crew[(i,j)] for i in range(num_flights) for j in range(num_persons)] + + # number of working persons + num_working = solver.IntVar(1, num_persons, 'num_working') + + + # + # constraints + # + + # number of working persons + solver.Add(num_working == solver.Sum( + [solver.IsGreaterOrEqualCstVar(solver.Sum([crew[(f,p)] + for f in range(num_flights)]), 1) + for p in range(num_persons) ]) ) + + for f in range(num_flights): + # size of crew + tmp = [crew[(f,i)] for i in range(num_persons)] + solver.Add(solver.Sum(tmp) == required_crew[f][0]) + + # attributes and requirements + for j in range(5): + tmp = [attributes[i][j]*crew[(f,i)] for i in range(num_persons) ] + solver.Add(solver.Sum(tmp) >= required_crew[f][j+1]) + + + # after a flight, break for at least two flights + for f in range(num_flights-2): + for i in range(num_persons): + solver.Add(crew[f,i] + crew[f+1,i] + crew[f+2,i] <= 1) + + # extra contraint: all must work at least two of the flights + # for i in range(num_persons): + # [solver.Add(solver.Sum([crew[f,i] for f in range(num_flights)]) >= 2) ] + + + # + # solution and search + # + solution = solver.Assignment() + solution.Add(crew_flat) + solution.Add(num_working) + + db = solver.Phase(crew_flat, + solver.CHOOSE_FIRST_UNBOUND, + solver.ASSIGN_MIN_VALUE) + + # + # result + # + solver.NewSearch(db) + num_solutions = 0 + while solver.NextSolution(): + num_solutions += 1 + print "Solution #%i" % num_solutions + print "Number working:", num_working.Value() + for i in range(num_flights): + for j in range(num_persons): + print crew[i,j].Value(), + print + print + + print "Flights:" + for flight in range(num_flights): + print "Flight", flight, "persons:", + for person in range(num_persons): + if crew[flight, person].Value() == 1: + print names[person], + print + print + + print "Crew:" + for person in range(num_persons): + print "%-10s flights" % names[person], + for flight in range(num_flights): + if crew[flight, person].Value() == 1: + print flight, + print + print + + if num_solutions >= sols: + break + solver.EndSearch() + + print + print "num_solutions:", num_solutions + print "failures:", solver.failures() + print "branches:", solver.branches() + print "wall_time:", solver.wall_time() + +num_solutions_to_show = 1 +if __name__ == '__main__': + if (len(sys.argv) > 1): + num_solutions_to_show = string.atoi(sys.argv[1]) + + main(num_solutions_to_show) diff --git a/python/divisible_by_9_through_1.py b/python/divisible_by_9_through_1.py new file mode 100644 index 0000000000..a935ea8d39 --- /dev/null +++ b/python/divisible_by_9_through_1.py @@ -0,0 +1,178 @@ +# Copyright 2010 Hakan Kjellerstrand hakank@bonetmail.com +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +""" + + Divisible by 9 through 1 puzzle in Google CP Solver. + + From http://msdn.microsoft.com/en-us/vcsharp/ee957404.aspx + ' Solving Combinatory Problems with LINQ' + ''' + Find a number consisting of 9 digits in which each of the digits + from 1 to 9 appears only once. This number must also satisfy these + divisibility requirements: + + 1. The number should be divisible by 9. + 2. If the rightmost digit is removed, the remaining number should + be divisible by 8. + 3. If the rightmost digit of the new number is removed, the remaining + number should be divisible by 7. + 4. And so on, until there's only one digit (which will necessarily + be divisible by 1). + ''' + + Also, see + 'Intel Parallel Studio: Great for Serial Code Too (Episode 1)' + http://software.intel.com/en-us/blogs/2009/12/07/intel-parallel-studio-great-for-serial-code-too-episode-1/ + + + This model is however generalized to handle any base, for reasonable limits. + The 'reasonable limit' for this model is that base must be between 2..16. + + Compare with the following models: + * MiniZinc: http://www.hakank.org/minizinc/divisible_by_9_through_1.mzn + * Comet : http://www.hakank.org/comet/divisible_by_9_through_1.co + * ECLiPSe : http://www.hakank.org/eclipse/divisible_by_9_through_1.ecl + * Gecode : http://www.hakank.org/gecode/divisible_by_9_through_1.cpp + + + This model was created by Hakan Kjellerstrand (hakank@bonetmail.com) + Also see my other Google CP Solver models: http://www.hakank.org/google_or_tools/ + +""" +import string, sys + +from constraint_solver import pywrapcp + +# +# Decomposition of modulo constraint +# +# This implementation is based on the ECLiPSe version +# mentioned in +# - A Modulo propagator for ECLiPSE' +# http://www.hakank.org/constraint_programming_blog/2010/05/a_modulo_propagator_for_eclips.html +# The ECLiPSe source code: +# http://www.hakank.org/eclipse/modulo_propagator.ecl +# +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(max(0,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) + + +# +# converts a number (s) <-> an array of integers (t) in the specific base. +# +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)])) + + +def main(base=10): + + # Create the solver. + solver = pywrapcp.Solver('Divisible by 9 through 1') + + # data + m = base**(base-1)-1 + n = base - 1 + + digits_str = "_0123456789ABCDEFGH" + + print "base:", base + + # declare variables + + # the digits + x = [solver.IntVar(1, base-1, 'x[%i]' % i) for i in range(n)] + + # the numbers, t[0] contains the answer + t = [solver.IntVar(0, m, 't[%i]' % i) for i in range(n)] + + # + # constraints + # + solver.Add(solver.AllDifferent(x, True)) + + for i in range(n): + mm = base-i-1 + toNum(solver, [x[j] for j in range(mm)], t[i], base) + my_mod(solver, t[i], mm, 0) + + + # + # solution and search + # + solution = solver.Assignment() + solution.Add(x) + solution.Add(t) + + db = solver.Phase(x, + solver.CHOOSE_FIRST_UNBOUND, + solver.ASSIGN_MIN_VALUE) + + + solver.NewSearch(db) + num_solutions = 0 + while solver.NextSolution(): + print "x: ", [x[i].Value() for i in range(n)] + print "t: ", [t[i].Value() for i in range(n)] + print "number base 10: %i base %i: %s" % (t[0].Value(),\ + base, + "".join([digits_str[x[i].Value()+1] for i in range(n)])) + print + num_solutions += 1 + solver.EndSearch() + + print "num_solutions:", num_solutions + print "failures:", solver.failures() + print "branches:", solver.branches() + print "wall_time:", solver.wall_time() + + +base = 10 +default_base = 10 +max_base = 16 +if __name__ == '__main__': + if len(sys.argv) > 1: + base = string.atoi(sys.argv[1]) + if base > max_base: + print "Sorry, max allowed base is %i. Setting base to %i..." % (max_base, default_base) + base = default_base + main(base) + + # for base in range(2, 17): + # print + # main(base) diff --git a/python/hidato.py b/python/hidato.py index a5dacf2960..de4b8d4998 100644 --- a/python/hidato.py +++ b/python/hidato.py @@ -34,7 +34,8 @@ * ECLiPSe: http://hakank.org/eclipse/hidato.ecl * SICStus: http://hakank.org/sicstus/hidato.pl - + Note: This model is very slow. Please see Laurent Perron's much faster + (and more elegant) model: hidato_table.py . This model was created by Hakan Kjellerstrand (hakank@bonetmail.com) Also see my other Google CP Solver models: http://www.hakank.org/google_or_tools/