import numpy as np
# from math import cos, sin, pi, sqrt, floor
import mpmath as mp
from mpmath import cos, sin, pi, sqrt, floor
from vectors import *

inf = float('inf')

def singleIneq(a, b): ### need to make it show nothing when |ai|=1
	s = ""
	n = len(a)
	found_nonzero = False
	for i in range(len(a)):
		if i == 0 and a[0] != 0:
			s += str(a[i]) + " " + chr(122-n+1) + " "
			found_nonzero = True
		else:
			ai = a[i]
			if ai > 0:
				if found_nonzero:
					s += "+ "
				else:
					found_nonzero = True
				s += str(ai) + " " + chr(122-n+1+i) + " "
			elif ai < 0:
				if not found_nonzero:
					s += str(ai) + " " + chr(122-n+1+i) + " "
					found_nonzero = True
				else:
					s += "- " + str(-ai) + " " + chr(122-n+1+i) + " "
	s += "<= " + str(b) + "\n"
	return s

class Interval:
	def __init__(self, a, b=None):
		self.L = -inf
		self.U = inf
		self.type = "Whole"
		if type(a) is tuple or type(a) is list:
			assert len(a) == 2, "not a valid tuple"
			self.L, self.U = a[0], a[1]
			if self.L == -inf and self.U == inf:
				self.type == "Whole"
			else:
				self.type = "Range"
		else:
			# assert (type(a) is float or type(a) is int or type(a) is np.int64 or type(a) is np.float64) and (type(b) is float or type(b) is int is type(b) is np.int64 or type(b) is np.float64), "not valid bounds"
			if a == 0 and b < 0:
				self.L, self.U = None, None
				self.type = "Empty"
			if a > 0:
				self.L, self.U = -inf, (1.*b)/a
				self.type = "Range"
			if a < 0:
				self.L, self.U = (1.*b)/a, inf
				self.type = "Range"
	
	def __str__(self):
		if self.type == "Empty":
			return "\emptyset"
		s = ""
		if self.L == -inf:
			s += "(-inf, "
		else:
			s += "[" + str(self.L) + ", "
		if self.U == inf:
			s += "inf)"
		else:
			s += str(self.U) + "]"
		return s
	
	def intersect(self, other):
		if self.type == "Whole" or other.type == "Empty":
			return other
		if other.type == "Whole" or self.type == "Empty":
			return self
		a, b, c, d = self.L, self.U, other.L, other.U
		if b < c or d < a:
			return Interval(0,-1)
		if a <= c <= b:
			return Interval((c,min(b,d)))
		if a <= d <= b:
			return Interval((max(a,c),d))
		return self
	
	def intersectZ(self, n=1):
		if self.type == "Whole":
			return list(range(n))
		if self.L == -inf:
			return list(range(int(floor(self.U-n)), int(floor(self.U))))
		if self.U == inf:
			return list(range(int(ceil(self.L)), int(ceil(self.L+n))))
		if self.type == "Empty":
			return []
		return list(range(int(ceil(self.L)), int(floor(self.U)+1)))

"""
t = [ Interval(0,1), Interval(0,-1), Interval((-1,1)), Interval((-1,1),2), Interval(2,5), Interval(-3,7) ]
for i in range(len(t)-1):
	i1, i2 = t[i], t[i+1]
	print str(i1), str(i2), str(i1.intersect(i2)), i1.intersectZ(7)
"""

def hand_coded_det(M):
	n = len(M)
	assert 1 <= n <= 2, "higher dimensions not supported yet (numerical issues with numpy/mpmath)"
	if n == 1:
		return M[0]
	else:
		return M[0][0]*M[1][1] - M[0][1]*M[1][0]

def hand_coded_adj(M):
	n = len(M)
	assert 1 <= n <= 2, "higher dimensions not supported yet (numerical issues with numpy/mpmath)"
	if n == 1:
		return [[mp.mpf(1)]]
	else:
		return np.array([[M[1][1], -M[0][1]], [-M[1][0], M[0][0]]])

def hand_coded_inv(M):
	d = hand_coded_det(M)
	assert d != 0, "singular matrix!"
	return hand_coded_adj(M)/(1.0*d)

"""
print hand_coded_inv([[1,2],[3,4]])
print hand_coded_inv([[mp.mpf(1),mp.mpf(2)],[mp.mpf(3),mp.mpf(4)]])
print hand_coded_inv([[mp.mpf(1),mp.mpf(2)],[mp.mpf(2),mp.mpf(4)]])
"""

