The Gurobi Solver
Gurobi offer a high performance solver with Linear Programming capabilities and a nice Python interface here. This page contains a number of examples of how it can be used to solve some of the puzzles published by the Sunday Times and the New Scientist in the UK.
Example 1 – Love Hearts
This puzzle was published in the Sunday Times in 2004 and is notable because the answer given was incorrect. The puzzle is repeated here for convenience:
Sunday Times Teaser Number 2161 – Love Hearts
by Hugh Bradley and Christopher Higgens
Bobby has made a Valentine card for his girlfriend Chin-We. He drew rows of identical touching circles in a triangular formation as shown below but with more rows.
He then placed a red heart in as many circles as possible without any of these forming equilateral triangles. He continued by placing pink hearts in some of the other circles and white circles in the remainder.
When he counted the number of hearts of different colours (in no particular order) he obtained three consecutive integers.
How many red hearts are on the card?
The following program solves this problem for n rows. It creates a binary variable for each circle and then seeks to maximise the number of such circles included subject to constraints that prevent the included circles forming any equilateral triangles.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 |
# Sunday Times Teaser 2161 -- Love Hearts # # by Hugh Bradley and Christopher Higgens # # O Bobby has made a Valentine card for his girlfriend # O O He has drawn rows of identical touching circles in # O O O a triangular formation (as shown on the left) with # O O O O the successive rows having 1, 2, 3 ... circles. He # has then placed a red heart in as many circles as # possible without any of them forming equilateral triangles. # He then put pink hearts in some of the remaining circles and # white hearts in the rest. When he counted the number of # hearts of different colours he obtained three consecutive # numbers. # # How many red hearts are on the card? # from __future__ import print_function import argparse from gurobipy import * # the circles are numbered as follows: # # 1 (1,1) # 2 3 (2,1) (2,2) # 4 5 6 (3,1) (3,2) (3,3) # .... def solve(r, solve=True, export=False, out_path=None): # intialise the model z = Model() # set to optimise for a maximum z.ModelSense = GRB.MAXIMIZE # tune the model z.setParam(GRB.Param.Symmetry, 1) # the dictionary of decision variables, one variable # for each circle with i in (1 .. r) as the row and # j in (1 .. i) as the position within the row a = {(i, j): z.addVar(vtype = GRB.BINARY, obj = 1.0, name = 'a[{:d}, {:d}]'.format(i, j)) for i in range(1, r + 1) for j in range(1, i + 1)} z.update() # the constraints - enumerate all equilateral triangles # and prevent any such triangles being formed by keeping # the number of included circles at its vertexes below 3 # for each row except the last for i in range(1, r): # for each position in this row for j in range(1, i + 1): # for each triangle of side length (k) with its upper vertex at # (i, j) and its sides parallel to those of the overall shape for k in range(1, r - i + 1): # the sets of 3 points at the same distances clockwise along the # sides of these triangles form k equilateral triangles for m in range(k): u, v, w = (i + m, j), (i + k, j + m), (i + k - m, j + k - m) z.addConstr(a[u] + a[v] + a[w] <= 2, name='a[{:s},{:s},{:s}]'.format(u, v, w)) if export and out_path: z.update() z.write(out_path + 'st2161.' + str(r) + '.mps') if solve: z.optimize() if out_path: z.write(out_path + 'st2161.' + str(r) + '.sol') path = 'C:\\Users\\Brian\\Documents\\Python\\' parser = argparse.ArgumentParser(prog='puzzle', description='Teaser 2161') parser.add_argument('-r', dest='rows', required=True, type=int, help='Puzzle size') args = parser.parse_args() solve(args.rows, out_path=path) |
The performance of the Gurobi solver allows us to find the maximum number of red hearts we can add without forming equilateral triangles for up to 14 rows.
rows | 9 | 10 | 11 | 12 | 13 | 14 |
red hearts | 17 | 20 | 22 | 25 | 28 | 31 |
variables | 45 | 55 | 66 | 78 | 91 | 105 |
constraints | 330 | 495 | 715 | 1001 | 1365 | 1820 |
The ‘consecutive integer’ condition on the three colours of hearts limits which of these solutions solves the teaser and gives the answer as 22 red hearts. The ‘official’ answer given by the Sunday Times was 16!
Example 2 – Quad Spoiling
This example was published in the New Scientist edition number 1225 on 30th October 1980. It is covered in some detail here.
New Scientist Enigma Number 82 – Quad Spoiling
by Stephen Ainley
There are a lot of quadrilaterals in the figure, 492 in fact, I reckon. I do not count ‘crossed quadrilaterals’, like ABCD, only simple ones, like ABED, ABFD, and so on. Now imagine the figure composed of 63 matches. My question is: what is the smallest number of matches you need to remove so as to spoil all 492 quadrilaterals? A quadrilateral is spoilt, of course, when one or more matches is missing from its boundary.
Here is a solution using the Python interface to the Gurobi solver. It creates a binary variable for each individual match, enumerates every quadrilateral and establishes a list of the matches on the boundary of each of them. The aim is to minimise the number variables that are true (missing matches) whilst ensuring that the sum of the variables for each quadrilateral is greater than zero (every quadrilateral has at least one missing match). This type of problem is often referred to as a Minimum Hitting Sum.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 |
# New Scientist Enigma Number 82 # # From New Scientist #1225, 30th October 19#80 # # 30th October 1980 by Stephen Ainley # # Quad-Spoiling # # /\ # / \ # A +----+ B # /\ / \ # / \ / \ # +----+-----+ # /\ / \ / \ # / \ / \ / \ # +----+-----+-----+ # /\ / \ / \ / \ # / \ / \ / \ / \ # +----+-----+-----+-----+ # /\ / \ / \ / \ / \ # / \ / \ / \ / \ / \ # +----+-----+-----+-----+-----+ E # /\ / \ / \ / \ / \ / \ # / \ / \ / \ / \ / \ / \ # +----+-----+-----+-----+-----+-----+ F # C D # # There are a lot of quadrilaterals in the figure, 492 in fact, I reckon. # I do not count 'crossed quadrilaterals', like ABCD, only simple ones, # like ABED, ABFD, and so on. # # Now imagine the figure composed of 63 matches. My question is: what is # the smallest number of matches you need to remove so as to spoil all # 492 quadrilaterals? A quadrilateral is spoilt, of course, when one or # more matches is missing from its boundary. from __future__ import print_function from itertools import combinations import argparse from gurobipy import * # enumerate the vertexes of the small triangles as: # # (0,0) 1 # (1,0) (1,1) 2 3 # (2,0) (2,1) (2,2) 4 5 6 # (3,0) (3,1) (3,2) (3,3) 7 8 9 10 # .... # the index for the vertex at row r, position in row c def ix(r, c): return r * (r + 1) // 2 + c + 1 # all the vertex points in a formation with 'r' rows def points(r): pts = [] for i in range(r + 1): for j in range(i + 1): pts.append((i, j)) return pts # four points are degenerate if three or more are on the # same line and cannot hence form a quadrilateral def degenerate(pts): for a, b, c in combinations(pts, 3): if (a[0] == b[0] == c[0] or a[1] == b[1] == c[1] or a[0] - a[1] == b[0] - b[1] == c[0] - c[1]): return True return False # check if two points are on a line that is aligned # with one of the three sides of the formation def aligned(p1, p2): return (p1[0] == p2[0] or p1[1] == p2[1] or p1[0] - p1[1] == p2[0] - p2[1]) # find all the quadrilaterals in the formation def quads(rows): # find all the vertex points p = points(rows) # now combine them four at a time # (note here that a < b < c < d) for t in combinations(p, 4): # check that they are not degenerate if not degenerate(t): a, b, c, d = t # ensure b is right of a if a[0] < b[0] and a[1] == b[1]: a, b = b, a # ensure d is right of c if c[0] < d[0] and c[1] == d[1]: c, d = d, c # now check that each pair of points lies # on a line parallel to one of the sides if (aligned(a, b) and aligned(a, c) and aligned(b, d) and aligned(c, d)): yield (a, b, d, c) # find the set of individual segments on the line between # two vertex points with segments identified by the indexes # of the two vertexes that they join def line_segs(a, b): if a[0] == b[0]: # they share a line parallel to the base return set(tuple(sorted((ix(a[0], i), ix(a[0], i + 1)))) for i in range(min(a[1], b[1]), max(a[1], b[1]))) elif a[1] == b[1]: # they share a line parallel to the left hand side return set(tuple(sorted((ix(i, a[1]), ix(i + 1, a[1])))) for i in range(min(a[0], b[0]), max(a[0], b[0]))) else: # they share a line parallel to the right hand side t = a[1] - a[0] return set(tuple(sorted((ix(i, t + i), ix(i + 1, t + i + 1)))) for i in range(min(a[0], b[0]), max(a[0], b[0]))) # list the segments on the perimeter of a quadrilateral def quad_segs(q): return ( line_segs(q[0], q[1]) | line_segs(q[1], q[2]) | line_segs(q[2], q[3]) | line_segs(q[3], q[0]) ) # solve for a given number of rows def solve(r, solve=True, export=False, out_path=None): # list the segments on all quadrilaterals within # a formation with this number of rows constraints = [quad_segs(q) for q in quads(r)] # compile a list of all segments all_segs = set() for s in constraints: all_segs |= s segs = tuple(sorted(all_segs)) # now map segment identifiers onto integers 1 .. size size = len(all_segs) inx = {s:(i + 1) for i, s in enumerate(segs)} xni = {(i + 1):str(s) for i, s in enumerate(segs)} # intialise the model z = Model() # set to optimise for a minimum z.ModelSense = GRB.MINIMIZE # tune the model z.setParam(GRB.Param.Symmetry, -1) # the dictionary of decision variables, one variable # for each circle with i in (0 .. r) as the row and # j in (0 .. i + 1) as the position within the row a = {(i, j): z.addVar(vtype = GRB.BINARY, obj = 1.0, name = 'a[{:d},{:d}]'.format(i, j)) for i, j in segs} z.update() # now add the constraints # for each quadrilateral c contains the segment identities for n, c in enumerate(constraints): # construct a linear expression constr = LinExpr() # for each segment for i, j in c: constr.add(a[(i, j)]) z.addConstr(constr >= 1, name=str(n + 1)) # export model in mps format if requested if export and out_path: z.update() z.write(out_path + 'ns0082.' + str(r) + '.mps') # solve the model and output the solution if requested if solve: z.optimize() if out_path: z.write(out_path + 'ns0082.' + str(r) + '.sol') out_path = 'C:\\Users\\Brian\\Documents\\lpsolve\\' parser = argparse.ArgumentParser(prog='puzzle', description='Teaser 2161') parser.add_argument('-r', dest='rows', required=True, type=int, help='Puzzle size') args = parser.parse_args() solve(args.rows) |
Here are the results for varying numbers of rows, the answer to the enigma itself being that 20 matches need to be removed.
rows | 4 | 5 | 6 | 7 | 8 | 9 |
matches | 9 | 14 | 20 | 27 | 35 | 44 |
variables | 30 | 45 | 63 | 84 | 108 | 135 |
constraints | 102 | 243 | 492 | 894 | 1500 | 3570 |
Jim Randell poostulates here that the number of matches needed to spoil all quadrilaterals follows A000096 at the The On-Line Encyclopedia of Integer Sequences (OEIS).