import numpy as np, random, scipy.linalg

"""
Explanation: Table of contents

Mathematical Background
The Null space and its orthogonal
Probabilistic interpretation of the curvature


Code documentation
A convenient basis
The Markov chain
Computation of the curvature


Appendix
Lemmas: the structure of zero modes



Mathematical Background
We investigate constant plaquettes of size d, with one plaquette per site. In order terms, we consider the map $$\begin{matrix}\mathbb{C}^{2^d}&\mapsto &(\mathbb{C}^2)^{\otimes N}\(C^{\delta_1\ldots \delta_d}) &\mapsto & \sum_{\epsilon}\prod_{i=1}^{N}C^{\epsilon_i\cdots \epsilon_{i+d-1}}|\epsilon\rangle,\end{matrix}$$ which to $C$ associates $|\phi(C)>$. This map is holomorphic.
We focus on the local geometry of this map. As we will see, it is never an immersion. The null tangential vectors are explicit, and we will also give a basis of their orthogonal space.
We also give a numerical scheme to compute, up to machine precision, the curvature, which is a symmetric 2-form on the holomorphic tangent space: $$G(\eta,\xi)=\overline{\partial}{\eta}\partial{\xi}\log(\langle\phi(C)|\phi(C)\rangle).$$ This form is positive, as the Levi form of a plurisubharmonic function, but it is definite positive only when restricted to the orthogonal of the null space above.
The Null space
In order to understand completely the geometry of the image set above, we will make two simplifying assumptions:
* N is a multiple of d(d-1);
* All coefficients $C^{\delta_1,\ldots,\delta_d}$ are non-zero.
The first hypothesis might seem strange, since we are mainly interested in the large $N$ case. In fact, when $N$ is a multiple of $d(d-1)$ the proofs are much simpler, and when it is not the case (or when we consider plaquettes on a chain, without periodicity), one must add terms to $G$, of order $N^{-1}$.
The second hypothesis is of course true in a generic configuration; it is also true along a (real) one-dimensional family of configurations. As we are interested in time evolution, one can argue that this hypothesis will remain valid along trajectories.
As all coefficients are non-zero, we make use of the projective properties of the map $\phi$ to consider a basis for variations of parameters which is more convenient than the natural one. If $\alpha \in (\mathbb{C}^2)^{\otimes d}$, we let $$\alpha \cdot |\phi(C)\rangle=\sum_{\delta}\alpha_{\delta}C^{\delta}\partial_{C^{\delta}}|\phi(C)\rangle.$$
This is motivated by the fact that $$C^{\delta}\partial_{C^{\delta}}|\phi(C)\rangle=\sum_{i=1}^N\sum_{\large{\epsilon\in A_i^{\delta}}}\prod_{j=1}^NC^{\epsilon_j\ldots \epsilon_{j+d-1}}|\epsilon\rangle.$$ Here, $$A_i^{\delta}={\epsilon\,|\,\epsilon_i=\delta_1,\ldots,\epsilon_{i+d-1}=\delta_d}.$$
The main result is
Theorem 1 $\alpha\cdot |\phi(C)\rangle=0$ if and only if $\alpha=(|0\rangle+|1\rangle)\otimes v-v\otimes(|0\rangle+|1\rangle)$ for some $v\in (\mathbb{C}^2)^{\otimes d-1}$.
Proof of Theorem 1 One inclusion proceeds by direct computation. Indeed, if $v=\sum v_{\delta}|\delta\rangle$, then
$$(|0\rangle+|1\rangle)\otimes v \cdot |\phi(C)\rangle = \sum_{\delta}\sum_{i=1}^Nv_{\delta}\sum_{\large{\epsilon\in A_i^{0\delta}\cup A_i^{1\delta}}}\prod_{j=1}^NC^{\epsilon_j\cdots \epsilon_{j+d-1}}|\epsilon\rangle,$$ and $$A_i^{0\delta}\cup A_i^{1\delta}=A_{i+1}^{\delta0}\cup A_{i+1}^{\delta_1},$$ so that $$(|0\rangle+|1\rangle)\otimes v-v\otimes(|0\rangle+|1\rangle) \cdot |\phi(C)\rangle=0.$$
The other inclusion is less direct and makes full uses of the hypotheses above.
Since all coefficients are non-zero, then $\alpha\cdot |\phi(C)\rangle=0$ if and only if, for every $\epsilon$, one has $\sum_{i=1}^N\alpha_{\epsilon_i\cdots \epsilon_{i+d-1}}=0.$
In particular, let $\delta\in {0,1}^d$ arbitrary, and let $\epsilon=\delta \delta \delta...$ a concatenation of size $N$. (We assumed that $N$ is a multiple of $d$). Then the equation above yields $\sum_{k=0}^{d-1}\alpha_{ls(d)}=0$, with $ls$ the operator of left shift.
In particular, $$\sum_{k=0}^{d-1}ls^k(\alpha)=0.$$
Using Lemma 1 in the Appendix, it follows that $\alpha=v-ls(v)$ for some $v$.
We will now prove that, for every word $w\in {a,b}^{d-2}$, one has $$\alpha_{0w0}+\alpha_{1w1}-\alpha_{0w1}-\alpha{1w0}=0.$$ Using Lemma 3, with $a=\frac{|0\rangle+|1\rangle}{\sqrt{2}}$ and $b=\frac{|0\rangle-|1\rangle}{\sqrt{2}}$, this will conclude the proof.
Recall that each $\epsilon \in {0,1}^N$ yields an equation on $\alpha$. Substracting the equation given by $\epsilon=...0w0w1w0w0...$ (recall N is a multiple of $d-1$) from the one given by $\epsilon=...0w0w0...$ yields $$\alpha_{0w1}+\ldots+\alpha_{1w0}=\alpha_{0w0}.$$
Now, substracting the equation given by $\epsilon=...0w0w1w1w0w0...$ from $\epsilon=...0w0w0...$ yields $$\alpha_{0w1}+\ldots+\alpha_{1w0}-\alpha_{1w0}+\alpha{1w1}-\alpha{0w1}+\alpha{0w1}+\ldots+\alpha{1w0}=\alpha{0w0}.$$
Hence,
$$\alpha_{1w1}+\alpha_{0w0}=\alpha_{0w1}-\alpha_{1w0}.$$
This concludes the proof.
Probabilistic interpretation of the curvature
The expression of the derivative of $|\phi(C)\rangle$ along a base vector was formulated as a sum with the same terms as the initial sum, but restricted on some configuration set $A^{\delta}_i$. These sums are ubiquitous in the computation of the curvature, which will allow a probabilistic reformulation. Since we are dealing with one-dimensional loops, the probabilistic model will be a Markov chain (restricted to identical final and initial states). If all coefficients are non-zero, this Markov chain is mixing, and this will allow a precise numerical treatment.
We want to compute $$G(\alpha,\beta)=\alpha\overline{\cdot}\beta\cdot \log(\langle \phi(C)|\phi(C)\rangle).$$ Here, $\overline{\cdot}$ stands for anti-holomorphic derivation. This is exactly $$G(\alpha,\beta)=\cfrac{\langle \alpha\cdot \phi(C)|\beta\cdot \phi(C)\rangle}{\langle \phi(C)|\phi(C)\rangle} - \cfrac{\langle \alpha\cdot \phi(C)|\phi(C)\rangle \langle \phi(C)|\beta\cdot \phi(C)\rangle}{(\langle \phi(C)|\phi(C)\rangle)^2}.$$
Now $$\langle \phi(C)|\phi(C)\rangle=\sum_{\epsilon}\prod_{j=1}^N|C^{\epsilon_j\ldots\epsilon_{j+d-1}}|^2.$$ Moreover, if $\alpha=C^{\delta}\partial_{C^{\delta}}$ and $\beta=C^{\delta'}\partial_{C^{\delta'}}$, then
$$\langle \alpha\cdot \phi(C)|\phi(C)\rangle=\sum_{i=1}^N\sum_{\large{\epsilon\in A^{\delta}i}}\prod_j|C^{\epsilon_j\ldots \epsilon{j+d-1}}|^2.$$ Finally,
$$\langle \alpha\cdot \phi(C)|\beta\cdot \phi(C)\rangle=\sum_{i_1=1}^N\sum_{i_2=1}^N\sum_{\large{\epsilon\in A^{\delta}{i_1}\cap A^{\delta'}{i_2}}}\prod_j|C^{\epsilon_j\ldots \epsilon_{j+d-1}}|^2.$$
Those expressions are quite clumsy. Let $({0,1}^N,\mathcal{P}({0,1}^N),\mathbb{P}C)$ be a probabilized space whose elementary events are words of length $N$ on ${0,1}$, with $$\mathbb{P}_C(\epsilon)=(\prod{j=1}^N|C^{\epsilon_j\cdots \epsilon_{j+d-1}}|^2)\langle \phi(C)|\phi(C)\rangle^{-1}.$$
Then the matrix elements of $G$ are simply $$G_{\delta,\delta'}=\sum_i\sum_j\mathbb{P}_C(A^{\delta}_i\cap A^{\delta'}_j)-\mathbb{P}_C(A^{\delta}_i)\mathbb{P}_C(A^{\delta}_j).$$
To evaluate those probabilities, we will reformulate this probabilized space with a Markov chain.
Consider the Markov chain, with ${0,1}^d$ as set of states, with transition matrix:
$$T(\delta_1\ldots \delta_d,\delta_2\ldots\delta_{d+1})=\cfrac{|C^{\delta_2\ldots\delta_{d+1}}|^2}{|C^{\delta_2\ldots \delta_d 0}|^2+|C^{\delta_2\ldots \delta_d 1}|^2}.$$
Consider realisations of this Markov chain, of length $N+1$, conditionned to having the same final and initial state. Since the $d-1$ last digits at one time coincide with the $d-1$ first digits at the following time, this realisation is a sequence of $0$ and $1$, of length $N+d$ (where the state, at step $j$, consists of the digits $j$ to $j+d-1$.). Moreover, the $N+1$ digits corresponds to the first one, and so on. Then the probability of this sequence $\epsilon$ is exactly $\mathbb{P}_C(\epsilon)$, as long as the initial state is chosen according to the law of $A_1$.
The matrix $T$ verifies the hypotheses of the Perron-Frobenius theorem: since all coefficients are non-zero, the matrix $T^d$ has only positive entries. In particular $T$ has 1 as only eigenvalue of modulus one (all other eigenvalues have strictly smaller modulus), and this eigenvalue is simple. From there we are able to identify the law of $A_i$. Indeed, on one hand the original probability space is rotation invariant so that $A_i$ and $A_j$ have the same law for any $i$ and $j$; on the other hand the law of $A_{i+1}$ (which is a $2^d$-dimensional vector) is $T$ times the law of $A_i$, plus an error of size $e^{-cN}$. Hence, if $\mu_C$ is the equilibrium measure for $T$, then $$\mathbb{P}_C(A^{\delta}_i)=\mu_C(\delta)+O(e^{-cN}).$$
Along the same lines, if $j$ is closer from $i$ to the left than to the right (meaning $j$ is between $i$ and $i+N/2$ on $\mathbb{Z}/N\mathbb{Z}$, then $$\mathbb{P}_C(A^{\delta'}_j|A^{\delta}_i)=T^{j-i}[\delta',\delta]+O(e^{-cN}).$$ In particular, $$\mathbb{P}_C(A^{\delta'}_j|A^{\delta}_i)=\mathbb{P}_C(A^{\delta'}_j+O(e^{-c|i-j|}).$$ The intuition of this estimate is that, if $i$ and $j$ are far apart (more than the reduced spectral radius of $T$ allows for), then $A_i$ and $A_j$ are independent up to an exponentially small error.
This offers a numerical scheme to compute $G$ up to machine precision. First, we compute and store the sequence of positive powers of $T$, until they converge to a matrix of rank one (at step $k_{M}$). Each line of this matrix is the invariant measure $\mu_C$. Then $$N^{-1}G_{\delta,\delta'}=\sum_{j=0}^{k_M}\mu_C(\delta)(T^j_{\delta',\delta}-\mu_C(\delta'))+\sum_{j=1}^{k_M}\mu_C(\delta')(T^j_{\delta,\delta'}-\mu_C(\delta'))+ O(\text{machine precision})+O(Ne^{-N}).$$
Code documentation
First and foremost, we need to import standard numerical libraries, as well as scipy.linalg to show the spectrum of the curvature form
End of explanation
"""


