reset()
from sage.structure.parent import Parent
from sage.structure.element import Element
from sage.structure.unique_representation import UniqueRepresentation
from sage.categories.rings import Rings
from sage.structure.richcmp import rich_to_bool, op_EQ, op_NE, op_LT, op_LE, op_GT, op_GE

#####################################################################
class Apln(UniqueRepresentation):
    
  def __init__(self,n):
      self.n = n
      ts = ['t'+str(t) for t in range(n+1)]
      ys = ['y'+str(t) for t in range(n+1)]
      degrees = tuple((n+1)*[(0,2)])+tuple((n+1)*[(1,2)])
      self.gca = GradedCommutativeAlgebra(QQ, names=ts+ys, degrees=degrees)
      t_relation = -1
      y_relation = 0
      #t_relation = sum(self.gca.gens()[:n])-1
      #y_relation = sum(self.gca.gens()[n:])
      for i in range(n+1):
          t_relation += self.gca.gens()[i]
          y_relation += self.gca.gens()[n+1+i]
      self.gca_with_relations = self.gca.quotient(self.gca.ideal(t_relation,y_relation),check=False)
      dfs = {self.gca.gens_dict()[ts[i]]:self.gca.gens_dict()[ys[i]] for i in range(n+1)}
      self.cdga = self.gca.cdg_algebra(dfs)
  
  def _repr_(self):
    return "Algebra Apln of dimension %s" % self.n
    
  def inject_variables(self):
      self.cdga.inject_variables()
      
  def generators(self):
      return self.cdga.gens()
      
  def reduce_with_relations(self,pol):
      return self.gca_with_relations(pol)
      
  def homogeneousPart(self,pol,degree):
      parts = [pol.homogeneous_parts()[k] for k in pol.homogeneous_parts() if list(k)[0]==degree]
      hompart = 0
      for p in parts:
          hompart += p
      return hompart

  def face(self,i):
      if(self.n==0):
          return 0
      domain = self.cdga
      codomain = self.cdga
      tImages = []
      yImages = []
      n = len(domain.gens())//2
      for k in range(n):
          if(k<i):
              tImages.append(codomain.gens()[k])
              yImages.append(codomain.gens()[n+k])
          elif(k==i):
              tImages.append(0)
              yImages.append(0)
          elif(k>i):
              tImages.append(codomain.gens()[k-1])
              yImages.append(codomain.gens()[n+k-1])
      images = tImages+yImages
      return domain.hom(codomain=codomain,im_gens=images)
      
  def _reduce_polynomial(self,p,n):
    ## Input:
    #  - p: polynomial in variables t_0,..t_n,y_0,...,y_n
    #  - n: dimension of the algebra (apl)_n where p is defined
    ## Output:
    #  - The polynomial after applying the corresponding relations to reduce the variables t_n, y_n
    #  - For practical reasons, we change the nomenclature in this way: t_i ==> xs_i, y_i ==> zs_i
    
    xs = ['x'+str(x) for x in range(n+1)]
    zs = ['z'+str(z) for z in range(n+1)]
    R = PolynomialRing(QQ,names=xs+zs)
    r_gens = R.gens()
    gca_gens = self.gca.gens()
    pol = R(0)
    for monomial_exponents in list(p.dict().keys()):
      monomial = R(1)
      coef = p.dict()[monomial_exponents]
      for i in range(n+1):
        monomial *= r_gens[i]**monomial_exponents[i]
        monomial *= r_gens[n+1+i]**monomial_exponents[self.n+1+i]
      pol += coef*monomial
    xn_relation = 1-sum(r_gens[:n])
    zn_relation = -sum(r_gens[n+1:-1])
    pol = pol.subs({r_gens[n]:xn_relation,r_gens[2*n+1]:zn_relation})
    return pol
    
  def _convert_integrable_polynomial(self,p,n):
    ## Input:
    #  - p: polynomial in variables t_0,..t_n,y_0,...,y_n
    #  - n: dimension of the algebra (apl)_n where p is defined
    ## Output:
    #  - A polynomial in the ring Q<xs_{0},..,xs_{n-1}>.
    #    (After reducing the variables t_{n},y_{n} by the corresponding relations, 
    #     the algorithm removes the variables y_{0},..,y_{n-1}). 

    xs = ['x'+str(x) for x in range(n)]
    zs = ['z'+str(z) for z in range(n)]
    R = PolynomialRing(QQ,names=xs+zs)
    pol = R(self._reduce_polynomial(p,n))
    r_gens = R.gens()
    relations = {}
    for i in range(n):
      relations[r_gens[n+i]]=1
    return pol.subs(relations)
        
  def integral(self,p,n):
    ## Input:
    #  - p: polynomial in variables t_0,..t_n,y_0,...,y_n
    #  - n: dimension of the algebra (apl)_n where p is defined
    ## Output:
    #  - the integral of the polynomial after reducing the variables t_n, y_n by the corresponding relations 
    #
    xs = ['x'+str(x) for x in range(n)]
    R = PolynomialRing(QQ,names=xs)
    r_gens = R.gens()
    pol = R(self._convert_integrable_polynomial(p,n))
    result = pol
    for i in range(n):
      result = result.integral(r_gens[n-1-i])
      result = result.subs({r_gens[n-1-i]:1-sum(r_gens[:n-1-i])})
    return R(((-1)**(n*(n-1)/2)))*result
            
