From: Antonio Jimenez Pastor Date: Mon, 18 Mar 2019 16:13:41 +0000 (+0100) Subject: Finished a first version of lazyDDRing X-Git-Url: http://git.risc.jku.at/gitweb/?a=commitdiff_plain;h=059eb4da6ec876aacd58828474a4e2a211402e18;p=ajpastor%2Fdiff_defined_functions.git Finished a first version of lazyDDRing Now it can compute derivatives, create the lazy ring for DD-finite functions and include elements that are also D-finite. Fixed some errors in lazy/conversion.py and lazy/lazyRing.py --- diff --git a/ajpastor/dd_functions/lazyDDRing.py b/ajpastor/dd_functions/lazyDDRing.py index 39d5bee..f14e0eb 100644 --- a/ajpastor/dd_functions/lazyDDRing.py +++ b/ajpastor/dd_functions/lazyDDRing.py @@ -38,8 +38,11 @@ class _LazyDDFunction(IntegralDomainElement): try: self.__poly = parent.poly_field()(el); except: - self.__raw = parent.base()(el); - self.__poly = self.parent().to_poly(self.__raw); + try: + self.__poly = self.parent().to_poly(el); + except: + self.__raw = parent.base()(el); + self.__poly = self.parent().to_poly(self.__raw); self.simplify(); # We try to simplify the element from the beginning @@ -67,6 +70,7 @@ class _LazyDDFunction(IntegralDomainElement): def simplify(self): self.__poly = self.parent().simplify(self.poly()); + return self; def variables(self): ''' @@ -176,8 +180,8 @@ class _LazyDDFunction(IntegralDomainElement): @cached_method def is_zero(self): result = (self.poly() == 0 ); - map_to_values = {repr(var) : self.parent().to_real(var)(0)}; - if((not result) and (self(**map_to_values) == 0 )): + map_to_values = {repr(var) : self.parent().to_real(var)(0) for var in self.poly().variables()}; + if((not result) and (self.poly()(**map_to_values) == 0)): if(not (self.__raw is None)): result = (self.raw() == 0); else: @@ -265,6 +269,7 @@ class LazyDDRing (UniqueRepresentation, ConversionSystem, IntegralDomain): self.__poly_field = self.__poly_ring.fraction_field(); self.__gen = self.__poly_ring.gens()[0]; self.__used = 0; + self.stop = False; ## Initializing the map of variables self.__map_of_vars = {}; # Map where each variable is associated with a 'base' object @@ -298,7 +303,9 @@ class LazyDDRing (UniqueRepresentation, ConversionSystem, IntegralDomain): return tuple(self.poly_ring()(key) for key in self.map_of_vars().keys()); def _to_poly_element(self, element): - if(not (element in self.base())): + if(element.parent() is self): + return element.poly(); + elif(not (element in self.base())): raise TypeError("Element is not in the base ring.\n\tExpected: %s\n\tGot: %s" %(element.parent(), self.base())); ## We try to cast to a polynomial. If we can do it, it means no further operations has to be done @@ -377,6 +384,37 @@ class LazyDDRing (UniqueRepresentation, ConversionSystem, IntegralDomain): ## Otherwise, we add the element to the variables and return it self.__add_function(element, gen=True); return self.__gen[self.__map_to_vars[element]]; + + def _relations(self): + ''' + Returns a Groebner Basis of the relations ideals known in this conversion system. + ''' + if(self._ConversionSystem__relations is None): + self._ConversionSystem__relations = []; + self._ConversionSystem__rel_ideal = ideal([Integer(0)]); + + return self._ConversionSystem__relations; + + def _groebner_basis(self): + if(len(self._ConversionSystem__relations) >= 1): + gens_in_relations = list(set(sum([[str(v) for v in el.variables()] for el in self._ConversionSystem__relations], []))); + Poly = PolynomialRing(self.poly_ring().base_ring(), gens_in_relations, order='deglex'); + rels = [Poly(el) for el in self._ConversionSystem__relations]; + self._ConversionSystem__rel_ideal = ideal(Poly, rels); + return self._ConversionSystem__rel_ideal.groebner_basis(); + return [self.poly_ring().base_ring().zero()]; + + def _simplify(self, poly): + gens = list(set([str(g) for g in self._ConversionSystem__rel_ideal.parent().ring().gens() if g != 1] + [str(g) for g in poly.variables()])); + Poly = PolynomialRing(self.poly_ring().base_ring(), gens); + try: + poly = Poly(poly); + except TypeError: + return poly; + try: + return self.poly_ring()(poly.reduce([Poly(el) for el in self._ConversionSystem__relations])); + except AttributeError: + return self.poly_ring()(poly); ################################################################################################ ### Other Methods for LazyRing @@ -401,10 +439,36 @@ class LazyDDRing (UniqueRepresentation, ConversionSystem, IntegralDomain): self.__map_of_derivatives = {}; def derivative(self, el, times=1): - raise NotImplementedError("method derivative not implemented"); + ''' + Method that computes the derivative of an element in the LazyDDRing. It performs a casting to 'self' before starting the algorithm. + ''' + el = self(el); + + if(times == 0): + return el; + elif(times > 1): + return self.derivative(self.derivative(el, times-1)); + + ## Splitting the derivative into the monomials + poly = el.poly(); + if(poly == 0): + return self.zero(); + coeffs = poly.coefficients(); + monomials = poly.monomials(); + + ## Base case (monomial case) + if(len(monomials) == 1): + gens = poly.variables(); + degrees = [poly.degree(g) for g in gens]; + + res = sum(prod(gens[j]**(degrees[j]-kronecker_delta(i,j)) for j in range(len(gens)))*degrees[i]*self.get_derivative(gens[i]) for i in range(len(gens))); + + return self(coeffs[0]*res); + else: ## Generic case (polynomial) + return self(sum(coeffs[i]*self.derivative(monomials[i]) for i in range(len(coeffs)))); def get_derivative(self, el): - raise NotImplementedError("method get_derivative not implemented"); + return self(self.__map_of_derivatives[el]).poly(); ################################################################################################ ### Other Integral Domain methods @@ -489,7 +553,7 @@ class LazyDDRing (UniqueRepresentation, ConversionSystem, IntegralDomain): ## For each variable in X.poly(), we get the new polynomial translate = {} for var in X.variables(): - translate[var] = self.to_poly(other.map_of_vars()[var]); + translate[var] = self.to_poly(other.map_of_vars()[str(var)]); ## We now plugin the expressions return _LazyDDFunction(self, pol(**translate)); @@ -518,7 +582,7 @@ class LazyDDRing (UniqueRepresentation, ConversionSystem, IntegralDomain): def __str__(self): final = "%s with %d variables\n{\n" %(repr(self),len(self.__gens)); for g in self.__gens: - final += "\t%s : %s,\n" %(g, repr(self.__map_of_vars[g])); + final += "\t%s : %s,\n" %(g, repr(self.__map_of_vars[str(g)])); final += "}"; return final @@ -568,9 +632,12 @@ class LazyDDRing (UniqueRepresentation, ConversionSystem, IntegralDomain): break; ## Create the tranformation into f + + # Companion matrix + companion = f.equation.companion(); + companion = Matrix(self, [[self(companion[i][j]) for j in range(companion.ncols())] for i in range(companion.nrows())]); if(order < f.getOrder()): # Relation found C = []; - companion = f.equation.companion(); for i in range(order): new_row = [companion[i][j] for j in range(order)]; if(i == relation[0]): @@ -579,7 +646,7 @@ class LazyDDRing (UniqueRepresentation, ConversionSystem, IntegralDomain): C = Matrix(self,C); self.__trans[f] = (relation[1][1],C); else: - self.__trans[f] = (self.poly_ring().zero(),Matrix(self, f.equation.companion())); + self.__trans[f] = (self.poly_ring().zero(),Matrix(self, companion)); ## Adding f as a vertex (edges must be added outside this function) self.__r_graph.add_vertex(f); @@ -587,12 +654,18 @@ class LazyDDRing (UniqueRepresentation, ConversionSystem, IntegralDomain): ## If f will become a generator, add it to the used variables if(gen): current = f; + before = self.__used; for i in range(order): + print "** Adding a new variable (%d)" %self.__used; self.__gens += [self.__gen[self.__used]]; - self.__map_of_vars[self.__gens[-1]] = current; + self.__map_of_vars[str(self.__gens[-1])] = current; self.__map_to_vars[current] = self.__used; + self.__map_of_derivatives[self.__gen[self.__used]] = self.__gen[self.__used+1]; self.__used += 1; current = current.derivative(); + + trans = self.__trans[f]; + self.__map_of_derivatives[self.__gen[before+order-1]] = sum([self.__gen[before+i]*trans[1][i][-1] for i in range(order)], trans[0]); return; diff --git a/ajpastor/lazy/conversion.py b/ajpastor/lazy/conversion.py index 2fdff3d..5a68890 100644 --- a/ajpastor/lazy/conversion.py +++ b/ajpastor/lazy/conversion.py @@ -8,6 +8,9 @@ from sage.rings.polynomial.multi_polynomial_ring import is_MPolynomialRing as is from sage.rings.polynomial.infinite_polynomial_ring import InfinitePolynomialRing_dense as isDenseIPolynomial; from sage.rings.polynomial.infinite_polynomial_ring import InfinitePolynomialRing_sparse as isSparseIPolynomial; +from sage.structure.element import is_Matrix; +from sage.structure.element import is_Vector; + class ConversionSystem(object): ## Main bulder def __init__(self, base): @@ -62,21 +65,33 @@ class ConversionSystem(object): ## Pulbic methods def add_relations(self, *relations): if(self.is_polynomial()): - try: - ## Adding the new relations - self._relations(); # We make sure the relations are initialized - - self.__add_relation(relations); - - ## Changing the ideal and computing a groebner basis - if(len(self.__relations) >= _sage_const_1 ): - self.__rel_ideal = ideal(self.poly_ring(), self.__relations); - self.__relations = self.__rel_ideal.groebner_basis(); - except TypeError as e: - raise e; + ## Adding the new relations + self._relations(); # We make sure the relations are initialized + + self.__add_relation(relations); + + ## Changing the ideal and computing a groebner basis + self.__relations = self._groebner_basis(); def clean_relations(self): self.__relations = []; + + def _add_relation(self, poly): + ''' + Auxiliar method to add a polynomial to the relations of the conversion system + ''' + reduced = self.simplify(poly); + if(reduced != _sage_const_0 ): + self.__relations += [reduced]; + + def _groebner_basis(self): + ''' + Auxiliar method to compute the groebner basis for the current set of relations on the conversion system + ''' + if(len(self.__relations) >= _sage_const_1 ): + self.__rel_ideal = ideal(self.poly_ring(), self.__relations); + return self.__rel_ideal.groebner_basis(); + return [self.poly_ring().zero()]; def to_poly(self, element): ''' @@ -100,12 +115,22 @@ class ConversionSystem(object): n = self.to_poly(element.numerator()); d = self.to_poly(element.denominator()); return self.poly_field()(n/d); - elif(isinstance(element, sage.matrix.matrix.Matrix)): + elif(isinstance(element.parent(), ConversionSystem)): + if(element.parent() is self): + return self._to_poly_element(element); + poly = element.parent().to_poly(element); + try: + n = self.to_poly(element.parent().to_real(poly.numerator())); + d = self.to_poly(element.parent().to_real(poly.denominator())); + return self.poly_field()(n/d); + except AttributeError: + return self.to_poly(element.parent().to_real(poly)); + elif(is_Matrix(element)): R = self.poly_ring(); if(element.parent().base().is_field()): R = self.poly_field(); return Matrix(R, [self.to_poly(row) for row in element]); - elif(isinstance(element, sage.modules.free_module_element.FreeModuleElement)): + elif(is_Vector(element)): R = self.poly_ring(); if(element.parent().base().is_field()): R = self.poly_field(); @@ -114,7 +139,7 @@ class ConversionSystem(object): return [self.to_poly(el) for el in element]; elif(isinstance(element, set)): return set([self.to_poly(el) for el in element]); - elif(isinstance(relation, tuple)): + elif(isinstance(element, tuple)): return tuple([self.to_poly(el) for el in element]); else: raise TypeError("Wrong type: Impossible to get polynomial value of element (%s)" %(element)); @@ -147,12 +172,12 @@ class ConversionSystem(object): n = self.to_real(poly.numerator()); d = self.to_real(poly.denominator()); return n/d; - elif(isinstance(poly, sage.matrix.matrix.Matrix)): + elif(is_Matrix(poly)): R = self.base(); if(poly.parent().base().is_field()): R = R.fraction_field(); return Matrix(R, [self.to_real(row) for row in poly]); - elif(isinstance(poly, sage.modules.free_module_element.FreeModuleElement)): + elif(is_Vector(poly)): R = self.base(); if(poly.parent().base().is_field()): R = R.fraction_field(); @@ -182,12 +207,12 @@ class ConversionSystem(object): if(element in self.poly_ring()): element = self.poly_ring()(element); try: # Weird case: fraction field fall in polynomial field - n = self.poly_ring()(element.numerator()).reduce(self._relations()); - d = self.poly_ring()(element.denominator()).reduce(self._relations()); + n = self._simplify(self.poly_ring()(element.numerator())); + d = self._simplify(self.poly_ring()(element.denominator())); return n/d; except AttributeError: try: - return self.poly_ring()(element).reduce(self._relations()); + return self._simplify(self.poly_ring()(element)); except AttributeError: return self.poly_ring()(element); elif(element in self.poly_field()): @@ -215,6 +240,15 @@ class ConversionSystem(object): return self.to_real(self.simplify(self.to_poly(element))); else: return element; + + def _simplify(self, poly): + ''' + Auxiliar method that make the simplification of an element in self.poly_ring(). + ''' + try: + return poly.reduce(self._relations()); + except AttributeError: + return poly; def mix_conversion(self, conversion): ''' @@ -266,12 +300,13 @@ class ConversionSystem(object): This method can be overwritten if needed. ''' + polynomial = self.poly_ring()(polynomial); variables = polynomial.variables(); if(len(variables) == 0): return self.base()(polynomial); - return polynomial(**{str(v) : self.map_of_vars()[str(v)] for v in variables}); + return self.base()(polynomial(**{str(v) : self.map_of_vars()[str(v)] for v in variables})); # multi = (len(variables) > _sage_const_1 ); # res = self.base().zero(); # for (k,v) in polynomial.dict().items(): @@ -311,13 +346,9 @@ class ConversionSystem(object): for el in relation: self.__add_relation(el); elif(relation in self.poly_ring()): - reduced = self.simplify(relation); - if(reduced != _sage_const_0 ): - self.__relations += [reduced]; + self._add_relation(self.poly_ring()(relation)); elif(relation in self.poly_field()): - reduced = self.simplify(relation.numerator()); - if(reduced != _sage_const_0 ): - self.__relations += [reduced]; + self._add_relation(self.poly_ring()(relation.numerator())); else: try: self.__add_relation(self.to_poly(relation)); diff --git a/ajpastor/lazy/lazyRing.py b/ajpastor/lazy/lazyRing.py index 8b0683e..9d7d4a9 100644 --- a/ajpastor/lazy/lazyRing.py +++ b/ajpastor/lazy/lazyRing.py @@ -318,7 +318,9 @@ class LazyRing (UniqueRepresentation, ConversionSystem, IntegralDomain): self.__create_poly_field(); def _to_poly_element(self, element): - if(not (element in self.base())): + if(element.parent() is self): + return element.poly(); + elif(not (element in self.base())): raise TypeError("Element is not in the base ring.\n\tExpected: %s\n\tGot: %s" %(element.parent(), self.base())); ## We try to cast to a polynomial. If we can do it, it means no further operations has to be done diff --git a/releases/diff_defined_functions.zip b/releases/diff_defined_functions.zip index 835704a..cbf4346 100644 Binary files a/releases/diff_defined_functions.zip and b/releases/diff_defined_functions.zip differ diff --git a/releases/diff_defined_functions__0.6.zip b/releases/diff_defined_functions__0.6.zip index 835704a..cbf4346 100644 Binary files a/releases/diff_defined_functions__0.6.zip and b/releases/diff_defined_functions__0.6.zip differ diff --git a/releases/old/diff_defined_functions__0.6__19.03.18_17:13:41.zip b/releases/old/diff_defined_functions__0.6__19.03.18_17:13:41.zip new file mode 100644 index 0000000..cbf4346 Binary files /dev/null and b/releases/old/diff_defined_functions__0.6__19.03.18_17:13:41.zip differ