Pulling x-values of spline integer y-values?

Question

I have the given sample data and interpolated spline:

import matplotlib.pyplot as plt 
import numpy as np
from scipy import interpolate

x = [0.1, 1.5, 2.3, 5.5, 6.7, 7]
y = [5, 4, 5, 2, 2, 3]
s = interpolate.UnivariateSpline(x, y, s=0)

xs = np.linspace(min(x), max(x), 10000) #values for x axis
ys = s(xs)

plt.figure()
plt.plot(xs, ys, label='spline')
plt.plot(x, y, 'x', label='collected data')
plt.legend()

I would like to pull the x values that correspond to the integer y values of the spline, but am not sure how to do this. I assume I will be using np.where() and have tried (to no avail):

root_index = np.where(ys == ys.round())

Does this answer your question? [Invert interpolation to give the variable associated with a desired interpolation function value](https://stackoverflow.com/questions/47275957/invert-interpolation-to-give-the-variable-associated-with-a-desired-interpolatio) — Sayandip Dutta, Feb 26 '21 at 20:09
Rather than generating the dense grid in the accepted answer, I would recommend the more analytical solution of calling `solve` on a `CubicSpline` object. — Mad Physicist, Feb 26 '21 at 21:23

Quang Hoang · Answer 1 · 2021-02-26T20:07:11.120

2

It's not a good idea comparing two floats with equal. Use:

# play with `thresh`
thresh = 1e-4

root_index = np.where(np.abs(ys - ys.round())<1e-4)

printt(xs[root_index])
print(ys[root_index])

Output:

array([0.1       , 0.44779478, 2.29993999, 4.83732373, 7.        ])
array([5.        , 4.00009501, 4.99994831, 3.00006626, 3.        ])

Numpy also has np.isclose. So something like:

root_index = np.isclose(ys, ys.round(), rtol=1e-5, atol=1e-4)
ys[root_index]

And you have the same output.

edited Feb 26 '21 at 20:07

answered Feb 26 '21 at 19:54

Quang Hoang

146,074
10
56
74

Hong, thank you for your solution, but why does this not retrieve the x-values for when the y-values are 2 and for all of the times y = 4 or 5? – Andrea Feb 26 '21 at 20:04
For the `x` values, use `xs[root_index]`. See updated answer. – Quang Hoang Feb 26 '21 at 20:06

JohanC · Accepted Answer · 2021-02-26T20:41:00.410

You could use the find_roots function from this post to find the exact interpolated x-values:

import matplotlib.pyplot as plt
import numpy as np
from scipy import interpolate

def find_roots(x, y):
    s = np.abs(np.diff(np.sign(y))).astype(bool)
    return x[:-1][s] + np.diff(x)[s] / (np.abs(y[1:][s] / y[:-1][s]) + 1)

x = [0.1, 1.5, 2.3, 5.5, 6.7, 7]
y = [5, 4, 5, 2, 2, 3]
s = interpolate.UnivariateSpline(x, y, s=0)

xs = np.linspace(min(x), max(x), 500)  # values for x axis
ys = s(xs)

plt.figure()
plt.plot(xs, ys, label='spline')

for y0 in range(0, 7):
    r = find_roots(xs, ys - y0)
    if len(r) > 0:
        plt.scatter(r, np.repeat(y0, len(r)), label=f'y0 = {y0}')
        for ri in r:
            plt.text(ri, y0, f'  {ri:.3f}', ha='left', va='center')
plt.legend()
plt.xlim(min(x) - 0.2, max(x) + 1)
plt.show()

PS: with much less points for the x's, e.g. xs = np.linspace(min(x), max(x), 50), the curve would look a bit bumpy, and the interpolation would be slightly different:

This looks like a linear interpolation of a sufficiently finely sampled grid. — Mad Physicist, Feb 26 '21 at 20:34

Mad Physicist · Answer 3 · 2021-02-27T15:29:27.523

The object returned by scipy.interpolate.UnivariateSpline is a wrapper about fitpack interpolation routines. You can get an identical interpolation using the CubicSpline class: replace s = interpolation.UnivariateSpline(x, y, s=0) with s = interpolation.CubicSpline(x, y). The results are identical, as you can see in the figure at the end.

The advantage to using CubicSpline is that it returns a PPoly object which has working roots and solve methods. You can use this to compute the roots with integer offsets.

To compute the range of possible integers, you can use the derivative method coupled with roots:

x_extrema = np.concatenate((x[:1], s.derivative().roots(), x[-1:]))
y_extrema = s(x_extrema)
i_min = np.ceil(y_extrema.min())
i_max = np.floor(y_extrema.max())

Now you can compute the roots of the offset spline:

x_ints = [s.solve(i) for i in np.arange(i_min, i_max + 1)]
x_ints = np.concatenate(x_ints)
x_ints.sort()

Here are some additional plotting commands and their output:

plt.figure()
plt.plot(xs, interpolate.UnivariateSpline(x, y, s=0)(xs), label='Original Spline')
plt.plot(xs, ys, 'r:', label='CubicSpline')
plt.plot(x, y, 'x', label = 'Collected Data')

plt.plot(x_extrema, y_extrema, 's', label='Extrema')
plt.hlines([i_min, i_max], xs[0], xs[-1], label='Integer Limits')
plt.plot(x_ints, s(x_ints), '^', label='Integer Y')

plt.legend()

You can verify the numbers:

>>> s(x_ints)
array([5., 4., 4., 5., 5., 4., 3., 2., 2., 3., 4., 5.])

Pulling x-values of spline integer y-values?

3 Answers3