#!/usr/bin/env python
# coding: utf-8

# In[1]:


get_ipython().run_line_magic('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()


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


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


# # Tiling a 3x3 square with Wang tiles

# In[6]:


tiling = T0.solver(3,3).solve(solver='Coin')
tiling.tikz()


# In[7]:


tiling = T0.solver(3,3).solve(solver='glucose')
tiling.tikz()


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


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


# In[22]:


time_finding_dominoes(T2, solver='picosat', verbose=True, timelimit=100)


# 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


# 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


# 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[ ]:





# In[29]:


df_dominoes_T3 = compare_solvers_for_dominoes(T3, solvers, radius_list=[1,2,3], timelimit=timelimit)
df_dominoes_T3


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


# In[ ]: