from qutip import *
import matplotlib.pyplot as plt
import numpy as np

val_s = ['site0','site1']
(H_cell_form,T_inter_cell_form,H_cell,T_inter_cell) = cell_structures( val_s)

H_cell_form

T_inter_cell_form

"""
Explanation: The Su-Schrieffer–Heeger (SSH) model
Saumya Biswas (saumyab@uoregon.edu)
The celebrated SSH model is analyzed with QuTiP's lattice module below.

The above figure shows a SSH model with 6 sites with periodic boundary condition. The same lattice with hardwall/aperiodic boundary condition would be the folloowing.

In the secod quantized formalism, the periodic lattice Hamiltonian can be written as
\begin{eqnarray}
H_{per} = -t_{intra} (c_{-2}^{\dagger} c_{-1} + c_{-1}^{\dagger} c_{-2} ) -t_{intra} (c_{0}^{\dagger} c_{1} + c_{1}^{\dagger} c_{0} ) -t_{intra} (c_{2}^{\dagger} c_{3} + c_{3}^{\dagger} c_{2} ) \nonumber \ 
 -t_{inter} (c_{-1}^{\dagger} c_{0} + c_{0}^{\dagger} c_{-1} ) -t_{inter} (c_{1}^{\dagger} c_{2} + c_{2}^{\dagger} c_{1} ) -t_{inter} (c_{3}^{\dagger} c_{-2} + c_{-2}^{\dagger} c_{3} ) \nonumber 
\end{eqnarray}
The aperiodic lattice Hamiltonian can be obtained by discarding the very last term.
\begin{eqnarray}
H_{aper} = -t_{intra} (c_{-2}^{\dagger} c_{-1} + c_{-1}^{\dagger} c_{-2} ) -t_{intra} (c_{0}^{\dagger} c_{1} + c_{1}^{\dagger} c_{0} ) -t_{intra} (c_{2}^{\dagger} c_{3} + c_{3}^{\dagger} c_{2} ) \nonumber \ 
 -t_{inter} (c_{-1}^{\dagger} c_{0} + c_{0}^{\dagger} c_{-1} ) -t_{inter} (c_{1}^{\dagger} c_{2} + c_{2}^{\dagger} c_{1} ) \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \nonumber 
\end{eqnarray}
The representation in terms of unit cell blocks become obvious once we resolve the terms into unit cell operators.
\begin{eqnarray}
H_{per}= \begin{bmatrix}
    c_{-2}^{\dagger} & c_{-1}^{\dagger}     
\end{bmatrix} 
\begin{bmatrix}
      o     & -t_1 \
    -t_1    & 0 
\end{bmatrix}
\begin{bmatrix}
    c_{-2}  \
    c_{-1}
\end{bmatrix} + 
\begin{bmatrix}
    c_{0}^{\dagger} & c_{1}^{\dagger}     
\end{bmatrix} 
\begin{bmatrix}
      o     & -t_1 \
    -t_1    & 0 
\end{bmatrix}
\begin{bmatrix}
    c_{0}  \
    c_{1}
\end{bmatrix} 
\nonumber \
+ \begin{bmatrix}
    c_{2}^{\dagger} & c_{3}^{\dagger}     
\end{bmatrix} 
\begin{bmatrix}
      o     & -t_1 \
    -t_1    & 0 
\end{bmatrix}
\begin{bmatrix}
    c_{2}  \
    c_{3}
\end{bmatrix} \nonumber \
+ \begin{bmatrix}
    c_{-2}^{\dagger} & c_{-1}^{\dagger}     
\end{bmatrix} 
\begin{bmatrix}
      o     & 0 \
    -t_2    & 0 
\end{bmatrix}
\begin{bmatrix}
    c_{0}  \
    c_{1}
\end{bmatrix} + 
\begin{bmatrix}
    c_{0}^{\dagger} & c_{1}^{\dagger}     
\end{bmatrix} 
\begin{bmatrix}
      o     & -t_2 \
      0     & 0 
\end{bmatrix}
\begin{bmatrix}
    c_{-2}  \
    c_{-1}
\end{bmatrix} 
\nonumber \
+ \begin{bmatrix}
    c_{0}^{\dagger} & c_{1}^{\dagger}     
\end{bmatrix} 
\begin{bmatrix}
      o     & 0 \
    -t_2    & 0 
\end{bmatrix}
\begin{bmatrix}
    c_{2}  \
    c_{3}
\end{bmatrix} + 
\begin{bmatrix}
    c_{2}^{\dagger} & c_{3}^{\dagger}     
\end{bmatrix} 
\begin{bmatrix}
      o     & -t_2 \
      0     & 0 
\end{bmatrix}
\begin{bmatrix}
    c_{0}  \
    c_{1}
\end{bmatrix} 
\nonumber \
+ \begin{bmatrix}
    c_{2}^{\dagger} & c_{3}^{\dagger}     
\end{bmatrix} 
\begin{bmatrix}
      o     & 0 \
    -t_2    & 0 
\end{bmatrix}
\begin{bmatrix}
    c_{-2}  \
    c_{-1}
\end{bmatrix} + 
\begin{bmatrix}
    c_{-2}^{\dagger} & c_{-1}^{\dagger}     
\end{bmatrix} 
\begin{bmatrix}
      o     & -t_2 \
      0     & 0 
\end{bmatrix}
\begin{bmatrix}
    c_{2}  \
    c_{3}
\end{bmatrix} \nonumber
\end{eqnarray}
Hence, $H_{TB}$ can be succinctly put in the form:
\begin{eqnarray}
H_{per/aper} = \sum_i \psi_i^{\dagger} D \psi_i + \sum_{i} \left( \psi_i^{\dagger} T \psi_{i+1} +  \psi_{i+1}^{\dagger} T^{\dagger} \psi_i  \right) \label{eq:TB_block} 
\end{eqnarray}
where $D = \begin{bmatrix}
    0 & -t_{intra}  \
    -t_{intra} & 0
\end{bmatrix}$ and $T = \begin{bmatrix}
    0 & 0  \
    -t_{inter} & 0
\end{bmatrix}$. The $\psi_i = \begin{bmatrix}
    c^0_{i}  \
    c^1_{i}
\end{bmatrix}$ is associated with $\vec{R}=\vec{R}i$ and $\psi{i+1}$ is associated with $\vec{R}=\vec{R}_i + \hat{x}$, with $\hat{x}$ being the unit vector along the x direction.
The equation above can be put into the alternate form (also changing the summation variables to m,n to distinguish them from the imaginary i):
\begin{eqnarray}
H_{per/aper} = \sum_{m,n} \psi_m^{\dagger} (D \delta_{m,n} + T \delta_{m,n-1} + T^{\dagger} \delta_{m,n+1} ) \psi_n \label{eq:TB_block_1} 
\end{eqnarray}
So far, we have see that a SSH lattice can be put into the unit cell representations and the lattice can be completely defied with two matrices D and T. In declaring an instance of Qutip.lattice.Lattice1d we only need to input the two matrices. 
In declaring an instance of Lattice1d class, the two arguments cell_Hamiltonian
and inter_hop are set to D and T respectively. The matrix structure of D and T can be
obtained from the aide function cell_structures().
End of explanation
"""


