Fix type stability issue with linear interpolation #113

[GlenHertz]

The fast path for linear interpolation when no interpolation is needed (time point is provided in vector) returned a different type which caused issues with types/performance with higher level code.

I ran some benchmarks and see no difference in performance.

[GlenHertz]

I'm not sure if I'm going about this the right way but it works.

[GlenHertz]

Note this changes how the derivative works. With the fast-path change the derivative of the LinearInterpolation at the exact time point is zero (from ForwardDiff) while previously it was the derivative of the current point and the next.

[andreasnoack]

Note this changes how the derivative works.

Good point. Hopefully people don't rely too much on the derivative here. If they did then they should probably consider a different interpolation. However, you could also consider restricting the branches to when isnan(ret). That was your original motivation anyway and will leave the derivative the way it used to be.

[GlenHertz]

The derivative at the exact point is undefined but I think it is better to stay consistent with before and if they depend on the derivative they should probably use a better interpolation scheme. I've seen other tools use the average slope from either side for this but that seems like a different/hack interpolation method to me.

[ChrisRackauckas]

Good point. Hopefully people don't rely too much on the derivative here.

They do more than they expect. For example, the ODE solver will differentiate at those points if using a stiff ODE solver, and it will step there if it's in a tstop. So this is actually pretty important. But, we already have frule and rrule overloads, it's just for ForwardDiff. This fixes potential ForwardDiff issues too, so it's a pretty good change.

[ChrisRackauckas]

I've seen other tools use the average slope from either side for this but that seems like a different/hack interpolation method to me.

The derivative on the right would be preferred due to the ODE solver context.

[ChrisRackauckas]

@GlenHertz can you do a patch bump?

[GlenHertz]

@ChrisRackauckas I'm not sure what a patch bump is but I added another commit for the derivative fix.

I'm trying to get my code to work with NaNs and I ran into this odd case:

t = 0:3; 
u = [0, -2, -1, -2]; 
A = LinearInterpolation(u, t); 
dA = t -> ForwardDiff.derivative(A, t)
dA(NaN) # returns -1.0

It is because if NaN is given then the idx will be the second last one and it will load real points from the vectors and then it has a derivative. I'll work on a patch.

[GlenHertz]

The tests fail right now because using a NaN input for t with Integer vector doesn't compute....I'll get back to this tomorrow.

[GlenHertz]

I have another update and this one looks good to me but please take a critical look. Originally, for input t=NaN searchsortedlast would find the last bracket which isn't correct and thus the derivative of NaN would be the derivative of the last bracket.
My fix looks a bit strange but it sets the inputs to the same values so the interpolation calculation is bad and therefore the derivative is bad too.

[ChrisRackauckas]

Test that the derivative at a point gives the derivative to the right?