Add derivatives for the splines

[avik-pal]

TODOs:

• Implement derivatives for all the interpolation methods
• Integrate with ChainRulesCore
• Verify correctness using FiniteDifferences

Follow up from SciML/DiffEqFlux.jl#380

[avik-pal]

I am finding it hard to get the BSpline derivative right. From my understanding the derivative of a BSpline is another BSpline of lower order. Reference

N = spline_coefficients(n, A.d - 1, A.k[2:(end - 1)], t)
ducum = zero(eltype(A.u))
for i = 1:(n - 1)
      ducum += N[i] * (A.c[i + 1] - A.c[i]) / (A.k[i + A.d + 1] - A.k[i + 1])
end
ducum * A.d

But the gradients seem to be way off.

[ChrisRackauckas]

That looks like it might introduce some numerical errors? I think @YingboMa might have comments on this.

[avik-pal]

@YingboMa any thoughts on this?

[ChrisRackauckas]

Is there a version where BSplines are left off that we can open an issue on to complete later?

[YingboMa]

Sorry, I missed this PR. I will review it today.

[YingboMa]

I can help with the BSpline derivative later today.

[YingboMa]

For the BSpline derivative, you should check pg. 72 of the NURBS book https://www.springer.com/gp/book/9783642973857. I am not familiar with BSplines, but I think a bulletproof way of differentiating BSplines is just to differentiate each basis function individually.

[avik-pal]

The current implementation was differentiating each basis function itself. I removed the BSpline implementation for the time being so that we can move forward with the PR, and can tackle that later.

[YingboMa]

I think you have an off-by-one error in your code, also the derivative of the coordinate transformation is missing.

This gives the correct derivative

using DataInterpolations, Calculus
ts = 0:4
data = sin.(ts)
A = BSplineInterpolation(data,float.(ts),3,:ArcLen,:Average)
n = length(A.t)
t = 0.1
idx = searchsortedlast(A.t,t)
idx == length(A.t) ? idx -= 1 : nothing
t_ = A.p[idx] + (t - A.t[idx])/(A.t[idx+1] - A.t[idx]) * (A.p[idx+1] - A.p[idx])

Nm1 = DataInterpolations.spline_coefficients(n, A.d-1, A.k, t_)
scale = (A.t[idx+1] - A.t[idx]) * (A.p[idx+1] - A.p[idx])
ducum = zero(eltype(A.u))
for i = 1:(n - 1)
    ducum += Nm1[i+1] * (A.c[i + 1] - A.c[i]) / (A.k[i + A.d + 1] - A.k[i + 1])
end
ducum * A.d * scale
Calculus.derivative(A, t)

julia> ducum * A.d * scale
1.0854013436178458

julia> Calculus.derivative(A, t)
1.0854013436149472

[ChrisRackauckas]

I rebased this branch back to current master since it was falling behind a bit.