I need to efficiently evaluate a bivariate spline on a B-spline basis. I have already calculated the knot positions and spline coefficients (independently of scipy classes/methods such as RectBivariateSpline which calculate these for you). Hence I only want to evaluate the spline function itself using the knots and coefficients I provide. This spline function needs to be evaluated on a cloud of points rather than a grid formed of a Cartesian product. I came across the following stackoverflow query which explains how to do it: Querying multiple points on a bivariate spline in the B-spline basis . For convenience, I reproduce the code which does what I want:
from scipy.interpolate import dfitpack
# Meaning of the different variables:
# x: x coordinates where to evaluate the bivariate spline function
# y: y coordinates where to evaluate the bivariate spline function
# xknots: the positions of the knots in the x direction
# yknots: the positions of the knots in the y direction
# coefs: spline coefficients in the bivariate spline function
# result: the bivariate spline function evaluated at x, y
#
# NOTE: I'm assuming x, y, xknots, yknots, and coefs have already been
# calculated. In the present example, x and y may be multi-
# dimentional arrays, but need to be transformed into a 1D array.
# coefs is a 2D array given that we are dealing with a bivariate
# spline function.
x = x.ravel(order="C")
y = y.ravel(order="C")
# NOTE: specify the order parameter to avoid leaving this to chance
# NOTE: access low-level fortran routine to carry out calculations (this
# runs much faster than any python equivalent, but the documentation
# is only available in the fortran routine or online)
result, ier = dfitpack.bispeu(
yknots, xknots, coefs.ravel(order="C"), degree, degree, y, x
)
if not ier == 0:
raise ValueError("Error in spline2grid. Error code from bispeu: %s" % ier)
However, accessing the low level Fortran routine bispeu directly no longer works in later versions of scipy. I get DeprecationWarning in scipy-1.14.1, and the program fails in later versions. Therefore I have resorted to doing this instead:
from scipy.interpolate import BivariateSpline
# This example uses the same naming convention and assumptions as
# in the previous code snippet.
# set up BivariateSpline object (to access bispeu indirectly)
spline_object = BivariateSpline() # is this legal ?
spline_object.degrees = (degree, degree)
spline_object.tck = (yknots, xknots, coefs.ravel(order="C"))
x = x.ravel(order="C")
y = y.ravel(order="C")
result = spline_object(y, x, grid=False)
The above solution works, but I'm wondering if it's legal or if there is a better way of doing this. For instance, the documentation to BivariateSpline state "This class is meant to be subclassed, not instantiated directly." However, if I try using a subclass such as RectBivariateSpline, it typically wants to recalculate the knots and coefficients (which is a waste of time and doesn't correspond to what I want), and may also try to force me to calculate the spline on a Cartesian grid rather than a cloud of points.
This spline function needs to be evaluated on a cloud of points rather than a grid formed of a Cartesian product.
You are looking for the "2D FITPACK splines / unstructured data" section of the interpolate docs.
Specifically, bisplrep and bisplev seem closest to your requirements. They differ from your requirements in only that they are based on SURFIT rather than BISPEU. (Edit: and that bisplev emits a gridded output)
NOTE: access low-level fortran routine to carry out calculations (this runs much faster than any python equivalent, but the documentation is only available in the fortran routine or online) result, ier = dfitpack.bispeu(
The question I would ask is if using low-level fitpack routines is really faster.
For example, BivariateSpline is implemented using bispeu. In fact, a lot of SciPy's high-level functions are just wrappers around Fortran/C++ code. Often it has additional error checking (e.g. did you pass an array? is the array an appropriate type? does the array have the correct memory ordering?) but I would be surprised if the overhead of those checks were significant enough to worry about.
However, accessing the low level Fortran routine bispeu directly no longer works in later versions of scipy. I get DeprecationWarning in scipy-1.14.1, and the program fails in later versions.
Right, that's part of SciPy's private API, so I would not be surprised if it breaks between versions.
The above solution works, but I'm wondering if it's legal or if there is a better way of doing this.
In my opinion, that would fall under part of SciPy's private API. It works, but it requires using undocumented details of the implementation that could change in a future release.
Edit: Example of how to use a BivariateSpline subclass without gridded input or output:
import scipy
import numpy as np
import matplotlib.pyplot as plt
x = np.random.uniform(-1, 1, size=100)
y = np.random.uniform(-1, 1, size=100)
def gaussfunc(x, y):
z = np.exp(-np.sqrt(x ** 2 + y ** 2))
return z
z = gaussfunc(x, y)
# plt.scatter(x, y, c=z)
# plt.colorbar()
evaluation_x = np.random.uniform(-1, 1, size=100)
evaluation_y = np.random.uniform(-1, 1, size=100)
# Fit spline
spline = scipy.interpolate.SmoothBivariateSpline(x, y, z, s=len(z) * 1e-3)
# Evaluate spline
evaluation_z = spline(evaluation_x, evaluation_y, grid=False)
plt.scatter(evaluation_x, evaluation_y, c=evaluation_z)
plt.colorbar()