diff --git a/aplK.py b/aplK.py index 4dbf70d..85c1494 100644 --- a/aplK.py +++ b/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 i in range(dim+1): - values.update(self._get_elements_of_dimension(i)) - S = SimplicialComplex(list(facets.keys()),is_mutable=False) - A = AlgebraAplK(S) + 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) + sol = M.solve_right(v) + return self._get_solution(sol, p_basis) + class MorphismFromQQToAplK(Morphism): def __init__(self, codomain):