2

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.

Improve this question
4
  • How do you define "clean pythonic way" and what's your metric for "a better way of doing this"? How can I compare two ways and measure which one is better? What do you mean with "legal"? If you want a code review of your correct and working code, you should post your question at Code Review SE Commented Jul 22 at 13:26
  • Apologies if my question is a bit vague. What I'm worried about is whether in some later version scipy or python, the above code snippet will stop working because I'm instantiating a class that is not meant to be instantiated directly. Also, is it standard practise to instantiate a class, and then manually overwrite its fields? Would it comply with PEP recommendations? Commented Jul 22 at 13:44
  • Also, one may wonder whether there's a function that does directly what I want to do, like the bispeu function, but which is not deprecated? I haven't found any in scipy, but perhaps I've missed something. Commented Jul 22 at 13:46
  • "whether in some later version scipy or python, the above code snippet will stop working" How could this anybody answer without a time machine? What is "standard practise"? Maybe also relevant here: Why is asking a question on "best practice" a bad thing? Commented Jul 22 at 13:52

1 Answer 1

Reset to default
2

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()
Improve this answer
Sign up to request clarification or add additional context in comments.

4 Comments

Thanks a lot for this answer. Yes, that's what I feared - the solution I propose is fragile and could break if implementation details change in later SciPy releases. I had seen the bisplev method, but it evaluates the spline over a Cartesian product of the x and y arrays. A workaround may involve looping through (x,y) pairs and applying bisplev to each pair individually - I will test this to see if it works reasonably fast. In answer to your question, the low-level fitpack routines are an order of magnitude faster than my naive python implementation of B-splines.
Ah, I didn't notice that about bisplev. One other note I found while researching this is that it is possible to use BivariateSpline subclasses without the output being gridded through the use of the grid=False parameter to BivariateSpline.__call__. See edit. This still doesn't address your concern about re-calculating coeffs - that will be done whenever you instantiate a new SmoothBivariateSpline.
I tried using bisplev within a loop but it turned out to be an order of magnitude slower. It may be an idea to use SmoothBivariateSpline given the grid=False parameter in the __call__ method, and it avoids instantiating BivariateSpline directly. One solution for the coefficients (and knot positions) may be to instantiate SmoothBivariateSpline with some minimal cost-efficient inputs then overwrite the relevant attributes with the appropriate coefficients and knot positions (like what I do with the BivariateSpline class above).
I think it would be more robust to stick with your original BivariateSpline idea than to selectively override SmoothBivariateSpline.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.