Refactor make_interp_spline to use a single solve

examples/1d-interp.py

@@ -14,7 +14,7 @@
 import ndsplines
 
 def gaussian(x_in):
-    z = norm.ppf(.995)
+    z = norm.ppf(.9995)
     x = z*(2*x_in-1)
     return norm.pdf(x)
 
@@ -28,44 +28,72 @@ def tanh(x_in):
 
 funcs = [gaussian, sin, tanh]
 
-x = np.linspace(0, 1, 7)
+x = np.linspace(0, 1, 9)
 xx = np.linspace(-.25, 1.25, 1024)
 k = 3
 
 
-test_bcs = np.array(list(itertools.chain(
-    itertools.product(["natural", "clamped"], repeat=2),
+scipy_test_bcs = np.array(list(itertools.chain(
+    itertools.product(["natural", "clamped", ], repeat=2),
     ((None,None),),
 )))
 
-NDspline_dict = {"natural": ndsplines.pinned, "clamped": ndsplines.clamped, None: 0}
+k0_bcs = np.array(list(
+itertools.product([(0,0), (0,-1)], repeat=2)
+)[:-1])[[-1,0,1]]
 
-for func in funcs:
-    fvals = func(x)
-    truef = func(xx)
-    plt.figure()
+NDspline_dict = {"natural": ndsplines.pinned, "clamped": ndsplines.clamped, "not-a-knot": ndsplines.notaknot}
 
-    plot_sel = slice(None)
+ndsplines_test_bcs = np.array([(NDspline_dict[item[0]], NDspline_dict[item[1]],)  for item in itertools.chain(
+    itertools.product(["natural", "clamped", ], repeat=2),
+    (("not-a-knot","not-a-knot"),),
+    itertools.product(["not-a-knot"],["natural", "clamped", ]),
+    itertools.product(["natural", "clamped", ], ["not-a-knot"]),
+)])
 
-    plt.gca().set_prop_cycle(None)
 
-    for test_bc in test_bcs[plot_sel,:]:
-        test_Bspline = interpolate.make_interp_spline(x, fvals, bc_type=list(test_bc))
-
-        splinef = test_Bspline(xx.copy(), extrapolate=True)
-        plt.plot(xx, splinef,'--', lw=3.0, label=str(test_bc) + ' (BSpline)')
-
-    plt.gca().set_prop_cycle(None)
-
-    for test_bc in test_bcs[plot_sel,:]:
-        test_NDBspline = ndsplines.make_interp_spline(x, fvals, bcs=(NDspline_dict[test_bc[0]], NDspline_dict[test_bc[1]]))
-        NDsplienf = test_NDBspline(xx.copy())
-        plt.plot(xx, NDsplienf[0], label=str(test_bc) + ' (ndspline)' )
+NDspline_bc_to_string = {tuple(v):k for k,v in NDspline_dict.items()}
+NDspline_bc_to_string[(0,-1)] = 'one-sided hold'
 
+for order in range(3,4):
+    for func in funcs:
+        fvals = func(x)
+        truef = func(xx)
+        plt.figure()
+
+        plot_sel = slice(None)
 
-    plt.plot(xx, truef, 'k--', label="True " + func.__name__)
+        plt.gca().set_prop_cycle(None)
 
-    plt.legend(loc='best')
-    plt.plot(x, fvals, 'o')
-    plt.show()
+        for test_bc in scipy_test_bcs[plot_sel,:]:
+            try:
+                test_Bspline = interpolate.make_interp_spline(x, fvals, bc_type=list(test_bc), k=order)
+            except ValueError:
+                continue
+            else:
+                splinef = test_Bspline(xx.copy(), extrapolate=True)
+                plt.plot(xx, splinef, '--', lw=3.0, label=str(test_bc) + ' (BSpline)')
+
+        plt.gca().set_prop_cycle(None)
+
+        if order == 0:
+            bc_iter = k0_bcs
+        else:
+            bc_iter = ndsplines_test_bcs
+        for test_bc in bc_iter[plot_sel,:]:
+            try:
+                test_NDBspline = ndsplines.make_interp_spline(x, fvals, bcs=test_bc, orders=order)
+            except ValueError:
+                continue
+            else:
+                NDsplienf = test_NDBspline(xx.copy())
+                plt.plot(xx, NDsplienf[0], label=', '.join([NDspline_bc_to_string[tuple(bc)] for bc in test_bc]) + ' (ndspline)' )
+
+
+        plt.plot(xx, truef, 'k--', label="True " + func.__name__)
+        plt.plot(x, fvals, 'ko')
+        plt.title('k=%d'%order)
+
+        plt.legend(loc='best')
+        plt.show()
 

examples/2d-interp.py

@@ -41,14 +41,14 @@ def func2d(x_in, y_in):
 funcs = [ wrap2d(*funcs_to_wrap) for funcs_to_wrap in itertools.combinations_with_replacement(funcs, r=2)]
 funcs.append(dist)
 
-x = np.linspace(0, 1, 7)
+x = np.linspace(0, 1, 5)
 y = np.linspace(0, 1, 7)
 
 xx = np.linspace(0,1,64) 
 yy = np.linspace(0,1,64)
 
-xx = np.linspace(-.25, 1.25, 64)
-yy = np.linspace(-.25, 1.25, 64)
+# xx = np.linspace(-.25, 1.25, 64)
+# yy = np.linspace(-.25, 1.25, 64)
 k = 3
 
 
@@ -62,13 +62,15 @@ def func2d(x_in, y_in):
 for func in funcs:
     fvals = func(meshx, meshy)
     truef = func(meshxx, meshyy)
-    test_NDBspline = ndsplines.make_interp_spline(gridxy, fvals,)
-    test_RectSpline = interpolate.RectBivariateSpline(x, y, fvals)
+    test_NDBspline = ndsplines.make_interp_spline(gridxy, fvals, orders=k)
+    test_RectSpline = interpolate.RectBivariateSpline(x, y, fvals, kx=k, ky=k)
 
     fig = plt.figure()
     ax = fig.add_subplot(111, projection='3d')
 
-    ax.plot_wireframe(meshxx, meshyy, truef, alpha=0.25, color='C0')
-    ax.plot_wireframe(meshxx, meshyy, test_NDBspline(gridxxyy)[0], color='C1')
-    ax.plot_wireframe(meshxx, meshyy, test_RectSpline(meshxx, meshyy, grid=False), color='C2')
+    ax.plot_wireframe(meshxx, meshyy, truef, alpha=0.25, color='C0', label='True')
+    ax.plot_wireframe(meshxx, meshyy, test_NDBspline(gridxxyy)[0], color='C1', label='ND')
+    ax.plot_wireframe(meshxx, meshyy, test_RectSpline(meshxx, meshyy, grid=False), color='C2', label='Rect')
+
+    plt.legend(loc='best')
     plt.show()