# this function creates a dictionary with keys prop(grf) and values list of graphs;
# so if prop( grf) is lambda g: g.charpoly(), the keys will be characteristic polynomials
# and the values will be the list of graphs from grfs with that particular characteristic polynomial
def prop_sets_sd( grfs, prop):
	dc = {}
	for G in grfs:
	    dc.setdefault( prop(G), []).append(G)
	return dc

# values of dc must be lists
def append_value(dc, key, value):
    # Check if key is in dict or not
    if key in dc:
        dc[key].append(value)
    else:
        dc[key] = [value]
    return dc

# this is equivalent to prop_sets_sd, but it is possible to understand the code
def prop_sets(grfs, prop):
    dc = {}
    for G in grfs: append_value(dc,prop(G),G)
    return dc
    
# given a dictionary, this creates a new dictionary consisting of the pairs (key,value)
# where the value satisfies the given predicate
#
# {k:v for (k,v) in dc.items() if pred(v)}
#
def dc_filter( dc, pred):
    dc2 = {}
    for k,v in dc.items():
        if pred(v): dc2[k] = v
    return dc2 
    
def coprop( grfs, prop):
    dc = prop_sets( grfs, prop)   
    dc = dc_filter( dc, lambda it: len(it)>1)
    return dc

	
###############################################################################

def is_asymmetric( X ):
	return (1 == X.automorphism_group( order=True, return_group=False) )
def orbits( X):
    return X.automorphism_group(return_group=False, orbits=True)
    
# returns order and orbits of vertex stabilizer
# NOTE: If G = X.automorphism_group() then G.orbits() is wrong in many cases.
def vx_stab( X, vx):
    vxs = X.vertices()
    vxs.remove( vx)
    return X.automorphism_group(return_group=False, orbits=True, order=True, \
      partition=[[vx],vxs])

# not tested much; might not be efficient on long lists of graphs
def is_arc_transitive( X):
    if X.is_vertex_transitive():
        vxs = X.vertices()
        vx = vxs[0]
        vxs.remove(vx)
        orbs = X.automorphism_group(partition=[[vx],vxs], return_group=False, orbits=True)
        nbrs = Set(X[vx])
        return nbrs in map( Set, orbs)
    else:
        return False

###############################################################################

# get walk matrix relative to V(X)
def walk_mat( G): 
	A = G.am()
	v = G.num_verts()
	jv = vector([1 for i in [0..v-1]]) # jv = vector(n*[1])
	W = A.iterates(jv, v, rows = False)
	return W

# get walk matrix relative to a subset of grf.vertices(), subset read in as a list	
def sub_walk( grf, vxs):
	n = grf.num_verts()
	b = vector([0 for ii in [1..n]])
	for vx in vxs: 
		b[ vx]=1
	A = grf.am()
	return A.iterates(b, n, rows = False)
	
# sub_walk should be extended to handle this case gracefully
def vx_walk( grf, vx):
	return sub_walk( grf, vx.list())


###############################################################################

def dist_to_set(X, v, vxs):
    ds = [X.distance(v, w) for w in vxs]
    return min(ds) 
    
def distance_partition(X,vxs):
    E = {}
    for v in X.vertices():
        dj = dist_to_set(X,v,vxs)
        E.setdefault(dj,[]).append(v)
    return E.values()

## this will fail under Python 3 - cannot use 'map'    	
def contract(H, e):
	G = H.copy()
	u, v = e[0], e[1]
	vedges = map(lambda x: [v,x], G.neighbors(u))
	G.delete_vertex(u)
	G.add_edges(vedges)
	return G
	
def subdiv(G):
	V = G.vertices()
	nv = len(V)
	E = G.edges(labels = False)
	ne = len(E)
	dc ={}
	for i in [0..nv-1]:
		for j in [0..ne-1]:
			if V[i] in E[j]:
				dc.setdefault(i,[]).append(j+nv)
	return Graph( dc)

# subdivide each edge k times, i.e., replace it by a path of length k+1   
def ksub( grf, k):
    def pedges( arc, k):
        op = []
        i = arc[0]
        j = arc[1]
        op = op + [(i,(arc,1))] + [(j,(arc,k))]
        if k>1:
            qq= [((arc,i),(arc,i+1)) for i in [1..k-1]]
            op = op + qq
        return op
    G = Graph()
    grf.relabel()
    arcs = [(i,j) for i in grf.vertices() for j in grf.vertices()\
      if i<j and (i,j) in grf.edges(labels=False)]
    pathvs = [(aa,ii) for aa in arcs for ii in [1..k]]
    G.add_vertices( grf.vertices()+pathvs)
    for arc in arcs:
        G.add_edges( pedges(arc, k))
    return G

def distance_matrix( G): # G.distance_matrix()
	v = G.num_verts()
	D = Matrix( v, v, 0) 
	H = G.copy()
	H.relabel()
	for i in [0..v-1]:
		for j in [0..v-1]:
			D[i,j] = H.distance(i,j)
	return D

# 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

# switching graph of grf1; grf1 and grf2 are switching equivalent
# if and only if swgrf( grf1).is_isomorphic( swgrf( grf2))		
def swgrf( grf):
	V = [(vx,0) for vx in grf.vertices()]+[(vx,1) for vx in grf.vertices()]
	def adj( va, vb):
		if va[0] in grf[vb[0]] and va[1]==vb[1]: return True
		if va[0] not in grf[vb[0]] and va[0] != vb[0] and va[1] != vb[1]:
			return True
		else: return False
	return Graph([V, lambda va, vb: adj( va, vb)], loops=False)
	
