Author: Erez Nesharim
Views : 137
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):
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)

[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)