# %load ../curv_basis.py
# a tentative basis

def lshift(c,d):
    return (2*c)%(2**d)+c/(2**(d-1))

def rtldigits(c,d):
    """Returns the reversed digits of c in base 2, of length d"""
    ccopy=c
    dig=[]
    for j in range(d):
        dig.append(ccopy%2)
        ccopy=ccopy/2
    return np.array(dig)

def rotation(d=3):
    """Rotates from the base along z, to the base along x."""
    U=np.zeros((2**d,2**d))
    for c in range(2**d):
        for k in range(2**d):
            U[c,k]=(-1)**(sum((1-rtldigits(c,d))*(1-rtldigits(k,d))))/2**(0.5*d)
    return U

def fullbasis(d=3):
    U=np.zeros((2**d,2**d))
    line=0
    #First we look at the zero vectors
    #they look like 1c-c1 for some c
    #but we have to give an orthogonal basis !
    for j in range(d-2):
        for c in range(2**j):
            #we consider the vectors from 1...10c0 to 0c01.....1 of size d-1
            #here they form a Toeplitz matrix of rank d-2-j,
            #which can easily be diagonalized
            for k in range(d-2-j):
                for l in range(d-2-j):
                    vect=2**(d-1)-2**(l+j+2)+c*2**(l+1)+2**l-1
                    val=np.sin(np.pi*(k+1)*(l+1)/(d-1-j))
                    U[line+k,2**(d-1)+vect]+=val
                    U[line+k,2*vect+1]-=val
                U[line+k]/=np.sqrt(np.dot(U[line+k],U[line+k]))
            line+=d-2-j
    #then we only forgot the vectors 1.....10 to 01......1
    for k in range(d-1):
        for l in range(d-1):
            vect=2**(d-1)-1-2**l
            val=np.sin(np.pi*(k+1)*(l+1)/d)
            U[line+k,2**(d-1)+vect]+=val
            U[line+k,2*vect+1]-=val
        U[line+k]/=np.sqrt(np.dot(U[line+k],U[line+k]))
    #The vector 1...1 is itself orthogonal
    line+=d-1
    U[line,2**d-1]=1.
    line+=1
    #There are several types of orthogonal (nonzero) vectors:
    #All 0c0 vectors
    for c in range(2**(d-2)):
        U[line+c,2*c]=1
    line+=2**(d-2)
    #All 10c0+0c01 vectors
    for c in range(2**(d-3)):
        U[line+c,2*c+2**(d-1)]=1./np.sqrt(2)
        U[line+c,4*c+1]=1./np.sqrt(2)
    line += 2**(d-3)
    #Some weird vectors
    for c in range(2**(d-3)):
        current=2*c+2**(d-1)+2**(d-2)
        shiftlist=[current]
        while current>=2**(d-1):
            current=lshift(current,d)
            shiftlist.append(current)
        for shifted in shiftlist:
            U[line+c,shifted]=1./np.sqrt(len(shiftlist))
    return U