t_intra = -0.5
t_inter = -0.5
H_cell[0,1] = t_intra
H_cell[1,0] = t_intra
T_inter_cell[1,0] = t_inter
H_cell = Qobj(H_cell)
T_inter_cell = Qobj(T_inter_cell)

"""
Explanation: Guided by cell_H_form and inter_cell_T_form, we can set values to cell_H and inter_cell_T which were initialized to all zero elements by cell_structures().
End of explanation
"""


H_cell = Qobj( np.array( [[ 0, t_intra ],[t_intra,0]] )  )
T_inter_cell =  Qobj( np.array( [[ 0, 0 ],[t_inter,0]] )  ) 

"""
Explanation: Using cell_structures() is completely optional. The user could equally well have defined cell_H and inter_cell_T directly.
End of explanation
"""


boundary_condition = "periodic"
cells = 3
cell_sites = 2
site_dof = [1]
SSH_lattice = Lattice1d(num_cell=cells, boundary = boundary_condition, cell_num_site = cell_sites, cell_site_dof = site_dof, Hamiltonian_of_cell = H_cell, inter_hop = T_inter_cell )

"""
Explanation: For our SSH lattice with 3 unit cells, 2 sites in each unit cell, and [1] degree of freedom per each site, we can initiate an instance of the Lattice1d class at this stage.
End of explanation
"""


