Author: Erez Nesharim
Views : 284
Description: Calculates the number wall of a sequence and plot
Compute Environment: Ubuntu 18.04 (Deprecated)
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 = [ * (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=*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
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+=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])
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
tetlen += 1
s = tetlen
quarter.append(False)
#if s==214:
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)
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)]
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=
merp =  + merp
#for i in range(0,7):
#   print(stab[i])
stab = [ [ merp[s] for s in t ] for t in stab ]
#for i in range(0,7):
#   print(stab[i])
print("final pass: states stabilised")
end = time.time()
print(end-start)

# 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)
end = time.time()
print(end-start)
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.43182682991 135.227772951 ('deficiencies: ', [1, 4609983, 6149634, 2305916, 0, 0, 0, 0, 0, 0]) 216.238996983 ('SquareTiling() initial pass: parameters = ', -38, 1602, -4800, 4800, 8, 8, 7.152557373046875e-06) ('first pass: segment tiled, state count', 201, ' ', 37.637633085250854) ('second pass: inflation consistent', 42.04904818534851) ('third pass: vertices complete, tetrad count = ', 1241, ' ', 60.32440400123596) final pass: states stabilised 72.2449820042 ('Pagoda construction: rows, cols = ', 1600, 9601) 288.48642087
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)
end = time.time()
print(end-start)
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):
"generates a rectangular portion of an automatic tiling by substitution and coding: Input 1. A k-substitution given by a list of arrays of size k^2. 2. An l-coding given by a list of arrays of size l^2.  3. A prolongable element s. 4. Number of rows m. 5. Number of columns n. Output is: A two dimensional array of size (m,n) which is a rectangular pattern with m rows and n column of the automatic tiling that is generated by the k-substitution applied to the prolongable element s, followed by the l-coding."
codeLeng = int(sqrt(len(codes)))
inflsLeng = int(sqrt(len(infls)))
m_infls = m//codeLeng + 1
n_infls = n//codeLeng + 1
infls_tab = [ * n_infls for t in range(0,m_infls)]
infls_tab = 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 = [ * 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)
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)

# Be sure to use a prolongable element. Currently inflsFromFile was generated by a version of SquareTiling with line 224 replaced with line 223, to make it equal to previous observations about the Pagoda wall. If inflsFromFile is generated using SquareTiling as is, call AutomaticNdTiling with the symbol 3 instead of 2
[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)
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)