class HyperplaneSolid:
	def __init__(self, n=0):
		assert n >= 0, "negative dimension doesn't exist"
		self.dim = n
		self.vecs = np.array([])
		self.bounds = np.array([])
		self.count = 0
		self.corners = []
		self.V = 0
	
	def __str__(self):
		assert self.dim <= 26, "i'm sorry Dave, I'm afraid I can't do that..."
		s = ""
		for i in range(self.count):
			s += singleIneq(self.vecs[i], self.bounds[i])
		return s[:-1]
	
	def enter(self, a, b): # want to add the constraint a.x <= b. for now, a will be a numpy array of floats
		assert len(a) == self.dim, "vector has the wrong dimension"
		if self.count == 0:
			self.vecs = np.array([a])
			self.bounds = np.array([b])
		else:
			self.vecs = np.append(self.vecs, [a], axis=0)
			self.bounds = np.append(self.bounds, [b])
		self.count += 1
		# print self.vecs, self.bounds, self.count
	
	def check(self, pt):
		assert len(pt) == self.dim, "point has the wrong dimension"
		for i in range(self.count):
			if np.dot(pt, self.vecs[i]) > self.bounds[i]:
				return False
		return True
	
	def select(self, entries): # given a list of n of the hyperplanes, returns the n-by-n matrix and corresponding bounds
		assert len(entries) == self.dim, "wrong number of hyperplanes for a potential corner"
		return (self.vecs[entries], self.bounds[entries])
	
	def findCorners(self):
		possibilities = []
		def combo(arr, n, r): # adapted from https://www.geeksforgeeks.org/print-subsets-given-size-set/
			def combo_aux(arr, n, r, index, data, i):
				if index == r:
					# print data
					data_copy = []
					for d in data:
						data_copy.append(d)
					possibilities.append(data_copy)
				elif i < n:
					data[index] = arr[i]
					combo_aux(arr, n, r, index+1, data, i+1), combo_aux(arr, n, r, index, data, i+1)
			d = list(range(r))
			combo_aux(arr, n, r, 0, d, 0)
		combo(list(range(self.count)), self.count, self.dim)
		# print possibilities, len(possibilities)
		verts = []
		for ind in possibilities:
			(M,B) = self.select(ind)
			if hand_coded_det(M) != 0:
				N = hand_coded_inv(M)
				pt = np.dot(N,B) # candidate: solution to M.x = B -> x = N.B
				bools = np.dot(self.vecs, pt) <= self.bounds
				bools[ind[0]] = True # for some reason there's a floating point error with the equality cases
				bools[ind[1]] = True
				if np.all(bools): # for each v.x <= b bounding the solid, is v.pt <= b?
					# return pt
					verts += [pt]
		# return verts
		self.corners = verts
		self.V = len(verts)
		# return None # empty intersection of half-spaces!
	
	def minimizer(self, a, b=0): # a is a np array and b is a float, representing a function a.x - b, but b is not actually used
		assert self.V > 0, "empty convex body!"
		v = self.corners[0]
		m = np.dot(v,a)
		for w in self.corners[1:]:
			n = np.dot(w,a)
			if n < m:
				v = w
				m = n
		return v
	"""
	def downDimension(self, basis, pt): # basis is a np array v1..n of a LLL basis of R^n, and pt is a np array of integers pt1..n representing pt1 v1 + ...
		assert np.shape(basis) == (self.dim, self.dim), "basis is for the wrong vector space"
		assert len(pt) == self.dim, "not a point in R^n"
		assert np.linalg.det(basis) != 0, "not a full-rank basis"
		tauinv = np.linalg.inv(np.transpose(basis))
		n = self.dim
		P = HyperplaneSolid(n-1)
		for i in range(len(self.vecs)):
			a = [ np.linalg.multi_dot([self.vecs[i], tauinv, basis[j]]) for j in range(n-1) ]
			b = self.bounds[i] - np.linalg.multi_dot([self.vecs[i], tauinv, pt])
			P.enter(a,b)
		return P
	"""
	def largeSimplex(self): # the underlying assumption here is that any n hyperplanes meet at a unique point if they meet at all
		n = self.dim
		assert self.V > n, "not enough vertices to form a proper simplex"

		"""----------- "hacky" version analogous to .findCorners() -----------"""
		possibilities = []
		def combo(arr, n, r): # adapted from https://www.geeksforgeeks.org/print-subsets-given-size-set/
			def combo_aux(arr, n, r, index, data, i):
				if index == r:
					# print data
					data_copy = []
					for d in data:
						data_copy.append(d)
					possibilities.append(data_copy)
				elif i < n:
					data[index] = arr[i]
					combo_aux(arr, n, r, index+1, data, i+1), combo_aux(arr, n, r, index, data, i+1)
			d = list(range(r))
			combo_aux(arr, n, r, 0, d, 0)
		combo(list(range(self.V)), self.V, n+1)
		max_volume = 0
		simplex = list(range(n+1))
		corners = np.array(self.corners)
		for e in possibilities:
			a = corners[e]
			a = (a-a[0])[1:]
			d = hand_coded_det(a)
			if d > max_volume:
				max_volume, simplex = d, e
		return simplex
		
		"""----------- version much closer to Lenstra's paper -----------"""
		"""simplex = range(n+1)
		if self.V == n+1:
			return simplex
		breakQ = False
		while not breakQ:
			breakQ = True
			for i in range(n+1):
				entries = simplex[:i] + simplex[(i+1):]
				m = np.array(self.corners)[entries]
				v = np.array([1,]*n)
				# print "~~~~~~~~~~~~~~~\n", entries, simplex, "\n", m, v,
				a = np.dot(hand_coded_inv(m), v)
				# print a
				for x in range(self.V):
					l = abs(np.dot(a,self.corners[x])-1)
					r = abs(np.dot(a,self.corners[simplex[i]])-1)
					if l > 1.5*r:
						simplex[i] = x
						# print "found a replacement! i =", i, ", x =", x, "and new simplex is", simplex
						breakQ = False
						break
				# don't need to add in the case of -g_i because |-g_i (x-v_j)| = |g_i (x-v_j)| and similalry for the RHS
		return simplex
		"""
	
	def ineqs(self):
		assert self.dim == 1, "can't make R-inequalities in high dimensions"
		# print self.vecs, self.bounds, type(self.vecs[0][0]), type(self.bounds[0])
		return [ Interval(self.vecs[k][0], self.bounds[k]) for k in range(self.count) ] ### need to test this!!!

	def Z1points(self):
		assert self.dim == 1, "wrong dimension"
		assert self.count > 0, "empty convex set"
		intervals = self.ineqs()
		i0 = intervals[0]
		# print "~~~~~~~~ running a Z1points instance with", len(intervals), "intervals"
		# print str(i0)
		for i in intervals[1:]:
			# print str(i),
			i0 = i0.intersect(i)
			# print str(i0)
		return [ [x] for x in i0.intersectZ() ]
	
	def Znpoints(self): # find a Z^self.dim-point, using Lenstra's algorithm
		assert self.dim > 0, "dimension too low"
		if self.dim == 1:
			return self.Z1points()
		n = self.dim
		self.findCorners()
		# print "corners are\n", self.corners
		if self.V <= n:
			return []
		simplex = self.largeSimplex()
		V = np.array([ [ self.corners[simplex[i]][j] for j in range(n) ] for i in range(n+1) ])
		# print "simplex is", simplex
		tauinv = np.transpose([ V[i] - V[n] for i in range(n) ])
		tau = hand_coded_inv(tauinv)
		# print "tauinv is\n", tauinv
		# print "tau is\n", tau
		p = [0.0]*n
		for i in range(n+1):
			p += V[i]
		p /= n+1.0
		badBasis = np.transpose(tau) # non-LLL basis, i.e. tau.d_i
		# print "badBasis is\n", badBasis
		basis = LLL(badBasis) # the LLL basis of {tau.d_i : i in [n]}
		# print "basis is\n", basis
		CoB_unrounded = np.dot(basis, hand_coded_inv(badBasis))
		CoB = [[round(y) for y in x] for x in CoB_unrounded]
		# print "CoB =\n", CoB
		q = np.dot(np.dot(hand_coded_inv(np.transpose(basis)), tau), p)
		pt = nearestInLattice(basis, q, q)
		# print "p is", p
		# print "q is", q
		# print "q is approximated in basis by", pt
		# print "corresponding badBasis entries are", np.dot(pt, CoB)
		v0 = np.dot(tau,V[n])
		AT = basis[:-1]
		A = np.transpose(AT)
		ATAinv = hand_coded_inv(np.dot(AT, A))
		h = np.linalg.norm(basis[-1]-np.linalg.multi_dot([A, ATAinv, AT, basis[-1]]))
		R = sqrt(n+0.0) # sqrt(n/(n+1.0)) ### figure out what the actual radius is
		# print "R =", R, "h =", h, "R/h =", R/h
		points = []
		for k in range(-int(floor(R/h)), int(floor(R/h))+1):
			new_pt = [0]*n
			for i in range(n-1):
				new_pt[i] = pt[i]
			new_pt[n-1] = pt[n-1] + k
			# print "~~~~~~~~~~~ trial new_pt is", new_pt, "~~~~~~~~~~~"
			# downDimension seems buggy. here is the "lazy" version: it only looks inside the simplex, instead of a translation of the HyperplaneSolid's constraints under tau.
			D = HyperplaneSolid(n-1)
			for i in range(n):
				w = [-1.0*basis[j][i] for j in range(n-1)]
				rhs1 = sum([pt[j]*basis[j][i] for j in range(n)])
				rhs = -1.0*v0[i] + rhs1 + k*basis[-1][i]
				D.enter(w,rhs)
			w_last = [sum(basis[j]) for j in range(n-1)]
			rhs_last = sum([pt[j]*sum(basis[j]) for j in range(n-1)])
			D.enter(w_last, 1 - rhs_last - (pt[n-1]+k)*sum(basis[-1]) + sum(v0))
			# print D
			lower_points = D.Znpoints()
			# print "  lower_points are", lower_points
			for p in lower_points:
				newZpoint = [0]*n
				for i in range(n):
					newZpoint[i] = new_pt[i]
				for i in range(n-1):
					newZpoint[i] += p[i]
				# print "   current newZpoint is", newZpoint
				points += [newZpoint]
		# print "points is", points
		## badPoints = [ np.linalg.multi_dot([p, CoB]) for p in points ]# convert the list points back into badBasis
		# badPoints is actually what we want!!! (i think)
		
		# print len(badPoints)
		# for p in badPoints:
		#	print p
		
		## return [[int(round(x)) for x in pt] for pt in badPoints]
		if len(points) == 0:
			return []
		return [[int(round(x)) for x in pt] for pt in np.dot(points, CoB) if self.check(pt)]

