SharedNumber wall.sagewsOpen in CoCalc
Author: Erez Nesharim
Views : 45
Description: 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.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) 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)