SharedNumber wall.sagewsOpen in CoCalc
Calculates the number wall of a sequence and plot
def Dragon(n):
    """
    This function returns the nth element of the dragon sequence
    """
    if n == 0:
        return 0
    else:
        while Mod(n, 2) == 0:
        #while n%2 == 0:
            n = n / 2
    return ((n-1)/2)%2

def Pagoda(n):
    """
    This function returns the nth element of the pagoda sequence
    """
    return Dragon(n+1) - Dragon(n-1)

def ThueMorse(n):
    """
    This function returns the nth element of the Thue-Morse sequence
    """
    if n == 0:
        return 0
    elif  n < 0:
        return 1-ThueMorse(-n-1)
    elif n%2==0: #n is odd
        return ThueMorse(n/2)
    else: #n is odd
        return 1-ThueMorse((n-1)/2)

def NumberWall(seq,mlo,mhi,nlo,nhi,char=0):
    """
    This function calculates a rectangular part of dimensions [mlo,mhi]*[nlo,nhi] of the number wall of a sequence mod p
    """
    leng=len(seq)
    off = (leng-(nhi-nlo+mhi+mhi+1))// 2 # focus offset
    wal_width = nhi+mhi+3 - (nlo-mhi-2)
    wal = [ [1]* (nhi+mhi+3 - (nlo-mhi-2)) for m in range(mlo,mhi+1)]
    rebase_j= -nlo+mhi+2
    #long wall of empty entries, typed as for sequence, subscript base = 1
    for m in range(mlo,mhi+1):  # assign each entry in matrix order
        i = m-mlo
        wal_nlo = nlo-mhi-2 + max(0,m)
        wal_nhi = nhi+mhi+3 - max(0,m)
        #for n in range(nlo-mhi-2,nhi+mhi+3):
        for n in range(wal_nlo,wal_nhi):
            j = n+rebase_j  # rebased subscripts
            if m < -1 or ((n < nlo-mhi-1 or n > nhi+mhi+1) and m != -1):
                wal[i][j] = Mod(0,char)
            elif m == -1 or ((n == nlo-mhi-1 or n == nhi+mhi+1) and m > -1):
                wal[i][j] = Mod(1,char) # initial sentinel windows
            elif m == 0 and n > nlo-mhi-1 and n < nhi+mhi+1:
                wal[i][j] = Mod(seq[j-2+off],char) # initial sequence
            # compute entry ( m gt 0 and n gt nlo-mhi-1 and n lt nhi+mhi+1 )
            elif wal[i-2][j] != 0: # zero-free case ( d = 0,1 )
                wal[i][j] = (wal[i-1][j]^2 - wal[i-1][j-1] * wal[i-1][j+1]) / wal[i-2][j]
            elif wal[i-1][j] == 0: # window interior ( wal[i-2,j] eq 0 )
                p = 0
                while wal[i-p-1][j] == 0:
                    p += 1
                q = 1
                while wal[i-p][j-q] == 0:
                    q += 1
                d = 1
                while wal[i-p][j+d-q] == 0:
                    d += 1
                # window defic.  d  and pane offsets  p, q
                if p < d-1:
                    wal[i][j] = Mod(0,char) # window pane
                else:  # inner frame
                    wal[i][j] = (1 - Integer(Mod((d-1)*(d-q),2))*2) * wal[i-q][j-q] * wal[i-d+q][j+d-q] / wal[i-d][j+d-q-q]
            else:  # window exterior ( wal[i-2,j] eq 0 and wal[i-1,j] ne 0 )
                d = 1
                while wal[i-1-d][j] == 0:
                    d += 1
                q = 1
                while wal[i-d][j-q] == 0:
                    q += 1 #window defic.  d  and pane offset  q

                P = wal[i-d-1][j+d-q-q] / wal[i-d-1][j+d-q-q-1] # frame ratios
                Q = wal[i-q-1][j-q] / wal[i-q-2][j-q]
                R = wal[i-d+q-1][j+d-q] / wal[i-d+q][j+d-q]
                S = wal[i-1][j] / wal[i-1][j+1]

                # outer frame
                wal[i][j] = wal[i-1][j] / R * ( Q * wal[i-d-2][j+d-q-q] / wal[i-d-1][j+d-q-q]
                                               + (1-Integer(Mod(d-q,2))*2) * ( P * wal[i-q-1][j-q-1] / wal[i-q-1][j-q]
                                                                              - S * wal[i-d+q-1][j+d-q+1] / wal[i-d+q-1][j+d-q]
                                                                             )
                                              )
    #    print(wal[i])
    # short wall: prune partial and sentinel columns from left and right
    wal = [ [ wal[i][j] for j in range(mhi+2,nhi-nlo+mhi+3) ] for i in range(0,mhi-mlo+1) ]

    return wal

