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

# $$
# \def\CC{\bf C}
# \def\QQ{\bf Q}
# \def\RR{\bf R}
# \def\ZZ{\bf Z}
# \def\NN{\bf N}
# $$
# # Rauzy induction of polygon partitions and toral $\mathbb{Z}^2$-rotations
# 
# This file contains the Sage code contained in the preprints:
# 
# > -   [arXiv:1906.01104v1](https://arxiv.org/abs/1906.01104v1), June 2019, 36 p.
# > -   [arXiv:1906.01104v2](https://arxiv.org/abs/1906.01104v2), May 2020, 40 p.
# > -   [arXiv:1906.01104v3](https://arxiv.org/abs/1906.01104v3), January 2021, 40 p.
# 
# The Jupyter notebook `arXiv_1906_01104.ipynb` is created from `arXiv_1906_01104.rst` with the command `sage -rst2ipynb arXiv_1906_01104.rst arXiv_1906_01104.ipynb`. The file `arXiv_1906_01104.rst` is in the `slabbe/demos` directory of the optional SageMath package [slabbe](https://pypi.org/project/slabbe/).
# 
# Running all examples below with `sage -t arXiv_1906_01104.rst` takes 10 seconds with sage-9.1 and slabbe-0.6.1 (or with sage-9.2 and slabbe-0.6.2).

# In[1]:


get_ipython().run_line_magic('display', 'latex           # not tested')


# First we construct the golden mean as a element of a quadratic number field because it is more efficient for arithmetic operations and comparisons:

# In[2]:


z = polygen(QQ, 'z')
K = NumberField(z**2-z-1, 'phi', embedding=RR(1.6))
phi = K.gen()


# We import the polygon partition $\mathcal{P}_0$ of $\mathbb{R}^2/\Gamma_0$ : which is predefined in `slabbe` :

# In[3]:


from slabbe.arXiv_1903_06137 import jeandel_rao_wang_shift_partition
P0 = jeandel_rao_wang_shift_partition()
P0


# In[4]:


P0.plot()                    # optional long


# We import polyhedron exchange transformations from the package:

# In[5]:


from slabbe import PolyhedronExchangeTransformation as PET


# We define the lattice $\Gamma_0$ and the maps $R_0^{e_1}$, $R_0^{e_2}$ which can be seen as a polygon exchange transformations on a rectangular fundamental domain of $\mathbb{R}^2/\Gamma_0$ :

# In[6]:


Gamma0 = matrix.column([(phi,0), (1,phi+3)])
fundamental_domain = polytopes.parallelotope([(phi,0), (0,phi+3)])
R0e1 = PET.toral_translation(Gamma0, vector((1,0)), fundamental_domain)
R0e2 = PET.toral_translation(Gamma0, vector((0,1)), fundamental_domain)


# The following allows to compute the induced partition $\mathcal{P}_1$ of $\mathbb{R}^2/\Gamma_1$, the substitution $\beta_0$ and the $\mathbb{Z}^2$-action $R_1$ on $\mathbb{R}^2/\Gamma_1$.

# In[7]:


y_le_1 = [1, 0, -1]   # syntax for the inequality y <= 1
P1,beta0 = R0e2.induced_partition(y_le_1, P0, substitution_type='column')
R1e1,_ = R0e1.induced_transformation(y_le_1)
R1e2,_ = R0e2.induced_transformation(y_le_1)


# In[8]:


R1e1


# In[9]:


R1e1.plot()              # optional long


# In[10]:


R1e2


# In[11]:


R1e2.plot()              # optional long


# In[12]:


P1


# In[13]:


P1.plot()                    # optional long


# In[14]:


beta0


# We keep $\mathcal{P}_2$, $\Gamma_2$ equal to $\mathcal{P}_1$, $\Gamma_1$ and we change the base of the action to get the $\mathbb{Z}^2$-action $R_2$ on $\mathbb{R}^2/\Gamma_2$.

# In[15]:


Gamma2 = Gamma1 = matrix.column([(phi,0), (0,1)])
P2 = P1
R2e1 = R1e1
R2e2 = (R1e1 * R1e2).merge_atoms_with_same_translation()


# The following confirms that the $R_2^{\boldsymbol{e}_2}$ is now a vertical rotation on the torus $\mathbb{R^2}/\Gamma_2$.