H = SSH_lattice.display_unit_cell(label_on = True)
T = SSH_lattice.display_lattice()

"""
Explanation: The model can be visualized with the display functions.
End of explanation
"""


SSH_Haml = SSH_lattice.Hamiltonian()
SSH_Haml

"""
Explanation: The Hamiltonian of the lattice can be obtained with the method Hamiltonian()
End of explanation
"""


chiral_op = SSH_lattice.distribute_operator(sigmaz())        
anti_commutator_chi_H = chiral_op * SSH_Haml + SSH_Haml * chiral_op
is_null = (np.abs(anti_commutator_chi_H) < 1E-10).all()
print(is_null)

"""
Explanation: Sublattice Projectors and Chiral Symmetry of the SSH model:
We explain now that, the Hamiltonian of the SSH model is bipartite. We can set the sites into two partitions. Each partition consists of the set of every other sites starting from site 0(or 1). Since, the Hamiltonian  does not enable any transition between sites of the same sublattice, only between them, the Hamiltonian is said to be bipartite.
\begin{eqnarray}
\hat{P}{A} = \sum\limits{m=1}^{N} |m,A\rangle \langle n, A|,\ \ \ \ \ \ \hat{P}{B} = \sum\limits{m=1}^{N} |m,B\rangle \langle n, B|, \ \ \ \ \ \
\hat{\Sigma}z = \hat{P}{A} - \hat{P}_{B}
\end{eqnarray}
The projectors $\hat{P}A$(or $\hat{P}_B$) project out a ket vector into its projection into sublattice A(or B). Because, the Hamiltonian is bipartite,
\begin{eqnarray}
\hat{P}{A} H \hat{P}{A} = \hat{P}{B} H \hat{P}_{B} = 0
\end{eqnarray}
The SSH Hamiltonian is said to have chiral symmetry, since
\begin{eqnarray}
\hat{\Sigma}_z H \hat{\Sigma}_z = -H
\end{eqnarray}
which is a generalization of the symmetry concept of a Hamiltonian, since, ordinarily, a Hamiltonian is said to have a symmetry if an unitary operator leaves it invariant.
\begin{eqnarray}
\hat{U} H \hat{U}^{\dagger} = H
\end{eqnarray}
Now, we form the operators $\hat{P}_A$,$\hat{P}_B$ and $\hat{\Sigma}_z$, and verify the properties of the SSH Hamiltonian.
End of explanation
"""


SSH_lattice.plot_dispersion()

"""
Explanation: Hence, it is verified that $\hat{\Sigma}_z$ and H indeed anticommute and the SSH Hamiltonian has chiral symmetry.
The dispersion relationship for the lattice can be obtained with plot_dispersion() method.
End of explanation
"""


[V,[S0,S1,S2,S3,S4,S5]]=SSH_Haml.eigenstates()

V

"""
Explanation: plot_dispersion() plots the 3(since number of unit cells in 3) points in k-space 
(the first and last one are the same) over the dispersion relation of an infinite
crystal.
End of explanation
"""


print(S0)
print(chiral_op*S0)
print(S5)

"""
Explanation: First, the eigen-values are the same as the ones obtained from the dispersion calculation. Second, they are symmetric about the value 0.
The second is a consequence of the chiral symmetry of the Hamiltonian, as we explain now.
\begin{eqnarray}
\hat{\bf{1}} = \hat{P}_A + \hat{P}_B, \ \ \ \ \, \hat{P}_A = \frac{1}{2}(\hat{\bf{1}}+\hat{\Sigma}_z), \ \ \ \ \, \hat{P}_B = \frac{1}{2}(\hat{\bf{1}}-\hat{\Sigma}_z) \nonumber \
H |\psi_n\rangle = E_n | \psi_n \rangle \implies H (\hat{\Sigma}_z | \psi_n \rangle) = - \hat{\Sigma}_z H | \psi_n \rangle = -E_n (\hat{\Sigma}_z|\psi_n\rangle)
\end{eqnarray}
So, if $|\psi_n\rangle$ is an eigenstate with eigenenergy $E_n$, $\hat{\Sigma}_z | \psi_n \rangle$ is also an eigenstate with energy $-E_n$ and the eigen-spectrum is symmetric about 0.
Here, S0 is the eigenvector with eigenvalue -1 and S5 is the eigenvector with eigenvalue +1. So, we can verify if S5 is the same eigenvector(withinn a phase factor) as ($\hat{\Sigma}_z*$S0).
End of explanation
"""