def WallDeficiencies(wal,mlo,mhi,nlo,nhi,d=10):
    """
    This function counts the deficiencies of a number wall in a rectangular part of dimensions [mlo,mhi]*[nlo,nhi]
    """
    defs=[0]*d

    for m in range (mlo,mhi): #for each top left corner in wall
        i=m-mlo
        for n in range (nlo,nhi):
            j=n-nlo
            if wal[i][j]!=0:
                if wal[i+1][j]==0 or wal[i][j+1]==0:
                    defs[1]+=1
                else: #nonzero deficiency
                    k=1
                    while m+k < mhi and wal[i+k][j+1]==0:
                        k=k+1
                    l=1
                    while n+l < nhi and wal[i+1][j+l]==0:
                        l=l+1
                    if m+k < mhi or n+l < nhi:
                        defs[max(k,l)+1]+=1
                    else:
                        defs[0]+=1 #broken window
    return defs



def SquareTiling(tab,mlo,mhi,nlo,nhi,tel,cid):
    """
    Return a tiling of the 2 dim table tab
    """
    import time
    start = time.time()
    #pass 0: initialise
    codes=[]
    tel2=tel//2
    cid2=cid//2
    end = time.time()
    print("SquareTiling() initial pass: parameters = ", mlo, mhi, nlo, nhi, tel, cid, end-start)
    #Pass 1: build codes list and mini-table of states
    stalen = 0  # current state count
    stab = [ [ 0 for n in range(nlo//cid,nhi//cid+1) ] for m in range(mlo//cid,mhi//cid+1)] # mini-table of states
    k=0
    for m in range(mlo+tel2,mhi-tel2+1):
        for n in range(nlo+tel2,nhi-tel2+1):
            if m%cid==cid2 and n%cid==cid2: # scan wall squares
                i = m - mlo
                j = n - nlo # rebase subscripts
                code = [ [ tab[i+k][j+l] for l in range(-tel2,tel2+1) ] for k in range(-tel2,tel2+1) ]
                if code in codes:
                    s=codes.index(code)+1 # state from code
                else:
                    codes.append(code) # set new code
                    stalen += 1
                    s=stalen
                i = m//cid - mlo//cid
                j = n//cid - nlo//cid # rebase subscripts
                stab[i][j] = s # update mini-table
    end = time.time()
    print("first pass: segment tiled, state count", stalen, " ", end-start)
    #print(stab)
    #Pass 2: build inflation table, scanning states in mini-table
    tetrads = [ [0, 0, 0, 0] for s in range(0,stalen)] # daughters indexed by state; later extendable to all tetrads
    for m in range(ceil(mlo/2)+tel2+cid2+1,floor(mhi/2)-tel2-cid2):
        for n in range(ceil(nlo/2)+tel2+cid2+1,floor(nhi/2)-tel2-cid2):
            if m%cid == cid2 and n%cid==cid2: #centre
                i = m//cid - mlo//cid
                j = n//cid - nlo//cid  # rebase subscripts
                s = stab[i][j]
                if s != 0: # for each set state
                    for h in range(0,4): # Update partial inflation or check consistency
                        k = h//2 + 2*m//cid - mlo//cid - 1
                        l = Integer(Mod(h,2)) + 2*n//cid - nlo//cid - 1 # rebase subscripts
                        t = stab[k][l]  # daughter
                        if t == 0:
                            print("missing tile", s, t, h, tetrads[s-1])
                            return [codes,stalen,stab,tetrads]
                        elif tetrads[s-1][h] == 0:
                            tetrads[s-1][h] = t
                        elif tetrads[s-1][h] !=t:
                            print("inconsistent", s, t, h, tetrads[s-1])
                            return [codes,stalen,stab,tetrads]
    end = time.time()
    print("second pass: inflation consistent", end-start)
    # Pass 3: build tetrad list, scanning states in mini-table
    quarter = [ False for s in range(0,stalen) ]  # tetrad in focal quarter
    tetlen = stalen
    for m in range(mlo+tel2,mhi-tel2-cid):
        for n in range(nlo+tel2,nhi-tel2-cid):
            if ( m%cid==cid2 and n%cid==cid2 ): # centre
                i = m // cid - mlo // cid
                j = n // cid - nlo // cid  # rebase subscripts
                tet = [ stab[i][j], stab[i][j+1], stab[i+1][j], stab[i+1][j+1] ]
                if 0 not in tet: # for each complete tetrad
                    if tet in tetrads:
                        s = tetrads.index(tet) + 1
                    else:             # new tetrad
                        tetrads.append(tet)
                        tetlen += 1
                        s = tetlen
                        quarter.append(False)
                    #if s==214:
                    #    print(tetrads,tet,tetlen,s,tetrads.index(tet))
                    #    return [codes,stalen,stab,tetrads]
                    if m+m < mhi-tel2-cid-1 and m+m > mlo+tel2 and n+n < nhi-tel2-cid-1 and n+n > nlo+tel2:
                        #print(s,len(quarter))
                        quarter[s-1] = True # within focal quarter
    if False in quarter:
        s = quarter.index(False)
        print("incomplete tetrad", s, tetrads[s]) # abort
        return [codes,stalen,stab,tetrads]
    end = time.time()
    print("third pass: vertices complete, tetrad count = ", tetlen, " ", end-start)
    # Pass 4: stabilise state indices
    b = 1.0*(mhi-mlo)
    c = 1.0*(nhi-nlo)
    dist = [ b+c for s in range(0,stalen) ] # minimum norm indexed by state
    for m in range(mlo+tel2,mhi-tel2+1):
        for n in range(nlo+tel2,nhi-tel2+1):
            if m%cid==cid2 and n%cid==cid2: # scan mini-table
                i = m//cid - mlo//cid
                j = n//cid - nlo//cid # rebase subscripts
                s = stab[i][j]
                if s!=0:
                    #dist[s-1] = min(dist[s-1], abs(m - cid) + abs(n) + 0.1*m/b + 0.1*n/b/c) # abs(m-cid) is used to adjust indices as in Magma for pagoda but not for dragon
                    dist[s-1] = min(dist[s-1], abs(m) + abs(n) + 0.1*m/b + 0.1*n/b/c) # update Manhatten distance from focus, row order, column order
    perm = [ s for s in range(1,stalen+1)]
    dist, perm = (list(t) for t in zip(*sorted(zip(dist, perm))))
    #ParallelSort(~dist, ~perm) #invert dist to perm
    codes = [codes[perm[i] - 1] for i in range(0,stalen)]
    tetrads1 = [[tetrads[perm[i] - 1][h] for h in range(0,4) ] for i in range(0,stalen)]
    tetrads1.extend([tetrads[s] for s in range(stalen,tetlen) ]) # permute state indices
    merp = [ s for s in range(1,stalen+1)]
    perm, merp = (list(t) for t in zip(*sorted(zip(perm,merp))))
    #ParallelSort(~perm, ~merp) # invert perm to merp
    tetrads = [ [ merp[tetrads1[t][h] - 1] for h in range(0,4)] for t in range(0,tetlen) ]# permute state values
    #efes=[0]
    merp = [0] + merp
    #for i in range(0,7):
     #   print(stab[i])
    stab = [ [ merp[s] for s in t ] for t in stab ]
    #print("made it")
    #for i in range(0,7):
     #   print(stab[i])
    print("final pass: states stabilised")
    end = time.time()
    print(end-start)
    return [codes,stalen,stab,tetrads]
# this run takes about 10min
import time
print("hello")
start = time.time()

char = 3
tel = 8
cid = 8     # tile edge length, centre distance
mhi = 1600
nd2 = 4800  # wall rows and cols/2
pad = floor((cid + tel)*5/2) #sentinel rows
print("Pagoda validation: rows, cols = ", mhi, 2*nd2+1)
pag = [Pagoda(n) for n in range (-(nd2+mhi),nd2+mhi+1)]
end = time.time()
print(end-start)
wal = NumberWall(pag,-pad,mhi,-nd2,nd2,char)
end = time.time()
print(end-start)
defs=WallDeficiencies(wal,-pad, mhi, -nd2, nd2)
print ("deficiencies: ", defs) # 1, 2 isolated zeros
end = time.time()
print(end-start)
[codes,stalen,stab,tetrads] = SquareTiling(wal, 2-pad, 2+mhi, -nd2, nd2, tel, cid) # edge overlap, focus 2 rows up
print("Pagoda construction: rows, cols = ", mhi, 2*nd2+1)
end = time.time()
print(end-start)

hello ('Pagoda validation: rows, cols = ', 1600, 9601) 0.453722000122 170.198223829 ('deficiencies: ', [1, 4609983, 6149634, 2305916, 0, 0, 0, 0, 0, 0]) 278.718138933 ('SquareTiling() initial pass: parameters = ', -38, 1602, -4800, 4800, 8, 8, 5.0067901611328125e-06) ('first pass: segment tiled, state count', 201, ' ', 44.72906804084778) ('second pass: inflation consistent', 53.48532009124756) ('third pass: vertices complete, tetrad count = ', 1241, ' ', 79.84820413589478) final pass: states stabilised 96.5376532078 ('Pagoda construction: rows, cols = ', 1600, 9601) 375.258188009
infls = [tetrads[s] for s in range(0,stalen)]
with open("infls_pag.txt", "w") as f:
    for s in infls:
        f.write(str(s) +"\n")

with open("codes_pag.txt", "w") as f:
    for s in codes:
        f.write(str(s) +"\n")

with open("prolongable_pag.txt", "w") as f:
    for s in [5,5,1,2]:
        f.write(str(s) +"\n")

# this run takes about 50min and might not be possible to run on the online server
import time
print("hello")
start = time.time()

char = 3
tel = 12
cid = 8     # tile edge length, centre distance
mhi = 6000
nd2 = 15000  # wall rows and cols/2
pad = floor((cid + tel)*5/2) #sentinel rows
print("Dragon validation: rows, cols = ", mhi, 2*nd2+1)
rpf = [Dragon(n) for n in range (-(nd2+mhi),nd2+mhi+1)]
end = time.time()
print(end-start)
wal = NumberWall(rpf,-pad,mhi,-nd2,nd2,char)
end = time.time()
print(end-start)
defs=WallDeficiencies(wal,-pad, mhi, -nd2, nd2)
print ("deficiencies: ", defs) # 1, 2, 3, 4 isolated zeros
end = time.time()
print(end-start)
[codes,stalen,stab,tetrads] = SquareTiling(wal, -pad, mhi, -nd2, nd2, tel, cid) # edge overlap, focus as a sequence
end = time.time()
print(end-start)
hello ('Dragon validation: rows, cols = ', 6000, 30001)
infls = [tetrads[s] for s in range(0,stalen)]
with open("infls_dragon.txt", "w") as f:
    for s in infls:
        f.write(str(s) +"\n")

with open("codes_dragon.txt", "w") as f:
    for s in codes:
        f.write(str(s) +"\n")

with open("prolongable_dragon.txt", "w") as f:
    for s in [1,2,3,4]:
        f.write(str(s) +"\n")
def AutomaticNdTiling(infls,codes,prolongable,m,n):
    codeLeng = int(sqrt(len(codes[0])))
    inflsLeng = int(sqrt(len(infls[0])))
    m_infls = m//codeLeng + 1
    n_infls = n//codeLeng + 1
    infls_tab = [[0] * n_infls for t in range(0,m_infls)]
    infls_tab[0][0] = prolongable
    for i in range(0,m_infls):
        for j in range(0,n_infls):
            infls_tab[i][j] = infls[infls_tab[i//inflsLeng][j//inflsLeng]-1][i%inflsLeng * inflsLeng + j%inflsLeng]
    coded_tab = [[0] * n for i in range(0,m)]
    for i in range(0,m):
        temp_infls = infls_tab[i//codeLeng]
        temp_row = i%codeLeng * codeLeng
        for j in range(0,n):
            coded_tab[i][j] = codes[temp_infls[j//codeLeng] - 1][temp_row + j%codeLeng]
    return [infls_tab,coded_tab]

inflsFromFile=[]
with open("infls_pag.txt", "r") as f:
    for line in f:
        line=line.replace(',','').replace('[','').replace(']','').split()
        inflsFromFile.append([int(n) for n in line])

codesFromFile=[]
with open("codes_pag.txt", "r") as f:
    for line in f:
        line=line.replace(',','').replace('[','').replace(']','').split()
        codesFromFile.append([int(n) for n in line])

prolongablesFromFile=[]
with open("prolongable_pag.txt", "r") as f:
    for line in f:
        line=line.replace(',','').replace('[','').replace(']','').split()
        prolongablesFromFile.append([int(n) for n in line])


overlap=1
codeLeng = len(codesFromFile[0])
sqrtCodeLeng = int(sqrt(codeLeng))
codesNoOverlap = [[codesFromFile[i][j] for j in range(0,codeLeng) if j%sqrtCodeLeng < sqrtCodeLeng - overlap] for i in range(0,len(codesFromFile))]

char = 3
tel = 8
cid = 8     # tile edge length, centre distance
mhi = 1600
nd = 4800  # wall rows and cols/2
pag = [Pagoda(n) for n in range (-mhi,nd+mhi+1)]
wal = NumberWall(pag,-2,mhi,0,nd,char)

[infls_tab,coded_tab] = AutomaticNdTiling(inflsFromFile,codesNoOverlap,2,2+mhi+1,nd+1)
wal==coded_tab
True
inflsFromFile=[]
with open("infls_dragon.txt", "r") as f:
    for line in f:
        line=line.replace(',','').replace('[','').replace(']','').split()
        inflsFromFile.append([int(n) for n in line])

codesFromFile=[]
with open("codes_dragon.txt", "r") as f:
    for line in f:
        line=line.replace(',','').replace('[','').replace(']','').split()
        codesFromFile.append([int(n) for n in line])

prolongablesFromFile=[]
with open("prolongable_dragon.txt", "r") as f:
    for line in f:
        line=line.replace(',','').replace('[','').replace(']','').split()
        prolongablesFromFile.append([int(n) for n in line])

char = 3
tel = 12
cid = 8     # tile edge length, centre distance
mhi = 600#0
nd = 1500#0  # wall rows and cols/2
#rpf = [Dragon(n) for n in range (-mhi,nd+mhi+1)]
#wal = NumberWall(rpf,-2,mhi,0,nd,char)

overlap=5
codeLeng = len(codesFromFile[0])
sqrtCodeLeng = int(sqrt(codeLeng))
codesNoOverlap = [[codesFromFile[i][j] for j in range(0,codeLeng) if j%sqrtCodeLeng < sqrtCodeLeng - overlap] for i in range(0,len(codesFromFile))]

[infls_tab_dragon,coded_tab_dragon] = AutomaticNdTiling(inflsFromFile,codesNoOverlap,4,2+mhi+1,nd+1)