#####################################################################
class AplK(UniqueRepresentation, Parent):
    
  def __init__(self,K, category=None):
    self.K = K
    self.n = K.dimension()
    Parent.__init__(self,category=category )
   
  def _repr_(self):
    return "Algebra AplK over %s" % self.K
    
  def _element_constructor_(self,values):
    return self.element_class(self,values)
    
  def dimension(self):
    return self.K.dimension()
    
  def simplicialComplex(self):
    return self.K
    
  def _n_cells_sorted(self,n):
    return self.K._n_cells_sorted(n)
    
#####################################################################
class AplKElement(Element):

  def __init__(self, parent,values):
    self.elements = values
    self.apl = Apln(parent.dimension())
    Element.__init__(self,parent)
    self.check()
    
  def _repr_(self):
    return str(self.elements)
    
  def dimension(self):
    return self.parent().dimension()
    
  def nullElement(self):
    C = self.__class__
    values = {}
    for k in self.elements.keys():
      values[k]=0
    return C(self.parent(),values)    
    
  def _add_(self,other):
    C = self.__class__
    if(self.elements.keys()!=other.elements.keys()):
      raise ValueError("error: keys do not match")
    newValues = {}
    for k in self.elements.keys():
      newValues[k] = self.elements[k]+other.elements[k]
    return C(self.parent(),newValues)
      
  def _sub_(self,other):
    C = self.__class__
    if(self.elements.keys()!=other.elements.keys()):
      raise ValueError("error: keys do not match")
    newValues = {}
    for k in self.elements.keys():
      newValues[k] = self.elements[k]-other.elements[k]
    return C(self.parent(),newValues)
      
  def _mul_(self,other):
    C = self.__class__
    if(self.elements.keys()!=other.elements.keys()):
      raise ValueError("error: keys do not match")
    newValues = {}
    for k in self.elements.keys():
      newValues[k] = self.elements[k]*other.elements[k]
    return C(self.parent(),newValues)
    
  def _richcmp_(self, other, op):
    if(op==op_EQ):
      if(self.elements.keys()!=other.elements.keys()):
        return False
      for k in self.elements.keys():
        if(self.elements[k]!=other.elements[k]):
          return False
      return True
    elif(op==op_NE):
      return not self._richcmp_(other,op_EQ)
    else:
      raise NotImplementedError("operation not implemented for this element class")
      
  def differential(self):
    C = self.__class__
    newValues = {}
    for k in self.elements.keys():
      newValues[k] = self.apl.cdga(self.elements[k]).differential()
    return C(self.parent(),newValues)
      
  def check(self):
    S = self.parent().simplicialComplex()
    dim = S.dimension()
    for n in range(1,dim+1):
      print("Checking dimension " + str(n))
      for simplex in S._n_cells_sorted(n):
        for f in range(n+1):
          sFace = tuple(S.face(simplex,f))
          aFace = self.apl.face(f)(self.elements[tuple(simplex)])
          if(self.elements[sFace] != aFace):
            raise ValueError("the assignation is not valid")
            
  def homogeneous_component(self,degree):
    C = self.__class__
    component = {}
    for k in self.elements.keys():
      component[k] = self.apl.homogeneousPart(self.elements[k],degree)
    return C(self.parent(),component)
    
  def homogeneous_components(self):
    components = {}
    for n in range(self.dimension()+1):
      components[n] = self.homogeneous_component(n)
    return components
               
  def is_homogeneous(self):
    isHomogeneous = True
    for k in self.elements.keys():
      isHomogeneous = isHomogeneous and self.apl.cdga(self.elements[k]).is_homogeneous()
    return isHomogeneous
  
  def is_in_kernel(self):
    diff = self.differential()
    if(diff == self.nullElement()):
      return True
    else:
      return False
      
  def _integrable_component(self):
    C = self.__class__
    component = {}
    for k in self.elements.keys():
      degree = len(k)-1
      component[k] = self.apl.homogeneousPart(self.elements[k],degree)
      #component[k] = self.homogeneous_component(dim).elements[k]
    return component
    
  def integral(self):
    component = self._integrable_component()
    cochain = {}
    #integration of each polynomial
    for k in component.keys():
      print("Integrating simplex " + str(k))
      n = len(k)
      pol = self.apl.cdga(component[k])
      cochain[k] = self.apl.integral(pol,n)
    return cochain

  
      
