`gradient` does not support vector args which is inconsistent with interpolation evaluation

[danmackinlay]

perhaps this is an understanding problem, but it seems from the documentation that there is a disanalogy between of interpolant evaluation and the interpolant gradient evaluation which is not reflected in the documentation regarding vector arguments. Specifically, if I evaluate an interpolant at a vector of locations, it returns a vector of values, no problems. If I ask for the gradient at a vector of locations, it throws an error. Is this behaviour expected?

julia> itp = interpolate([1.,1.,2], (BSpline(Cubic(Line(OnGrid())))))
julia> itp( [1.0, 2.0],)
[1.0, 1.0]
julia> Interpolations.gradient(itp, [1.0, 2.0],)
gradient of [1.0, 1.0, 2.0] not supported for position ([1.0, 2.0],)

Stacktrace:
 [1] error(::String) at ./error.jl:33
 [2] gradient(::Interpolations.BSplineInterpolation{Float64,1,OffsetArrays.OffsetArray{Float64,1,Array{Float64,1}},BSpline{Cubic{Line{OnGrid}}},Tuple{Base.OneTo{Int64}}}, ::Array{Float64,1}) at /home/me/.julia/packages/Interpolations/AiRSM/src/Interpolations.jl:369
 [3] top-level scope at In[84]:4

Of course I can call the gradient function multiple times to get multiple cvalues, but this seems to my limited understanding to construct the gradient weights and such multiple times, and so likely to be inefficient, as well as unnecessarily different to the normal interpolation evaluation. Am I missing something?

[danmackinlay]

I suppose there is also broadcasting; is that the recommended method?

julia> Interpolations.gradient.([itp], [1.0, 2.0],)
StaticArrays.SArray{Tuple{1},Float64,1,1}[[-0.25], [0.5]]

[timholy]

Broadcasting is the default option. I've contemplated removing it from evaluation of the value, but OTOH there are sometimes efficiencies to be gained by evaluating at many points, if they happen to be ordered: https://discourse.julialang.org/t/using-interpolations/23933/31?u=tim.holy

I therefore wouldn't be opposed to a PR that supported it. (I'm not planning on doing it anytime soon myself.) In the meantime, you can use Ref(itp) which is slightly more efficient than [itp].

[danmackinlay]

Well said @timholy ! And thanks. What I will do is my sadly classic MO: get my project protoype going with a hacked Fourier-based interpolation method, and return to it if time allows after the paper deadline. I give it odds 4:1 against but let's keep this ticket open for the happy chance I (or someone else) actually get allocated time to put something back in to this (incredibly useful) project.

If I actually get to hack on this, I will try to document a little too; it's very opaque to an outsider, the project internals at the moment, and whiles the code is very tidy it's not exactly clear what's happening. There is a bootstrap problem here ;)

[henry2004y]

I totally agree with @danmackinlay. At least there should an example of the gradient of vector in the document. Broadcast is the assumed usage, but I sometimes totally forget about it...