"""
Explanation: A convenient basis
The file curv_basis.py contains the decomposition along the null space and a complementary subspace, according to the discussion above.
End of explanation
"""


#We check that the change of basis U is indeed orthogonal
U=np.dot(fullbasis(4),rotation(4))
np.sum(np.abs(np.dot(U,U.T)-np.eye(2**4)))

"""
Explanation: The first step is the function rotation, which returns the matrix which exchanges the basis given by $(|0\rangle,|1\rangle)$ with the basis given by $(2^{-\frac 12}(|0\rangle+|1\rangle),2^{-\frac 12}(|0\rangle+|1\rangle))$, on $(\mathbb{C}^2)^{\otimes N}$.
The second step is the function fullbasis, which computes the decomposition of the space given by Lemma 2 (with $a=|0\rangle$ and $b=|1\rangle$. Here, the non-trivial part is to give an orthonormal basis of the space $E$.
Indeed, if $v$ is a base vector, $\langle a\otimes v-v\otimes a|a\otimes v'-v'\otimes a\rangle$ is equal to $2$ if $v=v'$ (unless $v=aa.....a$), and otherwise is equal to $-1$ if $v\otimes a=a\otimes v'$ (which happens only if $v$ begins by an $a$ and $v'=ls(v)$) or if $v'\otimes a=a\otimes v$ (which happens only if $v$ ends by an $a$ and $v=ls(v')$). The exceptional case when both $v\otimes a=a\otimes v'$ and $v'\otimes a=a\otimes v$ means that $v=awa$ and $v'=aw'a$ with $v'=ls(v)$ so that $v=a...a$, a case which was excluded a priori.
Hence, the structure of the matrix given by the previous scalar product is formed of tridiagonal Toeplitz matrices, with $2$ on the diagonal and $-1$ on the over- and underdiagonal; these matrices can be explicitly diagonalized.
The natural basis of $F$ and $G$ given in Lemma 2 is orthogonal (as every word appears only once, either in $F$, or in $G$), and requires only a normalization.
End of explanation
"""


