From: Antonio Jimenez Pastor Date: Tue, 29 Jan 2019 08:55:17 +0000 (+0100) Subject: Started the addition for a Lazy Conversion System for DDRings X-Git-Url: http://git.risc.jku.at/gitweb/?a=commitdiff_plain;h=318fac54d2a454f6e5c3072e313bc3e7a9cee554;p=ajpastor%2Fdiff_defined_functions.git Started the addition for a Lazy Conversion System for DDRings using all the ideas for simplification got in the diffalg_reduction method. --- diff --git a/ajpastor/dd_functions/ddExamples.py b/ajpastor/dd_functions/ddExamples.py index 1ea4f9b..3e5b8f6 100644 --- a/ajpastor/dd_functions/ddExamples.py +++ b/ajpastor/dd_functions/ddExamples.py @@ -1149,8 +1149,8 @@ def GenericHypergeometricFunction(num=[],den=[],init=_sage_const_1 ): if(not((numerator,denominator,initial) in __CACHED_HYPERGEOMETRIC)): ## Building differential operator get_op = lambda p : destiny_ring.element(p).equation; - op_num = x*prod(get_op([el,x]) for el in numerator); - op_den = x*get_op([_sage_const_0 ,_sage_const_1 ])*prod(get_op([el-_sage_const_1 ,x]) for el in denominator); + op_num = prod(get_op([el,x]) for el in numerator,x); + op_den = prod(get_op([el-_sage_const_1 ,x]) for el in denominator, get_op([_sage_const_0 ,x])); op = op_num - op_den; f = destiny_ring.element(op); diff --git a/ajpastor/dd_functions/lazyDDRing.py b/ajpastor/dd_functions/lazyDDRing.py new file mode 100644 index 0000000..be70a16 --- /dev/null +++ b/ajpastor/dd_functions/lazyDDRing.py @@ -0,0 +1,653 @@ + +# This file was *autogenerated* from the file ./lazyRing.sage +from sage.all_cmdline import * # import sage library + +_sage_const_2 = Integer(2); _sage_const_1 = Integer(1); _sage_const_0 = Integer(0); _sage_const_20 = Integer(20) +from sage.rings.ring import IntegralDomain; +from sage.structure.element import IntegralDomainElement; +from sage.categories.integral_domains import IntegralDomains; +from sage.categories.fields import Fields; +from sage.categories.pushout import ConstructionFunctor; + +from sage.rings.polynomial.polynomial_ring import is_PolynomialRing as isUniPolynomial; +from sage.rings.polynomial.multi_polynomial_ring import is_MPolynomialRing as isMPolynomial; + +from ajpastor.lazy.conversion import ConversionSystem; +from ajpastor.dd_functions.ddFunction import *; + +#################################################################################################### +#################################################################################################### +### ELEMENT CLASS +#################################################################################################### +#################################################################################################### +class _LazyDDFunction(IntegralDomainElement): + def __init__(self, parent, el): + if(not (isinstance(parent, LazyDDRing))): + raise TypeError("The parent for a _LazyDDFunction is not a LazyDDRing"); + + self.__raw = None; + self.__poly = None; + + IntegralDomainElement.__init__(self, parent); + + try: + self.__poly = parent.poly_ring()(el); + except: + try: + self.__poly = parent.poly_field()(str(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 + + ################################################################################################ + ### Methods for a LazyElement + ################################################################################################ + def raw(self): + ''' + Method that computes (if needed) and returns an element of `self.base()` that is equal to `self`. + ''' + if(self.__raw is None): + self.__raw = self.parent().to_real(self.__poly); + + return self.__raw; + + def poly(self): + ''' + Method that computes (if needed) and returns an polynomial such that the conversion using + self.parent() returns self.raw(). + ''' + if(self.__poly is None): + self.__poly = self.parent().to_poly(self.__raw); + + return self.__poly; + + def simplify(self): + self.__poly = self.parent().simplify(self.poly()); + + def variables(self): + ''' + Method that returns a tuple with the variables that appear in self.poly(). + + If such polynomial representation is a quotient of polynomials, it take the union of the variables in the numerator and denominator. + If such polynomial representation has no variables, returns an empty tuple. + ''' + if(self.poly() not in QQ): + if(self.poly() in self.parent().poly_ring()): + return self.parent().poly_ring()(self.poly()).variables(); + else: + var_n = self.parent().poly_ring()(self.poly().numerator()).variables(); + var_d = self.parent().poly_ring()(self.poly().denominator()).variables(); + + return tuple(set(var_n+var_d)); + return tuple(); + + def derivative(self, times = 1): + ''' + Method that computes the derivative of an element in the laziest way possible. + + This method compute the derivative of each of the variables in 'self.poly()' and then build + a new polynomial using this expressions. + + This method rely on the parent method 'get_derivative()' which saves the derivatives of the + variables as a dictionary + ''' + return self.parent().derivative(self, times); + + ################################################################################################ + ### Arithmetics methods + ################################################################################################ + def _add_(self, other): + if(isinstance(other, _LazyDDFunction)): + try: + return _LazyDDFunction(self.parent(), self.poly()+self.parent()(other).poly()); + except NotImplementedError as e: + pass; + + return NotImplemented; + + def _sub_(self,other): + return self.__add__(-other); + + def _neg_(self): + try: + return _LazyDDFunction(self.parent(), -self.poly()); + except NotImplementedError: + pass; + + return NotImplemented; + + def _mul_(self,other): + if(isinstance(other, _LazyDDFunction)): + try: + return _LazyDDFunction(self.parent(), self.poly()*self.parent()(other).poly()); + except NotImplementedError as e: + pass; + + return NotImplemented; + + def _div_(self,other): + if(isinstance(other, _LazyDDFunction)): + try: + return _LazyDDFunction(self.parent(), self.parent().poly_field()(self.poly())/self.parent().poly_field()(self.parent()(other).poly())); + except NotImplementedError as e: + pass; + + return NotImplemented; + + def _pow_(self,i): + try: + return _LazyDDFunction(self.parent(), self.poly()**i); + except NotImplementedError: + return NotImplemented; + + def __eq__(self, other): + return (self-other).is_zero(); + + def __call__(self, *input, **kwds): + if(self.__raw is None): + try: + return self.poly()(**{str(var):self.parent().to_real(var)(*input,**kwds) for var in self.variables()}); + except TypeError: + pass; + return self.raw()(*input,**kwds); + + ################################################################################################ + ### Non-trivial arithmetics methods + ################################################################################################ + def gcd(self,*input): + ''' + Method that a common divisor of 'self' and the input + ''' + if(len(input) > _sage_const_1): + return self.gcd(input); + + return _LazyDDFunction(self.parent(), gcd([self.poly()]+[self.parent()(el).poly() for el in input])); + + def lcm(self,*input): + ''' + Method that a common multiple of 'self' and the input + ''' + if(len(input) > _sage_const_1 ): + return self.gcd(input); + + return _LazyDDFunction(self.parent(), lcm([self.poly()]+[self.parent()(el).poly() for el in input])); + + def divides(self, other): + ''' + Method that returns True if 'other = a*self'. + + REMARK: If this methods return False does not mean we can not divide other by self in the level of 'base'. + ''' + return self.poly().divides(self.parent()(other).poly()); + + @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 )): + if(not (self.__raw is None)): + result = (self.raw() == 0); + else: + pol = None; + if(self.poly() in self.parent().poly_ring()): + pol = self.parent().poly_ring()(self.poly()); + else: + pol = self.parent().poly_ring()(self.poly().numerator()); + + factors = None; + try: + try: + factors = pol.factor(proof=True); + except NotImplementedError: + factors = pol.factor(proof=False); + except: + factors = [(pol,_sage_const_1 )]; + for factor in factors: + result = (self.parent().to_real(factor[_sage_const_0 ]) == _sage_const_0 ); + if(result): + self.parent().add_relations(factor[_sage_const_0 ]); + break; + + return result; + + @cached_method + def is_one(self): + # Computing self-1 and checking if it is zero + result = (self-1).is_zero(); + + # Doing a great simplification in __poly + if(result): + self.__poly = self.parent().poly_ring().one(); + return result; + + ################################################################################################ + ### Representation methods + ################################################################################################ + def __repr__(self): + return repr(self.poly()); + + def __str__(self): + if(self.__raw is None): + return "Lazy Element: %s" %(repr(self)); + else: + return "Lazy Element: %s\n%s" %(repr(self), str(self.raw())); + +#################################################################################################### +#################################################################################################### +### RING CLASS +#################################################################################################### +#################################################################################################### +class LazyDDRing (UniqueRepresentation, ConversionSystem, IntegralDomain): + + Element = _LazyDDFunction; + + def __init__(self, base, constants=QQ, var_name="z" category=None): + ''' + Implementation of a Covnersion System using InfinitePolynomialRing as a basic structure. + + Elements of 'base' will be represented in this ring using elements of + 'InfinitePolynomialRing(constants, names=[var_name])' and linear relations of the type + 'f = cg+d' where 'f in base', 'g in self.gens()' and 'c,d in constants'. + + This class captures all the functionality to work as a ConversionSystem and all the + data structures required for ompting lazily with those objects. + ''' + ## Checking the arguments + if(not (constants in Fields())): + raise TypeError("The argument 'constants' must be a Field"); + if(not (isinstance(base, Ring_w_Sequence))): + raise TypeError("The argument 'base' must be a Ring with Sequence"); + + ## Initializing the parent structures + ConversionSystem.__init__(self, base); + IntegralDomain.__init__(self, base, category); + + ## Initializing the attributes + self.__constants = constants; + + self.__poly_ring = InfinitePolynomialRing(constants, names=[var_name]); + self.__gen = self.__poly_ring.gens()[0]; + self.__used = 0; + + ## Initializing the map of variables + self.__map_of_vars = {}; # Map where each variable is associated with a 'base' object + self.__map_to_vars = {}; # Map where each variable is associated with a 'base' object + self.__r_graph = DiGraph(); # Graph of relations between the already casted elements + self.__trans = dict(); # Dictionary with the transformation to get into a node in 'r_graph' + + self.__gens = []; + + self.__map_of_derivatives = {}; # Map for each variable to its derivative (as a polynomial) + + ## Casting and Coercion system + self.base().register_conversion(LDRSimpleMorphism(self, self.base())); + + ################################################################################################ + ### Implementing methods of ConversionSystem + ################################################################################################ + def poly_ring(self): + return self.__poly_ring; + + def poly_field(self): + return self.__poly_field; + + def map_of_vars(self): + return self.__map_of_vars; + + def variables(self): + 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())): + 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 + try: + return self.poly_ring()(element); + except: + pass; + + ## Now we check if the element has a built method + try: + built = element.built(); + if(not (built is None)): + if(built[_sage_const_0 ] == "derivative"): + if(not(element in built[_sage_const_1 ])): + integral = self(built[_sage_const_1 ][_sage_const_0 ]); + if(not (integral.poly().is_monomial() and integral.poly().degree() == _sage_const_1 )): + return self(built[_sage_const_1 ][_sage_const_0 ]).derivative().poly(); + elif(built[_sage_const_0 ] == "polynomial"): + ## We check we do not have infinite recursion + if(not element in built[_sage_const_1 ][_sage_const_1 ].values()): + ## We have some building information + polys = {key : self.to_poly(built[_sage_const_1 ][_sage_const_1 ][key]) for key in built[_sage_const_1 ][_sage_const_1 ]}; + return built[_sage_const_1 ][_sage_const_0 ](**polys); + except AttributeError: + pass; + + ## Otherwise we look for a linear relation between the element and the variables + ## Breath-search for a linear relation + roots = [el for el in self.__r_graph.vertices() if self.__r_graph.in_degree() == 0];# Getting roots + from collections import deque; + to_search = deque(roots); # Initializing the queue + while(len(to_search) > 0): + # Updating the queue + current = to_search.popleft(); + for edge in self.__r_graph.outgoing_edges(current): + to_search.append(edge[1]); + + # Visiting the current node + relation = self.__find_relation(element, current); + if(not (relation is None)): + ## If it is not just the function, create a new node for future relations + if(relation[0] != 0): + self.__add_function(element); + self.__r_graph.add_edge((current, element, relation)); + + ## Creating the vector in the current level + v = [0 for i in range(self.__trans[current][1].nrows())]; + v[relation[0]] = relation[1]; + v = vector(self, v); + + ## Going to a root with the vector + v,dest = self.__pullup_vector(v, current); + + ## Returning the appropriate polynomial using the final vector + return self.__vector_to_poly(v, dest); + + ## At this point, no direct relation was found. We look for relation between element and the roots + for root in roots: + relation = self.__find_relation(root, element); + if(not (relation is None)): + # Adding the element to the graph + self.__add_function(element, gen=True); + self.__r_graph.add_edge((element, root, relation)); + + return self.__map_to_vars[element]; + + ## Otherwise, we add the element to the variables and return it + self.__add_function(element, gen=True); + return self.__map_to_vars[element]; + + ################################################################################################ + ### Other Methods for LazyRing + ################################################################################################ + def sequence(self, el, n): + if(not isinstance(el, _LazyDDFunction)): + return el.parent().sequence(el,n); + + return self.base().sequence(el.raw(),n); + + def clean_ring(self): + ## Cleaning the relations + self.clean_relations(); + + ## Cleaning all the variables for the laziness + self.__used = 0; + self.__map_of_vars = {}; + self.__map_to_vars = {}; + self.__r_graph = DiGraph(); + self.__trans = dict(); + self.__gens = []; + self.__map_of_derivatives = {}; + + def get_derivative(self, el): + if(el in self.__map_of_derivatives): + return self.__map_of_derivatives[el]; + + if(el in self.__constants): + return _sage_const_0 ; + else: + try: + el = self.poly_ring()(el); + if(el in self.poly_ring().gens()): + new_poly = self.to_poly(self.map_of_vars()[str(el)].derivative()); + self.__map_of_derivatives[el] = new_poly; + return new_poly; + except: + pass; + raise ValueError("Method 'get_derivative' can only be summoned with variables.\nGot: %s -- %s" %(el, el.parent())); + + ################################################################################################ + ### Other Integral Domain methods + ################################################################################################ + def base_ring(self): + return self.base().base_ring(); + + def characteristic(self): + return self.base().characteristic(); + + def _an_element_(self): + return self.one(); + + def element(self, X): + return self(X); + + ################################################################################################ + ### Other Ring methods + ################################################################################################ + def is_field(self): + return self.base().is_field(); + + def is_finite(self): + return self.base().is_finite(); + + def is_integrally_closed(self): + return self.base().is_integrally_closed(); + + def is_noetherian(self): + return self.base().is_noetherian(); + + ################################################################################################ + ### Coercion methods + ################################################################################################ + def _coerce_map_from_(self, S): + coer = None; + if(isinstance(S, LazyRing)): + coer = self.base()._coerce_map_from_(S.base()); + elif(S == self.base()): + coer = True; + else: + coer = self.base()._coerce_map_from_(S); + + if(not(coer is False) and not(coer is None)): + return True; + return None; + + def _has_coerce_map_from(self, S): + coer = self._coerce_map_from_(S); + return (not(coer is False) and not(coer is None)); + + def _element_constructor_(self, *args, **kwds): + if(len(args) < _sage_const_1 ): + print args + raise ValueError("Impossible to build an element without arguments"); + + i = _sage_const_0 ; + if(len(args) >= _sage_const_2 ): + if(not (args[_sage_const_0 ] is self)): + raise ValueError("RIOKO: What the .... are you sending to this method?"); + i = _sage_const_1 ; + X = args[i]; + + try: + if(not isinstance(X, _LazyDDFunction)): + ## If the element is not a LazyElement, then we try to create a new element with it + return _LazyDDFunction(self, X); + elif (X.parent() is self): + return X; + else: + ## Otherwise, X.parent() may have different variables + other = X.parent(); + pol = X.poly(); + + ## 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]); + + ## We now plugin the expressions + return _LazyDDFunction(self, pol(**translate)); + except TypeError: + raise TypeError("This element can not be casted to %s" %repr(self)); + + def construction(self): + return (LazyDDRingFunctor(), self.base()); + + def __contains__(self, X): + try: + return (X.parent() is self) or (self._has_coerce_map_from(X.parent())); + except AttributeError: + try: + self(X) + return True; + except Exception: + return False; + + ################################################################################################ + ### Representation methods + ################################################################################################ + def __repr__(self): + return "Lazy DD-Ring over (%s)" %(repr(self.base())); + + 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 += "}"; + return final + + ################################################################################################ + ### Private methods + ################################################################################################ + def __add_function(self, f, gen=False): + ## We look for relations with the derivatives of itself + derivatives = [f.derivative(times=i) for i in range(f.getOrder())]; + + order = f.getOrder(); + relation = None; + for i in range(f.getOrder(),0,-1): + n_relation = self.__find_relation(derivatives[i], derivatives[:i-1], _debug); + if(not (n_relation is None)): + order = i-1; + relation = n_relation; + else: + break; + + ## Create the tranformation into f + 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]): + new_row[-1] = relation[1][0]; + C += [[companion[i][j] for j in range(k-1)] + []]; + C = Matrix(self,C); + self.__trans[f] = (relation[1],C); + else: + self.__trans[f] = (self.zero(),Matrix(self, f.equation.companion())); + + ## Adding f as a vertex (edges must be added outside this function) + self.__r_graph.add_vertex(f); + + ## If f will become a generator, add it to the used variables + if(gen): + current = f; + for i in range(order): + self.__gens += [self.__gen[self.__used]]; + self.__map_of_vars[self.__gens[-1]] = current; + self.__map_to_vars[current] = self.__gens[-1]; + self.__used += 1; + current = current.derivative(); + + return; + + def __find_relation(self, g,f): + ''' + Method that get the relation between g and a f. If possible, it computes the derivatives up to + some order of f and check relations with them. + + It returns a tuple (k,res) where: + - k: the first index which we found a relation + - res: the relation (see __find_linear_relation to see how such relation is expressed). + + Then the following statement is always True: + g == f[k]*res[0] + res[1] + ''' + if(is_DDFunction(f)): + f = [f.derivative(times=i) for i in range(f.getOrder())]; + + for k in range(d): + res = __find_linear_relation(g,f[k]); + if(not (res is None)): + return (k,res); + + return None; + +def __find_linear_relation(f,g): + ''' + This method receives two DDFunctions and return two constants (c,d) such that f = cg+d. + None is return if those constants do not exist. + ''' + if(not (is_DDFunction(f) and is_DDFunction(g))): + raise TypeError("find_linear_relation: The functions are not DDFunctions"); + + ## Simplest case: some of the functions are constants is a constant + if(f.is_constant): + return (f.parent().zero(), f(0)); + elif(g.is_constant): + return None; + + try: + of = 1; og = 1; + while(f.getSequenceElement(of) == 0): of += 1; + while(g.getSequenceElement(og) == 0): og += 1; + + if(of == og): + c = f.getSequenceElement(of)/g.getSequenceElement(og); + d = f(0) - c*g(0); + + if(f == c*g + d): + return (c,d); + except Exception: + pass; + + return None; + +#################################################################################################### +#################################################################################################### +### Construction Functor for LazyRing +#################################################################################################### +#################################################################################################### +class LazyDDRingFunctor (ConstructionFunctor): + def __init__(self): + ID = IntegralDomains(); + self.rank = _sage_const_20 ; + super(LazyDDRingFunctor, self).__init__(ID,ID); + + ### Methods to implement + def _coerce_into_domain(self, x): + if(x not in self.domain()): + raise TypeError("The object [%s] is not an element of [%s]" %(x, self.domain())); + if(isinstance(x, LazyDDRing)): + return x.base(); + return x; + + def _apply_functor(self, x): + return LazyDDRing(x); + +#################################################################################################### +#################################################################################################### +### General Morphism for return to basic rings +#################################################################################################### +#################################################################################################### +class LDRSimpleMorphism (sage.categories.map.Map): + def __init__(self, domain, codomain): + super(LDRSimpleMorphism, self).__init__(domain, codomain); + + def _call_(self, p): + return self.codomain()(p.raw()); + diff --git a/ajpastor/dd_functions/toDiffAlgebraic.py b/ajpastor/dd_functions/toDiffAlgebraic.py index 98e7267..7ecfdbf 100644 --- a/ajpastor/dd_functions/toDiffAlgebraic.py +++ b/ajpastor/dd_functions/toDiffAlgebraic.py @@ -35,9 +35,6 @@ def get_InfinitePolynomialRingGen(parent, var, with_number = False): return gen; return None; -def get_InfinitePolynomialRingVaribale(parent, gen, number): - return parent(repr(gen).replace("*", str(number))); - #### TODO: Review this function def infinite_derivative(element, times=1, var=x, _infinite=True): ''' @@ -95,7 +92,7 @@ def infinite_derivative(element, times=1, var=x, _infinite=True): ## Case of degree 1 (add one to the index of the variable) if(len(element.variables()) == 1 and element.degree() == 1): g,n = get_InfinitePolynomialRingGen(parent, element, True); - return r_method(get_InfinitePolynomialRingVaribale(parent, g,n+1)); + return r_method(g[n+1]); ## Case of higher degree else: # Computing the variables in the monomial and their degrees (in order) @@ -208,8 +205,8 @@ def diffalg_reduction(poly, _infinite=False, _debug=False): maximal_elements = __digraph_roots(graph, _debug); ## Detecting case: all coefficients are already simpler - if(len(maximal_elements) == 1 and maximal_elements[0] == -1): - __dprint("Detected simplicity: returning the same polynomial over the basic ring"); + if(len(maximal_elements) == 1 and maximal_elements[0] == -1 and len(graph.outgoing_edges(-1)) == len(graph.edges())): + __dprint("Detected simplicity: returning the same polynomial over the basic ring", _debug); return r_method(InfinitePolynomialRing(R, names=["y"])(poly)); # Reducing the considered coefficients @@ -308,24 +305,29 @@ def __simplify_coefficients(base, coeffs, derivatives, _debug=False): R = []; # Checking the simplest cases (coeff[i] in R) + simpler = []; for i in range(n): if(coeffs[i] in base): if((-1) not in E): E += [-1]; R += [(-1,i,(0,0,coeffs[i]))]; + simpler += [i]; + __dprint("Simpler coefficiets (related with -1): %s" %simpler, _debug); # Starting the comparison for i in range(n): for j in range(i+1, n): # Checking (i,j) - rel = __find_relation(coeffs[j], derivatives[i], _debug); - if(not(rel is None)): - R += [(i,j,rel)]; - continue; + if(not(j in simpler)): + rel = __find_relation(coeffs[j], derivatives[i], _debug); + if(not(rel is None)): + R += [(i,j,rel)]; + continue; # Checking (j,i) - rel = __find_relation(coeffs[i], derivatives[j], _debug); - if(not(rel is None)): - R += [(j,i,rel)]; + if(not(i in simpler)): + rel = __find_relation(coeffs[i], derivatives[j], _debug); + if(not(rel is None)): + R += [(j,i,rel)]; # Checking if the basic case will be used because of the relations if((-1) not in E and any(e[2][1][1] != 0 for e in R)): @@ -350,12 +352,13 @@ def __simplify_derivatives(derivatives, _debug=False): ''' res = []; for i in range(len(derivatives)): + to_add = (None,None); for k in range(len(derivatives[i])-1,0,-1): relation = __find_relation(derivatives[i][k], derivatives[i][:k-1], _debug); if(not (relation is None)): - res += [(k,relation)]; + to_add = (k,relation); break; - res += [(None,None)]; + res += [to_add]; __dprint("Found relations with derivatives: %s" %res, _debug); return res; @@ -371,11 +374,14 @@ def __find_relation(g,df, _debug=False): Then the following statement is always True: g == df[k]*res[0] + res[1] ''' + __dprint("** Checking relations between %s and %s" %(repr(g), df), _debug); for k in range(len(df)): - res = __find_linear_relation(df[k], g); + res = __find_linear_relation(g,df[k]); if(not (res is None)): __dprint("Found relation:\n\t(%s) == [%s]*(%s) + [%s]" %(repr(g),repr(res[0]), repr(df[k]), repr(res[1])), _debug); + __dprint("*************************", _debug); return (k,res); + __dprint("Nothing found\n*************************", _debug); return None; @@ -387,6 +393,12 @@ def __find_linear_relation(f,g): if(not (is_DDFunction(f) and is_DDFunction(g))): raise TypeError("find_linear_relation: The functions are not DDFunctions"); + ## Simplest case: some of the functions are constants is a constant + if(f.is_constant): + return (f.parent().zero(), f(0)); + elif(g.is_constant): + return None; + try: of = 1; og = 1; while(f.getSequenceElement(of) == 0): of += 1; @@ -403,12 +415,6 @@ def __find_linear_relation(f,g): return None; - ## Simplest case: some of the functions are constants is a constant - if(f.is_constant): - return (f.parent().zero(), f(0)); - elif(g.is_constant): - return None; - def __min_span_tree(digraph, _debug=False): ''' Method that computes a forest from a directed graph containing all vertices and having the minimal depth. @@ -547,11 +553,11 @@ def __build_vector(coeffs, monomials, graph, drelations, cR, _debug=False): const_val = 0; # Doing a Tree Transversal in POstorder in the graph to pull up the vectors - stack = copy(maximal_elements); stack.reverse(); # Case for -1 if(constant): maximal_elements = [-1] + maximal_elements; - const_val = sum(cR(coeffs[e[1]])*monomials[coeffs[e[1]]] for e in graph.outgoing_edges(-1)); + # const_val = sum(cR(coeffs[e[1]])*monomials[coeffs[e[1]]] for e in graph.outgoing_edges(-1)); + stack = copy(maximal_elements); stack.reverse(); # Rest of the graph ready = []; @@ -563,31 +569,42 @@ def __build_vector(coeffs, monomials, graph, drelations, cR, _debug=False): # Visit the node: pull-up the vector edge = graph.incoming_edges(current)[0]; cu_vector = vectors[current]; - de_vector = vectors[edge[0]]; - relation = edge[2]; - C = trans[edge[0]][1]; - - __dprint("** Reducing node %d to %d" %(current, edge[0]), _debug); - __dprint("\t- Current vector: %s" %(cu_vector), _debug); - __dprint("\t- Destiny vector: %s" %(de_vector), _debug); - __dprint("\t- Relation: %s" %(str(relation)), _debug); - __dprint("\t- Matrix:\n\t\t%s" %(str(C).replace('\n','\n\t\t')), _debug); - __dprint("\t- Prev. constant: %s" %const_val, _debug); - - # Building vectors for all the required derivatives - ivectors = [vector(cR, [0 for i in range(relation[0])] + [1] + [0 for i in range(relation[0]+1,coeff_order[edge[0]])])]; - extra_cons = [relation[1][1]]; - - for i in range(len(cu_vector)): - extra_cons += [ivectors[-1][-1]*trans[edge[0]][0]]; - ivectors += [differential_movement(C, ivectors[-1], infinite_derivative)]; - vectors[edge[0]] = sum([vector(cR,[cu_vector[i]*ivectors[i][j] for j in range(len(ivectors[i]))]) for i in range(len(cu_vector))], de_vector); - const_val += sum([cu_vector[i]*extra_cons[i] for i in range(len(cu_vector))]); + ## Constant case; reducing to the constant term + if(edge[0] == -1): + __dprint("** Reducing node %d to constant" %(current), _debug); + __dprint("\t- Current vector: %s" %(cu_vector), _debug); + __dprint("\t- Prev. constant: %s" %const_val, _debug); + const_val += sum(cu_vector[i]*cR(coeffs[current].derivative(times=i)) for i in range(len(cu_vector))); + __dprint("\t- New constant: %s" %const_val, _debug); - __dprint("\t- New vector: %s" %vectors[edge[0]], _debug); - __dprint("\t- New constant: %s" %const_val, _debug); - __dprint("*************************", _debug); + ## General case + else: + de_vector = vectors[edge[0]]; + relation = edge[2]; + C = trans[edge[0]][1]; + + __dprint("** Reducing node %d to %d" %(current, edge[0]), _debug); + __dprint("\t- Current vector: %s" %(cu_vector), _debug); + __dprint("\t- Destiny vector: %s" %(de_vector), _debug); + __dprint("\t- Relation: %s" %(str(relation)), _debug); + __dprint("\t- Matrix:\n\t\t%s" %(str(C).replace('\n','\n\t\t')), _debug); + __dprint("\t- Prev. constant: %s" %const_val, _debug); + + # Building vectors for all the required derivatives + ivectors = [vector(cR, [0 for i in range(relation[0])] + [1] + [0 for i in range(relation[0]+1,coeff_order[edge[0]])])]; + extra_cons = [relation[1][1]]; + + for i in range(len(cu_vector)): + extra_cons += [ivectors[-1][-1]*trans[edge[0]][0]]; + ivectors += [differential_movement(C, ivectors[-1], infinite_derivative)]; + + vectors[edge[0]] = sum([vector(cR,[cu_vector[i]*ivectors[i][j] for j in range(len(ivectors[i]))]) for i in range(len(cu_vector))], de_vector); + const_val += sum([cu_vector[i]*extra_cons[i] for i in range(len(cu_vector))]); + + __dprint("\t- New vector: %s" %vectors[edge[0]], _debug); + __dprint("\t- New constant: %s" %const_val, _debug); + __dprint("*************************", _debug); # Getting out the element of the stack stack.pop(); @@ -1106,14 +1123,14 @@ def __dprint(obj, _debug): #################################################################################################### #### PACKAGE ENVIRONMENT VARIABLES #################################################################################################### -__all__ = ["is_InfinitePolynomialRing", "get_InfinitePolynomialRingGen", "get_InfinitePolynomialRingVaribale", "infinite_derivative", "diffalg_reduction", "toDifferentiallyAlgebraic_Below", "diff_to_diffalg", "inverse_DA", "func_inverse_DA", "guess_DA_DDfinite", "guess_homogeneous_DNfinite", "FaaDiBruno_polynomials", "Exponential_polynomials"]; +__all__ = ["is_InfinitePolynomialRing", "get_InfinitePolynomialRingGen", "infinite_derivative", "diffalg_reduction", "toDifferentiallyAlgebraic_Below", "diff_to_diffalg", "inverse_DA", "func_inverse_DA", "guess_DA_DDfinite", "guess_homogeneous_DNfinite", "FaaDiBruno_polynomials", "Exponential_polynomials"]; # Extra functions for debugging def simplify_coefficients(base, coeffs, derivatives, _debug=True): return __simplify_coefficients(base,coeffs,derivatives,_debug); def simplify_derivatives(derivatives, _debug=True): - return __simplify_derivatives(derivatives, _debug) + return __simplify_derivatives(derivatives, _debug); def find_relation(g,df, _debug=True): return __find_relation(g,df,_debug); diff --git a/releases/diff_defined_functions__0.6.zip b/releases/diff_defined_functions__0.6.zip index bff5f6b..eecea87 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.01.29_09:55:17.zip b/releases/old/diff_defined_functions__0.6__19.01.29_09:55:17.zip new file mode 100644 index 0000000..eecea87 Binary files /dev/null and b/releases/old/diff_defined_functions__0.6__19.01.29_09:55:17.zip differ