"""

A collection of graph recipes

"""

# Some strongly regular graphs

# I need the following because I need to control the vertex set for the Chang graphs	
def LKn(n):
	vxs = [(i,j) for i in [0..n-1] for j in [i+1..n-1]]
	def incident(p1,p2):
		return p1[0]==p2[0] or p1[0]==p2[1] or p1[1]==p2[0] or p1[1]==p2[1]
	return Graph([vxs, lambda p1,p2: incident(p1,p2)])

def Paley(q):
	K.<a> = GF(q)
	return Graph([K, lambda i,j: i != j and (i-j).is_square()])

# Paley49 = Paley(49)
K.<a> = GF(49)
pows = [a^(4*i) for i in [1..12]] + [a^(4*i+1) for i in [1..12]]
# Peisert graph: arc transitive and self complementary, not Paley
Peis49 = Graph([K, lambda i,j: i != j and i-j in pows])

def cubelike(vecls): # input a list of binary vectors of length d
	d = len(vecls[0])
	VS = VectorSpace(GF(2),d)
	return Graph([[0..(2^d-1)], lambda i,j: VS[i]-VS[j] in vecls])

def LatinSqGrf(n, fun):
	vxs = [tuple([i,j]) for i in [0..n-1] for j in [0..n-1]]
	return Graph([vxs, lambda t1, t2: fun(t1,t2,n)], loops=False)

def cyclicLSfun(tpl1,tpl2,n):
	return (tpl1[0] == tpl2[0]) or (tpl1[1] == tpl2[1])\
	 or ((tpl1[0]+tpl1[1]-tpl2[0]-tpl2[1]).mod(n) == 0)

Shrikande = LatinSqGrf(4, cyclicLSfun).complement()
LS25a = LatinSqGrf(5, cyclicLSfun)  # isomorphic to Paley graph on GF(25)
LS36a = LatinSqGrf(6, cyclicLSfun) 

# switch on the subset sub of the vertices of G
def switch(G,sub):
	H = G.copy()
	rest = Set(H.vertices()).difference(Set(sub))
	for u in sub:
		for v in rest:
			if H.has_edge(u,v): H.delete_edge(u,v)
			else: H.add_edge(u,v)
	return H

def Schlaefli():
	L8 = LK8.copy()
	S0 = switch(L8, L8.neighbors(L8.vertices()[0]))
	S0.delete_vertex((0,1))
	return S0

# Chang graphs on 28 vertices
# switch on C8, 2C4, C3+C5
#### This is not working properly ####

LK8 = (graphs.CompleteGraph(8)).line_graph()

C8 = LK8.subgraph(vertices = [LK8.vertices()[i] for i in [0, 7, 13, 18, 22, 25, 27, 6]])
ChangC8 = switch(LK8,C8)  

C44 = LK8.subgraph(vertices = [LK8.vertices()[i] for i in [0, 7, 13, 2, 22, 25, 27, 24]])
ChangC44 = switch(LK8,C44)

C35 = LK8.subgraph(vertices = [LK8.vertices()[i] for i in [0, 7, 1, 18, 22, 25, 27, 21]])
ChangC35 = switch(LK8,C35)

def HoffmanSingleton():
	hsvxs = [(x,y,i) for x in GF(5) for y in GF(5) for i in [0,1]]
	def hosi(s,t):
		if s[2] == t[2] == 0:
			return (s[0]==t[0]) and (s[1]-t[1] in [1,4])
		elif s[2] == t[2] == 1:
			return (s[0]==t[0]) and (s[1]-t[1] in [2,3])
		elif s[2] == 0:
			return s[1] == t[0]*s[0]+t[1]
		else:
			return t[1] == s[0]*t[0]+s[1]
	return Graph([hsvxs, hosi])

# Note: this calls HoffmanSingleton()
# vertices are cocliques of size 15 in HoffmanSingleton - perhaps they should be compacted?
def HigmanSims():
    hs = HoffmanSingleton()
    clqs = (hs.complement()).cliques_maximum()
    alqs = [Set(cl) for cl in clqs]
    return Graph([alqs, lambda s, t: len(s.intersection(t)) in [0,8]])

hisi = HigmanSims()
vxs_m22 = [v for v in hisi.vertices() if hisi.distance(v,hisi.vertices()[0])==2]	
vxs_gw = [v for v in hisi.vertices() if hisi.distance(v,hisi.vertices()[0])==2 and\
  hisi.distance(v,hisi.vertices()[2])==2]

# the M_22 graph on 77 vertices
M22 = hisi.subgraph( vertices = vxs_m22)

# the Gewirtz graph on 56 vertices
Gewirtz = hisi.subgraph( vertices=vxs_gw)

# another construction of M22
xs
#########################################################################################
## Now some distance-regular graphs:

# J(v,k,r), J(v,k), Kneser(v,k)
# the first version is much slower: too many calls to Set()

# def gen_johnson2(v,k,r):
# 	combs = Combinations([0..v-1], k)
# 	return Graph([[0..binomial(v,k)-1], lambda i,j:\
# 	  len(set(combs[i]).intersection(set(combs[j])))==r])
	
def gen_johnson(v,k,r):
    combs = [Set(it) for it in Combinations([0..v-1], k)]
    return Graph([[0..binomial(v,k)-1], lambda i,j: len((combs[i]).intersection(combs[j]))==r])

def Johnson(v,k):
	return gen_johnson(v,k,k-1)

def Kneser(v,k):
	return gen_johnson(v,k,0)

def Coxeter():
	kk = Kneser(7,3)
	clqs = (kk.complement()).cliques()
	dlqs = filter(lambda cl: len(cl)==7,clqs)
	cl0 = dlqs[0]
	sub_vs = [vx for vx in kk.vertices() if vx not in cl0]
	Cox = kk.subgraph(sub_vs)
	Cox.relabel()
	return(Cox)

