From d17f413ca195b9680f89f678b77de899765bcaac Mon Sep 17 00:00:00 2001 From: Hirish Chandrasekaran Date: Sun, 2 Jun 2024 17:27:30 -0700 Subject: [PATCH 01/20] Create recurrence.py - Feedback requested --- sumpy/recurrence.py | 296 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 296 insertions(+) create mode 100644 sumpy/recurrence.py diff --git a/sumpy/recurrence.py b/sumpy/recurrence.py new file mode 100644 index 000000000..1bf0519ff --- /dev/null +++ b/sumpy/recurrence.py @@ -0,0 +1,296 @@ +__copyright__ = """ +Copyright (C) 2024 Hirish Chandrasekaran +Copyright (C) 2024 Andreas Kloeckner +""" + +__license__ = """ +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in +all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +THE SOFTWARE. +""" + +from collections import namedtuple +from pyrsistent import pmap +from pytools import memoize +from sumpy.tools import add_mi +from itertools import accumulate +import sumpy.symbolic as sym +import logging +from typing import List +import sympy as sp +from sumpy.expansion.diff_op import LinearPDESystemOperator +from pytools.obj_array import make_obj_array + +#A similar function exists in sumpy.symbolic +def make_sympy_vec(name, n): + return make_obj_array([sp.Symbol(f"{name}{i}") for i in range(n)]) + + +__doc__ = """ +.. autoclass:: Recurrence +.. automodule:: sumpy.recurrence +""" + +#CREATE LAPLACE_3D +DerivativeIdentifier = namedtuple("DerivativeIdentifier", ["mi", "vec_idx"]) +partial2_x = DerivativeIdentifier((2,0,0), 0) +partial2_y = DerivativeIdentifier((0,2,0), 0) +partial2_z = DerivativeIdentifier((0,0,2), 0) +#Coefficients +list_pde_dict_3d = {partial2_x: 1, partial2_y: 1, partial2_z: 1} +laplace_3d = LinearPDESystemOperator(3,list_pde_dict_3d) + +#CREATE LAPLACE_2D +partial2_x = DerivativeIdentifier((2,0), 0) +partial2_y = DerivativeIdentifier((0,2), 0) +#Coefficients +list_pde_dict = {partial2_x: 1, partial2_y: 1} +laplace_2d = LinearPDESystemOperator(2,list_pde_dict) + +#CREATE HELMHOLTZ_2D +func_val = DerivativeIdentifier((0,0), 0) +#Coefficients +list_pde_dict = {partial2_x: 1, partial2_y: 1, func_val: 1} +helmholtz_2d = LinearPDESystemOperator(2,list_pde_dict) + +''' +get_pde_in_recurrence_form +Input: + - pde, a :class:`sumpy.expansion.diff_op.LinearSystemPDEOperator` pde such that assert(len(pde.eqs) == 1) + is true. +Output: + - ode_in_r, an ode in r which the POINT-POTENTIAL (has radial symmetry) satisfies away from the origin. + Note: to represent f, f_r, f_{rr}, we use the sympy variables f_{r0}, f_{r1}, .... So ode_in_r is a linear + combination of the sympy variables f_{r0}, f_{r1}, .... + - var, represents the variables for the input space: [x0, x1, ...] + - n_derivs, the order of the original PDE + 1, i.e. the number of derivatives of f that may be present + (the reason this is called n_derivs since if we have a second order PDE for example + then we might see f, f_{r}, f_{rr} in our ODE in r, which is technically 3 terms since we count + the 0th order derivative f as a "derivative." If this doesn't make sense just know that n_derivs + is the order the of the input sumpy PDE + 1) + +Description: We assume we are handed a system of 1 sumpy PDE (pde) and output the +pde in a way that allows us to easily replace derivatives with respect to r. In other words we output +a linear combination of sympy variables f_{r0}, f_{r1}, ... (which represents f, f_r, f_{rr} respectively) +to represent our ODE in r for the point potential. +''' +def get_pde_in_recurrence_form(laplace): + dim = laplace.dim + n_derivs = laplace.order + assert(len(laplace.eqs) == 1) + ops = len(laplace.eqs[0]) + derivs = [] + coeffs = [] + for i in laplace.eqs[0]: + derivs.append(i.mi) + coeffs.append(laplace.eqs[0][i]) + var = make_sympy_vec("x", dim) + r = sp.sqrt(sum(var**2)) + + eps = sp.symbols("epsilon") + rval = r + eps + + f = sp.Function("f") + f_derivs = [sp.diff(f(rval),eps,i) for i in range(n_derivs+1)] + + def compute_term(a, t): + term = a + for i in range(len(t)): + term = term.diff(var[i], t[i]) + return term + + pde = 0 + for i in range(ops): + pde += coeffs[i] * compute_term(f(rval), derivs[i]) + + n_derivs = len(f_derivs) + f_r_derivs = make_sympy_vec("f_r", n_derivs) + + for i in range(n_derivs): + pde = pde.subs(f_derivs[i], f_r_derivs[i]) + + return pde, var, n_derivs + + +ode_in_r, var, n_derivs = get_pde_in_recurrence_form(laplace_2d) + + +''' +generate_ND_derivative_relations +Input: + - var, a sympy vector of variables called [x0, x1, ...] + - n_derivs, the order of the original PDE + 1, i.e. the number of derivatives of f that may be present +Output: + - a vector that gives [f, f_r, f_{rr}, ...] in terms of f, f_x, f_{xx}, ... using the chain rule + (f, f_x, f_{xx}, ... in code is represented as f_{x0}, f_{x1}, f_{x2} and + f, f_r, f_{rr}, ... in code is represented as f_{r0}, f_{r1}, f_{r2}) + +Description: Using the chain rule outputs a vector that tells us how to write f, f_r, f_{rr}, ... as a linear +combination of f, f_x, f_{xx}, ... +''' +def generate_ND_derivative_relations(var, n_derivs): + f_r_derivs = make_sympy_vec("f_r", n_derivs) + f_x_derivs = make_sympy_vec("f_x", n_derivs) + f = sp.Function("f") + eps = sp.symbols("epsilon") + rval = sp.sqrt(sum(var**2)) + eps + f_derivs_x = [sp.diff(f(rval),var[0],i) for i in range(n_derivs)] + f_derivs = [sp.diff(f(rval),eps,i) for i in range(n_derivs)] + for i in range(len(f_derivs_x)): + for j in range(len(f_derivs)): + f_derivs_x[i] = f_derivs_x[i].subs(f_derivs[j], f_r_derivs[j]) + system = [f_x_derivs[i] - f_derivs_x[i] for i in range(n_derivs)] + + return sp.solve(system, *f_r_derivs, dict=True)[0] + + +''' +ode_in_r_to_x +Input: + - ode_in_r, a linear combination of f, f_r, f_{rr}, ... (in code represented as f_{r0}, f_{r1}, f_{r2}) + with coefficients as RATIONAL functions in var[0], var[1], ... + - var, array of sympy variables [x_0, x_1, ...] + - n_derivs, the order of the original PDE + 1, i.e. the number of derivatives of f that may be present +Output: + - ode_in_x, a linear combination of f, f_x, f_{xx}, ... with coefficients as rational + functions in var[0], var[1], ... + +Description: Translates an ode in the variable r into an ode in the variable x by substituting f, f_r, f_{rr}, ... + as a linear combination of f, f_x, f_{xx}, ... using the chain rule +''' +def ode_in_r_to_x(ode_in_r, var, n_derivs): + subme = generate_ND_derivative_relations(var, n_derivs) + ode_in_x = ode_in_r + f_r_derivs = make_sympy_vec("f_r", n_derivs) + for i in range(n_derivs): + ode_in_x = ode_in_x.subs(f_r_derivs[i], subme[f_r_derivs[i]]) + return ode_in_x + + +ode_in_x = ode_in_r_to_x(ode_in_r, var, n_derivs).simplify() +ode_in_x_cleared = (ode_in_x * var[0]**n_derivs).simplify() + +delta_x = sp.symbols("delta_x") +c_vec = make_sympy_vec("c", len(var)) + +''' +compute_poly_in_deriv +Input: + - ode_in_x_cleared, an ode in x, i.e. a linear combination of f, f_x, f_{xx}, ... + (in code represented as f_{x0}, f_{x1}, f_{x2}) with coefficients as POLYNOMIALS in var[0], var[1], ... + (i.e. not rational functions) + - n_derivs, the order of the original PDE + 1, i.e. the number of derivatives of f that may be present +Output: + - a polynomial in f, f_x, f_{xx}, ... (in code represented as f_{x0}, f_{x1}, f_{x2}) with coefficients + as polynomials in \delta_x where \delta_x = x_0 - c_0 that represents the ''shifted ODE'' - i.e. the ODE + where we substitute all occurences of \delta_x with x_0 - c_0 + +Description: Converts an ode in x, to a polynomial in f, f_x, f_{xx}, ..., with coefficients as polynomials +in \delta_x = x_0 - c_0. +''' +def compute_poly_in_deriv(ode_in_x_cleared, n_derivs): + #Note that generate_ND_derivative_relations will at worst put some power of $x_0^order$ in the denominator. To clear + #the denominator we can probably? just multiply by x_0^order. + ode_in_x_cleared = (ode_in_x * var[0]**n_derivs).simplify() + + ode_in_x_shifted = ode_in_x_cleared.subs(var[0], delta_x + c_vec[0]).simplify() + + f_x_derivs = make_sympy_vec("f_x", n_derivs) + poly = sp.Poly(ode_in_x_shifted, *f_x_derivs) + + return poly + +poly = compute_poly_in_deriv(ode_in_x, n_derivs) + +''' +compute_coefficients_of_poly +Input: + - poly, a polynomial in sympy variables f_{x0}, f_{x1}, ..., + (recall that this corresponds to f_0, f_x, f_{xx}, ...) with coefficients that are polynomials in \delta_x + where poly represents the ''shifted ODE''- i.e. we substitute all occurences of \delta_x with x_0 - c_0 +Output: + - a 2d array, each row giving the coefficient of f_0, f_x, f_{xx}, ..., + each entry in the row giving the coefficients of the polynomial in \delta_x + +Description: Takes in a polynomial in f_{x0}, f_{x1}, ..., w/coeffs that are polynomials in \delta_x +and outputs a 2d array for easy access to the coefficients based on their degree as a polynomial in \delta_x. +''' +def compute_coefficients_of_poly(poly, n_derivs): + #Returns coefficients in lexographic order. So lowest order first + def tup(i,n=n_derivs): + a = [] + for j in range(n): + if j != i: + a.append(0) + else: + a.append(1) + return tuple(a) + + coeffs = [] + for deriv_ind in range(n_derivs): + coeffs.append(sp.Poly(poly.coeff_monomial(tup(deriv_ind)), delta_x).all_coeffs()) + + return coeffs + +coeffs = compute_coefficients_of_poly(poly, n_derivs) + +i = sp.symbols("i") +s = sp.Function("s") + +''' +compute_recurrence_relation +Input: + - coeffs a 2d array that gives access to the coefficients of poly, where poly represents the coefficients of + the ''shifted ODE'' (''shifted ODE'' = we substitute all occurences of \delta_x with x_0 - c_0) + based on their degree as a polynomial in \delta_x + - n_derivs, the order of the original PDE + 1, i.e. the number of derivatives of f that may be present +Output: + - a recurrence statement that equals 0 where s(i) is the ith coefficient of the Taylor polynomial for + our point potential. + +Description: Takes in coeffs which represents our ``shifted ode in x" (i.e. ode_in_x with coefficients in \delta_x) +and outputs a recurrence relation for the point potential. +''' + +def compute_recurrence_relation(coeffs, n_derivs): + #Compute symbolic derivative + def hc_diff(i, n): + retMe = 1 + for j in range(n): + retMe *= (i-j) + return retMe + + #We are differentiating deriv_ind, which shifts down deriv_ind. Do this for one deriv_ind + r = 0 + for deriv_ind in range(n_derivs): + part_of_r = 0 + pow_delta = 0 + for j in range(len(coeffs[deriv_ind])-1, -1, -1): + shift = pow_delta - deriv_ind + 1 + pow_delta += 1 + temp = coeffs[deriv_ind][j] * s(i) * hc_diff(i, deriv_ind) + part_of_r += temp.subs(i, i-shift) + r += part_of_r + + for j in range(1, len(var)): + r = r.subs(var[j], c_vec[j]) + + return r.simplify() + +r = compute_recurrence_relation(coeffs, n_derivs) + + From 073426354b1c1a235840428b78360ff08be293c4 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Mon, 3 Jun 2024 16:07:17 -0500 Subject: [PATCH 02/20] Hackin with Andreas --- sumpy/recurrence.py | 111 +++++++++++++++++++------------------------- 1 file changed, 48 insertions(+), 63 deletions(-) diff --git a/sumpy/recurrence.py b/sumpy/recurrence.py index 1bf0519ff..7586cffa9 100644 --- a/sumpy/recurrence.py +++ b/sumpy/recurrence.py @@ -24,13 +24,6 @@ """ from collections import namedtuple -from pyrsistent import pmap -from pytools import memoize -from sumpy.tools import add_mi -from itertools import accumulate -import sumpy.symbolic as sym -import logging -from typing import List import sympy as sp from sumpy.expansion.diff_op import LinearPDESystemOperator from pytools.obj_array import make_obj_array @@ -45,50 +38,29 @@ def make_sympy_vec(name, n): .. automodule:: sumpy.recurrence """ -#CREATE LAPLACE_3D -DerivativeIdentifier = namedtuple("DerivativeIdentifier", ["mi", "vec_idx"]) -partial2_x = DerivativeIdentifier((2,0,0), 0) -partial2_y = DerivativeIdentifier((0,2,0), 0) -partial2_z = DerivativeIdentifier((0,0,2), 0) -#Coefficients -list_pde_dict_3d = {partial2_x: 1, partial2_y: 1, partial2_z: 1} -laplace_3d = LinearPDESystemOperator(3,list_pde_dict_3d) - -#CREATE LAPLACE_2D -partial2_x = DerivativeIdentifier((2,0), 0) -partial2_y = DerivativeIdentifier((0,2), 0) -#Coefficients -list_pde_dict = {partial2_x: 1, partial2_y: 1} -laplace_2d = LinearPDESystemOperator(2,list_pde_dict) - -#CREATE HELMHOLTZ_2D -func_val = DerivativeIdentifier((0,0), 0) -#Coefficients -list_pde_dict = {partial2_x: 1, partial2_y: 1, func_val: 1} -helmholtz_2d = LinearPDESystemOperator(2,list_pde_dict) -''' -get_pde_in_recurrence_form -Input: - - pde, a :class:`sumpy.expansion.diff_op.LinearSystemPDEOperator` pde such that assert(len(pde.eqs) == 1) - is true. -Output: - - ode_in_r, an ode in r which the POINT-POTENTIAL (has radial symmetry) satisfies away from the origin. - Note: to represent f, f_r, f_{rr}, we use the sympy variables f_{r0}, f_{r1}, .... So ode_in_r is a linear - combination of the sympy variables f_{r0}, f_{r1}, .... - - var, represents the variables for the input space: [x0, x1, ...] - - n_derivs, the order of the original PDE + 1, i.e. the number of derivatives of f that may be present - (the reason this is called n_derivs since if we have a second order PDE for example - then we might see f, f_{r}, f_{rr} in our ODE in r, which is technically 3 terms since we count - the 0th order derivative f as a "derivative." If this doesn't make sense just know that n_derivs - is the order the of the input sumpy PDE + 1) - -Description: We assume we are handed a system of 1 sumpy PDE (pde) and output the -pde in a way that allows us to easily replace derivatives with respect to r. In other words we output -a linear combination of sympy variables f_{r0}, f_{r1}, ... (which represents f, f_r, f_{rr} respectively) -to represent our ODE in r for the point potential. -''' def get_pde_in_recurrence_form(laplace): + ''' + get_pde_in_recurrence_form + Input: + - pde, a :class:`sumpy.expansion.diff_op.LinearSystemPDEOperator` pde such that assert(len(pde.eqs) == 1) + is true. + Output: + - ode_in_r, an ode in r which the POINT-POTENTIAL (has radial symmetry) satisfies away from the origin. + Note: to represent f, f_r, f_{rr}, we use the sympy variables f_{r0}, f_{r1}, .... So ode_in_r is a linear + combination of the sympy variables f_{r0}, f_{r1}, .... + - var, represents the variables for the input space: [x0, x1, ...] + - n_derivs, the order of the original PDE + 1, i.e. the number of derivatives of f that may be present + (the reason this is called n_derivs since if we have a second order PDE for example + then we might see f, f_{r}, f_{rr} in our ODE in r, which is technically 3 terms since we count + the 0th order derivative f as a "derivative." If this doesn't make sense just know that n_derivs + is the order the of the input sumpy PDE + 1) + + Description: We assume we are handed a system of 1 sumpy PDE (pde) and output the + pde in a way that allows us to easily replace derivatives with respect to r. In other words we output + a linear combination of sympy variables f_{r0}, f_{r1}, ... (which represents f, f_r, f_{rr} respectively) + to represent our ODE in r for the point potential. + ''' dim = laplace.dim n_derivs = laplace.order assert(len(laplace.eqs) == 1) @@ -113,20 +85,33 @@ def compute_term(a, t): term = term.diff(var[i], t[i]) return term - pde = 0 + ode_in_r = 0 for i in range(ops): - pde += coeffs[i] * compute_term(f(rval), derivs[i]) + ode_in_r += coeffs[i] * compute_term(f(rval), derivs[i]) n_derivs = len(f_derivs) f_r_derivs = make_sympy_vec("f_r", n_derivs) for i in range(n_derivs): - pde = pde.subs(f_derivs[i], f_r_derivs[i]) + ode_in_r = ode_in_r.subs(f_derivs[i], f_r_derivs[i]) - return pde, var, n_derivs + return ode_in_r, var, n_derivs + + +def test_recurrence_finder(): + from sumpy.expansion.diff_op import make_identity_diff_op, laplacian + w = make_identity_diff_op(2) + laplace2d = laplacian(w) + print(get_pde_in_recurrence_form(laplace2d)) + + assert 1 == 1 + + + + -ode_in_r, var, n_derivs = get_pde_in_recurrence_form(laplace_2d) +# ode_in_r, var, n_derivs = get_pde_in_recurrence_form(laplace_2d) ''' @@ -181,11 +166,11 @@ def ode_in_r_to_x(ode_in_r, var, n_derivs): return ode_in_x -ode_in_x = ode_in_r_to_x(ode_in_r, var, n_derivs).simplify() -ode_in_x_cleared = (ode_in_x * var[0]**n_derivs).simplify() +# ode_in_x = ode_in_r_to_x(ode_in_r, var, n_derivs).simplify() +# ode_in_x_cleared = (ode_in_x * var[0]**n_derivs).simplify() -delta_x = sp.symbols("delta_x") -c_vec = make_sympy_vec("c", len(var)) +# delta_x = sp.symbols("delta_x") +# c_vec = make_sympy_vec("c", len(var)) ''' compute_poly_in_deriv @@ -214,7 +199,7 @@ def compute_poly_in_deriv(ode_in_x_cleared, n_derivs): return poly -poly = compute_poly_in_deriv(ode_in_x, n_derivs) +# poly = compute_poly_in_deriv(ode_in_x, n_derivs) ''' compute_coefficients_of_poly @@ -246,10 +231,10 @@ def tup(i,n=n_derivs): return coeffs -coeffs = compute_coefficients_of_poly(poly, n_derivs) +# coeffs = compute_coefficients_of_poly(poly, n_derivs) -i = sp.symbols("i") -s = sp.Function("s") +# i = sp.symbols("i") +# s = sp.Function("s") ''' compute_recurrence_relation @@ -291,6 +276,6 @@ def hc_diff(i, n): return r.simplify() -r = compute_recurrence_relation(coeffs, n_derivs) +# r = compute_recurrence_relation(coeffs, n_derivs) From 7e26f2895080a9963bdd14affe54b8a3c42546a9 Mon Sep 17 00:00:00 2001 From: Hirish Chandrasekaran Date: Mon, 3 Jun 2024 20:01:16 -0700 Subject: [PATCH 03/20] Fix all flake8 issues --- sumpy/recurrence.py | 270 ++++++++++++++++++++++---------------------- 1 file changed, 133 insertions(+), 137 deletions(-) diff --git a/sumpy/recurrence.py b/sumpy/recurrence.py index 7586cffa9..227ada25f 100644 --- a/sumpy/recurrence.py +++ b/sumpy/recurrence.py @@ -23,12 +23,11 @@ THE SOFTWARE. """ -from collections import namedtuple import sympy as sp -from sumpy.expansion.diff_op import LinearPDESystemOperator from pytools.obj_array import make_obj_array -#A similar function exists in sumpy.symbolic + +#A similar function exists in sumpy.symbolic def make_sympy_vec(name, n): return make_obj_array([sp.Symbol(f"{name}{i}") for i in range(n)]) @@ -42,28 +41,33 @@ def make_sympy_vec(name, n): def get_pde_in_recurrence_form(laplace): ''' get_pde_in_recurrence_form - Input: - - pde, a :class:`sumpy.expansion.diff_op.LinearSystemPDEOperator` pde such that assert(len(pde.eqs) == 1) + Input: + - pde, a :class:`sumpy.expansion.diff_op.LinearSystemPDEOperator` pde such + that assert(len(pde.eqs) == 1) is true. - Output: - - ode_in_r, an ode in r which the POINT-POTENTIAL (has radial symmetry) satisfies away from the origin. - Note: to represent f, f_r, f_{rr}, we use the sympy variables f_{r0}, f_{r1}, .... So ode_in_r is a linear - combination of the sympy variables f_{r0}, f_{r1}, .... + Output: + - ode_in_r, an ode in r which the POINT-POTENTIAL (has radial symmetry) + satisfies away from the origin. + Note: to represent f, f_r, f_{rr}, we use the sympy variables + f_{r0}, f_{r1}, .... So ode_in_r is a linear combination of the sympy + variables f_{r0}, f_{r1}, .... - var, represents the variables for the input space: [x0, x1, ...] - - n_derivs, the order of the original PDE + 1, i.e. the number of derivatives of f that may be present - (the reason this is called n_derivs since if we have a second order PDE for example - then we might see f, f_{r}, f_{rr} in our ODE in r, which is technically 3 terms since we count - the 0th order derivative f as a "derivative." If this doesn't make sense just know that n_derivs - is the order the of the input sumpy PDE + 1) - - Description: We assume we are handed a system of 1 sumpy PDE (pde) and output the - pde in a way that allows us to easily replace derivatives with respect to r. In other words we output - a linear combination of sympy variables f_{r0}, f_{r1}, ... (which represents f, f_r, f_{rr} respectively) + - n_derivs, the order of the original PDE + 1, i.e. the number of + derivatives of f that may be present (the reason this is called n_derivs + since if we have a second order PDE for example then we might see f, f_{r}, + f_{rr} in our ODE in r, which is technically 3 terms since we count + the 0th order derivative f as a "derivative." If this doesn't make sense + just know that n_derivs is the order the of the input sumpy PDE + 1) + + Description: We assume we are handed a system of 1 sumpy PDE (pde) and output + the pde in a way that allows us to easily replace derivatives with respect to r. + In other words we output a linear combination of sympy variables + f_{r0}, f_{r1}, ... (which represents f, f_r, f_{rr} respectively) to represent our ODE in r for the point potential. ''' dim = laplace.dim n_derivs = laplace.order - assert(len(laplace.eqs) == 1) + assert (len(laplace.eqs) == 1) ops = len(laplace.eqs[0]) derivs = [] coeffs = [] @@ -75,10 +79,9 @@ def get_pde_in_recurrence_form(laplace): eps = sp.symbols("epsilon") rval = r + eps - f = sp.Function("f") - f_derivs = [sp.diff(f(rval),eps,i) for i in range(n_derivs+1)] - + f_derivs = [sp.diff(f(rval), eps, i) for i in range(n_derivs+1)] + def compute_term(a, t): term = a for i in range(len(t)): @@ -88,76 +91,63 @@ def compute_term(a, t): ode_in_r = 0 for i in range(ops): ode_in_r += coeffs[i] * compute_term(f(rval), derivs[i]) - n_derivs = len(f_derivs) f_r_derivs = make_sympy_vec("f_r", n_derivs) for i in range(n_derivs): ode_in_r = ode_in_r.subs(f_derivs[i], f_r_derivs[i]) - return ode_in_r, var, n_derivs -def test_recurrence_finder(): - from sumpy.expansion.diff_op import make_identity_diff_op, laplacian - w = make_identity_diff_op(2) - laplace2d = laplacian(w) - print(get_pde_in_recurrence_form(laplace2d)) - - assert 1 == 1 - - - - - - -# ode_in_r, var, n_derivs = get_pde_in_recurrence_form(laplace_2d) - - -''' -generate_ND_derivative_relations -Input: +def generate_ND_derivative_relations(var, n_derivs): + ''' + generate_ND_derivative_relations + Input: - var, a sympy vector of variables called [x0, x1, ...] - - n_derivs, the order of the original PDE + 1, i.e. the number of derivatives of f that may be present -Output: - - a vector that gives [f, f_r, f_{rr}, ...] in terms of f, f_x, f_{xx}, ... using the chain rule + - n_derivs, the order of the original PDE + 1, i.e. the number of derivatives of + f that may be present + Output: + - a vector that gives [f, f_r, f_{rr}, ...] in terms of f, f_x, f_{xx}, ... + using the chain rule (f, f_x, f_{xx}, ... in code is represented as f_{x0}, f_{x1}, f_{x2} and f, f_r, f_{rr}, ... in code is represented as f_{r0}, f_{r1}, f_{r2}) -Description: Using the chain rule outputs a vector that tells us how to write f, f_r, f_{rr}, ... as a linear -combination of f, f_x, f_{xx}, ... -''' -def generate_ND_derivative_relations(var, n_derivs): + Description: Using the chain rule outputs a vector that tells us how to + write f, f_r, f_{rr}, ... as a linear + combination of f, f_x, f_{xx}, ... + ''' f_r_derivs = make_sympy_vec("f_r", n_derivs) f_x_derivs = make_sympy_vec("f_x", n_derivs) f = sp.Function("f") eps = sp.symbols("epsilon") rval = sp.sqrt(sum(var**2)) + eps - f_derivs_x = [sp.diff(f(rval),var[0],i) for i in range(n_derivs)] - f_derivs = [sp.diff(f(rval),eps,i) for i in range(n_derivs)] + f_derivs_x = [sp.diff(f(rval), var[0], i) for i in range(n_derivs)] + f_derivs = [sp.diff(f(rval), eps, i) for i in range(n_derivs)] for i in range(len(f_derivs_x)): for j in range(len(f_derivs)): f_derivs_x[i] = f_derivs_x[i].subs(f_derivs[j], f_r_derivs[j]) system = [f_x_derivs[i] - f_derivs_x[i] for i in range(n_derivs)] - return sp.solve(system, *f_r_derivs, dict=True)[0] -''' -ode_in_r_to_x -Input: - - ode_in_r, a linear combination of f, f_r, f_{rr}, ... (in code represented as f_{r0}, f_{r1}, f_{r2}) +def ode_in_r_to_x(ode_in_r, var, n_derivs): + ''' + ode_in_r_to_x + Input: + - ode_in_r, a linear combination of f, f_r, f_{rr}, ... + (in code represented as f_{r0}, f_{r1}, f_{r2}) with coefficients as RATIONAL functions in var[0], var[1], ... - var, array of sympy variables [x_0, x_1, ...] - - n_derivs, the order of the original PDE + 1, i.e. the number of derivatives of f that may be present -Output: - - ode_in_x, a linear combination of f, f_x, f_{xx}, ... with coefficients as rational - functions in var[0], var[1], ... - -Description: Translates an ode in the variable r into an ode in the variable x by substituting f, f_r, f_{rr}, ... - as a linear combination of f, f_x, f_{xx}, ... using the chain rule -''' -def ode_in_r_to_x(ode_in_r, var, n_derivs): + - n_derivs, the order of the original PDE + 1, i.e. the number of derivatives of + f that may be present + Output: + - ode_in_x, a linear combination of f, f_x, f_{xx}, ... with coefficients as + rational functions in var[0], var[1], ... + + Description: Translates an ode in the variable r into an ode in the variable x + by substituting f, f_r, f_{rr}, ... as a linear combination of + f, f_x, f_{xx}, ... using the chain rule + ''' subme = generate_ND_derivative_relations(var, n_derivs) ode_in_x = ode_in_r f_r_derivs = make_sympy_vec("f_r", n_derivs) @@ -166,57 +156,54 @@ def ode_in_r_to_x(ode_in_r, var, n_derivs): return ode_in_x -# ode_in_x = ode_in_r_to_x(ode_in_r, var, n_derivs).simplify() -# ode_in_x_cleared = (ode_in_x * var[0]**n_derivs).simplify() - -# delta_x = sp.symbols("delta_x") -# c_vec = make_sympy_vec("c", len(var)) - -''' -compute_poly_in_deriv -Input: - - ode_in_x_cleared, an ode in x, i.e. a linear combination of f, f_x, f_{xx}, ... - (in code represented as f_{x0}, f_{x1}, f_{x2}) with coefficients as POLYNOMIALS in var[0], var[1], ... - (i.e. not rational functions) - - n_derivs, the order of the original PDE + 1, i.e. the number of derivatives of f that may be present -Output: - - a polynomial in f, f_x, f_{xx}, ... (in code represented as f_{x0}, f_{x1}, f_{x2}) with coefficients - as polynomials in \delta_x where \delta_x = x_0 - c_0 that represents the ''shifted ODE'' - i.e. the ODE - where we substitute all occurences of \delta_x with x_0 - c_0 - -Description: Converts an ode in x, to a polynomial in f, f_x, f_{xx}, ..., with coefficients as polynomials -in \delta_x = x_0 - c_0. -''' -def compute_poly_in_deriv(ode_in_x_cleared, n_derivs): - #Note that generate_ND_derivative_relations will at worst put some power of $x_0^order$ in the denominator. To clear +def compute_poly_in_deriv(ode_in_x, n_derivs, var): + ''' + compute_poly_in_deriv + Input: + - ode_in_x, a linear combination of f, f_x, f_{xx}, ... with coefficients as + rational functions in var[0], var[1], ... + - n_derivs, the order of the original PDE + 1, i.e. the number of derivatives + of f that may be present + Output: + - a polynomial in f, f_x, f_{xx}, ... (in code represented as f_{x0}, f_{x1}, + f_{x2}) with coefficients as polynomials in delta_x where delta_x = x_0 - c_0 + that represents the ''shifted ODE'' - i.e. the ODE where we substitute all + occurences of delta_x with x_0 - c_0 + + Description: Converts an ode in x, to a polynomial in f, f_x, f_{xx}, ..., + with coefficients as polynomials in delta_x = x_0 - c_0. + ''' + #Note that generate_ND_derivative_relations will at worst put some power of + #$x_0^order$ in the denominator. To clear #the denominator we can probably? just multiply by x_0^order. + delta_x = sp.symbols("delta_x") + c_vec = make_sympy_vec("c", len(var)) ode_in_x_cleared = (ode_in_x * var[0]**n_derivs).simplify() - ode_in_x_shifted = ode_in_x_cleared.subs(var[0], delta_x + c_vec[0]).simplify() - f_x_derivs = make_sympy_vec("f_x", n_derivs) poly = sp.Poly(ode_in_x_shifted, *f_x_derivs) - return poly -# poly = compute_poly_in_deriv(ode_in_x, n_derivs) - -''' -compute_coefficients_of_poly -Input: - - poly, a polynomial in sympy variables f_{x0}, f_{x1}, ..., - (recall that this corresponds to f_0, f_x, f_{xx}, ...) with coefficients that are polynomials in \delta_x - where poly represents the ''shifted ODE''- i.e. we substitute all occurences of \delta_x with x_0 - c_0 -Output: - - a 2d array, each row giving the coefficient of f_0, f_x, f_{xx}, ..., - each entry in the row giving the coefficients of the polynomial in \delta_x - -Description: Takes in a polynomial in f_{x0}, f_{x1}, ..., w/coeffs that are polynomials in \delta_x -and outputs a 2d array for easy access to the coefficients based on their degree as a polynomial in \delta_x. -''' + def compute_coefficients_of_poly(poly, n_derivs): + ''' + compute_coefficients_of_poly + Input: + - poly, a polynomial in sympy variables f_{x0}, f_{x1}, ..., + (recall that this corresponds to f_0, f_x, f_{xx}, ...) with coefficients + that are polynomials in delta_x where poly represents the ''shifted ODE'' + - i.e. we substitute all occurences of delta_x with x_0 - c_0 + Output: + - a 2d array, each row giving the coefficient of f_0, f_x, f_{xx}, ..., + each entry in the row giving the coefficients of the polynomial in delta_x + Description: Takes in a polynomial in f_{x0}, f_{x1}, ..., w/coeffs that are + polynomials in delta_x and outputs a 2d array for easy access to the + coefficients based on their degree as a polynomial in delta_x. + ''' + delta_x = sp.symbols("delta_x") + #Returns coefficients in lexographic order. So lowest order first - def tup(i,n=n_derivs): + def tup(i, n=n_derivs): a = [] for j in range(n): if j != i: @@ -224,42 +211,46 @@ def tup(i,n=n_derivs): else: a.append(1) return tuple(a) - + coeffs = [] for deriv_ind in range(n_derivs): - coeffs.append(sp.Poly(poly.coeff_monomial(tup(deriv_ind)), delta_x).all_coeffs()) - + coeffs.append( + sp.Poly(poly.coeff_monomial(tup(deriv_ind)), delta_x).all_coeffs()) + return coeffs -# coeffs = compute_coefficients_of_poly(poly, n_derivs) - -# i = sp.symbols("i") -# s = sp.Function("s") - -''' -compute_recurrence_relation -Input: - - coeffs a 2d array that gives access to the coefficients of poly, where poly represents the coefficients of - the ''shifted ODE'' (''shifted ODE'' = we substitute all occurences of \delta_x with x_0 - c_0) - based on their degree as a polynomial in \delta_x - - n_derivs, the order of the original PDE + 1, i.e. the number of derivatives of f that may be present -Output: - - a recurrence statement that equals 0 where s(i) is the ith coefficient of the Taylor polynomial for - our point potential. - -Description: Takes in coeffs which represents our ``shifted ode in x" (i.e. ode_in_x with coefficients in \delta_x) -and outputs a recurrence relation for the point potential. -''' - -def compute_recurrence_relation(coeffs, n_derivs): + +def compute_recurrence_relation(coeffs, n_derivs, var): + ''' + compute_recurrence_relation + Input: + - coeffs a 2d array that gives access to the coefficients of poly, where poly + represents the coefficients of the ''shifted ODE'' + (''shifted ODE'' = we substitute all occurences of delta_x with x_0 - c_0) + based on their degree as a polynomial in delta_x + - n_derivs, the order of the original PDE + 1, i.e. the number of derivatives + of f that may be present + Output: + - a recurrence statement that equals 0 where s(i) is the ith coefficient of + the Taylor polynomial for our point potential. + + Description: Takes in coeffs which represents our ``shifted ode in x" + (i.e. ode_in_x with coefficients in delta_x) and outputs a recurrence relation + for the point potential. + ''' + i = sp.symbols("i") + s = sp.Function("s") + c_vec = make_sympy_vec("c", len(var)) + #Compute symbolic derivative def hc_diff(i, n): retMe = 1 for j in range(n): retMe *= (i-j) return retMe - - #We are differentiating deriv_ind, which shifts down deriv_ind. Do this for one deriv_ind + + #We are differentiating deriv_ind, which shifts down deriv_ind. + #Do this for one deriv_ind r = 0 for deriv_ind in range(n_derivs): part_of_r = 0 @@ -270,12 +261,17 @@ def hc_diff(i, n): temp = coeffs[deriv_ind][j] * s(i) * hc_diff(i, deriv_ind) part_of_r += temp.subs(i, i-shift) r += part_of_r - + for j in range(1, len(var)): r = r.subs(var[j], c_vec[j]) - + return r.simplify() -# r = compute_recurrence_relation(coeffs, n_derivs) +def test_recurrence_finder(): + from sumpy.expansion.diff_op import make_identity_diff_op, laplacian + w = make_identity_diff_op(2) + laplace2d = laplacian(w) + print(get_pde_in_recurrence_form(laplace2d)) + assert 1 == 1 From f40dbe583924dcf1a6855fded63a7c4eb0636b69 Mon Sep 17 00:00:00 2001 From: Hirish Chandrasekaran Date: Tue, 4 Jun 2024 10:14:27 -0700 Subject: [PATCH 04/20] Flake8 Issues --- sumpy/recurrence.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/sumpy/recurrence.py b/sumpy/recurrence.py index 227ada25f..defcebf29 100644 --- a/sumpy/recurrence.py +++ b/sumpy/recurrence.py @@ -99,9 +99,9 @@ def compute_term(a, t): return ode_in_r, var, n_derivs -def generate_ND_derivative_relations(var, n_derivs): +def generate_nd_derivative_relations(var, n_derivs): ''' - generate_ND_derivative_relations + generate_nd_derivative_relations Input: - var, a sympy vector of variables called [x0, x1, ...] - n_derivs, the order of the original PDE + 1, i.e. the number of derivatives of @@ -148,7 +148,7 @@ def ode_in_r_to_x(ode_in_r, var, n_derivs): by substituting f, f_r, f_{rr}, ... as a linear combination of f, f_x, f_{xx}, ... using the chain rule ''' - subme = generate_ND_derivative_relations(var, n_derivs) + subme = generate_nd_derivative_relations(var, n_derivs) ode_in_x = ode_in_r f_r_derivs = make_sympy_vec("f_r", n_derivs) for i in range(n_derivs): @@ -173,7 +173,7 @@ def compute_poly_in_deriv(ode_in_x, n_derivs, var): Description: Converts an ode in x, to a polynomial in f, f_x, f_{xx}, ..., with coefficients as polynomials in delta_x = x_0 - c_0. ''' - #Note that generate_ND_derivative_relations will at worst put some power of + #Note that generate_nd_derivative_relations will at worst put some power of #$x_0^order$ in the denominator. To clear #the denominator we can probably? just multiply by x_0^order. delta_x = sp.symbols("delta_x") @@ -244,10 +244,10 @@ def compute_recurrence_relation(coeffs, n_derivs, var): #Compute symbolic derivative def hc_diff(i, n): - retMe = 1 + retme = 1 for j in range(n): - retMe *= (i-j) - return retMe + retme *= (i-j) + return retme #We are differentiating deriv_ind, which shifts down deriv_ind. #Do this for one deriv_ind From 66568194f3010e48b830aa729a8a1fc10cc6248e Mon Sep 17 00:00:00 2001 From: Hirish Chandrasekaran Date: Fri, 7 Jun 2024 11:58:31 -0700 Subject: [PATCH 05/20] Flake8 Docstring Issue --- sumpy/recurrence.py | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/sumpy/recurrence.py b/sumpy/recurrence.py index defcebf29..47b247bc4 100644 --- a/sumpy/recurrence.py +++ b/sumpy/recurrence.py @@ -39,7 +39,7 @@ def make_sympy_vec(name, n): def get_pde_in_recurrence_form(laplace): - ''' + """ get_pde_in_recurrence_form Input: - pde, a :class:`sumpy.expansion.diff_op.LinearSystemPDEOperator` pde such @@ -64,7 +64,7 @@ def get_pde_in_recurrence_form(laplace): In other words we output a linear combination of sympy variables f_{r0}, f_{r1}, ... (which represents f, f_r, f_{rr} respectively) to represent our ODE in r for the point potential. - ''' + """ dim = laplace.dim n_derivs = laplace.order assert (len(laplace.eqs) == 1) @@ -100,7 +100,7 @@ def compute_term(a, t): def generate_nd_derivative_relations(var, n_derivs): - ''' + """ generate_nd_derivative_relations Input: - var, a sympy vector of variables called [x0, x1, ...] @@ -115,7 +115,7 @@ def generate_nd_derivative_relations(var, n_derivs): Description: Using the chain rule outputs a vector that tells us how to write f, f_r, f_{rr}, ... as a linear combination of f, f_x, f_{xx}, ... - ''' + """ f_r_derivs = make_sympy_vec("f_r", n_derivs) f_x_derivs = make_sympy_vec("f_x", n_derivs) f = sp.Function("f") @@ -131,7 +131,7 @@ def generate_nd_derivative_relations(var, n_derivs): def ode_in_r_to_x(ode_in_r, var, n_derivs): - ''' + """ ode_in_r_to_x Input: - ode_in_r, a linear combination of f, f_r, f_{rr}, ... @@ -147,7 +147,7 @@ def ode_in_r_to_x(ode_in_r, var, n_derivs): Description: Translates an ode in the variable r into an ode in the variable x by substituting f, f_r, f_{rr}, ... as a linear combination of f, f_x, f_{xx}, ... using the chain rule - ''' + """ subme = generate_nd_derivative_relations(var, n_derivs) ode_in_x = ode_in_r f_r_derivs = make_sympy_vec("f_r", n_derivs) @@ -157,7 +157,7 @@ def ode_in_r_to_x(ode_in_r, var, n_derivs): def compute_poly_in_deriv(ode_in_x, n_derivs, var): - ''' + """ compute_poly_in_deriv Input: - ode_in_x, a linear combination of f, f_x, f_{xx}, ... with coefficients as @@ -172,7 +172,7 @@ def compute_poly_in_deriv(ode_in_x, n_derivs, var): Description: Converts an ode in x, to a polynomial in f, f_x, f_{xx}, ..., with coefficients as polynomials in delta_x = x_0 - c_0. - ''' + """ #Note that generate_nd_derivative_relations will at worst put some power of #$x_0^order$ in the denominator. To clear #the denominator we can probably? just multiply by x_0^order. @@ -186,7 +186,7 @@ def compute_poly_in_deriv(ode_in_x, n_derivs, var): def compute_coefficients_of_poly(poly, n_derivs): - ''' + """ compute_coefficients_of_poly Input: - poly, a polynomial in sympy variables f_{x0}, f_{x1}, ..., @@ -199,7 +199,7 @@ def compute_coefficients_of_poly(poly, n_derivs): Description: Takes in a polynomial in f_{x0}, f_{x1}, ..., w/coeffs that are polynomials in delta_x and outputs a 2d array for easy access to the coefficients based on their degree as a polynomial in delta_x. - ''' + """ delta_x = sp.symbols("delta_x") #Returns coefficients in lexographic order. So lowest order first @@ -221,7 +221,7 @@ def tup(i, n=n_derivs): def compute_recurrence_relation(coeffs, n_derivs, var): - ''' + """ compute_recurrence_relation Input: - coeffs a 2d array that gives access to the coefficients of poly, where poly @@ -237,7 +237,7 @@ def compute_recurrence_relation(coeffs, n_derivs, var): Description: Takes in coeffs which represents our ``shifted ode in x" (i.e. ode_in_x with coefficients in delta_x) and outputs a recurrence relation for the point potential. - ''' + """ i = sp.symbols("i") s = sp.Function("s") c_vec = make_sympy_vec("c", len(var)) From 633cb417193e66204dd0657167cc53eccebfe4d3 Mon Sep 17 00:00:00 2001 From: Hirish Chandrasekaran Date: Fri, 7 Jun 2024 12:59:46 -0700 Subject: [PATCH 06/20] Fix pylint issues --- sumpy/recurrence.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/sumpy/recurrence.py b/sumpy/recurrence.py index 47b247bc4..7befcb76b 100644 --- a/sumpy/recurrence.py +++ b/sumpy/recurrence.py @@ -80,6 +80,7 @@ def get_pde_in_recurrence_form(laplace): eps = sp.symbols("epsilon") rval = r + eps f = sp.Function("f") + # pylint: disable=not-callable f_derivs = [sp.diff(f(rval), eps, i) for i in range(n_derivs+1)] def compute_term(a, t): @@ -121,8 +122,10 @@ def generate_nd_derivative_relations(var, n_derivs): f = sp.Function("f") eps = sp.symbols("epsilon") rval = sp.sqrt(sum(var**2)) + eps + # pylint: disable=not-callable f_derivs_x = [sp.diff(f(rval), var[0], i) for i in range(n_derivs)] f_derivs = [sp.diff(f(rval), eps, i) for i in range(n_derivs)] + # pylint: disable=not-callable for i in range(len(f_derivs_x)): for j in range(len(f_derivs)): f_derivs_x[i] = f_derivs_x[i].subs(f_derivs[j], f_r_derivs[j]) @@ -258,6 +261,7 @@ def hc_diff(i, n): for j in range(len(coeffs[deriv_ind])-1, -1, -1): shift = pow_delta - deriv_ind + 1 pow_delta += 1 + # pylint: disable=not-callable temp = coeffs[deriv_ind][j] * s(i) * hc_diff(i, deriv_ind) part_of_r += temp.subs(i, i-shift) r += part_of_r From e39e0f58dc263af57d1bed8bb7ce9013f115f8cf Mon Sep 17 00:00:00 2001 From: Hirish Chandrasekaran Date: Mon, 17 Jun 2024 15:41:09 -0700 Subject: [PATCH 07/20] Add test for Laplace2D --- sumpy/recurrence.py | 34 ++++++++++++++++++++++++++++------ 1 file changed, 28 insertions(+), 6 deletions(-) diff --git a/sumpy/recurrence.py b/sumpy/recurrence.py index 7befcb76b..f42208e86 100644 --- a/sumpy/recurrence.py +++ b/sumpy/recurrence.py @@ -25,7 +25,8 @@ import sympy as sp from pytools.obj_array import make_obj_array - +from sumpy.expansion.diff_op import make_identity_diff_op, laplacian +import math #A similar function exists in sumpy.symbolic def make_sympy_vec(name, n): @@ -272,10 +273,31 @@ def hc_diff(i, n): return r.simplify() -def test_recurrence_finder(): - from sumpy.expansion.diff_op import make_identity_diff_op, laplacian +def test_recurrence_finder_laplace(): + """ + test_recurrence_finder_laplace + Description: Checks that the recurrence finder works correctly for the Laplace + 2D point potential. + """ w = make_identity_diff_op(2) laplace2d = laplacian(w) - print(get_pde_in_recurrence_form(laplace2d)) - - assert 1 == 1 + ode_in_r, var, n_derivs = get_pde_in_recurrence_form(laplace2d) + ode_in_x = ode_in_r_to_x(ode_in_r, var, n_derivs).simplify() + poly = compute_poly_in_deriv(ode_in_x, n_derivs, var) + coeffs = compute_coefficients_of_poly(poly, n_derivs) + i = sp.symbols("i") + s = sp.Function("s") + r = compute_recurrence_relation(coeffs, n_derivs, var) + + def coeff_laplace(i): + x, y = sp.symbols("x,y") + c_vec = make_sympy_vec("c", 2) + true_f = sp.log(sp.sqrt(x**2 + y**2)) + return sp.diff(true_f, x, i).subs(x, c_vec[0]).subs( + y, c_vec[1])/math.factorial(i) + d = 6 + # pylint: disable=not-callable + val = r.subs(i, d).subs(s(d+1),coeff_laplace(d+1)).subs( + s(d), coeff_laplace(d)).subs(s(d-1), coeff_laplace(d-1)).subs( + s(d-2), coeff_laplace(d-2)).simplify() + assert val == 0 From b987dbed4dd24c059107f407d1aaee19d67cada9 Mon Sep 17 00:00:00 2001 From: Hirish Chandrasekaran Date: Tue, 18 Jun 2024 08:54:51 -0700 Subject: [PATCH 08/20] Add test Laplace3D --- sumpy/recurrence.py | 34 ++++++++++++++++++++++++++++++++++ 1 file changed, 34 insertions(+) diff --git a/sumpy/recurrence.py b/sumpy/recurrence.py index f42208e86..347737aaf 100644 --- a/sumpy/recurrence.py +++ b/sumpy/recurrence.py @@ -301,3 +301,37 @@ def coeff_laplace(i): s(d), coeff_laplace(d)).subs(s(d-1), coeff_laplace(d-1)).subs( s(d-2), coeff_laplace(d-2)).simplify() assert val == 0 + + +def test_recurrence_finder_laplace_three_d(): + """ + test_recurrence_finder_laplace_three_d + Description: Checks that the recurrence finder works correctly for the Laplace + 3D point potential. + """ + w = make_identity_diff_op(3) + laplace3d = laplacian(w) + print(laplace3d) + ode_in_r, var, n_derivs = get_pde_in_recurrence_form(laplace3d) + ode_in_x = ode_in_r_to_x(ode_in_r, var, n_derivs).simplify() + poly = compute_poly_in_deriv(ode_in_x, n_derivs, var) + coeffs = compute_coefficients_of_poly(poly, n_derivs) + i = sp.symbols("i") + s = sp.Function("s") + r = compute_recurrence_relation(coeffs, n_derivs, var) + + def coeff_laplace_three_d(i): + x, y, z = sp.symbols("x,y,z") + c_vec = make_sympy_vec("c", 3) + true_f = 1/(sp.sqrt(x**2 + y**2 + z**2)) + return sp.diff(true_f, x, i).subs(x, c_vec[0]).subs( + y, c_vec[1]).subs(z, c_vec[2])/math.factorial(i) + + + d = 6 + # pylint: disable=not-callable + val = r.subs(i, d).subs(s(d+1),coeff_laplace_three_d(d+1)).subs( + s(d), coeff_laplace_three_d(d)).subs(s(d-1), coeff_laplace_three_d(d-1)).subs( + s(d-2), coeff_laplace_three_d(d-2)).simplify() + + assert val == 0 \ No newline at end of file From 10dbb926a4695a1fc7dc2dbc2ec159c7214975a9 Mon Sep 17 00:00:00 2001 From: Hirish Chandrasekaran Date: Tue, 18 Jun 2024 13:55:44 -0700 Subject: [PATCH 09/20] Flake8 --- sumpy/recurrence.py | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/sumpy/recurrence.py b/sumpy/recurrence.py index 347737aaf..7ae522c25 100644 --- a/sumpy/recurrence.py +++ b/sumpy/recurrence.py @@ -23,10 +23,11 @@ THE SOFTWARE. """ +import math import sympy as sp from pytools.obj_array import make_obj_array from sumpy.expansion.diff_op import make_identity_diff_op, laplacian -import math + #A similar function exists in sumpy.symbolic def make_sympy_vec(name, n): @@ -297,7 +298,7 @@ def coeff_laplace(i): y, c_vec[1])/math.factorial(i) d = 6 # pylint: disable=not-callable - val = r.subs(i, d).subs(s(d+1),coeff_laplace(d+1)).subs( + val = r.subs(i, d).subs(s(d+1), coeff_laplace(d+1)).subs( s(d), coeff_laplace(d)).subs(s(d-1), coeff_laplace(d-1)).subs( s(d-2), coeff_laplace(d-2)).simplify() assert val == 0 @@ -327,11 +328,11 @@ def coeff_laplace_three_d(i): return sp.diff(true_f, x, i).subs(x, c_vec[0]).subs( y, c_vec[1]).subs(z, c_vec[2])/math.factorial(i) - d = 6 # pylint: disable=not-callable - val = r.subs(i, d).subs(s(d+1),coeff_laplace_three_d(d+1)).subs( - s(d), coeff_laplace_three_d(d)).subs(s(d-1), coeff_laplace_three_d(d-1)).subs( + val = r.subs(i, d).subs(s(d+1), coeff_laplace_three_d(d+1)).subs( + s(d), coeff_laplace_three_d(d)).subs(s(d-1), + coeff_laplace_three_d(d-1)).subs( s(d-2), coeff_laplace_three_d(d-2)).simplify() assert val == 0 \ No newline at end of file From 558114c059fe02e082ddebdcf39367ec273b7ee8 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Mon, 24 Jun 2024 15:34:51 -0500 Subject: [PATCH 10/20] Work on improving docs --- doc/expansion.rst | 5 +++ sumpy/recurrence.py | 74 +++++++++++++++++++++++++-------------------- 2 files changed, 47 insertions(+), 32 deletions(-) diff --git a/doc/expansion.rst b/doc/expansion.rst index 5d72d735a..ea2680340 100644 --- a/doc/expansion.rst +++ b/doc/expansion.rst @@ -27,3 +27,8 @@ Estimating Expansion Orders --------------------------- .. automodule:: sumpy.expansion.level_to_order + +Recurrences +----------- + +.. automodule:: sumpy.recurrence diff --git a/sumpy/recurrence.py b/sumpy/recurrence.py index 7ae522c25..7e0b8636a 100644 --- a/sumpy/recurrence.py +++ b/sumpy/recurrence.py @@ -1,3 +1,12 @@ +""" +.. autofunction:: get_pde_in_recurrence_form +.. autofunction:: generate_nd_derivative_relations +.. autofunction:: ode_in_r_to_x +.. autofunction:: compute_poly_in_deriv +.. autofunction:: compute_coefficients_of_poly +.. autofunction:: compute_recurrence_relation +""" + __copyright__ = """ Copyright (C) 2024 Hirish Chandrasekaran Copyright (C) 2024 Andreas Kloeckner @@ -23,35 +32,32 @@ THE SOFTWARE. """ +import numpy as np import math import sympy as sp +from typing import Tuple from pytools.obj_array import make_obj_array -from sumpy.expansion.diff_op import make_identity_diff_op, laplacian +from sumpy.expansion.diff_op import ( + make_identity_diff_op, laplacian,LinearPDESystemOperator) -#A similar function exists in sumpy.symbolic -def make_sympy_vec(name, n): +# similar to make_sym_vector in sumpy.symbolic, but returns an object array +# instead of a sympy.Matrix. +def _make_sympy_vec(name, n): return make_obj_array([sp.Symbol(f"{name}{i}") for i in range(n)]) -__doc__ = """ -.. autoclass:: Recurrence -.. automodule:: sumpy.recurrence -""" - - -def get_pde_in_recurrence_form(laplace): +def get_pde_in_recurrence_form(pde: LinearPDESystemOperator) -> Tuple[ + sp.Expr, np.ndarray, int + ]: """ - get_pde_in_recurrence_form Input: - - pde, a :class:`sumpy.expansion.diff_op.LinearSystemPDEOperator` pde such - that assert(len(pde.eqs) == 1) - is true. + - *pde*, representing a scalar PDE. Output: - ode_in_r, an ode in r which the POINT-POTENTIAL (has radial symmetry) satisfies away from the origin. Note: to represent f, f_r, f_{rr}, we use the sympy variables - f_{r0}, f_{r1}, .... So ode_in_r is a linear combination of the sympy + :math:`f_{r0}`, f_{r1}, .... So ode_in_r is a linear combination of the sympy variables f_{r0}, f_{r1}, .... - var, represents the variables for the input space: [x0, x1, ...] - n_derivs, the order of the original PDE + 1, i.e. the number of @@ -67,16 +73,19 @@ def get_pde_in_recurrence_form(laplace): f_{r0}, f_{r1}, ... (which represents f, f_r, f_{rr} respectively) to represent our ODE in r for the point potential. """ - dim = laplace.dim - n_derivs = laplace.order - assert (len(laplace.eqs) == 1) - ops = len(laplace.eqs[0]) + if len(pde.eqs) != 1: + raise ValueError("PDE must be scalar") + + dim = pde.dim + n_derivs = pde.order + assert (len(pde.eqs) == 1) + ops = len(pde.eqs[0]) derivs = [] coeffs = [] - for i in laplace.eqs[0]: + for i in pde.eqs[0]: derivs.append(i.mi) - coeffs.append(laplace.eqs[0][i]) - var = make_sympy_vec("x", dim) + coeffs.append(pde.eqs[0][i]) + var = _make_sympy_vec("x", dim) r = sp.sqrt(sum(var**2)) eps = sp.symbols("epsilon") @@ -95,7 +104,7 @@ def compute_term(a, t): for i in range(ops): ode_in_r += coeffs[i] * compute_term(f(rval), derivs[i]) n_derivs = len(f_derivs) - f_r_derivs = make_sympy_vec("f_r", n_derivs) + f_r_derivs = _make_sympy_vec("f_r", n_derivs) for i in range(n_derivs): ode_in_r = ode_in_r.subs(f_derivs[i], f_r_derivs[i]) @@ -109,6 +118,7 @@ def generate_nd_derivative_relations(var, n_derivs): - var, a sympy vector of variables called [x0, x1, ...] - n_derivs, the order of the original PDE + 1, i.e. the number of derivatives of f that may be present + Output: - a vector that gives [f, f_r, f_{rr}, ...] in terms of f, f_x, f_{xx}, ... using the chain rule @@ -119,8 +129,8 @@ def generate_nd_derivative_relations(var, n_derivs): write f, f_r, f_{rr}, ... as a linear combination of f, f_x, f_{xx}, ... """ - f_r_derivs = make_sympy_vec("f_r", n_derivs) - f_x_derivs = make_sympy_vec("f_x", n_derivs) + f_r_derivs = _make_sympy_vec("f_r", n_derivs) + f_x_derivs = _make_sympy_vec("f_x", n_derivs) f = sp.Function("f") eps = sp.symbols("epsilon") rval = sp.sqrt(sum(var**2)) + eps @@ -155,7 +165,7 @@ def ode_in_r_to_x(ode_in_r, var, n_derivs): """ subme = generate_nd_derivative_relations(var, n_derivs) ode_in_x = ode_in_r - f_r_derivs = make_sympy_vec("f_r", n_derivs) + f_r_derivs = _make_sympy_vec("f_r", n_derivs) for i in range(n_derivs): ode_in_x = ode_in_x.subs(f_r_derivs[i], subme[f_r_derivs[i]]) return ode_in_x @@ -182,10 +192,10 @@ def compute_poly_in_deriv(ode_in_x, n_derivs, var): #$x_0^order$ in the denominator. To clear #the denominator we can probably? just multiply by x_0^order. delta_x = sp.symbols("delta_x") - c_vec = make_sympy_vec("c", len(var)) + c_vec = _make_sympy_vec("c", len(var)) ode_in_x_cleared = (ode_in_x * var[0]**n_derivs).simplify() ode_in_x_shifted = ode_in_x_cleared.subs(var[0], delta_x + c_vec[0]).simplify() - f_x_derivs = make_sympy_vec("f_x", n_derivs) + f_x_derivs = _make_sympy_vec("f_x", n_derivs) poly = sp.Poly(ode_in_x_shifted, *f_x_derivs) return poly @@ -245,7 +255,7 @@ def compute_recurrence_relation(coeffs, n_derivs, var): """ i = sp.symbols("i") s = sp.Function("s") - c_vec = make_sympy_vec("c", len(var)) + c_vec = _make_sympy_vec("c", len(var)) #Compute symbolic derivative def hc_diff(i, n): @@ -292,7 +302,7 @@ def test_recurrence_finder_laplace(): def coeff_laplace(i): x, y = sp.symbols("x,y") - c_vec = make_sympy_vec("c", 2) + c_vec = _make_sympy_vec("c", 2) true_f = sp.log(sp.sqrt(x**2 + y**2)) return sp.diff(true_f, x, i).subs(x, c_vec[0]).subs( y, c_vec[1])/math.factorial(i) @@ -323,7 +333,7 @@ def test_recurrence_finder_laplace_three_d(): def coeff_laplace_three_d(i): x, y, z = sp.symbols("x,y,z") - c_vec = make_sympy_vec("c", 3) + c_vec = _make_sympy_vec("c", 3) true_f = 1/(sp.sqrt(x**2 + y**2 + z**2)) return sp.diff(true_f, x, i).subs(x, c_vec[0]).subs( y, c_vec[1]).subs(z, c_vec[2])/math.factorial(i) @@ -335,4 +345,4 @@ def coeff_laplace_three_d(i): coeff_laplace_three_d(d-1)).subs( s(d-2), coeff_laplace_three_d(d-2)).simplify() - assert val == 0 \ No newline at end of file + assert val == 0 From 6afa2ea23faf52bffba492d3ed88f4bf9d9115af Mon Sep 17 00:00:00 2001 From: Hirish Chandrasekaran Date: Tue, 25 Jun 2024 13:48:13 -0700 Subject: [PATCH 11/20] Make function for producing recurrence from pde --- sumpy/recurrence.py | 30 ++++++++++++++++++++---------- 1 file changed, 20 insertions(+), 10 deletions(-) diff --git a/sumpy/recurrence.py b/sumpy/recurrence.py index 7e0b8636a..476f34434 100644 --- a/sumpy/recurrence.py +++ b/sumpy/recurrence.py @@ -284,6 +284,24 @@ def hc_diff(i, n): return r.simplify() +def get_recurrence_from_pde(pde): + """ + Input: + - *pde*, representing a scalar PDE. + + Output: + - r, a recurrence relation for a Line-Taylor expansion. + + Description: Takes in a pde, outputs a recurrence. + """ + ode_in_r, var, n_derivs = get_pde_in_recurrence_form(pde) + ode_in_x = ode_in_r_to_x(ode_in_r, var, n_derivs).simplify() + poly = compute_poly_in_deriv(ode_in_x, n_derivs, var) + coeffs = compute_coefficients_of_poly(poly, n_derivs) + r = compute_recurrence_relation(coeffs, n_derivs, var) + return r + + def test_recurrence_finder_laplace(): """ test_recurrence_finder_laplace @@ -292,13 +310,9 @@ def test_recurrence_finder_laplace(): """ w = make_identity_diff_op(2) laplace2d = laplacian(w) - ode_in_r, var, n_derivs = get_pde_in_recurrence_form(laplace2d) - ode_in_x = ode_in_r_to_x(ode_in_r, var, n_derivs).simplify() - poly = compute_poly_in_deriv(ode_in_x, n_derivs, var) - coeffs = compute_coefficients_of_poly(poly, n_derivs) + r = get_recurrence_from_pde(laplace2d) i = sp.symbols("i") s = sp.Function("s") - r = compute_recurrence_relation(coeffs, n_derivs, var) def coeff_laplace(i): x, y = sp.symbols("x,y") @@ -323,13 +337,9 @@ def test_recurrence_finder_laplace_three_d(): w = make_identity_diff_op(3) laplace3d = laplacian(w) print(laplace3d) - ode_in_r, var, n_derivs = get_pde_in_recurrence_form(laplace3d) - ode_in_x = ode_in_r_to_x(ode_in_r, var, n_derivs).simplify() - poly = compute_poly_in_deriv(ode_in_x, n_derivs, var) - coeffs = compute_coefficients_of_poly(poly, n_derivs) + r = get_recurrence_from_pde(laplace3d) i = sp.symbols("i") s = sp.Function("s") - r = compute_recurrence_relation(coeffs, n_derivs, var) def coeff_laplace_three_d(i): x, y, z = sp.symbols("x,y,z") From bd42c961fbaadbc1414b26e3e8ad472fc1e97764 Mon Sep 17 00:00:00 2001 From: Hirish Chandrasekaran Date: Wed, 26 Jun 2024 16:00:24 -0700 Subject: [PATCH 12/20] Update documentation --- sumpy/recurrence.py | 137 +++++++++++++++++++++++--------------------- 1 file changed, 72 insertions(+), 65 deletions(-) diff --git a/sumpy/recurrence.py b/sumpy/recurrence.py index 476f34434..3aefab1f0 100644 --- a/sumpy/recurrence.py +++ b/sumpy/recurrence.py @@ -31,14 +31,13 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. """ - -import numpy as np import math -import sympy as sp from typing import Tuple +import numpy as np +import sympy as sp from pytools.obj_array import make_obj_array from sumpy.expansion.diff_op import ( - make_identity_diff_op, laplacian,LinearPDESystemOperator) + make_identity_diff_op, laplacian, LinearPDESystemOperator) # similar to make_sym_vector in sumpy.symbolic, but returns an object array @@ -49,23 +48,26 @@ def _make_sympy_vec(name, n): def get_pde_in_recurrence_form(pde: LinearPDESystemOperator) -> Tuple[ sp.Expr, np.ndarray, int - ]: +]: """ Input: - *pde*, representing a scalar PDE. + Output: - - ode_in_r, an ode in r which the POINT-POTENTIAL (has radial symmetry) - satisfies away from the origin. - Note: to represent f, f_r, f_{rr}, we use the sympy variables - :math:`f_{r0}`, f_{r1}, .... So ode_in_r is a linear combination of the sympy - variables f_{r0}, f_{r1}, .... + - ode_in_r, an ode in r which the point-potential corresponding to the PDE + satisfies away from the origin. We assume that the point-potential has + radial symmetry. + Note: to represent :math:`f, f_r, f_{rr}`, we use the sympy variables + f_{r0}, f_{r1}, ... So ode_in_r is a linear combination of the + sympy variables f_{r0}, f_{r1}, ... - var, represents the variables for the input space: [x0, x1, ...] - n_derivs, the order of the original PDE + 1, i.e. the number of derivatives of f that may be present (the reason this is called n_derivs - since if we have a second order PDE for example then we might see f, f_{r}, - f_{rr} in our ODE in r, which is technically 3 terms since we count - the 0th order derivative f as a "derivative." If this doesn't make sense - just know that n_derivs is the order the of the input sumpy PDE + 1) + since if we have a second order PDE for example then we might see + :math:`f, f_{r}, f_{rr}` in our ODE in r, which is technically 3 terms + since we count the 0th order derivative f as a "derivative." If this + doesn't make sense just know that n_derivs is the order the of the input + sumpy PDE + 1) Description: We assume we are handed a system of 1 sumpy PDE (pde) and output the pde in a way that allows us to easily replace derivatives with respect to r. @@ -111,23 +113,21 @@ def compute_term(a, t): return ode_in_r, var, n_derivs -def generate_nd_derivative_relations(var, n_derivs): +def generate_nd_derivative_relations(var: np.ndarray, n_derivs: int) -> dict: """ - generate_nd_derivative_relations Input: - - var, a sympy vector of variables called [x0, x1, ...] - - n_derivs, the order of the original PDE + 1, i.e. the number of derivatives of - f that may be present + - *var*, a sympy vector of variables called [x0, x1, ...] + - *n_derivs*, the order of the original PDE + 1, i.e. the number of + derivatives of f that may be present Output: - - a vector that gives [f, f_r, f_{rr}, ...] in terms of f, f_x, f_{xx}, ... - using the chain rule - (f, f_x, f_{xx}, ... in code is represented as f_{x0}, f_{x1}, f_{x2} and - f, f_r, f_{rr}, ... in code is represented as f_{r0}, f_{r1}, f_{r2}) + - a vector that gives [f, f_r, f_{rr}, ...] in terms of f, f_x, f_{xx}, ... + using the chain rule + (f, f_x, f_{xx}, ... in code is represented as f_{x0}, f_{x1}, f_{x2} and + f, f_r, f_{rr}, ... in code is represented as f_{r0}, f_{r1}, f_{r2}) Description: Using the chain rule outputs a vector that tells us how to - write f, f_r, f_{rr}, ... as a linear - combination of f, f_x, f_{xx}, ... + write f, f_r, f_{rr}, ... as a linear combination of f, f_x, f_{xx}, ... """ f_r_derivs = _make_sympy_vec("f_r", n_derivs) f_x_derivs = _make_sympy_vec("f_x", n_derivs) @@ -145,19 +145,19 @@ def generate_nd_derivative_relations(var, n_derivs): return sp.solve(system, *f_r_derivs, dict=True)[0] -def ode_in_r_to_x(ode_in_r, var, n_derivs): +def ode_in_r_to_x(ode_in_r: sp.Expr, var: np.ndarray, n_derivs: int) -> sp.Expr: """ - ode_in_r_to_x Input: - - ode_in_r, a linear combination of f, f_r, f_{rr}, ... - (in code represented as f_{r0}, f_{r1}, f_{r2}) - with coefficients as RATIONAL functions in var[0], var[1], ... - - var, array of sympy variables [x_0, x_1, ...] - - n_derivs, the order of the original PDE + 1, i.e. the number of derivatives of - f that may be present + - *ode_in_r*, a linear combination of f, f_r, f_{rr}, ... + (in code represented as f_{r0}, f_{r1}, f_{r2}) + with coefficients as RATIONAL functions in var[0], var[1], ... + - *var*, array of sympy variables [x_0, x_1, ...] + - *n_derivs*, the order of the original PDE + 1, i.e. the number of + derivatives of f that may be present + Output: - - ode_in_x, a linear combination of f, f_x, f_{xx}, ... with coefficients as - rational functions in var[0], var[1], ... + - ode_in_x, a linear combination of f, f_x, f_{xx}, ... with coefficients as + rational functions in var[0], var[1], ... Description: Translates an ode in the variable r into an ode in the variable x by substituting f, f_r, f_{rr}, ... as a linear combination of @@ -171,19 +171,22 @@ def ode_in_r_to_x(ode_in_r, var, n_derivs): return ode_in_x -def compute_poly_in_deriv(ode_in_x, n_derivs, var): +def compute_poly_in_deriv(ode_in_x: sp.Expr, n_derivs: int, var: + np.ndarray) -> sp.polys.polytools.Poly: """ - compute_poly_in_deriv Input: - - ode_in_x, a linear combination of f, f_x, f_{xx}, ... with coefficients as - rational functions in var[0], var[1], ... - - n_derivs, the order of the original PDE + 1, i.e. the number of derivatives - of f that may be present + - *ode_in_x*, a linear combination of f, f_x, f_{xx}, ... with coefficients + as rational functions in var[0], var[1], ... + - *n_derivs*, the order of the original PDE + 1, i.e. the number of + derivatives of f that may be present + Output: - - a polynomial in f, f_x, f_{xx}, ... (in code represented as f_{x0}, f_{x1}, - f_{x2}) with coefficients as polynomials in delta_x where delta_x = x_0 - c_0 - that represents the ''shifted ODE'' - i.e. the ODE where we substitute all - occurences of delta_x with x_0 - c_0 + - a polynomial in math:`f, f_x, f_{xx}, ...` (in code represented as + math:`f_{x0}, f_{x1}, f_{x2}`) with coefficients as polynomials in + math:`\\delta_x` where + ammath:`delta_x = x_0 - c_0` that represents the ''shifted ODE'' - i.e. + the ODE where we substitute all occurences of math:`delta_x` with + math:`x_0 - c_0` Description: Converts an ode in x, to a polynomial in f, f_x, f_{xx}, ..., with coefficients as polynomials in delta_x = x_0 - c_0. @@ -200,17 +203,21 @@ def compute_poly_in_deriv(ode_in_x, n_derivs, var): return poly -def compute_coefficients_of_poly(poly, n_derivs): +def compute_coefficients_of_poly(poly: sp.polys.polytools.Poly, + n_derivs: int) -> list: """ - compute_coefficients_of_poly Input: - - poly, a polynomial in sympy variables f_{x0}, f_{x1}, ..., - (recall that this corresponds to f_0, f_x, f_{xx}, ...) with coefficients - that are polynomials in delta_x where poly represents the ''shifted ODE'' - - i.e. we substitute all occurences of delta_x with x_0 - c_0 + - *poly*, a polynomial in sympy variables math:`f_{x0}, f_{x1}, ...`, + (recall that this corresponds to math:`f_0, f_x, f_{xx}, ...`) with + coefficients that are polynomials in delta_x where poly represents the + ''shifted ODE'' i.e. we substitute all occurences of math:`\\delta_x` + with math:`x_0 - c_0` + Output: - - a 2d array, each row giving the coefficient of f_0, f_x, f_{xx}, ..., - each entry in the row giving the coefficients of the polynomial in delta_x + - coeffs, a 2d array, each row giving the coefficient of + math:`f_0, f_x, f_{xx}, ...`, each entry in the row giving the + coefficients of the polynomial in math:`\\delta_x` + Description: Takes in a polynomial in f_{x0}, f_{x1}, ..., w/coeffs that are polynomials in delta_x and outputs a 2d array for easy access to the coefficients based on their degree as a polynomial in delta_x. @@ -237,17 +244,17 @@ def tup(i, n=n_derivs): def compute_recurrence_relation(coeffs, n_derivs, var): """ - compute_recurrence_relation Input: - - coeffs a 2d array that gives access to the coefficients of poly, where poly - represents the coefficients of the ''shifted ODE'' - (''shifted ODE'' = we substitute all occurences of delta_x with x_0 - c_0) - based on their degree as a polynomial in delta_x - - n_derivs, the order of the original PDE + 1, i.e. the number of derivatives - of f that may be present + - *coeffs* a 2d array that gives access to the coefficients of poly, where + poly represents the coefficients of the ''shifted ODE'' + (''shifted ODE'' = we substitute all occurences of delta_x with x_0 - c_0) + based on their degree as a polynomial in delta_x) + - *n_derivs*, the order of the original PDE + 1, i.e. the number of + derivatives of f that may be present + Output: - - a recurrence statement that equals 0 where s(i) is the ith coefficient of - the Taylor polynomial for our point potential. + - r, a recurrence statement that equals 0 where s(i) is the ith coefficient + of the Taylor polynomial for our point potential. Description: Takes in coeffs which represents our ``shifted ode in x" (i.e. ode_in_x with coefficients in delta_x) and outputs a recurrence relation @@ -290,7 +297,8 @@ def get_recurrence_from_pde(pde): - *pde*, representing a scalar PDE. Output: - - r, a recurrence relation for a Line-Taylor expansion. + - r, a recurrence relation for a coefficients of a Line-Taylor expansion of + the point potential. Description: Takes in a pde, outputs a recurrence. """ @@ -336,7 +344,6 @@ def test_recurrence_finder_laplace_three_d(): """ w = make_identity_diff_op(3) laplace3d = laplacian(w) - print(laplace3d) r = get_recurrence_from_pde(laplace3d) i = sp.symbols("i") s = sp.Function("s") @@ -355,4 +362,4 @@ def coeff_laplace_three_d(i): coeff_laplace_three_d(d-1)).subs( s(d-2), coeff_laplace_three_d(d-2)).simplify() - assert val == 0 + assert val == 0 \ No newline at end of file From 411eb6acee5889cd6af59a852f0266e97388586b Mon Sep 17 00:00:00 2001 From: Hirish Chandrasekaran Date: Thu, 27 Jun 2024 13:17:15 -0700 Subject: [PATCH 13/20] List comprehension --- sumpy/recurrence.py | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/sumpy/recurrence.py b/sumpy/recurrence.py index 3aefab1f0..f50a106aa 100644 --- a/sumpy/recurrence.py +++ b/sumpy/recurrence.py @@ -234,10 +234,8 @@ def tup(i, n=n_derivs): a.append(1) return tuple(a) - coeffs = [] - for deriv_ind in range(n_derivs): - coeffs.append( - sp.Poly(poly.coeff_monomial(tup(deriv_ind)), delta_x).all_coeffs()) + coeffs = [sp.Poly(poly.coeff_monomial(tup(deriv_ind)), delta_x).all_coeffs() + for deriv_ind in range(n_derivs)] return coeffs From 1a5c539c273ce0f00838fc61c045cb31e39aaf37 Mon Sep 17 00:00:00 2001 From: Hirish Chandrasekaran Date: Thu, 27 Jun 2024 14:02:08 -0700 Subject: [PATCH 14/20] From hardcode to loop --- sumpy/recurrence.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/sumpy/recurrence.py b/sumpy/recurrence.py index f50a106aa..c4fe6a739 100644 --- a/sumpy/recurrence.py +++ b/sumpy/recurrence.py @@ -355,9 +355,9 @@ def coeff_laplace_three_d(i): d = 6 # pylint: disable=not-callable - val = r.subs(i, d).subs(s(d+1), coeff_laplace_three_d(d+1)).subs( - s(d), coeff_laplace_three_d(d)).subs(s(d-1), - coeff_laplace_three_d(d-1)).subs( - s(d-2), coeff_laplace_three_d(d-2)).simplify() + r_sub = r.subs(i, d) + for i in range(d-2, d+2): + r_sub = r_sub.subs(s(i), coeff_laplace_three_d(i)) + r_sub = r_sub.simplify() - assert val == 0 \ No newline at end of file + assert r_sub == 0 From 80d8eff04c26f41d152c021cc482531e11f57ba1 Mon Sep 17 00:00:00 2001 From: Hirish Chandrasekaran Date: Sun, 30 Jun 2024 17:31:00 -0700 Subject: [PATCH 15/20] Added get_recurrence_order --- sumpy/recurrence.py | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/sumpy/recurrence.py b/sumpy/recurrence.py index c4fe6a739..3b58017f3 100644 --- a/sumpy/recurrence.py +++ b/sumpy/recurrence.py @@ -289,6 +289,23 @@ def hc_diff(i, n): return r.simplify() +def get_recurrence_order(coeffs): + """ + Input: + - *coeffs*, represents coefficients of a scalar ODE. + + Output: + - true_order, the order of the recurrence relation that will be produced. + """ + orders = [] + for i in range(len(coeffs)): + for j in range(len(coeffs[i])): + if coeffs[i][j] != 0: + orders.append(i - j) + true_order = (max(orders)-min(orders)+1) + return true_order + + def get_recurrence_from_pde(pde): """ Input: From b284e475024a958964c2dffa4c8df3a72d5ce450 Mon Sep 17 00:00:00 2001 From: Hirish Chandrasekaran Date: Mon, 1 Jul 2024 10:49:35 -0700 Subject: [PATCH 16/20] Unit test for order --- sumpy/recurrence.py | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/sumpy/recurrence.py b/sumpy/recurrence.py index 3b58017f3..05a21de63 100644 --- a/sumpy/recurrence.py +++ b/sumpy/recurrence.py @@ -322,7 +322,7 @@ def get_recurrence_from_pde(pde): poly = compute_poly_in_deriv(ode_in_x, n_derivs, var) coeffs = compute_coefficients_of_poly(poly, n_derivs) r = compute_recurrence_relation(coeffs, n_derivs, var) - return r + return r, get_recurrence_order(coeffs) def test_recurrence_finder_laplace(): @@ -333,7 +333,7 @@ def test_recurrence_finder_laplace(): """ w = make_identity_diff_op(2) laplace2d = laplacian(w) - r = get_recurrence_from_pde(laplace2d) + r, order = get_recurrence_from_pde(laplace2d) i = sp.symbols("i") s = sp.Function("s") @@ -348,6 +348,7 @@ def coeff_laplace(i): val = r.subs(i, d).subs(s(d+1), coeff_laplace(d+1)).subs( s(d), coeff_laplace(d)).subs(s(d-1), coeff_laplace(d-1)).subs( s(d-2), coeff_laplace(d-2)).simplify() + assert order == 4 assert val == 0 @@ -359,7 +360,7 @@ def test_recurrence_finder_laplace_three_d(): """ w = make_identity_diff_op(3) laplace3d = laplacian(w) - r = get_recurrence_from_pde(laplace3d) + r, order = get_recurrence_from_pde(laplace3d) i = sp.symbols("i") s = sp.Function("s") @@ -376,5 +377,5 @@ def coeff_laplace_three_d(i): for i in range(d-2, d+2): r_sub = r_sub.subs(s(i), coeff_laplace_three_d(i)) r_sub = r_sub.simplify() - - assert r_sub == 0 + assert order == 4 + assert r_sub == 0 \ No newline at end of file From 5bbcfc96269cff15f6c9a58c474f1579c72f76b2 Mon Sep 17 00:00:00 2001 From: Hirish Chandrasekaran Date: Mon, 1 Jul 2024 10:54:22 -0700 Subject: [PATCH 17/20] Flake8 --- sumpy/recurrence.py | 22 ++++++++++------------ 1 file changed, 10 insertions(+), 12 deletions(-) diff --git a/sumpy/recurrence.py b/sumpy/recurrence.py index 05a21de63..ccf280953 100644 --- a/sumpy/recurrence.py +++ b/sumpy/recurrence.py @@ -117,7 +117,7 @@ def generate_nd_derivative_relations(var: np.ndarray, n_derivs: int) -> dict: """ Input: - *var*, a sympy vector of variables called [x0, x1, ...] - - *n_derivs*, the order of the original PDE + 1, i.e. the number of + - *n_derivs*, the order of the original PDE + 1, i.e. the number of derivatives of f that may be present Output: @@ -214,8 +214,8 @@ def compute_coefficients_of_poly(poly: sp.polys.polytools.Poly, with math:`x_0 - c_0` Output: - - coeffs, a 2d array, each row giving the coefficient of - math:`f_0, f_x, f_{xx}, ...`, each entry in the row giving the + - coeffs, a 2d array, each row giving the coefficient of + math:`f_0, f_x, f_{xx}, ...`, each entry in the row giving the coefficients of the polynomial in math:`\\delta_x` Description: Takes in a polynomial in f_{x0}, f_{x1}, ..., w/coeffs that are @@ -243,15 +243,15 @@ def tup(i, n=n_derivs): def compute_recurrence_relation(coeffs, n_derivs, var): """ Input: - - *coeffs* a 2d array that gives access to the coefficients of poly, where + - *coeffs* a 2d array that gives access to the coefficients of poly, where poly represents the coefficients of the ''shifted ODE'' (''shifted ODE'' = we substitute all occurences of delta_x with x_0 - c_0) based on their degree as a polynomial in delta_x) - - *n_derivs*, the order of the original PDE + 1, i.e. the number of + - *n_derivs*, the order of the original PDE + 1, i.e. the number of derivatives of f that may be present Output: - - r, a recurrence statement that equals 0 where s(i) is the ith coefficient + - r, a recurrence statement that equals 0 where s(i) is the ith coefficient of the Taylor polynomial for our point potential. Description: Takes in coeffs which represents our ``shifted ode in x" @@ -312,7 +312,7 @@ def get_recurrence_from_pde(pde): - *pde*, representing a scalar PDE. Output: - - r, a recurrence relation for a coefficients of a Line-Taylor expansion of + - r, a recurrence relation for a coefficients of a Line-Taylor expansion of the point potential. Description: Takes in a pde, outputs a recurrence. @@ -327,9 +327,8 @@ def get_recurrence_from_pde(pde): def test_recurrence_finder_laplace(): """ - test_recurrence_finder_laplace - Description: Checks that the recurrence finder works correctly for the Laplace - 2D point potential. + Description: Test the recurrence relation produced for the Laplace 2D point + potential. """ w = make_identity_diff_op(2) laplace2d = laplacian(w) @@ -354,7 +353,6 @@ def coeff_laplace(i): def test_recurrence_finder_laplace_three_d(): """ - test_recurrence_finder_laplace_three_d Description: Checks that the recurrence finder works correctly for the Laplace 3D point potential. """ @@ -378,4 +376,4 @@ def coeff_laplace_three_d(i): r_sub = r_sub.subs(s(i), coeff_laplace_three_d(i)) r_sub = r_sub.simplify() assert order == 4 - assert r_sub == 0 \ No newline at end of file + assert r_sub == 0 From c1eb7dd8e6444a1bdbb697fe0dcc3e6f0b6a3608 Mon Sep 17 00:00:00 2001 From: Hirish Chandrasekaran Date: Mon, 1 Jul 2024 14:07:09 -0700 Subject: [PATCH 18/20] Update get_recurrence_order --- sumpy/recurrence.py | 25 ++++++++++++++----------- 1 file changed, 14 insertions(+), 11 deletions(-) diff --git a/sumpy/recurrence.py b/sumpy/recurrence.py index ccf280953..c77fec574 100644 --- a/sumpy/recurrence.py +++ b/sumpy/recurrence.py @@ -32,7 +32,7 @@ THE SOFTWARE. """ import math -from typing import Tuple +from typing import Sequence, Tuple import numpy as np import sympy as sp from pytools.obj_array import make_obj_array @@ -289,21 +289,24 @@ def hc_diff(i, n): return r.simplify() -def get_recurrence_order(coeffs): +def get_recurrence_order(coeffs: Sequence[Sequence[sp.Expr]]) -> int: """ Input: - - *coeffs*, represents coefficients of a scalar ODE. - + - *coeffs*, represents coefficients of the normalized, + center-shifted ODE (se above) + with the outer sequence reflecting the order of the derivative, + and the second expansion reflecting expansion in the shift + $\delta_x$ Output: - true_order, the order of the recurrence relation that will be produced. """ - orders = [] - for i in range(len(coeffs)): - for j in range(len(coeffs[i])): - if coeffs[i][j] != 0: - orders.append(i - j) - true_order = (max(orders)-min(orders)+1) - return true_order + orders = { + i - j + for i, deriv_order_coeff in enumerate(coeffs) + for j, shift_coeff in enumerate(deriv_order_coeff) + if shift_coeff + } + return max(orders)-min(orders)+1 def get_recurrence_from_pde(pde): From c18a399b9038eb46db2eab4bfaf446a97ddfd7bc Mon Sep 17 00:00:00 2001 From: Hirish Chandrasekaran Date: Tue, 9 Jul 2024 21:10:30 -0700 Subject: [PATCH 19/20] Add parametric recurrence finder --- sumpy/recurrence.py | 107 ++++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 103 insertions(+), 4 deletions(-) diff --git a/sumpy/recurrence.py b/sumpy/recurrence.py index c77fec574..d4ce2d7ad 100644 --- a/sumpy/recurrence.py +++ b/sumpy/recurrence.py @@ -5,6 +5,10 @@ .. autofunction:: compute_poly_in_deriv .. autofunction:: compute_coefficients_of_poly .. autofunction:: compute_recurrence_relation +.. autofunction:: get_recurrence_parametric_from_pde +.. autofunction:: get_recurrence_parametric_from_coeffs +.. autofunction:: auto_product_rule_single_term +.. autofunction:: compute_coefficients_of_poly_parametric """ __copyright__ = """ @@ -293,10 +297,10 @@ def get_recurrence_order(coeffs: Sequence[Sequence[sp.Expr]]) -> int: """ Input: - *coeffs*, represents coefficients of the normalized, - center-shifted ODE (se above) + center-shifted ODE (see above) with the outer sequence reflecting the order of the derivative, and the second expansion reflecting expansion in the shift - $\delta_x$ + math:`\\delta_x` Output: - true_order, the order of the recurrence relation that will be produced. """ @@ -328,6 +332,101 @@ def get_recurrence_from_pde(pde): return r, get_recurrence_order(coeffs) +def compute_coefficients_of_poly_parametric(poly, n_derivs, var): + """ + Input: + - *poly*, a polynomial in sympy variables math:`f_{x0}, f_{x1}, ...`, + (recall that this corresponds to math:`f_0, f_x, f_{xx}, ...`) with + coefficients that are polynomials in math:`x_0` where poly represents the + TRUE ODE. + - *n_derivs*, the order of the original PDE + 1, i.e. the number of + derivatives of f that may be present + - *var*, array of sympy variables [x_0, x_1, ...] + + Output: + - coeffs, a 2d array, each row giving the coefficient of + math:`f_0, f_x, f_{xx}, ...`, each entry in the row giving the + coefficients of the polynomial in math:`x_0` + + Description: Takes in a polynomial in f_{x0}, f_{x1}, ..., w/coeffs that are + polynomials in math:`x_0` and outputs a 2d array for easy access to the + coefficients based on their degree as a polynomial in math:`x_0`. + """ + def tup(i, n=n_derivs): + a = [] + for j in range(n): + if j != i: + a.append(0) + else: + a.append(1) + return tuple(a) + + coeffs = [] + for deriv_ind in range(n_derivs): + coeffs.append(sp.Poly(poly.coeff_monomial(tup(deriv_ind)), + var[0]).all_coeffs()[::-1]) + + return coeffs + + +def auto_product_rule_single_term(p, m, var): + """ + Input: + - *p*, degree of monomial + - *m*, order of derivative + + Output: + - recurrence relation for ODE math:`x_0^p f^(m)(x_0)` + """ + n = sp.symbols("n") + s = sp.Function("s") + result = 0 + for i in range(p+1): + temp = 1 + for j in range(i): + temp *= (n - j) + # pylint: disable=not-callable + temp *= math.comb(p, i) * s(n-i+m) * var[0]**(p-i) + result += temp + return result + + +def get_recurrence_parametric_from_coeffs(coeffs, var): + """ + Input: + - *coeffs* + + Output: + - recurrence relation for full ODE + """ + final_recurrence = 0 + for m, _ in enumerate(coeffs): + for p, _ in enumerate(coeffs[m]): + final_recurrence += coeffs[m][p] * auto_product_rule_single_term(p, + m, var) + return final_recurrence + + +def get_recurrence_parametric_from_pde(pde): + """ + Input: + - *pde*, representing a scalar PDE. + + Output: + - r, a recurrence relation for a coefficients of a Line-Taylor expansion of + the point potential. + + Description: Takes in a pde, outputs a recurrence. + """ + ode_in_r, var, n_derivs = get_pde_in_recurrence_form(pde) + ode_in_x = ode_in_r_to_x(ode_in_r, var, n_derivs).simplify() + ode_in_x_cleared = (ode_in_x * var[0]**n_derivs).simplify() + f_x_derivs = _make_sympy_vec("f_x", n_derivs) + poly = sp.Poly(ode_in_x_cleared, *f_x_derivs) + coeffs = compute_coefficients_of_poly_parametric(poly, n_derivs, var) + return get_recurrence_parametric_from_coeffs(coeffs, var) + + def test_recurrence_finder_laplace(): """ Description: Test the recurrence relation produced for the Laplace 2D point @@ -361,7 +460,7 @@ def test_recurrence_finder_laplace_three_d(): """ w = make_identity_diff_op(3) laplace3d = laplacian(w) - r, order = get_recurrence_from_pde(laplace3d) + r, _ = get_recurrence_from_pde(laplace3d) i = sp.symbols("i") s = sp.Function("s") @@ -378,5 +477,5 @@ def coeff_laplace_three_d(i): for i in range(d-2, d+2): r_sub = r_sub.subs(s(i), coeff_laplace_three_d(i)) r_sub = r_sub.simplify() - assert order == 4 + #assert order == 4 assert r_sub == 0 From 671df6813da9953504125c8f1012ef31ab34ce6b Mon Sep 17 00:00:00 2001 From: Hirish Chandrasekaran Date: Wed, 10 Jul 2024 09:30:45 -0700 Subject: [PATCH 20/20] Update recurrence.py --- sumpy/recurrence.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/sumpy/recurrence.py b/sumpy/recurrence.py index d4ce2d7ad..f4173da8a 100644 --- a/sumpy/recurrence.py +++ b/sumpy/recurrence.py @@ -377,6 +377,7 @@ def auto_product_rule_single_term(p, m, var): Output: - recurrence relation for ODE math:`x_0^p f^(m)(x_0)` + s(i) """ n = sp.symbols("n") s = sp.Function("s") @@ -394,12 +395,14 @@ def auto_product_rule_single_term(p, m, var): def get_recurrence_parametric_from_coeffs(coeffs, var): """ Input: - - *coeffs* + - *coeffs*, take the ODE Output: - recurrence relation for full ODE """ final_recurrence = 0 + #Outer loop is derivative direction + #Inner is polynomial order of x_0 for m, _ in enumerate(coeffs): for p, _ in enumerate(coeffs[m]): final_recurrence += coeffs[m][p] * auto_product_rule_single_term(p,