# %load ../curvature.py
import numpy as np

def generate_C(d,m=1,M=2):
    """
    Generates a bunch of randomly chosen coefficients, far from zero
    """
    C=[]
    for k in range(2**d):
        C.append(random.random()*(M-m)+m)
    return C

def sparse_mult(mat,C):
    """returns an efficient multiplication of mat by the transfer matrix.
    Runs in quadratic time.
    """
    prod=np.zeros((len(C),len(C)),dtype=float)
    for i in range(len(C)):
        for j in range(len(C)):
            prod[i,j]=mat[i,(2*j)%len(C)]*C[(2*j)%len(C)]
            prod[i,j]+=mat[i,(2*j+1)%len(C)]*C[(2*j+1)%len(C)]
            prod[i,j]/=C[(2*j)%len(C)]+C[(2*j+1)%len(C)]
    return prod

def Markov_powers(C,ERROR=1.0e-14):
    """Computes the powers of the Markov chain
    for constant plaquettes"""
    powerseq=[np.eye(len(C))]
    converged=False
    while not converged:
        current=powerseq[-1]
        new=sparse_mult(current,C)
        powerseq.append(new)
        dist=0.
        for k in range(len(C)/2):
            dist+=abs(new[k,0]-new[k,1])
        if dist<ERROR:
            converged=True
        #print dist
    print "Markov chain convergence after",len(powerseq),"steps"
    return powerseq

