--- /dev/null
+
+# 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());
+
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):
'''
## 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)
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
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)):
'''
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;
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;
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;
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.
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 = [];
# 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();
####################################################################################################
#### 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);
projects / ajpastor / diff_defined_functions.git / commitdiff