Repository: treverhines/RBF

Issue:
Removal of lambdify
Hello, cool package you have here! The PDE weight matrix functionality in particular has been rather useful. Just wondering, was there any reason why the ```set_symbolic_to_numeric_method``` toggle to use ```lambdify``` needed to be removed entirely? Admittedly, usecases involving calculations >15 dimensions that don't cause most computers to burst into flames through their sheer memory requirements alone seem quite scarce, but at least having the option to do so was nice even if it came at a performance cost. This is with reference to commit [3cad702](https://github.com/treverhines/RBF/commit/3cad7024f8da4e19682f3cbecc4e92f2ae1f234c).

Discussion:

Patch:
diff --git a/rbf/basis.py b/rbf/basis.py
index ede9778..41a75b9 100644
--- a/rbf/basis.py
+++ b/rbf/basis.py
@@ -61,6 +61,7 @@ import numpy as np
 from scipy.sparse import csc_matrix
 from scipy.spatial import cKDTree
 from sympy.utilities.autowrap import ufuncify
+from sympy import lambdify
 
 from rbf.utils import assert_shape
 
@@ -278,7 +279,15 @@ class RBF(object):
 
         Notes
         -----
-        The number of spatial dimensions currently cannot exceed 15.
+        The default method for converting the symbolic RBF to a numeric
+        function limits the number of spatial dimensions `D` to 15. If there
+        are more than 15 dimensions, the `lambdify` method is automatically
+        used instead. However, functions produced by this method are unable
+        to perform in-place operations. Hence, if both `D` > 15 and `out` 
+        are specified, in-place operations are emulated for compatibility.
+        This preserves syntactic consistency with code in 15 dimensions or 
+        less, but incurs memory and performance overhead costs due to the
+        creation of intermediate arrays.
 
         The derivative order can be arbitrarily high, but some RBFs, such as
         Wendland and Matern, become numerically unstable when the derivative
@@ -315,7 +324,14 @@ class RBF(object):
         # reshape eps from (..., m) to (..., 1, m)
         eps = eps[..., None, :]
         # evaluate the cached function for the given `x`, `c`, and `eps`
-        out = self._cache[diff](*x, *c, eps, out=out)
+        func = self._cache[diff]
+        if out is None:
+            out = func(*x, *c, eps)
+        elif isinstance(func, np.ufunc):
+            #dim <= 15 ufunc, supports in-place ops
+            out = func(*x, *c, eps, out=out)
+        else:
+            out[...] = func(*x, *c, eps)
         return out
 
     def center_value(self, eps=1.0, diff=(0,)):
@@ -403,9 +419,13 @@ class RBF(object):
             # <= supp` and 0 otherwise.
             expr = sympy.Piecewise((expr, r_sym <= self.supp), (0, True))
 
-        func = ufuncify(
-            x_sym + c_sym + (EPS,), expr, backend='numpy', tempdir=tempdir
-            )
+        if dim <= 15:
+            func = ufuncify(
+                x_sym + c_sym + (EPS,), expr, backend='numpy', tempdir=tempdir
+                )
+
+        else:
+            func = lambdify(x_sym + c_sym + (EPS,), expr, modules=['numpy'])
 
         self._cache[diff] = func
         logger.debug('The numeric function has been created and cached.')


Test patch:
diff --git a/test/test_basis.py b/test/test_basis.py
index fcac012..75e781f 100644
--- a/test/test_basis.py
+++ b/test/test_basis.py
@@ -212,4 +212,17 @@ class Test(unittest.TestCase):
     rbf.basis.phs2(x, y, out=out5)
     self.assertTrue(np.allclose(out1, out5))
 
+    # test dimensions > 15
+    x = np.random.random((5, 16))
+    y = np.random.random((4, 16))
+    out6 = rbf.basis.phs2(x, y)
+
+    out7 = np.empty((5, 4))
+    rbf.basis.phs2(x, y, out=out7)
+    self.assertTrue(np.allclose(out6, out7))
+
+    out8 = np.empty((4, 5)).T
+    rbf.basis.phs2(x, y, out=out8)
+    self.assertTrue(np.allclose(out6, out8))
+
 #unittest.main()