"""

Copyright (c) 2011, E. R. Vaughan. All rights reserved.

Redistribution and use in source and binary forms, with or without modification,
are permitted provided that the following conditions are met:

1) Redistributions of source code must retain the above copyright notice, this
list of conditions and the following disclaimer.

2) Redistributions in binary form must reproduce the above copyright notice,
this list of conditions and the following disclaimer in the documentation and/or
other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR
ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

"""

import sys
import fractions
import getopt
import json
import itertools

try:
	from sage.all import *
	using_sage = True
except ImportError:
	using_sage = False

should_print_help = False

try:
	opts, args = getopt.gnu_getopt(sys.argv[1:], "", ["help", "admissible-graphs",
		"flags", "r-matrices", "qdash-matrices",
		"pair-densities", "q-matrices", "verify-bound",
		"sharp-graphs", "flag-algebra-coefficients"])

except getopt.GetoptError:
	should_print_help = True
	opts, args = ((), ())

if len(args) != 1:
	should_print_help = True

action = ""

for o, a in opts:
	if o == "--help":
		should_print_help = True
	elif o == "--admissible-graphs":
		action = "print admissible graphs"
	elif o == "--flags":
		action = "print flags"
	elif o == "--r-matrices":
		action = "print r_matrices"
	elif o == "--qdash-matrices":
		action = "print qdash_matrices"
	elif o == "--pair-densities":
		action = "print pair densities"
	elif o == "--q-matrices":
		action = "print q_matrices"
	elif o == "--verify-bound":
		action = "verify bound"
	elif o == "--sharp-graphs":
		action = "print sharp graphs"
	elif o == "--flag-algebra-coefficients":
		action = "print flag algebra coefficients"

if should_print_help:
	sys.stdout.write("""Usage: inspect_certificate.py CERTIFICATE OPTIONS
Possible options:
--help                       Display this message.
--admissible-graphs          Display the admissible graphs.
--flags                      Display the types and flags.
--r-matrices                 Display the R matrices.
--qdash-matrices             Display the Q' matrices.
--q-matrices                 Compute and display the Q matrices.
--pair-densities             Compute and display the flag pair densities.
--verify-bound               Verify the bound.
--sharp-graphs               Display the admissible graphs that are sharp.
--flag-algebra-coefficients  Display each admissible graph's flag algebra coefficient.
"""
	)
	sys.exit(0)

certificate_filename = args[0]

try:
	if certificate_filename[-3:] == ".gz":
		import gzip
		certf = gzip.open(certificate_filename)
	elif certificate_filename[-4:] == ".bz2":
		import bz2
		certf = bz2.BZ2File(certificate_filename)
	else:
		certf = open(certificate_filename)
except IOError:
	sys.stdout.write("Could not open certificate.\n")
	sys.exit(1)

certificate = json.load(certf)

sys.stdout.write("Problem is \"%s\".\n" % certificate["description"])
sys.stdout.write("Claimed bound is %s.\n" % certificate["bound"])

minimize = "minimize" in certificate["description"]

if using_sage:
	x = polygen(QQ)
else:
	x = None

if "field" in certificate.keys():

	if certificate["field"] != "RationalField()":
	
		sys.stdout.write("Field used is \"%s\".\n" % certificate["field"])
	
		if not using_sage:
			sys.stdout.write("This script must be run using Sage's python, as it does not use the rational field.\n")
			sys.exit(1)
		
		try:
			K = sage_eval(certificate["field"], locals={'x': x})
		except AttributeError:
			K = RationalField()
		x = K.gen()
		
if action == "":
	sys.exit(0)

if "3-graph" in certificate["description"]:
	edge_size = 3
elif "2-graph" in certificate["description"]:
	edge_size = 2
else:
	sys.stdout.write("Unsupported graph kind.\n")
	sys.exit(1)

oriented = "oriented" in certificate["description"]

def induced_subgraph (g, S):

	global oriented
	good_edges = tuple(e for e in g[1] if all(x in S for x in e))
	p = [0 for i in range(g[0] + 1)]
	for i in range(len(S)):
		p[S[i]] = i + 1
	if oriented:
		edges = tuple(sorted([tuple([p[x] for x in e]) for e in good_edges]))
	else:
		edges = tuple(sorted([tuple(sorted([p[x] for x in e])) for e in good_edges]))
	return (len(S), edges)

