diff --git a/aplK.py b/aplK.py index b035dfd..dea1778 100644 --- a/aplK.py +++ b/aplK.py @@ -83,6 +83,58 @@ class AplN(DifferentialGCAlgebra_multigraded, UniqueRepresentation): return value*factorial(self.n)*prod(self.gens()[self.n+2:]) 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)] + V = VectorSpace(GF(tdegree),N-1) + listCoefs = [tuple([0]+list(v)+N*[0]) for v in V] + 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): """ @@ -179,21 +231,6 @@ class AplN(DifferentialGCAlgebra_multigraded, UniqueRepresentation): # 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): """" @@ -210,8 +247,8 @@ class AplN(DifferentialGCAlgebra_multigraded, UniqueRepresentation): 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) + 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: @@ -270,7 +307,7 @@ class AplN(DifferentialGCAlgebra_multigraded, UniqueRepresentation): primitive[idx] += 1 new_coef = coef/primitive[idx] a = monomial.parent() - return self._get_term(new_coef, primitive, a.gens()) + return a.get_term(a,new_coef, primitive, a.gens()) def _leading_primitive(self): """" @@ -283,11 +320,12 @@ class AplN(DifferentialGCAlgebra_multigraded, UniqueRepresentation): # 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() + apl = polynomial.parent() + aGens = apl.gens() primitive = 0 for exp in polynomial.dict().keys(): coef = polynomial.dict()[exp] - monomial = self._get_term(coef, exp, aGens) + 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()]) @@ -302,6 +340,9 @@ class AplN(DifferentialGCAlgebra_multigraded, UniqueRepresentation): """ 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): @@ -503,23 +544,6 @@ class SimplicialAplN(AplN, UniqueRepresentation): 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): @@ -531,47 +555,28 @@ class SimplicialAplN(AplN, UniqueRepresentation): new_elements.append(tuple(monomial)) return new_elements """ - - def findsubsests(S,m): - return set(itertools.combinations(S,m)) - - def _add_kernel_terms_to_primitive_basis(self, tdegree, ydegree): - apl = self.parent() - N = len(apl.gens())//2 - ygens = apl.gens()[N+1:] #N+1 because we remove y0 generator - listProds = [prod(s) for s in findsubsets(ygens,ydegree)] - V = VectorSpace(GF(tdegree),N-1) - listCoefs = [tuple([0]+list(v)+N*[0]) for v in V] - listPolinomials = [self._get_term(1, c, apl.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] - def _get_primitive_basis(self): """ Returns a basis of the primitive of 'self' """ - if not self: - return [] basis = [] p = self.reduce() + apl = self.parent() for monomial in p.dict().keys(): - basis += p._get_primitive_basis_of_a_monomial(monomial) + 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 = p.total_degree()#p.degree()[0] - #basis += p._add_kernel_terms_to_primitive_basis(t_degree, y_degree-1) + y_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(self._get_term(1, b, apl.gens())) + terms.append(apl.get_term(apl,1, b, apl.gens())) terms.sort(reverse=False) return terms @@ -579,14 +584,15 @@ class SimplicialAplN(AplN, UniqueRepresentation): """ 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 = self._get_term( - 1, monomial, self.parent().gens()) + term = apl.get_term( + apl,1, monomial, apl.gens()) differential_basis.append(term) differential_basis = list(set(differential_basis)) differential_basis.sort(reverse=False) @@ -617,16 +623,18 @@ class SimplicialAplN(AplN, UniqueRepresentation): #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]) + 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 + 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 @@ -933,7 +941,7 @@ class AlgebraAplK(UniqueRepresentation, Parent): def __init__(self, parent, values): Element.__init__(self, parent) self.elements = self._reduce(values) - self.check() + #self.check() def _repr_(self): """ @@ -1041,6 +1049,7 @@ class AlgebraAplK(UniqueRepresentation, Parent): #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) def differential(self): """ @@ -1122,19 +1131,21 @@ class AlgebraAplK(UniqueRepresentation, Parent): 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)) * \ + cochain[k] = ((-1) ** (((dim-1) * (dim-2)) // 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 = face_pol._get_term( - 1, monomial, face_pol.parent().gens()) + term = apl.get_term( + apl,1, monomial, apl.gens()) face_basis.append(term) face_basis = list(set(face_basis)) face_basis.sort(reverse=False) @@ -1147,7 +1158,11 @@ class AlgebraAplK(UniqueRepresentation, Parent): dim = S.dimension() min_dim = self.degree() for n in reversed(range(min_dim, dim+1)): - for simplex in S._n_cells_sorted(n): + #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 @@ -1170,10 +1185,11 @@ class AlgebraAplK(UniqueRepresentation, Parent): 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 = b_diff._get_term( - 1, monomial, b_diff.parent().gens()) + term = apl.get_term( + apl,1, monomial, apl.gens()) diff_basis.append(term) diff_basis = list(set(diff_basis)) diff_basis.sort(reverse=False) @@ -1196,12 +1212,15 @@ class AlgebraAplK(UniqueRepresentation, Parent): 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) + 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) + #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(): @@ -1209,10 +1228,11 @@ class AlgebraAplK(UniqueRepresentation, Parent): 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 = dif_p._get_term( - 1, monomial, dif_p.parent().gens()) + term = apl.get_term( + apl,1, monomial, apl.gens()) if term in differentials: row_index = differentials.index( term) @@ -1235,20 +1255,23 @@ class AlgebraAplK(UniqueRepresentation, Parent): 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 = pol._get_term( - 1, monomial, pol.parent().gens()) + 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): @@ -1272,19 +1295,25 @@ class AlgebraAplK(UniqueRepresentation, Parent): def _get_row_face_matrix(self, primitives, face_idx, face_basis, offset_simplex, offset_face, len_row): - M = zero_matrix(QQ, 1, 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) + 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 = face_term._get_term( - 1, monomial, face_term.parent().gens()) + 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) @@ -1294,15 +1323,19 @@ class AlgebraAplK(UniqueRepresentation, Parent): 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) + 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)): - for simplex in S._n_cells_sorted(n): + #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] @@ -1310,6 +1343,7 @@ class AlgebraAplK(UniqueRepresentation, Parent): 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): @@ -1337,17 +1371,22 @@ class AlgebraAplK(UniqueRepresentation, Parent): return A._element_constructor_(values) def solve_primitive(self): + print("start _get_primitive_basis()",flush=True) p_basis = self._get_primitive_basis() - print("lenght primitive basis: ",len(p_basis)) - d_basis = self._get_differential_basis(p_basis) + print("lenght primitive basis: ",len(p_basis),flush=True) + d_basis = self._get_differential_basis(p_basis) + print("lenght 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()) + 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)