drcubic.py

import networkx as nx
import sympy

# intersection numbers 'iota'
K4 = [[3], [1]]
K33 = [[3, 2], [1, 3]]
Petersen = [[3, 2], [1, 1]]
Cube = [[3, 2, 1], [1, 2, 3]]
Heawood = [[3, 2, 2], [1, 1, 3]]
Pappus = [[3, 2, 2, 1], [1, 1, 2, 3]]
Coxeter = [[3, 2, 2, 1], [1, 1, 1, 2]]
Tutte = [[3, 2, 2, 2], [1, 1, 1, 3]]
Dodecahedron = [[3, 2, 1, 1, 1], [1, 1, 1, 2, 3]]
Desargues = [[3, 2, 2, 1, 1], [1, 1, 2, 2, 3]]
Biggs = [[3, 2, 2, 2, 1, 1, 1], [1, 1, 1, 1, 1, 1, 3]]
Foster = [[3, 2, 2, 2, 2, 1, 1, 1], [1, 1, 1, 1, 2, 2, 2, 3]]

smith1971 = [K4, K33, Petersen, Cube, Heawood, Pappus, Coxeter, Tutte, Dodecahedron, Desargues, Biggs, Foster]

"""
From Wolfram:
    Given a distance-regular graph G with integers  b_i,c_i,i=0,...,d such that for any two vertices
    x,y in G at distance i=d(x,y), there are exactly c_i neighbors of y in G_(i-1)(x) and b_i neighbors of y in G_(i+1)(x), the sequence
        iota(gamma)={b_0,b_1,...,b_(d-1);c_1,...,c_d}
    is called the intersection array of G.
See also https://en.wikipedia.org/wiki/Distance-regular_graph.
"""
def intersection_array_to_R(iota):
    diam = len(iota[0])
    b = lambda i: iota[0][i] if i in range(diam) else 0
    c = lambda i: iota[1][i-1] if i-1 in range(diam) else 0
    a = lambda i: b(0) - b(i) - c(i)

    B = sympy.zeros(diam+1)
    B[0,:2] = [[a(0), b(0)]]
    for i in range(1, diam):
        B[i,i-1:i+2] = [[c(i), a(i), b(i)]]
    B[-1,-2:] = [[c(diam), a(diam)]]

    R1 = sympy.Matrix([[B[i,j]/sum(B[i,:])  # normalized
        for j in range(diam+1)] for i in range(diam+1)])
    R = [R1**i for i in range(diam+1)]
    for i in range(diam+1):  # gaussian elimination of the "top layer"
        R[i] /= R[i][0,i]
        for j in range(i+1, diam+1):
            R[j] -= R[j][0,i]*R[i]

    return R

def intersection_array_to_L(iota):
    R = intersection_array_to_R(iota)
    L = [A.T for A in R]
    return L

def demo(which=range(12)):
    for i in which:
        print('intersection array', smith1971[i], sep='\n')
        print('left representation', intersection_array_to_L(smith1971[i]), sep='\n')
        print('irreducible characters', repr(characters(intersection_array_to_R(smith1971[i]))), sep='\n')
        print()

def characters(R):
    I = tuple(range(len(R)))
    x = sympy.var(['x'+str(i) for i in I])
    eqns = [x[i]*x[j]-sum(R[i][j,k]*x[k] for k in I)
        for i in I for j in I] + [x[0]-1]
    chars = sympy.solve(eqns, x)
    return sympy.Matrix(chars)

if __name__ == '__main__':
    demo()