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

# # Substitutive structure of Jeandel-Rao aperiodic tilings
# 
# Author: Sébastien Labbé, August 23th 2018.
# 
# ## Prerequisites

# Sage version 8.3 or above should work:

# In[1]:


version()


# The installation of the [optional Sage package slabbe v0.4.3](https://pypi.org/project/slabbe/0.4.3) is a prerequisite. Remove the dash sign (`#`) and evaluate the following cell to install slabbe optional package:

# In[2]:


#!sage -pip install slabbe==0.4.3


# The installation of [graphviz](https://graphviz.org/) and of the optional Sage package dot2tex is a prerequisite to construct the graphs at the end of the notebook.

# In[3]:


#!sage -pip install dot2tex


# ## Choosing a solver
# 
# The time it takes to evaluate all cells of this notebook depend on the chosen available solver.
# 
# Solver | Time to run this notebook
# --- | ---
# *Gurobi* | about 4 minutes
# *dancing_links* | about 69 minutes
# *other MILP solver* | unknown
# 
# The above timing were computed on a 2017 computer with 8 available cpus.
# 
# Gurobi is a proprietary software, but it is free for researchers and students. See this [tutorial to make Gurobi available through Sage](https://doc.sagemath.org/html/en/thematic_tutorials/linear_programming.html#using-cplex-or-gurobi-through-sage).

# In[4]:


from sage.doctest.external import has_gurobi
if has_gurobi():
    solver = 'gurobi'
else:
    solver = 'dancing_links'
print("We are using solver = '{}'".format(solver))


# ## Colors setup

# In[5]:


from collections import defaultdict
color = defaultdict(lambda : 'white')
color.update({0:'white', 1:'red', 2:'cyan', 3:'green', 4:'lightgray'})
color.update({str(k):v for k,v in color.items()})


# ## Definition of Jeandel-Rao Tiles $\mathcal{T}_0$

# In[6]:


from slabbe import WangTileSet, Substitution2d
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)]
tiles = [map(str,t) for t in tiles]
T0 = WangTileSet(tiles)
T0.tikz(font=r'\small', size=1, ncolumns=11, color=color)


# ## Computing $\mathcal{T}_1$

# In[7]:


get_ipython().run_line_magic('time', 'T0.find_markers(i=2, radius=1, solver=solver) # 47s with dancing_links, 1.4s with Gurobi')


# In[8]:


M0 = [0,1]
T1,s0,_ = T0.find_substitution(M=M0, i=2, side='left', solver=solver)


# In[9]:


T1.tikz(font=r'\small', size=1.2)


# In[10]:


show(s0)


# ## Computing $\mathcal{T}_2$

# In[11]:


get_ipython().run_line_magic('time', 'T1.find_markers(i=2, radius=1, solver=solver) #1min 16s with dancing_links, 2s with Gurobi')


# In[12]:


M1 = [8, 9, 10, 11, 12]
T2,s1,_ = T1.find_substitution(M=M1, i=2, side='left', solver=solver)


# In[13]:


T2.tikz(font=r'\small', size=1.2)


# In[14]:


show(s1)


# ## Computing $\mathcal{T}_3$

# In[15]:


get_ipython().run_line_magic('time', 'T2.find_markers(i=2, radius=2, solver=solver) # 5min 25s with dancing_links, 15s with Gurobi')


# In[16]:


M2 = [8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19]
T3,s2,_ = T2.find_substitution(M=M2, i=2, side='left', radius=2, solver=solver)


# In[17]:


T3.tikz(font=r'\small', size=1.2)


# In[18]:


show(s2)


# ## Computing $\mathcal{T}_4'$

# In[19]:


get_ipython().run_line_magic('time', 'T3.find_markers(i=2, radius=3, solver=solver) # 10min 12s with dancing_links, 59s with Gurobi')


# In[20]:


M3 = [0, 1, 2, 3]
T4p,s3p,_ = T3.find_substitution(M=M3, i=2, side='right', radius=3, solver=solver)


# In[21]:


