In [1]:
%matplotlib inline
%run ../python-libs/bits.py
%run ../python-libs/timing.py
%run ../python-libs/graycodes.py
%run ../python-libs/symbols.py
In [2]:
import sys
sys.setrecursionlimit(10000000)
The following code implements a backtracking function that yields a sequence of tiling of a given board, possibly containing some forbidden cells, according to a set of available shapes.
In [3]:
import math, random
def polyominoes(dim, shapes, availables='ones', max_depth_reached=None,
forbidden=[], simulated_annealing={'enabled':False, 'energy':lambda s: 0, 'scheduler':None},
nogood=lambda shape_placement, fewer_positions, iso, availables: False,
pruning=lambda coord, positions, shapes: False):
"""
Returns an iterable of arrangements for the Polyominos problem.
"""
rows, cols = dim
sol = []
if not availables or availables == 'ones':
availables = {s.name:1 for s in shapes} #[1]*len(shapes)
elif availables == 'inf':
availables = {s.name:-1 for s in shapes} #[-1]*len(shapes) # trick: represent ∞ piece availability by decreasing negative numbers
def place(S, positions):
for r, c in positions:
S = clear_bit(S, r + rows*c)
return S
initial_temperature = sum(availables[s.name] for s in shapes)
schedulers = {'cauchy': lambda t: initial_temperature/(t),
'standard':lambda t: initial_temperature*.95**t,
'log':lambda t: initial_temperature/math.log(t),}
def simulated_annealing_choice(step, competitors, temperature_scheduler, boltzman_const=1):
candidate, outsider = competitors
energy = simulated_annealing['energy']
delta = energy(outsider) - energy(candidate) # in this way it is always positive
if delta > 0: return outsider
elif not delta: return candidate if random.random() > 0.5 else outsider
# therefore: delta < 0
prob = math.exp(delta/(boltzman_const * temperature_scheduler(step)))
return outsider if random.random() < prob else candidate
"""
freq = int(prob*accuracy)
population = ([o]*freq) + ([s]*(accuracy-freq))
return random.choice(population)
"""
def gen(positions, attempts):
p = low_bit(positions)
c, r = divmod(p, rows)
if pruning((r,c), positions, {s for s in shapes if availables[s.name]}):
raise StopIteration()
for i, s in enumerate(shapes):
if not availables[s.name]: continue
#________________________________________________________________
if simulated_annealing['enabled']:
left_shapes = {s for s in shapes if availables[s.name]} - {s}
annealing_shape = left_shapes.pop() if left_shapes else s
s = simulated_annealing_choice(step=sum(availables[s.name] for s in shapes),
competitors=(s, annealing_shape),
temperature_scheduler=schedulers[simulated_annealing['scheduler']])
#________________________________________________________________
for j, iso in enumerate(s.isomorphisms(r, c)):
if all(0 <= rr < rows and 0 <= cc < cols and is_on(positions, rr + rows*cc) for rr, cc in iso):
fewer_positions = place(positions, iso)
if nogood((s,p), fewer_positions, iso, availables): continue
availables[s.name] -= 1
sol.append((s, positions, (r,c), iso),)
yield from [sol] if not (fewer_positions and attempts) else gen(fewer_positions, attempts-1)
sol.pop()
availables[s.name] += 1
return gen(place(set_all(rows*cols), forbidden), max_depth_reached or -1)
In [127]:
import math, random
def polyominoes_bidirectional(dim, shapes, availables='ones', max_depth_reached=None,
forbidden=[], simulated_annealing={'enabled':False, 'energy':lambda s: 0, 'scheduler':None},
nogood=lambda shape_placement, fewer_positions, iso, availables: False,
pruning=lambda coord, placement, shapes: False):
"""
Returns an iterable of arrangements for the Polyominos problem.
"""
rows, cols = dim
sol = []
if not availables or availables == 'ones':
availables = {s.name:1 for s in shapes} #[1]*len(shapes)
elif availables == 'inf':
availables = {s.name:-1 for s in shapes} #[-1]*len(shapes) # trick: represent ∞ piece availability by decreasing negative numbers
def coord_in_alternative(current_cell):
r, c = current_cell
return rows-1-r, cols-1-c
def place(positions, iso):
current_positions, alternative_positions = positions
for cell in iso:
r, c = cell
current_positions = clear_bit(current_positions, r + rows*c)
ar, ac = coord_in_alternative(cell)
alternative_positions = clear_bit(alternative_positions, ar + rows*ac)
return current_positions, alternative_positions
initial_temperature = sum(availables[s.name] for s in shapes)
schedulers = {'cauchy': lambda t: initial_temperature/(t),
'standard':lambda t: initial_temperature*.95**t,
'log':lambda t: initial_temperature/math.log(t),}
def simulated_annealing_choice(step, competitors, temperature_scheduler, boltzman_const=1):
candidate, outsider = competitors
energy = simulated_annealing['energy']
delta = energy(outsider) - energy(candidate) # in this way it is always positive
if delta > 0: return outsider
elif not delta: return candidate if random.random() > 0.5 else outsider
# therefore: delta < 0
prob = math.exp(delta/(boltzman_const * temperature_scheduler(step)))
return outsider if random.random() < prob else candidate
"""
freq = int(prob*accuracy)
population = ([o]*freq) + ([s]*(accuracy-freq))
return random.choice(population)
"""
def gen(positions, attempts, insertion_side):
current_positions, alternative_positions = positions
p = low_bit(current_positions)
c, r = divmod(p, rows)
if pruning((r,c), positions, {s for s in shapes if availables[s.name]}):
raise StopIteration()
for i, s in enumerate(shapes):
if not availables[s.name]: continue
#________________________________________________________________
if simulated_annealing['enabled']:
left_shapes = {s for s in shapes if availables[s.name]} - {s}
annealing_shape = left_shapes.pop() if left_shapes else s
s = simulated_annealing_choice(step=sum(availables[s.name] for s in shapes),
competitors=(s, annealing_shape),
temperature_scheduler=schedulers[simulated_annealing['scheduler']])
#________________________________________________________________
for j, iso in enumerate(s.isomorphisms(r, c)):
can_be_place_in_current = all(0 <= rr < rows
and 0 <= cc < cols
and is_on(current_positions, rr + rows*cc)
for rr, cc in iso)
can_be_place_in_alternative = all(0 <= rr < rows
and 0 <= cc < cols
and is_on(alternative_positions, rr + rows*cc)
for rr, cc in map(coord_in_alternative, iso))
if can_be_place_in_current and can_be_place_in_alternative:
fewer_positions = place(positions, iso) # `place` consumes and produces a pair of bitmasks
fewer_current_positions, fewer_alternative_positions = fewer_positions
if nogood((s,p), fewer_current_positions, iso, availables): continue
availables[s.name] -= 1
placed_iso = iso if not insertion_side else list(map(coord_in_alternative, iso))
sol.append((s, positions, max(placed_iso), placed_iso),)
swapped_fewer_positions = fewer_alternative_positions, fewer_current_positions
holes_left = fewer_current_positions and fewer_alternative_positions
yield from [sol] if not (holes_left and attempts) else gen(swapped_fewer_positions,
attempts-1,
1-insertion_side)
sol.pop()
availables[s.name] += 1
current_positions = set_all(rows*cols)
for r, c in forbidden:
current_positions = clear_bit(current_positions, r + rows*c)
initial_positions = current_positions, set_all(rows*cols)
return gen(positions=initial_positions,
attempts=max_depth_reached or -1,
insertion_side=0)
In the following cell we define a function that generates all parallelogram polyominoes with semiperimeter sp.
In [33]:
def parallelogram_polyominoes(sp):
"""
Returns the set of shapes denoting all *parallelogram polyominoes* of semiperimeter `sp`.
The relation to the main problem is $sp = n + 2$, where the board edge is $2^{n}$.
"""
import itertools
steps = [(1,0), # go to cell below
(0,1),] # go to cell on the right
# we assume the canonical ordering inside a board, namely ascending
# order from top to bottom, and from left to right, as in cells above.
# first of all build Catalan paths that have no intersection point,
# using only vertical and horizontal steps (here we do not distinguish
# among upward/downward vertical steps, what is important is to remain
# consistent to board ordering).
initial_prefix = [(0,0)] # always start placing cell at anchor coordinate (r, c)
prefixes = {tuple(initial_prefix)}
candidates = set()
n = sp - 2
while prefixes:
prefix = prefixes.pop()
if len(prefix) > n:
candidates.add(prefix)
continue
last_row_offset, last_col_offset = prefix[-1]
dominating_prefixes = [(last_row_offset + ro, last_col_offset + co) for ro, co in steps]
dominated_prefixes = [(last_row_offset + ro, last_col_offset + co) for ro, co in steps]
for sub, dom in [(sub,dom) for sub, dom in zip(dominated_prefixes, dominating_prefixes) if sub <= dom]:
prefixes.add(prefix + (sub,))
prefixes.add(prefix + (dom,))
parallelograms = {(tuple(dom), tuple(sub)) for sub, dom in itertools.product(candidates, repeat=2)
if sub[-1] == dom[-1]}
polyominoes = {frozenset(fst) | frozenset(snd) for fst, snd in parallelograms}
def fill(polyomino):
r_max, c_max = max(r for r,_ in polyomino), max(c for _,c in polyomino)
filled = set()
for coord in itertools.product(range(1, r_max), range(1, c_max)):
if coord in polyomino: continue
row = [(rr, cc) for rr, cc in polyomino if coord[0] == rr]
column = [(rr, cc) for rr, cc in polyomino if coord[1] == cc]
if (min(column, default=coord) < coord < max(column, default=coord) and
min(row, default=coord) < coord < max(row, default=coord)):
filled.add(coord)
#print("filled: ", polyomino, polyomino | filled)
return polyomino | filled
return {fill(p) for p in polyominoes}
In [34]:
def rotate_clockwise(shape):
"""
Returns a shape rotated clockwise by π/2 radians.
"""
clockwise = [(c, -r) for r, c in shape] # by matrix multiplication for rotations
#print(clockwise)
co, ro = min((c, r) for r, c in clockwise)
#print(ro, co)
return frozenset((r-ro, c-co) for r, c in clockwise)
def rotate_clockwise_isomorphisms(isos):
rotations = set()
for i in isos:
rotating = i
rotating = rotate_clockwise(rotating)
rotations.add(rotating)
rotating = rotate_clockwise(rotating)
rotations.add(rotating)
rotating = rotate_clockwise(rotating)
rotations.add(rotating)
#print(isos | frozenset(rotations))
return isos | rotations
#______________________________________________________________________________
def vertical_mirror(shape):
m = max([c for r, c in shape])
return frozenset((r, m-c) for r, c in shape)
def vertical_isomorphism(isos):
return isos | frozenset(vertical_mirror(i) for i in isos)
#______________________________________________________________________________
def make_shapes(primitive_polyos, isomorphisms=lambda isos: isos):
prefix = "polyo"
def o(j, isos):
return shape_spec(name='{}{}'.format(prefix, j),
isomorphisms=lambda r, c: [ [(r+ro, c+co) for ro, co in iso] for iso in isos])
return [o(j, isomorphisms({p})) for j, p in enumerate(primitive_polyos)]
Sandbox to test:
In [35]:
c = frozenset({(0,0), (1,0), (1,1), (1,2), (1,3)})
print(c)
c = rotate_clockwise(c)
print(c)
c = rotate_clockwise(c)
print(c)
c = rotate_clockwise(c)
print(c)
In [36]:
def describe(semiperimeter, polyominoes):
from IPython.display import Markdown
n, catalan = semiperimeter - 2, semiperimeter - 1
size=2**n
area = 0
for p in polyominoes:
area += len(p)
theoretical_area = 4**n
code='''
**Problem instance**:
- using {pp} parallelogram polyominoes, which is the catalan number $c_{catalan}$ (0-based indexing)
- each polyomino has semiperimeter $sp$ equals to {sp}, so let $n=sp-2={n}$
- board edges have size $2^{n}={size}$ each
- theoretically, area is known to be $4^{n}={fa}$, which *is {pred}* equal to counted area {ca}'''.format(
sp=semiperimeter, n=n, size=size, pp=len(polyominoes), fa=theoretical_area,
pred='not' if theoretical_area != area else '', ca=area, catalan=catalan)
return Markdown(code), n, catalan, polyominoes, theoretical_area, area, size
Using function describe above, we state the problem instance clearly with:
In [37]:
semiperimeter = 6
description, *facts = describe(semiperimeter, parallelogram_polyominoes(semiperimeter))
n, catalan, polyos, theoretical_area, area, size = facts
description
Out[37]:
The order in which shapes are provided is very important since it affects the quality of backtracking:
In [38]:
#shapes = make_shapes(pp, isomorphisms=vertical_isomorphism)#, isomorphisms=lambda isos: vertical_isomorphism(rotate_clockwise_isomorphisms(isos)))
shapes = make_shapes(polyos)#, isomorphisms=lambda isos: vertical_isomorphism(rotate_clockwise_isomorphisms(isos)))
def convex_hull(iso):
lower_row, lower_col = min(iso)
greatest_row, greatest_col = max(iso) # find the "convex hull" respect to the cell at the bottom-right most location
return [(r, c)
for r in range(lower_row, greatest_row+1)
for c in range(lower_col, greatest_col+1)] # convex hull area
def area_key(s, weights={0:0, 1:0, 2:0, 3:0, 4:1}):
import random
isos = s.isomorphisms(0,0) # place the shape somewhere, the origin is pretty good
iso = isos.pop() # all isomorphisms, of course, have the same are, therefore pick the first one
filled_area = len(iso) # since a piece is represent as an iterable of coordinates `(r, c)`
convex_hull_area = len(convex_hull(iso))
#key = filled_area*(weights['filled'] + weights['convex']/convex_hull_area) + weights['bias']*random.random() # return the ratio of the filled area respect the convex hull
#key = 0*max(iso)[1] + filled_area * convex_hull_area / max_convex_hull
def key(x):
max_r, max_c = max(iso)
max_convex_hull = (semiperimeter-1)**2
#return max_r + max_c*x + (filled_area * convex_hull_area / max_convex_hull)*x**2
return (weights[0]*random.random()*x**0 +
weights[1]*max_r*x**1 +
weights[2]*max_c*x**2 +
weights[3]*filled_area*x**3 +
weights[4]*(convex_hull_area / max_convex_hull)*x**4)
return key(2)
shapes = sorted(shapes, key=area_key, reverse=True)
Represent all available parallelogram polyominoes:
In [39]:
polyominoes_boards = [symbols_pretty([(shape, None, (0,0), iso)],
(semiperimeter, semiperimeter),
{shape.name:'▢'},
joiner=lambda board: board,
axis=True)
for i, shape in enumerate(shapes) for iso in shape.isomorphisms(0,0)]
group=5
for i in range(0, len(polyominoes_boards), group):
for k in range(semiperimeter+2):
grouped_board = []
for j in range(i, i + group):
if j >= len(polyominoes_boards): break
grouped_board.append(polyominoes_boards[j][k])
print(' '.join(grouped_board))
Thoughts:
available array holds (for instance, that every piece has been used so that no pieces
left for each solution)
In [204]:
############### DIRTY ######################################
dim = (size, size)
def nogood_predicate(placed_shape, fewer_positions, iso, fringe, availables):
if not fewer_positions:
return False
shape, p = placed_shape
#print("given fringe ", fringe)
rows=size
c,r = divmod(p, rows)
#if r == rows-1 or c == size-1: return True
def cell_fringe(cell):
r, c = cell
return {(r, c-1)}#, (r+1, c)}
#return {(r, c-1),(r-1, c),(r+1, c),(r, c+1),} # this is the complete fringe
max_iso_cell_r, max_iso_cell_c = max(iso)
fringe_iso = {(rr, cc)
for cell in iso
for (rr,cc) in cell_fringe(cell)
if 0 <= r <= rr < max_iso_cell_r + 1
and 0 <= c-1 <= cc < max_iso_cell_c
and is_on(fewer_positions, rr + rows*cc)}
#if max_iso_cell_r+1 == rows-1 and is_on(fewer_positions, max_iso_cell_r+1 + rows*(max_iso_cell_c-1)):
# fringe_iso.add((max_iso_cell_r+1, max_iso_cell_c-1),)
if False:
for rr, cc in iso:
if cc == size-1:
if rr-1 >= 0 and is_on(fewer_positions, rr-1 + rows*cc):
fringe_iso.add((rr-1, cc),)
if rr+1 < rows and is_on(fewer_positions, rr+1 + rows*cc):
fringe_iso.add((rr+1, cc),)
#print("standard fringe_iso ", fringe_iso)
hole = set()
for rr, cc in fringe_iso:
back_cc = cc
while back_cc >= c and is_on(fewer_positions, rr + rows*back_cc):
hole.add((rr, back_cc),)
back_cc -= 1
if hole:
#print(hole)
min_c, min_r = min((c, r) for r, c in hole)
chances = {s for s in shapes
if s != shape
and availables[s.name]
and all(is_on(fewer_positions, rr + rows*cc)
for rr, cc in s.isomorphisms(min_r, min_c).pop())}
#and len(set(s.isomorphisms(0,0).pop())) <= len(hole)}
return not chances
else:
return False
if not fringe_iso: return False # `iso` placement fills a hole exactly
min_c, min_r = min((c, r) for r, c in fringe_iso)
min_cell = min_r, min_c
#print(min_cell)
fringe_iso = {(rr,cc) for (rr,cc) in fringe_iso if rr != min_r} | {min_cell} # discard any cell occurring more than one on the row where the min cell lies
#print('doubles-free fringe ', fringe_iso)
class unbounded_area(Exception): pass
def ant(cell, remaining_steps):
r, c = cell
free = is_on(fewer_positions, r + rows*c)
if not remaining_steps:
return set()
return {cell} if free else set()
covered_bottom = ant((r+1, c), remaining_steps-1) if r < rows-1 else set()
covered_right = ant((r, c+1), remaining_steps-1) if c < rows-1 else set()
return {cell} | covered_bottom | covered_right
for cell in fringe | fringe_iso:
rr, cc = cell
covered = ant(cell, semiperimeter-1)
#print("covered ", covered)
filling_shapes = {s for s in shapes
if s != shape
and availables[s.name]
and set(s.isomorphisms(rr, cc).pop()) <= covered}
if not filling_shapes:
#print("not filling")
return True
fringe.update(fringe_iso)
#print("updated fringe ", fringe)
return False
rows=size
#if r == rows-1: return True
"""
latest_placed_greater_row, latest_placed_greater_col, = max(iso)
if latest_placed_greater_row == rows-1: # the latest shape has its greatest cell on the last row
for q in range(0, semiperimeter-1):
if latest_placed_greater_col-q > 0 and is_on(fewer_positions, rows*(latest_placed_greater_col-q)-1):
return True
"""
columns_ahead = 1
for ahead in range(p+1, rows*(c+1+columns_ahead)): # look ahead up to the end of the column where last shape has been placed
if not is_on(fewer_positions, ahead): continue
next_free_col, next_free_row, = divmod(ahead, rows)
if next_free_row == rows-1 and next_free_col > columns_ahead:
return True
"""if next_free_row == rows-1:
for q in range(2, semiperimeter):
if not is_on(fewer_positions, rows*(next_free_col+q)-2):
return True
"""
if not (is_on(fewer_positions, (next_free_row+1) + rows*next_free_col) or
is_on(fewer_positions, next_free_row + rows*(next_free_col+1))):
return True
return False
In [40]:
dim = (size, size)
def nogood_predicate(placed_shape, fewer_positions, iso, availables):
if not fewer_positions:
return False
shape, p = placed_shape
rows=size
c,r = divmod(p, rows)
max_r, max_c = max(iso)
if False and (max_r == rows-2 or max_c == size-2):
return True
for (rr, cc) in iso:
closing_downward = cc and is_on(fewer_positions, rr+rows*(cc-1)) and (rr+1==rows or not is_on(fewer_positions, rr+1+rows*(cc-1)))
closing_upwards = rr and is_on(fewer_positions, rr-1+rows*cc) and (cc+1==size or not is_on(fewer_positions, rr-1+rows*(cc+1)))
if closing_downward or closing_upwards:
return True
#return False
for i in range(size):
if not fewer_positions:
break
ahead, fewer_positions = low_bit_and_clear(fewer_positions)
next_free_col, next_free_row, = divmod(ahead, rows)
if next_free_row >= r and next_free_col > c:
break
#if next_free_row == rows-1 and not is_on(fewer_positions, (next_free_row-1) + rows*next_free_col):
# return True
#if next_free_row+1 == rows or next_free_col+1 == size:
# return True
if not (is_on(fewer_positions, (next_free_row+1) + rows*next_free_col)
or is_on(fewer_positions, next_free_row + rows*(next_free_col+1))):
return True
#continue
if all(not is_on(fewer_positions, next_free_row+1+rows*cc)
for cc in range(next_free_col, next_free_col+3) if cc < size):
return True
return False
In [41]:
energy=lambda item: area_key(item, weights={0:0, 1:0, 2:0, 3:1, 4:0, })
def not_insertion_on_edges(coord, positions, shapes):
r, c = coord
if r >= size-1 or c >= size-1:
return True
return False
for s in shapes:
iso = s.isomorphisms(r, c).pop()
max_r, max_c = max(iso)
if max_r < size and max_c < size:
break
else:
return True
return False
# 114 in 1'56'' for sp=7
polys_sols = polyominoes(dim, shapes, availables="ones", forbidden=[], max_depth_reached=38,
simulated_annealing={'enabled':False, 'energy':energy, 'scheduler':'cauchy'},
nogood=nogood_predicate, pruning=not_insertion_on_edges)
pretty_tilings = markdown_pretty(polys_sols, dim, shapes, raw_text=True)
In [54]:
print(next(pretty_tilings))
and pretty print some experiments:
In [ ]:
symbols = lower_greek_symbols() + capital_greek_symbols() + latin_symbols()
symbols_map = {s.name:symbols[i] for i,s in enumerate(shapes)}
count = 0
with open('sp_6.txt', 'a') as f:
for s in polys_sols:
f.write('sol index {}:\n{}\n'.format(count, ppretty(s, dim, symbols_map)))
count += 1
count
Maybe quadtrees should be helpful:
<span xmlns:dct="http://purl.org/dc/terms/" property="dct:title">Backtracking tutorial</span> by <a xmlns:cc="http://creativecommons.org/ns#" href="massimo.nocentini@unifi.it" property="cc:attributionName" rel="cc:attributionURL">Massimo Nocentini</a> is licensed under a Creative Commons Attribution-NonCommercial 4.0 International License.
Based on a work at <a xmlns:dct="http://purl.org/dc/terms/" href="https://github.com/massimo-nocentini/competitive-programming/blob/master/tutorials/backtrack.ipynb" rel="dct:source">https://github.com/massimo-nocentini/competitive-programming/blob/master/tutorials/backtrack.ipynb</a>.