对称之美：李代数 $E_8$ 的根系图

#coding=utf-8


#----------------------------------
#This script draws the 2d projection of E8 root system
#into its Coxeter plane. For the math behind this script,
#see Humphreys‘s book "Reflection groups and Coxeter groups"
#section 17, chapter 3
#----------------------------------


import numpy as np
import matplotlib.pyplot as plt
from itertools import product, combinations

roots = []
#roots of the form (+-1,+-1,0,0,0,0,0,0)
for i,j in combinations(range(8),2):
    for x,y in product([-1,1],repeat=2):
        v = np.zeros(8)
        v[i] = x
        v[j] = y
        roots.append(v)

#roots of the form 1/2 * (+-1, +-1, ..., +-1)
for v in product([-1,1],repeat=8):
    if sum(v) % 4 == 0:
        roots.append(np.array(v)*0.5)

#get all edges
edges = []
for i in range(240):
    for j in range(i+1,240):
        if abs(np.sum((roots[i]-roots[j])**2)-2)<1e-8:
            edges.append([i,j])

print("There are %03d vertices and %04d edges"%(len(roots),len(edges)))
# choose a complete set of simple roots. They are copied from

#    https://en.wikipedia.org/wiki/E8_(mathematics)#Cartan_matrix

# These simple roots are given by the rows of the matrix delta.
delta = np.eye(8)
for i in range(5):
    delta[i,i+1] = -1
delta[5,6] = 1
delta[6] = -0.5
delta[7,5] = 1
delta[7,6] = -1
delta[7,7] = 0

print("The simple roots:")
print delta
"""
To costruct a bais of the coxeter plane, we fist split the simple roots
into two disjoint sets I and J. The roots in I are pairwise orthorgonal (and so is J).
By our choice they are I={0,2,4,6} and J={1,3,5,7}, which means for any pair
k,l in I (or in J), delta[k] and delta[l] are orthogonal.

The cartan matrix of E8 has a unique maximal positive eigenvalue with eigenvector c. Then
u = sum (c[i] * delta[i]) for i in I
v = sum (c[j] * delta[j]) for j in J
form a basis of the Coexter plane. Gram-Schimdt them, and project the root system to this plane. i.e. calculate the inner product <u,x>,<v,x> for all x in the root system, this is the 2d point we want. 

The Dynkin graph:

     1---2---3---4---5---6---7
                     |
                     8
"""
I = [0,2,4,6]
J = [1,3,5,7]

#cartan matrix of the E8 root system
cartan = np.dot(delta,delta.transpose())
print "The cartan matrix:"
print cartan

eigenvals,eigenvecs= np.linalg.eigh(cartan)

#the eigenvalues returned by eigh() are in ascending order,
#and the eigenvectors are listed by columns.
#c is the eigenvector for the max-eigenvalue of the catan matrix. 
c = eigenvecs[:,7]
print "The eigenvector for the maximal eigenvalue of the Cartan matrix:"
print c

#u, v form a basis of the coxeter plane.
u = np.sum([c[i]*delta[i] for i in I],axis=0)
v = np.sum([c[j]*delta[j] for j in J],axis=0)
print "A basis of the Coxeter plane:"
print u, v

# Gram-Schimdt. Then u,v are orthogonal unit vectors
u /= np.linalg.norm(u)
v = v - np.dot(u,v)*u
v /= np.linalg.norm(v)
print "u,v after Gram-Schimdt and normalized to unit vectors:"
print u,v

fig = plt.figure(figsize=(8,8),dpi=100)
ax = fig.add_axes([0,0,1,1],aspect=1)
ax.axis([-1.3,1.3,-1.3,1.3])
ax.axis(‘off‘)

def edge_color_map(x):
    if x < 0.3:
       return (0.75, 0.5, 0.25)  
    elif x < 0.4:
        return (0, 0, 1)
    elif x < 0.5:
        return (1, 0, 1)
    elif x < 0.6:
        return (0, 1, 1)
    elif x < 0.75:
        return (1, 1, 0)
    elif x < 0.8:
        return (0, 1, 0)
    elif x < 1.0:
        return (0.5, 0.25, 0.75)
    else:
        return (1, 0, 0)
    
roots_2d = np.zeros((240,2))
modulus = np.zeros(240)
colors = np.zeros((240,3))
#project all roots to 2d by u,v
#and assign a color to each 2d point
for i in range(240):
    x = np.dot(roots[i],u)
    y = np.dot(roots[i],v)
    roots_2d[i] = [x,y]
    modulus[i] = np.linalg.norm([x,y])
    colors[i] = edge_color_map(modulus[i])

LineWidth = max(modulus)/8.0

#draw edges:
for e in edges:
    x = roots_2d[e[0]]
    y = roots_2d[e[1]]
    a, b = x
    c, d = y
    if modulus[e[0]] > modulus[e[1]]:
        color = colors[e[0]]
    else:
        color = colors[e[1]]
    ax.plot([a,c], [b,d],color=color,lw=LineWidth)

#draw vertices:
for i in range(240):
    x,y = roots_2d[i]
    ax.plot([x],[y],‘o‘,color=colors[i],ms=4)

plt.savefig("E8_root_system.png")