Function call syntax support

[chiyahn]

Here are a couple of changes I made in the code for function call syntax support to address issues raised in #198. Unit tests (function-call-syntax.jl) have been added to ensure that this extra feature does not cause conflicts with a previous version.

[jlperla]

It is worth seeing if the v0.7 failure there is similar to the one for Interpolation.jl as a whole, or if you need something more. If there is an issue, then maybe look into Compat.jl

[jlperla]

@Nosferican @sglyon @ChrisRackauckas

Hi guys, can I get you to take a look at this PR (and maybe see if you know someone who could merge it if you are happy with it)? The hope is that we can add in enough convenience wrappers to ensure that we can teach economics lectures and the notation not be absurdly inconvenient compared to alternatives.

[ChrisRackauckas]

The implementation looks good to me. It's really just forwarding function calls to getindex and adding tests.

[jlperla]

Another reason to consider this. I don't think that broadcasting is possible with the [] access, or it is highly non-standard... but it definitely is possible now. i.e.

using Interpolations
x = 0:0.1:1
y = x.^2;
itp = interpolate((x,), y, Gridded(Linear()))
itp[0.52]
function (itp::Interpolations.GriddedInterpolation{T,N,TCoefs,IT,K,pad})(args...)where {T,N,TCoefs,IT,K,pad}
    itp[args...]
end
itp.([0.52 0.53]) #Broadcast!

[sglyon]

Implementation is also good.

I have no hesitation about merging this implementation of the feature (good work @chiyahn ). The only question in my mind is whether or not we want to merge this feature.

I'd like to wait to hear from either @timholy or @tlycken before we make the final call on this one

[ChrisRackauckas]

That is correct that you cannot broadcast indexing, at least without making it a getindex function call.

[pkofod]

sorry if I'm being dense, but do you really need to specify all the parameters? wouldn't

(itp::TheITPType)(args...) = itp[args...]

work?

[tomasaschan]

A lot of water has passed under the bridges since I was last really involved in the community, and it seems to me that while at that time array indexing was the thing that got a lot of love, function calling has gotten at least as much love since. I don't feel that I have the insight to make a final call on this one, but I have no really strong feelings against it. Both ideas (array indexing vs function calling) derive from valid, albeit different, mental models of what interpolation is and how it works, so neither is inherently unsuitable.

What I do feel, though, is that it's probably a bad idea to support both [] and () for any extended period of time. The expected semantics of them are pretty different, and I wouldn't be surprised if merging this without making any other changes would soon result in issues from users who find the behavior of the library confusing in some cases, because it doesn't behave the way they expect callable objects to behave. OTOH, accommodating those cases could make the array indexing behavior confusing instead; we don't expect them to behave exactly the same.

Therefore, I'm reluctant to recommend merging this feature without a plan for how to pivot the entire library to function call semantics instead of array indexing semantics. This can, naturally, not be done overnight (and there would probably have to be two major version bumps involved - one after which we deprecate [], and one more before we remove it) so I don't expect it to be implemented before this merges, but I would like there to be some idea of a plan for how (and perhaps when and/or who)...

[ChrisRackauckas]

I honestly don't see how having both is confusing. Mention it in one spot in the docs and it's fine without any more.

[tomasaschan]

@ChrisRackauckas As I said, I don't want my words to be a final call in any direction. I don't want to stand in the way of the community here - but my opinion was asked for, so I gave it. Feel free to ignore me :)

But to further explain my reasoning: if you see A(x, y) and B[x, y] in a script, do you expect them to follow the same semantics? Are there any cases where they do not? Do you expect to be able to refactor one in ways you don't expect to be able to refactor the other? If there are no such cases, then all my points above are moot :)

[jlperla]

But to further explain my reasoning: if you see A(x, y) and B[x, y] in a script, do you expect them to follow the same semantics? Are there any cases where they do not?

Semantics matter, and to tell you the truth I am not sure I see any examples I would use that have the B[x,y] semantics (in the economics domain, at least) unless x and y were always countable. So going back in time I would say that () is the right semantics in all cases where you can call the function with a non-integer, as f[2.2] doesn't make any semantic sense to me but f(2.2) does.

Since we can't go back in time, I would say that the choices are (1) allow both in all cases; (2) allow both in any case where the argument can be a real; (3) remove the [] for any case where the argument can be a real so that () becomes the only semantics; and (4) do nothing. I think the best solution is probably (3), but it is a breaking change and a lot of work. (1) is good enough for now and can be revisted.

For, (4), which that means is that you are effectively saying that this library cannot be the default interpolation library for a large problem domain with incoming users with Matlab-level programming experience. To give a few concrete examples:

• I recently suggested someone swap out for the Interpolations.jl library, but the only way it was possible was to tell them to add in crazy boilerplate code that would scare off someone coming from matlab (see https://discourse.julialang.org/t/how-to-call-a-function-that-returns-a-function-efficiently/11851/22?u=jlperla )
• Try to explain https://lectures.quantecon.org/jl/julia_libraries.html#univariate-with-irregular-grid where you call all functions with (), every other interpolation library with () but then you have to use [] for Interpolations.jl.... Makes generic programming very difficult, since you can't pass around interpolation objects as "functions".
• To add the functions yourselves, it is unequivocal type piracy of the worst kind. This could cause confusing grief if other libraries they use do the same thing by just copying the () code around, and perhaps mess up the generic type parameters...

[jlperla]

More generally, in order for Interpolations.jl to be a solid foundation for economic applications, I have identified the following issues:

• Unlike every other library, you cannot call as a function. (this fixes it).
• Performance issues compared to other implementations. For example, the person who tried Interpolations.jl found that Spencer's LinInterp in QuantEcon was faster. Hopefully that can be fixed.
• Convenience constructors are necessary, or else the code is extremely confusing. See https://lectures.quantecon.org/jl/julia_libraries.html#univariate-interpolation
• We really need an implementation of cubic splines on irregular grids. See Gridded Cubic Splines (WIP) #193 for the start of the implementation.
• I tried it and suspect Interpolations.jl either doesn't support static arrays but I think it is an important thing to implement efficiently since many cases have interpolation over relatively small arrays (e.g. discretized markov chains)

With those, I am happy to tell my students and all economists I know to use the library. I also really don't like having the split implementation in of LinInterp in QuantEcon, and think it should be deprecated in favor of a cleaner Interpolations.jl implementation and am willing to argue with spencer, john, et al with that goal. I suspect spencer already smelled this coming :-)

But, more importantly, I am willing to finance grad students and RAs in order to make these things possible with those goals in mind for teaching this fall.

[tomasaschan]

I think the best solution is probably (3), but it is a breaking change and a lot of work.

I agree; this is the gist of what I was trying to say. Most of the rest of it was arguing against choosing (1) (or, IMO even worse, (2)), as a long-term plan.

BUT - since it's a lot of work (I know) and I'm not in a place where I have the time to put in said work (I wish...) I don't think it's reasonable to expect all this work to be in place just to open up all the possibilities the () syntax would allow. Therefore, accepting (1) for now with an explicitly stated goal of moving toward (3) is probably good enough.

A simple enough plan for how to do that:

1. Merge this PR, and tag a new minor release. This can be done now (more or less - see comment at the end...)

2. Deprecate [] in favor of (), and tag a new major release. Migth need some extra thought with deprecation messages to suggest broadcasting instead of array indexing, etc, since such things are currently handled by AbstractArray stuff in base.

3. Re-write the library to not support [], but only use () instead. Might be as easy as replacing the getindex definitions with corresponding function call implementations, but might also require some thought about how to handle array indexing etc (maybe broadcasting takes care of that?).

4. Tag a new major release, with the deprecations removed.

Since we haven't even tagged 1.0 yet, we might want to do that first. This library is IMO mature enough that it should probably be 1.0 already. (Maybe Julia 0.7 compat issues, see e.g. #204 should be solved first?)

[jlperla]

I agree completely.

For tagging, I was hoping to sneak some stuff into the library over the next month or so, a few thoughts on what could happen

1. Convenience constructors: In order to increase the simplicity of using the library, add in a few helpers to make simple things simple (i.e. univariate linear interpolation with regular and irregular grids could have a very clean interface, etc.) At this point the construction is very general and powerful, but beyond the cognitive scope of a matlab user for simple things (myself included... I can't figure out most of it :-) ). We would then add documentation and an example/docs
2. Performance of Univariate Linear Interpolation: Make the linear interpolation comparable in performance to QuantEcon's LinInterp, or else it is tough to convince Spencer and John to use it.
3. Get v0.7 supported
4. Get Gridded cubic splines working, if perhaps not full optimized yet.

If you agree (in principle, if not in the details) to all of this, then I can get @chiyahn to start on executing these with the hope that we have an easy to teach, higher performance library for v0.7 this fall.

[rdeits]

I've been using this package for a couple of years, and I would certainly support the use of () notation. As posted earlier, the biggest benefit of doing so would be integration with Julia's really excellent broadcasting, but I think where this really shines is with broadcast fusion, so that something like:

2 .* itp.(x .+ 1)

would both (a) work for an arbitrary array x and (b) perform exactly one pass over the data x and allocate no intermediate arrays. Getting the same result with the current [] notation requires something like:

2 .* getindex.(Ref(itp), x .+ 1)

[stevengj]

I agree with @tlycken that there should be one syntax, not two — this isn't Perl 😉. Why not @deprecate the getindex syntax in favor of the call syntax, which is arguably more natural here?

[jlperla]

@stevengj While I agree in principle, I don't have time to support RAs doing a full overhaul on this deprecation right now, as it looks like it would be an enormous change... I would much rather use the existing patch with the plan to deprecate later so that we can add in other essential features for the library to be usable (e.g. convenience constructors, performance improvement)

However, what if we:

• Change the https://github.com/JuliaMath/Interpolations.jl/blob/master/README.md to only use the () notation
• Put a comment at the top of that saying [] will be deprecated in the future

Is that enough to get through? There is just too much other work to prepare for teaching with v0.7 this fall otherwise...

[sglyon]

I am on board with merging this, under two conditions:

1. at least one other maintainer gives the green light.
2. We create an issue that describes what/how/when should be depcreated in light of the discussion above

[jlperla]

@sglyon Sounds great. It looked like other maintainers gave the "thumbs up" to the plan in #206 (comment) which @chiyahn will execute. We will also add in a "Deprecate []" as a separate issue describing roughly would be needed for the full deprecation. When this is done, we will ping this PR to get everyone's thoughts.

[jlperla]

@sglyon OK, I believe that @chiyahn has finished everything you and others have asked of him (assuming that the thumbsup in #206 (comment) was an endoresement of the approach)

As you requested, he also added in the #212 as a reminder (and with some notes) for the full deprecation. Doublecheck and merge if you agree! Next he is working on some simple, QuantEcon compatible, convenience constructors.

[tomasaschan]

Before this is merged, I think #213 brings up a counter-argument to making this transition in the way we talked about. In short, ForwardDiff.jl requires, among other things, that

The types of array inputs must be subtypes of AbstractArray . Non-AbstractArray array-like types are not officially supported.

By removing getindex we also effectively make interpolation objects not-arrays, meaning less compatibility with stuff that people seem to use this library for. Maybe we need to think a little harder on what the API should be like?

[stevengj]

By removing getindex we also effectively make interpolation objects not-arrays, meaning less compatibility with stuff that people seem to use this library for.

The question here is about supporting getindex with a continuous "index", which doesn't behave like an AbstractArray and probably won't work with other code expecting ordinary indexed containers anyway.

If you want to keep getindex only for a genuine array-like interface to the underlying discrete data, that's a different matter.

[jlperla]

By removing getindex we also effectively make interpolation objects not-arrays, meaning less compatibility with stuff that people seem to use this library for

For univariate functions, ForwardDiff.jl would work perfectly with the () interface (whereas before you would need to create a function to forward the [] to () on your own.

I think you are misreading the ForwardDiff.jl semantics. For multivariate functions, it auto-differentiates functions of AbstractArrays, not AbstractArrays themselves. This is a standard issue with any users of AD (or optimizers) that require things in vectors. The standard trick is to do something like (untested)

f(z1,z2) = z1^2 - z2
ForwardDiff.gradient(z -> f(z[1], z[2], [0.1, 0.2])

Where in the case of an interpolation, the f would just be the interpolation object.

Of course, down the road you could make this sort of thing easier by adding in a () overload for AbstractArrays in the interpolation, which does something to the effect of (interp::ConcreteInterpolationType)(x::AbstractArray) = interp(x...) I know that splatting isn't that easy, but hopefully that gets the point across. My feeling, though is that it isn't worth it until we see a direct need, but the new move to () makes it possible.

Hopefully with Capstan.jl it will be able to start using AD on functions with multiple arguments to prevent that stacking, but it has nothing to do with interpolation.

---

But all of this comes down to a simple statement: in julia, you call functions with () (stacked in an array or otherwise) and the sematnics [] is reserved for array-like interfaces where you pass in discrete arguments.

[tomasaschan]

I think you are misreading the ForwardDiff.jl semantics.

That might very well be true; if so, carry on and don't mind the confused guy on the sidelines :)

[sglyon]

Ok, I am starting the timer and will merge within the next 12-24 hours unless I hear back on this thread that it shouldn't happen.

[timholy]

Supporting () natively isn't really done yet, nor is it particularly simple, because currently we rely on Base's array-indexing machinery to support things like

itp[1:3, 5:7]
itp[CartesianIndex(2,3), 4, CartesianIndex(5,6)]    # expands to itp[2,3,4,5,6]

I will dig a little and see how much we can co-opt that machinery without having to duplicate it.

[sglyon]

Thanks @timholy I'm appreciative that you are working on this. I'll readily admit that you understand all the benefits/tradeoffs of relying on getindex machinery better than I do!

[timholy]

np, you made the right call by merging it.