From: Antonio Jimenez Pastor Date: Mon, 27 Aug 2018 14:49:06 +0000 (+0200) Subject: Added the file with the D^k --> DA functions. X-Git-Url: http://git.risc.jku.at/gitweb/?a=commitdiff_plain;h=6ff1fe5504da4c7ed1311486ff8bdfc43164eb72;p=ajpastor%2Fdiff_defined_functions.git Added the file with the D^k --> DA functions. Modified the ddFunction file in order to compile properly. --- diff --git a/ajpastor/dd_functions/__init__.py b/ajpastor/dd_functions/__init__.py index b793718..4bb1f3e 100644 --- a/ajpastor/dd_functions/__init__.py +++ b/ajpastor/dd_functions/__init__.py @@ -10,6 +10,10 @@ try: from .symbolic import *; except Exception: print "Error loading module dd_functions.symbolic"; +try: + from .toDiffAlgebraic import *; +except Exception: + print "Error loading module dd_functions.toDiffAlgebraic"; from pkgutil import extend_path; __path__ = extend_path(__path__, __name__); diff --git a/ajpastor/dd_functions/ddFunction.py b/ajpastor/dd_functions/ddFunction.py index 3cd5b17..470446b 100644 --- a/ajpastor/dd_functions/ddFunction.py +++ b/ajpastor/dd_functions/ddFunction.py @@ -2422,7 +2422,7 @@ def _get_derivation_matrix(poly, dic, simpl = False): coeffs = poly.coefficients(); sols = [_get_derivation_matrix(mon, dic, simpl=True) for mon in poly.monomials()]; D = Matrix.block_diagonal(*[sol[1] for sol in sols]); - v = vector(D[0][0].parent(), sum([[coeffs[i]*el for el in sols[i][2]] for i in range(len(sols))], []); + v = vector(D[0][0].parent(), sum([[coeffs[i]*el for el in sols[i][2]] for i in range(len(sols))], [])) return D,v; ## More than one variable diff --git a/ajpastor/dd_functions/toDiffAlgebraic.py b/ajpastor/dd_functions/toDiffAlgebraic.py new file mode 100644 index 0000000..bb76b27 --- /dev/null +++ b/ajpastor/dd_functions/toDiffAlgebraic.py @@ -0,0 +1,316 @@ + +from sage.all_cmdline import * # import sage library + +################################################################################ +### +### Second approach: use an infinite polynomial field +### +################################################################################ +def is_InfinitePolynomialRing(parent): + from sage.rings.polynomial.infinite_polynomial_ring import InfinitePolynomialRing_dense as dense; + from sage.rings.polynomial.infinite_polynomial_ring import InfinitePolynomialRing_sparse as sparse; + + return isinstance(parent, dense) or isinstance(parent, sparse); + +def get_InfinitePolynomialRingGen(parent, var, with_number = False): + for gen in parent.gens(): + spl = str(var).split(repr(gen)[:-1]); + if(len(spl) == 2): + if(with_number): + return (gen, int(spl[1])); + else: + return gen; + return None; + +def get_InfinitePolynomialRingVaribale(parent, gen, number): + return parent(repr(gen).replace("*", str(number))); + +def infinite_derivative(element, times=1): + if(times in ZZ and times > 1): + return infinite_derivative(infinite_derivative(element, times-1)) + if(not is_InfinitePolynomialRing(element.parent())): + try: + return element.derivative(); + except AttributeError: + return element.parent().zero(); + if(element.is_monomial()): + if(len(element.variables()) == 1): + parent = element.parent() + return parent(repr(parent.gen()).replace("*", str(int(str(element).split("_")[-1])+1))); + else: + degrees = element.degrees(); + variables = element.variables(); + part = prod(variables[i]^(degrees[i]-1) for i in range(len(degrees))); + return sum([degrees[i]*infinite_derivative(variables[i])*prod(variables[j] for j in range(len(variables)) if i != j) for i in range(len(variables))]); + else: + coefficients = element.coefficients(); + monomials = element.monomials(); + return sum([coefficients[i].derivative()*monomials[i] + coefficients[i]*infinite_derivative(monomials[i]) for i in range(len(monomials))]); + +################################################################################ +### +### Special cases for diff_to_diffalg +### +################################################################################ +## Polynomial in x case +def diff_to_diffalg_poly(func): + base, is_x = __get_base_field(func.parent()); + final_ring = PolynomialRing(base, ["x", "y0"], order='deglex'); + x = final_ring.gens()[0]; + y = final_ring.gens()[1]; + return y - final_ring(func); + +## D-finite case +def diff_to_diffalg_dfinite(func): + if(not (isinstance(func.parent(), DDRing) and func.parent().depth()==1)): + raise TypeError("diff_to_diffalg_dfinite called with a non-dfinite function"); + if(func.is_constant): + return diff_to_diffalg(func(0)); + order = func.getOrder(); + base, is_x = __get_base_field(func.parent().base()); + + name_vars = []; + if(is_x): + name_vars += ["x"]; + name_vars += ["y%d" %i for i in range(order+1)]; + final_ring = PolynomialRing(base, name_vars, order='deglex'); + if(is_x): + x = final_ring.gens()[0]; + y = final_ring.gens()[1:]; + else: + y = final_ring.gens(); + + return sum(y[i]*final_ring(func.equation[i]) for i in range(order+1)); + +## DD-finite case +def diff_to_diffalg_ddfinite(func): + if(not (isinstance(func.parent(), DDRing) and func.parent().depth()==2)): + raise TypeError("diff_to_diffalg_ddfinite called with a non-ddfinite function"); + + base, is_x = __get_base_field(func.parent().base().base()); + + polys = [diff_to_diffalg(coeff) for coeff in func.equation.getCoefficients()]; + orders = [coeff.getOrder() for coeff in func.equation.getCoefficients()]; + S = sum(orders); + name_vars = ["y%d" %i for i in range(S+func.getOrder())]; + if(is_x): + name_vars = ["x"] + name_vars; + final_ring = PolynomialRing(base, name_vars, order='deglex'); + x = None; + if(is_x): + x = final_ring.gens()[0]; + y = final_ring.gens()[1:]; + else: + y = final_ring.gens(); + first_row = []; + for i in range(len(orders)): + if(func.equation[i].getOrder() > 1 or polys[i].degree() > 1): + first_row += [[final_ring.one()] + [final_ring.zero() for j in range(1,orders[i])]]; + first_row += [final_ring.zero()]; + rows = [first_row]; + + + for i in range(S): + # Each element is r[-1][a]'+r[-1][a-1]+r[-1][-1]*eq + # Part: r[-1][a]' + new_row = [[__get_derivative(el, [],[],y,x) for el in list] for list in rows[-1][:-1]] + [__get_derivative(rows[-1][-1], [],[],y,x)]; + + # Part: r[-1][a-1] + for i in range(len(rows[:-1])): + for j in range(1, len(rows[-1][:-1][i])): + new_row[i][j] += rows[-1][i][j-1]; + + # Part: r[-1][-1]*eq + for i in range(len(new_row)-1): + for j in range(len(new_row[i])): + coeff = polys[i].coefficient({y[-1]:1}); + new_row[i][j] -= polys[i].coefficient({y[j]:1})/coeff; + new_row[-1] -= polys[i].coefficients()[-1]/coeff; + + dens = sum([[el.denominator() for el in list] for list in new_row[:-1]],[]); + dens += [new_row[-1].denominator()]; + new_lcm = lcm(dens); + fin_new_row = [[new_lcm*el for el in list] for list in new_row[:-1]]; + fin_new_row += [new_lcm*new_row[-1]]; + rows += [fin_new_row]; + + mat = [sum(rows[i][:-1], []) + [rows[i][-1]] for i in range(len(rows))]; + M = Matrix(mat); + return M.determinant(); + +@cached_function +def diff_to_diffalg(func): + ## Particular imports for this method + from sage.categories.pushout import pushout; + + ## Dificult case: diff. definable function + if(isinstance(func, DDRing.Element)): + order = func.getOrder(); + base = func.parent().base(); + + ######################################################################## + ## Recursive-case + if(func.parent().depth() > 2): + # Main recursive step + coeff_poly = [diff_to_diffalg(coeff) for coeff in func.equation.getCoefficients()]; + + # Building some data for getting the final polynomial ring + base = reduce(pushout, [poly.parent().base() for poly in coeff_poly]); + coeff_order = []; + is_x = False; + for i in range(order+1): + try: + ## Case x in poly.parent() + coeff_poly[i].parent()('x'); + coeff_order += [coeff_poly[i].parent().ngens()-2]; + is_x = True; + except: + coeff_order += [coeff_poly[i].parent().ngens()-1]; + + # Getting the variables needed for the middle step + name_vars = ["z%d_%d" %(i,j) for i in range(order+1) for j in range(coeff_order[i],-1,-1)]; + if(is_x): + name_vars += ["x"]; + + exe = 0; + eqs = []; + while(True): + max_order = 1+exe;#max(coeff_order)+exe; + + # Building the middle polynomial ring + mon_order = [TermOrder('deglex', coeff_order[i]+1) for i in range(order+1)]; + if(is_x): + mon_order += [TermOrder('deglex', 1)]; + mon_order += [TermOrder('lex', order+max_order+1)]; + mon_order = reduce(lambda p,q : p + q, mon_order); + middle_ring = PolynomialRing(base, name_vars_f, order=mon_order); + + # Recovering the variables of the polynomial ring + gens = middle_ring.gens(); + z = []; + j = 0; + for i in range(order+1): + z += [list(gens[j:j+coeff_order[i]+1])] + z[-1].reverse(); + j += coeff_order[i]+1; + if(is_x): + x = gens[j]; + j += 1; + else: + x = None; + y = list(gens[j:]); y.reverse(); + + # Computing all the required polynomials for the reduction + # We start with the equations from the coefficients + eqs = [middle_ring(el) for el in eqs]; + for i in range(len(eqs),order+1): + args = {"y%d" %j : z[i][j] for j in range(coeff_order[i]+1)}; + eqs += [coeff_poly[i](**args)]; + + # We keep computing the derivatives from the equation of func + if(len(eqs) < order+2): + eqs += [sum(y[i]*z[i][0] for i in range(order+1))]; + for i in range(len(eqs), order+2+max_order): + eqs += [__get_derivative(eqs[-1], eqs[:order+1], z,y,x)]; + + # Now we compute a Groebner-basis to eliminate all possible variables + # It is important to remark that the order of the variables were chosen + # so all possible variables are eliminated and the smallest possible + # polynomial equation remains. + gb = ideal(eqs).groebner_basis(); + + # We get the final ring for the equation + final_vars = [str(y[i]) for i in range(len(y))]; + if(is_x): + final_vars = ["x"] + final_vars; + final_ring = PolynomialRing(base, final_vars, order='deglex'); + try: + return final_ring(str(gb[-1])); + except TypeError: + print "The polynomial obtained has too many variables:\n\t- Expected: %s\n\t- Got: %s" %(final_ring.gens(), gb[-1]); + + exe += 1; + + ######################################################################## + ## Basic cases + elif(func.parent().depth() == 2): + return diff_to_diffalg_ddfinite(func); + else: + return diff_to_diffalg_dfinite(func); + + ############################################################################ + ## Base case: polynomial or non-diff_definable functions + return diff_to_diffalg_poly(func); + +def __get_base_field(base): + ## Particular imports for this method + from sage.rings.polynomial.polynomial_ring import is_PolynomialRing; + from sage.misc.functional import is_field; + + ## Checking the case is R(x) + if(is_field(base) and not (base.base() == base)): + return __get_base_field(base.base()); + + ## Checking the case is R[x] + is_x = False; + if(is_PolynomialRing(base) and str(base.gens()[0]) == 'x'): + base = base.base(); + is_x = True; + + if(not base.is_field()): + base = base.fraction_field(); + + return base, is_x; + +__CACHE_OF_DICS = {}; + +def __get_derivative(poly, equations, z,y,x): + if(poly.is_constant()): + return 0; + global __CACHE_OF_DICS; + key = tuple([tuple(equations), tuple([tuple(row) for row in z]), tuple(y)]); + dic_of_derivatives = __CACHE_OF_DICS.get(key, {}); + + coeffs = poly.coefficients(); monomials = poly.monomials(); + gens = poly.parent().gens(); + der_mon = []; + for mon in monomials: + degrees = mon.degrees(); + basic = prod(gens[i]**max(0,degrees[i]-1) for i in range(len(gens))); + extra = poly.parent().zero(); + for i in range(len(gens)): + if(degrees[i] > 0): + extra += degrees[i]*__get_variable_derivative(dic_of_derivatives, gens[i], equations, z, y ,x)*prod(gens[j] for j in range(len(gens)) if (j != i and degrees[j] > 0)); + der_mon += [basic*extra]; + + __CACHE_OF_DICS[key] = dic_of_derivatives; + + element = sum(coeffs[i]*der_mon[i] for i in range(len(coeffs))); + if(hasattr(element, "numerator")): + element = element.numerator(); + + return element; + +def __get_variable_derivative(cache, var, equations, z,y,x): + if(var not in cache): + if(var == x): + cache[var] = 1; + if(var in y): + cache[var] = y[y.index(var)+1]; + else: + for i in range(len(z)): + if(var in z[i]): + j = z[i].index(var); + if(j < len(z[i])-1): + cache[var] = z[i][j+1]; + else: + # Computing derivative of z[i][-1] + d = equations[i].degree(z[i][j]); + alpha = [equations[i].coefficient({z[i][j]: k}) for k in range(d+1)]; + cache[var] = -sum(__get_derivative(alpha[i], equations,z,y,x)*z[i][-1]^i for i in range(d+1))/sum(alpha[i]*i*z[i][-1]^(i-1) for i in range(1,d+1)); + return cache[var]; + +################################################################################################### +### PACKAGE ENVIRONMENT VARIABLES +################################################################################################### +__all__ = ["is_InfinitePolynomialRing", "get_InfinitePolynomialRingGen", "get_InfinitePolynomialRingVaribale", "infinite_derivative"]; diff --git a/releases/diff_defined_functions__0.5.zip b/releases/diff_defined_functions__0.5.zip index 246fc51..e7b8dd5 100644 Binary files a/releases/diff_defined_functions__0.5.zip and b/releases/diff_defined_functions__0.5.zip differ diff --git a/releases/old/diff_defined_functions__0.5__18.08.27_16:49:05.zip b/releases/old/diff_defined_functions__0.5__18.08.27_16:49:05.zip new file mode 100644 index 0000000..e7b8dd5 Binary files /dev/null and b/releases/old/diff_defined_functions__0.5__18.08.27_16:49:05.zip differ