#######################################################################################
nbrs=[[1,2,3,11], [0,2,4,5], [0,1,3,9], [0,2,6,7], [1,5,7,8], [1,4,6,8],\
  [3,5,7,10],[3,4,6,10],[4,5,9,11],[2,8,10,11],[6,7,9,11],[0,8,9,10]]	
d = dict(zip([0..11],nbrs))

# probably the smallest graph that is walk regular but not vertex transitive
G = Graph(dict(zip([0..11],nbrs))) 
H = (graphs.CubeGraph(3)).line_graph()
# G and H are cospectral:
# H.characteristic_polynomial().factor()
#   (x - 4) * (x - 2)^3 * x^3 * (x + 2)^5
# G.characteristic_polynomial().factor()
#   (x - 4) * (x - 2)^3 * x^3 * (x + 2)^5
#
	
# construct a Cayley graph for Sym(4)
# P  = PermutationGroup([[(1,2)], [(2,3)], [(3,4)]])
# gp = P.cayley_graph()
# gq = gp.to_undirected()
# gq.spectrum()

ds1 = [(1,0,0), (4,0,0),\
 (0,1,1), (0,4,1), (2,0,1), (4,2,1), (4,3,1),\
 (0,1,2), (0,4,2), (1,2,2), (1,3,2), (2,0,2), (3,1,2), (3,2,2),\
 (3,3,2), (3,4,2), (4,0,2),\
 (0,1,3), (0,4,3), (1,0,3), (2,2,3), (2,3,3)]

H1 = PermutationGroup([[(1,2,3,4,5)],[(6,7,8,9,10)],[(2,3,5,4)]])
y,z,x = H1.gens()
C = [ x^it[0]*y^it[1]*z^it[2] for it in ds1]
HS = H1.cayley_graph(generators=C).to_undirected() # HiSi again

H2 = PermutationGroup([[(1,2,3,4,5)],[(6,7,8,9,10)],[(2,3,5,4),(7,10),(9,8)]])
b,c,a = H2.gens()
D = [ a^it[0]*b^it[1]*c^it[2] for it in ds1]

H3 = PermutationGroup([[(1,2,3,4,5)],[(6,7,8,9,10)],[(2,4,5,3),(7,9,10,8)]])
y,z,x = H3.gens()

ds3 = [(1,2,0),(2,1,0),(2,2,0),(3,3,0),(3,4,0),(4,3,0),\
 (0,4,1),(1,1,1),(2,2,1),(2,3,1),(2,4,1),(3,2,1),(3,3,1),(4,0,1),(4,2,1),(4,4,1),\
 (0,1,2),(0,2,2),(0,3,2),(1,0,2),(1,3,2),(2,0,2),(2,2,2),(3,0,2),(3,1,2),(3,3,2),\
 (0,3,3),(1,1,3),(1,4,3),(2,2,3),(3,0,3),(3,3,3),(3,4,3),(4,1,3),(4,3,3),(4,4,3)]
E = [ x^it[0]*y^it[1]*z^it[2] for it in ds3]
HW = H3.cayley_graph(generators=E)# Hall-Wales

##
# generalized quadrangle
##
def beta(u,v):
    return u[0]*v[2]+u[1]*v[3] -u[2]*v[0] -u[3]*v[1]
V = VectorSpace( GF(3), 4)
Points = [u[1] for u in V.subspaces(1)]
Lines = [sub for sub in V.subspaces(2)\
	  if beta(sub.matrix()[0],sub.matrix()[1])==0]

W3 = Graph([range(len(Points)),\
   lambda i,j: beta(Points[i],Points[j])==0], loops=False)
W3d = Graph( [range(len(Lines)), lambda i,j: i != j\
    and det( Lines[i].matrix().stack(Lines[j].matrix())) == 0])

##
# McLaughlin graph
##
def McL():
    C = codes.ExtendedBinaryGolayCode()
    D = C.punctured([0])
    words = [ Set(it.support()) for it in D if hamming_weight(it)==7]
    MG = Graph( [words, lambda a,b: len(a.intersection(b))==1])
    MM = MG.copy()
    MM.add_vertices([0..22])
    edges = [ (i,a) for i in [0..22] for a in words if i not in a]
    MM.add_edges( edges)
    McL = switch( MM, MM[0])
    McL.delete_vertex(0)
    return McL
    
#############
## cfi graphs
#############

def cfi( grf):
    VA = [(v,w,i) for v in grf.vertices() for w in grf[v] for i in [0,1]]
    VM = [(v,evns) for v in grf.vertices() for evns in Set(grf[v]).subsets()\
      if len(evns) % 2 == 0]
    VAVA = [(x,y) for x in VA for y in VA if x[2]==y[2] and x[0]==y[1] and x[1]==y[0]]
    VAVM = [(x,z) for x in VA for z in VM if x[0]==z[0]\
      and ((x[2]==0 and x[1] in z[1]) or (x[2]==1 and x[1] not in z[1]))]
    G = Graph()
    G.add_vertices(VA+VM)
    G.add_edges(VAVA+VAVM)
    return G
    
def twist( cfigrf, i, j):
    H = cfigrf.copy()
    H.delete_edges( [((i,j,0),(j,i,0)), ((i,j,1),(j,i,1))] )
    H.add_edges( [((i,j,0),(j,i,1)), ((i,j,1),(j,i,0))] )
    return H
    
######

# get the symmetric square of a graph
def symm_sq( grf):
    P=list(Subsets(grf.vertices(),2))
    n=len(P)
    return Graph([P, lambda pa,pb: len(pa.symmetric_difference(pb))==2 and grf.has_edge(pa.symmetric_difference(pb))])