From d70728cae36e855b82da2a299f653c4924b53616 Mon Sep 17 00:00:00 2001 From: Antonio Jimenez Pastor Date: Tue, 3 Jul 2018 11:18:29 +0200 Subject: [PATCH] Added a .gitignore so no .pyc file is tracked. Small changes in fles ddFunction.py and matrix.py --- .gitignore | 2 + ajpastor/dd_functions/ddFunction.py | 127 +++++++++++++++++++++++++++++++++++- ajpastor/misc/matrix.py | 18 +++++ 3 files changed, 145 insertions(+), 2 deletions(-) create mode 100644 .gitignore diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..fd20fdd --- /dev/null +++ b/.gitignore @@ -0,0 +1,2 @@ + +*.pyc diff --git a/ajpastor/dd_functions/ddFunction.py b/ajpastor/dd_functions/ddFunction.py index 14e13db..3cd5b17 100644 --- a/ajpastor/dd_functions/ddFunction.py +++ b/ajpastor/dd_functions/ddFunction.py @@ -2338,8 +2338,130 @@ def _m_replace(string, to_replace, escape=("(",")")): res += string[current_index:] return res; -def _get_equation_poly(poly, dic): - raise DevelopementError("_get_equation_poly"); +def _get_equation_poly(poly, dic): + ''' + Method that computes the differential equation that the evaluation + of a polynomial over some functions satisfies. + + This procedure has 5 main steps: + (1) Compute a vector space $W$ generated by a vector of generators $(g1,...,gk)$ + (2) Compute a derivation matrix $D$ for $W$ w.r.t. the generators + (3) Compute a initial vector $v_0$ that represents poly in $W$. + (4) Compute the ansatz system $S$ + (5) Solve the ansatz system + + The steps (1), (2), (3) can be done simultaneously in the method _get_derivation_matrix(poly,dic). + The step (4) use the inductive formula $v_{i+1} = Dv_i + \partial v_i$. Then $S = (v0 | v1 | ...)$. + The step (5) will be done using Bareiss Algorithm. + ''' + ### Steps (1)-(3) + D, v0 = _get_derivation_matrix(poly, dic); + + ### Step 4 + R = D[0][0].parent(); + if(sage.rings.fraction_field.is_FractionField(R)): + R = R.base(); + F = R.fraction_field(); + + v = [v0]; + for i in range(D.nrows()): + v += [D*v[-1] + vector(F, [el.derivative() for el in v[-1]])]; + + S = Matrix(F, v).transpose(); + + ### Step 5 + ## TODO We rely on the system solver of the operator of the functions we are + ## plug-in to the polynomial. + solution = dic.keys()[0].equation._get_element_nullspace(S); + + ## We create the final operator + return dic.keys()[0].parent().element([el for el in solution]).equation; + +def _get_derivation_matrix(poly, dic, simpl = False): + ''' + Method that computes the derivation matrix of a vector space that + contains the vector space generated by the evaluation of poly + using the dictionary dic and a vector that represent poly into + that vector space. + + + This method has a recursive structure: + - If for any key in dic, dic[key][i] appears in poly for i > 0, + then we simplify the polynomial so it only contains the variables + dic[key][0]. We return the matrix obtained with that polyomial and + the associated vector. + - If poly has multiple terms, we return the direct sum of the matrices + for each term and the direct sum of the vectors. + - If poly is a monomial with multiple variables, we return the Kronecker + sum of the matrices of the result and the tensor product of the vectors. + - If poly is a monomial with just one variable, we compute the + Kronecker sum of the companion matrix associated with the function + related with that variable degree times and return the vector (1,0,...,0). + + WARNING: + - The argument simpl must never be called with True but within this method + - The very base case can be improve using the fact that we are computing + the power of an element. + ''' + ## Simplification step + #if(not simpl): + # non_simple = sum([dic[key][1:] for key in dic], []); + # if(any(var in poly.variables() for var in non_simple)): + # to_eval = {} + # for key in dic: + # current = dic[key]; + # for i in range(1,len(current)): + # to_eval[str(current[i])] = current[0]; + # D, v = _get_derivation_matrix(poly(**to_eval, dic, simpl=True); + # ## TODO + # raise DevelopementError("_get_derivation_matrix"); + + ## From this point on, we can assume only simplified polynomials appears + ## Case non-monomial + if(not poly.is_monomial()): + 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))], []); + + return D,v; + ## More than one variable + elif(len(poly.variables()) > 1): + var = poly.variables()[0]; + deg = poly.degree(var); + oth = poly.parent()(poly/(var**deg)); + + D1,v1 = _get_derivation_matrix(var**deg, dic); + D2,v2 = _get_derivation_matrix(oth, dic); + + D = ajpastor.misc.Matrix.kronecker_sum(D1, D2); + v = vector(v1.parent().base(), sum(([el for el in row] for row in v1.tensor_product(v2)), [])); + + return D,v; + ## A power of a variable (using Faa di Bruno's formula) + elif(len(poly.variables()) == 1): + var = poly.variables()[0]; + deg = poly.degree(); + ## Getting the function and the offset for the initial values + func = None; + for key in dic.keys(): + if(var in dic[key]): + func = key; + offset = dic[key].index(var); + break; + ## Checking the variable can be easily represent + if(offset > func.getOrder()): + raise DevelopementError("_get_derivation_matrix"); + ## TODO Go on from here + + + raise DevelopementError("_get_derivation_matrix"); + +## TODO Saved function from the console +def get_sequences(a,b): + M = [[tensor_power(get_vector_prime(n+1),i+1) for i in range(a)] for n in range(b)] + return [[[list(M[i][j]).index(k^(j+1))+1 for k in [1] + [el for el in primes_first_n(i)]] for j in range(5)] for i in range(5)]; + __CACHED_INIT_POLY = {}; def _get_initial_poly(poly, dic, m): @@ -2396,6 +2518,7 @@ def _get_initial_poly(poly, dic, m): if(var in dic[key]): func = key; offset = dic[key].index(var); + break; init = lambda n : func.getInitialValue(offset+n); diff --git a/ajpastor/misc/matrix.py b/ajpastor/misc/matrix.py index 1e66a8d..5822b73 100644 --- a/ajpastor/misc/matrix.py +++ b/ajpastor/misc/matrix.py @@ -339,6 +339,24 @@ def simplify_matrix(M): #################################################################################### ### +### MATRIX OPERATIONS +### +### In this section we include functionality to compute some operations of +### matrices. +### +#################################################################################### +def direct_sum(*matrices): + return Matrix.block_diagonal(*matrices); + +def kronecker_sum(M,N): + if(not (M.is_square() and N.is_square())): + raise TypeError("Only square matrices for the Kronecker sum"); + + n1 = M.nrows(); n2 = N.nrows(); + return M.tensor_product(identity_matrix(n2)) + identity_matrix(n1).tensor_product(N); + +#################################################################################### +### ### MATRICIAL D-MOVE ### ### In this section we include functionality to compute a differential movement -- 2.1.4