calquezar/AplK
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

Update 'aplK.py'

Browse Source
This commit is contained in:
calquezar 2024-02-26 12:12:47 +01:00
parent 13857b8aca
commit ca4194694c
1 changed files with 127 additions and 88 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

215
aplK.py
View File

@ -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)