r"""
<Short one-line summary that ends with no period>

<Paragraph description>

EXAMPLES::

<Lots and lots of examples>

AUTHORS:

- YOUR NAME (2005-01-03): initial version

- person (date in ISO year-month-day format): short desc

"""

# ****************************************************************************
#       Copyright (C) 2021 YOUR NAME <your email>
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 2 of the License, or
# (at your option) any later version.
#                  https://www.gnu.org/licenses/
# ****************************************************************************
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
from sage.categories.morphism import Morphism
from sage.rings.rational_field import QQ
from sage.algebras.commutative_dga import GradedCommutativeAlgebra
from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
from sage.functions.other import factorial

#####################################################################
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)
        #self.inject_variables()
  
    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
        
    def integral2(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 = 0
        for k in pol.dict().keys(): # k is a monomial
            numerator = 1
            denominator = n
            for e in k: # e is an exponents
                numerator *= factorial(e)
                denominator += e
            coef = pol.dict()[k]
            result += coef*numerator/factorial(denominator)
        return R(((-1)**(n*(n-1)/2)))*result
            
#####################################################################

class MorphismFromQQToAplK(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)
    
#####################################################################

class AplK(UniqueRepresentation, Parent):
    
    def __init__(self,K, category=None):
        self.K = K
        self.n = K.dimension()
        Parent.__init__(self,category=category )
        f = MorphismFromQQToAplK(self)
        self.register_coercion(f)
   
    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):
        r"""
        Return sum of the elements

        INPUT:

        - ``other`` -- An AplKEelement


        OUTPUT: An AplKElement equal to the sum of self and the element passed as argument

        EXAMPLES:

        This example illustrates ... ::

            sage: S = SimplicialComplex([[0,1],[0,2],[1,2]],is_mutable=False) 
            sage: apl = Apln(S.dimension())
            sage: (t0, t1, y0, y1) = apl.generators()  
            sage: params1 = {(0,): t0, (1,): t0, (2,):3*t0, (0,1): t0+t1,(0,2): t0+3*t1, (1,2): t0+3*t1}
            sage: params2 = {(0,): 3*t0, (1,): 4*t0, (2,):5*t0, (0,1): 3*t0+4*t1,(0,2): 3*t0+5*t1, (1,2): 4*t0+5*t1}
            sage: A = AplK(S)
            sage: e1 = AplKElement(A,params1)
            sage: e2 = AplKElement(A,params2)  
            sage: e1+e2
            {(0,): 4*t0, (1,): 5*t0, (2,): 8*t0, (0, 1): 4*t0 + 5*t1, (0, 2): 4*t0 + 8*t1, (1, 2): 5*t0 + 8*t1}
          
        """
        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)-1
            pol = self.apl.cdga(component[k])
            cochain[k] = self.apl.integral(pol,n)
        return cochain