Calling with function syntax?

[jlperla]

I think your perspective changes somewhat depending on what you are using it for, but I strongly dislike having to call the interpolated function with array indexing as opposed to functions. For example, given the following interpolation of a function?

using Interpolations
x = 0:0.1:1
y = x.^2;
itp = interpolate((x,), y, Gridded(Linear()))
itp[0.52]

Now add in...

function (itp::Interpolations.GriddedInterpolation{T,N,TCoefs,IT,K,pad})(args...)where {T,N,TCoefs,IT,K,pad}
    itp[args...]
end
itp(0.52) # works now!
itp.([0.52 0.53]) #Broadcast!

Is there any reason not to add in this function (at least for gridded?)

[tomasaschan]

@jlperla IIUC, what you're suggesting is to implement itp(0.52) and itp.([0.52 0.53]) as a complement to itp[0.52] and itp[[0.52 0.53]] - is that correctly understood? Is there a reason for wanting this change besides liking one syntax over the other?

[jlperla]

It is very practical, not just philosophical. But to give some perspective on how I think about this.

In Julia, functions that map integers, i.e. f : N -> R are usually called by the square brackets, f[2] etc. When you interpolate those sorts of functions, you are filling in some of the integers. Interpolation spits out something you can use as f[3] with any integer holes filled in. Now, you can use any higher-order functions or write code for those. There is certainly precedent in math for this, for example difference equations are often written f[t] or f_t in math.

In Julia, functions that map general parameters to parameters (e.g. f : R -> R) are usually called as f(1.1). If I generate an interpolation, then I have a new function which works across a larger set of R than I defined it at. But conceptually, it is still a function and users expect to call it as such, which means they expect ( )

What this means is that if I am writing generic code on functions, unless the interpolation supports f(2.0) notation, it won't work. For example, if I want to use quadrature on the interpolated function with quadgk then I have to go quadgk( x-> f[x], 0.0, 1.0) which doesn't make much sense if we think of it already as a function. map also won't work, and the fusing operations are a little odd

# get f from an interpolation
h(x) = 2x;
g(x) = x.^2 + h.(x) + f.[x] #Why the asymmetry?

... all of this is to say: I think you should simply support the ( ) calling of the interpolation function anywhere you support the [ ] calling and in the exact same way. You could only turn it on for gridded interpolations, but my guess it is better to just make them synonymous.

[tomasaschan]

In Julia, functions that map integers, i.e. f : N -> R are usually called by the square brackets, f[2] etc.

This is, I think, not exactly the common interpretation of f[2] in Julia-land. Rather, square brackets are used to indicate the operation of indexing into a matrix-like thing. It would be more idiomatic to talk about A[i] than of
f[n]. (I realize that, from a mathematical viewpoint matrix indexing is still a function from N^n -> T, with T being the element type of the matrix, but in computer science, and in Julia in particular, those are regarded as different concepts.)

This package was built with the interpretation of an interpolation as an extension of the concept of an Array - an N-dimensional, matrix-like object which happens to support indexing with all real numbers, and not just integers. Thus, the square brackets make more sense than regular braces.

That said, it's probably trivial to extend interpolation objects with function call syntax outside of the package - just define it so it forwards to getindex and everything should work as expected.

[jlperla]

Yes, I understand when A[i] has indexing like behavior, but a number of features of it have nothing to do with filling in integer indices. For example, https://github.com/JuliaMath/Interpolations.jl#gridded-interpolation calls itp[4,1.2] That doesn't make much sense for array indexing? It does, however, make plenty of sense to me as itp(4, 1.2)

Also, the notation in the examples

tfine = 0:.01:1
xs, ys = [itp[t,1] for t in tfine], [itp[t,2] for t in tfine]

seems extremely odd to someone used to working with functions. how can you index 0.02 into an array?

Furthermore, it isn't really trivial to extend the function call outside of the package do the complexity of generic programming. Certainly not for a typical user (like the students of mine I would like to use this package). In every other function approximation library I can find, you call function-like things as functions, broadcasting as functions, etc. and it is very confusing that this one doesn't support that.

Is there harm in just adding in the correct, no-overhead forwarding of itp[1.1, 2.1] to itp(1.1, 2.1) for those who think of interpolation as creating an interpolated function?

[jlperla]

For example of all the other libraries I can find that interpolate and use functions:

• http://juliaapproximation.github.io/ApproxFun.jl/latest/faq.html#How-do-I-interpolate-a-function-at-a-specified-grid?-1
• http://quantecon.github.io/QuantEcon.jl/latest/api/QuantEcon.html#QuantEcon.LinInterp
• https://github.com/kbarbary/Dierckx.jl
• https://github.com/EconForge/splines.jl/blob/master/cubic_splines_julia.ipynb
  etc.

