Contact
CoCalc Logo Icon
StoreFeaturesDocsShareSupport News Sign UpSign In
| Download

Calculates the number wall of a sequence and plot

Views: 377
Image: ubuntu2204
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): "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[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) # 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[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)