# In[16]:


R2e2


# In[17]:


R2e2.plot()              # optional long


# The following allows to compute the induced partition $\mathcal{P}_3$ of $\mathbb{R}^2/\Gamma_3$, the substitution $\beta_2$ and the $\mathbb{Z}^2$-action $R_3$ on $\mathbb{R}^2/\Gamma_3$.

# In[18]:


x_le_1 = [1, -1, 0]   # syntax for x <= 1
P3,beta2 = R2e1.induced_partition(x_le_1, P2, substitution_type='row')
R3e1,_ = R2e1.induced_transformation(x_le_1)
R3e2,_ = R2e2.induced_transformation(x_le_1)
R3e1


# In[19]:


P3


# In[20]:


P3.plot()                    # optional long


# In[21]:


beta2


# The following allows to compute the induced partition $\mathcal{P}_4$ of $\mathbb{R}^2/\Gamma_4$, the substitution $\beta_3$ and the $\mathbb{Z}^2$-action $R_4$ on $\mathbb{R}^2/\Gamma_4$

# In[22]:


x_le_phi_inv = [phi^-1, -1, 0]    # syntax for x <= phi^-1
P4,beta3 = R3e1.induced_partition(x_le_phi_inv, P3, substitution_type='row')
R4e1,_ = R3e1.induced_transformation(x_le_phi_inv)
R4e2,_ = R3e2.induced_transformation(x_le_phi_inv)
R4e2


# In[23]:


P4


# In[24]:


P4.plot()                    # optional long


# In[25]:


beta3


# The following allows to compute the induced partition $\mathcal{P}_5$ of $\mathbb{R}^2/\Gamma_5$, the substitution $\beta_4$ and the $\mathbb{Z}^2$-action $R_5$ on $\mathbb{R}^2/\Gamma_5$.

# In[26]:


y_le_phi_inv = [phi^-1, 0, -1]    # syntax for y <= phi^-1
P5,beta4 = R4e2.induced_partition(y_le_phi_inv, P4, substitution_type='column')
R5e1,_ = R4e1.induced_transformation(y_le_phi_inv)
R5e2,_ = R4e2.induced_transformation(y_le_phi_inv)
P5


# In[27]:


P5.plot()                    # optional long


# In[28]:


beta4


# We rescale the partition $\mathcal{P}_5$ :

# In[29]:


P5_scaled = (-phi*P5).translate((1,1))
R5e1_scaled = (-phi*R5e1).translate_domain((1,1))
R5e2_scaled = (-phi*R5e2).translate_domain((1,1))


# In[30]:


P5_scaled.plot()                  # optional long


# The following allows to compute the induced partition $\mathcal{P}_6$ of $\mathbb{R}^2/\Gamma_6$, the substitution $\beta_5$ and the $\mathbb{Z}^2$-action $R_6$ on $\mathbb{R}^2/\Gamma_6$.

# In[31]:


P6,beta5 = R5e1_scaled.induced_partition(x_le_phi_inv, P5_scaled, substitution_type='row')
R6e1,_ = R5e1_scaled.induced_transformation(x_le_phi_inv)
R6e2,_ = R5e2_scaled.induced_transformation(x_le_phi_inv)
P6


# In[32]:


P6.plot()                    # optional long


# In[33]:


beta5


# The following allows to compute the induced partition $\mathcal{P}_7$ of $\mathbb{R}^2/\Gamma_7$, the substitution $\beta_6$ and the $\mathbb{Z}^2$-action $R_7$ on $\mathbb{R}^2/\Gamma_7$.

# In[34]:


P7,beta6 = R6e2.induced_partition(y_le_phi_inv, P6, substitution_type='column')
R7e1,_ = R6e1.induced_transformation(y_le_phi_inv)
R7e2,_ = R6e2.induced_transformation(y_le_phi_inv)
P7


# In[35]:


P7.plot()                    # optional long


# In[36]:


beta6


# We rescale the partition $\mathcal{P}_7$ :

# In[37]:


P7_scaled = (-phi*P7).translate((1,1))
R7e1_scaled = (-phi*R7e1).translate_domain((1,1))
R7e2_scaled = (-phi*R7e2).translate_domain((1,1))


# In[38]:


P7_scaled.plot()                  # optional long


