diff --git a/sumpy/expansion/__init__.py b/sumpy/expansion/__init__.py
index 108bc67db749a50b17ac36a302675b4486df2787..9bd5b2d200bfed58397f3fffa9db51e5797eb486 100644
--- a/sumpy/expansion/__init__.py
+++ b/sumpy/expansion/__init__.py
@@ -372,32 +372,81 @@ class LinearRecurrenceBasedDerivativeWrangler(DerivativeWrangler):
return LinearRecurrenceBasedMiDerivativeTaker(expr, var_list, self)
-class LaplaceDerivativeWrangler(LinearRecurrenceBasedDerivativeWrangler):
+class NewLinearRecurrenceBasedDerivativeWrangler \
+ (LinearRecurrenceBasedDerivativeWrangler):
- def try_get_recurrence_for_derivative(self, deriv, in_terms_of):
- deriv = np.array(deriv, dtype=int)
+ @memoize_method
+ def _get_stored_ids_and_coeff_mat(self):
+ from scipy.linalg.interpolative import interp_decomp
+ from six import iteritems
+ from sumpy.tools import nullspace
- for dim in np.where(2 <= deriv)[0]:
- # Check if we can reduce this dimension in terms of the other
- # dimensions.
+ tol = 1e-13
+ stored_identifiers = []
- reduced_deriv = deriv.copy()
- reduced_deriv[dim] -= 2
- coeffs = {}
+ import time
+ start_time = time.time()
+ logger.debug("computing recurrence for Taylor coefficients: start")
- for other_dim in range(self.dim):
- if other_dim == dim:
- continue
- needed_deriv = reduced_deriv.copy()
- needed_deriv[other_dim] += 2
- needed_deriv = tuple(needed_deriv)
+ pde_dict = self.get_pde_dict()
+ P = []
+
+ mis = self.get_full_coefficient_identifiers()
+ coeff_ident_enumerate_dict = dict((tuple(mi), i) for
+ (i, mi) in enumerate(mis))
+
+ for mi in mis:
+ for pde_mi, coeff in iteritems(pde_dict):
+ diff = np.array(mi, dtype=int) - pde_mi
+ if (diff >= 0).all():
+ eq = [0]*len(mis)
+ for pde_mi2, coeff2 in iteritems(pde_dict):
+ c = tuple(pde_mi2 + diff)
+ eq[coeff_ident_enumerate_dict[c]] = 1
+ P.append(eq)
+
+ P = np.array(P, dtype=np.float64)
+ N = nullspace(P, atol=tol)
+ k, idx, _ = interp_decomp(N.T, tol)
+ idx = idx[:k]
+ S = np.linalg.solve(N[idx, :].T, N.T).T
+ stored_identifiers = [mis[i] for i in idx]
+
+ coeff_matrix = defaultdict(lambda: [])
+ for i in range(S.shape[0]):
+ for j in range(S.shape[1]):
+ if np.abs(S[i][j]) > tol:
+ coeff_matrix[i].append((j, S[i][j]))
- if needed_deriv not in in_terms_of:
- break
+ logger.debug("computing recurrence for Taylor coefficients: "
+ "done after {dur:.2f} seconds"
+ .format(dur=time.time() - start_time))
- coeffs[needed_deriv] = -1
- else:
- return coeffs
+ logger.debug("number of Taylor coefficients was reduced from {orig} to {red}"
+ .format(orig=len(self.get_full_coefficient_identifiers()),
+ red=len(stored_identifiers)))
+
+ return stored_identifiers, coeff_matrix
+
+ @memoize_method
+ def get_coefficient_matrix(self, rscale):
+ _, coeff_matrix = self._get_stored_ids_and_coeff_mat()
+ return coeff_matrix
+
+ def get_derivative_taker(self, expr, var_list):
+ from sumpy.tools import NewLinearRecurrenceBasedMiDerivativeTaker
+ return NewLinearRecurrenceBasedMiDerivativeTaker(expr, var_list, self)
+
+
+class LaplaceDerivativeWrangler(NewLinearRecurrenceBasedDerivativeWrangler):
+
+ def get_pde_dict(self):
+ pde_dict = {}
+ for d in range(self.dim):
+ mi = [0]*self.dim
+ mi[d] = 2
+ pde_dict[tuple(mi)] = 1
+ return pde_dict
class HelmholtzDerivativeWrangler(LinearRecurrenceBasedDerivativeWrangler):
diff --git a/sumpy/tools.py b/sumpy/tools.py
index 964c64cc440c807c02b3f8811c6a3f0de6c84b04..e2f6f86b36ad632809c72e0b51a8bdca926349f1 100644
--- a/sumpy/tools.py
+++ b/sumpy/tools.py
@@ -139,6 +139,37 @@ class LinearRecurrenceBasedMiDerivativeTaker(MiDerivativeTaker):
return expr
+
+class NewLinearRecurrenceBasedMiDerivativeTaker(MiDerivativeTaker):
+ """
+ The derivative taker for expansions that use
+ :class:`sumpy.expansion.LinearRecurrenceBasedDerivativeWrangler`
+ """
+
+ def __init__(self, expr, var_list, wrangler):
+ super(NewLinearRecurrenceBasedMiDerivativeTaker, self).__init__(
+ expr, var_list)
+ self.wrangler = wrangler
+
+ @memoize_method
+ def diff(self, mi):
+ """
+ :arg mi: a multi-index (tuple) indicating how many x/y derivatives are
+ to be taken.
+ """
+ try:
+ expr = self.cache_by_mi[mi]
+ except KeyError:
+ closest_mi = self.get_closest_cached_mi(mi)
+ expr = self.cache_by_mi[closest_mi]
+
+ for next_deriv, next_mi in (
+ self.get_derivative_taking_sequence(closest_mi, mi)):
+ expr = expr.diff(next_deriv)
+ self.cache_by_mi[next_mi] = expr
+
+ return expr
+
# }}}
@@ -448,4 +479,48 @@ def my_syntactic_subs(expr, subst_dict):
return expr
+# Source: http://scipy-cookbook.readthedocs.io/items/RankNullspace.html
+# Author: SciPy Developers
+# License: BSD-3-Clause
+
+def nullspace(A, atol=1e-13, rtol=0):
+ """Compute an approximate basis for the nullspace of A.
+
+ The algorithm used by this function is based on the singular value
+ decomposition of `A`.
+
+ Parameters
+ ----------
+ A : ndarray
+ A should be at most 2-D. A 1-D array with length k will be treated
+ as a 2-D with shape (1, k)
+ atol : float
+ The absolute tolerance for a zero singular value. Singular values
+ smaller than `atol` are considered to be zero.
+ rtol : float
+ The relative tolerance. Singular values less than rtol*smax are
+ considered to be zero, where smax is the largest singular value.
+
+ If both `atol` and `rtol` are positive, the combined tolerance is the
+ maximum of the two; that is::
+ tol = max(atol, rtol * smax)
+ Singular values smaller than `tol` are considered to be zero.
+
+ Return value
+ ------------
+ ns : ndarray
+ If `A` is an array with shape (m, k), then `ns` will be an array
+ with shape (k, n), where n is the estimated dimension of the
+ nullspace of `A`. The columns of `ns` are a basis for the
+ nullspace; each element in numpy.dot(A, ns) will be approximately
+ zero.
+ """
+ from numpy.linalg import svd
+ A = np.atleast_2d(A)
+ u, s, vh = svd(A)
+ tol = max(atol, rtol * s[0])
+ nnz = (s >= tol).sum()
+ ns = vh[nnz:].conj().T
+ return ns
+
# vim: fdm=marker