is_null = (np.abs(S0.dag()*S5) < 1E-10).all()
print(is_null) # Are S0 and S5 orthogonal?

dimH = chiral_op.dims
identity_H = Qobj(np.identity(6), dims=dimH)
identity_H
PA = 0.5*( identity_H + chiral_op)
PB = 0.5*( identity_H - chiral_op)
support_S5_A = S5.dag() * PA * S5
support_S5_B = S5.dag() * PB * S5
print(support_S5_A == support_S5_B) # Does S5 have equal support on A and B?

"""
Explanation: Clearly, S5 is the same eigenvector as ($\hat{\Sigma}_z*$S0).
Since, $\hat{\Sigma}_z | \psi_n \rangle$ and $| \psi_n \rangle$ are eigenvectors of a Hermitian oerator,H with distinct eigenvalues, they must be orthogonal.
\begin{eqnarray}
E_n \ne 0 \implies 0 = \langle \phi_n| \hat{\Sigma}_z | \phi_n \rangle = \langle \phi_n \hat{P}_A |\phi_n \rangle - \langle \phi_n| \hat{P}_B |\phi_n \rangle
\end{eqnarray}
i.e., an eigenstate with nonzero eigenvalue has equal support on both sublattices. We, now check this for S5.
End of explanation
"""


(knxA, qH_ks) = SSH_lattice.bulk_Hamiltonians()
qH_ks

"""
Explanation: We discuss the implications of $E_n = 0$ for an eigenvector later in the context of edge states later.
Unsurprisingly, diagonalizing the Hamiltonian gives the same spectrum of eigen-values
as the one obtained from the plot_dispersion() function. We shall soon illustrate, translational symmetry is a very useful computational hack. Here, we see how they can produce the eigen-values and eigen-vectors of the Hamiltonian ($6\times6$, in our example) from diagonalizing a $2\times2$ matrix. The reduction in size by a factor of 3 comes from the fact that the lattice of 3 cells repeats itself infinitely on both ends.
Using Translational Symmetry:
Any periodic lattice Hamiltonian can be diagonalized easier exploiting translational symmetry, this feature is not specific to SSH model alone.
\begin{eqnarray}
|\psi_n(k) \rangle = |k \rangle \otimes | u_{n}(k) \rangle   \nonumber \
| u_{n}(k) \rangle = a_n(k)|a\rangle  + b_n(k)|b\rangle \nonumber \
\end{eqnarray}
The vectors $| u_{nk}(r) \rangle \in H_{internal}$ are the eigenstates of the bulk momentum space Hamiltonian $H(k)$ defined as
\begin{eqnarray}
\langle k | H_{bulk} | k \rangle = \sum\limits_{\alpha,\beta \in {A,B}} \langle k, \alpha | H_{bulk} | k, \beta \rangle | \alpha \rangle \langle \beta |  \nonumber \
H(k)|u_{n}(k) \rangle = E(k)| u_{n}(k) \rangle
\end{eqnarray}
In a lattice with N unitcells, since $|\psi_n(k) \rangle$ is required to be periodic over N cells, it needs to be invariant with a translation by N cells, and the valid/good quantum number for k are $0,1,2,...,\frac{2\pi}{N}$.
The dispersion, $E(k)$ can be obtained from get_dispersion() method. A list of $H(k)$ at the valid quantum numbers can be produced from the bulk_Hamiltonian_array() method.
End of explanation
"""


(knxA, vec_kns) = SSH_lattice.cell_periodic_parts()
vec_kns

"""
Explanation: The array of $|u_n(k) \rangle$(in terms of its expansion in {a(k),b(k)}) at the good quantum numbers, k can be produced with the method array_of_unk().
End of explanation
"""


knxA

"""
Explanation: In both cases, knxA is simply an array containing the valid values of k, in units of $2\pi/a$, a being the length of the unit cell i.e. 1. a is always 1 in all methods().
End of explanation
"""


