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

<Paragraph description>

EXAMPLES::

<Lots and lots of examples>

AUTHORS:

- Carlos Alquézar Baeta (2023-02-07): initial version
- Miguel Ángel Marco (2023-02-07): initial version
- Jorge Martín Morales (2023-02-07): initial version

"""

# ****************************************************************************
#       Copyright (C) 2023 Carlos Alquézar Baeta <alquezar@unizar.es>
#
# 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.algebras.commutative_dga import *
from sage.structure.parent import Parent
from sage.structure.element import Element
from sage.structure.richcmp import rich_to_bool, op_EQ, op_NE, op_LT, op_LE, op_GT, op_GE
from sage.rings.polynomial.polydict import ETuple
from copy import deepcopy
import itertools

class AplN(DifferentialGCAlgebra_multigraded, UniqueRepresentation):
    """
      A differential multi-graded commutative algebra

      This class extends the DifferentialGCAlgebra_multigraded class
      to create the algebra of polynomial differential forms with coefficients in QQ

    INPUT:

    - ``A`` -- a multi-graded commutative algebra
    - ``differential`` -- the differential rules

    """

    def __init__(self, A, differential):
        """
        Initialize ``self``.

        INPUT:

        - ``A`` -- a multi-graded commutative algebra
        - ``differential`` -- a differential
        """

        self.n = len(A._names) // 2 - 1
        DifferentialGCAlgebra_multigraded.__init__(self, A, differential)
        self.t_relation = sum(self.gens()[:self.n+1])-1
        self.y_relation = sum(self.gens()[self.n+1:])

    def _base_repr(self):
        """
        Return the base string representation of ``self``.
        """

        s = DifferentialGCAlgebra._base_repr(self)
        old = '{}'.format(self._degrees)
        new = '{}'.format([i[0] for i in self._degrees_multi])
        return s.replace(old, new)

    def integrand(self, value):
        """
        Returns a polynomial whose definite integral is equal to 'value'

        OUTPUT:
            - A polynomial
        """
        if (value in QQ) or ((value in A) and (value.degree().is_zero())):
            #return value*factorial(self.n)*prod(self.gens()[self.n+2:])
            
            gens = self.gens()
            tExp = [0*randint(0,4) for i in range(self.n+1)] # Random pol in variables t_i, 1<=i<=n
            yExp = [0]+[1]*self.n
            exp =  tuple(tExp + yExp)
            #print(value, " - ", exp, flush=True)
            numerator = factorial(sum(tExp)+self.n)
            denominator = prod([factorial(e) for e in tExp])
            coef = value*QQ(numerator)/QQ(denominator)
            return coef*self.get_term(self, 1, exp, gens)
            
        else:
            raise ValueError("only field elements can be integrated.")
     
    @cached_function       
    def findsubsets(self,S,m):
        return set(itertools.combinations(S,m))    
            
    @cached_function
    def get_term(self, coef, exp, gens):
        """"

        Returns the term specified by its coefficient, exponents and generators

        INPUT::
            - coef: The coefficient of the term
            - exp: A tuple with the exponents of each variable
            - gens: The generators

        OUTPUT::
            - A polynomial
        """
        return coef*prod(list(map(lambda exp, gen: gen**exp, exp, gens)))
            
    @cached_function
    def get_primitive_basis_of_a_monomial(self, monomial):
        """
        Returns a basis of the primitive of a given monomial

        Example:

        """
        n_vars = len(monomial) // 2
        combinations = []
        for i in range(n_vars):
            mono = list(deepcopy(monomial))
            if mono[n_vars+i]:
                mono[i] += 1
                mono[n_vars+i] = 0
                combinations.append(tuple(mono))
        return list(set(combinations))
        
    @cached_function
    def add_kernel_terms_to_primitive_basis(self, tdegree, ydegree):
        N = len(self.gens())//2
        ygens = self.gens()[N+1:] #N+1 because we remove y0 generator
        listProds = [prod(s) for s in self.findsubsets(self, ygens,ydegree)]
        listCoefs = [tuple(2*N*[0])]
        if tdegree>1:
            #V = VectorSpace(GF(tdegree),N-1)
            #listCoefs = [tuple([0]+list(v)+N*[0]) for v in V]
            l = list(range(tdegree+1))
            listCoefs = [tuple([0]+ list(i)) for i in cartesian_product((N-1)*[l]).list()]
        listPolinomials = [self.get_term(self, 1, c, self.gens()) for c in listCoefs]
        cp = cartesian_product([listProds,listPolinomials])
        elements = [prod(c) for c in cp]
#            for e in elements:
#                print(e,flush=True)
        return [tuple(e.dict().keys())[0] for e in elements]

    class Element(GCAlgebra_multigraded.Element, DifferentialGCAlgebra.Element):
        """
        Element class of the commutative differential multi-graded algebra (Apl)n.
        """

        def is_homogeneous(self):
            """
            Return "True" if "self" is homogeneous and "False" otherwise.
            """
            if self.is_zero():
                return True
            degrees = set([k[0] for k in self.homogeneous_parts()])
            return len(degrees) == 1

        def total_degree(self):
            """
            Return the total degree of this element.
            """
            homParts = self.homogeneous_parts()
            if not homParts: # the polynomial is equal to zero
                return self.parent().n
            degrees = set([k[0] for k in homParts])
            return max(degrees)

        def homogeneous_part(self, degree):
            """
            Returns the homogeneous part of the element for a given degree.
            """

            if (degree < 0):
                raise ValueError('`degree` must be a non-negative integer')
            parts = [self.homogeneous_parts()[k]
                     for k in self.homogeneous_parts() if list(k)[0] == degree]
            return self.parent()(sum(parts))

        #@cached_function
        def reduce(self, relations=[]):
            """"
            The method reduces the polynomial by the Ideal <sum(t_i)-1, sum(y_j)>

            INPUT::
                - p: polynomial in variables t_0,..t_n,y_0,...,y_n
                - relations (optional): if relations are given, the polynomial is reduced following
                    this rules. If no relations are given, the polynomial is reduced taking into
                    account the relations of the parent algebra.

            OUTPUT::
                - The polynomial after applying the corresponding relations
            """
            apl = self.parent()
            if relations:
                return DifferentialGCAlgebra.Element.reduce(self, relations)
            else:
                return DifferentialGCAlgebra.Element.reduce(self, [apl.t_relation, apl.y_relation])

        def _integrable_part(self):
            """
            The method returns the coefficient in Q[t_1,...,t_n] of the term
            of degree `n`

            OUTPUT:
                - A polynomial in the ring Q<t_1,..t_n>.
                  (After reducing the variables t_{0},y_{0} by the corresponding relations,
                  the algorithm removes the variables y_{1},..,y_{n}).
            """

            apl = self.parent()
            # needed in case reduce() returns an integer
            pol = apl(self.reduce())
            relations = []
            for i in range(apl.n+1):
                relations.append(apl.gens()[apl.n+1+i]-1)
            pol = apl(pol.homogeneous_part(degree=apl.n))
            return apl(pol.reduce(relations=relations))

        def integral(self):
            """
            Returns the definite integral of the term of degree `n`

            OUTPUT:
                - The rational number obtained by the definite integral
                  of the term of degree `n`

            """

            apl = self.parent()
            pol = self._integrable_part()
            result = 0
            for k in pol.dict().keys():  # k is a monomial
                numerator = 1
                denominator = apl.n
                for e in k:  # e is an exponents
                    numerator *= factorial(e)
                    denominator += e
                coef = pol.dict()[k]
                result += coef*numerator/factorial(denominator)
            # return ((-1) ** ((apl.n*(apl.n-1)) // 2)) * result
            return result

        def _get_terms(self):
            """"

            Returns self seen as a collection of terms in the variables y_j with coefficients in QQ<t_1,..,t_n>

            OUTPUT::
                - A dictionary with all terms

            """
            a = self.parent()
            n = len(a.gens())//2
            ts_only_gens = a.gens()[:n]+tuple(n*[1])
            ys_only_gens = tuple(n*[1])+a.gens()[n:]
            terms = {}
            pol = self.dict()
            for exp in pol.keys():
                key = a.get_term(a,1, exp, ys_only_gens)
                coef = a.get_term(a,pol[exp], exp, ts_only_gens)
                if key not in terms.keys():
                    terms[key] = coef
                else:
                    terms[key] += coef
            return terms

        def _get_leading_pol(self):
            """"

            Returns the leading polynomial

            OUTPUT::
                - A polynomial
            """
            terms = self._get_terms()
            keys = list(terms.keys())
            keys.sort(reverse=True)
            return terms[keys[0]]*keys[0]

        def _get_idx(self, exp):
            """"

            Returns the idx 'i' of the first generator 'y_i'  that appears in exp

            INPUT::
                - exp: A tuple with the exponents

            OUTPUT::
                - A non-negative number lower or equal to 'n'
            """
            n = len(exp)//2
            filter_ts = n*[0]+n*[1]
            diff_forms = list(map(lambda v, f: v*f, exp, filter_ts))
            return next((i for i, x in enumerate(list(diff_forms)) if x), None)-n

        def _monomial_primitive(self, monomial, idx):
            """"

            Returns the primitive of the given monomial

            INPUT::
                - monomial: A monomial
                - idx: the index of the variable y_j that will be used to calculate the pseudo-primitive

            OUTPUT::
                - A polynomial
            """
            exp = list(monomial.dict().keys())[0]
            coef = monomial.dict()[exp]
            n = len(exp)//2
            if (exp[n+idx] == 0):
                raise IndexError(
                    "It is not possible to get a primitive using the given idx.")
            primitive = list(deepcopy(exp))
            primitive[n+idx] = 0
            primitive[idx] += 1
            new_coef = coef/primitive[idx]
            a = monomial.parent()
            return a.get_term(a,new_coef, primitive, a.gens())

        def _leading_primitive(self):
            """"
            Returns the primitive of the leading term of self

            OUTPUT::
                - A polynomial
            """
            polynomial = self._get_leading_pol()
            # necessary condition: all monomials have apl.gens[n+idx] factor
            exp = list(polynomial.dict().keys())[0]
            idx = self._get_idx(exp)
            apl = polynomial.parent()
            aGens = apl.gens()
            primitive = 0
            for exp in polynomial.dict().keys():
                coef = polynomial.dict()[exp]
                monomial = apl.get_term(apl,coef, exp, aGens)
                primitive += self._monomial_primitive(monomial, idx)
            terms = primitive._get_terms()
            result = sum([terms[key]*key for key in terms.keys()])
            return result

        def primitive(self):
            """"
            Retuns a polynomial whose differential is self

            OUTPUT::
                - A polynomial
            """
            if self.differential() != 0:
                raise ValueError("This polynomial is not closed.")
            if not self: # it is the constant 0
                a = self.parent()
                return a(0)
            w = self._leading_primitive()
            theta = self - w.differential()
            if (theta == 0):
                return w
            else:
                return w+theta.primitive()

@cached_function
def AlgebraPL(n):
    """
    The graded commutative differential cochain algebra Apl of dimension ``n``

    INPUT:

    - ``n`` -- A non-negative integer corresponding to the number of generators of degree 0.
    """

    if (n < 0):
        raise ValueError('`n` must be a non-negative integer')
    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)])
    gca = GradedCommutativeAlgebra(QQ, names=ts+ys, degrees=degrees)
    differentials = {gca.gens_dict()[ts[i]]: gca.gens_dict()[
        ys[i]] for i in range(n+1)}
    return AplN(gca, differentials)

#############################################
# TEST
# print("Test AlgebraPL")
#A = AlgebraPL(5)
# A.inject_variables()
#p = t0**2 + t1 + y3 + y1*y2
#q = p.differential()
#############################################


class SimplicialAplN(AplN, UniqueRepresentation):
    """

      The simplicial commutative cochain algebra Apl

    """

    def _base_repr(self):
        """
        Return the base string representation of ``self``.
        """

        s = DifferentialGCAlgebra._base_repr(self)
        old = '{}'.format(self._degrees)
        new = '{}'.format([i[0] for i in self._degrees_multi])
        return "Simplicial " + s.replace(old, new)

    class Element(AplN.Element):
        """
        Element class of the simplicial commutative differential algebra (Apl)n.
        """
        @cached_method
        def face_morphism(self, i):
            """"

            Returns the morphism face 'i' of self

            INPUT::

                - i: A non-negative integer lower or equal to 'n'. Represents the face index

            OUTPUT::

                - A morphism of cochain algebras
            """
            n = self.parent().n
            if(n == 0):
                return 0
            domain = self.parent()
            codomain = SimplicialAlgebraplN(n-1)
            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-1+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-2])
            images = tImages+yImages
            return domain.hom(codomain=codomain, im_gens=images)
        
        @cached_method
        def face(self, i):
            """"

            Returns the face 'i' of self

            INPUT::

                - i: A non-negative integer lower or equal to 'n'. Represents the face index

            OUTPUT::

                - A polynomial that is the face 'i' of self
            """

            return self.face_morphism(i)(self)

        @cached_method
        def degeneration_morphism(self, j):
            """"

            Returns the morphism degeneration 'j' of self

            INPUT::
                - j: A non-negative integer lower or equal to 'n'. Represents the degeneration index

            OUTPUT::
                - A morphism of cochain algebras
            """

            n = self.parent().n
            domain = self.parent()
            codomain = SimplicialAlgebraplN(n+1)
            if(j > n):
                print(
                    "the degeneration index must be a non-negative number lower or equal to "+str(n))
                return False
            elif(j < 0):
                print("the degeneration index must be non-negative")
                return False
            else:
                tImages = []
                yImages = []
                n = len(domain.gens())//2
                for k in range(n):
                    if(k < j):
                        tImages.append(codomain.gens()[k])
                        yImages.append(codomain.gens()[n+1+k])
                    elif (k == j):
                        tImages.append(codomain.gens()[k]+codomain.gens()[k+1])
                        yImages.append(codomain.gens()[
                                       n+1+k]+codomain.gens()[n+1+k+1])
                    elif (k > j):
                        tImages.append(codomain.gens()[k+1])
                        yImages.append(codomain.gens()[n+1+k+1])
                images = tImages+yImages
                return domain.hom(codomain=codomain, im_gens=images)

        def degeneration(self, j):
            """"

            Returns the degeneration 'j' of self

            INPUT::
                - j: A non-negative integer lower or equal to 'n'. Represents the degeneration index

            OUTPUT::
                - A polynomial that is the degeneration 'j' of self
            """
            return self.degeneration_morphism(j)(self)
        
        #@cached_function
        def propagate(self, face_list, dest_dim):
            """
            Returns a polynomial in dimension 'dest_dim' that has 'self' as a face

            INPUT:
                - face_list: indexes of the elements that form the face of the simplex we want to propagate
                - dest_dim: Dimension of the codomain cochain algebra
                - EXAMPLE:  s = SimplicialAlgebraplN(2)
                            s.inject_variables()
                            p = (t1^2+t2^2)*y1*y2
                            p.propagate2([0,1,2],3)
                            t1^2*t2*y0*y1 + t2^3*y0*y1 - t1^3*y0*y2 - t1*t2^2*y0*y2 + t0*t1^2*y1*y2 + t0*t2^2*y1*y2

            """
            if not (self.is_homogeneous() and self.total_degree() == self.parent().n):
                raise ValueError(
                    "The polynomial is not homogeneous of maximum degree.")
            self_dim = self.parent().n
            indexes = list(range(self_dim+1))
            apl = SimplicialAlgebraplN(dest_dim)
            y_gens = [apl.gens()[dest_dim+1+i] for i in face_list]
            t_gens = [apl.gens()[i] for i in face_list]
            graded_part = apl(0)
            for i in reversed(indexes):
                t = t_gens[i]
                y_prod = deepcopy(y_gens)
                del y_prod[i]
                y_prod = prod(y_prod)
                # ((-1)**(simplex_dim+i))*t*y_prod
                graded_part += ((-1)**i)*t*y_prod
            terms = self.dict()
            coeff_in_t = apl(0)
            for exp in terms.keys():
                coef = terms[exp]
                ts_exp = exp[:self_dim+1]
                coeff_in_t += coef * \
                    prod(list(map(lambda t, e: t**e, t_gens, ts_exp)))
            return coeff_in_t*graded_part
        
        """
        def _add_kernel_terms_to_primitive_basis(self, basis):
            new_elements = []
            for e in basis:
                monomial = list(deepcopy(e))
                n = len(e) // 2
                monomial[:n] = n*[0]
                new_elements.append(tuple(monomial))
            return new_elements
        """
        
        
        def _get_primitive_basis(self):
            """
            Returns a basis of the primitive of 'self'
            """
            basis = []
            p = self.reduce()
            apl = self.parent()
            for monomial in p.dict().keys():
                basis += apl.get_primitive_basis_of_a_monomial(apl,monomial)
            basis = list(set(basis))
            #basis += p._add_kernel_terms_to_primitive_basis(basis)
            t_degree = 2
            y_degree = self.total_degree() #self.parent().n #p.total_degree()#p.degree()[0]
            basis += apl.add_kernel_terms_to_primitive_basis(apl,t_degree, y_degree-1)
            basis = list(set(basis))
            terms = []
            apl = self.parent()
            for b in basis:
                #dim = len(b)//2 - 1
                terms.append(apl.get_term(apl,1, b, apl.gens()))
            terms.sort(reverse=False)
            return terms

        def _get_differential_basis(self):
            """
            Returns a basis of the differential of 'self'
            """
            apl = self.parent()
            primitive_basis = self._get_primitive_basis()
            differential_basis = []
            for b in primitive_basis:
                diff = b.differential().reduce()
                if diff:  # as we want a basis, we don't consider the coeffient of the differential
                    for monomial in diff.dict().keys():
                        term = apl.get_term(
                            apl,1, monomial, apl.gens())
                        differential_basis.append(term)
            differential_basis = list(set(differential_basis))
            differential_basis.sort(reverse=False)
            return differential_basis

        def _get_coordinate_vector_differential(self, diff_basis):
            """
            Returns the coordinates of 'self'
            in terms of the differential basis
            """
            p_reduced = self.reduce()
            #dic = p.dict()
            #terms = [p_reduced.parent()({t: dic[t]}) for t in dic.keys()]
            monomials = list(p_reduced.dict().keys())
            v = vector(QQ, len(diff_basis)*[0])
            for i in range(len(diff_basis)):
                t = diff_basis[i]
                term = list(t.dict().keys())[0]  # there is only one monomial
                if term in monomials:
                    v[i] = p_reduced.dict()[term]
            return v

        def _differential_matrix(self, primitive_basis, differential_basis):
            """
            Returns the matrix of the differential map
            """
            #primitive_basis = get_primitive_basis_of_a_polynomial(p.reduce())
            #differential_basis = get_differential_basis_of_a_polynomial(p.reduce())
            nrows = len(differential_basis)
            ncols = len(primitive_basis)
            M = zero_matrix(QQ, nrows, ncols, sparse=True)
            #print("zero_matrix created",flush=True)
            for j in range(ncols):
                #print("ncol :"+str(j),flush=True)
                primitive = primitive_basis[j]
                #diff_dict = self._tuple_differential(primitive)
                diff_primitive = primitive.differential().reduce()
                monomials = list(diff_primitive.dict().keys())
                for i in range(nrows):
                    #print("nrow :"+str(i),flush=True)
                    diff_basis_element = differential_basis[i]
                    term = list(diff_basis_element.dict().keys())[0]  # there is only one monomial
                    if term in monomials:
                        M[i, j] = diff_primitive.dict()[term]
            return M

        def solve_primitive(self):
            p = self.reduce()
            primitive_basis = p._get_primitive_basis()
            diff_basis = p._get_differential_basis()
            v = p._get_coordinate_vector_differential(diff_basis)
            Mdif = p._differential_matrix(primitive_basis, diff_basis)
            sol = Mdif.solve_right(v)
            return sum(map(lambda x, y: x*y, sol, primitive_basis))

@cached_function
def SimplicialAlgebraplN(n):
    """
    Returns the graded commutative differential cochain algebra (Apl)_n as part of
    the simplicial cochain algebra A_{PL}

    INPUT:

    - ``n`` -- A non-negative integer corresponding to the number of generators of degree 0.
    """

    if (n < 0):
        raise ValueError('`n` must be a non-negative integer')
    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)])
    gca = GradedCommutativeAlgebra(QQ, names=ts+ys, degrees=degrees)
    differentials = {gca.gens_dict()[ts[i]]: gca.gens_dict()[
        ys[i]] for i in range(n+1)}
    return SimplicialAplN(gca, differentials)


#############################################
# TEST
#print("TEST SimplicialAlgebraN")
# s = SimplicialAlgebraplN(5)
# s.inject_variables()
# ## face and degeneration maps
# print("Test face and degeneration maps")
# p = t1+y2
# p.degeneration(0).face(1)==p
# m1 = p.degeneration(1).face_morphism(1)
# m2 = p.degeneration_morphism(1)
# print((m1*m2).is_identity()) # true
# ## Propagation
# print("Test progatation")
# s = SimplicialAlgebraplN(2)
# s.inject_variables()
# p = (t1**2+t2**2)*y1*y2
# q = p.propagate([1,5,7],10)
# print(q)
# print(q.face(10).face(9).face(8).face(6).face(4).face(3).face(2).face(0).reduce() == p)
#############################################

class AlgebraAplK(UniqueRepresentation, Parent):
    """

      The graded commutative differential cochain algebra Apl(K)

    """

    def __init__(self, simplicial_complex, category=None):
        """
        The constructor of Apl(K)

        INPUT:
            - simplicial_complex: The simplicial complex used to build the algebra

        """
        self.K = simplicial_complex
        self.n = simplicial_complex.dimension()
        Parent.__init__(self, category=category)
        f = MorphismFromQQToAplK(self)
        self.register_coercion(f)

    def _repr_(self):
        """
        Return the base string representation of ``self``.
        """
        return "The graded commutative differential cochain algebra Apl(K) of %s" % self.K

    def _element_constructor_(self, values):
        """
        The constructor of an element of this class
        """
        return self.element_class(self, values)

    def dimension(self):
        """
        The dimension of the algebra. It is equivalent to the maximum degree possible that a non-null element of this algebra can take.
        """
        return self.n

    def simplicialComplex(self):
        """
        The simplicial complex of this algebra.
        """
        return self.K

    def _n_cells_sorted(self, n):
        """
        Returns all simplices of dimension 'n', as a sorted list.
        """
        return self.K._n_cells_sorted(n)

    def _faces_of_length(self, simplex, length):
        """
        Returns all faces of dimension 'length' for a given simplex 's'
        """
        return [tuple(l) for l in list(powerset(simplex)) if len(l)==length]
        #return [tuple(e) for e in Set(simplex).subsets(length)]

    @cached_function
    def zero(self):
        """
        Returns the zero element of the algebra
        """
        values = {}
        for i in range(self.K.dimension()+1):
            s = SimplicialAlgebraplN(i)
            for f in self.K._n_cells_sorted(i):
                values[tuple(f)] = s(0)
        return self._element_constructor_(values)

    def one(self):
        """
        Returns the one element of the algebra
        """
        values = {}
        for i in range(self.K.dimension()+1):
            s = SimplicialAlgebraplN(i)
            for f in self.K._n_cells_sorted(i):
                values[tuple(f)] = s(1)
        return self._element_constructor_(values)

    null_element = zero

    def lift_cochain(self, cochain):
        """
        Returns an element of the algebra Apl(K) whose integral is the cochain given as input.
        """
        # first we find polynomials compatible with the integration value
        values = {}
        max_dim = 0
        for k in cochain.keys():
            dim = len(k)-1
            if dim > max_dim:
                max_dim = dim
            apl = SimplicialAlgebraplN(dim)
            if dim > 0:
                values[k] = ((1) ** ((dim * (dim-1)) // 2)) * \
                    apl.integrand(cochain[k])
            else:
                values[k] = apl(cochain[k])
        
        # second, we propagate the simplicial information
        simplicial_structure = {}
        for n in reversed(range(max_dim+1)):
            keys = [k for k in values.keys() if len(
                k) == n+1]  # keys of the n-simplices
            for k in keys:
                simplicial_structure[k] = values[k] + \
                    self._propagate(cochain, k, values)
        # print(simplicial_structure)

        return AlgebraAplK.Element(self, simplicial_structure)

    def _propagate(self, cochain, simplex, values): 
        """
        Returns the part of the polynomial assigned to 'simplex' related to the
        propagation of the simplicial information of its faces
        """
        dim = len(simplex)-1
        apl = SimplicialAlgebraplN(dim)
        propagation = apl(0)
        for n in reversed(range(dim)):
            faces = self._faces_of_length(simplex, n+1)
            apl = SimplicialAlgebraplN(n)
            # we propagate the polynomial (y_1...y_n)
            #pol = apl(prod(apl.gens()[n+2:]))
            dict_set_keys = {c:set(c) for c in values.keys()}
            for f in faces: #CUIDADO AQUÍ CON CÓMO ESTÁN ORDENADOS LOS IDENTIFICADORES DE LOS SIMPLICES                
                if set(f) in dict_set_keys.values():
                    #print("Dim "+str(n), "cochain: ",cochain[f],"pol: ", pol, flush=True)
                    pol = values[f]
                    face = [k for k, v in dict_set_keys.items() if v == set(f)][0] # get_keys_from_value(dict, value)
                    #value = cochain[face]
                    face_indexes = [simplex.index(i) for i in face]
                    #print("dim: ",n," pol: ",pol," face_idx: ",face_indexes)
                    #propagation += value * \
                    #    ((1) ** ((n * (n-1)) // 2)) * factorial(n) * \
                    #    pol.propagate(face_indexes, dim)
                    propagation += pol.propagate(face_indexes, dim)
        return propagation

    @cached_function
    def cohomology_gens(self, degree):
        """
        Returns the cohomolosy generators of a given degree
        """
        print("Cohomology_gens, degree " + str(degree), flush=True)
        gens = {}
        basis = self.K.cohomology_ring(QQ).basis()
        for k in basis.keys():
            if k[0] == degree:
                gens[k] = basis[k]
        return gens

    @cached_function
    def cohomology_group(self, n):
        print("Cohomology_group, degree " + str(n), flush=True)
        chainComplex = self.K.chain_complex(base_ring=QQ)
        cochainComplex = chainComplex.dual()
        Hn = cochainComplex.differential(n).right_kernel(
        )/cochainComplex.differential(n-1).column_space()
        return Hn

    def cochain_to_vector(self, cochain, dim):
        """
        Converts a cochain, given as a dictionary, to a vector
        """
        basis = self.K.n_chains(dim, base_ring=QQ, cochains=True).basis()
        v = vector(QQ, len(basis)*[0])
        i = 0
        for simplex in basis.keys():
            v[i] = cochain[tuple(simplex)]
            i = i+1
        return v

    def vector_to_cochain(self, v, dim):
        """
        Converts a cochain, given as a vector, to a dictionary
        """
        cochain = {}
        for d in range(self.K.dimension()+1):
            basis = self.K.n_chains(d, base_ring=QQ, cochains=True).basis()
            if d != dim:
                for simplex in basis.keys():
                    cochain[tuple(simplex)] = 0
            else:
                i = 0
                for simplex in basis.keys():
                    cochain[tuple(simplex)] = v[i]
                    i = i+1
        return cochain

    def lift_cochain_complex(self, cochain, dim):
        """
            (REVISAR)
        Returns an element of the preimage by the differential of a given cochain

        INPUT:
            - cochain: a cochain as a vector
            - dim: dimension of the cochain
        """
        # S =  the Simplicial Complex
        # cochain = the cochain you want to find a preimage by the differential
        # dim = the dimension of the cochain
        cochain_vector = self.cochain_to_vector(cochain, dim)
        cochainComplex = self.K.chain_complex(base_ring=QQ).dual()
        diff = cochainComplex.differential(dim-1)
        #show(diff.echelon_form())
        result = diff.solve_right(cochain_vector)
        return self.vector_to_cochain(result)

    def lift_cohomology_generator(self, g):
        """
        Returns a representative cochain associated to a given cohomology generator.
        """
        cycle = g.to_cycle()
        cochain = {}
        for n in range(self.K.dimension()+1):
            for simplex in self.K._n_cells_sorted(n):
                cochain[simplex.tuple()] = 0
        for i in cycle:
            # print(type(i[0].tuple()))
            cochain[i[0].tuple()] = i[1]
        return cochain

    @cached_function
    def get_cohomology_aplk(self, degree):
        """
        Returns a set of representative elements of Apl(K), one for each cohomology generator of a given degree
        """
        #print("Cohomology_aplk, degree " + str(degree), flush=True)
        representative_aplk_elements = {}
        dic_gens = self.cohomology_gens(self,degree)
        for key in dic_gens.keys():
            gen = dic_gens[key]
            cochain = self.lift_cohomology_generator(gen)
            representative_aplk_elements[key] = self.lift_cochain(cochain)
        return representative_aplk_elements

    def first_cohomology_group(self):
        """
        Returns the degree >= 1 of the first non-null cohomology generator.
        """
        basis = self.K.cohomology_ring(base_ring=QQ).basis()
        indexes = list(basis.keys())
        indexes.sort(reverse=False)
#         print(indexes)
        for n in range(1, len(indexes)):
            if indexes[n][0] > 0:
                return indexes[n][0]
        return -1

    class Element(Element):
        """
        Element class of the simplicial commutative differential algebra (Apl)n.
        """

        def __init__(self, parent, values):
            Element.__init__(self, parent)
            self.elements = self._reduce(values)
            #self.check()

        def _repr_(self):
            """
            Return the base string representation of ``self``.
            """
            return str({a: b for a, b in self.elements.items() if not b.is_zero()})
            # return str(self.elements)

        def dimension(self):
            """
            Returns the dimension where 'self' lives
            """
            return self.parent().dimension()

        def _add_(self, other):
            """
            Addition operation
            """
            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):
            """
            Substration operation
            """
            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):
            """
            Multiplication operation
            """
            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 __pow__(self, n):
            """
            Return ``self`` raised to the power ``n``.
            """
            C = self.__class__
            newValues = {}
            for k in self.elements.keys():
                newValues[k] = self.elements[k]**n
            return C(self.parent(), newValues)

        def _richcmp_(self, other, op):
            """
            Comparison operations. Only implemented equality operations (EQ and NE)
            """
            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 _reduce(self, values):
            """
                Returns the polynomials given as input, after applying the corresponding relations
                of the algebra where each polynomial belongs.

                INPUT:
                    - values: a dictionary where the values are polynomials in (Apl)n
            """
            reduced_values = {}
            for k in values.keys():
                n = len(k)-1
                reduced_values[k] = values[k].reduce()
            return reduced_values

        def check(self):
            """
            Checks if the initial assignation respects the simplicial structure. For every simplex S, the algorithm checks
            if the polynomial assigned to each face of S is equal to the same face of the polynomial assigned to S.
            """
            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.elements[simplex.tuple()].face(f).reduce()
                        #print("face: " + str(f),"simplex: " + str(simplex),"sFace: " + str(self.elements[sFace]), "aFace: " + str(aFace))
                        if(self.elements[sFace] != aFace):
                            raise ValueError("the assignation is not valid")
            #print("Check ok", flush=True)
            return True

        def differential(self):
            """
            Returns the differential of 'self'
            """
            C = self.__class__
            newValues = {}
            for k in self.elements.keys():
                newValues[k] = self.elements[k].differential()
            return C(self.parent(), newValues)

        def is_closed(self):
            """
            Returns 'True' if 'self' is closed, i.e., its differential is zero.
            """
            diff = self.differential()
            if(diff == self.parent().null_element()):
                return True
            else:
                return False

        def homogeneous_part(self, degree):
            """
            Returns the homogeneous component of grade 'degree' of 'self'
            """
            C = self.__class__
            part = {}
            for k in self.elements.keys():
                s = SimplicialAlgebraplN(len(k)-1)
                part[k] = s(self.elements[k].homogeneous_part(degree=degree))
            return C(self.parent(), part)

        def homogeneous_parts(self):
            """
            Returns all graded components of 'self'
            """
            parts = {}
            for n in range(self.dimension()+1):
                hom_part = self.homogeneous_part(degree=n)
                if hom_part != self.parent().null_element():
                    parts[n] = hom_part
            return parts

        def is_homogeneous(self):
            """
            Returns 'True' if self is composed by homogeneous polynomials of the same degree
            """
            isHomogeneous = True
            S = self.parent().simplicialComplex()
            dim = S.dimension()
            for n in range(dim+1):
                for simplex in S._n_cells_sorted(n):
                    pol = self.elements[simplex.tuple()]
                    if pol:
                        isHomogeneous = isHomogeneous and pol.is_homogeneous()
                        if not isHomogeneous:
                            break
                if not isHomogeneous:
                    break
            return isHomogeneous

        def degree(self):
            if not self.is_homogeneous():
                return -1
            else:
                degree = -1
                for e in self.elements.values():
                    if e:
                        degree = e.total_degree()#e.degree()[0]
                        break
                return degree

        def integral(self):
            """
            Returns the cochain obtained as integration of 'self'
            """
            cochain = {}
            # integration of each polynomial
            for k in self.elements.keys():
                dim = len(k)-1
                # the sign is due to stoke's theorem
                cochain[k] = ((1) ** ((dim*(dim-1)) // 2)) * \
                    self.elements[k].integral()
            return cochain

        def _face_basis(self, basis, face_idx):
            face_basis = []
            for b in basis:
                face_pol = b.face(face_idx).reduce()
                apl = face_pol.parent()
                if face_pol:
                    for monomial in face_pol.dict().keys():
                        # remove coefficients
                        term = apl.get_term(
                            apl,1, monomial, apl.gens())
                        face_basis.append(term)
            face_basis = list(set(face_basis))
            face_basis.sort(reverse=False)
            return face_basis

        def _get_primitive_basis(self):

            primitive_basis = {}
            S = self.parent().simplicialComplex()
            dim = S.dimension()
            min_dim = self.degree()
            for n in reversed(range(min_dim, dim+1)):
                print("Primitive Basis: Working on dimension " + str(n), flush=True)
                cells_sorted = S._n_cells_sorted(n)
                for simplex in cells_sorted:
                    #print("Primitive Basis: Working on simplex " + str(simplex), flush=True)
                    #print("simplex "+str(simplex)+ " is at position " +str(cells_sorted.index(simplex)) + " out of " + str(len(cells_sorted)), flush=True)
                    if simplex.tuple() not in primitive_basis.keys():
                        primitive_basis[simplex.tuple()] = []
                    # add primitive elements related to the assigned polynomial
                    pol = self.elements[simplex.tuple()]
                    primitive_basis[simplex.tuple()] += pol._get_primitive_basis()
                    primitive_basis[simplex.tuple()] = list(set(primitive_basis[simplex.tuple()]))
                    primitive_basis[simplex.tuple()].sort(reverse=False)
                    #print(primitive_basis[simplex.tuple()], flush=True)
                    # propagate the primitive basis through its faces
                    for f in range(n+1):
                        simplex_face = S.face(simplex, f).tuple()
                        if simplex_face not in primitive_basis.keys():
                            primitive_basis[simplex_face] = []
                        primitive_basis[simplex_face] += self._face_basis(primitive_basis[simplex.tuple()], f)
                        primitive_basis[simplex_face] = list(set(primitive_basis[simplex_face]))
                        primitive_basis[simplex_face].sort(reverse=False)
            return primitive_basis

        def _diff_basis(self, basis):
            diff_basis = []
            for b in basis:
                b_diff = b.differential().reduce()
                if b_diff:
                    apl = b_diff.parent()
                    for monomial in b_diff.dict().keys():
                        # remove coefficients
                        term = apl.get_term(
                            apl,1, monomial, apl.gens())
                        diff_basis.append(term)
            diff_basis = list(set(diff_basis))
            diff_basis.sort(reverse=False)
            return diff_basis

        def _get_differential_basis(self, primitive_basis):
            differential_basis = {}
            S = self.parent().simplicialComplex()
            dim = S.dimension()
            min_dim = self.degree()
            for n in reversed(range(min_dim, dim+1)):
                for simplex in S._n_cells_sorted(n):
                    prim_diff_basis = primitive_basis[simplex.tuple()]
                    differential_basis[simplex.tuple()] = self._diff_basis(prim_diff_basis)
            return differential_basis

        def _get_differential_matrix(self, p_basis, d_basis):
            S = self.parent().simplicialComplex()
            dim = S.dimension()
            min_dim = self.degree()
            ncols = sum([len(l) for l in p_basis.values()])
            nrows = sum([len(l) for l in d_basis.values()])
            print("nrows: "+str(nrows),flush=True)
            print("ncols: "+str(ncols),flush=True)
            M = zero_matrix(QQ, nrows, ncols, sparse=True)
            row_offset = 0
            col_index = 0
            for n in reversed(range(min_dim-1, dim+1)):
                print("_get_differential_matrix: dimension "+str(n),flush=True)
                for simplex in S._n_cells_sorted(n):
                    #print(simplex.tuple(), row_offset, col_index,flush=True)
                    if simplex.tuple() in p_basis.keys():
                        primitives = p_basis[simplex.tuple()]
                        if simplex.tuple() in d_basis.keys():
                            differentials = d_basis[simplex.tuple()]
                            if primitives and differentials:
                                for p in primitives:
                                    dif_p = p.differential().reduce()
                                    apl = dif_p.parent()
                                    dic = dif_p.dict()
                                    for monomial in dic.keys():
                                        term = apl.get_term(
                                            apl,1, monomial, apl.gens())
                                        if term in differentials:
                                            row_index = differentials.index(
                                                term)
                                            M[row_offset + row_index,
                                                col_index] = dic[monomial]
                                    col_index += 1
                                row_offset += len(differentials)
                            elif primitives:
                                col_index += len(primitives)
                        elif primitives:
                            col_index += len(primitives)

            return M

        def _get_vector_diff_coordinates(self, d_basis):
            len_v = sum([len(l) for l in d_basis.values()])
            v = zero_vector(QQ, len_v)
            offset = 0
            S = self.parent().simplicialComplex()
            dim = S.dimension()
            min_dim = self.degree()
            for n in reversed(range(min_dim, dim+1)):
                print("Vector diff coordinates: dimension " + str(n), flush=True)
                for s in S._n_cells_sorted(n):
                    simplex = s.tuple()
                    if simplex in d_basis.keys():
                        differentials = d_basis[simplex]
                        pol = self.elements[simplex].reduce()
                        apl = pol.parent()
                        if pol and differentials:
                            dic = pol.dict()
                            for monomial in dic.keys():
                                term = apl.get_term(
                                    apl,1, monomial, apl.gens())
                                if term in differentials:
                                    v_index = differentials.index(term)
                                    v[offset + v_index] = dic[monomial]
                        offset += len(differentials)
                        #print("Vector diff coordinates: offset " + str(offset), flush=True)
            return v

        def _get_offsets_matrix(self, p_basis):
            offsets = {}
            S = self.parent().simplicialComplex()
            dim = S.dimension()
            min_dim = self.degree()
            offset = 0
            first_simplex = True
            for n in reversed(range(min_dim-1, dim+1)):
                for simplex in S._n_cells_sorted(n):
                    if simplex.tuple() in p_basis.keys():
                        primitives = p_basis[simplex.tuple()]
                        if first_simplex:
                            offsets[simplex.tuple()] = 0
                            first_simplex = False
                        else:
                            offsets[simplex.tuple()] = offset
                        offset += len(primitives)
            return offsets

        def _get_row_face_matrix(self, primitives, face_idx, face_basis,
                                 offset_simplex, offset_face, len_row):
            M = zero_matrix(QQ, 1, len_row, sparse=True)
            #print("_get_row_face_matrix",flush=True)
            #print("len(face_basis): "+str(len(face_basis)),flush=True)
            for fb in face_basis:
                #print("fb: "+str(fb),flush=True)
                fb_idx = face_basis.index(fb)
                row = zero_matrix(QQ, 1, len_row, sparse=True)
                row[0, offset_face+fb_idx] = -1
                for term in primitives:
                    #print("term: "+str(term),flush=True)
                    face_term = term.face(face_idx).reduce()
                    apl = face_term.parent()
                    dic = face_term.dict()
                    i = primitives.index(term)
                    for monomial in dic:
                        t = apl.get_term(
                            apl,1, monomial, apl.gens())
                        if t == fb:
                            #print("fb = "+str(t),flush=True)
                            coef = dic[monomial]
                            row[0, offset_simplex+i] = coef
                M = M.stack(row)
            return M[1:]

        def _get_face_matrix(self, p_basis):

            offsets = self._get_offsets_matrix(p_basis)
            len_row = sum([len(l) for l in p_basis.values()])
            M = zero_matrix(QQ, 1, len_row,sparse=True)
            S = self.parent().simplicialComplex()
            dim = S.dimension()
            min_dim = self.degree()
            for n in reversed(range(min_dim, dim+1)):
                print("Face matrix: Working on dimension " + str(n), flush=True)
                cells_sorted = S._n_cells_sorted(n)
                for simplex in cells_sorted:
                    #print("simplex "+str(simplex)+ " is at position " +str(cells_sorted.index(simplex)) + " out of " + str(len(cells_sorted)), flush=True)
                    offset_simplex = offsets[simplex.tuple()]
                    primitives = p_basis[simplex.tuple()]
                    for f in range(n+1):
                        #print("Face matrix: Working on dimension " + str(n)+" and face "+str(f), flush=True)
                        simplex_face = S.face(simplex, f).tuple()
                        face_primitives = p_basis[simplex_face]
                        offset_face = offsets[simplex_face]
                        rows = self._get_row_face_matrix(primitives, f, face_primitives,
                                                         offset_simplex, offset_face,
                                                         len_row)
                        M = M.stack(rows)
                        #print("Face Matrix Dimensions " + str(M.dimensions()), flush=True)
            return M[1:]

        def _get_solution(self, sol, p_basis):
            S = self.parent().simplicialComplex()
            dim = S.dimension()
            min_dim = self.degree()-1
            offset = 0
            values = {}
            for n in reversed(range(dim+1)):
                apl = SimplicialAlgebraplN(n)
                for s in S._n_cells_sorted(n):
                    simplex = s.tuple()
                    pol = apl(0)
                    if simplex in p_basis.keys():
                        primitives = p_basis[simplex]
                        if primitives:
                            for p in primitives:
                                i = primitives.index(p)
                                coef = sol[offset+i]
                                pol += coef*p
                            offset += len(primitives)
                    values[simplex] = pol
            A = self.parent()
#             return values
            return A._element_constructor_(values)

        def solve_primitive(self):
            print("start _get_primitive_basis()",flush=True)
            p_basis = self._get_primitive_basis()
            print("length primitive basis: ",len(p_basis),flush=True)
            d_basis = self._get_differential_basis(p_basis)  
            print("length differential basis: ",len(d_basis),flush=True)
            diff_matrix = self._get_differential_matrix(p_basis, d_basis)
            print("_get_vector_diff_coordinates()", flush=True)
            v_dif = self._get_vector_diff_coordinates(d_basis)
            print("_get_face_matrix()", flush=True)
            face_matrix = self._get_face_matrix(p_basis)
            nrows = face_matrix.nrows()
            v = vector(QQ, list(v_dif)+nrows*[0])
            M = diff_matrix.stack(face_matrix)
            print("Matrix size: ", M.nrows(), M.ncols(),flush=True)
            sol = M.solve_right(v)
            print("Solution found: building element",flush=True)
            return self._get_solution(sol, p_basis)


class MorphismFromQQToAplK(Morphism):

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

    def _call_(self, x):
        S = self.simplicialComplex
        params = {}
        for d in range(S.dimension()+1):
            s = SimplicialAlgebraplN(d)
            for c in S._n_cells_sorted(d):
                params[tuple(c)] = s(x)
        A = AlgebraAplK(S)
        return AlgebraAplK.Element(A, params)