#!/usr/bin/env python
# encoding: utf-8
"""
SmallCaseChecker.py

Created by Jamie Radcliffe on 2011-9-11.

Verifies small cases of the results in "Extremal Graphs for Homomorphisms II", Jonathan Cutler and A.J. Radcliffe,
submitted to the Journal of Graph Theory.
"""

def binomial(n,k):
	"""binomial coefficient"""
	if n<0 or k<0:
		raise ValueError
	if k > n:
		return 0
	if n==0:
		return 1
	if k == 0 or n == k:
		return 1
	if k == 1:
		return n
	if k == 2:
		return n*(n-1)/2
	return binomial(n-1,k)+binomial(n-1,k-1)

def code_add(n,code,domq):
	"""Returns the code for the graph K_1 union/join T(n-1,code); domq == True if it's a join"""
	return code + (2**(n-1) if domq else 0)
		
def H_add(n,H,domq):
	"""Returns the number of homs when an isloated/dominating vertex is added to an n-1 vertex graph with H homs;
	domq == True if it's dominating"""
	return H + 2**(n-1)+1 if domq else 3*H

class Threshold(object):
	"""Describes a threshold graph"""
	def __init__(self,n,code):
		super(Threshold, self).__init__()
		if code >= 2**n or code < 0:
			raise ValueError, "Code out of range in Threshold.init"
		self.code = code
		self.n = n
		self._H = None
		self._e = None
	
	@property
	def H(self):
		"""Computes the number of homomorphisms from Threshold(seq) to the fox"""
		if self._H is not None:
			return self._H
		homs = 1
		code = self.code
		# homs accumulates the number of homomorphisms from the graph T(i,code') where code' is the tail
		# end of self.code. 
		for i in xrange(self.n):
			if code % 2:
				homs = homs + 2**i + 1
			else:
				homs = 3*homs
			code //= 2
		self._H = homs
		return homs
		
	@property
	def e(self):
		"""Computes the number of edges in Threshold(n,code)"""
		sofar = 0
		code = self.code // 2
		# sofar accumulates the number of edges in T(i,code') where code' is the tail
		# end of self.code
		for i in xrange(1,self.n):
			if code % 2:
				sofar += i
			code //= 2
		return sofar
		
	def __eq__(self,other):
		"""Test for equality of graphs"""
		if not isinstance(other,Threshold) or self.n != other.n:
			return False
		return self.code // 2 == other.code // 2
		
	def __ne__(self,other):
		"""Test for inequality of graphs"""
		return not (self == other)
		
	def __str__(self):
		"""Return printable form of self"""
		return "{1:0{0}b}".format(self.n,self.code)
		
	def addiso(self):
		"""Adds an isolated vertex"""
		return Threshold(self.n+1,code_add(self.n+1,self.code,False))
		
	def adddom(self):
		"""Adds a dominating vertex"""
		return Threshold(self.n+1,code_add(self.n+1,self.code,True))
		
		
def split_graph_code(n,k):
	"""Returns the code for K_k v E_(n-k) with the complete side on the left"""
	return 2**(n-k) * (2**(k)-1)
	
def lex(n):
	"""Generator for codes of lex graphs by number of edges"""
	# done = No. of vertices of full degree; same as vertex number currently being added to
	# joinedto =  No. of vertices to the right that vertex done is currently joined to
	yield 0
	for done in xrange(n):
		for joinedto in xrange(1,n - done):
			yield split_graph_code(n,done) + 2**joinedto

def colex(n):
	"""Generator for codes of colex grapns by numbers of edges"""
	# done = No. of vertices in a clique already
	# joinedto = No of vertices in the clique that vertex done is joined to
	yield 0
	for done in xrange(1,n):
		for joinedto in xrange(1,done+1):
			yield 2**(done+1) - 1 - (2**(done-joinedto) + 1 if joinedto < done else 1)
		
def anti_triangles(n):
	"""Generator for codes of anti-triangle graphs by size of complete side"""
	for q in xrange(4,n+1):
		yield 2**q-8

def kjes(n):
	"""Generator for codes of KJE's"""
	# done = Number of vertices in the K side
	# joinedto = number from done that the extra vertex is joined to
	for done in xrange(2,n-2):
		for joinedto in xrange(1,done):
			yield split_graph_code(n,done+1) - 2**(n-joinedto-1)

def kjeis(n):
	"""Generator for codes of KJEIs"""
	for code in kjes(n-1):
		yield code
	for k in xrange(1,n-2):
		yield split_graph_code(n-1,k)
		