eigen_states = SSH_lattice.bloch_wave_functions()
eigen_states

"""
Explanation: bloch_wave_functions() yields an ordered array for the eigenvalues and eigenvectors(which are bloch wave functions) for the Hailtonian of the lattice.
End of explanation
"""


t_intra = -0.5
t_inter = -0.35
H_cell = Qobj( np.array( [[ 0, t_intra ],[t_intra,0]] )  )
T_inter_cell =  Qobj( np.array( [[ 0, 0 ],[t_inter,0]] )  ) 
SSH_lattice_TrI = Lattice1d(num_cell=100, boundary = "periodic", cell_num_site = 2, cell_site_dof = [1], Hamiltonian_of_cell = H_cell, inter_hop = T_inter_cell )

SSH_lattice_TrI.plot_dispersion()

SSH_H_t = SSH_lattice_TrI.Hamiltonian()
D = SSH_H_t.eigenenergies()

plt.plot(D,'ro')
plt.xlabel('index of eigen values')
plt.ylabel('eigen values')
plt.show()
plt.close()

"""
Explanation: Knowing the cell periodic part of a bloch wavefunction, $u_n(k)$ suffices to calculate it. The translational symmetry enables calculation of the eigenstates of a $6\times6$ matrix through diagonalizing a $2\times2$ matrix.
Topology: Winding Number:
Due to the chiral symmetry, the bulk momentum space Hamiltonian,$H(k)$ can be written in terms of $\sigma_x$ and $\sigma_y$ components alone.
\begin{eqnarray}
H(k)= h_x(k)\sigma_x + h_y(k)\sigma_y
\end{eqnarray}
For this specific model, where $\bf{h}(k)$ moves about in a 2d plane($h_x-h_y$), winding number is a topological invariant characterizing the topology of the model. It enumerates the number of times, $\bf{h}(k)$ traverses around the origin in a positive sense as k is varied from 0 to $2\pi$.
The method winding_number() evaluates the following integral to determine it as well as plots the trajectory of $\bf{h}(k)$ in the $h_x-h_y$ plane.
\begin{eqnarray}
\nu = \frac{1}{2\pi i}\int\limits_{-\pi}^{\pi} dk \frac{d}{dk}Log(h(k)) \ \ \ \text{,where} \ \ \ \ \ \ \ \ h(k) = h_x(k)- ih_y(k)
\end{eqnarray}
The trivial Insulator:
The latice become a gapped system when $t_{inter}$ does not equal $t_{intra}$. For the case of $|t_{inter}| < |t_{intra}|$, the lattice is a topologically trivial insulator.
End of explanation
"""


SSH_lattice_TrI.winding_number()

"""
Explanation: Unsurprisingly, we see a gap in the dispersion relationship and also in the spectrum of
eigen values of the Hamiltonian.
However, a calculation of the winding number yields 0, since it has trivial topology.
End of explanation
"""


t_intra = -0.5
t_inter = -0.65
H_cell = Qobj( np.array( [[ 0, t_intra ],[t_intra,0]] )  )
T_inter_cell =  Qobj( np.array( [[ 0, 0 ],[t_inter,0]] )  ) 
pSSH_lattice_nTrI = Lattice1d(num_cell=100, boundary = "periodic", cell_num_site = 2, cell_site_dof = [1], Hamiltonian_of_cell = H_cell, inter_hop = T_inter_cell )

pSSH_lattice_nTrI.plot_dispersion()

pSSH_H_nt = pSSH_lattice_nTrI.Hamiltonian()
nD = pSSH_H_nt.eigenenergies()

plt.plot(nD,'ro')
plt.ylabel('eigen values')
plt.show()
plt.close()

"""
Explanation: The Topologically Nontrivial Insulator
The other limit ($|t_{inter}| < |t_{intra}| $) is the interesting one, although the same
calculations on a lattice with periodic boundary condition reveal nothing interesting.
End of explanation
"""


t_intra = -0.5
t_inter = -0.65
cell_H = Qobj( np.array( [[ 0, t_intra ],[t_intra,0]] )  )
inter_cell_T =  Qobj( np.array( [[ 0, 0 ],[t_inter,0]] )  ) 
apSSH_lattice_nTrI = Lattice1d(num_cell=100, boundary = "aperiodic", cell_num_site = 2, cell_site_dof = [1], cell_Hamiltonian = cell_H, inter_hop = inter_cell_T )

