From 28749dd6a7761c12c79f125a01efa202bc07a748 Mon Sep 17 00:00:00 2001
From: ben <ben@sixpearls.com>
Date: Tue, 11 Jun 2019 16:21:53 -0700
Subject: [PATCH 01/10] minimal port of scipy make_interp_spline

---
 ndsplines/ndsplines.py | 174 ++++++++++++++++++++++++++++++++++++++---
 1 file changed, 163 insertions(+), 11 deletions(-)

diff --git a/ndsplines/ndsplines.py b/ndsplines/ndsplines.py
index 13f6f69..66b754c 100644
--- a/ndsplines/ndsplines.py
+++ b/ndsplines/ndsplines.py
@@ -13,6 +13,11 @@
 __all__ = ['pinned', 'clamped', 'extrap', 'periodic', 'BSplineNDInterpolator',
            'make_interp_spline', 'make_lsq_spline', 'default_implementation']
 
+import operator
+from scipy._lib.six import string_types
+from scipy.linalg import (get_lapack_funcs, LinAlgError,
+                          cholesky_banded, cho_solve_banded)
+
 """
 TODOs:
 
@@ -211,6 +216,9 @@ def make_lsq_spline(x, y, knots, orders, w=None, check_finite=True):
     return temp_spline
 
 
+from scipy.interpolate._bsplines import _as_float_array, _get_dtype, _augknt, _convert_string_aliases, _process_deriv_spec, _bspl, prod, _not_a_knot
+
+
 def make_interp_spline(x, y, bcs=0, orders=3):
     """
     Construct an interpolating B-spline.
@@ -266,13 +274,36 @@ def make_interp_spline(x, y, bcs=0, orders=3):
     knots = []
     coefficients = np.pad(y, np.r_[np.c_[0,0], deriv_specs], 'constant')
 
+    axis=0
+    check_finite=True
+
     for i in np.arange(ndim)+1:
         all_other_ax_shape = np.asarray(np.r_[coefficients.shape[1:i],
             y.shape[i+1:]], dtype=np.int)
         x_line_sel = ((i-1,) + (0,)*(i-1) + (slice(None,None),) +
             (0,)*(ndim-i))
-        xp = x[x_line_sel]
-        order = orders[i-1]
+        x_slice = x[x_line_sel]
+        k = orders[i-1]
+
+        bc_type=(bc_map[(bcs[i-1,0])],
+                 bc_map[(bcs[i-1,1])])
+
+        #=======================================
+        # START FROM SCIPY
+        #=======================================
+        if bc_type is None or bc_type == 'not-a-knot':
+            deriv_l, deriv_r = None, None
+        elif isinstance(bc_type, string_types):
+            deriv_l, deriv_r = bc_type, bc_type
+        else:
+            try:
+                deriv_l, deriv_r = bc_type
+            except TypeError:
+                raise ValueError("Unknown boundary condition: %s" % bc_type)
+        #=======================================
+        # END FROM SCIPY
+        #=======================================
+
         for idx in np.ndindex(*all_other_ax_shape):
             offset_axes_remaining_sel = (tuple(idx[i-1:] + 
                 deriv_specs[i:,0]))
@@ -281,15 +312,136 @@ def make_interp_spline(x, y, bcs=0, orders=3):
                 offset_axes_remaining_sel)
             coeff_line_sel = ((Ellipsis,) + idx[:i-1] + (slice(None,None),)
                 + offset_axes_remaining_sel)
-            line_spline = interpolate.make_interp_spline(xp,
-                coefficients[y_line_sel].T,
-                k = order,
-                bc_type=(bc_map[(bcs[i-1,0])],
-                         bc_map[(bcs[i-1,1])]),
-
-            )
-            coefficients[coeff_line_sel] = line_spline.c.T
-        knots.append(line_spline.t)
+
+            y_slice = coefficients[y_line_sel].T
+
+            t = None
+
+            #=======================================
+            # START FROM SCIPY
+            #=======================================
+
+            if not -y_slice.ndim <= axis < y_slice.ndim:
+                raise ValueError("axis {} is out of bounds".format(axis))
+            if axis < 0:
+                axis += y_slice.ndim
+
+            # special-case k=0 right away
+            if k == 0:
+                if any(_ is not None for _ in (t, deriv_l, deriv_r)):
+                    raise ValueError("Too much info for k=0: t and bc_type can only "
+                                     "be None.")
+                x_slice = _as_float_array(x_slice, check_finite)
+                t = np.r_[x_slice, x_slice[-1]]
+                c = np.asarray(y_slice)
+                c = np.rollaxis(c, axis)
+                c = np.ascontiguousarray(c, dtype=_get_dtype(c.dtype))
+                return BSpline.construct_fast(t, c, k, axis=axis)
+
+            # special-case k=1 (e.g., Lyche and Morken, Eq.(2.16))
+            if k == 1 and t is None:
+                if not (deriv_l is None and deriv_r is None):
+                    raise ValueError("Too much info for k=1: bc_type can only be None.")
+                x_slice = _as_float_array(x_slice, check_finite)
+                t = np.r_[x_slice[0], x_slice, x_slice[-1]]
+                c = np.asarray(y_slice)
+                c = np.rollaxis(c, axis)
+                c = np.ascontiguousarray(c, dtype=_get_dtype(c.dtype))
+                return BSpline.construct_fast(t, c, k, axis=axis)
+
+            x_slice = _as_float_array(x_slice, check_finite)
+            y_slice = _as_float_array(y_slice, check_finite)
+            k = operator.index(k)
+
+            # come up with a sensible knot vector, if needed
+            if t is None:
+                if deriv_l is None and deriv_r is None:
+                    if k == 2:
+                        # OK, it's a bit ad hoc: Greville sites + omit
+                        # 2nd and 2nd-to-last points, a la not-a-knot
+                        t = (x_slice[1:] + x_slice[:-1]) / 2.
+                        t = np.r_[(x_slice[0],)*(k+1),
+                                   t[1:-1],
+                                   (x_slice[-1],)*(k+1)]
+                    else:
+                        t = _not_a_knot(x_slice, k)
+                else:
+                    t = _augknt(x_slice, k)
+
+            t = _as_float_array(t, check_finite)
+
+            y_slice = np.rollaxis(y_slice, axis)    # now internally interp axis is zero
+
+            if x_slice.ndim != 1 or np.any(x_slice[1:] <= x_slice[:-1]):
+                raise ValueError("Expect x_slice to be a 1-D sorted array_like.")
+            if k < 0:
+                raise ValueError("Expect non-negative k.")
+            if t.ndim != 1 or np.any(t[1:] < t[:-1]):
+                raise ValueError("Expect t to be a 1-D sorted array_like.")
+            if x_slice.size != y_slice.shape[0]:
+                raise ValueError('x_slice and y_slice are incompatible.')
+            if t.size < x_slice.size + k + 1:
+                raise ValueError('Got %d knots, need at least %d.' %
+                                 (t.size, x_slice.size + k + 1))
+            if (x_slice[0] < t[k]) or (x_slice[-1] > t[-k]):
+                raise ValueError('Out of bounds w/ x_slice = %s.' % x_slice)
+
+            # Here : deriv_l, r = [(nu, value), ...]
+            deriv_l = _convert_string_aliases(deriv_l, y_slice.shape[1:])
+            deriv_l_ords, deriv_l_vals = _process_deriv_spec(deriv_l)
+            nleft = deriv_l_ords.shape[0]
+
+            deriv_r = _convert_string_aliases(deriv_r, y_slice.shape[1:])
+            deriv_r_ords, deriv_r_vals = _process_deriv_spec(deriv_r)
+            nright = deriv_r_ords.shape[0]
+
+            # have `n` conditions for `nt` coefficients; need nt-n derivatives
+            n = x_slice.size
+            nt = t.size - k - 1
+
+            if nt - n != nleft + nright:
+                raise ValueError("The number of derivatives at boundaries does not "
+                                 "match: expected %s, got %s+%s" % (nt-n, nleft, nright))
+
+            # set up the LHS: the collocation matrix + derivatives at boundaries
+            kl = ku = k
+            ab = np.zeros((2*kl + ku + 1, nt), dtype=np.float_, order='F')
+            _bspl._colloc(x_slice, t, k, ab, offset=nleft)
+            if nleft > 0:
+                _bspl._handle_lhs_derivatives(t, k, x_slice[0], ab, kl, ku, deriv_l_ords)
+            if nright > 0:
+                _bspl._handle_lhs_derivatives(t, k, x_slice[-1], ab, kl, ku, deriv_r_ords,
+                                        offset=nt-nright)
+
+            # set up the RHS: values to interpolate (+ derivative values, if any)
+            extradim = prod(y_slice.shape[1:])
+            rhs = np.empty((nt, extradim), dtype=y_slice.dtype)
+            if nleft > 0:
+                rhs[:nleft] = deriv_l_vals.reshape(-1, extradim)
+            rhs[nleft:nt - nright] = y_slice.reshape(-1, extradim)
+            if nright > 0:
+                rhs[nt - nright:] = deriv_r_vals.reshape(-1, extradim)
+
+            # solve Ab @ x_slice = rhs; this is the relevant part of linalg.solve_banded
+            if check_finite:
+                ab, rhs = map(np.asarray_chkfinite, (ab, rhs))
+            gbsv, = get_lapack_funcs(('gbsv',), (ab, rhs))
+            lu, piv, c, info = gbsv(kl, ku, ab, rhs,
+                    overwrite_ab=True, overwrite_b=True)
+
+            if info > 0:
+                raise LinAlgError("Collocation matix is singular.")
+            elif info < 0:
+                raise ValueError('illegal value in %d-th argument of internal gbsv' % -info)
+
+            c = np.ascontiguousarray(c.reshape((nt,) + y_slice.shape[1:]))
+
+            #=======================================
+            # START FROM SCIPY
+            #=======================================
+
+            coefficients[coeff_line_sel] = c.T
+        knots.append(t)
     return BSplineNDInterpolator(knots, coefficients, orders, 
         np.all(bcs==periodic, axis=1),
         (bcs >= 0))

From 1c3609699a566de90471e60ed827898b09780818 Mon Sep 17 00:00:00 2001
From: ben <ben@sixpearls.com>
Date: Wed, 12 Jun 2019 08:31:09 -0700
Subject: [PATCH 02/10] re-arranged so knots are only calculated once per
 dimension