def sw_mat( grf):
    A = grf.am()
    v = grf.num_verts()
    K = matrix(QQ,n,n,n^2*[1])-1 # B = ones_matrix(n)-1-A
    B = K-A
    C = matrix([A,B,B,A])
    return C

def cone( grf):
	cg = grf.complement()
	cg.add_vertex()
	return cg.complement()
	
def onesum(X,i,Y,u):
	n = X.num_verts()
	X1 = X.copy()
	X1.relabel()
	Y1 = Y.copy()
	Y1.relabel()
	dc = dict([(w,w+n) for w in Y.vertices()])
	dc[u] = i 
	Y1.relabel(dc)
	return X1.union(Y1)
	
def twosum(X,i,j, Y,u,v): 
	n = X.num_verts()
	X1 = X.copy()
	X1.relabel()
	Y1 = Y.copy()
	Y1.relabel()
	dc = dict([(w,w+n) for w in Y.vertices()])
	dc[u] = i; dc[v] = j 
	Y1.relabel(dc)
	return X1.union(Y1)


##################################################################################
##  Matrices
##################################################################################
	
# delete i-th row and j-th column from matrix M
def dropij( M,i,j):
    return M.delete_rows([i]).delete_columns([j])

# delete i-th row and column
def dropi( M,i):
    return dropij(M,i,i)

# delete first row and column	
def drop1( M):
    M.subdivide(1,1)
    return M.subdivision(1,1)
            
def psinv( mm):
	# mm had better be a real symmetric matrix
	eps = 1e-6
	n = mm.nrows()
	nn = Matrix(n,n,0)
	sp = spcdec( mm)
	for pr in sp:
		if abs(pr[0]) > eps:
			nn += (1/pr[0])*pr[1]
	return Matrix(nn)
   
##################################################################################

def psiuv(X, u, v):
    g = X.characteristic_polynomial()
    Y = X.copy()
    Y.delete_vertex(u)
    gu = Y.characteristic_polynomial()
    Y = X.copy()
    Y.delete_vertex(v)
    gv = Y.characteristic_polynomial()
    Y = X.copy()
    Y.delete_vertices([u,v])
    guv = Y.characteristic_polynomial()
    phiphi = gu*gv - g*guv
    return phiphi.sqrt()
	
# g is X.characteristic_polynomial() and f = psiuv(X,u,v) and la is an eigenvalue
# then f(la)/g'(la) is (E_la)_{u,v}

##################################################################################
'''
This code is intended to read and write binary matrices.
In the file each matrix consists of an integer n (number of rows)
followed by n lines of equal length, each of which is a non-empty 01-string. 
Blank lines between matrices should cause no problem. Lines between matrices 
that start with '#' are also ignored (and so could contain comments).
'''

# it does not like blank lines at the end of the file though

# convert 01-vector to a binary string and back
def bv_to_str( vec):
	return ''.join( [num.str() for num in vec])
def str_to_bv( ss):
	return vector( map( int, list(ss)))

# h = open('teds.txt', 'r')
# lines = h.readlines()
# h.close()
	
def get_matrices( lines):
	mats = []
	while lines != []:
		line = lines.pop(0).strip()
		while line == '' or line[0] == '#':
			line = lines.pop(0).strip()
		n = int( line)
		A = matrix( map( lambda row: str_to_bv( row.strip()), lines[:n]))
		mats += [A]
		del lines[:n]
	return mats

# mats is a list of binary matrices, op will be a list of strings
def strmats( mats):
	op = []
	while mats != []:
		M = mats.pop(0)
		nr = M.nrows()
		mat_strings = map( bv_to_str, M.rows())
		op = op + [str(nr)] + mat_strings
	return op

def write_strs( string_list, fh):
	for it in string_list:
		h.write( it+'\n')
        
##################################################################################
	
def is_walkreg( grf):
    cpd = grf.charpoly().derivative()
    vxs = grf.vertices()
    n = grf.num_verts()
    for i in vxs:
        gi = grf.copy()
        gi.delete_vertex(i)
        cpi = gi.charpoly()
        if n*cpi != cpd: return False
    return True
    
# get all maximal paths in grf starting at vx
def find_paths( grf, vx):
    paths=[]
    wrk=[[vx]]
    while wrk!=[]:
        pth = wrk.pop()
        nbrs = []
        for vx in grf[pth[0]]:
            if vx not in pth:
                nbrs += [vx]
        if nbrs==[]: paths += [pth]
        else:
            for vx in nbrs: wrk += [[vx]+pth]
    return paths
    
def path_tree( grf, vx):
    paths=[]
    edges=[]
    wrk=[[vx]]
    while wrk!=[]:
        pth = wrk.pop()
        nbrs = []
        for vx in grf[pth[0]]:
            if vx not in pth:
                nbrs += [vx]
        if nbrs==[]: paths += [pth]
        else:
            for vx in nbrs: 
                new_path = [vx]+pth
                wrk += [new_path]
                edges += [(pth,new_path)]
    return Graph(edges)

def line_dg( grf):
    return DiGraph([arcs(grf), lambda va,vb: va[1]==vb[0]])

# strict line digraph
def sline_dg( grf):
    return DiGraph([arcs(grf), lambda va,vb: va[1]==vb[0] and vb[1]!=va[0]])