apSSH_lattice_nTrI.winding_number()

apSSH_H_nt = apSSH_lattice_nTrI.Hamiltonian()
[n_D,Vx] = apSSH_H_nt.eigenstates()

plt.plot(n_D,'ro')
plt.ylabel('eigen values')
plt.show()
plt.close()

"""
Explanation: Topologically Nontrivial Insulator: Hardwall Boundary Condition
To reveal the topological nontriviality, we have to put the insulator next to a 
boundary, i.e. use a hardwall boundary condition and the spectrum of eigenvalues of
the Hamiltonian shows the conspicuous edge states which plotted as a eigenvector
are dense towards the edge considerably. However, the dispersion looks completely the
same.
End of explanation
"""


xA = [i for i in range(200)]
Es0 = np.abs(Vx[99])
Es1 = np.abs(Vx[100])
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 4), sharey=True)
ax1.plot(xA, Es0, label="State0")
ax2.plot(xA, Es1, label="State1")
ax1.legend()
ax2.legend()
plt.show()
fig.suptitle('Mid-gap Edge states')
plt.close()

"""
Explanation: We observe the much celebrated mid-gap edge states, who appear at mid-gaps. The nontriviality of the topological insulator is revealed with their placement next to a boundary. The edge states are the consequences of the bulk-bundary correspondence. Plotting the absolute values of the eigen-functions reveal that they are localized at the edges. We plot the absolute values of the eigenfunctions corresponding to the two eigenvalues found at mid-gaps.
End of explanation
"""


n_D[96:104]

"""
Explanation: It is quite obvious, the edge states are localized towards the edges.
The eigenstates are not unique, they are only unique within a phase factor. The routine used in QuTiP only yields eigenstates that are one of the possibilities within the phase factor.
Let us check out the two edge state eigenvalues as well as the three preceding and three succeeding ones.
End of explanation
"""


chiral_op_nTrI = apSSH_lattice_nTrI.distribute_operator(sigmaz())
dimH_nTrI = chiral_op_nTrI.dims
identity_H_nTrI = Qobj(np.identity(200), dims=dimH_nTrI)
identity_H_nTrI
PA_200 = 0.5*( identity_H_nTrI + chiral_op_nTrI)
PB_200 = 0.5*( identity_H_nTrI - chiral_op_nTrI)

xA = [i for i in range(200)]
Es0 = np.abs(PA_200*Vx[99])
Es1 = np.abs(PB_200*Vx[99])
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 4), sharey=True)
ax1.plot(xA, Es0, label="State0")
ax2.plot(xA, Es1, label="State1")
ax1.legend()
ax2.legend()
plt.show()
fig.suptitle('Mid-gap Edge states')
plt.close()

"""
Explanation: It is obvious that, n_D[99] and n_D[100] are values 0 within numerical precision. There is an interesting property of the eigenstates of a chiral Hamiltonian with eigenvalue 0.
\begin{eqnarray}
H|\psi_n \rangle = 0 \implies H(\hat{P}_{A/B}| \psi_n\rangle) = \frac{1}{2}H(\hat{\bf{1}} \pm \Sigma_z )| \psi_n \rangle = \frac{1}{2}(\hat{\bf{1}} \mp \Sigma_z )H| \psi_n \rangle = 0
\end{eqnarray}
Therefore, $\hat{P}_{A/B} |\psi_n\rangle$ is an eigenstate with eigenvalue 0 as well. As a result, the edge-states can be chosen to have support in only one sublattice.
End of explanation
"""


is_null = (np.abs(apSSH_H_nt*PA_200*Vx[99]) < 1E-10).all()
print(is_null)

qutip.about()

qutip.cite()

"""
Explanation: $\hat{P}{A} |\psi_n\rangle$ and $\hat{P}{B} |\psi_n\rangle$ are orthogonal eigenstates with eigenvalue 0, and we can choose the two edge states to be localized at two edges of the lattice. Just to be sure, we also check if $H\hat{P}_{A/B} |\psi_n\rangle $ is indeed 0.
End of explanation
"""