---
 ndsplines/ndsplines.py | 95 ++++++++++++++++++++++--------------------
 1 file changed, 50 insertions(+), 45 deletions(-)

diff --git a/ndsplines/ndsplines.py b/ndsplines/ndsplines.py
index 66b754c..e8f7647 100644
--- a/ndsplines/ndsplines.py
+++ b/ndsplines/ndsplines.py
@@ -1,5 +1,12 @@
 import numpy as np
-from scipy import interpolate
+
+import operator
+from scipy._lib.six import string_types
+from scipy.linalg import (get_lapack_funcs, LinAlgError,
+                          cholesky_banded, cho_solve_banded)
+from scipy.interpolate._bsplines import (_as_float_array, _get_dtype, _augknt, 
+    _convert_string_aliases, _process_deriv_spec, _bspl as _sci_bspl, prod, 
+    _not_a_knot)
 
 from ndsplines import _npy_bspl
 
@@ -13,10 +20,7 @@
 __all__ = ['pinned', 'clamped', 'extrap', 'periodic', 'BSplineNDInterpolator',
            'make_interp_spline', 'make_lsq_spline', 'default_implementation']
 
-import operator
-from scipy._lib.six import string_types
-from scipy.linalg import (get_lapack_funcs, LinAlgError,
-                          cholesky_banded, cho_solve_banded)
+
 
 """
 TODOs:
@@ -216,9 +220,6 @@ def make_lsq_spline(x, y, knots, orders, w=None, check_finite=True):
     return temp_spline
 
 
-from scipy.interpolate._bsplines import _as_float_array, _get_dtype, _augknt, _convert_string_aliases, _process_deriv_spec, _bspl, prod, _not_a_knot
-
-
 def make_interp_spline(x, y, bcs=0, orders=3):
     """
     Construct an interpolating B-spline.
@@ -288,6 +289,8 @@ def make_interp_spline(x, y, bcs=0, orders=3):
         bc_type=(bc_map[(bcs[i-1,0])],
                  bc_map[(bcs[i-1,1])])
 
+        t = None
+
         #=======================================
         # START FROM SCIPY
         #=======================================
@@ -300,10 +303,44 @@ def make_interp_spline(x, y, bcs=0, orders=3):
                 deriv_l, deriv_r = bc_type
             except TypeError:
                 raise ValueError("Unknown boundary condition: %s" % bc_type)
+
+        if k == 0:
+            if any(_ is not None for _ in (t, deriv_l, deriv_r)):
+                raise ValueError("Too much info for k=0: t and bc_type can only "
+                                 "be None.")
+            x_slice = _as_float_array(x_slice, check_finite)
+            t = np.r_[x_slice, x_slice[-1]]
+
+        if k == 1 and t is None:
+            if not (deriv_l is None and deriv_r is None):
+                raise ValueError("Too much info for k=1: bc_type can only be None.")
+            x_slice = _as_float_array(x_slice, check_finite)
+            t = np.r_[x_slice[0], x_slice, x_slice[-1]]
+
+        # come up with a sensible knot vector, if needed
+        if t is None:
+            if deriv_l is None and deriv_r is None:
+                if k == 2:
+                    # OK, it's a bit ad hoc: Greville sites + omit
+                    # 2nd and 2nd-to-last points, a la not-a-knot
+                    t = (x_slice[1:] + x_slice[:-1]) / 2.
+                    t = np.r_[(x_slice[0],)*(k+1),
+                               t[1:-1],
+                               (x_slice[-1],)*(k+1)]
+                else:
+                    t = _not_a_knot(x_slice, k)
+            else:
+                t = _augknt(x_slice, k)
+
+        t = _as_float_array(t, check_finite)
+
         #=======================================
         # END FROM SCIPY
         #=======================================
 
+
+
+
         for idx in np.ndindex(*all_other_ax_shape):
             offset_axes_remaining_sel = (tuple(idx[i-1:] + 
                 deriv_specs[i:,0]))
@@ -315,62 +352,29 @@ def make_interp_spline(x, y, bcs=0, orders=3):
 
             y_slice = coefficients[y_line_sel].T
 
-            t = None
+
 
             #=======================================
             # START FROM SCIPY
             #=======================================
 
-            if not -y_slice.ndim <= axis < y_slice.ndim:
-                raise ValueError("axis {} is out of bounds".format(axis))
-            if axis < 0:
-                axis += y_slice.ndim
-
             # special-case k=0 right away
             if k == 0:
                 if any(_ is not None for _ in (t, deriv_l, deriv_r)):
                     raise ValueError("Too much info for k=0: t and bc_type can only "
                                      "be None.")
-                x_slice = _as_float_array(x_slice, check_finite)
-                t = np.r_[x_slice, x_slice[-1]]
                 c = np.asarray(y_slice)
-                c = np.rollaxis(c, axis)
-                c = np.ascontiguousarray(c, dtype=_get_dtype(c.dtype))
-                return BSpline.construct_fast(t, c, k, axis=axis)
 
             # special-case k=1 (e.g., Lyche and Morken, Eq.(2.16))
             if k == 1 and t is None:
                 if not (deriv_l is None and deriv_r is None):
                     raise ValueError("Too much info for k=1: bc_type can only be None.")
-                x_slice = _as_float_array(x_slice, check_finite)
-                t = np.r_[x_slice[0], x_slice, x_slice[-1]]
                 c = np.asarray(y_slice)
-                c = np.rollaxis(c, axis)
-                c = np.ascontiguousarray(c, dtype=_get_dtype(c.dtype))
-                return BSpline.construct_fast(t, c, k, axis=axis)
 
             x_slice = _as_float_array(x_slice, check_finite)
             y_slice = _as_float_array(y_slice, check_finite)
             k = operator.index(k)
 
-            # come up with a sensible knot vector, if needed
-            if t is None:
-                if deriv_l is None and deriv_r is None:
-                    if k == 2:
-                        # OK, it's a bit ad hoc: Greville sites + omit
-                        # 2nd and 2nd-to-last points, a la not-a-knot
-                        t = (x_slice[1:] + x_slice[:-1]) / 2.
-                        t = np.r_[(x_slice[0],)*(k+1),
-                                   t[1:-1],
-                                   (x_slice[-1],)*(k+1)]
-                    else:
-                        t = _not_a_knot(x_slice, k)
-                else:
-                    t = _augknt(x_slice, k)
-
-            t = _as_float_array(t, check_finite)
-
-            y_slice = np.rollaxis(y_slice, axis)    # now internally interp axis is zero
 
             if x_slice.ndim != 1 or np.any(x_slice[1:] <= x_slice[:-1]):
                 raise ValueError("Expect x_slice to be a 1-D sorted array_like.")
@@ -403,14 +407,15 @@ def make_interp_spline(x, y, bcs=0, orders=3):
                 raise ValueError("The number of derivatives at boundaries does not "
                                  "match: expected %s, got %s+%s" % (nt-n, nleft, nright))
 
+
             # set up the LHS: the collocation matrix + derivatives at boundaries
             kl = ku = k
             ab = np.zeros((2*kl + ku + 1, nt), dtype=np.float_, order='F')
-            _bspl._colloc(x_slice, t, k, ab, offset=nleft)
+            _sci_bspl._colloc(x_slice, t, k, ab, offset=nleft)
             if nleft > 0:
-                _bspl._handle_lhs_derivatives(t, k, x_slice[0], ab, kl, ku, deriv_l_ords)
+                _sci_bspl._handle_lhs_derivatives(t, k, x_slice[0], ab, kl, ku, deriv_l_ords)
             if nright > 0:
-                _bspl._handle_lhs_derivatives(t, k, x_slice[-1], ab, kl, ku, deriv_r_ords,
+                _sci_bspl._handle_lhs_derivatives(t, k, x_slice[-1], ab, kl, ku, deriv_r_ords,
                                         offset=nt-nright)
 
             # set up the RHS: values to interpolate (+ derivative values, if any)

From 5f0f134818f1bae5bc52bbb6f612a7e5b7dc48ca Mon Sep 17 00:00:00 2001
From: ben <ben@sixpearls.com>
Date: Wed, 12 Jun 2019 09:31:55 -0700
Subject: [PATCH 03/10] append knots after processing knots

---
 ndsplines/ndsplines.py | 6 +-----
 1 file changed, 1 insertion(+), 5 deletions(-)

diff --git a/ndsplines/ndsplines.py b/ndsplines/ndsplines.py
index e8f7647..09d3e1e 100644
--- a/ndsplines/ndsplines.py
+++ b/ndsplines/ndsplines.py
@@ -29,9 +29,6 @@
 create wrapper with callback to allow for creating anti-derivative splines, etc
 (use 1D operations that can be iterated over )
 
-break out some of the knot normalization stuff (in interpolate.make_interp_spline)
-to make it easier to make lsq splines? would also be useful for making the
-make_interp_spline more efficient
 
 make sure these can be pickled (maybe store knots, coeffs, orders, etc to matfile?
 
@@ -338,7 +335,7 @@ def make_interp_spline(x, y, bcs=0, orders=3):
         # END FROM SCIPY
         #=======================================
 
-
+        knots.append(t)
 
 
         for idx in np.ndindex(*all_other_ax_shape):
@@ -446,7 +443,6 @@ def make_interp_spline(x, y, bcs=0, orders=3):
             #=======================================
 
             coefficients[coeff_line_sel] = c.T
-        knots.append(t)
     return BSplineNDInterpolator(knots, coefficients, orders, 
         np.all(bcs==periodic, axis=1),
         (bcs >= 0))

From 94de3151fca40e3b9f52414b92b6851db522c3b8 Mon Sep 17 00:00:00 2001
From: ben <ben@sixpearls.com>
Date: Sat, 15 Jun 2019 09:42:49 -0700
Subject: [PATCH 04/10] change ndim --> xdim, mdim --> ydim. will probably
 change that, but now it will be easier since no conflict with ndarray.ndim

---
 ndsplines/ndsplines.py | 106 ++++++++++++++++++++---------------------
 1 file changed, 53 insertions(+), 53 deletions(-)

diff --git a/ndsplines/ndsplines.py b/ndsplines/ndsplines.py
index 09d3e1e..a66098f 100644
--- a/ndsplines/ndsplines.py
+++ b/ndsplines/ndsplines.py
@@ -47,29 +47,29 @@ class BSplineNDInterpolator(object):
     Parameters
     ----------
     knots : list of ndarrays,
-        shapes=[n_1+orders[i-1]+1, ..., n_ndim+orders[-1]+1], dtype=np.float_
-    coefficients : ndarray, shape=(mdim, n_1, n_2, ..., n_ndim), dtype=np.float_
-    orders : ndarray, shape=(ndim,), dtype=np.int_
-    periodic : ndarray, shape=(ndim,), dtype=np.bool_
-    extrapolate : ndarray, shape=(ndim,2), dtype=np.bool_
+        shapes=[n_1+orders[i-1]+1, ..., n_xdim+orders[-1]+1], dtype=np.float_
+    coefficients : ndarray, shape=(ydim, n_1, n_2, ..., n_xdim), dtype=np.float_
+    orders : ndarray, shape=(xdim,), dtype=np.int_
+    periodic : ndarray, shape=(xdim,), dtype=np.bool_
+    extrapolate : ndarray, shape=(xdim,2), dtype=np.bool_
     """
     def __init__(self, knots, coefficients, orders, periodic=False, extrapolate=True, implementation=default_implementation):
         self.implementation = implementation
         self.knots = knots
         self.coefficients = coefficients
-        self.ndim = len(knots) # dimension of knots
+        self.xdim = len(knots) # dimension of knots
         self.nis = np.array([len(knot) for knot in knots])
-        self.mdim = coefficients.shape[0] # dimension of coefficeints
-        self.orders = np.broadcast_to(orders, (self.ndim,))
-        self.periodic = np.broadcast_to(periodic, (self.ndim,))
-        self.extrapolate = np.broadcast_to(extrapolate, (self.ndim,2))
+        self.ydim = coefficients.shape[0] # dimension of coefficeints
+        self.orders = np.broadcast_to(orders, (self.xdim,))
+        self.periodic = np.broadcast_to(periodic, (self.xdim,))
+        self.extrapolate = np.broadcast_to(extrapolate, (self.xdim,2))
 
-        self.coefficient_op = [0,] + list(i for i in range(2,self.ndim+2)) + [1,]
-        self.u_ops = [[1, i+2] for i in range(self.ndim)]
+        self.coefficient_op = [0,] + list(i for i in range(2,self.xdim+2)) + [1,]
+        self.u_ops = [[1, i+2] for i in range(self.xdim)]
         self.output_op = [0,1]
 
         self.coefficient_selector_base = np.meshgrid(*[np.arange(order+1) for order in self.orders], indexing='ij')
-        self.coefficient_shape_base = (self.ndim,)+tuple(self.orders+1)
+        self.coefficient_shape_base = (self.xdim,)+tuple(self.orders+1)
 
         self.current_max_num_points = 0
         self.allocate_workspace_arrays(1)
@@ -80,15 +80,15 @@ def compute_basis_coefficient_selector(self, x, nus=0):
         """
         Parameters
         ----------
-        x : ndarray, shape=(self.ndim, s) dtype=np.float_
-        nus : int or ndarray, shape=(self.ndim,) dtype=np.int_
+        x : ndarray, shape=(self.xdim, s) dtype=np.float_
+        nus : int or ndarray, shape=(self.xdim,) dtype=np.int_
         """
         num_points = x.shape[-1]
 
         if not isinstance(nus, np.ndarray):
             nu = nus
 
-        for i in range(self.ndim):
+        for i in range(self.xdim):
             t = self.knots[i]
             k = self.orders[i]
             if isinstance(nus, np.ndarray):
@@ -119,31 +119,31 @@ def allocate_workspace_arrays(self, num_points):
         if self.current_max_num_points < num_points:
             self.current_max_num_points = num_points
             self.basis_workspace = np.empty((
-                self.ndim, 
+                self.xdim, 
                 self.current_max_num_points,
                 2*np.max(self.orders)+3
             ), dtype=np.float_)
-            self.interval_workspace = np.empty((self.ndim, self.current_max_num_points, ), dtype=np.intc)
+            self.interval_workspace = np.empty((self.xdim, self.current_max_num_points, ), dtype=np.intc)
             self.coefficient_selector = np.empty(self.coefficient_shape_base + (self.current_max_num_points,), dtype=np.int_)
 
     def __call__(self, x, nus=0):
         """
         Parameters
         ----------
-        x : ndarray, shape=(self.ndim, ...) dtype=np.float_
-            Point(s) to evaluate spline on. Output will be (self.mdim,...)
-        nus : ndarray, broadcastable to shape=(self.ndim,) dtype=np.int_
+        x : ndarray, shape=(self.xdim, ...) dtype=np.float_
+            Point(s) to evaluate spline on. Output will be (self.ydim,...)
+        nus : ndarray, broadcastable to shape=(self.xdim,) dtype=np.int_
             Order of derivative(s) for each dimension to evaluate
             
         """
         if not isinstance(x, np.ndarray):
             x = np.array(x)
         x_shape, x_ndim = x.shape, x.ndim
-        x = np.ascontiguousarray(x.reshape((self.ndim, -1)), dtype=np.float_)
+        x = np.ascontiguousarray(x.reshape((self.xdim, -1)), dtype=np.float_)
         num_points = x.shape[-1]
 
         if isinstance(nus, np.ndarray):
-            if nus.ndim != 1 or nus.size != self.ndim:
+            if nus.ndim != 1 or nus.size != self.xdim:
                 raise ValueError("nus is wrong shape")
 
         self.allocate_workspace_arrays(x.shape[-1])
@@ -153,8 +153,8 @@ def __call__(self, x, nus=0):
         y_out = np.einsum(self.coefficients[coefficient_selector], self.coefficient_op, 
             *self.u_arg, 
             self.output_op)
-        if (x_ndim > 1) and ((self.ndim > 1) or (size[0] != self.ndim)):
-            y_out = y_out.reshape((self.mdim,) + x_shape[1:])
+        if (x_ndim > 1) and ((self.xdim > 1) or (size[0] != self.xdim)):
+            y_out = y_out.reshape((self.ydim,) + x_shape[1:])
         return y_out
 
 def make_lsq_spline(x, y, knots, orders, w=None, check_finite=True):
@@ -163,13 +163,13 @@ def make_lsq_spline(x, y, knots, orders, w=None, check_finite=True):
 
     Parameters
     ----------
-    x : array_like, shape (ndim, num_points)
+    x : array_like, shape (xdim, num_points)
         Abscissas.
-    y : array_like, shape (mdim, num_points)
+    y : array_like, shape (ydim, num_points)
         Ordinates.
-    knots : iterable of array_like, shape (n_1 + orders[0] + 1,), ... (n_ndim, + orders[-1] + 1)
+    knots : iterable of array_like, shape (n_1 + orders[0] + 1,), ... (n_xdim, + orders[-1] + 1)
         Knots and data points must satisfy Schoenberg-Whitney conditions.
-    orders : ndarray, shape=(ndim,), dtype=np.int_
+    orders : ndarray, shape=(xdim,), dtype=np.int_
     w : array_like, shape (num_points,), optional
         Weights for spline fitting. Must be positive. If ``None``,
         then weights are all equal.
@@ -185,8 +185,8 @@ def make_lsq_spline(x, y, knots, orders, w=None, check_finite=True):
 
 
     """
-    ndim = x.shape[0]
-    mdim = y.shape[0]
+    xdim = x.shape[0]
+    ydim = y.shape[0]
     num_points = x.shape[1]
     assert x.shape[1] == y.shape[1]
     assert x.ndim == 2
@@ -197,7 +197,7 @@ def make_lsq_spline(x, y, knots, orders, w=None, check_finite=True):
 
     knot_shapes = tuple(knot.size - order - 1 for knot, order in zip(knots, orders))
     
-    temp_spline = BSplineNDInterpolator(knots, np.empty(mdim), orders)
+    temp_spline = BSplineNDInterpolator(knots, np.empty(ydim), orders)
     temp_spline.allocate_workspace_arrays(num_points)
     temp_spline.compute_basis_coefficient_selector(x)
 
@@ -210,7 +210,7 @@ def make_lsq_spline(x, y, knots, orders, w=None, check_finite=True):
     # TODO: implemnet weighting matrix, which I think is just matrix multiply by diag(w) on left for both observation matrix and output.
 
     lsq_coefficients, lsq_residuals, rank, singular_values = np.linalg.lstsq(observation_matrix, y.T)
-    temp_spline.coefficients = lsq_coefficients.T.reshape((mdim,) + knot_shapes )
+    temp_spline.coefficients = lsq_coefficients.T.reshape((ydim,) + knot_shapes )
 
     # TODO: I think people will want this matrix, is there a better way to give this to a user?
     temp_spline.observation_matrix = observation_matrix
@@ -223,15 +223,15 @@ def make_interp_spline(x, y, bcs=0, orders=3):
 
     Parameters
     ----------
-    x : array_like, shape (ndim, n_1, n_2, ..., n_ndim) or arguments to np.meshgrid to construct same
-        Abscissas.
-    y : array_like, shape (mdim, n_1, n_2, ..., n_ndim)
-        Ordinates.
+    x : array_like broadcastable to (xdim, n_1, n_2, ..., n_xdim) or arguments to np.meshgrid to construct same
+        Abscissas. 
+    y : array_like, 
+        Ordinates. shape (ydim, n_1, n_2, ..., n_xdim)
     bcs : (list of) 2-tuples or None
-    orders : ndarray, shape=(ndim,), dtype=np.int_
+    orders : ndarray, shape=(xdim,), dtype=np.int_
         Degree of interpolant for each axis (or broadcastable)
-    periodic : ndarray, shape=(ndim,), dtype=np.bool_
-    extrapolate : ndarray, shape=(ndim,), dtype=np.bool_
+    periodic : ndarray, shape=(xdim,), dtype=np.bool_
+    extrapolate : ndarray, shape=(xdim,), dtype=np.bool_
 
     Notes
     -----
@@ -241,32 +241,32 @@ def make_interp_spline(x, y, bcs=0, orders=3):
     """
     if isinstance(x, np.ndarray): # mesh
         if x.ndim == 1:
-            ndim = 1
+            xdim = 1
             x = x[None,...]
         else:
-            ndim = x.ndim - 1
+            xdim = x.ndim - 1
 
     elif not isinstance(x, str) and len(x): # vectors, or try
-        ndim = len(x)
+        xdim = len(x)
         x = np.meshgrid(x, indexing='ij')
         
     else:
         raise ValueError("Don't know how to interpret x")
     
-    if y.ndim == ndim:
+    if y.ndim == xdim:
         # how can you tell if y needs a [None, ...] or [...] ?
-        # same question as to whether ndim is y.ndim or y.ndim-1
+        # same question as to whether xdim is y.ndim or y.ndim-1
         # I think this is the right answer. 
         y = y[None, ...]
-    elif y.ndim == ndim+1:
+    elif y.ndim == xdim+1:
         pass
     else:
         raise ValueError("incompatible dimension size")
-    mdim = y.shape[0]
-    # generally, x.shape = (ndim, n1, n2, ..., n_ndim)
-    # and y.sahpe = (mdim, n1, n2, ..., n_ndim)
-    bcs = np.broadcast_to(bcs, (ndim,2))
-    orders = np.broadcast_to(orders, (ndim,))
+    ydim = y.shape[0]
+    # generally, x.shape = (xdim, n1, n2, ..., n_xdim)
+    # and y.sahpe = (ydim, n1, n2, ..., n_xdim)
+    bcs = np.broadcast_to(bcs, (xdim,2))
+    orders = np.broadcast_to(orders, (xdim,))
     deriv_specs = np.asarray((bcs[:,:]>0),dtype=np.int)
 
     knots = []
@@ -275,11 +275,11 @@ def make_interp_spline(x, y, bcs=0, orders=3):
     axis=0
     check_finite=True
 
-    for i in np.arange(ndim)+1:
+    for i in np.arange(xdim)+1:
         all_other_ax_shape = np.asarray(np.r_[coefficients.shape[1:i],
             y.shape[i+1:]], dtype=np.int)
         x_line_sel = ((i-1,) + (0,)*(i-1) + (slice(None,None),) +
-            (0,)*(ndim-i))
+            (0,)*(xdim-i))
         x_slice = x[x_line_sel]
         k = orders[i-1]
 

From d66a435ddd5d79c027542a2cef12c654a7a94f58 Mon Sep 17 00:00:00 2001
From: ben <ben@sixpearls.com>
Date: Sun, 16 Jun 2019 11:00:48 -0700
Subject: [PATCH 05/10] breaking up algorithm, separating out BCs, etc

---
 ndsplines/ndsplines.py | 101 ++++++++++++++++++++---------------------
 1 file changed, 49 insertions(+), 52 deletions(-)

diff --git a/ndsplines/ndsplines.py b/ndsplines/ndsplines.py
index a66098f..db14492 100644
--- a/ndsplines/ndsplines.py
+++ b/ndsplines/ndsplines.py
@@ -17,7 +17,7 @@
 else:
     default_implementation = _bspl
 
-__all__ = ['pinned', 'clamped', 'extrap', 'periodic', 'BSplineNDInterpolator',
+__all__ = ['pinned', 'clamped', 'notaknot', 'BSplineNDInterpolator',
            'make_interp_spline', 'make_lsq_spline', 'default_implementation']
 
 
@@ -34,12 +34,10 @@
 
 
 """
-pinned = 2
-clamped = 1
-extrap = 0
-periodic = -1
+clamped = np.array([1,0.0])
+pinned = np.array([2,0.0])
+notaknot = np.array([0,0.0])
 
-bc_map =  {clamped: "clamped", pinned: "natural", extrap: None, periodic: None}
 
 
 class BSplineNDInterpolator(object):
@@ -265,9 +263,11 @@ def make_interp_spline(x, y, bcs=0, orders=3):
     ydim = y.shape[0]
     # generally, x.shape = (xdim, n1, n2, ..., n_xdim)
     # and y.sahpe = (ydim, n1, n2, ..., n_xdim)
-    bcs = np.broadcast_to(bcs, (xdim,2))
     orders = np.broadcast_to(orders, (xdim,))
-    deriv_specs = np.asarray((bcs[:,:]>0),dtype=np.int)
+
+    bcs = np.broadcast_to(bcs, (xdim,2,2))
+    deriv_specs = np.asarray((bcs[:,:,0]>0),dtype=np.int)
+    nak_spec = np.asarray((bcs[:,:,0]==0),dtype=np.bool)
 
     knots = []
     coefficients = np.pad(y, np.r_[np.c_[0,0], deriv_specs], 'constant')
@@ -283,23 +283,11 @@ def make_interp_spline(x, y, bcs=0, orders=3):
         x_slice = x[x_line_sel]
         k = orders[i-1]
 
-        bc_type=(bc_map[(bcs[i-1,0])],
-                 bc_map[(bcs[i-1,1])])
-
         t = None
 
         #=======================================
         # START FROM SCIPY
         #=======================================
-        if bc_type is None or bc_type == 'not-a-knot':
-            deriv_l, deriv_r = None, None
-        elif isinstance(bc_type, string_types):
-            deriv_l, deriv_r = bc_type, bc_type
-        else:
-            try:
-                deriv_l, deriv_r = bc_type
-            except TypeError:
-                raise ValueError("Unknown boundary condition: %s" % bc_type)
 
         if k == 0:
             if any(_ is not None for _ in (t, deriv_l, deriv_r)):
@@ -316,7 +304,7 @@ def make_interp_spline(x, y, bcs=0, orders=3):
 
         # come up with a sensible knot vector, if needed
         if t is None:
-            if deriv_l is None and deriv_r is None:
+            if np.all(nak_spec[i-1, :]):
                 if k == 2:
                     # OK, it's a bit ad hoc: Greville sites + omit
                     # 2nd and 2nd-to-last points, a la not-a-knot
@@ -331,6 +319,45 @@ def make_interp_spline(x, y, bcs=0, orders=3):
 
         t = _as_float_array(t, check_finite)
 
+        # Here : deriv_l, r = [(nu, value), ...]
+        deriv_l_ords, deriv_r_ords = bcs[i-1, :, 0].astype(np.int_)
+    
+        if deriv_l_ords == 0:
+            deriv_l_ords = np.array([])
+            deriv_l_vals = np.array([])
+            nleft = 0
+        else:
+            deriv_l_ords = np.array([deriv_l_ords])
+            deriv_l_vals = np.broadcast_to(bcs[i-1, 0, 1], ydim)
+            nleft = 1
+
+        if deriv_r_ords == 0:
+            deriv_r_ords = np.array([])
+            deriv_r_vals = np.array([])
+            nright = 0
+        else:
+            deriv_r_ords = np.array([deriv_r_ords])
+            deriv_r_vals = np.broadcast_to(bcs[i-1, 1, 1], ydim)
+            nright = 1
+
+        # have `n` conditions for `nt` coefficients; need nt-n derivatives
+        n = x_slice.size
+        nt = t.size - k - 1
+
+        if nt - n != nleft + nright:
+            raise ValueError("The number of derivatives at boundaries does not "
+                             "match: expected %s, got %s+%s" % (nt-n, nleft, nright))
+
+        # set up the LHS: the collocation matrix + derivatives at boundaries
+        kl = ku = k
+        ab = np.zeros((2*kl + ku + 1, nt), dtype=np.float_, order='F')
+        _sci_bspl._colloc(x_slice, t, k, ab, offset=nleft)
+        if nleft > 0:
+            _sci_bspl._handle_lhs_derivatives(t, k, x_slice[0], ab, kl, ku, deriv_l_ords)
+        if nright > 0:
+            _sci_bspl._handle_lhs_derivatives(t, k, x_slice[-1], ab, kl, ku, deriv_r_ords,
+                                    offset=nt-nright)
+
         #=======================================
         # END FROM SCIPY
         #=======================================
@@ -387,34 +414,6 @@ def make_interp_spline(x, y, bcs=0, orders=3):
             if (x_slice[0] < t[k]) or (x_slice[-1] > t[-k]):
                 raise ValueError('Out of bounds w/ x_slice = %s.' % x_slice)
 
-            # Here : deriv_l, r = [(nu, value), ...]
-            deriv_l = _convert_string_aliases(deriv_l, y_slice.shape[1:])
-            deriv_l_ords, deriv_l_vals = _process_deriv_spec(deriv_l)
-            nleft = deriv_l_ords.shape[0]
-
-            deriv_r = _convert_string_aliases(deriv_r, y_slice.shape[1:])
-            deriv_r_ords, deriv_r_vals = _process_deriv_spec(deriv_r)
-            nright = deriv_r_ords.shape[0]
-
-            # have `n` conditions for `nt` coefficients; need nt-n derivatives
-            n = x_slice.size
-            nt = t.size - k - 1
-
-            if nt - n != nleft + nright:
-                raise ValueError("The number of derivatives at boundaries does not "
-                                 "match: expected %s, got %s+%s" % (nt-n, nleft, nright))
-
-
-            # set up the LHS: the collocation matrix + derivatives at boundaries
-            kl = ku = k
-            ab = np.zeros((2*kl + ku + 1, nt), dtype=np.float_, order='F')
-            _sci_bspl._colloc(x_slice, t, k, ab, offset=nleft)
-            if nleft > 0:
-                _sci_bspl._handle_lhs_derivatives(t, k, x_slice[0], ab, kl, ku, deriv_l_ords)
-            if nright > 0:
-                _sci_bspl._handle_lhs_derivatives(t, k, x_slice[-1], ab, kl, ku, deriv_r_ords,
-                                        offset=nt-nright)
-
             # set up the RHS: values to interpolate (+ derivative values, if any)
             extradim = prod(y_slice.shape[1:])
             rhs = np.empty((nt, extradim), dtype=y_slice.dtype)
@@ -443,6 +442,4 @@ def make_interp_spline(x, y, bcs=0, orders=3):
             #=======================================
 
             coefficients[coeff_line_sel] = c.T
-    return BSplineNDInterpolator(knots, coefficients, orders, 
-        np.all(bcs==periodic, axis=1),
-        (bcs >= 0))
+    return BSplineNDInterpolator(knots, coefficients, orders,)

From fd02edc77b437891caa3e5355df2c1fde0cd11fe Mon Sep 17 00:00:00 2001
From: ben <ben@sixpearls.com>
Date: Mon, 17 Jun 2019 20:04:56 -0700
Subject: [PATCH 06/10] this seems to give same results as scipy.bspline for
 k=0,1,2,3.

---
 examples/1d-interp.py  |  71 ++++++++------
 ndsplines/ndsplines.py | 207 ++++++++++++++++++++++++++---------------
 2 files changed, 174 insertions(+), 104 deletions(-)

diff --git a/examples/1d-interp.py b/examples/1d-interp.py
index ce6cac6..ade8a9e 100644
--- a/examples/1d-interp.py
+++ b/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,57 @@ 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),
+    itertools.product(["natural", "clamped", ], repeat=2),
     ((None,None),),
 )))
 
-NDspline_dict = {"natural": ndsplines.pinned, "clamped": ndsplines.clamped, None: 0}
+# test_bcs = (None, (None,None,), "not-a-knot", ("not-a-knot", "not-a-knot"),)
 
-for func in funcs:
-    fvals = func(x)
-    truef = func(xx)
-    plt.figure()
-
-    plot_sel = slice(None)
-
-    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_dict = {"natural": ndsplines.pinned, "clamped": ndsplines.clamped, None: ndsplines.notaknot}
 
+for order in range(4):
+    for func in funcs:
+        fvals = func(x)
+        truef = func(xx)
+        plt.figure()
     
-    plt.plot(xx, truef, 'k--', label="True " + func.__name__)
+        plot_sel = slice(None)
     
-    plt.legend(loc='best')
-    plt.plot(x, fvals, 'o')
-    plt.show()
+        plt.gca().set_prop_cycle(None)
+    
+        for test_bc in 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)')
+            # break
+    
+        plt.gca().set_prop_cycle(None)
+    
+        for test_bc in test_bcs[plot_sel,:]:
+            try:
+                test_NDBspline = ndsplines.make_interp_spline(x, fvals, bcs=(NDspline_dict[test_bc[0]], NDspline_dict[test_bc[1]]), orders=order)
+            except ValueError:
+                continue
+            else:
+                NDsplienf = test_NDBspline(xx.copy())
+                plt.plot(xx, NDsplienf[0], label=str(test_bc) + ' (ndspline)' )
+            # break
+    
+        
+        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()
     
diff --git a/ndsplines/ndsplines.py b/ndsplines/ndsplines.py
index db14492..8e519aa 100644
--- a/ndsplines/ndsplines.py
+++ b/ndsplines/ndsplines.py
@@ -214,6 +214,27 @@ def make_lsq_spline(x, y, knots, orders, w=None, check_finite=True):
     temp_spline.observation_matrix = observation_matrix
     return temp_spline
 
+def _not_a_knot(x, k, left=True, right=True):
+    if k > 2 and k % 2 == 0:
+        raise ValueError("Odd degree for now only. Got %s." % k)
+    if k==0: # special case for k = 0, only apply to right
+        left = False
+    t = x
+    if left:
+        t = np.r_[(t[0],)*(k+1), t[(k -1) //2 +1:]]
+    if right:
+        t = np.r_[t[:-(k-1)//2 -1 or None], (t[-1],)*(k+1)]
+    return t
+
+def _augknt(x, k, left=True, right=True):
+    t = x
+    if left:
+        t = np.r_[(t[0],)*k, t]
+    if right:
+        t = np.r_[t, (t[-1],)*k]
+    return t
+
+
 
 def make_interp_spline(x, y, bcs=0, orders=3):
     """
@@ -226,6 +247,10 @@ def make_interp_spline(x, y, bcs=0, orders=3):
     y : array_like, 
         Ordinates. shape (ydim, n_1, n_2, ..., n_xdim)
     bcs : (list of) 2-tuples or None
+        2-tuple defines for each side a 2-tuple of (deriv_spec, spec_value)
+        Use deriv_spec == 0 for not-a-knot boundary condition
+        For k=0, both spec_values = 0 implements nearest-neighbor,
+        a single side with spec_value = 0 uses zero-order-hold from that direction
     orders : ndarray, shape=(xdim,), dtype=np.int_
         Degree of interpolant for each axis (or broadcastable)
     periodic : ndarray, shape=(xdim,), dtype=np.bool_
@@ -233,9 +258,7 @@ def make_interp_spline(x, y, bcs=0, orders=3):
 
     Notes
     -----
-    TODO: use scipy source to implement this more efficiently?
-    i.e., do knot computation once for each dimension, then coefficients for
-    each line 
+    Special case boundary condition - for k=0, 
     """
     if isinstance(x, np.ndarray): # mesh
         if x.ndim == 1:
@@ -267,7 +290,7 @@ def make_interp_spline(x, y, bcs=0, orders=3):
 
     bcs = np.broadcast_to(bcs, (xdim,2,2))
     deriv_specs = np.asarray((bcs[:,:,0]>0),dtype=np.int)
-    nak_spec = np.asarray((bcs[:,:,0]==0),dtype=np.bool)
+    nak_spec = np.asarray((bcs[:,:,0]<=0),dtype=np.bool)
 
     knots = []
     coefficients = np.pad(y, np.r_[np.c_[0,0], deriv_specs], 'constant')
@@ -283,44 +306,76 @@ def make_interp_spline(x, y, bcs=0, orders=3):
         x_slice = x[x_line_sel]
         k = orders[i-1]
 
+        left_nak, right_nak = nak_spec[i-1, :]
+        both_nak = left_nak and right_nak
+
+        # Here : deriv_l, r = [(nu, value), ...]
+        deriv_l_ords, deriv_r_ords = bcs[i-1, :, 0].astype(np.int_)
+
         t = None
 
         #=======================================
         # START FROM SCIPY
         #=======================================
 
+        # should there be a general check for k <= deriv_spec ?
+
         if k == 0:
-            if any(_ is not None for _ in (t, deriv_l, deriv_r)):
+            # all derivatives are fully defined, can only set 0-th derivative,
+            # aka not-a-knot boundary conditions to both sides
+            if not both_nak:
                 raise ValueError("Too much info for k=0: t and bc_type can only "
-                                 "be None.")
+                                 "be notaknot.")
             x_slice = _as_float_array(x_slice, check_finite)
-            t = np.r_[x_slice, x_slice[-1]]
 
-        if k == 1 and t is None:
-            if not (deriv_l is None and deriv_r is None):
-                raise ValueError("Too much info for k=1: bc_type can only be None.")
+            # not the same as not a knot... _not_a_knot would apply to both ends,
+            # this only applies to right. added as special case to 1-sided
+            # _not_a_knot implementation.
+            left_zero, right_zero = (bcs[i-1, :, 1]==0)
+
+            if left_zero and right_zero:
+                t = np.r_[x_slice[0], (x_slice[:-1] + x_slice[1:])/2., x_slice[-1]]
+            elif not left_zero and right_zero:
+                t = np.r_[x_slice, x_slice[-1]]
+            elif left_zero == 0 and not right_zero:
+                t = np.r_[x_slice[0], x_slice]
+            else:
+                raise ValueError("Set deriv_spec = 0, with up to one side = -1 for k=0")
+            
+
+            # actually, I suspect which side it is applied to determines
+            # whether it is causal or anti-causal zero-order hold.
+            # greville sites + x[0], x[-1] may give nearest neighbor?
+
+        if k == 1:
+            # all derivatives are fully defined, can only set 0-th derivative,
+            # aka not-a-knot boundary conditions to both sides
+            if not both_nak:
+                raise ValueError("Too much info for k=1: bc_type can only be notaknot.")
             x_slice = _as_float_array(x_slice, check_finite)
-            t = np.r_[x_slice[0], x_slice, x_slice[-1]]
-
-        # come up with a sensible knot vector, if needed
-        if t is None:
-            if np.all(nak_spec[i-1, :]):
-                if k == 2:
-                    # OK, it's a bit ad hoc: Greville sites + omit
-                    # 2nd and 2nd-to-last points, a la not-a-knot
-                    t = (x_slice[1:] + x_slice[:-1]) / 2.
-                    t = np.r_[(x_slice[0],)*(k+1),
-                               t[1:-1],
-                               (x_slice[-1],)*(k+1)]
-                else:
-                    t = _not_a_knot(x_slice, k)
+            # apply _not_a_knot
+
+        if k==2:
+
+            if both_nak:
+                # special, ad-hoc case using greville sites
+                t = (x_slice[1:] + x_slice[:-1]) / 2.
+                t = np.r_[(x_slice[0],)*(k+1),
+                           t[1:-1],
+                           (x_slice[-1],)*(k+1)]
+
+            elif left_nak or right_nak:
+                raise ValueError("For k=2, can set both sides or neither side to notaknot")
             else:
-                t = _augknt(x_slice, k)
+                t = x_slice
 
-        t = _as_float_array(t, check_finite)
+        elif k != 0:
+            t = _not_a_knot(x_slice, k, left_nak, right_nak)
 
-        # Here : deriv_l, r = [(nu, value), ...]
-        deriv_l_ords, deriv_r_ords = bcs[i-1, :, 0].astype(np.int_)
+        t = _augknt(t, k, not left_nak, not right_nak)
+
+        t = _as_float_array(t, check_finite)
+        
     
         if deriv_l_ords == 0:
             deriv_l_ords = np.array([])
@@ -343,8 +398,9 @@ def make_interp_spline(x, y, bcs=0, orders=3):
         # have `n` conditions for `nt` coefficients; need nt-n derivatives
         n = x_slice.size
         nt = t.size - k - 1
-
-        if nt - n != nleft + nright:
+        
+        # this also catches if deriv_spec > k-1, possibly?
+        if np.clip(nt - n, 0, np.inf).astype(int) != nleft + nright:
             raise ValueError("The number of derivatives at boundaries does not "
                              "match: expected %s, got %s+%s" % (nt-n, nleft, nright))
 
@@ -384,58 +440,59 @@ def make_interp_spline(x, y, bcs=0, orders=3):
 
             # special-case k=0 right away
             if k == 0:
-                if any(_ is not None for _ in (t, deriv_l, deriv_r)):
+                if not both_nak:
                     raise ValueError("Too much info for k=0: t and bc_type can only "
-                                     "be None.")
+                                     "be notaknot.")
                 c = np.asarray(y_slice)
 
             # special-case k=1 (e.g., Lyche and Morken, Eq.(2.16))
-            if k == 1 and t is None:
-                if not (deriv_l is None and deriv_r is None):
+            elif k == 1:
+                if not both_nak:
                     raise ValueError("Too much info for k=1: bc_type can only be None.")
                 c = np.asarray(y_slice)
 
-            x_slice = _as_float_array(x_slice, check_finite)
-            y_slice = _as_float_array(y_slice, check_finite)
-            k = operator.index(k)
-
-
-            if x_slice.ndim != 1 or np.any(x_slice[1:] <= x_slice[:-1]):
-                raise ValueError("Expect x_slice to be a 1-D sorted array_like.")
-            if k < 0:
-                raise ValueError("Expect non-negative k.")
-            if t.ndim != 1 or np.any(t[1:] < t[:-1]):
-                raise ValueError("Expect t to be a 1-D sorted array_like.")
-            if x_slice.size != y_slice.shape[0]:
-                raise ValueError('x_slice and y_slice are incompatible.')
-            if t.size < x_slice.size + k + 1:
-                raise ValueError('Got %d knots, need at least %d.' %
-                                 (t.size, x_slice.size + k + 1))
-            if (x_slice[0] < t[k]) or (x_slice[-1] > t[-k]):
-                raise ValueError('Out of bounds w/ x_slice = %s.' % x_slice)
-
-            # set up the RHS: values to interpolate (+ derivative values, if any)
-            extradim = prod(y_slice.shape[1:])
-            rhs = np.empty((nt, extradim), dtype=y_slice.dtype)
-            if nleft > 0:
-                rhs[:nleft] = deriv_l_vals.reshape(-1, extradim)
-            rhs[nleft:nt - nright] = y_slice.reshape(-1, extradim)
-            if nright > 0:
-                rhs[nt - nright:] = deriv_r_vals.reshape(-1, extradim)
-
-            # solve Ab @ x_slice = rhs; this is the relevant part of linalg.solve_banded
-            if check_finite:
-                ab, rhs = map(np.asarray_chkfinite, (ab, rhs))
-            gbsv, = get_lapack_funcs(('gbsv',), (ab, rhs))
-            lu, piv, c, info = gbsv(kl, ku, ab, rhs,
-                    overwrite_ab=True, overwrite_b=True)
-
-            if info > 0:
-                raise LinAlgError("Collocation matix is singular.")
-            elif info < 0:
-                raise ValueError('illegal value in %d-th argument of internal gbsv' % -info)
-
-            c = np.ascontiguousarray(c.reshape((nt,) + y_slice.shape[1:]))
+            else:
+                x_slice = _as_float_array(x_slice, check_finite)
+                y_slice = _as_float_array(y_slice, check_finite)
+                k = operator.index(k)
+
+
+                if x_slice.ndim != 1 or np.any(x_slice[1:] <= x_slice[:-1]):
+                    raise ValueError("Expect x_slice to be a 1-D sorted array_like.")
+                if k < 0:
+                    raise ValueError("Expect non-negative k.")
+                if t.ndim != 1 or np.any(t[1:] < t[:-1]):
+                    raise ValueError("Expect t to be a 1-D sorted array_like.")
+                if x_slice.size != y_slice.shape[0]:
+                    raise ValueError('x_slice and y_slice are incompatible.')
+                if t.size < x_slice.size + k + 1:
+                    raise ValueError('Got %d knots, need at least %d.' %
+                                     (t.size, x_slice.size + k + 1))
+                if (x_slice[0] < t[k]) or (x_slice[-1] > t[-k]):
+                    raise ValueError('Out of bounds w/ x_slice = %s.' % x_slice)
+
+                # set up the RHS: values to interpolate (+ derivative values, if any)
+                extradim = prod(y_slice.shape[1:])
+                rhs = np.empty((nt, extradim), dtype=y_slice.dtype)
+                if nleft > 0:
+                    rhs[:nleft] = deriv_l_vals.reshape(-1, extradim)
+                rhs[nleft:nt - nright] = y_slice.reshape(-1, extradim)
+                if nright > 0:
+                    rhs[nt - nright:] = deriv_r_vals.reshape(-1, extradim)
+
+                # solve Ab @ x_slice = rhs; this is the relevant part of linalg.solve_banded
+                if check_finite:
+                    ab, rhs = map(np.asarray_chkfinite, (ab, rhs))
+                gbsv, = get_lapack_funcs(('gbsv',), (ab, rhs))
+                lu, piv, c, info = gbsv(kl, ku, ab, rhs,
+                        overwrite_ab=True, overwrite_b=True)
+
+                if info > 0:
+                    raise LinAlgError("Collocation matix is singular.")
+                elif info < 0:
+                    raise ValueError('illegal value in %d-th argument of internal gbsv' % -info)
+
+                c = np.ascontiguousarray(c.reshape((nt,) + y_slice.shape[1:]))
 
             #=======================================
             # START FROM SCIPY

From a11235f07701eeb5cd61e04f95e03c94a30581e1 Mon Sep 17 00:00:00 2001
From: ben <ben@sixpearls.com>
Date: Mon, 17 Jun 2019 22:43:29 -0700
Subject: [PATCH 07/10] able to use not-a-knot on each side for k>2, k odd.
 decent API for nearest-neighbor vs left zoh vs right zoh

---
 examples/1d-interp.py  | 38 ++++++++++++++++++++++++++++----------
 ndsplines/ndsplines.py |  8 ++++----
 2 files changed, 32 insertions(+), 14 deletions(-)

diff --git a/examples/1d-interp.py b/examples/1d-interp.py
index ade8a9e..6d3ba62 100644
--- a/examples/1d-interp.py
+++ b/examples/1d-interp.py
@@ -33,14 +33,30 @@ def tanh(x_in):
 k = 3
 
 
-test_bcs = np.array(list(itertools.chain(
+scipy_test_bcs = np.array(list(itertools.chain(
     itertools.product(["natural", "clamped", ], repeat=2),
     ((None,None),),
 )))
 
+k0_bcs = np.array(list(
+itertools.product([(0,0), (0,-1)], repeat=2)
+)[:-1])[[-1,0,1]]
+
+NDspline_dict = {"natural": ndsplines.pinned, "clamped": ndsplines.clamped, "not-a-knot": ndsplines.notaknot}
+
+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"]),
+)])
+
+
+NDspline_bc_to_string = {tuple(v):k for k,v in NDspline_dict.items()}
+NDspline_bc_to_string[(0,-1)] = 'one-sided hold'
+
 # test_bcs = (None, (None,None,), "not-a-knot", ("not-a-knot", "not-a-knot"),)
 
-NDspline_dict = {"natural": ndsplines.pinned, "clamped": ndsplines.clamped, None: ndsplines.notaknot}
 
 for order in range(4):
     for func in funcs:
@@ -52,7 +68,7 @@ def tanh(x_in):
     
         plt.gca().set_prop_cycle(None)
     
-        for test_bc in test_bcs[plot_sel,:]:
+        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:
@@ -60,19 +76,21 @@ def tanh(x_in):
             else:
                 splinef = test_Bspline(xx.copy(), extrapolate=True)
                 plt.plot(xx, splinef, '--', lw=3.0, label=str(test_bc) + ' (BSpline)')
-            # break
-    
+
         plt.gca().set_prop_cycle(None)
-    
-        for test_bc in test_bcs[plot_sel,:]:
+        
+        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=(NDspline_dict[test_bc[0]], NDspline_dict[test_bc[1]]), orders=order)
+                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=str(test_bc) + ' (ndspline)' )
-            # break
+                plt.plot(xx, NDsplienf[0], label=', '.join(*tuple([NDspline_bc_to_string[tuple(bc)] for bc in test_bc])) + ' (ndspline)' )
     
         
         plt.plot(xx, truef, 'k--', label="True " + func.__name__)
diff --git a/ndsplines/ndsplines.py b/ndsplines/ndsplines.py
index 8e519aa..416d3e4 100644
--- a/ndsplines/ndsplines.py
+++ b/ndsplines/ndsplines.py
@@ -337,7 +337,7 @@ def make_interp_spline(x, y, bcs=0, orders=3):
                 t = np.r_[x_slice[0], (x_slice[:-1] + x_slice[1:])/2., x_slice[-1]]
             elif not left_zero and right_zero:
                 t = np.r_[x_slice, x_slice[-1]]
-            elif left_zero == 0 and not right_zero:
+            elif left_zero and not right_zero:
                 t = np.r_[x_slice[0], x_slice]
             else:
                 raise ValueError("Set deriv_spec = 0, with up to one side = -1 for k=0")
@@ -377,7 +377,7 @@ def make_interp_spline(x, y, bcs=0, orders=3):
         t = _as_float_array(t, check_finite)
         
     
-        if deriv_l_ords == 0:
+        if nak_spec[i-1,0]:
             deriv_l_ords = np.array([])
             deriv_l_vals = np.array([])
             nleft = 0
@@ -386,7 +386,7 @@ def make_interp_spline(x, y, bcs=0, orders=3):
             deriv_l_vals = np.broadcast_to(bcs[i-1, 0, 1], ydim)
             nleft = 1
 
-        if deriv_r_ords == 0:
+        if nak_spec[i-1,1]:
             deriv_r_ords = np.array([])
             deriv_r_vals = np.array([])
             nright = 0
@@ -425,7 +425,7 @@ def make_interp_spline(x, y, bcs=0, orders=3):
             offset_axes_remaining_sel = (tuple(idx[i-1:] + 
                 deriv_specs[i:,0]))
             y_line_sel = ((Ellipsis,) + idx[:i-1] + 
-                (slice(deriv_specs[i-1,0],-deriv_specs[i-1,0] or None),) +
+                (slice(deriv_specs[i-1,0],-deriv_specs[i-1,1] or None),) +
                 offset_axes_remaining_sel)
             coeff_line_sel = ((Ellipsis,) + idx[:i-1] + (slice(None,None),)
                 + offset_axes_remaining_sel)

From 31f328d1f485ae5a9b20e914f2c297ca38da8c2f Mon Sep 17 00:00:00 2001
From: ben <ben@sixpearls.com>
Date: Tue, 18 Jun 2019 08:31:31 -0700
Subject: [PATCH 08/10] cleaning up

---
 examples/1d-interp.py  |  5 +--
 ndsplines/ndsplines.py | 82 +++++++++++++-----------------------------
 2 files changed, 25 insertions(+), 62 deletions(-)

diff --git a/examples/1d-interp.py b/examples/1d-interp.py
index 6d3ba62..8349264 100644
--- a/examples/1d-interp.py
+++ b/examples/1d-interp.py
@@ -55,9 +55,6 @@ def tanh(x_in):
 NDspline_bc_to_string = {tuple(v):k for k,v in NDspline_dict.items()}
 NDspline_bc_to_string[(0,-1)] = 'one-sided hold'
 
-# test_bcs = (None, (None,None,), "not-a-knot", ("not-a-knot", "not-a-knot"),)
-
-
 for order in range(4):
     for func in funcs:
         fvals = func(x)
@@ -90,7 +87,7 @@ def tanh(x_in):
                 continue
             else:
                 NDsplienf = test_NDBspline(xx.copy())
-                plt.plot(xx, NDsplienf[0], label=', '.join(*tuple([NDspline_bc_to_string[tuple(bc)] for bc in test_bc])) + ' (ndspline)' )
+                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__)
diff --git a/ndsplines/ndsplines.py b/ndsplines/ndsplines.py
index 416d3e4..cfd6e1b 100644
--- a/ndsplines/ndsplines.py
+++ b/ndsplines/ndsplines.py
@@ -4,9 +4,8 @@
 from scipy._lib.six import string_types
 from scipy.linalg import (get_lapack_funcs, LinAlgError,
                           cholesky_banded, cho_solve_banded)
-from scipy.interpolate._bsplines import (_as_float_array, _get_dtype, _augknt, 
-    _convert_string_aliases, _process_deriv_spec, _bspl as _sci_bspl, prod, 
-    _not_a_knot)
+from scipy.interpolate._bsplines import (_as_float_array, _bspl as _sci_bspl, 
+    prod,)
 
 from ndsplines import _npy_bspl
 
@@ -312,25 +311,17 @@ def make_interp_spline(x, y, bcs=0, orders=3):
         # Here : deriv_l, r = [(nu, value), ...]
         deriv_l_ords, deriv_r_ords = bcs[i-1, :, 0].astype(np.int_)
 
-        t = None
-
-        #=======================================
-        # START FROM SCIPY
-        #=======================================
-
+        x_slice = _as_float_array(x_slice, check_finite)
         # should there be a general check for k <= deriv_spec ?
 
         if k == 0:
             # all derivatives are fully defined, can only set 0-th derivative,
-            # aka not-a-knot boundary conditions to both sides
+            # special case for nearest-neighbor, causal/anti-causal zero-order
+            # hold
             if not both_nak:
                 raise ValueError("Too much info for k=0: t and bc_type can only "
                                  "be notaknot.")
-            x_slice = _as_float_array(x_slice, check_finite)
 
-            # not the same as not a knot... _not_a_knot would apply to both ends,
-            # this only applies to right. added as special case to 1-sided
-            # _not_a_knot implementation.
             left_zero, right_zero = (bcs[i-1, :, 1]==0)
 
             if left_zero and right_zero:
@@ -341,22 +332,16 @@ def make_interp_spline(x, y, bcs=0, orders=3):
                 t = np.r_[x_slice[0], x_slice]
             else:
                 raise ValueError("Set deriv_spec = 0, with up to one side = -1 for k=0")
-            
-
-            # actually, I suspect which side it is applied to determines
-            # whether it is causal or anti-causal zero-order hold.
-            # greville sites + x[0], x[-1] may give nearest neighbor?
 
+        # special-case k=1 (e.g., Lyche and Morken, Eq.(2.16))
         if k == 1:
             # all derivatives are fully defined, can only set 0-th derivative,
             # aka not-a-knot boundary conditions to both sides
             if not both_nak:
                 raise ValueError("Too much info for k=1: bc_type can only be notaknot.")
-            x_slice = _as_float_array(x_slice, check_finite)
-            # apply _not_a_knot
 
         if k==2:
-
+            # it's possible this may be the best behavior for all even k > 0
             if both_nak:
                 # special, ad-hoc case using greville sites
                 t = (x_slice[1:] + x_slice[:-1]) / 2.
@@ -377,21 +362,21 @@ def make_interp_spline(x, y, bcs=0, orders=3):
         t = _as_float_array(t, check_finite)
         
     
-        if nak_spec[i-1,0]:
+        if left_nak:
             deriv_l_ords = np.array([])
             deriv_l_vals = np.array([])
             nleft = 0
         else:
-            deriv_l_ords = np.array([deriv_l_ords])
+            deriv_l_ords = np.array([bcs[i-1, 0, 0]], dtype=np.int_)
             deriv_l_vals = np.broadcast_to(bcs[i-1, 0, 1], ydim)
             nleft = 1
 
-        if nak_spec[i-1,1]:
+        if right_nak:
             deriv_r_ords = np.array([])
             deriv_r_vals = np.array([])
             nright = 0
         else:
-            deriv_r_ords = np.array([deriv_r_ords])
+            deriv_r_ords = np.array([bcs[i-1, 1, 0]], dtype=np.int_)
             deriv_r_vals = np.broadcast_to(bcs[i-1, 1, 1], ydim)
             nright = 1
 
@@ -414,12 +399,21 @@ def make_interp_spline(x, y, bcs=0, orders=3):
             _sci_bspl._handle_lhs_derivatives(t, k, x_slice[-1], ab, kl, ku, deriv_r_ords,
                                     offset=nt-nright)
 
-        #=======================================
-        # END FROM SCIPY
-        #=======================================
-
         knots.append(t)
 
+        if k >= 2:
+            if x_slice.ndim != 1 or np.any(x_slice[1:] <= x_slice[:-1]):
+                raise ValueError("Expect x_slice to be a 1-D sorted array_like.")
+            if k < 0:
+                raise ValueError("Expect non-negative k.")
+            if t.ndim != 1 or np.any(t[1:] < t[:-1]):
+                raise ValueError("Expect t to be a 1-D sorted array_like.")
+            if t.size < x_slice.size + k + 1:
+                raise ValueError('Got %d knots, need at least %d.' %
+                                 (t.size, x_slice.size + k + 1))
+            if (x_slice[0] < t[k]) or (x_slice[-1] > t[-k]):
+                raise ValueError('Out of bounds w/ x_slice = %s.' % x_slice)
+
 
         for idx in np.ndindex(*all_other_ax_shape):
             offset_axes_remaining_sel = (tuple(idx[i-1:] + 
@@ -432,44 +426,20 @@ def make_interp_spline(x, y, bcs=0, orders=3):
 
             y_slice = coefficients[y_line_sel].T
 
-
-
-            #=======================================
-            # START FROM SCIPY
-            #=======================================
-
             # special-case k=0 right away
             if k == 0:
-                if not both_nak:
-                    raise ValueError("Too much info for k=0: t and bc_type can only "
-                                     "be notaknot.")
                 c = np.asarray(y_slice)
 
             # special-case k=1 (e.g., Lyche and Morken, Eq.(2.16))
             elif k == 1:
-                if not both_nak:
-                    raise ValueError("Too much info for k=1: bc_type can only be None.")
                 c = np.asarray(y_slice)
 
             else:
-                x_slice = _as_float_array(x_slice, check_finite)
                 y_slice = _as_float_array(y_slice, check_finite)
                 k = operator.index(k)
 
-
-                if x_slice.ndim != 1 or np.any(x_slice[1:] <= x_slice[:-1]):
-                    raise ValueError("Expect x_slice to be a 1-D sorted array_like.")
-                if k < 0:
-                    raise ValueError("Expect non-negative k.")
-                if t.ndim != 1 or np.any(t[1:] < t[:-1]):
-                    raise ValueError("Expect t to be a 1-D sorted array_like.")
                 if x_slice.size != y_slice.shape[0]:
                     raise ValueError('x_slice and y_slice are incompatible.')
-                if t.size < x_slice.size + k + 1:
-                    raise ValueError('Got %d knots, need at least %d.' %
-                                     (t.size, x_slice.size + k + 1))
-                if (x_slice[0] < t[k]) or (x_slice[-1] > t[-k]):
-                    raise ValueError('Out of bounds w/ x_slice = %s.' % x_slice)
 
                 # set up the RHS: values to interpolate (+ derivative values, if any)
                 extradim = prod(y_slice.shape[1:])
@@ -494,9 +464,5 @@ def make_interp_spline(x, y, bcs=0, orders=3):
 
                 c = np.ascontiguousarray(c.reshape((nt,) + y_slice.shape[1:]))
 
-            #=======================================
-            # START FROM SCIPY
-            #=======================================
-
             coefficients[coeff_line_sel] = c.T
     return BSplineNDInterpolator(knots, coefficients, orders,)

From 958d633611c5cb75553fa935c895b7a9b35c1073 Mon Sep 17 00:00:00 2001
From: ben <ben@sixpearls.com>
Date: Tue, 18 Jun 2019 08:50:01 -0700
Subject: [PATCH 09/10] dont over-write solution matrix for this case

---
 ndsplines/ndsplines.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/ndsplines/ndsplines.py b/ndsplines/ndsplines.py
index cfd6e1b..574e56c 100644
--- a/ndsplines/ndsplines.py
+++ b/ndsplines/ndsplines.py
@@ -455,7 +455,7 @@ def make_interp_spline(x, y, bcs=0, orders=3):
                     ab, rhs = map(np.asarray_chkfinite, (ab, rhs))
                 gbsv, = get_lapack_funcs(('gbsv',), (ab, rhs))
                 lu, piv, c, info = gbsv(kl, ku, ab, rhs,
-                        overwrite_ab=True, overwrite_b=True)
+                        overwrite_ab=False, overwrite_b=True)
 
                 if info > 0:
                     raise LinAlgError("Collocation matix is singular.")

From 522d4ffa82523b7b6991c86fd33d449513e10e42 Mon Sep 17 00:00:00 2001
From: ben <ben@sixpearls.com>
Date: Sat, 22 Jun 2019 20:22:23 -0700
Subject: [PATCH 10/10] single solve fo coefficient

---
 examples/1d-interp.py  |   2 +-
 examples/2d-interp.py  |  18 ++---
 ndsplines/ndsplines.py | 158 ++++++++++++++++++++---------------------
 3 files changed, 87 insertions(+), 91 deletions(-)

diff --git a/examples/1d-interp.py b/examples/1d-interp.py
index 8349264..5aa68f3 100644
--- a/examples/1d-interp.py
+++ b/examples/1d-interp.py
@@ -55,7 +55,7 @@ def tanh(x_in):
 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(4):
+for order in range(3,4):
     for func in funcs:
         fvals = func(x)
         truef = func(xx)
diff --git a/examples/2d-interp.py b/examples/2d-interp.py
index dc55d19..d3fcb01 100644
--- a/examples/2d-interp.py
+++ b/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()
diff --git a/ndsplines/ndsplines.py b/ndsplines/ndsplines.py
index 243c45b..35371f8 100644
--- a/ndsplines/ndsplines.py
+++ b/ndsplines/ndsplines.py
@@ -208,15 +208,16 @@ def make_lsq_spline(x, y, knots, orders, w=None, check_finite=True):
     return temp_spline
 
 def _not_a_knot(x, k, left=True, right=True):
-    if k > 2 and k % 2 == 0:
-        raise ValueError("Odd degree for now only. Got %s." % k)
     if k==0: # special case for k = 0, only apply to right
         left = False
-    t = x
+    if k % 2 == 0: # use greville sights if even
+        t = (x[1:] + x[:-1])/2.
+    else:
+        t = x
     if left:
-        t = np.r_[(t[0],)*(k+1), t[(k -1) //2 +1:]]
+        t = np.r_[(x[0],)*(k+1), t[(k -1) //2 +1:]]
     if right:
-        t = np.r_[t[:-(k-1)//2 -1 or None], (t[-1],)*(k+1)]
+        t = np.r_[t[:-((k-1)//2) -1 or None], (x[-1],)*(k+1)]
     return t
 
 def _augknt(x, k, left=True, right=True):
@@ -244,6 +245,8 @@ def make_interp_spline(x, y, bcs=0, orders=3):
         Use deriv_spec == 0 for not-a-knot boundary condition
         For k=0, both spec_values = 0 implements nearest-neighbor,
         a single side with spec_value = 0 uses zero-order-hold from that direction
+        For k>0, k even, can only apply not-a-knot BC to both sides or to neither.
+        In this case, not-a-knot uses Greville sites for knot locations.
     orders : ndarray, shape=(xdim,), dtype=np.intc
         Degree of interpolant for each axis (or broadcastable)
     periodic : ndarray, shape=(xdim,), dtype=np.bool_
@@ -307,6 +310,22 @@ def make_interp_spline(x, y, bcs=0, orders=3):
 
         x_slice = _as_float_array(x_slice, check_finite)
         # should there be a general check for k <= deriv_spec ?
+        y_deriv_spec_slice = slice(deriv_specs[i-1,0],-deriv_specs[i-1,1] or None)
+        swapped_view = np.swapaxes(coefficients, 0, i)
+        y_swapped_view = swapped_view[y_deriv_spec_slice, ...]
+        y_swapped_view_shape = y_swapped_view.shape
+        y_lines = _as_float_array(y_swapped_view.reshape(x_slice.shape + (-1,)), check_finite)
+
+        """
+        for idx in np.ndindex(*all_other_ax_shape):
+            offset_axes_remaining_sel = (tuple(idx[i-1:] + 
+                deriv_specs[i:,0]))
+            y_line_sel = ((Ellipsis,) + idx[:i-1] + 
+                (slice(deriv_specs[i-1,0],-deriv_specs[i-1,1] or None),) +
+                offset_axes_remaining_sel)
+            coeff_line_sel = ((Ellipsis,) + idx[:i-1] + (slice(None,None),)
+                + offset_axes_remaining_sel)
+        # """
 
         if k == 0:
             # all derivatives are fully defined, can only set 0-th derivative,
@@ -327,34 +346,32 @@ def make_interp_spline(x, y, bcs=0, orders=3):
             else:
                 raise ValueError("Set deriv_spec = 0, with up to one side = -1 for k=0")
 
-        # special-case k=1 (e.g., Lyche and Morken, Eq.(2.16))
+
+            knots.append(t)
+            # coefficients don't change
+            continue
+
         if k == 1:
             # all derivatives are fully defined, can only set 0-th derivative,
             # aka not-a-knot boundary conditions to both sides
             if not both_nak:
                 raise ValueError("Too much info for k=1: bc_type can only be notaknot.")
 
-        if k==2:
-            # it's possible this may be the best behavior for all even k > 0
-            if both_nak:
-                # special, ad-hoc case using greville sites
-                t = (x_slice[1:] + x_slice[:-1]) / 2.
-                t = np.r_[(x_slice[0],)*(k+1),
-                           t[1:-1],
-                           (x_slice[-1],)*(k+1)]
-
-            elif left_nak or right_nak:
-                raise ValueError("For k=2, can set both sides or neither side to notaknot")
-            else:
-                t = x_slice
+            t = _not_a_knot(x_slice, k)
+            knots.append(t)
+            # coefficients don't change
+            continue
 
-        elif k != 0:
-            t = _not_a_knot(x_slice, k, left_nak, right_nak)
+        if k%2 == 0:
+            if left_nak != right_nak:
+                raise ValueError("For k>0, k even, can set both sides or neither side to notaknot")
 
+        t = _not_a_knot(x_slice, k, left_nak, right_nak)
         t = _augknt(t, k, not left_nak, not right_nak)
 
         t = _as_float_array(t, check_finite)
-        
+        knots.append(t)
+
     
         if left_nak:
             deriv_l_ords = np.array([])
@@ -393,70 +410,47 @@ def make_interp_spline(x, y, bcs=0, orders=3):
             _sci_bspl._handle_lhs_derivatives(t, k, x_slice[-1], ab, kl, ku, deriv_r_ords,
                                     offset=nt-nright)
 
-        knots.append(t)
-
-        if k >= 2:
-            if x_slice.ndim != 1 or np.any(x_slice[1:] <= x_slice[:-1]):
-                raise ValueError("Expect x_slice to be a 1-D sorted array_like.")
-            if k < 0:
-                raise ValueError("Expect non-negative k.")
-            if t.ndim != 1 or np.any(t[1:] < t[:-1]):
-                raise ValueError("Expect t to be a 1-D sorted array_like.")
-            if t.size < x_slice.size + k + 1:
-                raise ValueError('Got %d knots, need at least %d.' %
-                                 (t.size, x_slice.size + k + 1))
-            if (x_slice[0] < t[k]) or (x_slice[-1] > t[-k]):
-                raise ValueError('Out of bounds w/ x_slice = %s.' % x_slice)
+        if x_slice.ndim != 1 or np.any(x_slice[1:] <= x_slice[:-1]):
+            raise ValueError("Expect x_slice to be a 1-D sorted array_like.")
+        if k < 0:
+            raise ValueError("Expect non-negative k.")
+        if t.ndim != 1 or np.any(t[1:] < t[:-1]):
+            raise ValueError("Expect t to be a 1-D sorted array_like.")
+        if t.size < x_slice.size + k + 1:
+            raise ValueError('Got %d knots, need at least %d.' %
+                             (t.size, x_slice.size + k + 1))
+        if (x_slice[0] < t[k]) or (x_slice[-1] > t[-k]):
+            raise ValueError('Out of bounds w/ x_slice = %s.' % x_slice)
+
+        k = operator.index(k)
+
+        if x_slice.size != y_lines.shape[0]:
+            raise ValueError('x_slice and y_lines are incompatible.')
+
+        # set up the RHS: values to interpolate (+ derivative values, if any)
+        extradim = y_lines.shape[1]
+        rhs = np.empty((nt, extradim), dtype=y_lines.dtype)
+        if nleft > 0:
+            rhs[:nleft] = deriv_l_vals.reshape(-1, extradim)
+        rhs[nleft:nt - nright] = y_lines.reshape(-1, extradim)
+        if nright > 0:
+            rhs[nt - nright:] = deriv_r_vals.reshape(-1, extradim)
 
+        # solve Ab @ x_slice = rhs; this is the relevant part of linalg.solve_banded
+        if check_finite:
+            ab, rhs = map(np.asarray_chkfinite, (ab, rhs))
+        gbsv, = get_lapack_funcs(('gbsv',), (ab, rhs))
 
-        for idx in np.ndindex(*all_other_ax_shape):
-            offset_axes_remaining_sel = (tuple(idx[i-1:] + 
-                deriv_specs[i:,0]))
-            y_line_sel = ((Ellipsis,) + idx[:i-1] + 
-                (slice(deriv_specs[i-1,0],-deriv_specs[i-1,1] or None),) +
-                offset_axes_remaining_sel)
-            coeff_line_sel = ((Ellipsis,) + idx[:i-1] + (slice(None,None),)
-                + offset_axes_remaining_sel)
+        lu, piv, c, info = gbsv(kl, ku, ab, rhs,
+                overwrite_ab=True, overwrite_b=True)
 
-            y_slice = coefficients[y_line_sel].T
+        if info > 0:
+            raise LinAlgError("Collocation matix is singular.")
+        elif info < 0:
+            raise ValueError('illegal value in %d-th argument of internal gbsv' % -info)
 
-            # special-case k=0 right away
-            if k == 0:
-                c = np.asarray(y_slice)
+        swapped_view[...] = c.reshape((-1,) +  y_swapped_view_shape[1:])
+            
 
-            # special-case k=1 (e.g., Lyche and Morken, Eq.(2.16))
-            elif k == 1:
-                c = np.asarray(y_slice)
 
-            else:
-                y_slice = _as_float_array(y_slice, check_finite)
-                k = operator.index(k)
-
-                if x_slice.size != y_slice.shape[0]:
-                    raise ValueError('x_slice and y_slice are incompatible.')
-
-                # set up the RHS: values to interpolate (+ derivative values, if any)
-                extradim = prod(y_slice.shape[1:])
-                rhs = np.empty((nt, extradim), dtype=y_slice.dtype)
-                if nleft > 0:
-                    rhs[:nleft] = deriv_l_vals.reshape(-1, extradim)
-                rhs[nleft:nt - nright] = y_slice.reshape(-1, extradim)
-                if nright > 0:
-                    rhs[nt - nright:] = deriv_r_vals.reshape(-1, extradim)
-
-                # solve Ab @ x_slice = rhs; this is the relevant part of linalg.solve_banded
-                if check_finite:
-                    ab, rhs = map(np.asarray_chkfinite, (ab, rhs))
-                gbsv, = get_lapack_funcs(('gbsv',), (ab, rhs))
-                lu, piv, c, info = gbsv(kl, ku, ab, rhs,
-                        overwrite_ab=False, overwrite_b=True)
-
-                if info > 0:
-                    raise LinAlgError("Collocation matix is singular.")
-                elif info < 0:
-                    raise ValueError('illegal value in %d-th argument of internal gbsv' % -info)
-
-                c = np.ascontiguousarray(c.reshape((nt,) + y_slice.shape[1:]))
-
-            coefficients[coeff_line_sel] = c.T
     return BSplineNDInterpolator(knots, coefficients, orders,)