"""
a = HyperplaneSolid(2)
a.enter([-1,1],2)
a.enter([1,-1],2)
a.enter([1,1],2)
a.enter([-1,-1],2)
print(str(a))
a.findCorners()
print a.corners
print a.minimizer(np.array([3,2]))
print a.minimizer([2,3])
"""

"""
print "============ TESTING .Z1points() ============"
seg = HyperplaneSolid(1)
seg.enter([3],100.1)
seg.enter([-2],14)
print seg.Z1points() # it gets this right!
"""

"""
n = 7
ngon = HyperplaneSolid(2)
for k in range(n):
	pt = [cos(2*pi*k/n),sin(2*pi*k/n)]
	ngon.enter(pt,1)
print str(ngon)
ngon.findCorners()
print ngon.corners, ngon.V
print "============ TESTING .largeSimplex() ============"
print ngon.largeSimplex()
print np.array(ngon.corners)[ngon.largeSimplex()] # it gets this right!
"""


"""
n = 5
shift = [3,4]
ngon = HyperplaneSolid(2)
for k in range(n):
	pt = [0.5*cos(2*pi*k/n),0.5*sin(2*pi*k/n)]
	ngon.enter(pt,1+np.dot(shift,[0.5*cos(2*pi*k/n),0.5*sin(2*pi*k/n)]))
# print ngon
ngon.findCorners()
# print ngon.corners, ngon.V

print "============ TESTING .Znpoints() ============"
ngon.Znpoints()
# looks good!
"""