def minimal_typed_isomorph (g, t):

	global oriented
	s = t[0]
	n = g[0]
	min_edges = g[1]
	for p in itertools.permutations(range(s + 1, n + 1), n - s):
		pt = tuple(range(1, s + 1)) + p
		if oriented:
			edges = tuple(sorted([tuple([pt[x - 1] for x in e]) for e in g[1]]))
		else:
			edges = tuple(sorted([tuple(sorted([pt[x - 1] for x in e])) for e in g[1]]))
		if edges < min_edges:
			min_edges = edges
	return (n, min_edges)

def graph_to_string (g):
	s = "%d:" % g[0]
	for e in g[1]:
		for x in e:
			s += "%d" % x
	return s

def string_to_graph (s):
	global edge_size
	n = int(s[0])
	if n < 0 or n > 9 or s[1] != ":":
		return None
	e = len(s) - 2
	if e == 0:
		return (n, ())
	if e % edge_size != 0:
		return None
	m = e / edge_size
	edges = []
	for i in range(0, m):
		# N.B. +2 because of n: header
		edge = [int(s[(edge_size * i) + 2 + j]) for j in range(edge_size)]
		es = set(edge)
		if len(es) != edge_size or not es <= set(range(1, n + 1)):
			return None
		edges.append(tuple(edge))
	
	return (n, tuple(edges))

admissible_graphs = [string_to_graph(s) for s in certificate["admissible_graphs"]]
nag = len(admissible_graphs)
types = [string_to_graph(s) for s in certificate["types"]]
flags = [[string_to_graph(s) for s in f] for f in certificate["flags"]]

if action == "print admissible graphs":
	
	sys.stdout.write("There are %d admissible graphs:\n" % nag)
	for i in range(nag):
		sys.stdout.write("%d. %s\n" % (i + 1, graph_to_string(admissible_graphs[i])))
	sys.exit(0)

if action == "print flags":
	for t in range(len(types)):
		nf = len(flags[t])
		sys.stdout.write("Type %d (%s) has %d flags:\n" % (t + 1, graph_to_string(types[t]), nf))
		for i in range(nf):
			sys.stdout.write("   %d. %s\n" % (i + 1, graph_to_string(flags[t][i])))
	sys.exit(0)

def to_string (a):
	if a == None:
		return "Not used."
	if isinstance(a, list):
		return map(to_string, a)
	s = str(a)
	try:
		n = int(s)
		return n
	except ValueError:
		return s

if action == "print r_matrices":

	for t in range(len(types)):
		sys.stdout.write("R matrix for type %d (%s):\n" % (t + 1, graph_to_string(types[t])))
		sys.stdout.write("%s\n" % to_string(certificate["r_matrices"][t]))
	sys.exit(0)

if action == "print qdash_matrices":

	for t in range(len(types)):
		sys.stdout.write("Q' matrix for type %d (%s):\n" % (t + 1, graph_to_string(types[t])))
		sys.stdout.write("%s\n" % to_string(certificate["qdash_matrices"][t]))
	sys.exit(0)

n = certificate["order_of_admissible_graphs"]
vertices = tuple(range(1, n + 1))

sys.stdout.write("Computing Q matrices...\n")

Qs = []

for t in range(len(types)):

	if certificate["qdash_matrices"][t] == None:
		Qs.append(None)
		continue

	if using_sage:
		QD = [[sage_eval(str(s), locals={'x': x}) for s in row] for row in certificate["qdash_matrices"][t]]	
	else:
		QD = [[fractions.Fraction(s) for s in row] for row in certificate["qdash_matrices"][t]]	
	
	if certificate["r_matrices"][t] == None:
		Qs.append(QD)
		continue

	if using_sage:
		R = [[sage_eval(str(s), locals={'x': x}) for s in row] for row in certificate["r_matrices"][t]]
	else:
		R = [[fractions.Fraction(s) for s in row] for row in certificate["r_matrices"][t]]
	
	nq = len(QD)
	nf = len(flags[t])
	
	Q = [[0 for j in range(i, nf)] for i in range(nf)]
	for l in range(nq):
		for k in range(l, nq):
			qlk = QD[l][k - l]
			if qlk != 0:
				for i in range(nf):
					for j in range(i, nf):
						Q[i][j - i] += R[i][l] * R[j][k] * qlk
						if k != l:
							Q[i][j - i] += R[i][k] * R[j][l] * qlk
	Qs.append(Q)
	

if action == "print q_matrices":

	for t in range(len(types)):
		sys.stdout.write("Q matrix for type %d (%s):\n" % (t + 1, graph_to_string(types[t])))
		sys.stdout.write("%s\n" % to_string(Qs[t]))
	sys.exit(0)

