In [1]:
%matplotlib inline
Dependencies¶
This notebook depends on the most recent development version of (slabbe)[https://pypi.python.org/pypi/slabbe] Optional Sage Package. The latest release can be installed with sage -pip install slabbe.
This notebook depends on the use of Glucose as a Sage optional package: https://trac.sagemath.org/ticket/26361
In [2]:
from slabbe import WangTileSet
Jeandel-Rao's 11 Wang tiles¶
We construct Jeandel-Rao's set of Wang tiles $T_0$ as well at tile sets $T_2$ and $T_3$ from https://arxiv.org/abs/1808.07768
In [3]:
tiles = [(2,4,2,1), (2,2,2,0), (1,1,3,1), (1,2,3,2), (3,1,3,3),
(0,1,3,1), (0,0,0,1), (3,1,0,2), (0,2,1,2), (1,2,1,4), (3,3,1,2)]
T0 = WangTileSet([map(str,t) for t in tiles])
T0.tikz()
Out[3]:
In [4]:
tiles = [(0, 0, 0, 1), (0, 1, 3, 1), (0, 2, 1, 2), (1, 1, 3, 1), (1, 2, 3, 2), (3, 1, 0, 2), (3, 1, 3, 3), (3, 3, 1, 2),
(203, 1, 210, 0), (203, 3, 211, 0), (210, 2, 211, 1), (210, 2, 231, 0), (211, 2, 213, 1),
(211, 2, 233, 0), (213, 1, 210, 1), (213, 1, 230, 0), (213, 3, 211, 1), (230, 1, 203, 0),
(231, 1, 203, 0), (233, 1, 213, 0)]
T2 = WangTileSet([map(str,t) for t in tiles])
T2.tikz()
Out[4]:
In [5]:
tiles = [(0, 0, 0, 1), (0, 1, 3, 1), (3, 1, 0, 2), (3, 1, 3, 3), (2030, 1, 2103, 0), (2033, 1, 2113, 0),
(2103, 1, 2110, 1), (2103, 1, 2310, 0), (2103, 3, 2111, 1), (2110, 2, 2131, 1), (2110, 2, 2331, 0), (2111, 2, 2133, 1),
(2113, 1, 2130, 1), (2113, 1, 2330, 0), (2130, 0, 2300, 0), (2130, 1, 2103, 1), (2131, 1, 2103, 1), (2133, 1, 2113, 1),
(2300, 0, 2030, 0), (2300, 1, 2033, 0), (2310, 1, 2033, 0), (2330, 0, 2130, 0), (2330, 1, 2133, 0), (2331, 1, 2133, 0)]
T3 = WangTileSet([map(str,t) for t in tiles])
T3.tikz()
Out[5]:
Tiling a 3x3 square with Wang tiles¶
In [6]:
tiling = T0.solver(3,3).solve(solver='Coin')
tiling.tikz()
Out[6]:
In [7]:
tiling = T0.solver(3,3).solve(solver='glucose')
tiling.tikz()
Out[7]:
Reduction to SAT:¶
The reduction to SAT of the tiling of a $n\times n$ square by $t$ Wang tiles needs $t\times n^2$ variables and at most $2n^2 t^2$ constraints.
In [8]:
for k in range(10,70,10):
print k, T0.solver(k,k).sat_solver('cryptominisat')
Reduction to MILP¶
The reduction to MILP of the tiling of a $n\times n$ square needs $11\times n^2$ variables and $2n(n-1)+n^2$ constraints.
In [9]:
for k in range(10,70,10):
print k, T0.solver(k,k).milp()
Finding a nxn valid Wang square tiling : comparing solvers¶
In [10]:
import time,signal
def time_solving_squares(wang_tile_set, solver='gurobi', maxsize=1000, step=1, timelimit=5, verbose=False):
r"""
Return a list of tuples (n,t) where t is the time spent (in seconds) finding
a valid Wang tiling of the n times n square.
EXAMPLES::
sage: time_solving_squares(solver='gurobi', timelimit=1)
[(1, 0.0019371509552001953),
(2, 0.004050016403198242),
(3, 0.01121377944946289),
(4, 0.015249013900756836),
(5, 0.0812530517578125),
(6, 0.08473992347717285),
(7, 0.23214101791381836),
(8, 0.33846092224121094),
(9, 0.6193099021911621),
(10, 1)]
"""
L = []
for i in range(1, maxsize, step):
try:
signal.alarm(timelimit)
start = time.time()
solution = wang_tile_set.solver(i,i).solve(solver=solver, ncpus=8)
end = time.time()
except (AlarmInterrupt, ValueError): # dancing_links catches the alarm and returns valueError
L.append( (i,timelimit) )
break
else:
signal.alarm(0) # cancels the alarm
duration = end - start
L.append( (i,duration) )
if verbose:
print( (i,duration) )
return L
Testing:
In [11]:
time_solving_squares(T0, solver='Gurobi', timelimit=1)
Out[11]:
In [12]:
import pandas as pd
def compare_solvers(wang_tile_set, solvers, maxsize=1000, step=1, timelimit=5):
wang_tile_set = WangTileSet(list(wang_tile_set)) # reinitialize cached values
series_dict = {}
for solver in solvers:
print('Testing solver {} ...'.format(solver))
L = time_solving_squares(wang_tile_set, solver=solver, maxsize=maxsize, step=step, timelimit=timelimit)
index,data = zip(*L)
series = pd.Series(data=data, index=index, dtype=float)
series_dict[solver] = series
print(' -> solved {} instances'.format(len(series)))
sorted_solvers = sorted(solvers, key=lambda solver:len(series_dict[solver]))
return pd.DataFrame(series_dict, columns=sorted_solvers)
Define the list of solvers:
In [13]:
LPsolvers = ['Gurobi', 'LP']
slowLPsolvers = ['GLPK', 'PPL']
SATsolvers = ['cryptominisat', 'picosat', 'glucose']
uninstalled = ['Coin', 'GLPK/exact', 'CPLEX', 'CVXOPT', 'rsat'] # to be added to the comparison one day
solvers = slowLPsolvers + SATsolvers + LPsolvers + ['dancing_links']
In [14]:
df_square_T0 = compare_solvers(T0, solvers, timelimit=10)
In [41]:
ax = df_square_T0.plot(figsize=(15,10), marker='o')
ax.set(xlabel="n, the square size", ylabel="time (seconds)",
title='Time to find a valid Wang Tiling of a n x n square using Jeandel Rao 11 tile set $T_0$')
ax.get_figure().savefig('T0_square_tilings.svg')
In [16]:
df_square_T2 = compare_solvers(T2, solvers, timelimit=10)
In [42]:
ax = df_square_T2.plot(figsize=(15,10), marker='o')
ax.set(xlabel="n, the square size", ylabel="time (seconds)",
title='Time to find a valid Wang Tiling of a n x n square using the 20 Wang tile set $T_2$')
ax.get_figure().savefig('T2_square_tilings.svg')
In [18]:
df_square_T3 = compare_solvers(T3, solvers, timelimit=10)
In [43]:
ax = df_square_T3.plot(figsize=(15,10), marker='o')
ax.set(xlabel="n, the square size", ylabel="time (seconds)",
title='Time to find a valid Wang Tiling of a n x n square using the 24 Wang tile set $T_3$')
ax.get_figure().savefig('T3_square_tilings.svg')
Finding all dominoes admitting a surrounding : comparing solvers¶
In [20]:
import time,signal
def time_finding_dominoes(wang_tile_set, solver='gurobi', radius_list=[1,2,3], timelimit=5, verbose=False):
r"""
EXAMPLES::
sage: time_finding_dominoes(solver='gurobi', timelimit=30)
[(1, 2.8133392333984375e-05),
(2, 1.2874603271484375e-05),
(3, 5.506983995437622),
(4, 10.366616010665894),
(5, 20.101346015930176),
(6, 30)]
"""
L = []
for radius in radius_list:
try:
signal.alarm(timelimit)
start = time.time()
dominoes = {}
for direction in [1, 2]:
dominoes[direction] = wang_tile_set.dominoes_with_surrounding(radius=radius, i=direction, solver=solver, ncpus=1)
if verbose:
print('{} computes {} dominoes in the direction e{} with'
' surrounding radius {}'.format(solver, len(dominoes[direction]), direction, radius))
end = time.time()
except (AlarmInterrupt, ValueError): # dancing_links catches the alarm and returns valueError
#L.append( (radius,timelimit, None, None) )
break
else:
signal.alarm(0) # cancels the alarm
duration = end - start
L.append( (radius,duration,len(dominoes[1]), len(dominoes[2]) ) )
if verbose:
print('Computation time for radius {}: {} seconds'.format(radius,duration))
return L
In [21]:
time_finding_dominoes(T2, solver='dancing_links', verbose=True, timelimit=10)
Out[21]:
In [22]:
time_finding_dominoes(T2, solver='picosat', verbose=True, timelimit=100)
Out[22]:
In [23]:
import pandas as pd
def compare_solvers_for_dominoes(wang_tile_set, solvers, radius_list=[1,2,3], timelimit=5, verbose=False):
wang_tile_set = WangTileSet(list(wang_tile_set)) # reinitialize cached values
series_dict = {}
for solver in solvers:
print('Testing solver {} ...'.format(solver))
L = time_finding_dominoes(wang_tile_set, solver=solver, radius_list=radius_list,
timelimit=timelimit, verbose=verbose)
if L:
index,data,dir1,dir2 = zip(*L)
print(' -> Found {} and {} dominoes in the direction'
' 1 and 2 resp for radius in {}'.format(dir1, dir2, radius_list))
print(' -> time data: {}'.format(data))
else:
index = []
data = []
series = pd.Series(data=data, index=index, dtype=float)
series_dict[solver] = series
sorted_solvers = sorted(solvers, key=lambda solver:len(series_dict[solver]))
return pd.DataFrame(series_dict, columns=sorted_solvers)
In [24]:
solvers = ['dancing_links'] + SATsolvers + LPsolvers # + slowLPsolvers
timelimit = 1000
In [25]:
df_dominoes_T0 = compare_solvers_for_dominoes(T0, solvers, radius_list=[1,2,3], timelimit=timelimit)
df_dominoes_T0
Out[25]:
In [44]:
ax = df_dominoes_T0.plot(figsize=(15,10), marker='o')
ax.set(xlabel="surrounding radius", ylabel="time (seconds)",
title='Finding the dominoes allowing a surrounding of given radius '
'for the Jeandel Rao 11 tile set $T_0$ in less than {} seconds'.format(timelimit))
ax.get_figure().savefig('T0_dominoes_surrounding.svg')
In [27]:
df_dominoes_T2 = compare_solvers_for_dominoes(T2, solvers, radius_list=[1,2,3], timelimit=timelimit)
df_dominoes_T2
Out[27]:
In [45]:
ax = df_dominoes_T2.plot(figsize=(15,10), marker='o')
ax.set(xlabel="surrounding radius", ylabel="time (seconds)",
title='Finding all dominoes allowing a surrounding of given radius '
'for the 20 Wang tile set $T_2$ in less than {} seconds'.format(timelimit))
ax.get_figure().savefig('T2_dominoes_surrounding.svg')
In [29]:
df_dominoes_T3 = compare_solvers_for_dominoes(T3, solvers, radius_list=[1,2,3], timelimit=timelimit)
df_dominoes_T3
Out[29]:
In [46]:
ax = df_dominoes_T3.plot(figsize=(15,10), marker='o')
ax.set(xlabel="surrounding radius", ylabel="time (seconds)",
title='Finding all dominoes allowing a surrounding of given radius '
'for the 24 Wang tile set $T_3$ in less than {} seconds'.format(timelimit))
ax.get_figure().savefig('T3_dominoes_surrounding.svg')