T4p.tikz(font=r'\small', size=1.2)


# In[22]:


show(s3p)


# ## Check that we can remove 2 tiles from $\mathcal{T}_4'$
# 
# Only Gurobi can solve this in a reasonable amount of time (6s)

# In[23]:


def check_2_tiles_to_remove_from_T4p():
    assert T4p[18] == ('21103', '1', '23310', '0')
    S = T4p.solver(width=71, height=9, preassigned_tiles={(35,4):18})
    return S.has_solution(solver='Gurobi')
get_ipython().run_line_magic('time', 'check_2_tiles_to_remove_from_T4p()')


# ## Computing $\mathcal{T}_4$

# In[24]:


forbidden = ('21103', '1', '23310', '0'), ('23310', '0', '21330', '0')
id_tiles = [(i,t) for (i,t) in enumerate(T4p) if t not in forbidden]
indices,tiles = zip(*id_tiles)
d = dict(enumerate(indices))
iota = Substitution2d.from_permutation(d)
T4 = WangTileSet(tiles)
T4.tikz(font=r'\small', size=1.2)


# In[25]:


show(iota)


# In[26]:


s3 = s3p * iota
show(s3)


# ## Defining $\mathcal{T}_5$

# In[27]:


tiles5 = [('2030', '1', '2103', '0'),
 ('2033', '1', '2113', '0'),
 ('2103', '5', '2310', '0'),
 ('2113', '5', '2130', '1'),
 ('2113', '5', '2330', '0'),
 ('2130', '6', '2300', '0'),
 ('2130', '1', '2103', '5'),
 ('2133', '1', '2113', '1'),
 ('2300', '0', '2030', '6'),
 ('2300', '1', '2033', '6'),
 ('2310', '1', '2033', '6'),
 ('2330', '0', '2130', '6'),
 ('2330', '1', '2133', '6'),
 ('20330', '0', '21130', '0'),
 ('21030', '6', '23100', '0'),
 ('21030', '1', '21103', '1'),
 ('21033', '1', '21113', '1'),
 ('21103', '5', '21310', '1'),
 ('21113', '5', '21330', '1'),
 ('21130', '6', '21300', '1'),
 ('21130', '6', '23300', '0'),
 ('21300', '0', '21030', '5'),
 ('21300', '1', '21033', '5'),
 ('21310', '0', '21030', '5'),
 ('21310', '1', '21033', '5'),
 ('21330', '0', '21130', '1'),
 ('21330', '0', '21130', '5'),
 ('23100', '0', '20330', '6'),
 ('23300', '0', '21330', '6')]
T5 = WangTileSet(tiles5)
T5.tikz(font=r'\small', size=1.2)


# ## The morphism $\pi:\mathcal{T}_5\to\mathcal{T}_4$

# In[28]:


def project_56(t, p56='10'):
    L = []
    d = dict(zip('56', p56))
    for w in t:
        for a in '56':
            w = w.replace(a, d[a])
        L.append(w)
    return tuple(L)


# In[29]:


T4_tiles_list = T4.tiles()
d = {}
for i,t in enumerate(T5):
    pt = project_56(t, p56='10')
    j = T4_tiles_list.index(pt)
    d[i] = j
s4_pi = Substitution2d.from_permutation(d)
show(s4_pi)


# ## Computing $\mathcal{T}_6$ and the shear topological conjugacy $\eta:\Omega_6\to\Omega_5$

# In[30]:


get_ipython().run_line_magic('time', 'T6,s5_sheer = T5.shear(radius=2) # 20s with dancing_links, 22s with Gurobi')
T6.tikz(font=r'\small', size=1.2)


# In[31]:


show(s5_sheer)


# ## Computing $\mathcal{T}_7$

# In[32]:


get_ipython().run_line_magic('time', 'T6.find_markers(i=1, radius=1, solver=solver) #14min 31s with dancing_links, 17s with Gurobi')


# In[33]:


M6 = [6, 8, 9, 10, 11, 12, 21, 22, 23, 24, 26, 27, 28]
T7,s6,_ = T6.find_substitution(M=M6, i=1, radius=1, side='left', solver=solver)


