Feature/interpolation

[mkgessen]

This PR implements following features:

• Linear and Cubic 2D and 3D interpolators
  • with nearest neighbour and linear extrapolators
• Quadratic extrapolators in 1D
• Documentation at docs/build/html/api_reference/core/functions.html#interpolators

[jacklovell]

Thanks for this @mkgessen: I'm really keen to see this merged. It will help to fix https://github.com/cherab/core/issues/233 if we can get Cherab to use these interpolators rather than its existing ones. @CnlPepper is this going to make it into raysect 0.7.x, or will we have to wait for the 1.x refactor to take advantage of it?

[CnlPepper]

Thanks for this @mkgessen: I'm really keen to see this merged. It will help to fix cherab/core#233 if we can get Cherab to use these interpolators rather than its existing ones. @CnlPepper is this going to make it into raysect 0.7.x, or will we have to wait for the 1.x refactor to take advantage of it?

Hi Jack, we are developing this for 0.7.1. Once it is merged I'll make a Raysect release. Please be aware that these are not a direct drop in replacement for the cherab interpolators, there are some API differences. The code and API's are cleaner in the new interpolators, the unnecessary solve dependency is removed and they are faster, so all round should be worth the effort to move to these. I'd recommend deprecating the cherab interpolators for a release before they are removed.

[vsnever]

Hi All, I started testing these new interpolators in Cherab, and found undefined behaviour when interpolators are initialised with arrays of shape 1 in any dimension. Namely, this

f1d = Interpolator1DArray([1.], [2.], 'cubic', 'nearest', np.inf)
f2d = Interpolator2DArray([1.], [2.], [[3.]], 'cubic', 'nearest', np.inf, np.inf)
f3d = Interpolator2DArray([1.], [1., 2.], [[3., 4.]], 'cubic', 'nearest', np.inf, np.inf)

will pass all the checks in __init__, but may cause memory errors. Note that the checks like if (np.diff(x) <= 0).any() will pass fine because (np.array([]) <= 0).any() returns False. I will illustrate the problem on 1D case, but 2D and 3D interpolators have similar issues.

Calling f1d(1) leads to undefined behaviour because in this case find_index() returns 0 and the following branch is executed:

    elif px == self.x[self._last_index]:
        return self._interpolator.evaluate(px, index - 1)

So, the negative index goes to self._interpolator.evaluate(), which is decorated with @cython.wraparound(False), which means that negative indexing is not allowed.

Also, in _Interpolator1DCubic.__init__() the memoryviews _mask_a and _a will point to zero-shaped arrays, which is also a problem.

I think that trying to initialise these interpolators with arrays of shape 1 in any dimension should throw an error.

[Sorchard1]

Just an update on what I have changed since the last few comments:

I think that trying to initialise these interpolators with arrays of shape 1 in any dimension should throw an error.

Agreed this should be checked and I have added this check in init in 1D, 2D and 3D and updated the unit tests in 1D, 2D and 3D.

Looking a lot better. I did notice an issue however that needs resolving relating to an anti-pattern in a section of the code that needs to be optimised. I've expanded on the original comment.

I have refactored this, then noticed another similar case (which also made the calculation of derivatives in an array more confusing) so I have refactored in 2D and 3D where relevant.

[CnlPepper]

@Sorchard1 If you can address the last two queries (which probably apply to both 2D and 3D), I'll merge the code.