# Here we define functions that allow us to compute eigenvalues of subshifts and to plot Rauzy fractals
# The code to plot Rauzy fractal is inspired by the one of Timo Jolivet,
# but it correct some bugs and generalize to any infinite word and any projection.
#
# Exemples are below

from badic import *

def ReturnSubstitution (s, u, getrw=False, verb=False):
    """
    Give the return substitution of s on the return word u.
    u is assumed to be a prefix of an infinite fixed point.

    INPUT:

    - ``s`` - a WordMorphism

    - ``u`` - a word

    - ``getrw`` - Boolean (default: ``False``) - if True, return also the set of return words.

    - ``verb`` - Boolean (default: ``False``) - if True, print more informations.

    OUTPUT:

    A WordMorphism, and also a dictionnary of return words if getrw is True.

    EXAMPLES::
    
        sage: s = WordMorphism('1->213,2->4,3->5,4->1,5->21')
        sage: ReturnSubstitution(s^3, '1')
        WordMorphism: 0->01012, 1->0102012, 2->0102

    """
    A = list(s.domain().alphabet())
    x = s.fixed_point(u[0])
    m = s.incidence_matrix()
    R = dict() # dict of return words
    nR = 0 # number of return words
    if verb:
        print(A)
    a = dag.AnyWord(A,A).concat(dag.Word(u,A=A))
    A = a.alphabet
    v = vector([0 for i in range(len(A))])
    rv = copy(v)
    dA = dict()
    if verb:
        print(A)
    for i,t in enumerate(A):
        dA[t] = i
    vu = copy(v) # abelianized of u
    for t in u:
        vu[dA[t]] += 1
    lu = len(u)
    lv1 = [] # lv1[i] = abelian vector of starting of s(m[i]) where m[i] is the ith return word
    lv2 = [] # lv2[i] = abelian vector of ending of s(m[i]) where m[i] is the ith return word
    lim = dict() # lim[i] = image if m[i] by s as return words
    cim = 0 # current return word for which we compute the image
    lcim = [] # image of the current return word
    lstop = []
    if verb:
        print(a)
    i = 0
    i1 = -1 # d?but du mot de retour courant
    crw = -1 # current return word
    e = a.initial_state
    while True:
        if verb:
            print("i=",i)
        l = dA[x[i]]
        e = a.succ(e, l)
        if verb:
            print('l={}, v={}'.format(l,v))
        v[l] += 1 # update the abelian vector
        #print(e)
        if a.is_final(e):
            if i1 < 0:
                if i != lu-1:
                    raise ValueError("u={} should be a prefix of a fixed point !".format(u))
                i1 = 0
            else:
                # return word found
                i2 = i+1-lu
                r = x[i1:i2]
                i1 = i2
                #print("i1={}".format(i1))
                if r not in R:
                    R[r] = nR
                    #print(r)
                    if verb:
                        print("r={}".format(r))
                        print("rv={}, v-vu={}".format(rv, v-vu))
                    lv1.append(m*rv)
                    lv2.append(m*(v-vu))
                    if verb:
                        print(lv1, lv2)
                    lstop.append(sum(lv2[nR]))
                    nR += 1
                rv = copy(v-vu)
                lcim.append(R[r])
                if i1 >= lstop[cim]:
                    if i1 > lstop[cim]:
                        raise RuntimeError("D?passement !!!")
                    # we read all the image of the current return word
                    if verb:
                        print("image de {} : {}".format(cim, lcim))
                    lim[cim] = lcim
                    lcim = []
                    cim += 1
                    if verb:
                        print("\ncim = {}".format(cim))
                    if cim == nR:
                        # end
                        if getrw:
                            return WordMorphism(lim), R
                        return WordMorphism(lim)
                    # change the starting point in order to read the image of the current return word
                    if verb:
                        print("v={}".format(v))
                    rv = lv1[cim]
                    v = copy(lv1[cim])+vu
                    if verb:
                        print("v={}".format(v))
                        print("stop={}".format(lstop[cim]))
                    i1 = sum(lv1[cim])
                    i = i1 + lu-1
                    if verb:
                        print("i1={}".format(i1))
                        print(lv1)
                    e = a.final_states[0]
                if verb:
                    print(r, R[r])
        i += 1
        #if i > 1000:
        #    break
    if verb:
        print("nR = {}, cim={}".format(nR, cim))
        print(lim)
        print("Return words : {}".format(R))

def Prop8 (su, R, verb=False):
    """
    Return the substitution of Proposition 8 in [Durand-Petite].
    
    INPUT:
    
    - ``su`` - a return substitution
    
    - ``R`` - the set of return words
    
    - ``verb`` - Boolean (default: ``False``) if True, print additionnal informations.
    
    """
    i = 1
    while(True):
        su2 = su^i
        for k,j in enumerate(R):
            if len(j) > len(su2(k)):
                break
        else:
            break
        i += 1
    if verb:
        print("su ?l?v? ? la puissance {}".format(i))
    su = su^i

    ln = [len(j) for j in R]
    lm = [len(su(j)) for j in range(len(R))]

    def psi (lj):
        r = []
        for j in lj:
            r += [(j,i) for i in range(ln[j])]
        return r
    if verb:
        print(R)
    B = [(R.index(j),p) for j in R for p in range(len(j))]
    xi = dict()
    for (j,p) in B:
        if p < ln[j]-1:
            xi[(j,p)] = psi([su(j)[p]])
        else:
            xi[(j,p)] = psi(su(j)[ln[j]-1: lm[j]])
    if verb:
        print(xi)
    xi2 = dict()
    for (j,p) in xi:
        xi2[B.index((j,p))] = [B.index(t) for t in xi[(j,p)]]
    return WordMorphism(xi2)

def initials_period (s):
    """
    Return the initials period of s.
    Starting letters of s^n are eventually periodic,
    with a period which is the initials period.
    """
    A = s.domain().alphabet()
    a = DetAutomaton([(a, 0, s(a)[0]) for a in A], avoidDiGraph=True)
    return lcm([len(c) for c in a.strongly_connected_components()])

def proprify (s, u=None, verb=False, check_non_trivial_coboundaries=False):
    """
    Return a substitution whose a power is left-proper.
    If check_non_trivial_coboundaries is set to True, then it checks if there is no non-trivial coboundary
    and if it is the case, then it return the substitution s^k where k is the initials period.
    """
    # iterate to have a fixed point
    for k in range(1,2000):
        fp = (s^k).fixed_points()
        if len(fp) > 0:
            break
    if verb:
        print("power {}".format(k))
    s = s^k
    if u is None:
        u = s.fixed_points()[0][:1] # choose the first letter of a fixed point
    if verb:
        print("u={}".format(u))
    # compute the return substitution on u
    su,R = ReturnSubstitution(s, u, True)
    R = list(R.keys())
    if verb:
        print("return subst: {}".format(su))
        print("return words: {}".format(R))
    if check_non_trivial_coboundaries:
        # check if there are non-trivial coboundaries
        mr = matrix([w.abelian_vector() for w in R])
        mr2 = []
        for l in mr.BKZ():
            if not l.is_zero():
                mr2.append(l)
        mr2 = matrix(mr2)
        if mr2.is_square():
            d = abs(det(mr2))
            if d == 1:
                if verb:
                    print("No non-trivial coboundary.")
                ip = initials_period(s)
                if verb:
                    print("Initials period = {}".format(ip))
                return s^ip
            else:
                if verb:
                    print("Potential coboundary with den {}".format(d))
    return Prop8(su, R)

def pseudodet (m):
    """
    Return the pseudo-determinant of the matrix.
    """
    return product([t for t in m.eigenvalues() if t != 0])

def generalized_eigenvectors (m, verb=False):
    """
    Return a concatenation of basis of generalized eigenspaces for eigenvalues >= 1.
    """
    pi = m.charpoly()
    I = identity_matrix(m.nrows())
    r = []
    rvp = []
    for pi, k in m.charpoly().factor():
        if verb:
            print(pi, k)
        lr = pi.roots(ring=QQbar)
        for vp, _ in lr:
            if verb:
                print(vp)
            if abs(vp) >= 1:
                K.<b> = NumberField(vp.minpoly(), embedding=vp)
                m2 = matrix(m, ring=K)
                if verb:
                    print("vp = {}, multiplicity = {}".format(vp, k))
                r += matrix((m-vp*I)^k, ring=K).right_kernel().basis()
                #r += ((m2-vp*I)^k).right_kernel().basis()
                rvp.append(vp)
    return r, rvp

def diff(v1, v2):
    """
    Return the difference of two vectors in two different NumberField.
    """
    K = v1[0].parent()
    R.<x> = PolynomialRing(K)
    pi = R(v2[0].parent().polynomial())
    c = v2[0].parent().gen()
    for f,_ in pi.factor():
        if f(c) == 0:
            L = K.extension(f, "c")
            break
    K = L.absolute_field('a')
    pi = K.polynomial()
    lr = K.polynomial().roots(ring=QQbar)
    for r,_ in lr:
        try:
            K.<b> = NumberField(pi, embedding=r)
            return vector([K(t) for t in v1]) - vector([K(t) for t in v2])
        except:
            continue

def get_w_matrix (m, verb=False):
    """
    Return an integer matrix whose left-kernel is the set of w satisfying the condition of Theorem 1.1 in [Me].
    """
    global K
    if abs(pseudodet(m))-1 != 0:
        raise ValueError("The matrix is not pseudo-unimodular.")
    b = max(m.eigenvalues())
    if verb:
        print("Perron eigenvalue: {}".format(b))
    lv, lvp = generalized_eigenvectors(m, verb=verb)
    if verb:
        print("Basis of generalized eigenspaces: {}".format(list(zip(lv,lvp))))
    # compute a basis of the orthogonal of w
    lvo = []
    nn = 0
    for v in lv:
        if sum(v) == 0:
            lvo.append(v)
        else:
            v /= sum(v)
            if nn == 0:
                v0 = v
            else:
                lvo.append(diff(v,v0))
            nn += 1
    if verb:
        print("Basis of the orthogonal of w: {}".format(lvo))
    # convert each vector in the splitting field
    lm = []
    for v in lvo:
        lm.append(matrix([list(t) for t in v]))
    wm = block_matrix([lm])
    den = 1
    for v in wm:
        for t in v:
            den = lcm(den, t.denom())
    wm = matrix(wm*den, ring=ZZ)
    if verb:
        print("Matrix whose integer left kernel is the set of possible w: {}".format(wm))
    return wm

def Perron_eig (m):
    """
    Return the Perron eigenvector of sum 1 with coefficients in a NumberField.
    """
    b = max(m.eigenvalues())
    K.<b> = NumberField(b.minpoly(), embedding=b)
    m = matrix(m, ring=K)
    v0 = (m - b*identity_matrix(K, m.nrows())).right_kernel().basis()[0] # compute the Perron eigenvector
    v0 /= sum(v0) # of sum 1
    return v0

def eigenvalues_proper (m, verb=False):
    """
    Compute the eigenvalues of the subshift from the incidence matrix of the substitution.
    We assume that eigenvalues are exp(2 i pi a), with a(1,...,1)m^n ---> 0 mod 1,
    which is the case if the substitution is proper.
    """
    v0 = Perron_eig(m)
    wm = get_w_matrix(m, verb=verb)
    # compute eigenvalues
    lvp = []
    if wm.nrows() == 0:
        lk = identity_matrix(ZZ, m.nrows())
    else:
        lk = wm.kernel().basis() # compute solutions in ZZ
    if verb:
        print("Basis of the Z-module where w lives: {}".format(lk))
    for k in lk:
        vp = k*v0
        lvp.append(vp)
    if verb:
        print("List of (non-independant) eigenvalues: {}".format(lvp))
    # convert in the number field of b
    b = max(m.eigenvalues())
    K.<b> = NumberField(b.minpoly(), embedding=b)
    lvp = [K(t) for t in lvp]
    # convert to a matrix
    mvp = matrix([list(t) for t in lvp])
    if verb:
        print("Corresponding matrix: {}".format(mvp))
    den = 1
    # compute denominator
    for v in mvp:
        for t in v:
            den = lcm(den, t.denom())
    if verb:
        print("den = {}".format(den))
    # convert to a matrix in ZZ
    mvp2 = matrix(mvp*den, ring=ZZ)
    # compute a basis of the lattice
    mvp3 = mvp2.BKZ()
    if verb:
        print("matrix after BKZ: {}".format(mvp3))
    # return the basis in a number field
    lvp = [K(t)/den for t in mvp3 if K(t) != 0]
    b = max(m.eigenvalues())
    if b.degree() == len(lvp):
        if verb:
            print("Express the result in the NumberField of the Perron eigenvalue.")
        # express the result in the NumberField of b
        return lvp
    else:
        if verb:
            print("Find a small NumberField to express the result.")
        # find a small NumberField to express the result
        pi = product([t.minpoly() for t in lvp])
        K.<b> = pi.splitting_field()
        if verb:
            print("Splitting field: {}".format(K))
        pi = K.polynomial()
        lr = pi.roots(ring=QQbar)
        for r, _ in lr: # try to find a correct embedding
            try:
                L.<b> = NumberField(pi, embedding=r)
                return [L(t) for t in lvp]
                break
            except:
                continue
        raise RuntimeError("Could not find a field to express the result !")

def eigenvalues(s, verb=False, check_non_trivial_coboundaries=False):
    """
    Compute eigenvalues of the subshift of the substitution s.
    """
    sp = proprify(s, None, False, check_non_trivial_coboundaries)
    if verb:
        print("proprified substitution: {}".format(sp))
    m = sp.incidence_matrix()
    if verb:
        print(m.charpoly().factor())
    lvp = eigenvalues_proper(m, verb=verb)
    m = s.incidence_matrix()
    b = max(m.eigenvalues())
    if len(lvp) == b.degree():
        if verb:
            print("convert in another NumberField")
        # express the result in the NumberField of b
        K.<b> = NumberField(m.charpoly(), embedding=b)
        return [K(t) for t in lvp]
    return lvp

def pseudoinv (m):
    """
    compute the pseudo-inverse of m
    im(m) and ker(m) are assumed to be supplementary
    """
    # basis of kernel
    mk = matrix(m.right_kernel().basis()).transpose()
    # basis of image
    mi = matrix(m.transpose().image().basis()).transpose()
    # diagonalize by block
    p = matrix(mk.columns()+mi.columns()).transpose()
    md = p^(-1)*m*p
    # restriction of m to the image
    N = md.submatrix(mk.ncols(),mk.ncols(),mi.ncols(),mi.ncols())
    # projection to the image along the kernel
    Vi = (p^(-1)).submatrix(mk.ncols(),0,mi.ncols(),m.nrows())
    return mi*N^(-1)*Vi

def usual_projection (s):
    """
    Return the usual projection, giving the usual Rauzy fractal.
    """
    A = s.domain().alphabet()
    d = s.rauzy_fractal_projection()
    r = []
    for a in A:
        r.append(d[a])
    return matrix(r).transpose()

def conformize (s, V):
    """
    Compute a matrix M such that the Rauzy fractal is conformal for the projection M*V.
    """
    m = s.incidence_matrix()
    pi = m.charpoly()
    b = max(m.eigenvalues())
    pib = b.minpoly()
    P = pi/pib
    # basis of kernel of pib(m)
    mb = matrix(pib(m).right_kernel().basis()).transpose()
    # basis of kernel of P(m)
    mP = matrix(P(m).right_kernel().basis()).transpose()
    # diagonalize by block
    p = matrix(mb.columns()+mP.columns()).transpose()
    # projection on ker(pib(m)) along ker(P(m))
    PI = p*diagonal_matrix([1 for _ in range(mb.ncols())]+[0 for _ in range(mP.ncols())])*p^(-1)
    # compute the usual projection associated to V
    Vu = matrix(V*PI, ring=RR)
    # usual projection
    V0 = usual_projection(s)
    # matrix of conformization
    k = V.nrows()
    return (V0*p).submatrix(0,0,k,k)*(Vu*p).submatrix(0,0,k,k)^(-1)

def Vp (s, lvp, conform=True, prec=52, verb=False):
    """
    Compute the projection V' giving a extension of a torus translation
    for the susbtitution s and for independant eigenvalues lvp of the subshift.
    """
    m = s.incidence_matrix()
    # compute order of vp 0
    pi = m.charpoly()
    k = 0
    x = pi.variables()[0]
    while pi(0) == 0:
        pi /= x
        k += 1
    if verb:
        print("eigenvalue 0 of order {}".format(k))
    if k == 0:
        k = 1
    # compute line vectors w
    l1 = vector([1 for _ in range(m.ncols())])
    lv = []
    RR = RealField(20000)
    lgv, _ = generalized_eigenvectors(m)
    for vp in lvp:
        for i in range(100):
            va = vector([round(RR(t)) for t in vp*l1*m^(k*i)])
            if i > 0:
                if va == rva*m^k:
                    v = va - vp*l1*m^(k*i)
                    # check that v is orthogonal to every generalized eigenvector for eigenvalues >=1
                    ok = True
                    for v2 in lgv:
                        if v*v2 != 0:
                            ok = False
                            break
                    if ok:
                        break
            rva = va
        else:
            raise RuntimeError("Could not find n such that a(n+k) = a(n) m^k and such that (vp(1,...,1)-a(n))m^N ---> 0.")
        # compute an inverse of va by m^(k*i)
        mi = pseudoinv(m^k)
        w = va*(mi^i)
        # check it is in Z
        for t in w:
            if t not in ZZ:
                raise RuntimeError("w is not in Z !")
        lv.append(vp*l1 - w)
    V = matrix(lv)
    if conform:
        M = conformize(s, V)
        RR = RealField(prec)
        return M*matrix(V, ring=RR)
    return V

def rauzy_fractal_points (u, V, n=200000, prec=52, exchange=False, translate=None, modZ=False, verb=False):
    """
    Compute points of the Rauzy fractal of an infinite word u for the projection V.
    """
    d = dict()
    #ind = dict()
    proj = dict()
    if translate is None:
        translate = zero_vector(V.nrows())
    else:
        try:
            translate = [vector(v) for v in translate]
        except:
            raise ValueError("translate argument must be a list of translation vectors")
    A = u.parent().alphabet()
    if verb:
        print("alphabet: {}".format(A))
    I = identity_matrix(len(A))
    RR = RealField(prec)
    for i,a in enumerate(A):
        #ind[a] = i
        proj[a] = vector([RR(t) for t in V*I[i]])
        d[a] = []
    u = iter(u)
    ra = next(u)
    #V = matrix(V, ring=RR)
    S = zero_vector(V.nrows())
    #S = zero_vector(len(A))
    for t in translate:
        for _ in range(n):
            S += proj[ra]
            if modZ:
                for i in range(len(S)):
                    S[i] -= round(S[i])
            #S[ind[ra]] += 1
            a = next(u)
            if exchange:
                d[ra].append(S+t) #V*S)
            else:
                d[a].append(S+t) #V*S)
            ra = a
    RR = RealField(52)
    return d

def rauzy_fractal_plot(u, V, n=200000, prec=52, exchange=False, translate=None, modZ=False, M=None, verb=False):
    """
    Plot the generalized Rauzy fractal of the substitution s for the projection V.
    
    INPUT:
    
    - ``u`` - an infinite word, or a substitution
    
    - ``V`` - a projection map
    
    - ``n`` - number of plotted points
    
    - ``prec`` - int (default: ``52``) precision of real numbers used for the computation
    
    - ``exchange`` - bool (default: ``False``) if True, draw the domains after exchange
    
    - ``translate`` - list (default: ``None``) list of vectors of translation
    
    - ``modZ`` - bool (default: ``False``) if True, compute modulo Z
    
    - ``M`` - matrix (default: ``None``) change of basis (used only if modZ is True)
    
    """
    if type(u) == WordMorphism:
        u = u.periodic_points()[0][0]
    if verb:
        print("u = {}".format(u))
    d = rauzy_fractal_points(u, V, n, prec, exchange, translate, modZ, verb)
    from matplotlib import cm
    colormap = cm.__dict__['hsv']
    col_dict = {}
    A = u.parent().alphabet()
    for i, a in enumerate(A):
        col_dict[a] = colormap(float(i)/float(len(A)))[:3]
    G = Graphics()
    if V.nrows() == 2 or V.nrows() == 3:
        for a in A:
            if modZ and M is not None:
                G += points([M*v for v in d[a]], color=col_dict[a], size=1)
            else:
                G += points(d[a], color=col_dict[a], size=1)
    elif V.nrows() == 1:
        for a in A:
            G += plot([x[0] for x in d[a]], color=col_dict[a], thickness=1)
    else:
        raise ValueError("Can't plot in dimension {}.".format(V.nrows()))
    G.set_aspect_ratio(1)
    return G

def print_info (s):
    """
    Print informations about the substitution s:
    - primitivity
    - incidence matrix
    - characteristic polynomial
    - eigenvalues of the incidence matrix
    - pseudo-determinant
    - left-proprifed subtitution
    - basis of the Z-module of alpha such that exp(2 i pi alpha) is an eigenvalue of the subshift
    - indicate if it is weakly mixing, or if it finds that it is a finite extension of a torus translation.
    """
    print(s)
    print("primitive? {}".format(s.is_primitive()))
    m = s.incidence_matrix()
    print("incidence matrix:")
    print(m)
    print("charpoly: {}".format(m.charpoly().factor()))
    print("eigenvalues: {}".format(m.eigenvalues()))
    print("pseudo-det: {}".format(ZZ(pseudodet(m))))
    sp = proprify(s, None, False, False)
    print("left-proprified substitution:")
    if len(sp.domain().alphabet()) > 40:
        print("substitution on an alphabet of {} letters.".format(len(sp.domain().alphabet())))
    else:
        print(sp)
    lvp = eigenvalues_proper(sp.incidence_matrix(), verb=False)
    m = s.incidence_matrix()
    b = max(m.eigenvalues())
    if b.degree() == len(lvp):
        K.<b> = NumberField(b.minpoly(), embedding=b)
        lvp = [K(t) for t in lvp]
    print("Eigenvalues of the subshift:")
    print(lvp)
    print("where {}={} is root of {}.".format(lvp[0].parent().gen(), lvp[0].parent().gen().n(), lvp[0].parent().gen().minpoly()))
    if len(lvp) == 1 and abs(lvp[0]) == 1:
        print("The subshift is weakly mixing.")
    if len(lvp) == b.minpoly().degree():
        # test if condition of Thm 1.1 is satisfied
        lv, lvp = generalized_eigenvectors(sp.incidence_matrix(), verb=False)
        if len([v for v in lv if sum(v) != 0]) == 1:
            print("The subshift is a finite extension of a minimal translation on the torus T^{}.".format(b.minpoly().degree()-1))

# Computation of proprification

# Exemple of substitution due to Timo Jolivet
s = WordMorphism('1->213,2->4,3->5,4->1,5->21')
show(s)
print("cube of the substitution:")
s = s^3
show(s)
print("fixed point:")
show(s.fixed_point('1'))
su, R = ReturnSubstitution(s, '1', getrw=True)
print("set of return words on 1, with the correspondance: {}".format(R))
print("return substitution:")
show(su)
print("proprified substitution:")
show(proprify(s))

cube of the substitution:

fixed point:

set of return words on 1, with the correspondance: {word: 142: 0, word: 1352: 1, word: 13: 2}
return substitution:

proprified substitution:

# Computation of eigenvalues of the subshift

# Exemple of Ferenczi-Mauduit-Nogueira (Example 2 in [FMN])
s = WordMorphism('a->abdd,b->bc,c->d,d->a')

show(s)
m = s.incidence_matrix()
print("incidence matrix:")
show(m)
print("eigenvalues of the incidence matrix:")
show(m.eigenvalues())
print("Computation of eigenvalues of the subshift")
eigenvalues_proper(m, True)

incidence matrix:

eigenvalues of the incidence matrix:

Computation of eigenvalues of the subshift
Perron eigenvalue: 2.132241882311901?
x^4 - 2*x^3 - x^2 + 2*x - 1 1
-1.132241882311901?
vp = -1.132241882311901?, multiplicity = 1
2.132241882311901?
vp = 2.132241882311901?, multiplicity = 1
0.500000000000000? - 0.4052327261871813?*I
0.500000000000000? + 0.4052327261871813?*I
Basis of generalized eigenspaces: [((1, b^3 - b^2 - 2*b, b^2 - b - 2, b - 1), -1.132241882311901?), ((1, b^3 - b^2 - 2*b, b^2 - b - 2, b - 1), 2.132241882311901?)]
Basis of the orthogonal of w: [(-6/41*b^3 + 9/41*b^2 - 17/41*b + 7/41, 18/41*b^3 - 27/41*b^2 - 31/41*b + 20/41, -14/41*b^3 + 21/41*b^2 + 15/41*b - 11/41, 2/41*b^3 - 3/41*b^2 + 33/41*b - 16/41)]
Matrix whose integer left kernel is the set of possible w: [  7 -17   9  -6]
[ 20 -31 -27  18]
[-11  15  21 -14]
[-16  33  -3   2]
Basis of the Z-module where w lives: [
(1, 1, 1, 1),
(0, 3, 4, 1)
]
List of (non-independant) eigenvalues: [1, -b^2 + b + 4]
Corresponding matrix: [ 1  0  0  0]
[ 4  1 -1  0]
den = 1
matrix after BKZ: [ 1  0  0  0]
[ 0  1 -1  0]
Find a small NumberField to express the result.
Splitting field: Number Field in b with defining polynomial x^2 + 2*x - 1

[1, b]

# do the computation in the NumberField of the Perron eigenvalue
lv = m.eigenvectors_right()
print(lv[0][0])
v1 = lv[0][1][0]
v1 = v1/sum(v1)
print(lv[1][0])
v0 = lv[1][1][0]
v0 = v0/sum(v0)
vd = v1-v0
b = max(m.eigenvalues())
pi = b.minpoly()
K.<b> = NumberField(pi, embedding=b)
v = vector([K(t) for t in vd])
v0 = vector([K(t) for t in v0])
41*v

-1.132241882311901?
2.132241882311901?

(-6*b^3 + 9*b^2 - 17*b + 7, 18*b^3 - 27*b^2 - 31*b + 20, -14*b^3 + 21*b^2 + 15*b - 11, 2*b^3 - 3*b^2 + 33*b - 16)

ms = 41*matrix([list(t) for t in v])
ms = matrix(ms, ring=ZZ)
ms

[  7 -17   9  -6]
[ 20 -31 -27  18]
[-11  15  21 -14]
[-16  33  -3   2]

ms.left_kernel()

Free module of degree 4 and rank 2 over Integer Ring
Echelon basis matrix:
[1 1 1 1]
[0 3 4 1]

la = [w*v0 for w in ms.kernel().basis()]
la

[1, -b^2 + b + 4]

(la[1]-3).minpoly()

x^2 - 2

# Exemple 9.1 due to Ferenczi-Mauduit-Nogueira (Example 1 in [FMN])
s = WordMorphism('a->abbbccccccccccdddddddd,b->bccc,c->d,d->a')
print_info(s)

a->abbbccccccccccdddddddd, b->bccc, c->d, d->a
primitive? True
incidence matrix:
[ 1  0  0  1]
[ 3  1  0  0]
[10  3  0  0]
[ 8  0  1  0]
charpoly: x^4 - 2*x^3 - 7*x^2 - 2*x + 1
eigenvalues: [-1.492066037647537?, -0.6702116225208424?, 0.2559980601477473?, 3.906279600020632?]
pseudo-det: 1
left-proprified substitution:
substitution on an alphabet of 126 letters.
Eigenvalues of the subshift:
[1, b]
where b=-2.41421356237309 is root of x^2 + 2*x - 1.

# Exemple 9.2
s = WordMorphism('1->2,2->3,3->14,4->5,5->1425')
print_info(s)

1->2, 2->3, 3->14, 4->5, 5->1425
primitive? True
incidence matrix:
[0 0 1 0 1]
[1 0 0 0 1]
[0 1 0 0 0]
[0 0 1 0 1]
[0 0 0 1 1]
charpoly: x^2 * (x^3 - x^2 - x - 1)
eigenvalues: [0, 0, 1.839286755214161?, -0.4196433776070806? - 0.6062907292071993?*I, -0.4196433776070806? + 0.6062907292071993?*I]
pseudo-det: 1
left-proprified substitution:
0->01234, 1->5,6,7,8,9,10,11, 10->12,13,14,15, 11->0,1,2,3,4,0,1,2,3,4,12,13,14,15, 12->01234, 13->5,6,7,8,9,10,11, 14->12,13,14,15, 15->0,1,2,3,4,12,13,14,15, 2->12,13,14,15, 3->01234, 4->0,1,2,3,4,12,13,14,15, 5->01234, 6->5,6,7,8,9,10,11, 7->12,13,14,15, 8->01234, 9->5,6,7,8,9,10,11
Eigenvalues of the subshift:
[15/11*b^2 + 13/11*b + 10/11, 17/11*b^2 + 14/11*b + 4/11, -2/11*b^2 - 1/11*b - 5/11]
where b=1.83928675521416 is root of x^3 - x^2 - x - 1.
The subshift is a finite extension of a minimal translation on the torus T^2.

# proprify
sp = proprify(s)
# compute eigenvalues of the subshift
m = sp.incidence_matrix()
lvp = eigenvalues_proper(m)
lvp

[2/11*b^2 + 1/11*b + 1/11, 3/11*b^2 - 4/11*b - 4/11, -1/11*b^2 + 5/11*b - 6/11]

# choose independant eigenvalues
ivp = vector([lvp[0], lvp[1]])
ivp

(2/11*b^2 + 1/11*b + 1/11, 3/11*b^2 - 4/11*b - 4/11)

# Compute the projection V' associated to eigenvalues ivp (and conformize in order to do a nicer plot)
V = Vp(sp, ivp, conform=True)
#print(V)

# plot the Rauzy fractal R'
g = rauzy_fractal_plot(sp, V, n=200000)
g.show(figsize=15)

# plot the Rauzy fractal R' in a fundamental domain of T^2

V = Vp(sp, ivp, conform=False)
M = conformize(sp, V)
P = matrix([[1,0],[1,1]])*matrix([[1,9],[0,1]])*matrix([[1,0],[-1,1]]) # change of basis in SL(2,Z)
print(M*P)
V = matrix(V, ring=RR)
g = rauzy_fractal_plot(sp, P^(-1)*V, modZ=True, M=M*P, n=100000)
g

[ -2.19148788395506 -0.228155493655208]
[ 0.508851778832337  -1.11514250804019]

# translate
N = 1
tr = [(i,j) for i in range(-N,N+1) for j in range(-N,N+1)]
g = rauzy_fractal_plot(sp, P^(-1)*V, n=200000, translate=tr, modZ=True, M=M*P)
g.show(figsize=25)

# plot the domain exchange
rauzy_fractal_plot(sp, usual_projection(sp))

# plot the domain exchange after exchange
rauzy_fractal_plot(sp, usual_projection(sp), exchange=True)

# plot the usual Rauzy fractal R in a fundamental domain of T^2
# we see that the domain exchange does not give directly a tile
# the map psi is non-trivial
V = usual_projection(sp)
g = rauzy_fractal_plot(sp, P^(-1)*M^(-1)*V, modZ=True, M=M*P, n=100000)
g

# plot the Rauzy fractal R' directly for the substitution s
# it seems to give directly a domain exchange
lvp = eigenvalues_proper(s.incidence_matrix())
ivp = lvp[:2]
V = Vp(s, ivp)
rauzy_fractal_plot(s, V)

# Exemple 9.3
s = WordMorphism('a->Ab,b->A,A->aB,B->a')
print_info(s)

A->aB, B->a, a->Ab, b->A
primitive? True
incidence matrix:
[0 0 1 1]
[1 0 0 0]
[1 1 0 0]
[0 0 1 0]
charpoly: (x^2 - x - 1) * (x^2 + x + 1)
eigenvalues: [-0.618033988749895?, 1.618033988749895?, -0.50000000000000000? - 0.866025403784439?*I, -0.50000000000000000? + 0.866025403784439?*I]
pseudo-det: -1
left-proprified substitution:
0->0123456, 1->789, 10->0,1,2,3,4,5,6,10, 11->0123456, 12->789, 13->10, 14->11,12,13,14,10,15,16,22,23,24,17,18,19,20,21, 15->0123456, 16->7,8,9,10,15,16, 17->0123456, 18->789, 19->10, 2->10, 20->11,12,13,14, 21->10,15,16,17,18,19,20,21,10,15,16,22,23,24,17,18,19,20,21, 22->0123456, 23->789, 24->10,11,12,13,14,10,15,16,22,23,24, 3->11,12,13,14, 4->10, 5->15,16, 6->17,18,19,20,21,10,15,16,22,23,24,0,1,2,3,4,5,6,10,15,16,22,23,24,17,18,19,20,21, 7->0123456, 8->789, 9->10,15,16,22,23,24,17,18,19,20,21
Eigenvalues of the subshift:
[b + 1, -b]
where b=1.61803398874989 is root of x^2 - x - 1.
The subshift is a finite extension of a minimal translation on the torus T^1.

sp = proprify(s)
sp

WordMorphism: 0->0123456, 1->789, 10->0,1,2,3,4,5,6,10, 11->0123456, 12->789, 13->10, 14->11,12,13,14,10,15,16,22,23,24,17,18,19,20,21, 15->0123456, 16->7,8,9,10,15,16, 17->0123456, 18->789, 19->10, 2->10, 20->11,12,13,14, 21->10,15,16,17,18,19,20,21,10,15,16,22,23,24,17,18,19,20,21, 22->0123456, 23->789, 24->10,11,12,13,14,10,15,16,22,23,24, 3->11,12,13,14, 4->10, 5->15,16, 6->17,18,19,20,21,10,15,16,22,23,24,0,1,2,3,4,5,6,10,15,16,22,23,24,17,18,19,20,21, 7->0123456, 8->789, 9->10,15,16,22,23,24,17,18,19,20,21

V2

[-0.170820393249936  0.723606797749977   1.17082039324993  0.276393202250020   1.17082039324993  0.723606797749977  -3.74852915724959 -0.170820393249936  0.723606797749977 -0.618033988749893 -0.170820393249936 -0.170820393249936  0.723606797749977   1.17082039324993  -1.51246117974981 -0.170820393249936  0.276393202250020 -0.170820393249936  0.723606797749977   1.17082039324993  0.276393202250020  -1.95967477524976 -0.170820393249936  0.723606797749977 -0.618033988749893]

# plot psi for sp
# we see it is a translation by pieces
lvp = eigenvalues_proper(sp.incidence_matrix())
print(lvp, lvp[0].parent())
V1 = usual_projection(sp)
V2 = Vp(sp, [QQbar(lvp[0])])
V = matrix([V1[0],V2[0]], ring=RR)
rauzy_fractal_plot(sp, V, 5000)

[1/3*b + 1/3, -1/3*b + 2/3] Number Field in b with defining polynomial x^2 - 7*x + 1 with b = 6.854101966249684?

# plot a domain exchange conjugate to the subshift
V = usual_projection(sp)
rauzy_fractal_plot(sp, V, 1000)

# ... and after exchange
rauzy_fractal_plot(sp, V, 1000, exchange=True)

# plot R'
V = Vp(sp, [QQbar(lvp[0])], conform=False)
g = rauzy_fractal_plot(sp, V, 1000)
g.show(figsize=25)

# plot R' in T^2
# we see that it covers two times

M = matrix([[1]])
g = rauzy_fractal_plot(sp, V, modZ=True, M=M, n=500)
g.show(figsize=25)

g = rauzy_fractal_plot(sp, V, modZ=True, M=M, n=5000) # , n=500000
g.show(figsize=25) #figsize

# Exemple 9.4
s = WordMorphism("1->11116,2->1,3->1111112,4->1111113,5->466,6->566")
print_info(s)

1->11116, 2->1, 3->1111112, 4->1111113, 5->466, 6->566
primitive? True
incidence matrix:
[4 1 6 6 0 0]
[0 0 1 0 0 0]
[0 0 0 1 0 0]
[0 0 0 0 1 0]
[0 0 0 0 0 1]
[1 0 0 0 2 2]
charpoly: (x^2 - 4*x - 1) * (x^2 - x - 1)^2
eigenvalues: [-0.2360679774997897?, 4.236067977499789?, -0.618033988749895?, -0.618033988749895?, 1.618033988749895?, 1.618033988749895?]
pseudo-det: -1
left-proprified substitution:
substitution on an alphabet of 68 letters.
Eigenvalues of the subshift:
[1, -45/4*b - 11/4]
where b=4.23606797749979 is root of x^2 - 4*x - 1.

sp = proprify(s)
sp

WordMorphism: 0->00012000120001200034567, 1->0, 10->0, 11->12, 12->0, 13->0, 14->0, 15->12, 16->0, 17->0, 18->0, 19->12, 2->0,0,1,2,0,0,0,1,2,0,0,0,1,2,0,0,0,8,9,10,11,12,13,14,15,16,17,18,19,20,21, 20->0, 21->0,0,8,9,10,11,12,13,14,15,16,17,18,19,20,21,0,0,0,0,0,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,..., 22->0, 23->0, 24->0, 25->12, 26->0, 27->0, 28->0, 29->12, 3->0, 30->0, 31->0, 32->0, 33->12, 34->0, 35->0, 36->0, 37->34567, 38->0, 39->0, 4->0, 40->0, 41->12, 42->0, 43->0, 44->0, 45->12, 46->0, 47->0,0,1,2,0,0,0,1,2,0,0,0,1,2,0,0,0,1,2,0,0,0,0,0,0,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,..., 48->0, 49->0, 5->0, 50->0, 51->12, 52->0, 53->0, 54->0, 55->12, 56->0, 57->0, 58->0, 59->12, 6->12, 60->0, 61->0, 62->0, 63->34567, 64->0, 65->0, 66->0, 67->1,2,0,0,0,1,2,0,0,0,1,2,0,0,0,1,2,0,0,0,1,2,0,0,0,1,2,0,0,0,1,2,0,0,0,0,0,48,49,50,..., 7->0,0,0,1,2,0,0,0,1,2,0,0,0,8,9,10,11,12,13,14,15,16,17,18,19,20,21,0,0,0,0,0,22,23,24,25,26,27,28,29,..., 8->0, 9->0

# plot psi for sp
lvp = eigenvalues_proper(sp.incidence_matrix())
print(lvp, lvp[0].parent())
V1 = usual_projection(sp)
V2 = Vp(sp, [QQbar(lvp[1])])
V = matrix([V1[0],V2[0]], ring=RR)
rauzy_fractal_plot(sp, V, 5000)

[1, -45/16*b + 1/16] Number Field in b with defining polynomial x^2 - 18*x + 1 with b = 17.94427190999916?

# Exemple 9.5
s = WordMorphism('1->1116,2->1,3->2,4->3,5->1146,6->566')
print_info(s)

1->1116, 2->1, 3->2, 4->3, 5->1146, 6->566
primitive? True
incidence matrix:
[3 1 0 0 2 0]
[0 0 1 0 0 0]
[0 0 0 1 0 0]
[0 0 0 0 1 0]
[0 0 0 0 0 1]
[1 0 0 0 1 2]
charpoly: (x^2 - x - 1) * (x^4 - 4*x^3 + 2*x^2 - x + 1)
eigenvalues: [-0.618033988749895?, 1.618033988749895?, 0.7512735895760328?, 3.484794387923757?, -0.1180339887498949? - 0.6066580492747911?*I, -0.1180339887498949? + 0.6066580492747911?*I]
pseudo-det: -1
left-proprified substitution:
0->001200120034567, 1->0, 10->12, 11->0, 12->0, 13->12, 14->0, 15->0, 16->23,24,25,26,27,28,0,8,9,10,11,12,13,14,15,16,0,0,1,2,0,0,17,18,19,20,21,22,0,8,9,10,11,12,13,14,15,16,0,8,..., 17->0, 18->0, 19->12, 2->0,1,2,0,0,1,2,0,0,3,4,5,6,7,0,8,9,10,11,12,13,14,15,16, 20->0, 21->0, 22->1,2,0,0,3,4,5,6,7,0,8,9,10,11,12,13,14,15,16,0,0,0,1,2,0,0,17,18,19,20,21,22,0,8,9,10,11,12,13,14,..., 23->0, 24->0, 25->12, 26->0, 27->0, 28->1,2,0,0,3,4,5,6,7,0,8,9,10,11,12,13,14,15,16,0,0,1,2,0,0,17,18,19,20,21,22,0,8,9,10,11,12,13,14,15,..., 3->0, 4->0, 5->12, 6->0, 7->0,1,2,0,0,3,4,5,6,7,0,8,9,10,11,12,13,14,15,16,0,0,1,2,0,0,17,18,19,20,21,22,0,8,9,10,11,12,13,14,..., 8->0, 9->0
Eigenvalues of the subshift:
[1, b]
where b=-105.227937358744 is root of x^2 + 83*x - 2339.

# convert the eigenvalues into the NumberField of the golden number
pi = x^2-x-1
b = max(pi.roots(ring=AA))[0]
K.<b> = NumberField(pi, embedding=b)
lvp = eigenvalues(s)
[K(t) for t in lvp]

[1, -57*b - 13]

# Exemple 9.6
s = WordMorphism('1->15,2->2122,3->122,4->13,5->14122')
print_info(s)

1->15, 2->2122, 3->122, 4->13, 5->14122
primitive? True
incidence matrix:
[1 1 1 1 2]
[0 3 2 0 2]
[0 0 0 1 0]
[0 0 0 0 1]
[1 0 0 0 0]
charpoly: (x^2 - x - 1) * (x^3 - 3*x^2 - x - 1)
eigenvalues: [-0.618033988749895?, 1.618033988749895?, 3.382975767906238?, -0.1914878839531188? - 0.5088517788327380?*I, -0.1914878839531188? + 0.5088517788327380?*I]
pseudo-det: -1
left-proprified substitution:
0->01, 1->23456, 10->23, 11->12,13,14,15,4,5,6, 12->9,10,11, 13->12,13,14,15, 14->12,13,14,15, 15->456, 2->01, 3->78, 4->9,10,11, 5->12,13,14,15, 6->456, 7->01, 8->456, 9->01
Eigenvalues of the subshift:
[1]
where 0=0.000000000000000 is root of x.
The subshift is weakly mixing.

# it is not clear if the substitution satisfies the strong coincidence...
u = (s^8)('1')
v = (s^8)('2')
for i in range(len(u)):
    if u[:i].abelian_vector() == v[:i].abelian_vector() and u[i] == v[i]:
        print("ok")
        break

# ... but the substitution seems to define directly a domain exchange
rauzy_fractal_plot(s, usual_projection(s))

rauzy_fractal_plot(s, usual_projection(s), exchange=True)

# proprify
sp = proprify(s)
sp

WordMorphism: 0->01, 1->23456, 10->23, 11->12,13,14,15,4,5,6, 12->9,10,11, 13->12,13,14,15, 14->12,13,14,15, 15->456, 2->01, 3->78, 4->9,10,11, 5->12,13,14,15, 6->456, 7->01, 8->456, 9->01

# domain exchange
g = rauzy_fractal_plot(sp, usual_projection(sp)) #, 300000)
g.show() #figsize=15)

g = rauzy_fractal_plot(sp, usual_projection(sp), exchange=True) #700000
g.show() #figsize=15)

# Exemple 9.7
s = WordMorphism('1->16,2->122,3->12,4->3,5->124,6->15')
print_info(s)

1->16, 2->122, 3->12, 4->3, 5->124, 6->15
primitive? True
incidence matrix:
[1 1 1 0 1 1]
[0 2 1 0 1 0]
[0 0 0 1 0 0]
[0 0 0 0 1 0]
[0 0 0 0 0 1]
[1 0 0 0 0 0]
charpoly: (x^3 - 3*x^2 + 2*x - 1) * (x^3 - x - 1)
eigenvalues: [2.324717957244746?, 0.3376410213776270? - 0.5622795120623013?*I, 0.3376410213776270? + 0.5622795120623013?*I, 1.324717957244746?, -0.6623589786223730? - 0.5622795120623013?*I, -0.6623589786223730? + 0.5622795120623013?*I]
pseudo-det: 1
left-proprified substitution:
0->01, 1->2301456, 10->11,12,13,11,12,13,0,1,11,12,13,11,12,13,0,1,11,12,13, 11->01, 12->23, 13->0,1,11,12,13,11,12,13,0,1,11,12,13,11,12,13, 14->01, 15->2,3,0,1,11,12,13,11,12,13, 2->01, 3->2,3,0,1,7,8,9,10, 4->01, 5->23, 6->0,1,11,12,13,11,12,13,14,15, 7->01, 8->23, 9->01
Eigenvalues of the subshift:
[-1, -3*b^2 + 3*b - 2, -9*b^2 + 6*b - 4]
where b=2.32471795724475 is root of x^3 - 3*x^2 + 2*x - 1.

# the square of s is left-proper, so there is no need to proprify
lvp = eigenvalues_proper(s.incidence_matrix())
lvp, lvp[0].parent()

([-1, -3*b, -3*b^2],
 Number Field in b with defining polynomial x^3 - 3*x^2 + 2*x - 1 with b = 2.324717957244746?)

# plot the Rauzy fractal R'
V = Vp(s, lvp[1:])
g = rauzy_fractal_plot(s, V, 2000000)
g.show() #figsize=20)

# plot the domain exchange before...
rauzy_fractal_plot(s, usual_projection(s))

# ... and after
rauzy_fractal_plot(s, usual_projection(s), exchange=True)

# Exemple 9.8
s = WordMorphism('1->114,2->122,3->2,4->13')
print_info(s)

1->114, 2->122, 3->2, 4->13
primitive? True
incidence matrix:
[2 1 0 1]
[0 2 1 0]
[0 0 0 1]
[1 0 0 0]
charpoly: (x^2 - 3*x + 1) * (x^2 - x - 1)
eigenvalues: [0.3819660112501051?, 2.618033988749895?, -0.618033988749895?, 1.618033988749895?]
pseudo-det: -1
left-proprified substitution:
0->012, 1->0, 10->8,9,10,8,9,10, 2->1234, 3->0, 4->567, 5->0, 6->12, 7->3,4,8,9,10, 8->0, 9->12
Eigenvalues of the subshift:
[-1, -2*b]
where b=2.61803398874989 is root of x^2 - 3*x + 1.

# there is no need to proprify since the square of s is proper
print("eigenvalues of the incidence matrix:")
print(s.incidence_matrix().eigenvalues())
print("check condition of Thm 1.1")
for vp, lv, _ in s.incidence_matrix().eigenvectors_right():
    if abs(vp) < 1:
        continue
    print()
    print(vp)
    for v in lv:
        print(v)
        print(sum(v))

eigenvalues of the incidence matrix:
[0.3819660112501051?, 2.618033988749895?, -0.618033988749895?, 1.618033988749895?]
check condition of Thm 1.1

2.618033988749895?
(1, 0.2360679774997897?, 0.1458980337503155?, 0.3819660112501051?)
1.763932022500211?

1.618033988749895?
(1, -1, 0.3819660112501051?, 0.618033988749895?)
1

# plot psi
# we see it does not seems to be finite-to-one
lvp = eigenvalues_proper(s.incidence_matrix())
print(lvp, lvp[0].parent())
V1 = usual_projection(s)
V2 = Vp(s, [QQbar(lvp[1])])
V = matrix([V1[0],V2[0]], ring=RR)
rauzy_fractal_plot(s, V, 5000)

[1, -2*b] Number Field in b with defining polynomial x^2 - 3*x + 1 with b = 2.618033988749895?

# Rauzy fractal R'
# we see that it covers many times (but a finite number of times)
V2 = Vp(s, [QQbar(lvp[1])], conform=False)
g = rauzy_fractal_plot(s, V2, 2000)
g.show(figsize=15)

# domain exchange
rauzy_fractal_plot(s, V1, 1000)

rauzy_fractal_plot(s, V1, 1000, exchange=True)

# Exemple 9.9 (Timo Jolivet)
s = WordMorphism('1->213,2->4,3->5,4->1,5->21')
print_info(s)

1->213, 2->4, 3->5, 4->1, 5->21
primitive? True
incidence matrix:
[1 0 0 1 1]
[1 0 0 0 1]
[1 0 0 0 0]
[0 1 0 0 0]
[0 0 1 0 0]
charpoly: (x^2 + x + 1) * (x^3 - 2*x^2 + x - 1)
eigenvalues: [-0.50000000000000000? - 0.866025403784439?*I, -0.50000000000000000? + 0.866025403784439?*I, 1.754877666246693?, 0.1225611668766537? - 0.744861766619745?*I, 0.1225611668766537? + 0.744861766619745?*I]
pseudo-det: 1
left-proprified substitution:
0->012, 1->3456, 2->012345678, 3->012, 4->3456, 5->012, 6->78012345678, 7->012, 8->345601278
Eigenvalues of the subshift:
[4*b^2 - b + 2, -b^2 - 1, b^2]
where b=1.75487766624669 is root of x^3 - 2*x^2 + x - 1.
The subshift is a finite extension of a minimal translation on the torus T^2.

# Rauzy fractal of s with overlaps
op = .3
s.rauzy_fractal_plot(opacity={'1':op,'2':op,'3':op,'4':op,'5':op}, n=500000)

sp = proprify(s)
sp

WordMorphism: 0->012, 1->3456, 2->012345678, 3->012, 4->3456, 5->012, 6->78012345678, 7->012, 8->345601278

# domain exchange obtained with a proprification
g = rauzy_fractal_plot(sp, usual_projection(sp)) #, 800000)
g.show() #figsize=15)

g = rauzy_fractal_plot(sp, usual_projection(sp), exchange=True) #, 800000
g.show() #figsize=15)

# Rauzy fractal R'
lvp = eigenvalues_proper(sp.incidence_matrix())
V = Vp(sp, [lvp[0],lvp[1]])
g = rauzy_fractal_plot(sp, V)
g.show(figsize=20)

# Rauzy fractal R' in T^2

# choose a basis of the lattice
V = Vp(sp, [lvp[0],lvp[1]], conform=False)
M = conformize(sp, V)
P = matrix([[-7,-4],[2,1]]) # change of basis in SL(2,Z)
print(M)
V = matrix(V, ring=RR)
g = rauzy_fractal_plot(sp, P^(-1)*V, modZ=True, M=M*P) # , n=500000
g.show() #figsize=20)

[0.569840290998021  2.49444101849304]
[ 2.61428255736394  9.14998895077367]

# translate
N = 1
tr = [(i,j) for i in range(-N,N+1) for j in range(-N,N+1)]
g = rauzy_fractal_plot(sp, P^(-1)*V, n=200000, translate=tr, modZ=True, M=M*P)
g.show(figsize=25)

# Rauzy fractal R in T^2
# we see that the domain exchange does not give directly a tiling
# the map psi is non trivial

# choose a basis of the lattice
V = Vp(sp, [lvp[0],lvp[1]], conform=False)
M = conformize(sp, V)
P = matrix([[-7,-4],[2,1]]) # change of basis in SL(2,Z)
print(M)
V = matrix(V, ring=RR)
g = rauzy_fractal_plot(sp, P^(-1)*M^(-1)*usual_projection(sp), modZ=True, M=M*P) # , n=500000
g.show() #figsize=20)

[0.569840290998021  2.49444101849304]
[ 2.61428255736394  9.14998895077367]

# Example of Ferenczi-Mauduit-Nogueira
s = WordMorphism('a->abdd,b->bc,c->d,d->a')
print_info(s)

a->abdd, b->bc, c->d, d->a
primitive? True
incidence matrix:
[1 0 0 1]
[1 1 0 0]
[0 1 0 0]
[2 0 1 0]
charpoly: x^4 - 2*x^3 - x^2 + 2*x - 1
eigenvalues: [-1.132241882311901?, 2.132241882311901?, 0.500000000000000? - 0.4052327261871813?*I, 0.500000000000000? + 0.4052327261871813?*I]
pseudo-det: -1
left-proprified substitution:
0->012345, 1->6, 10->6,11,12,13,14, 11->012345, 12->6, 13->7,8,9,10, 14->11,12,13,14,11,12,13,14, 15->012345, 16->6, 17->7,8,9,10, 18->11,12,13,14, 19->15,16,17,18,19,20,21, 2->7,8,9,10, 20->6, 21->11,12,13,14, 3->11,12,13,14, 4->15,16,17,18,19,20,21, 5->6, 6->01234566, 7->012345, 8->6, 9->7,8,9,10
Eigenvalues of the subshift:
[1, -b + 2]
where b=-0.414213562373095 is root of x^2 - 2*x - 1.

# Rauzy fractal with overlaps
s = WordMorphism('a->Ab,b->Ac,c->A,A->aB,B->aC,C->a')
s.rauzy_fractal_plot()

sp = proprify(s)
sp

WordMorphism: 0->012345, 1->678, 10->6,7,8,9,10,11,12,13,14,9,10,15,16,17, 11->012345, 12->678, 13->9,10, 14->11,12,13,14,9,10,15,16,17,0,1,2,3,4,5,9,10,15,16,17,0,1,2,3,4,5,0,1,2,3,4,5,6,7,8,9,10, 15->012345, 16->678, 17->9,10,11,12,13,14,9,10,15,16,17,0,1,2,3,4,5,9,10,15,16,17,0,1,2,3,4,5, 2->9,10, 3->11,12,13,14, 4->9,10, 5->15,16,17,0,1,2,3,4,5,9,10,15,16,17,0,1,2,3,4,5,0,1,2,3,4,5,6,7,8,9,10,15,16,17,0,1,2,3,4,5,..., 6->012345, 7->678, 8->9,10,11,12,13,14,9,10,15,16,17,0,1,2,3,4,5,0,1,2,3,4,5,6,7,8,9,10, 9->012345

# domain exchange conjugate to the subshift
rauzy_fractal_plot(sp, usual_projection(sp))

rauzy_fractal_plot(sp, usual_projection(sp), exchange=True)

# eigenvalues
lvp = eigenvalues_proper(sp.incidence_matrix())
lvp

[1/8*b^2 - 1/4*b + 1/8, 1/4*b^2 - 1/4, -3/8*b^2 - 1/4*b - 3/8]

# Rauzy fractal R'
V = Vp(sp, [lvp[0],lvp[1]])
g = rauzy_fractal_plot(sp, V)
g.show(figsize=20)

# Rauzy fractal R' in T^2
# we see that it tiles
# the subshift is a measurably conjugate to a translation on T^2

# choose a basis of the lattice
V = Vp(sp, [lvp[0],lvp[1]], conform=False)
M = conformize(sp, V)
P = matrix([[1,0],[2,1]])*matrix([[1,2],[0,1]])*matrix([[1,0],[3,1]]) # change of basis in SL(2,Z)
print(M*P)
V = matrix(V, ring=RR)
g = rauzy_fractal_plot(sp, P^(-1)*V, modZ=True, M=M*P, n=200000) # , n=500000
g.show(figsize=15)

[  -0.999999999998863    0.647798871253414]
[2.01794136955868e-12     1.72143323723311]

sp.rauzy_fractal_plot()