# In[34]:


T7.tikz(font=r'\small', size=1.2)


# In[35]:


show(s6)


# ## Computing $\mathcal{T}_8$

# In[36]:


get_ipython().run_line_magic('time', 'T7.find_markers(i=1, radius=1, solver=solver) #4min 44s with dancing_links, 6s with Gurobi')


# In[37]:


M7 = [0, 1, 2, 9, 10, 11, 12]
T8,s7,_ = T7.find_substitution(M=M7, i=1, radius=1, solver=solver)


# In[38]:


T8.tikz(font=r'\small', size=1.2)


# In[39]:


show(s7)


# ## Computing $\mathcal{T}_9$

# In[40]:


get_ipython().run_line_magic('time', 'T8.find_markers(i=2, radius=2, solver=solver) # 5min 9s with dancing_links, 15s with Gurobi')


# In[41]:


M8 = [0, 1, 2, 3, 4, 5, 6, 7]
T9,s8,_ = T8.find_substitution(M=M8, i=2, radius=2, solver=solver)


# In[42]:


T9.table()


# In[43]:


show(s8)


# ## Computing $\mathcal{T}_{10}$

# In[44]:


get_ipython().run_line_magic('time', 'T9.find_markers(i=1, radius=1, solver=solver) # 6min 23s with dancing_links, 8s with Gurobi')


# In[45]:


M9 = [0, 1, 2, 3, 9, 10, 11, 12, 13]
T10,s9,_ = T9.find_substitution(M=M9, i=1, radius=1, solver=solver)


# In[46]:


T10.table()


# In[47]:


show(s9)


# ## Computing $\mathcal{T}_{11}$

# In[48]:


get_ipython().run_line_magic('time', 'T10.find_markers(i=2, radius=2, solver=solver) # 3min 37s with dancing_links, 11s with Gurobi')


# In[49]:


M10 = [0, 1, 2, 3, 4, 5]
T11,s10,_ = T10.find_substitution(M=M10, i=2, radius=2, solver=solver)


# In[50]:


T11.table()


# In[51]:


show(s10)


# ## Computing $\mathcal{T}_{12}$

# In[52]:


get_ipython().run_line_magic('time', 'T11.find_markers(i=1, radius=1, solver=solver) # 5min 25s with dancing_links, 7s with Gurobi')


# In[53]:


M11 = [0, 1, 2, 9, 10, 11]
T12,s11,_ = T11.find_substitution(M=M11, i=1, radius=1, solver=solver)


# In[54]:


T12.table()


# In[55]:


show(s11)


# ## Definition of $\mathcal{U}$

# In[56]:


tilesU = ['FOJO','FOHL','JMFP','DMFK','HPJP','HPHN','HKFP','HKDP','BOIO','GLEO','GLCL',
'ALIO','EPGP','EPIP','IPGK','IPIK','IKBM','IKAK','CNIP',]
tilesU = map(tuple, tilesU)
U = WangTileSet(tilesU)
U.tikz(size=1.2)


# ## Equivalence of $\mathcal{U}$ and $\mathcal{T}_{12}$

# In[57]:


U.is_equivalent(T12, certificate=True)


# ## Computing $\mathcal{T}_{13}$

# In[58]:


get_ipython().run_line_magic('time', 'U.find_markers(i=2, radius=2, solver=solver) # 4min 39s with dancing_links, 13s with Gurobi')


# In[59]:


M12 = [0, 1, 2, 3, 4, 5, 6, 7]
T13,s12,_ = U.find_substitution(M=M12, i=2, radius=2, solver=solver)


# In[60]:


T13.tikz(size=1.2)


# In[61]:


show(s12)


# ## Computing $\mathcal{T}_{14}$

# In[62]:


get_ipython().run_line_magic('time', 'T13.find_markers(i=1, radius=1, solver=solver) # 5min 34s with dancing_links, 7s with Gurobi')


# In[63]:


