###############################################################################
##   Characteristic cone given by equation module 						     ##
##                                                                           ##
##   Author: Ayan Sengupta, Mousumi Mukherjee and Debasattam Pal             ##
##   Version: 1.0                                                            ##
###############################################################################


reset()
import re
import sys


#code edited on 17 Jan 2017
# Written in SageMath version 7.4 and Python 2.7

''' 
This program takes the generators of the cone and the partial difference equations of discrete 
autonomous nD system as inputs and gives a 'yes/no' certificate according to whether the given cone 
is a characteristic cone for the given system or not.

First the dimension of the system, a kernel representation matrix and the number of cone 
generating vectors are taken as input. Then the matrix and the vectors are given by
the user.

As output the program will give if the given cone is a characteristic cone or not.

The dimension of the system determines the number of free variables of the system. According to the
input given by the user, the free variables are defined in x, as x1,x2,x3... . So the variables 
that are valid during matrix input are x1, x1^-1, x2, x2^-1, x3, x3^-1 ... . To make the input 
more readable, and reduce typo errors, user can use y variables instead of x inverses i.e. 
user can give input in form of x1, y1, x2, y2, x3, y3 ..., where y1, y2, y3... are considered as
x1^-1, x2^-1, x3^-1...

Transposed column vectors are taken as input for generating the cone. For example
 -c1 = 0 1 0 is actually the vector [0 1 0]^T 

A few examples are shown below: 

    EXAMPLE 1:

   - Dimension of the system: 3
   - Dimension of the kernel representation matrix: (nxm): 3x2
   - Number of vectors for generating the cone: 3

   -Please enter the Kernel Representation matrix:
   -->  x3, x2^3
   -->  x2, y2*y3 -1 
   -->  y1*y3, x1+x2

   -please enter the transposed column vectors to define the cone:
   -c1= 1 0 0
   -c2= 0 1 0
   -c3= 0 0 1

   -This is a characteristic cone    
	
	Now in the next example we will use inverse of x, in places of y.

     EXAMPLE 2:
   
   - Dimension of the system: 3
   - Dimension of the kernel representation matrix: (nxm): 3x2
   - Number of vectors for generating the cone: 3

   -Please enter the Kernel Representation matrix:
   -->  x2, x2^-1*x3^-1 -1
   -->  x3, x2^3 
   -->  x1^-1*x3^-1, x1+x2 

   -please enter the transposed column vectors to define the cone:
   -c1= 1 0 0
   -c2= 0 1 0
   -c3= 0 0 1

   -This is a characteristic cone 

CAUTION:

As the implementation of the code depends on computing the groebner basis of the module, there might
be cases where the groebner basis is too massive. In such cases computation gets too heavy and, SAGE or particularly, SINGULAR 
crashes after some time. Such an example is shown below:

EXAMPLE:

Dimension of the system: 3
Size of the kernel representation matrix: (nxm): 3x2
Number of vectors for generating the cone: 3

Please enter your Kernel Representation matrix:
Example: -> x1, x2, x3 

->  x1^5*x3+x2^6+x3^7*x2, x3*y1 -1
->  x2, y2*y3 -1
->  y1*y3, x1+x2

please enter the transposed column vectors of the cone:
c1= 1 0 0
c2= 0 1 0
c3= 0 0 1

HOW TO RUN THE CODE IN SAGE:

First open SAGE from the terminal. Then you can load a file in a Sage session with the command "load".

example: 

if your file is in the current directory where you started sage then: 
 load("your_file.sage")

if your file is in another directory then :
  load("~/the_file_location/your_file.sage")


'''

def variableNames(num, v_num,matrix_r,matrix_c):
	global x_var,y_var,t_var,r_vector,p_vector,v_vector,rg_module,z_module
	x_var = [var('x%s'  %p) for p in range(1, num + 1)]
	y_var = [var('y%s' %p) for p in range(1, num + 1)]
	t_var = [var('t%s' %p) for p in range(1, v_num + 1)]
	r_vector = [var('r%s' %p) for p in range(1, matrix_r + 1)]
	p_vector = [var('p%s' %p) for p in range(1, num*matrix_c + 1)]
	v_vector = [var('v%s' %p) for p in range(1, v_num*matrix_c + 1)]
	rg_module = [var('rg%s' %p) for p in range(1, num_var*2*matrix_c + 1)]
	z_module = [var('z%s' %p) for p in range(1, num_var*2*matrix_c + 1)]


def polynomialRing():
	

	variables = str(tuple(x_var+y_var+t_var))
	P = singular.ring(0, variables,'(lp,c)')


	return P

def matrixInput(ring,rows,columns):

	print("\nPlease enter your Kernel Representation matrix:")
	print("Example: -> x1, x2, x3 \n")
	

	matrix =[]
	for i in range(rows):
		matrix.append([x for x in raw_input("->  ").split(",")])
		
	for i in range(rows):
		for j in range(columns):
			matrix[i][j]= ''.join( c for c in matrix[i][j] if  c not in '{(})[]' )
			replace_positions = [m.start() for m in re.finditer('\^-', matrix[i][j])]
			
			a = list(matrix[i][j])
			
			if len(replace_positions)>0:
				for p in xrange(len(replace_positions)):       	
					a[replace_positions[p]-2] = 'y'
					a[replace_positions[p]+1] = ''
			a2 = ''.join(a)
			a3 = re.sub('\^','**',a2)
			
			matrix[i][j] = eval(a3)
	return matrix

