Get x-values corresponding to y-value on non-monotonic curves

Question

First, I'm not sure SO is the good place to ask this kind of question.
If not please tell me where to ask and I'll remove it.

However, since I not found any post that address my problem, my question and my code below may be useful to other people ...

The code part is a bit long, although I tried to keep only the minimum. However, this part should only serves to show that I have done some research and that I am looking for something better.

The question:

From a x-values and y-values lists, I want to find (or guess) the x-values which corresponds to a given y-value.

When the curve defined by X and Y values is monotonic, I can simply interpolate the function x = f(y) :

import numpy as np
from scipy import interpolate

class AxisCam:
    def __init__(self, x=None, y=None):
        self.x = x if x else []
        self.y = y if y else []

        if len(self.x):
            self.xMin = min(self.x)
            self.xMax = max(self.x)
        else:
            self.xMin = None
            self.xMax = None

        if len(self.y):
            self.yMin = min(self.y)
            self.yMax = max(self.y)
        else:
            self.yMin = None
            self.yMax = None

        self._interpolX, self._interpolY = self.setInterpolator()

    def setInterpolator(self, interpolator=interpolate.interp1d):
        """
        Define the interpolator to use to approximate the axis cam positions
        :param interpolator: interpolator function to use, default is scipy.interpolate.interp1d
        :return: a tuple with the interpolator functions for x and y values
        """
        if len(self.x) <= 0 or len(self.y) <= 0:
            return None, None
        with np.errstate(divide='ignore', invalid='ignore'):  # silent the warnings caused by the interpolator
            self._interpolX = interpolator(self.y, self.x)  # x = f(y)
            self._interpolY = interpolator(self.x, self.y)  # y = f(x)
        return self._interpolX, self._interpolY

    def getX(self, yValue):
        """
        Return x-value corresponding to a y-value using the interpolator
        :param yValue: y-value we want to know the corresponding x-value
        :return: x-value corresponding to the given y-value
        """
        if yValue < self.yMin:
            raise ValueError("value should be greater than the minimum y-value")
        elif yValue > self.yMax:
            raise ValueError("value should be lesser than the maximum y-value")
        return float(self._interpolX(yValue))

    def getY(self, value):
        """
        Return a y-value corresponding to a x-value using the interpolator
        :param value: x-value we want to know the corresponding y-value
        :return: the y-value corresponding to the given x-value
        """
        if value < self.xMin:
            raise ValueError("value should be greater than the minimum x-value")
        elif value > self.xMax:
            raise ValueError("value should be lesser than the maximum x-value")
        return float(self._interpolY(value))

x = [0, 0.351906, 0.703812, 1.055718]  # The 1024 values for X and Y can be retrieved here : https://pastebin.com/5eHsRjZ3
y = [0.0, 0.000306, 0.002419, 0.008111]
ac = AxisCam(x, y)
print(ac.getX(100))  # returns 30.124163768271398

However, when the curve is non-monotonic, I can't. An exception is raised

ValueError: x must be strictly increasing

So, for now, I splitted the curve into monotonics part using getMonotonicParts method below and I can interpolate the function x = f(y) on each monotonic parts.

import numpy as np
from scipy import interpolate

class AxisCam:
    def __init__(self, x=None, y=None):
        self.x = x if x else []
        self.y = y if y else []

        if len(self.y):
            self.yMin = min(self.y)
            self.yMax = max(self.y)
        else:
            self.yMin = None
            self.yMax = None

        self._monotonicParts = self.getMonotonicParts()

    def getMonotonicParts(self, interpolator=interpolate.interp1d):
        parts = []
        prevY = None  # will store the previous value of y to compare with at each iteration
        startIdx = None  # will store the index of self.x and self.y where the monotonic part start from
        direction = 0  # 0: Unknown - 1 : constant - 2: ascending - 3: descending
        lenY = len(self.y)
        for i, (x, y) in enumerate(zip(self.x, self.y)):
            if prevY is None:
                prevY = y
            if startIdx is None:
                startIdx = i

            prevDir = direction
            direction = 1 if y == prevY else 2 if y > prevY else 3
            if prevDir != 0 and prevDir != direction:  # Direction has changed => we have a new monotonic part
                endIdx = i - 1
                if direction == 3:  # y values are increasing => we can interpolate on it
                    interp_func = interpolator(self.y[startIdx:endIdx], self.x[startIdx:endIdx])
                elif direction == 1:  # y values are decreasing => we need to reverse it to interpolate on it
                    xValues = self.x[startIdx:endIdx]
                    xValues.reverse()
                    yValues = self.y[startIdx:endIdx]
                    yValues.reverse()
                    interp_func = interpolator(yValues, xValues)
                else:  # y values are the same on the range => return one of these
                    def interp_func(value): return self.y[startIdx]
                parts.append({'start': startIdx,
                              'end': endIdx,
                              'x0': self.x[startIdx],
                              'y0': self.y[startIdx],
                              'x1': self.x[endIdx],
                              'y1': self.y[endIdx],
                              'interp': interp_func})
                startIdx = i
            elif i == lenY - 1:  # Add element on the last iteration
                endIdx = i
                if direction == 2:
                    interp = interpolator(self.y[startIdx:endIdx], self.x[startIdx:endIdx])
                else:
                    interp = None
                parts.append({'start': startIdx,
                              'end': endIdx,
                              'x0': self.x[startIdx],
                              'y0': self.y[startIdx],
                              'x1': self.x[endIdx],
                              'y1': self.y[endIdx],
                              'interp': interp})
            prevY = y
        return parts

    def getX(self, yValue):
        """
        Return a list of x-values corresponding to a y-value using the interpolator
        :param yValue: y-value we want to know the corresponding x-value
        :return: a list of x-values corresponding to the given y-value
        """
        if yValue < self.yMin:
            raise ValueError("value should be greater than the minimum y-value")
        elif yValue > self.yMax:
            raise ValueError("value should be lesser than the maximum y-value")
        xValues = []
        for part in self._monotonicParts:
            if part['y0'] <= yValue <= part['y1'] or part['y0'] >= yValue >= part['y1']:
                xValues.append(float(part['interp'](yValue)))
        return xValues

x = []  # The 1024 values for X and Y can be retrieved here : https://pastebin.com/SL9RYYxY
y = []  # /!\ It is not the same values that the previous example /!\
ac = AxisCam(x, y)
print(ac.getX(100))  # returns [122.96996037206237, 207.6239552142487]

My solution works quite well but seems a bit far-fetched to me and I'm wondering if there is another better methods to do this.

Answer 1

I don't know of any standard routine to do this kind of multi-interpolation. But if you want this to scale, you should refactor your code a bit to use everything that numpy offers. For example, you could do something like this:

import numpy as np
from scipy.interpolate import interp1d

# convert data lists to arrays
x, y = np.array(x), np.array(y)

# sort x and y by x value
order = np.argsort(x)
xsort, ysort = x[order], y[order]

# compute indices of points where y changes direction
ydirection = np.sign(np.diff(ysort))
changepoints = 1 + np.where(np.diff(ydirection) != 0)[0]

# find groups of x and y within which y is monotonic
xgroups = np.split(xsort, changepoints)
ygroups = np.split(ysort, changepoints)
interps = [interp1d(y, x, bounds_error=False) for y, x in zip(ygroups, xgroups)]

# interpolate all y values
yval = 100
xvals = np.array([interp(yval) for interp in interps])

print(xvals)
# array([          nan,  122.96996037,  207.62395521,           nan])

Here nan indicates a value that is out of range (this algorithm treats repeated values as a separate group).