"""
print "============ TESTING .Z1points() ============"
line = HyperplaneSolid(1)
line.enter([1],5)
line.enter([-2],3)
line.enter([3],8)
line.enter([-1],0.6)
print str(line)
print line.Z1points()
# looks okay!
"""





"""
a = HyperplaneSolid(7)
a.enter([-3,-5,4,2.6,0,0,-3.99],2)
a.enter([0,-1,2,3,4,5,6],7)
print(str(a))
# print(singleIneq([3,-5,4,2.6,0,0,-3.99],2))
"""


""" # scratch code to get the list generation method for corner working
def combo2(arr, n, r, index, data, i):
	if index == r:
		# return data
		print data
	elif i < n:
		data[index] = arr[i]
		# return [combo2(arr, n, r, index+1, data, i+1), combo2(arr, n, r, index, data, i+1)]
		combo2(arr, n, r, index+1, data, i+1), combo2(arr, n, r, index, data, i+1)
	return []
def combo(arr, n, r):
	data = range(r)
	combo2(arr, n, r, 0, data, 0)
	# return combo2(arr, n, r, 0, data, 0)
# print combo(range(5), 5, 3)

def combo3(arr, n, r, index, data, i, outputs):
	if index == r:
		return outputs + [data]
	elif i < n:
		data[index] = arr[i]
		o = combo3(arr, n, r, index+1, data, i+1, outputs)
		return combo3(arr, n, r, index, data, i+1, o)
	return outputs

# print combo3(range(5), 5, 3, 0, range(5), 0, [])

def f(x):
	def g(y):
		return 5
	return x+g(x)
# print f(0), f(1)
"""