def curvature(C,ERROR=1.0e-14):
    """Computes, via a Markov chain, the curvature tensor
    for constant plaquettes up to some prescribed error."""
    #compute the sequence (M**i) until convergence
    powerseq=Markov_powers(C,ERROR)
    #the lines of the last matrix are the equilibium measure
    eq_m=powerseq[-1].T[1]
    print "Equilibrium measure: ", eq_m
    G=np.zeros((len(C),len(C)),dtype=float)
    M=powerseq[1]
    for i in range(len(C)):
        for j in range(len(C)):
            for k in range(len(powerseq)):
                current=powerseq[k]
                G[i,j]+=eq_m[j]*(current[i,j]-eq_m[i])
                if k !=0:
                    G[i,j]+=eq_m[i]*(current[j,i]-eq_m[j])
    return G


"""
Explanation: Computation of the curvature
The file curvature.py encodes the Markov chain to compute the curvature.
End of explanation
"""


C=generate_C(4)

"""
Explanation: We first have to generate the square norm of the coefficients $C$; here they are chosen as uniform independent random variables, shifted away from zero (as we want to ensure that all coefficients are non-zero)
End of explanation
"""


G=curvature(C)

"""
Explanation: Then we can compute the curvature (which, for the moment, is expressed in the canonical basis).
End of explanation
"""


G

"""
Explanation: As a control tool, we printed the step at which the Markov chain converges, and the equilibrium measure.
The matrix G is indeed quite complicated at this step.
End of explanation
"""


np.dot(U,np.dot(G,U.T))

"""
Explanation: However, in the base given by the matrix $U$, the quadratic form is simpler.
End of explanation
"""


Gred=np.dot(U,np.dot(G,U.T))[8:15].T[8:15]
scipy.linalg.eigvalsh(Gred)

"""
Explanation: As proved in Theorem 1, only the lower right quadrant of this matrix is non-zero; moreover this quadrant is positive.
End of explanation
"""


import matplotlib.pyplot as plt
C=generate_C(9)
G=curvature(C)
U=np.dot(fullbasis(9),rotation(9))
Gred=np.dot(U,np.dot(G,U.T))[2**8:2**9-1].T[2**8:2**9-1]
vals=scipy.linalg.eigvalsh(Gred)
plt.plot(np.log(vals),'ro')
plt.show()

"""
Explanation: This code has exponential complexity in $d$ (both in time and in memory). There is numerical evidence that the spectral gap of the Markov chain is uniform in $d$ (as long as the coefficients $C$ are uniformly away from zero and infinity), and the code works on a laptop up to $d=9$ or $10$.
End of explanation
"""


