#Este codigo lleva 20 dias corriendo en cluster y me pregunto cuando terminara?? o esta en un bucle infinito???
from sympy import *
from sympy.physics.quantum import TensorProduct
u=symbols('u'); v=u-1

umatrix=Matrix([
   [u,0,0,0],
   [0,u,0,0],
   [0,0,1,0],
   [0,0,0,1]
])

rmatrix=Matrix([
   [u,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
   [0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0],
   [0,0,0,0,0,0,0,0,u,0,0,0,0,0,0,0],
   [0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0],
   [0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
   [0,0,0,0,0,u,0,0,0,0,0,0,0,0,0,0],
   [0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0],
   [0,0,0,0,0,0,0,0,0,0,0,0,0,u,0,0],
   [0,0,1,0,0,0,0,0,v,0,0,0,0,0,0,0],
   [0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0],
   [0,0,0,0,0,0,0,0,0,0,u,0,0,0,0,0],
   [0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0],
   [0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0],
   [0,0,0,0,0,0,0,1,0,0,0,0,0,v,0,0],
   [0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0],
   [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,u]
])

def num_strands(list_braid):
    """ computes the size of the braid """
    res=abs(list_braid[0])+1
    for elem in list_braid[1:]:
        if abs(elem)>res-1: res=abs(elem)+1
    return res

def tensor_identities(num):
    res={1:Matrix([
   [1,0,0,0],
   [0,1,0,0],
   [0,0,1,0],
   [0,0,0,1]
    ])}
    for ind in range(2,num+1): res[ind]=eye(4**(ind))
    return res

def tensor_phi(strands,tensor_id_dic,indexp):
    res=rmatrix**(indexp/abs(indexp))
    if abs(indexp)-1: res=TensorProduct(tensor_id_dic[abs(indexp)-1],res)
    if strands-abs(indexp)-1: res=TensorProduct(res,tensor_id_dic[strands-abs(indexp)-1])
    return res

def tensor_phies(strands,tensor_id_dic,list_braid):
    res=dict()
    for indexp in list_braid:
        print(indexp)
        if indexp not in res: res[indexp]=tensor_phi(strands,tensor_id_dic,indexp)
    return res

def mu_matrix(num):
    res=umatrix
    for index in range(num-1): res=TensorProduct(res,umatrix)
    return res

def braid_knot_inv(list_braid):
    strands=num_strands(list_braid)
    tensor_id_dic=tensor_identities(strands-2)
    tensor_ph_dic=tensor_phies(strands,tensor_id_dic,list_braid)
    phi_braid=tensor_ph_dic[list_braid[0]]
    for elem in list_braid[1:]:
        print(7)
        phi_braid=phi_braid*tensor_ph_dic[elem]  
    final_product=phi_braid*mu_matrix(strands)
    wrt=0
    for elem in list_braid: wrt=wrt+(elem/abs(elem))
    print('u^('+str((1-(3*wrt)-strands)/2)+')*('+str(final_product.trace())+')/('+str(2*(1+u))+')')
    #print('u^('+str(int((1-(3*wrt)-strands)/2))+')*('+str(final_product.trace())+')/('+str(2*(1+u))+')')

#braid_knot_inv([1])
#braid_knot_inv([-1, -2, 3, -2, -3, 2, -1, -3, -3, 2, -3])
braid_knot_inv([1,-2,-3,-4,3,3,-5,4,-3,2,-1,-3,-2,-4,3,-2,-2,-2,5,4,-3])
#braid_knot_inv([1,-2])