scipy/scipy — pull request #3260 by bmcage (2014-01-31 13:55:35+00:00)
Title: This adds root finding to pchip as this is possible to do optimal because pchip is a monotonic interpolation.

Alternative would be fsolve and friends, but those would not use the known structure of pchip.

A typical use for pchip is an interpolation that can be inverted. This adds a root finding to pchip that uses the pchip form to invert fast.

I am not sure this is the best way. It might be better to use phcip to obtain yi=f(xi) and yi' = f'(xi), then create the cubic monotonic interpolation through yi with derivatives yi' as self.pchipinv. This should be the correct f^-1.
Then when asking roots, one can use this self.pchipinv.

The algorithm proposed here works differently: search the cubic for the y value given, calculate the 3 cubic roots. As it is a monotonic function only 1 of the roots will be in the correct interval.

Some problems:
1. computing cubic roots gives imaginary parts which need to be discarded.
See the np.allclose(ires, 0., atol=1e-10) to discard this. -10 is hard coded
2. return value can be numerically close to the bounding edges.
See the np.allclose(x1-rres, 0., atol=1e-10) to recognize this. -10 is hard coded
3. Use of newaxis to dump the results: x[c] = ress[:, np.newaxis]
I believe this will be correct for the use cases possible