M13 = [0, 1, 3, 8, 9, 14, 15]
T14,s13,_ = T13.find_substitution(M=M13, i=1, radius=1, solver=solver)


# In[64]:


T14.tikz(size=1.2)


# In[65]:


show(s13)


# ## The self-similarity $\omega_{12}:\Omega_{12}\to\Omega_{12}$

# In[66]:


is_equiv,f,g,permUtoT14 = U.is_equivalent(T14, certificate=True)
assert is_equiv
s12to12 = s12 * s13 * permUtoT14


# In[67]:


show(s12to12)


# In[68]:


s12to12.wang_tikz(U, U, font=r'\scriptsize', scale=1, codomain_color=color, ncolumns=4, 
            direction='right', extra_space=1.5)


# ## The morphism $\omega_0\omega_1\omega_2\omega_3:\Omega_{4}\to\Omega_{0}$

# In[69]:


s0to4 = s0*s1*s2*s3p*iota
s0to4.wang_tikz(T4, T0, font=r'\scriptsize', scale=.8, codomain_color=color, ncolumns=10, direction='down',
            extra_space=1.5)


# ## The morphism $\eta\omega_6\omega_7\omega_8\omega_9\omega_{10}\omega_{11}:\Omega_{12}\to\Omega_{5}$

# In[70]:


s5to12 = s5_sheer*s6*s7*s8*s9*s10*s11
M = matrix(2, [1,1,0,1])
s5to12sheer = s5to12.apply_matrix_transformation(M)
s5to12sheer.wang_tikz(U, T5, font=r'\scriptsize', scale=.9, ncolumns=2, extra_space=1.5)


# ## The morphism $\pi\eta\omega_6\omega_7\omega_8\omega_9\omega_{10}\omega_{11}:\Omega_{12}\to\Omega_{4}$

# In[71]:


s4to12 = s4_pi*s5to12
M = matrix(2, [1,1,0,1])
s4to12sheer = s4to12.apply_matrix_transformation(M)
s4to12sheer.wang_tikz(U, T4, font=r'\scriptsize', scale=.9, ncolumns=2, extra_space=1.5)


# ## The morphism $\omega_0\omega_1\omega_2\omega_3\pi\eta\omega_6\omega_7\omega_8\omega_9\omega_{10}\omega_{11}:\Omega_{12}\to\Omega_{0}$

# In[72]:


s0to12 = s0to4*s4to12sheer
s0to12.wang_tikz(U, T0, font=r'\scriptsize', scale=.8, codomain_color=color, id=None, 
                 ncolumns=4, direction='right', extra_space=1.5)


# ## Computing the eight fixed point of $\omega_{12}:\Omega_{12}\to\Omega_{12}$

# In[73]:


from slabbe import TikzPicture
G = s12to12.prolongable_origins()
TikzPicture.from_graph(G, rankdir='right')


# ## Horizontal Rauzy graph for the tiles $\mathcal{T}_{12}$

# In[74]:


Gh = DiGraph(T12.dominoes_with_surrounding(i=1,radius=2), format='list_of_edges')
TikzPicture.from_graph(Gh)


# ## Vertical Rauzy graph for the tiles $\mathcal{T}_{12}$

# In[75]:


Gv = DiGraph(T12.dominoes_with_surrounding(i=2,radius=2), format='list_of_edges')
TikzPicture.from_graph(Gv)


# ## Frequencies of tiles

# Inverse of frequencies of tiles in $\Omega_4$:

# In[76]:


M = matrix(s12to12)
z = polygen(QQ, 'z')
K.<phi> = NumberField(z**2-z-1, 'phi', embedding=RR(1.6))
MK = M.change_ring(K)
v = MK.eigenvectors_right()[0][1][0]
v4 = matrix(s4to12) * v
sorted((1/a, i) for (i,a) in enumerate(v4/sum(v4)))


# Inverse of frequencies of tiles in $\Omega_0$:

# In[77]:


v0 = matrix(s0to4)*matrix(s4to12) * v
sorted((1/a, i) for (i,a) in enumerate(v0/sum(v0)))