ENH: implement Akima interpolation in 1D #3112

Closed
wants to merge 2 commits into
from

Conversation

Projects
None yet
5 participants
Contributor

andreas-h commented Dec 3, 2013

As promised in #2885, here's my PR for Akima interpolation in 1D. While there are several improvements thinkable, I believe it's already useful as-is and would appreciate inclusion in the upcoming 0.14 release.

Owner

pv commented Dec 3, 2013

fill_value doesn't do anything. Would be better to remove both it and the periodic keywords.

Otherwise, looks OK.

Coverage Status

Changes Unknown when pulling 06f435c9a5a4dfd77463a7b7f81f34036e2c5fce on andreas-h:Akima1D into * on scipy:master*.

Coverage Status

Coverage remained the same when pulling dfcbcb9469f816bffc27e92fb86e46358b57d945 on andreas-h:Akima1D into f86894e on scipy:master.

Coverage Status

Coverage remained the same when pulling dfcbcb9469f816bffc27e92fb86e46358b57d945 on andreas-h:Akima1D into f86894e on scipy:master.

@ev-br ev-br commented on an outdated diff Dec 3, 2013

scipy/interpolate/interpolate.py
+ # These are the indices where the the slope at breakpoint is defined:
+ id_ = np.nonzero(f12 > 1e-9 * np.max(f12))[0]
+ # set the slope at breakpoint
+ t[id_] = (f1[id_] * m[id_ + 1] + f2[id_] * m[id_ + 2]) / f12[id_]
+
+ # calculate the higher order coefficients
+ c = (3. * m[2:-2] - 2. * t[:-1] - t[1:]) / dx
+ d = (t[:-1] + t[1:] - 2. * m[2:-2]) / dx ** 2
+
+ coeff = np.zeros((4, x.size - 1))
+ coeff[3] = y[:-1]
+ coeff[2] = t[:-1]
+ coeff[1] = c
+ coeff[0] = d
+
+ PPoly.__init__(self, coeff, x, extrapolate=False)
@ev-br

ev-br Dec 3, 2013

Member

a super call would play better with subclassing

@ev-br

ev-br Dec 3, 2013

Member

oh, mea culpa
I think you want super(Akima1D, self).__(coeff, x)

Member

ev-br commented Dec 3, 2013

Periodic boundary conditions: https://gist.github.com/ev-br/7110020

Contributor

andreas-h commented Dec 3, 2013

Thanks for the pointer to your code. If you want me to, I can prepare a seperate PR to include periodic boundary conditions for PPoly.

Contributor

andreas-h commented Dec 3, 2013

Updated subclassing to use super() call.

Member

ev-br commented Dec 3, 2013

Feel free to do whatever you feel like with these lines of code.
We had a discussion with Pauli about what to do with out-of-bounds behavior in his PPoly PR. I think the decision was to leave them to a user.
Given that periodic boundary conditions are literally five lines of code, I'm not sure it's worth inclusion into scipy. And if we do have PBCs, why not also antiperiodic (DMFT people use antiperiodic Akimas), mirror-symmetric (b-splines in signal) etc.

An optional enhancement would be to clarify the relation of the Akima splines and pchip.

Irrespective of this, I'm +1 on this PR.

Owner

pv commented Dec 3, 2013

One thing needs to be added: it should support multidimensional y-arrays, not only 1D. This feature is very often needed.

Coverage Status

Coverage remained the same when pulling f91f7c2 on andreas-h:Akima1D into 10d88c4 on scipy:master.

Contributor

andreas-h commented Dec 9, 2013

I agree multidimensional y arrays would be useful. The necessary machinery is already there in the _Interpolator1D class. However, when I try to inherit from both PPoly and _Interpolator1D, I get into some inheritance problems I've never encountered before:

TypeError: Error when calling the metaclass bases
multiple bases have instance lay-out conflict

I'm not sure where this comes from, but it might be because the __slots__ are defined in both parent classes. Any ideas how I could go about this?

As an alternative approach to this, I could go a different route: implement the Akima algorithm in a helper class _Akima1D and then extent the interp1d class to support kind=akima. What do you think?

Owner

pv commented Dec 9, 2013

It would probably be enough to support n-d so that the interpolation axis is the first one. The PPoly/BPoly can handle this, so in Akima the only thing to do is to make sure the coefficient computation is vectorized accordingly.

Owner

pv commented Dec 13, 2013

The current code in __init__ almost handles the n-d case, too. It seems the only thing needed is to change the coeff array shape to (4,) + y.shape, and to add some n-d tests.

