EscapeOfMassSelecta.ipynb

Kernel: SageMath 10.7

#define routines related to 1-dim and 2-dim substitutions, the computation of a number wall and 
#symmetries of the thue-morse number wall 
def SubstitutiveSequence(infls,prolong,leng):
    """
    This function returns the substitution sequence of length leng that is
    generated by repeated applications of the substitution infls on a 
    prolongable element prolong
    """
    infls_len = [len(infls[i]) for i in range(0,len(infls))]
    old_seq = [prolong]
    new_seq = infls[prolong]
    old_len = 1
    new_len = len(new_seq)
    next_len = new_len
    while next_len<leng:
        for n in range(old_len,new_len):
            next_len = next_len + infls_len[new_seq[n]]
        old_seq = new_seq
        new_seq = [0]*next_len
        new_seq[0:new_len] = old_seq
        cur_len = new_len
        for n in range(old_len,new_len):
            new_cur_len = cur_len + infls_len[new_seq[n]]
            new_seq[cur_len:new_cur_len] = infls[new_seq[n]]
            cur_len = new_cur_len
        old_len = new_len
        new_len = next_len
    return new_seq[0:leng]

def AutomaticGenerator(infls,codes,prolong,leng):
    """
    This function returns the automatic sequence of length leng that is
    generated by repeated applications of the substitution infls on a 
    prolongable element prolong followed by
    an application of the coding codes
    """
    seq = SubstitutiveSequence(infls,prolong,leng)
    len_seq = len(seq)
    coded_seq = [0]*len_seq
    coded_seq = [codes[seq[n]] for n in range(0,len_seq)]
    return coded_seq

def AutomaticGenerator2d(infls,codes,prolong,m,n):
    """
    This function returns the substitution array of size m times n leng that is
    generated by repeated applications of the substitution infls on a 
    prolongable element prolong. 
    The input codes is a list, and each symbol is coded into a symbol
    """
    inflsLeng = len(infls[0])
    temp_infls_row = [0] * n
    infls_tab = [temp_infls_row]*m
    
    infls_tab[0][0] = prolong
    for i in range(0,m):
        prev_row_index = i//inflsLeng
        row_index = i%inflsLeng
        for j in range(0,n):
            temp_infls_row[j] = infls[infls_tab[prev_row_index][j//inflsLeng]][row_index][j%inflsLeng]
        infls_tab[i] = copy(temp_infls_row)
    
    coded_tab = [[0] * n for i in range(0,m)]
    for i in range(0,m):
        temp_infls = infls_tab[i]
        for j in range(0,n):
            coded_tab[i][j] = codes[temp_infls[j]]
    return [infls_tab,coded_tab]

def SubAndCodGenerator2d(infls,codes,prolongable,m,n):
    """
    This function returns the substitution array of size m times n leng that is
    generated by repeated applications of the substitution infls on a 
    prolongable element prolong. 
    The input codes is a list of square arrays, and each symbol is coded into a square of symbols
    """
    codeLeng = len(codes[0])
    inflsLeng = len(infls[0])
    m_infls = ceil(m/codeLeng)
    n_infls = ceil(n/codeLeng)
    temp_infls_row = [0] * n_infls
    infls_tab = [temp_infls_row]*m_infls
    
    infls_tab[0][0] = prolongable
    for i in range(0,m_infls):
        prev_row_index = i//inflsLeng
        row_index = i%inflsLeng
        for j in range(0,n_infls):
            temp_infls_row[j] = infls[infls_tab[prev_row_index][j//inflsLeng]][row_index][j%inflsLeng]
        infls_tab[i] = copy(temp_infls_row)
    
    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
        for j in range(0,n):
            coded_tab[i][j] = codes[temp_infls[j//codeLeng]][temp_row][j%codeLeng]
    return [infls_tab,coded_tab]

def SimpleWall(seq,char=0):
    """
    For a sequence over a field with char elements, this function caluculates its number wall. 
    If char == 0 the number wall is calculated over the rationals.
    """
    # setting up the field
    if char!=0: 
        s=GF(char)
    else:
        s=QQ
    leng=len(seq)
    wal_width = leng + 4
    wal_height = floor(leng/2) + 2
    wal = [ [s(1)]* wal_width for m in range(0,wal_height)]
    for i in range(0,wal_height):  # assign each entry in matrix order
        for j in range(0,wal_width):
            if i == 0 or (j < i or j > wal_width - 1 - i):
                wal[i][j] = s(0)            # padding zeros on the first row and first and last columns 
            elif i == 1:
                wal[i][j] = s(1)            # the second row of the number wall consists ones 
            elif i == 2:
                wal[i][j] = s(seq[j-2])     # the third row of the number wall is the sequence itself 
            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] == 0 )
                p = 2
                while wal[i-p-1][j] == 0:
                    p += 1
                q = 1
                while wal[i-p][j-q] == 0 and p!=q:
                    q += 1
                d = q+1
                while wal[i-p][j+d-q] == 0 and p!=d-1:
                    d += 1                   # window defic.  d  and pane offsets  p, q

                if p < d-1:
                    wal[i][j] = s(0)         # 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] == 0 and wal[i-1,j] != 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]
                                                                             )
                                              )
    
    wal = [[wal[i][j] for j in range(2,wal_width - 2)] for i in range(2,wal_height)]
    
    return wal

def thm_rotation(s):
    """
    This function implements the rotation symmetry on a simple set of symbols that automtically generate the Thue-Morse 
    number-wall.
    """
    if s==0:
        return 2
    elif s==1:
        return 3
    elif s==2:
        return 4
    elif s==3:
        return 5
    elif s==4:
        return 6
    elif s==5:
        return 7
    elif s==6:
        return 0
    elif s==7:
        return 1
    elif s==10:
        return 11
    elif s==11:
        return 12 
    elif s==12:
        return 13
    elif s==13:
        return 10
    else:
        return s
    
def thm_reflection(s):
    """
    This function implements the reflection along the anti-diagonal symmetry on a simple set of symbols that automtically generate the Thue-Morse 
    number-wall.
    """
    if s==0:
        return 1
    elif s==1:
        return 0
    elif s==2:
        return 7
    elif s==3:
        return 6
    elif s==4:
        return 5
    elif s==5:
        return 4
    elif s==6:
        return 3
    elif s==7:
        return 2
    else:
        return s

def thm_inversion(s):
    """
    This function implements the inversion symmetry on a simple set of symbols that automtically generate the Thue-Morse 
    number-wall.
    """
    if s==0:
        return 1
    elif s==1:
        return 0
    elif s==2:
        return 3
    elif s==3:
        return 2
    elif s==4:
        return 5
    elif s==5:
        return 4
    elif s==6:
        return 7
    elif s==7:
        return 6
    if s==8:
        return 9
    elif s==9:
        return 8
    else:
        return s

def thm_matrix_rotation(M):
    """
    This function implements the rotation symmetry on matrices of size 2 over a simple set of symbols that automtically generate the Thue-Morse 
    number-wall.
    """
    temp = [[0,0],[0,0]]
    temp[0][0] = thm_inversion(thm_rotation(M[1][0]))
    temp[0][1] = thm_inversion(thm_rotation(M[0][0]))
    temp[1][0] = thm_inversion(thm_rotation(M[1][1]))
    temp[1][1] = thm_inversion(thm_rotation(M[0][1]))
    return temp
    
def thm_matrix_reflection(M):
    """
    This function implements the reflection symmetry on matrices of size 2 over a simple set of symbols that automtically generate the Thue-Morse 
    number-wall.
    """
    temp = [[0,0],[0,0]]
    temp[0][0] = thm_inversion(thm_reflection(M[1][1]))
    temp[0][1] = thm_inversion(thm_reflection(M[0][1]))
    temp[1][0] = thm_inversion(thm_reflection(M[1][0]))
    temp[1][1] = thm_inversion(thm_reflection(M[0][0]))
    return temp

def thm_matrix_inversion(M):
    """
    This function implements the inversion symmetry on matrices of size 2 over a simple set of symbols that automtically generate the Thue-Morse 
    number-wall.
    """
    temp = [[0,0],[0,0]]
    temp[0][0] = thm_inversion(M[0][0])
    temp[0][1] = thm_inversion(M[0][1])
    temp[1][0] = thm_inversion(M[1][0])
    temp[1][1] = thm_inversion(M[1][1])
    return temp

#compute the diagonally-aligned number wall of the Thue-Morse sequence
import time
start = time.time()
sequence_length=128
char=2
fieldSize = char
if char!=0:
    s = GF(fieldSize)
else:
    s=QQ

infls = [[0, 1], [1, 0]]
codes = [0,1]
prolongable = 0
seq = AutomaticGenerator(infls,codes,prolongable,sequence_length)
wal=SimpleWall(seq,fieldSize)

#align wal diagonally
diagonally_aligned_wal = [ [s(0)] * (sequence_length+2) for m in range(0,sequence_length+2)]

for m in srange(2,sequence_length+2):
    temp = m%2
    for n in srange(0,sequence_length+2):
        if temp==n%2 and n<m:
            diagonally_aligned_wal[m][n] = wal[(m-n)/2-1][(m+n)/2-1]

for m in range(0,sequence_length+2):
    diagonally_aligned_wal[m][m] = s(1)

#define a substitution and coding that generate the Thue-Morse diagonally-aligned number wall
infls = [
    [
        [0, 10], 
        [8, 1]], 
    [
        [1, 10], 
        [9, 0]], 
    [
        [9, 3], 
        [2, 11]], 
    [
        [8, 2], 
        [3, 11]],     
    [
        [5, 8], 
        [12, 4]], 
    [
        [4, 9], 
        [12, 5]], 
    [
        [13, 6], 
        [7, 9]], 
    [
        [13, 7], 
        [6, 8]], 
    [
        [3, 4], 
        [0, 7]], 
    [
        [2, 5], 
        [1, 6]], 
    [
        [14, 14], 
        [10, 14]], 
    [
        [11, 14], 
        [14, 14]], 
    [
        [14, 12], 
        [14, 14]], 
    [
        [14, 14], 
        [14, 13]],
    [
        [14,14],
        [14,14]]
]

codes_with_overlap = [
    [
        [0, 1, 0, 0, 0],
        [0, 0, 1, 0, 0],
        [0, 1, 0, 1, 0],
        [1, 0, 1, 0, 1],
        [0, 1, 0, 0, 0]],
    [
        [0, 1, 0, 0, 0],
        [1, 0, 1, 0, 0],
        [0, 0, 0, 1, 0],
        [0, 0, 0, 0, 1],
        [0, 0, 0, 1, 0]],
    [
        [0, 1, 0, 0, 0],
        [1, 0, 1, 0, 1],
        [0, 1, 0, 1, 0],
        [0, 0, 1, 0, 0],
        [0, 1, 0, 0, 0]],
    [
        [0, 0, 0, 1, 0],
        [0, 0, 0, 0, 1],
        [0, 0, 0, 1, 0],
        [1, 0, 1, 0, 0],
        [0, 1, 0, 0, 0]],
    [
        [0, 0, 0, 1, 0],
        [1, 0, 1, 0, 1],
        [0, 1, 0, 1, 0],
        [0, 0, 1, 0, 0],
        [0, 0, 0, 1, 0]],
    [
        [0, 1, 0, 0, 0],
        [1, 0, 0, 0, 0],
        [0, 1, 0, 0, 0],
        [0, 0, 1, 0, 1],
        [0, 0, 0, 1, 0]],
    [
        [0, 0, 0, 1, 0],
        [0, 0, 1, 0, 0],
        [0, 1, 0, 1, 0],
        [1, 0, 1, 0, 1],
        [0, 0, 0, 1, 0]],
    [
        [0, 0, 0, 1, 0],
        [0, 0, 1, 0, 1],
        [0, 1, 0, 0, 0],
        [1, 0, 0, 0, 0],
        [0, 1, 0, 0, 0]],
    [
        [0, 1, 0, 0, 0],
        [0, 0, 1, 0, 1],
        [0, 1, 0, 1, 0],
        [1, 0, 1, 0, 0],
        [0, 0, 0, 1, 0]],
    [
        [0, 0, 0, 1, 0],
        [1, 0, 1, 0, 0],
        [0, 1, 0, 1, 0],
        [0, 0, 1, 0, 1],
        [0, 1, 0, 0, 0]],
    [
        [0, 0, 0, 0, 0],
        [0, 0, 0, 0, 0],
        [0, 0, 0, 0, 0],
        [1, 0, 0, 0, 0],
        [0, 1, 0, 0, 0]],
    [
        [0, 1, 0, 0, 0],
        [1, 0, 0, 0, 0],
        [0, 0, 0, 0, 0],
        [0, 0, 0, 0, 0],
        [0, 0, 0, 0, 0]],
    [
        [0, 0, 0, 1, 0],
        [0, 0, 0, 0, 1],
        [0, 0, 0, 0, 0],
        [0, 0, 0, 0, 0],
        [0, 0, 0, 0, 0]],
    [
        [0, 0, 0, 0, 0],
        [0, 0, 0, 0, 0],
        [0, 0, 0, 0, 0],
        [0, 0, 0, 0, 1],
        [0, 0, 0, 1, 0]],
    [
        [0, 0, 0, 0, 0],
        [0, 0, 0, 0, 0],
        [0, 0, 0, 0, 0],
        [0, 0, 0, 0, 0],
        [0, 0, 0, 0, 0]]
]


codes = [[code[i][0:4] for i in range(0,4)] for code in codes_with_overlap]

prolongable = 0
l=sequence_length+1
width=l
leng=l+1
infls_tab,coded_tab=SubAndCodGenerator2d(infls,codes,prolongable,width,leng)

#check if diagonally_aligned_wall without the zeroth row equals coded_tab 
print(diagonally_aligned_wal[1:] == coded_tab)

True

#plot a portion of infls_tab
pl=32
for i in range(0,pl):
    temp = []
    for j in range(0,pl):
        if infls_tab[i][j]==0:
            temp += 'a0 '
        elif infls_tab[i][j]==2:
            temp += str('a1 ')
        elif infls_tab[i][j]==4:
            temp += str('a2 ')
        elif infls_tab[i][j]==6:
            temp += str('a3 ')
        elif infls_tab[i][j]==1:
            temp += str('b0 ')
        elif infls_tab[i][j]==3:
            temp += str('b1 ')
        elif infls_tab[i][j]==5:
            temp += str('b2 ')
        elif infls_tab[i][j]==7:
            temp += str('b3 ')
        elif infls_tab[i][j]==8:
            temp += str('da ')
        elif infls_tab[i][j]==9:
            temp += str('db ')
        elif infls_tab[i][j]==10:
            temp += str('c0 ')
        elif infls_tab[i][j]==11:
            temp += str('c1 ')
        elif infls_tab[i][j]==12:
            temp += str('c2 ')
        elif infls_tab[i][j]==13:
            temp += str('c3 ')
        elif infls_tab[i][j]==14:
            temp += str(' o ')
    print(*temp,sep = '')

a0 c0  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o 
da b0 c0  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o 
b1 a2 b0 c0  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o 
a0 b3 db a0 c0  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o 
da a1 b2 da b0 c0  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o 
b1 c1 c2 a2 db a0 c0  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o 
a0 c0 c3 b3 a1 b2 a0 c0  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o 
da b0 a3 da b0 a3 da b0 c0  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o 
b1 a2 db b1 a2 db b1 a2 b0 c0  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o 
a0 b3 a1 c1 c2 b2 a0 b3 db a0 c0  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o 
da a1 c1  o  o c2 b2 da a1 b2 a0 c0  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o 
b1 c1  o  o  o  o c2 a2 b0 a3 da b0 c0  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o 
a0 c0  o  o  o  o c3 b3 db b1 a2 db a0 c0  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o 
da b0 c0  o  o c3 a3 da a1 c1 c2 b2 da b0 c0  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o 
b1 a2 b0 c0 c3 a3 b1 a2 b0 c0 c3 a3 b1 a2 b0 c0  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o 
a0 b3 db a0 b3 db a0 b3 db a0 b3 db a0 b3 db a0 c0  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o 
da a1 b2 da a1 b2 da a1 b2 da a1 b2 da a1 b2 da b0 c0  o  o  o  o  o  o  o  o  o  o  o  o  o  o 
b1 c1 c2 a2 b0 a3 b1 c1 c2 a2 b0 a3 b1 c1 c2 a2 db a0 c0  o  o  o  o  o  o  o  o  o  o  o  o  o 
a0 c0 c3 b3 db b1 c1  o  o c2 a2 db a0 c0 c3 b3 a1 b2 a0 c0  o  o  o  o  o  o  o  o  o  o  o  o 
da b0 a3 da a1 c1  o  o  o  o c2 b2 da b0 a3 da b0 a3 da b0 c0  o  o  o  o  o  o  o  o  o  o  o 
b1 a2 db b1 c1  o  o  o  o  o  o c2 a2 db b1 a2 db b1 a2 db a0 c0  o  o  o  o  o  o  o  o  o  o 
a0 b3 a1 c1  o  o  o  o  o  o  o  o c2 b2 a0 b3 a1 c1 c2 b2 da b0 c0  o  o  o  o  o  o  o  o  o 
da a1 c1  o  o  o  o  o  o  o  o  o  o c2 b2 da b0 c0 c3 a3 b1 a2 b0 c0  o  o  o  o  o  o  o  o 
b1 c1  o  o  o  o  o  o  o  o  o  o  o  o c2 a2 db a0 b3 db a0 b3 db a0 c0  o  o  o  o  o  o  o 
a0 c0  o  o  o  o  o  o  o  o  o  o  o  o c3 b3 a1 b2 da a1 b2 da a1 b2 a0 c0  o  o  o  o  o  o 
da b0 c0  o  o  o  o  o  o  o  o  o  o c3 a3 da b0 a3 b1 c1 c2 a2 b0 a3 da b0 c0  o  o  o  o  o 
b1 a2 b0 c0  o  o  o  o  o  o  o  o c3 a3 b1 a2 db b1 c1  o  o c2 a2 db b1 a2 b0 c0  o  o  o  o 
a0 b3 db a0 c0  o  o  o  o  o  o c3 b3 db a0 b3 a1 c1  o  o  o  o c2 b2 a0 b3 db a0 c0  o  o  o 
da a1 b2 da b0 c0  o  o  o  o c3 a3 da a1 b2 da b0 c0  o  o  o  o c3 a3 da a1 b2 da b0 c0  o  o 
b1 c1 c2 a2 db a0 c0  o  o c3 b3 db b1 c1 c2 a2 db a0 c0  o  o c3 b3 db b1 c1 c2 a2 db a0 c0  o 
a0 c0 c3 b3 a1 b2 a0 c0 c3 b3 a1 b2 a0 c0 c3 b3 a1 b2 a0 c0 c3 b3 a1 b2 a0 c0 c3 b3 a1 b2 a0 c0 
da b0 a3 da b0 a3 da b0 a3 da b0 a3 da b0 a3 da b0 a3 da b0 a3 da b0 a3 da b0 a3 da b0 a3 da b0 

#plot a portion of coded_tab
print(*([1]+[0]*(pl-1)))
for i in range(0,pl-1):
    print(*coded_tab[i][0:pl])

1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
1 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 1 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 1 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
1 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 1 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 1 0 1 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 1 0 1 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
1 0 1 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 1 0 0 0 0 0 1 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 1 0 0 0 1 0 1 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 1 0 1 0 1 0 0 0 1 0 1 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
1 0 1 0 1 0 0 0 0 0 1 0 1 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 1 0 1 0 1 0 1 0 0 0 0 0 1 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0
0 1 0 1 0 1 0 1 0 1 0 0 0 1 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0
1 0 1 0 0 0 1 0 0 0 1 0 1 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0
0 0 0 1 0 1 0 0 0 0 0 1 0 0 0 1 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 0
0 0 0 0 1 0 0 0 0 0 0 0 1 0 1 0 1 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0
0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 1 0 1 0 1 0 1 0 1 0 0 0 0 0 0 0 0
1 0 1 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 1 0 1 0 1 0 0 0 0 0 0 0
0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 1 0 0 0 1 0 0 0 0 0 0
0 0 1 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 1 0 1 0 0 0 0 0 1 0 0 0 0 0
0 1 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 1 0 1 0 0 0 1 0 1 0 0 0 0
1 0 1 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 1 0 1 0 1 0 1 0 0 0
0 1 0 0 0 1 0 0 0 0 0 1 0 1 0 0 0 1 0 0 0 0 0 1 0 1 0 0 0 1 0 0
0 0 1 0 1 0 1 0 0 0 1 0 0 0 1 0 1 0 1 0 0 0 1 0 0 0 1 0 1 0 1 0
0 1 0 1 0 0 0 1 0 1 0 1 0 1 0 1 0 0 0 1 0 1 0 1 0 1 0 1 0 0 0 1

from PIL import Image
import numpy as np
import time
start = time.time()

dawal = np.zeros([width+1, leng], dtype=np.uint8)
dawal[0][0]=0
for n in srange(1,leng):
    dawal[0][n]=255
for m in srange(1,width+1):
    for n in srange(0,leng):
        dawal[m][n]=(1-coded_tab[m-1][n]/(char-1))*255

# Convert Numpy array to palette image
pi = Image.fromarray(dawal)

# Display and save
pi.save('thm_dawal_'+str(l-1)+'.png')

end = time.time()
print(end-start)

0.6398732662200928

# convert flat to wall
start = time.time()

l1=ceil(l/2)
wal = np.ones([l1,l+1], dtype=np.uint8)
for m in srange(0,l1):
    for n in srange(0,l+1):
        if m<=n & n<=l-m:
            #print(m-l0,n,m-l0+n+1,n-m+l0,flat[m-l0+n+1][n-m+l0])
            wal[m][n] = dawal[m+n][n-m]
        else:
            wal[m][n] = 255

# Convert Numpy array to palette image
pi = Image.fromarray(wal)

# Display and save
pi.save('thm_wal_'+str(l-1)+'.png')

end = time.time()
print(end-start)

0.026792287826538086

#make a list of all the size two squares contained in infls_tab
tetrads=[]
size = floor(l/4)
for i in range(0,size-1):
    for j in range(0,size-1):
        tetrad = [[infls_tab[i][j],infls_tab[i][j+1]],[infls_tab[i+1][j],infls_tab[i+1][j+1]]]
        #tetrad = [[[coded_tab[i][j]],[coded_tab[i][j+1]]],[[coded_tab[i+1][j]],[coded_tab[i+1][j+1]]]]
        
        if tetrad not in tetrads:
            tetrads.append(tetrad)
print(len(tetrads))

64

#verify that all the size two squares are consistent with respect to an overlap of size 1
for tetrad in tetrads:
    if not all([[codes_with_overlap[tetrad[0][0]][j][4] for j in range(0,5)]==[codes_with_overlap[tetrad[0][1]][j][0] for j in range(0,5)],
    [codes_with_overlap[tetrad[1][0]][j][4] for j in range(0,5)]==[codes_with_overlap[tetrad[1][1]][j][0] for j in range(0,5)],
    codes_with_overlap[tetrad[0][0]][4]==codes_with_overlap[tetrad[1][0]][0],
    codes_with_overlap[tetrad[0][1]][4]==codes_with_overlap[tetrad[1][1]][0]]):
        print(false)
print(true)        

True

#creates two lists of all the size two squares contained in infls_tab up to the actions of H and G, respectively

tetradsH=[]
tetradsG=[]

size = floor(l/4)
for k in range(0,64):
    tetrad = [[tetrads[k][0][0],tetrads[k][0][1]],[tetrads[k][1][0],tetrads[k][1][1]]]
    
    temp1 = thm_matrix_rotation(tetrad)
    temp2 = thm_matrix_rotation(temp1)
    temp3 = thm_matrix_rotation(temp2)
    
    tetradHOrbit = [tetrad,
                   temp1,
                   temp2,
                   temp3,
                   thm_matrix_reflection(tetrad),
                   thm_matrix_reflection(temp3),
                   thm_matrix_reflection(temp2),
                   thm_matrix_reflection(temp1)]
    
    #print(tetradHOrbit)
    
    if all([tet not in tetradsH for tet in tetradHOrbit]):
        tetradsH.append(tetrad)
    
    tetradHCoset = [thm_matrix_inversion(g) for g in tetradHOrbit]
    tetradGOrbit = tetradHOrbit + tetradHCoset

    if all([tet not in tetradsG for tet in tetradGOrbit]):
        tetradsG.append(tetrad)
print(len(tetradsH))
print(len(tetradsG))

16
11

#make a list of all the size two squares contained in infls_tab
pairs=[]
size = floor(l/4)
for i in range(0,size-1):
    for j in range(0,size-1):
        pair = [infls_tab[i][j],infls_tab[i][j+1]]
        
        if pair not in pairs:
            pairs.append(pair)
print(len(pairs))

31

#creates two lists of all the 1*2 patterns contained in infls_tab up to the actions of H2={1,\eta\rho,\rho\eta,\rho^2} and 
#G2={1,\eta\rho,\rho\eta,\rho^2}*{\iota}, the index 2 subgroups of H and G that preserve the rows, respectively

pairsH=[]
pairsG=[]

size = floor(l/4)
for k in range(0,31):
    pair = [pairs[k][0],pairs[k][1]]
    
    temp1 = [thm_rotation(thm_rotation(pair[1])),thm_rotation(thm_rotation(pair[0]))]
    temp2 = [thm_reflection(thm_rotation(pair[1])),thm_reflection(thm_rotation(pair[0]))]
    temp3 = [thm_rotation(thm_reflection(pair[0])),thm_rotation(thm_reflection(pair[1]))]
    
    pairHOrbit = [pair,
                   temp1,
                   temp2,
                   temp3]
        
    if all([pair not in pairsH for pair in pairHOrbit]):
        pairsH.append(pair)
    
    pairHCoset = [[thm_inversion(g[0]),thm_inversion(g[1])] for g in pairHOrbit]
    pairGOrbit = pairHOrbit + pairHCoset

    if all([pair not in pairsG for pair in pairGOrbit]):
        pairsG.append(pair)
print(len(pairsH))
print(len(pairsG))

13
9

for i in range(0,9):
    print(pairsG[i])

[0, 10]
[10, 14]
[14, 14]
[8, 1]
[3, 4]
[4, 1]
[3, 11]
[11, 12]
[11, 14]

for k in range(0,16):
    for i in range(0,2):
        temp = []
        for j in range(0,2):
            if tetradsH[k][i][j]==0:
                temp += 'a0 '
            elif tetradsH[k][i][j]==2:
                temp += str('a1 ')
            elif tetradsH[k][i][j]==4:
                temp += str('a2 ')
            elif tetradsH[k][i][j]==6:
                temp += str('a3 ')
            elif tetradsH[k][i][j]==1:
                temp += str('b0 ')
            elif tetradsH[k][i][j]==3:
                temp += str('b1 ')
            elif tetradsH[k][i][j]==5:
                temp += str('b2 ')
            elif tetradsH[k][i][j]==7:
                temp += str('b3 ')
            elif tetradsH[k][i][j]==8:
                temp += str('da ')
            elif tetradsH[k][i][j]==9:
                temp += str('db ')
            elif tetradsH[k][i][j]==10:
                temp += str('c0 ')
            elif tetradsH[k][i][j]==11:
                temp += str('c1 ')
            elif tetradsH[k][i][j]==12:
                temp += str('c2 ')
            elif tetradsH[k][i][j]==13:
                temp += str('c3 ')
            elif tetradsH[k][i][j]==14:
                temp += str(' o ')
        print(*temp,sep = '')
    print('')

a0 c0 
da b0 

c0  o 
b0 c0 

 o  o 
c0  o 

 o  o 
 o  o 

da b0 
b1 a2 

b0 c0 
a2 b0 

b1 a2 
a0 b3 

c0  o 
a0 c0 

b3 db 
a1 b2 

db a0 
b2 da 

a1 b2 
c1 c2 

b1 c1 
a0 c0 

c1 c2 
c0 c3 

a0 c0 
b2 a0 

a1 b2 
b0 a3 

c1 c2 
 o  o 

for k in range(0,11):
    for i in range(0,2):
        temp = []
        for j in range(0,2):
            if tetradsG[k][i][j]==0:
                temp += 'a0 '
            elif tetradsG[k][i][j]==2:
                temp += str('a1 ')
            elif tetradsG[k][i][j]==4:
                temp += str('a2 ')
            elif tetradsG[k][i][j]==6:
                temp += str('a3 ')
            elif tetradsG[k][i][j]==1:
                temp += str('b0 ')
            elif tetradsG[k][i][j]==3:
                temp += str('b1 ')
            elif tetradsG[k][i][j]==5:
                temp += str('b2 ')
            elif tetradsG[k][i][j]==7:
                temp += str('b3 ')
            elif tetradsG[k][i][j]==8:
                temp += str('da ')
            elif tetradsG[k][i][j]==9:
                temp += str('db ')
            elif tetradsG[k][i][j]==10:
                temp += str('c0 ')
            elif tetradsG[k][i][j]==11:
                temp += str('c1 ')
            elif tetradsG[k][i][j]==12:
                temp += str('c2 ')
            elif tetradsG[k][i][j]==13:
                temp += str('c3 ')
            elif tetradsG[k][i][j]==14:
                temp += str(' o ')
        print(*temp,sep = '')
    print('')

a0 c0 
da b0 

c0  o 
b0 c0 

 o  o 
c0  o 

 o  o 
 o  o 

da b0 
b1 a2 

b0 c0 
a2 b0 

b1 a2 
a0 b3 

db a0 
b2 da 

a1 b2 
c1 c2 

c1 c2 
c0 c3 

c1 c2 
 o  o