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

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:])
        else:
            raise ValueError("only field elements can be integrated.")

    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.
            """

            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.
            """

            degrees = set([k[0] for k in self.homogeneous_parts()])
            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))

        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()
            pol = apl(self.reduce()) #needed in case reduce() returns an integer
            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_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)))

        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 =  self._get_term(1, exp,ys_only_gens)
                coef = self._get_term(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 self._get_term(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)
            aGens = polynomial.parent().gens()
            primitive = 0
            for exp in polynomial.dict().keys():
                coef = polynomial.dict()[exp]
                monomial = self._get_term(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.")
            w = self._leading_primitive()
            theta = self - w.differential()
            if (theta == 0):
                return w
            else:
                return w+theta.primitive()

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)

        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)

        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)
                graded_part += ((-1)**i)*t*y_prod#((-1)**(simplex_dim+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 _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))

        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()
            for monomial in p.dict().keys():
                basis += p._get_primitive_basis_of_a_monomial(monomial)
            basis = list(set(basis))
            basis += p._add_kernel_terms_to_primitive_basis(basis)
            basis = list(set(basis))
            terms = []
            apl = self.parent()
            for b in basis:
                dim = len(monomial)//2 - 1
                terms.append(self._get_term(1, b, apl.gens()))
            terms.sort(reverse=False)
            return terms

        def _get_differential_basis(self):
            """
            Returns a basis of the differential of 'self'
            """
            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 = self._get_term(1, monomial, self.parent().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 = matrix(QQ, nrows, ncols, nrows*ncols*[0])
            for j in range(ncols):
                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):
                    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))

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
        """
        return [tuple(e) for e in Set(simplex).subsets(length)]

    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)
        #print(simplicial_structure)
        return AlgebraAplK.Element(self,simplicial_structure)

    def _propagate(self, cochain, simplex):
        """
        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)
            pol = apl(prod(apl.gens()[n+2:])) # we propagate the polynomial (y_1...y_n)
            for f in faces:
                if f in cochain.keys():
                    value = cochain[f]
                    face_indexes = [simplex.index(i) for i in f]
                    #print("dim: ",n," pol: ",pol," face_idx: ",face_indexes," dim: ", n)
                    propagation += value*((-1) ** ((n * (n-1)) // 2)) * factorial(n) * pol.propagate(face_indexes, dim)
        return propagation

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

    def cohomology_group(self, n):
        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

    def get_cohomology_aplk(self,degree):
        """
        Returns a set of representative elements of Apl(K), one for each cohomology generator of a given degree
        """
        representative_aplk_elements = {}
        dic_gens = self.cohomology_gens(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")

        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.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()
                if face_pol:
                    for monomial in face_pol.dict().keys():
                        #remove coefficients
                        term = face_pol._get_term(1, monomial, face_pol.parent().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)):
                for simplex in S._n_cells_sorted(n):
                    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)
                    # 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:
                    for monomial in b_diff.dict().keys():
                        #remove coefficients
                        term = b_diff._get_term(1, monomial, b_diff.parent().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()])
            M = zero_matrix(QQ, nrows, ncols)
            row_offset = 0
            col_index = 0
            for n in reversed(range(min_dim-1,dim+1)):
                for simplex in S._n_cells_sorted(n):
#                     print(simplex.tuple(), row_offset, col_index)
                    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()
                                    dic = dif_p.dict()
                                    for monomial in dic.keys():
                                        term  = dif_p._get_term(1, monomial, dif_p.parent().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)):
                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()
                        if pol and differentials:
                            dic = pol.dict()
                            for monomial in dic.keys():
                                term  = pol._get_term(1, monomial, pol.parent().gens())
                                if term in differentials:
                                    v_index = differentials.index(term)
                                    v[offset + v_index] = dic[monomial]
                        offset += len(differentials)
            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)
            for fb in face_basis:
                fb_idx = face_basis.index(fb)
                row = zero_matrix(QQ,1, len_row)
                row[0,offset_face+fb_idx] = -1
                for term in primitives:
                    face_term = term.face(face_idx).reduce()
                    dic = face_term.dict()
                    i = primitives.index(term)
                    for monomial in dic:
                        t  = face_term._get_term(1, monomial, face_term.parent().gens())
                        if t == fb:
                            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)
            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):
                    offset_simplex = offsets[simplex.tuple()]
                    primitives = p_basis[simplex.tuple()]
                    for f in range(n+1):
                        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)
            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):
            p_basis = self._get_primitive_basis()
            d_basis =  self._get_differential_basis(p_basis)
            diff_matrix = self._get_differential_matrix(p_basis, d_basis)
            v_dif = self._get_vector_diff_coordinates(d_basis)
            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())
            sol = M.solve_right(v)
            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)