# The following allows to compute the induced partition $\mathcal{P}_8$ of $\mathbb{R}^2/\Gamma_8$, the substitution $\beta_7$ and the $\mathbb{Z}^2$-action $R_8$ on $\mathbb{R}^2/\Gamma_8$.

# In[39]:


P8,beta7 = R7e1_scaled.induced_partition(x_le_phi_inv, P7_scaled, substitution_type='row')
R8e1,_ = R7e1_scaled.induced_transformation(x_le_phi_inv)
R8e2,_ = R7e2_scaled.induced_transformation(x_le_phi_inv)
P8


# In[40]:


P8.plot()                    # optional long


# In[41]:


beta7


# The following allows to compute the induced partition $\mathcal{P}_9$ of $\mathbb{R}^2/\Gamma_9$, the substitution $\beta_8$ and the $\mathbb{Z}^2$-action $R_9$ on $\mathbb{R}^2/\Gamma_9$.

# In[42]:


P9,beta8 = R8e2.induced_partition(y_le_phi_inv, P8, substitution_type='column')
R9e1,_ = R8e1.induced_transformation(y_le_phi_inv)
R9e2,_ = R8e2.induced_transformation(y_le_phi_inv)
P9


# In[43]:


P9.plot()                    # optional long


# In[44]:


beta8


# We rescale the partition $\mathcal{P}_9$ :

# In[45]:


P9_scaled = (-phi*P9).translate((1,1))
R9e1_scaled = (-phi*R9e1).translate_domain((1,1))
R9e2_scaled = (-phi*R9e2).translate_domain((1,1))


# In[46]:


P9_scaled.plot()                  # optional long


# The following allows to compute the induced partition $\mathcal{P}_{10}$ of $\mathbb{R}^2/\Gamma_{10}$, the substitution $\beta_9$ and the $\mathbb{Z}^2$-action $R_{10}$ on $\mathbb{R}^2/\Gamma_{10}$.

# In[47]:


P10,beta9 = R9e1_scaled.induced_partition(x_le_phi_inv, P9_scaled, substitution_type='row')
R10e1,_ = R9e1_scaled.induced_transformation(x_le_phi_inv)
R10e2,_ = R9e2_scaled.induced_transformation(x_le_phi_inv)
P10


# In[48]:


P10.plot()                    # optional long


# In[49]:


beta9


# We show that $\mathcal{P}_8$ and $\mathcal{P}_{10}$ are equivalent:

# In[50]:


P8.is_equal_up_to_relabeling(P10)


# In[51]:


from slabbe import Substitution2d
tau = Substitution2d.from_permutation(P8.keys_permutation(P10))
tau


# In[52]:


beta8*beta9*tau    # the self-similarity for P8


# We may check that the self-similarity for $\mathcal{P}_8$ satisfies $\zeta^{-1}\beta_8\beta_9\tau\zeta=\beta_{\mathcal{U}}$.

# In[53]:


zeta = Substitution2d.from_permutation({0:0, 1:1, 2:9, 3:7, 4:8, 5:11, 6:10, 
7:6, 8:2, 9:4, 10:5, 11:3, 12:18, 13:14, 14:16, 15:13, 16:12, 17:17, 18:15})
betaU = Substitution2d({0: [[17]], 1: [[16]], 2: [[15], [11]], 3: [[13], [9]], 4: [[17], [8]], 5: [[16], [8]], 6: [[15], [8]], 7: [[14], [8]], 8: [[14, 6]], 9: [[17, 3]], 10: [[16, 3]], 11: [[14, 2]], 12: [[15, 7], [11, 1]], 13: [[14, 6], [11, 1]], 14: [[13, 7], [9, 1]], 15: [[12, 6], [9, 1]], 16: [[18, 5], [10, 1]], 17: [[13, 4], [9, 1]], 18: [[14, 2], [8, 0]]})
zeta.inverse()*beta8*beta9*tau*zeta  == betaU


# In[54]:


betaU


# Observe that the $19 \times 19$ incidence matrix of $\beta_8\beta_9\tau$ is not hyperbolic but, as shown by its characteristic polynomial, it is hyperbolic on a 8-dimensional subspace:

# In[55]:


(beta8*beta9*tau).incidence_matrix().charpoly().factor()