[tomasaschan]

Replacing the array-indexing style of evaluating interpolation objects would be a crazily breaking change, and it would require a much better motivation than "it feels more natural to me" (because that argument can be, and is, used by both sides). So, to be frank, that is not going to happen anytime soon, and quite possibly never.

What is still on the table is extending the current behavior with function call syntax in a backwards-compatible manner. It's been a while since i wrote some serious Julia code, but from my interpretation of this section of the manual it shouldn't require many lines of code to write something that treats interpolation objects like callable functors, and basically defines itp(args...) = getindex(itp, args...), punting all of the actual implementation to the existing library code.

Due to the extensible nature of Julia, this is just as easy to do outside of the package as inside it, and my guess is that it would, in total, only be a handful of lines of code. My suggestion is that you try to make that work in your own code - it shouldn't require any deep knowledge of the library, just that you want to make the function call forward to getindex, which you already know now. If that turns out to solve your problems, a PR for including that code in the package (maybe as a submodule that is optional to import) definitely has a good chance of being accepted.

[ChrisRackauckas]

Due to the extensible nature of Julia, this is just as easy to do outside of the package as inside it, and my guess is that it would, in total, only be a handful of lines of code. My suggestion is that you try to make that work in your own code - it shouldn't require any deep knowledge of the library, just that you want to make the function call forward to getindex, which you already know now. If that turns out to solve your problems, a PR for including that code in the package (maybe as a submodule that is optional to import) definitely has a good chance of being accepted.

It shouldn't be outside of this package because that could be considered a form of type-piracy and would (at least should) be blocked in METADATA. And IMO it's weird that the call doesn't work because it means that any things which you would use an interpolation of data for, like for approximating integrals via quadgk, don't work without putting anonymous functions over it to exactly do (x)->interp[x].

I think it would be quite a usability enhancement to just define the call. It would just be:

(b::BSplineInterpolation)(args...) = b[args...]

etc. for each concrete interpolation (it cannot be done on an abstract type).

[tomasaschan]

It shouldn't be outside of this package because that could be considered a form of type-piracy and would (at least should) be blocked in METADATA.

Again, it's been a while since I was actively involved with the Julia community, but this is very contrary to the sentiment that was common back when (as in, up to a couple of years ago). What changed? Are there other instances of this discussion, that I could read up on? Other packages or situations where this has been discussed and/or discouraged?

Also, note that I'm absolutely not saying it's never going to be in the package (lots of negations there, sorry...). What I'm saying is that

1. it's not difficult to define whatever you need yourself, since it requires so little code.
2. there's no technical limitation that says it has to be inside the package, so you don't have to wait for it to be implemented here before you can make use of it.
3. I think there's value in trying the concept out outside of the package, before deciding if it should be included, especially since it's so easy to try the concept out outside of the package
4. it has significantly larger chances of actually being included if someone files a PR that implements it
5. item 4 is especially true if the sentiment of the issue filed is more like "your library doesn't make sense! fix it!" rather than "wouldn't it be nice if...?"

[ChrisRackauckas]

Again, it's been a while since I was actively involved with the Julia community, but this is very contrary to the sentiment that was common back when (as in, up to a couple of years ago). What changed? Are there other instances of this discussion, that I could read up on? Other packages or situations where this has been discussed and/or discouraged?

No, type-piracy has never been encouraged. Write your own functions on other people's types, you can extend other people's functions with your own types, and mixtures like that are fine. It's very much discouraged to extend functions (or calls) that you don't own with types that you don't own, since that form of type-piracy can break other people's code without warning (there are some weird examples you can cook up). It's the same as in other languages where it's generally not a good idea to monkey patch some other person's package.

there's no technical limitation that says it has to be inside the package, so you don't have to wait for it to be implemented here before you can make use of it.

There is a technical limitation that if two packages do this outside of the package, and do it differently, then they will overwrite each other, which is why this kind of thing is discouraged / not allowed with registered packages.

[tomasaschan]

OK, that makes sense - thanks for explaining.

There is a technical limitation that if two packages do this outside of the package, and do it differently, then they will overwrite each other, which is why this kind of thing is discouraged / not allowed with registered packages.

That's not what I had in mind - I say you (or @jlperla) could do this in your own scripts, to see if it solves the problems you hope it will. If it does, feel free to file a PR. But there's no need to include it in any package - this or another - to be able to use it in a program.