Write polynomial interpolation class

[rmjarvis]

@joezuntz has already started working on this on branch joe_polynomial. (See also branch polynomial_interp that @danielgruen started work on, but he said he thinks all of that work is superseded by what Joe has done.)

Not much to say here. Should take an arbitrary order for each coordinate. Maybe options for the kind of polynomial? Cheby, Lagrange, etc. That might not be important though.

[joezuntz]

Will do.

[joezuntz]

@rmjarvis I just committed a big chunk of this with a test collection.

It turns out to be trivial to do all the other polynomials, as numpy has a whole collection of modules with very similar interfaces making Hermite, Chebyshev, Laguerre, and Legendre polynomials alongside the standard module.

I've put a collection of tests in - they're a bit disorganized and there's definitely functionality that's not tested yet. They helped pick up a bunch of errors! Any advice or comments v welcome.

It's currently pegged to use the same polynomial order for all the parameters and for the u and v directions. I could make it more flexible easily if there is a reason for this - my instinct is that u and v should have the same order but not necessarily all the parameters, but do you have any thoughts?

There's a few other TODO's in there also:

• I'm using numpy arrays for the positions - you're using galsim positions elsewhere I think?
• Probably want to normalize the coordinate ranges to keep all the coefficients the same order of magnitude.
• Could refactor to pull out all the polynomial-specific bits to a subclass and the generic model-fitting bits in the base.

One idea: a nice milestone for Piff would be being able to run something fairly PSFEx-like on a single DES exposure file: polynomial interpolation and a pixel-wise star model. It's a bit boring but it would be a nice test of the machinery. Does that sound good? I can add it as a github milestone if so.

[rmjarvis]

Great! Thanks Joe, I'll take a look later this weekend. I'm traveling today and tomorrow, and I won't have much time until tomorrow night probably.

my instinct is that u and v should have the same order but not necessarily all the parameters,

I think we should allow each parameter to have arbitrary order. With our 2k x 4k chips, we may want to have different order for u and v when doing one chip-at-a-time interpolation. So we should probably enable that option in general. And more generally, we may have other parameters that we might want a different order on. Like linear in color for instance.

a nice milestone for Piff would be being able to run something fairly PSFEx-like on a single DES exposure file: polynomial interpolation and a pixel-wise star model.

Agreed! Maybe even worth an official release when we get there. :)

[joezuntz]

I've enabled using different polynomial orders for different parameters. Using different ones for different position vectors will require a little more thought and work as currently we are specifying the max total power (i.e. if p(x,y) = sum c_{ij} x^i y^j we use all the values where i+j<=order, so we don't have a separate polynomial for the two of them).

[rmjarvis]

Anyone know what PSFEx does when you tell it to use different orders for the two directions?

One option that might make sense is to use uorder as the maximum order in u, vorder as the maximum order in v, but also impose max(uorder,vorder) as the maximum combined i+j of any term.

So if uorder = 2, vorder = 3, you would get { 1, u, v, u^2, uv, v^2, u^2 v, u v^2, v^3 }.

[joezuntz]

I think this is probably equation 14 of the PSFEx manual. I'm not sure I understand it as it's not clear if a distortion group is the same as a context vector, but I think it's saying that you have groups of parameters (e.g., u and v might be one group and the colours might be another) then within a group it just uses a maximum total exponent, but then between groups it just multiplies the polynomials together.

So I think the answer is "either, depending on what settings you use".

[joezuntz]

Ah yes, and below that it says that in the default configuration the two parameters (presumably location) are by default in the same group, so that it would use a max total exponent (the default max is 2).

[rmjarvis]

I see. So it doesn't have any way to get this hybrid thing that I suggested.

How about for now, we follow suit and just go with what you have (position components use a single order with other parameters having their own separate order), and we can discuss what extension to this we might want down the line in a separate issue.

[joezuntz]

okay, cool