######To use this code, change '.txt' to '.sage' and  use "sage PGFAN.sage" command.
import sys


maple.read("#/path to PGGW.txt/# PGGW.txt")
import re,time
R.< x,y,z,u,a,b,c > = PolynomialRing(QQ)
fixedOrder = maple.plex(x,y,z,u) 
parOrder = maple.plex(a,b,c)
variables = [x,y,z,u] 
par = [a,b,c]
coneVariables = []
index = 0
#************************************* unll lists
#*************************************
branchList = []
PGfanList = [] 
allBranches = []



#************************************* classes
#*************************************
class GB:
	def __init__(self, N, W, G, facetL, ineqL, vectL, id ):
		self.N = N
		self.W = W
		self.G = G
		self.facetList = facetL
		self.ineqList = ineqL
		self.vectorsList = vectL
		self.id = id

class chain:
	def __init__(self, flag, id):
		self.flag = flag
		self.id = id


#************************************* functions
#*************************************
def Gnew(G):
	
	N = []
	W = []
	gg = []

	for j in range(1,maple.nops(G[2])+1):
				g=G[2][j]
				list=[]
				list.append(g[1])
				list.append(g[2])
				gg.append(list)

	
	for l in range(1,maple.nops(G[0])+1):
		N.append(G[0][l])
	
	for m in range(1,maple.nops(G[1])+1):
		W.append(G[1][m])
	
	return [N,W,gg]
##############
id=identity_matrix(len(variables))
identity=[]

for i in range(0,len(variables)):
	identity.append(list(id[i]))

def identity_vector(v): 
	for i in range(0,len(identity)):
		if identity[i]==v:
			return true
	return false
##############
def search(vectors,v):
	
	for i in range(0,len(vectors)):
		a = vectors[i][0]-vectors[i][1]
		b=[]
		#print maple.igcd(maple.op(a))
		gc = maple.igcd(maple.op(a))
		for j in range(1,maple.nops(a)+1):
			b.append(a[j]/gc)
		
		vv=[]
		for j in range(1,maple.nops(v)+1):
			vv.append(v[j])
		if b==vv:
			return(i)
##############
def GroebnerCone(G,n):
	
	p = MixedIntegerLinearProgram()
	v = p.new_variable(real=True, nonnegative=True)	
	if index == 0 :
		for i in range(0,len(n)):
			coneVariables.append(v[i])
	#print "conev",coneVariables
	p.set_objective(None)
	vec=[]
	ineq=[]			
	facet=[]
	vect=[]
	vv=list(zero_vector(len(n)))
	
	for j in range(0,len(G)):
		f=G[j][1]
		lt=G[j][0]
		p1 = maple.support(f,n)
		p2 = maple.support(lt,n)
		
		
		
		

 		for i in range(1,maple.nops(p1)+1):
 				
				q = p2[1]-p1[i]
				flag=False
				for j in range(1,maple.nops(n)+1):
					
					if q[j]<> 0:
						flag=true


				
				
				if flag==true:
					vec.append([p2[1],p1[i]])
					m = 0
					for i in range(1,maple.nops(n)+1):
						
						m = m+q[i]*v[i-1]
					p.add_constraint(m >= 0)
			
	po = p.polyhedron()	
	ine = po.inequalities_list()
	m = 0
	for i in range(0,len(n)):
		m = m+v[i]

	p.add_constraint(m == 1)
	po1 = p.polyhedron()
	
	
	
	for i in range(0,len(ine)):
		Mn = 0;
		
		ii= maple.elim(ine[i],1)
		w=[]
		for j in range(1,len(n)+1):
			w.append(ii[j])
		
		if identity_vector(w)==False:
			k=search(vec,ii)
			
			vect.append(vec[k])
			facet.append(ii)
		
		
		
		for j in range(0,len(n)):
			Mn = Mn+ii[j+1]*v[j]
		
		ineq.append(Mn>=0)
	return(facet,ineq,vect)
##############
def CheckFinalResult(list):
	for i in range(0,len(list)):
		
		if list[i].flag==0 :
			return 'false'
	return 'true'
##############
def find_interiorPoint(facet,f):
	
	p = MixedIntegerLinearProgram()
	v = p.new_variable(real=True, nonnegative=True)	
	p.set_objective(None)
	
	for i in range(0,len(facet)):
		Mn = 0
		
		for j in range(1,maple.nops(variables)+1):			
			Mn = Mn+facet[i][j]*v[j];
	    	if facet[i]==f:
	    		p.add_constraint(Mn == 0);
	    	else:
				p.add_constraint(Mn >= 0);
			
				
	m = 0
	for i in range(1,maple.nops(variables)+1):
		m = m+v[i]

	p.add_constraint(m == 1)
	po = p.polyhedron()	
	vv=po.center()
	point=[]
	for i in range(0,len(vv)):
		point.append(vv[i]);
	return(point)