#####################################################################
#S = SimplicialComplex([[0,1,2]],is_mutable=False)  
S = SimplicialComplex([[0,1],[0,2],[1,2]],is_mutable=False)                                                                                 
apl = Apln(S.dimension())
apl.inject_variables()

params1 = {}
params1[(0,)] = t0
params1[(1,)] = t0
params1[(2,)] = 3*t0
params1[(0,1)] = t0+t1
params1[(0,2)] = t0+3*t1
params1[(1,2)] = t0+3*t1

params2 = {}
params2[(0,)] = 3*t0
params2[(1,)] = 4*t0
params2[(2,)] = 5*t0
params2[(0,1)] = 3*t0+4*t1
params2[(0,2)] = 3*t0+5*t1
params2[(1,2)] = 4*t0+5*t1

params3 = {}
for d in range(S.dimension()+1):
  for c in S._n_cells_sorted(d):
    params3[tuple(c)]=randint(1,100)*t0
    
params4 = {}
params4[(0,)] = 3*t0+y0
params4[(1,)] = 4*t0
params4[(2,)] = 5*t0
params4[(0,1)] = 3*t0+4*t1+y0
params4[(0,2)] = 3*t0+5*t1+y0
params4[(1,2)] = 4*t0+5*t1

A = AplK(S)
e1 = AplKElement(A,params1)
e2 = AplKElement(A,params2)
#e3 = AplKElement(A,params3)
e4 = AplKElement(A,params4)


class morphism_from_QQ_to_AplK(Morphism):

  def __init__(self, codomain):
    self.simplicialComplex= codomain.K
    Morphism.__init__(self, QQ, codomain)

  def _call_(self, x):
    S = self.simplicialComplex
    params1 = {}
    for d in range(S.dimension()+1):
      for c in S._n_cells_sorted(d):
        params1[tuple(c)]=x
    A = AplK(S)
    return AplKElement(A,params1)
    
f = morphism_from_QQ_to_AplK(A)
A.register_coercion(f)

#p = t0*t1**2*y0+t0*y1


S = SimplicialComplex([[0,1,2,3,4,5,6]],is_mutable=False)                                                                                 
apl = Apln(S.dimension())
apl.inject_variables()
p = y0*y1*y2*y3*y4*y5+y0*y1*y2*y3*y6+t4*t5

initial_assignation = {}
initial_assignation[(0,1,2,3,4,5,6)] = p

def generateElement(S,initial_assignation):
  apl = Apln(S.dimension())
  aplk_structure = initial_assignation.copy()
  
  for n in reversed(range(1,S.dimension()+1)):
    print("Generating dimension "+ str(n))
    for simplex in S._n_cells_sorted(n):
      p = aplk_structure[tuple(simplex)]
      for f in range(n+1):
        sFace = tuple(S.face(simplex,f))
        aFace = apl.face(f)(p)
        if(sFace not in aplk_structure.keys()):
          aplk_structure[tuple(sFace)] = aFace
        elif(aplk_structure[tuple(sFace)] != aFace):
          raise ValueError("Error during the checking. The assignation is not valid")
  A = AplK(S)
  return AplKElement(A,aplk_structure)
          
          
g = generateElement(S,initial_assignation)
#g.check()
#cochain = g.integral()
  
#n = 1
#a = Apln(n)
#a.inject_variables()
#p = t0*y1
##print(p)
#print(a._reduce_polynomial(p,n))
#print(a._convert_integrable_polynomial(p,n))
#print(a.integral(p,n))