Enigmatic Code
Programming Enigma Puzzles
Enigma 144: Spoiling triangles and quads
Posted by on 11 November 2013
From New Scientist #1289, 21st January 1982 [link] [link]
A is a grid of 63 matchsticks. B shows one way of removing 21 matches so as to “spoil” all the triangles in the grid (there are 78 of various sizes) by removing one or more matches from the boundary of each. And C shows a way (an improvement due to Miss Gregory and Mr Buckle on the previously published solution to Enigma 82) of removing 20 matches so as to spoil all the 492 quadrilaterals. The problem now is: what is the smallest number of matches you must remove so as to spoil all the triangles and all the quadrilaterals in A?
Note that C actually has 21 matches removed, although it is possible to solve Enigma 82 by removing only 20 matches (see the diagram attached to my solution).
[enigma144]
Enigmatic Code 
As noted this problem is an extension of Enigma 82. I augmented the program I wrote for Enigma 82 to generate a suitable data file for the GLPK solver, although I found running glpsol on the output was taking too much time and memory (much more than for the equivalent problem for Enigma 82), so I used lp_solve which seemed to be both faster and use less memory (and can take the same format data). But even then it took 8h57m to complete an exhaustive search (although it finds the solution after only 5m26s, so if you already know the answer you can use the [[
-o]] option to find a solution of the appropriate size). A commercial ILP solver, that’s able to use multiple CPU cores would no doubt improve on this.Solution: The smallest number of matches required to spoil all quadrilaterals and triangles is 29.
Here’s the Python program I used to generate the data for the ILP solver.
from itertools import (product, combinations) from enigma import (irange, flatten, is_pairwise_distinct, Accumulator, printf) # consider the grid to be made of points (x, y) for 0 < x + y < N + 1 # lines along x=n, y=n, x=y. # points reachable from p along each axis def points(p, n): (x0, y0) = p px = set((x, y0) for x in irange(0, n - y0) if x != x0) py = set((x0, y) for y in irange(0, n - x0) if y != y0) m = x0 + y0 pz = set((x, m - x) for x in irange(0, m) if x != x0) return (px, py, pz) # points on the line between p0 and p1 def line(p0, p1): if p0[0] == p1[0]: return set((p0[0], y) for y in irange(min(p0[1], p1[1]), max(p0[1], p1[1]))) if p0[1] == p1[1]: return set((x, p0[1]) for x in irange(min(p0[0], p1[0]), max(p0[0], p1[0]))) m = sum(p0) if m == sum(p1): return set((x, m - x) for x in irange(min(p0[0], p1[0]), max(p0[0], p1[0]))) # are the points collinear? def is_collinear(p0, p1, p2): return ( p0[0] == p1[0] == p2[0] or p0[1] == p1[1] == p2[1] or sum(p0) == sum(p1) == sum(p2) ) # do the points identify a non-crossing quad? def is_quad(p0, p1, p2, p3): return ( is_pairwise_distinct(p0, p1, p2, p3) and # the lines (p0, p2) and (p1, p3) must not cross (not line(p0, p2).intersection(line(p1, p3))) and (not line(p0, p1).intersection(line(p2, p3))) and # check the none of the points are collinear (not is_collinear(p0, p1, p2)) and (not is_collinear(p0, p1, p3)) and (not is_collinear(p0, p2, p3)) and (not is_collinear(p1, p2, p3)) ) # generate all quads on the n x n triangular array def quads(n): seen = set() # consider each point in the lattice for y in irange(0, n): for x in irange(0, n - y): p0 = (x, y) # extend the quad along two of the vertices for (i, j) in combinations(points(p0, n), 2): # consider points along these axes for (p1, p2) in product(i, j): # and the final point completes the quad for p3 in set(flatten(points(p1, n))).intersection(flatten(points(p2, n))): if not is_quad(p0, p1, p2, p3): continue q = tuple(sorted((p0, p1, p2, p3))) if q in seen: continue yield (p0, p1, p2, p3) seen.add(q) # return the set of edges on this line def edges(p0, p1): if p0[0] == p1[0]: return set((2 * p0[0], 2 * y + 1) for y in range(min(p0[1], p1[1]), max(p0[1], p1[1]))) if p0[1] == p1[1]: return set((2 * x + 1, 2 * p0[1]) for x in range(min(p0[0], p1[0]), max(p0[0], p1[0]))) m = sum(p0) if m == sum(p1): return set((2 * x + 1, 2 * (m - x) - 1) for x in range(min(p0[0], p1[0]), max(p0[0], p1[0]))) return set() # return the set of edges in this quad def quad_edges(q): (p0, p1, p2, p3) = q return edges(p0, p1).union(edges(p0, p2), edges(p1, p3), edges(p2, p3)) # generate all triangles on the triangular array with side length n def tris(n): # consider each point in the lattice for y in irange(0, n): for x in irange(0, n - y): # triangle edge length for e in irange(1, n - (x + y)): yield ((x, y), (x + e, y), (x, y + e)) if not (e > y): yield ((x, y), (x + e, y), (x + e, y - e)) # return the set of edges in this triangle def tri_edges(t): (p0, p1, p2) = t return edges(p0, p1).union(edges(p0, p2), edges(p1, p2)) import sys N = 6 if len(sys.argv) < 2 else int(sys.argv[1]) # find all the quads on an NxN triangular array qs = list() for q in quads(N): # label the matchsticks (with an integer) m = set(e[1] + 100 * e[0] for e in quad_edges(q)) qs.append(m) printf("[quad {n}] {q} {m}", n=len(qs)) # find all the triangles on an NxN triangular array ts = list() for t in tris(N): m = set(e[1] + 100 * e[0] for e in tri_edges(t)) ts.append(m) printf("[tri {n}] {t} {m}", n=len(ts)) # find a minimum hitting set for the quads # output a file suitable for processing by GLPK # glpsol -m enigma82.mod -d enigma82.data from enigma import sprintf name = 'enigma144' fn = name + '.data' printf("generating {fn}") with open(fn, 'w') as f: f.write("data;\n") ss = qs + ts ms = sorted(set().union(*ss)) f.write(sprintf("set MATCHES := {m};\n", m=' '.join(str(x) for x in ms))) f.write(sprintf("set SHAPES := {s};\n", s=' '.join(str(x) for (x, s) in enumerate(ss)))) f.write(sprintf("param a: {m} :=\n", m=' '.join(str(x) for x in ms))) for (x, s) in enumerate(ss): f.write(sprintf(" {x} {s}\n", s=' '.join(str(int(x in s)) for x in ms))) f.write(";\n") f.write("end;\n")And the LP model:
/* sets - the matches and the shapes */ set MATCHES; set SHAPES; /* parameters - the inclusion matrix */ param a {i in SHAPES, j in MATCHES}; /* decision variables for the matches */ var x {j in MATCHES} binary >= 0; /* objective - minimum hitting set */ minimize size: sum{j in MATCHES} x[j]; /* constraint - the hitting set spoils each quad */ s.t. hit{i in SHAPES}: sum{j in MATCHES} a[i,j] * x[j] >= 1; end;I would then run the solver like this:
For N=6, lp_solve 5.5.2.0 finds a minimal solution in 8h57m. If [[
-o 29]] is applied the solution is found in 5m26s.For N=5, lp_solve 5.5.2.0 finds a minimal solution in 24.5s. If [[
-o 20]] is applied the solution is found in 122ms.As noted, this problem is an extension of Enigma 82, so here is an extension to my code that uses the Minimum Hitting Set implementation in MiniZinc that I used to solve Enigma 82, along with neater code for generating the triangles.
The size 6 grid is solved in 5.3 seconds, the size 7 grid in 22 seconds, the size 8 grid in 44 seconds, the size 9 grid in 4 minutes, and the size 10 grid in 12 hours (using the
highssolver).from enigma import (namedtuple, irange, subsets, ordered, cproduct, repeat, union, arg, printf) from utils_mzn import hitting_set # edge length of the triangular grid N = arg(6, 0, int) # a point has x and y values P = namedtuple("P", "x y") # points are non-negative (x, y) s.t. x + y <= N pts = set() for y in irange(0, N): pts.update(P(x, y) for x in irange(0, N - y)) # make a set of matchsticks from quad edges # each matchstick is of the form ((ax, ay), (bx, by)) def sticks(edges): ss = list() for (a, b) in edges: if a.x == b.x: ss.extend((P(a.x, y), P(a.x, y + 1)) for y in irange(a.y, b.y - 1)) elif a.y == b.y: ss.extend((P(x, a.y), P(x + 1, a.y)) for x in irange(a.x, b.x - 1)) elif a.x + a.y == b.x + b.y: ss.extend((P(a.x - k - 1, a.y + k + 1), P(a.x - k, a.y + k)) for k in irange(0, a.x - b.x - 1)) return ss # rotate a tuple of sticks rot = lambda ss: tuple(ordered(P(a.y, N - a.y - a.x), P(b.y, N - b.y - b.x)) for (a, b) in ss) # find quads quads = set() # choose a base level for the quadrilateral for y1 in irange(0, N - 2): # choose positions for the base points for (ax, bx) in subsets(irange(0, N - y1), size=2): # choose a top level for y2 in irange(y1 + 1, N - 1): k = y2 - y1 # choose position for the top points for (px, qx) in cproduct([[ax - k, ax], [bx - k, bx]]): if px < qx and P(px, y2) in pts and P(qx, y2) in pts: # find the matchsticks involved from the edges of the quadrilateral quad = sticks([ (P(ax, y1), P(bx, y1)), # base (P(px, y2), P(qx, y2)), # top (P(ax, y1), P(px, y2)), # left (P(bx, y1), P(qx, y2)), # right ]) # record all 3 rotations of the quadrilateral for qs in repeat(rot, quad, 2): quads.add(ordered(*qs)) # find triangles tris = set() # choose a base level for the triangle for y in irange(0, N - 1): # choose positions for the base points for (ax, bx) in subsets(irange(0, N - y), size=2): k = bx - ax # find matchsticks involved from the edges of the triangle tri = sticks([ (P(ax, y), P(bx, y)), # base (P(ax, y), P(ax, y + k)), # left (P(bx, y), P(bx - k, y + k)), # right ]) tris.add(ordered(*tri)) if not (k > y): # and the inverted triangle tri = sticks([ (P(ax, y), P(bx, y)), # top (P(ax + k, y - k), P(ax, y)), # left (P(bx, y - k), P(bx, y)), # right ]) tris.add(ordered(*tri)) printf("[N={N} -> {t} tris, {q} quads]", t=len(tris), q=len(quads)) # find a minimum hitting set hs = hitting_set(union([tris, quads]), solver="minizinc --solver highs") printf("minimum hitting set size {n}", n=len(hs))Solution: The smallest number of matches required to spoil all quadrilaterals and triangles is 29.