##############
def Valid(inequality,facet,v,var):
	
	negative = []
	for j in range(1,maple.nops(facet)+1):
		negative.append(-1*facet[j])

	for i in range(0,len(inequality)):
		idd = inequality[i].id
		
		ineq = allBranches[idd].ineqList
		f = allBranches[idd].facetList
		
		
		
		if str(maple.IsInSystem(ineq,v,var))=='true':
			
			
			flag ='false'
			for j in range(0,len(f)):
				if str(maple.vec(negative,f[j]))=='true':
					flag='true'
			
			if flag=='true':
				return 'false'

	return 'true'

##############
def checkValid(localI,facet,v,var):
	
	for i in range(0,len(localI)):
		a = localI[i]
		if Valid(a , facet , v , var)=='true':
			return 'true'
			
	return 'false'
##############
def AddAndUpdate(G,local):
	
	newlocal = []
	
	for i in range(0,len(G)):
		idd = G[i].id
		
		for j in range(0,len(local)):
			a = local[j]
			b = []
			for k in range(0,len(a)):
				b.append(a[k])
			
			#if str(maple.MemberList(G[i],aa))=='false' :

			ch = chain(0,idd)
			b.append(ch)
			
			newlocal.append(b)
			

	
	return newlocal
##############	
def MergeLists(l):
	N = []
	W = []
	for i in range(0, len(l)):
		idd = l[i].id
		N += allBranches[idd].N
		W += allBranches[idd].W
	can = maple.ccc(N,W,maple.plex(maple.op(par)))
	#print "can",can
  	if str(can[1]) == 'true':
  		return 'true'
  	else:
  		return 'false'
  
   
##############
def convert(list):
	global ALL_Lists,GroebnerList,InequalityList,index
	print "START convert ************************************************************************"
	#print "inputList"
	#print list
	for i in range(0,len(list)):
		print allBranches[list[i].id].N , allBranches[list[i].id].W, allBranches[list[i].id].G
  	
	r = []
	if CheckFinalResult(list) == 'true':

		print "Done! :))))))))))))))))))))))))))))))))))))))))))))))))))))))) "
		Result = []
		for i in range(0,len(list)):
			idd = list[i].id
			Result.append([allBranches[idd].N , allBranches[idd].W, allBranches[idd].G  ])
		PGfanList.append(Result)
		



	else:
		
		
		newResult = []
		input = []
		for i in range(0,len(list)):
		    ch = chain(1,list[i].id)
		    input.append(ch)
			
		newResult.append(input)
		
		
		
		
		for i in range(0,len(list)):
			
			if list[i].flag==0:
				idd = list[i].id 
				G =[allBranches[idd].N , allBranches[idd].W ,allBranches[idd].G ]
				print "G"
				print G
				
				
				
				facet=allBranches[idd].facetList
				print "facets"
				print facet
				
				vectors=allBranches[idd].vectorsList
				
				
				for j in range(0,len(facet)):
					

					v=find_interiorPoint(facet,facet[j])
					
					
					valid = checkValid(newResult,facet[j],v,coneVariables)
					
					
					if (valid =='true') and (str(maple.isValid(vectors[j][0],vectors[j][1],coneVariables,fixedOrder))=='true'):
						print "isvalid for",facet[j],"true"
						
						GNew = maple.flip(G,v ,facet[j], variables,par)
 						print "GNew'"
 						print GNew
						
						

						gForApdate = []
						for k in range(1,maple.nops(GNew)+1):

							h = GNew[k]
							

							N = []
							W = []
							gb = []

							for l in range(1,maple.nops(h[3])+1):
								g=h[3][l]
								li=[]
								li.append(g[1])
								li.append(g[2])
								gb.append(li)

							
							for l in range(1,maple.nops(h[1])+1):
								N.append(h[1][l])
	
							for m in range(1,maple.nops(h[2])+1):
								W.append(h[2][m])
							q = [N,W,gb]
							cone = GroebnerCone(gb , variables)
							new = GB(N,W,gb,cone[0],cone[1], cone[2],index)
							allBranches.append(new)
							
							gForApdate.append(new)
							index += 1
							
							
							
							
							
							
							
						
						newResult = AddAndUpdate(gForApdate,newResult )
						
                        
						
						
						
		print "new Result"
		for i in range(0,len(newResult)):
			bb = newResult[i]
			print "list",i
			for j in range(0,len(bb)):
				print allBranches[bb[j].id].N , allBranches[bb[j].id].W, allBranches[bb[j].id].G
			if MergeLists(bb)=='true':
				branchList.append(bb)


		
		
		print "finish*******************************************************************"


#************************************* main
#*************************************


ideal = [a*x + b*y*x, c*u*y  + z -u*2]
firstGS = maple.PGBMAIN2([],[1],ideal, parOrder , fixedOrder )