def coneModule(num,num_var,matrix_c):
	'''This function takes the cone generating vectors as input. The vectors provided must be transposed 
    column vectors.  '''
	print("\nplease enter the transposed column vectors of the cone:")
	
	vector_poly = []
	for i in xrange(num):
		x_vector =1
		y_vector =1
		vector_monomial =1
		print 'c' + str(i + 1) + '=',
		
		temp = raw_input().strip(" ").split(" ")
		vector = [int(num) for num in temp]
		
		if len(vector) != num_var:
			print("Wrong Vector length")
			break
		for p in xrange(len(vector)):

			if vector[p] > 0:
				x_vector = (x_var[p])^vector[p]
				
			if vector[p] < 0:
				y_vector = (y_var[p])^abs(vector[p])	
				
			vector_monomial = x_vector*y_vector*vector_monomial
			x_vector =1
			y_vector =1  
			
		poly_temp = t_var[i] - vector_monomial
		
		for k in xrange(matrix_c):
			list_temp =[0]*matrix_c
			list_temp[k] = poly_temp
			vector_poly.append(list_temp)
	return vector_poly


	
def varModule():
	var_matrix =[]
	all_var =x_var+y_var
	for i in xrange(len(all_var)):
		for k in xrange(matrix_c):
			list_temp =[0]*matrix_c
			list_temp[k] = all_var[i]
			var_matrix.append(list_temp)
	return var_matrix		

def laurentModule(matrix_c):
	laurent_matrix = []
	for i in xrange(num_var):
		laurent_temp = x_var[i]*y_var[i] - 1
		for k in xrange(matrix_c):
			list_temp =[0]*matrix_c
			list_temp[k] = laurent_temp
			laurent_matrix.append(list_temp)
			
	return laurent_matrix		

def checkRemainder(remainder):
	remainder = eval('+'.join(remainder.replace(" ",'').replace("^",'**').split(',\n')))
	try:
		remainder = remainder.variables()
	except:
		return True
	if not set(remainder)<=set(t_var):
		return False   
	else:
		return True


def singularCode(B_matrix,C_matrix,L_matrix,V_matrix):
	for i in xrange(len(r_vector)):
		singular.set('vector',str(r_vector[i]),str(B_matrix[i]) )
	singular.set('module','input',str(r_vector).replace('[','').replace(']',''))

	for i in xrange(len(v_vector)):
		singular.set('vector',str(v_vector[i]),str(C_matrix[i]) )
	singular.set('module','cone',str(v_vector).replace('[','').replace(']',''))

	for i in xrange(len(p_vector)):
		singular.set('vector',str(p_vector[i]),str(L_matrix[i]) )
	singular.set('module','laurent',str(p_vector).replace('[','').replace(']',''))
	
	for i in xrange(len(z_module)):
		singular.set('module',str(z_module[i]),str(V_matrix[i]) )
	

	singular.set('module','added_module','input+cone+laurent')
	singular.set('module','G','groebner(added_module)')

	for i in xrange(num_var*matrix_c*2):
		
		characteristic_cone = 1
		
		
		singular.set('module',str(rg_module[i]),str('reduce('+str(z_module[i])+',G)'))
		remainder_temp = singular.get(str(rg_module[i]))
		if 'rg' in str(remainder_temp):
			for k in range(1,matrix_c+1):
				rem = singular.get(str(rg_module[i])+'['+str(k)+',1]')
				result = checkRemainder(rem)
				if result == False:
					characteristic_cone =0
					break

		else:	
		    result = checkRemainder(remainder_temp)
		    if result == False:
		    	characteristic_cone =0
		    	break
	if characteristic_cone == 0:
		return False
	else:
		return True


    
    	
def output(logic):
	if logic == True:
		print("\nThis is a characteristic cone")
	else:
		print("\nThis is NOT a characteristic cone")



#################### MAIN #######################
test_again = True
num_var = int(input("Dimension of the system: "))
matrix_r, matrix_c = raw_input("Size of the kernel representation matrix: (nxm): ").split('x')
matrix_r = int(matrix_r)
matrix_c = int(matrix_c)
if matrix_r < matrix_c:
	sys.exit("Wrong dimension of kernel representation matrix")

v_num = input("Number of vectors for generating the cone: ")
variableNames(num_var,v_num,matrix_r,matrix_c)
P_ring = polynomialRing() 
input_module = matrixInput(P_ring,matrix_r,matrix_c)
cone_module = coneModule(v_num, num_var, matrix_c)
laurent_module = laurentModule(matrix_c)
var_module = varModule()
singular_setup = singularCode(input_module,cone_module,laurent_module,var_module)
output(singular_setup)
	

while test_again ==True:
	bool_again = raw_input("\nDo you want to test it with another set of cone vectors? (y/n): ")
	if bool_again == 'y' or bool_again =='Y':
		v_num = input("Number of vectors for generating the cone: ")
		cone_module = coneModule(v_num, num_var, matrix_c)
		singular_setup = singularCode(input_module,cone_module,laurent_module,var_module)
		output(singular_setup)


	else:
		test_again = False