sys.stdout.write("Computing pair densities...\n")
pair_densities = {}

for t in range(len(types)):
	s = types[t][0]
	m = flags[t][0][0] # assumes all flags are same order
	nf = len(flags[t])
	for g in admissible_graphs:
		pairs_found = [[0 for j in range(k, nf)] for k in range(nf)]
		total = 0
		for p in itertools.permutations(vertices):
			if p[s] > p[m]:
				continue
			is_good_permutation = True
			for c in range(s, n):
				if c == s or c == m:
					continue
				if p[c - 1] > p[c]:
					is_good_permutation = False
					break
			if not is_good_permutation:
				continue
			total += 1
			
			it = induced_subgraph(g, p[:s])
			if it != types[t]:
				continue
			
			def induced_flag(g, tv, ov):
				type_vertices = list(tv)
				other_vertices = list(ov)
				edges = []
				for e in g[1]:
					if any(not (x in tv or x in ov) for x in e):
						continue
					edge = []
					for x in e:
						if x in tv:
							edge.append(1 + type_vertices.index(x))
						else:
							edge.append(len(tv) + 1 + other_vertices.index(x))
					edges.append(tuple(edge))
				return (len(tv + ov), tuple(edges))		
		
			f1 = minimal_typed_isomorph(induced_flag(g, p[:s], p[s:m]), types[t])
			f2 = minimal_typed_isomorph(induced_flag(g, p[:s], p[m:2*m-s]), types[t])

			if1 = flags[t].index(f1)
			if2 = flags[t].index(f2)
			
			if if1 <= if2:
				pairs_found[if1][if2 - if1] += 1
			else:
				pairs_found[if2][if1 - if2] += 1

		for k in range(nf):
			for j in range(k, nf):
				pf = pairs_found[k][j - k]
				if pf > 0:
					if j == k:
						pf *= 2
					if using_sage:
						pair_densities[(t, admissible_graphs.index(g), k, j)] = Integer(pf) / (total * 2)
					else:
						pair_densities[(t, admissible_graphs.index(g), k, j)] = fractions.Fraction(pf, total * 2)

if action == "print pair densities":

	for i in range(nag):
		sys.stdout.write("Pair densities for admissible graph %d (%s):\n" % (i + 1, graph_to_string(admissible_graphs[i])))

		for t in range(len(types)):
			sys.stdout.write("   Non-zero densities for type %d (%s):\n" % (t + 1, graph_to_string(types[t])))
		
			for key, d in pair_densities.iteritems():
				if key[0] != t or key[1] != i:
					continue
				sys.stdout.write("      Flags %d and %d (%s and %s): %s\n" % (key[2] + 1, key[3] + 1,
					graph_to_string(flags[t][key[2]]), graph_to_string(flags[t][key[3]]), d))
	sys.exit(0)

sys.stdout.write("Computing bound...\n")

if using_sage:
	bounds = [sage_eval(str(s), locals={'x': x}) for s in certificate["admissible_graph_densities"]]
else:
	bounds = [fractions.Fraction(s) for s in certificate["admissible_graph_densities"]]

for key in pair_densities.keys():
	t, i, j, k = key
	if Qs[t] == None: # check that type is used
		continue
	d = pair_densities[key]
	v = Qs[t][j][k-j]
	if minimize:
		v *= -1
	if j == k:
		bounds[i] += d * v
	else:
		bounds[i] += d * v * 2

if minimize:
	bound = min(bounds)
else:
	# sage seems to be unable to return max of a list containing 0!
	nonzero_bounds = [b for b in bounds if b != 0]
	
	if len(nonzero_bounds) == 0:
		bound = 0
	else:
		bound = max(nonzero_bounds)

sys.stdout.write("Bound is %s.\n" % bound)

if action == "print sharp graphs":

	sys.stdout.write("Sharp graphs are:\n")
	for i in range(nag):
		if bounds[i] == bound:
			sys.stdout.write("%d. %s\n" % (i + 1, graph_to_string(admissible_graphs[i])))
	sys.exit(0)

if action == "print flag algebra coefficients":

	sys.stdout.write("There are %d admissible graphs:\n" % nag)
	for i in range(nag):
		sys.stdout.write("%d. (%s) : %s\n" % (i + 1, graph_to_string(admissible_graphs[i]), bounds[i]))
	sys.exit(0)