If the axis transposition is a feature that could be useful, it should be added in PPoly. There's no need to inherit from _Interpolator1D here, as it's only a helper class whose main purpose is to make axis transpositions etc. easier.

Contributor

andreas-h commented Feb 8, 2014

Done; Akima now doesn't inherit _Interpolator1D any more, and supports n-d y.

Is there a more elegant way to do the "broadcasting" I'm doing in l. 1433-1435? It feels hackish what I'm doing there.

Contributor

andreas-h commented Feb 8, 2014

Sorry, just realized I amended my changes to the last commit and did a --force push. Seems like I'm still doing stupid things with git ... :-/

@pv pv commented on an outdated diff Feb 8, 2014

scipy/interpolate/interpolate.py
+ def __init__(self, x, y):
+ # Original implementation in MATLAB by N. Shamsundar (BSD licensed), see
+ # http://www.mathworks.de/matlabcentral/fileexchange/1814-akima-interpolation
+ if np.any(np.diff(x) < 0.):
+ raise ValueError("x must be strictly ascending")
+ if x.ndim != 1:
+ raise ValueError("x must be 1-dimensional")
+ if x.size < 2:
+ raise ValueError("at least 2 breakpoints are needed")
+ if x.size != y.shape[0]:
+ raise ValueError("x.shape must equal y.shape[0]")
+
+ # determine slopes between breakpoints
+ m = np.empty((x.size + 3, ) + y.shape[1:])
+ dx = np.diff(x)
+ for i in range(y.ndim - 1):
@pv

pv Feb 8, 2014

Owner

Equivalent to dx = dx[(slice(None),) + (None,)*(y.ndim-1)] ?
If so, better not use as_strided

Owner

pv commented Feb 8, 2014

Force-pushing to pull requests is also OK

Contributor

andreas-h commented Feb 18, 2014

If this is okay, can we get it into 0.14?

Owner

pv commented Feb 18, 2014

I'll (i) rebase, (ii) rename to AkimaInterpolator, and (iii) move the implementation to _monotone.py together with pchip

Member

ev-br commented Feb 18, 2014

Could you also add a 'see also' crosslinks to both docstrings.

Owner

pv commented Feb 18, 2014

Ok, actually this has to wait for tomorrow from my side, got to go...
But looks ready to go to me.

@ev-br ev-br and 1 other commented on an outdated diff Feb 18, 2014

scipy/interpolate/interpolate.py
+ """
+ Akima interpolator
+
+ Fit piecewise cubic polynomials, given vectors x and y. The interpolation
+ method by Akima uses a continuously differentiable sub-spline built from
+ piecewise cubic polynomials. The resultant curve passes through the given
+ data points and will appear smooth and natural.
+
+ Parameters
+ ----------
+ x : ndarray, shape (m, )
+ 1-D array of monotonically increasing real values.
+ y : ndarray, shape (m, ...)
+ N-D array of real values. The length of *y* along the first axis must
+ be equal to the length of *x*.
+ axis : int, optional
@ev-br

ev-br Feb 18, 2014

Member

left over? The init below does not seem to have an axis argument (I realize I'm a little late to the party...)

@andreas-h

andreas-h Feb 18, 2014

Contributor

True ... will force-push fix

Member

ev-br commented Feb 18, 2014

(this is by no means intended to block merging)

Can you comment on this: https://gist.github.com/ev-br/9081651

I've to admit I did not read Akima's paper, but I somehow had an impression that akima splines were supposed to not overshoot?

Contributor

andreas-h commented Feb 18, 2014

interesting. I'll need to check, hopefully tomorrow. I know Akima's are supposed not to overshoot, but I don't remember if they're guaranteed not to do so. Have to go now, but will come back to this very soon.

Member

ev-br commented Feb 19, 2014

Handling of arbitrary axis can be done in a way similar to pchip: basically a single private helper is enough to make the behavior consistent with that of polyint-based interpolators. Which can be done now or left as an optional enhancement.

pv added the PR label Feb 19, 2014

Owner

pv commented Feb 19, 2014

I think the implementation here is correct. I get the same overshooting curve from other Akima implementations, so it's a property of the algorithm.

Owner

pv commented Feb 19, 2014

Thanks a lot Andreas! PR merged in 358d1aa.

pv closed this Feb 19, 2014

pv added this to the 0.14.0 milestone Feb 19, 2014

andreas-h deleted the andreas-h:Akima1D branch Feb 24, 2014

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment