adding classes and methods

aplK.py

@@ -2,6 +2,7 @@ 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
 import copy
 
 class AplN(DifferentialGCAlgebra_multigraded, UniqueRepresentation):
@@ -463,6 +464,114 @@ class SimplicialAplN(AplN, UniqueRepresentation):
                 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(copy.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(copy.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):
     """
@@ -574,6 +683,156 @@ class AlgebraAplK(UniqueRepresentation, Parent):
                 values[tuple(f)]=s(0)
         return self._element_constructor_(values)
 
+    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.
@@ -761,16 +1020,233 @@ class AlgebraAplK(UniqueRepresentation, Parent):
                         break
                 return degree
             
-        def get_skeleton(self, dim):
+        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)
+            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 = {}
-            facets = self._get_elements_of_dimension(dim)
+            for n in reversed(range(dim+1)):
-            for i in range(dim+1):
+                apl = SimplicialAlgebraplN(n)
-                values.update(self._get_elements_of_dimension(i))
+                for s in S._n_cells_sorted(n):
-            S = SimplicialComplex(list(facets.keys()),is_mutable=False)
+                    simplex = s.tuple()
-            A = AlgebraAplK(S)
+                    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)
+            sol = M.solve_right(v)
+            return self._get_solution(sol, p_basis)
     
+
 class MorphismFromQQToAplK(Morphism):
 
     def __init__(self, codomain):