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
p2161.1

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.


# 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)

# 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.


# 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)

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).

The Gurobi Solver

Example 1 – Love Hearts

Example 2 – Quad Spoiling

Recent Posts

Recent Comments

Archives

Categories

Meta