for i in range(1,maple.nops(firstGS)+1):
	consist = maple.canspec(firstGS[i][1],firstGS[i][2],maple.plex(maple.op(par)))
	if str(consist[1])=='true':
		if str(consist[3])==str([]):
			c3 = firstGS[i][2]
		else:
			c3 = consist[3]
		
		sim = maple.nonZeroTerms(firstGS[i][3],variables,par,consist[2])
		Rh = maple.reducing(sim,fixedOrder)
		Gmark = maple.mark(Rh,fixedOrder)
		G = Gnew([consist[2],c3,Gmark])
		print G
        Gcone = GroebnerCone(G[2],variables)
        gb = GB(G[0], G[1], G[2], Gcone[0], Gcone[1], Gcone[2], index)
       
        allBranches.append(gb)
     
        ch = chain(0,index)

        branchList.append([ch])
        index+=1
       


while branchList != [] :
	l = branchList.pop()
	convert(l)


#********************************************
#********************************************

file1 = open('finalResult','w')
for i in PGfanList:
	file1.write(str(i))
	file1.write('\n')
	file1.write('\n')
	file1.write("...........................")
	file1.write('\n')
file1.write('\n')
file1.write('***********************************************************')
file1.write('***********************************************************')
file1.write('\n')

Simplifyoutput = maple.SimplifyOutput(PGfanList,par)
for i in range(1,maple.nops(Simplifyoutput)+1):
	file1.write(str(Simplifyoutput[i]))
	file1.write('\n')
	file1.write('\n')
	file1.write("...........................")
	file1.write('\n')

#************************************** correctness test
#**************************************
from random import randint
def generatePoint(a,b,num,par):
	
	ind = len(par)
	re = []
	for i in range(0,num):
		l = []
		
		for j in range(0,ind):
			c = randint(a,b)
			l.append(c)
   

		re.append(l)

		
	return re
	
################
f4 = open('EqualityOfGFan','w')
def EqualityOfGFan(myResult,sageResult,rd):
	print "............................."
	print "start EqualityOfGFan"
	if maple.nops(myResult) != len(sageResult):
		f4.write('Myresult')
		f4.write('\n')
		f4.write(str(myResult))
		f4.write('\n')
		f4.write('......................................................')
		f4.write('\n')
		f4.write('sage')
		f4.write('\n')
		f4.write(str(sageResult))
		f4.write('\n')
		f4.write('********************************************************')
        f4.write('\n') 



	for i in range(1,maple.nops(myResult)+1):
		flag = 'false'
		ele = maple.getIdeal(myResult[i],rd,par)
		L=[]
		for k in range(1,maple.nops(ele)+1):
			L.append(ele[k])
		t1=maple.sim(L,variables)
		for j in range(0,len(sageResult)):
			t2=maple.sim(sageResult[j],variables);
			
			if t1==t2:
				flag = 'true'
				sageResult.pop(j)
				break
		if flag == 'false':
			
			temp=[]
			for k in range(i,maple.nops(myResult)+1):
				temp.append(myResult[k])
			

	
	if sageResult==[]:
		return 'true'
	else:
		return 'false'
	print "finish.........................."	
################
rd = generatePoint(-10,10,30,par)
for i in range(0,len(rd)):
	print i 
	if str(maple.CheckIdeal(ideal,rd[i],par,fixedOrder))=='true':
		ind= maple.checkPoint(Simplifyoutput,rd[i],par)
		if str(ind)=='false':
			print "false",rd[i]
		else:
			print "true"

file2 = open('correctnessTest','w')
RD = maple.getPoint2(-10,10,20,par)
for i in range(1,maple.nops(RD)+1):
	print i
	print RD[i]
	

	if str(maple.CheckIdeal(ideal,RD[i],par,fixedOrder))=='true':
		file2.write("************************** start test")
		file2.write(str(i))
		file2.write("**************************")
		file2.write('\n')
		
	
		ind= maple.checkPoint(Simplifyoutput,RD[i],par)
	
		if str(ind)=='false':
			file2.write("parameter condition false in point   ")
			file2.write(str(RD[i]))
			file2.write('\n')
			
			
			
		else:
		
		
			id=maple.getIdeal(ideal,RD[i],par)
			L=[]
			for j in range(1,maple.nops(id)+1):
				L.append(id[j])
			I = R.ideal(L)
			GF = I.groebner_fan()
			sageR=GF.reduced_groebner_bases()

			e = EqualityOfGFan(Simplifyoutput[ind][3],sageR,RD[i])
			if e=='false':
			
				
				file2.write('\n')
				file2.write('\n')
				file2.write("information:")
				file2.write('\n')
				file2.write("we are testing in point   ")
				file2.write(str(RD[i]))
				file2.write("which is in myResult    ")
				file2.write(str(Simplifyoutput[ind]))
				file2.write('\n')
				file2.write("..........................")
				file2.write('\n')
				file2.write("sageResult:   ")
				file2.write('\n')
				file2.write(str(sageR))
				file2.write('\n')
				file2.write("..........................")
				file2.write('\n')
			
				
		file2.write("************************** finish test ****************************************")
		file2.write('\n')