def optimal_pair(n,e,oldpairs):
	"""Generates the optimal (code,H) pair for an n vertex graph with e edges, based on the complete list of optimal pairs 
	for n-1."""
	if n < 1:
		raise ValueError
	if e < n-1:
		domq = False
		h0 = H_add(n,oldpairs[e][1],False)
	elif e > binomial(n-1,2):
		domq = True
		h1 = H_add(n,oldpairs[e-n+1][1],True)
	else:
		h0 = H_add(n,oldpairs[e][1],False)
		h1 = H_add(n,oldpairs[e-n+1][1],True)
		if h1 == h0:
			print "Tie between {0} and {1}".format(Threshold(n,code_add(n,oldpairs[e-n+1][0],True)),Threshold(n,code_add(n,oldpairs[e][0],False)))
		domq = (h1>h0)
	if domq:
		return (code_add(n,oldpairs[e-n+1][0],True),h1)
	else:
		return (code_add(n,oldpairs[e][0],False),h0) 

def optimal_list_from_previous(n,oldpairs):
	"""Generates the list of optimal (code,H) pairs, for each value of e, using the information in oldpairs"""
	newpairs = tuple(optimal_pair(n,e,oldpairs) for e in xrange(binomial(n,2)+1))
	# print "Got newpairs={0}".format(newpairs)
	return newpairs	
	
def optimal_list(n):
	"""Returns a tuple of pairs consisting of a wr-optimal code and the corresponding H value, one for each
	possible value of e."""
	oldpairs = ((0,3),)
	for i in xrange(2,n+1):
		oldpairs = optimal_list_from_previous(i,oldpairs)
	return oldpairs

def legal_optimal_code(e, code, lexs, colexs, anticodes, kjecodes, kjeicodes):
	"""Determines whether this is an expected optimal code"""
	if code == lexs[e]:
		return True
	if code == colexs[e]:
		return True
	if e in anticodes and code == anticodes[e]:
		return True
	if e in kjecodes and code == kjecodes[e]:
		return True
	if e in kjeicodes and code == kjeicodes[e]:
		return True
	return False
	
def test_optimals(n):
	"""Checks that all the optimal (code,H) pairs for orders at most n pass test."""
	oldpairs = ((0,3),)
	failed = False
	for i in xrange(2,n+1):
		oldpairs = optimal_list_from_previous(i,oldpairs)
		anticodes = dict((Threshold(i,code).e,code) for code in anti_triangles(i))
		kjecodes = dict((Threshold(i,code).e,code) for code in kjes(i))
		kjeicodes = dict((Threshold(i,code).e,code) for code in kjeis(i))
		lexs = [code for code in lex(i)]
		colexs = [code for code in colex(i)]
		for e, (code, H) in enumerate(oldpairs):
			if not legal_optimal_code(e,code,lexs,colexs,anticodes,kjecodes,kjeicodes):
				print "T({0},{1}) failed test".format(i,Threshold(i,code))
				failed = True
	return not failed

def verify_optimals(n):
	"""Checks that the optimal (code,H) pairs returned by optimal_list(k) for k<=n agree with those given but the """
	oldpairs = ((0,3),)
	failed = False
	for i in xrange(2,n+1):
		oldpairs = optimal_list_from_previous(i,oldpairs)
		if not oldpairs == crude_optimal_list(i):
			print "Found discrepancy for i = {0}".format(i)
			print oldpairs
			print crude_optimal_list(i)
			failed = True
	return not failed
	
def crude_optimal_pair(n,e):
	"""Does the brute force optimization to find the optimal code and H for graphs with n vertices and e edges"""
	Ts = [Threshold(n,code) for code in xrange(2**n)]
	pairs = [(T.code,T.H) for T in Ts if T.e == e]
	return max(pairs,key=lambda p: p[1])
	
def crude_optimal_list(n):
	"""Returns the list of optimal (code,H) pairs determined crudely"""
	return tuple(crude_optimal_pair(n,e) for e in xrange(binomial(n,2)+1))
		
def main():
	
	# n = 15
	# if verify_optimals(n):
	# 	print "Verified optimal graphs are in fact optimal for n <= {0}".format(n)
	# else:
	# 	print "Failure to verify that optimal graphs are in fact optimal for n <={0}".format(n)
	
	n = 250
	if test_optimals(n):
		print "Verified all optimal graphs are in the five classes for n <= {0}".format(n)
	else:
		print "Failure to verify that all optimal graphs are in the five classes for n <= {0}".format(n)
		
if __name__ == '__main__':
	main()