DM-13065: SipApproximation

[TallJimbo]

Best to start with the docs for the afw/geom/SipApproximation class and the afw/math/polynomials namespace (in afw/math/polynomials.h) to get the big picture.

[r-owen]

I'm going to postpone the rest of the review until I can see DM-13156 and DM-13157 in JIRA (which is down at the moment).

[r-owen]

every -> ever

[r-owen]

I have never seen this done before. Clearly you could make this a real class by using template parameters. I assume you are trying to avoid inheritance, as per your documentation in polynomials.h ("the classes... are not polymorphic (no virtual functions)"). A few more words of explanation here would be appreciated.

[TallJimbo]

Concepts are very much a part of C++ - they're not in the language yet, but they've been a part of how people think about the language since the inception of the STL (see, e.g. http://en.cppreference.com/w/cpp/concept). Sometimes people try to implement those with a templated base class (i.e. CRTP, as in the Eigen class hierarchies), and I tried that a bit here, but I ultimately rejected that approach because I felt the associated boilerplate obfuscated things more than it helped.

But this way of documenting a Concept by telling Doxygen that it's a class is something new; it seemed the best way to provide a separate documentation page with the right sections that had more-or-less the right relationship to the concrete classes that implement the Concept. I'm pretty happy with the result, but documentation is best judged by someone other than the author and I'm open to other ways of documenting this. Given how focused it is on C++, it's a bit strange that Doxygen doesn't include something native for Concepts.

[r-owen]

A few words in the documentation would be helpful -- at least a link to the cppreference page you site. I had no idea the word "concept" was formal in this note.

[TallJimbo]

Good idea. Will do. I'll also capitalize Concept in the docs, which I should have done before.

[r-owen]

I strongly suspect you should put this code (including useful functions based on it) in a new package, rather than afw, since it seems potentially generally useful and has a minimum of dependencies. I'm thinking along the lines of your RFC for using sphgeom and I'd like to avoid a change like that in the future.

[TallJimbo]

Working on this was part of the inspiration for that RFC, but because some of this code does depend on afw.geom it can't be moved out of afw yet, and I'd rather not split it up based on that dependency problem; I'd rather address the dependency problem by moving both the geometry primitives and this code to another package. But since that's a lot more work, I can't make a good case that it's what LSST should spend its time on, and hence the RFC.

[r-owen]

More comments. Still working my way through this.

[r-owen]

I know this is picky, but I suggest "can support" -> "supports". I first misunderstood this as the maximum order you can request for this class, partly because I think I recall reading about other classes in this ticket that are vaguely similar but whose size can be increased if needed.

Also, if the workspace may contain data for a smaller order than the maximum, is it reasonable to call this getMaxOrder in hopes of perhaps eventually having getOrder return the current order.

[r-owen]

void -> avoid

[r-owen]

If I understand correctly, this gives a linear way of going through x and y coefficients. I know of two standard orders for doing this: x0y0 x1y0 x0y1 x2y0 x1y1 x0y2.... and x0y0 x0y1 x1y0 x0y2 x1y1 x2y0.... I believe we have different bits of code in our stack that uses each of the two orders. The point is: please document which order you are using here, and anywhere else this ambiguity comes up. I hope whatever it is will match the existing FunctionLibrary.h code and will become our de-facto standard.

[TallJimbo]

This is documented in the classes in PackedIndex.h; I've added a reference to that documentation here. The convention is the latter of the two examples you gave. I'm not sure what one FunctionLibrary uses, but this is the one used by lsst.shapelet and (internally) by ChebyshevBoundedField.

[r-owen]

Given that pre-increment is pretty much guaranteed to be faster (especially since post-increment calls pre-increment), is it necessary to define post-increment? I'm not complaining about this code (which is nice and simple); it is more of a C++ style question.

[TallJimbo]

It is necessary to meet the formal requirements of being a conformant STL InputIterator.

[r-owen]

These aren't Chebyshev polynomials, right? I suggest looking for all instances of Chebyshev in this file and the 1d file.

[r-owen]

I encourage you to try to find a less ambiguous name, such as unscaled. (Note that Transform supports simplify for a very different purpose.)

It is a pity it is necessary, but I imagine it was needed in order to compute the final SIP coefficients.

[TallJimbo]

@r-owen, I think this is ready for another look.

I actually took some time to look more closely into why unscaled involved a loss of precision, including adding a "safe sum" algorithm that should avoid any loss of precision beyond the truncation imposed by the actual result type (you still can't compute more digits than fit in a double, of course). The results were a bit surprising: the loss of precision in unscaled really isn't that bad, even without the safe algorithm. The safe algorithm is slightly better, so I've kept it around (especially as I'm much more confident in its worst-case performance), but the bottom line is that all of the loss of precision in the high-level test in afw happens in the fit of the unscaled polynomials (which we can't use the safe sum on, because it's a matrix inversion, not a simple sum), rather than the conversion to unscaled polynomials. The 50*DEFAULT_RTOL seen in the tests is scarier than it looks - note that this is still an accuracy of ~1E-15. As a result, I don't actually see a problem with renaming unscaled to simplified, and that's what I've done here.

Also, note that there is an afw PR for this, with mostly obsolete comments dating back to when the polynomial was in afw. Please take a look at that, too.

[r-owen]

I have some suggestions to improve interoperability but overall this looks very nice.

Do you anticipate any syntactic sugar for making a WCS from one of these (e.g. a new function)?

[r-owen]

Outstanding documentation. Very clear and helpful.

[r-owen]

Please state that crpix is 0-based (I sure hope it is!). In case it helps SkyWcs::getPixelOrigin() returns this value and its documentation string is:
``
/**
* Get the pixel origin, in pixels, using the LSST convention
*
* This is CRPIX1 - 1, CRPIX2 -1 in FITS terminology
*/

[r-owen]

Consider "set up" -> "generate" (or "build" or "construct" or...)

[r-owen]

I'm uneasy about this. It's trivial to make a SkyWcs from SIP coefficients (call makeTanSipWcs, though if you want a version that handles arbitrary projections we would have to add that, it should be trivial). Instead of having this code I suggest you make a WCS that way, extract the pixelToIwc portion by calling getPixelToIntermediateWorldCoords, and call the regular fitter with the crpix, cd matrix and SIP coefficients you used to make the WCS.

[TallJimbo]

The purpose of this is to be able to test that SipApproximation internal evaluate identical transformations the same way that SkyWcs does (which is something that needs to be verified before testing that it evaluates approximation transforms similarly), so I can't really use SkyWcs as a replacement for this.

[r-owen]

"a point" -> "an array of points"

[r-owen]

Please document behavior for invalid inputs, e.g. order < 0 (if 0 is sensible), gridDimensions <= 1 (especially <= 0) and an empty bbox. I'm hoping you raise an exception, but I suppose an assert would do.

[r-owen]

The units "scaled IWC" seem a bit odd. However, I assume this is naturally what comes out of the computation, and if so, I agree it's probably the best choice (especially since in the standard case it'll be in reasonable units).

[TallJimbo]

Yup, that's the natural result of the computation, though putting all of the deviations in pixel units was actually the motivation.

[r-owen]

makeSkyWcs and makeTanSipWcs use Eigen::Matrix2d const &cdMatrix and the latter uses the same type for the SIP coefficients. Please do the same unless there is a very strong reason not to. It would improve interoperability. For instance the user could pass in the output of wcs.getCdMatrix() as the cd (or, preferably, cdMatrix) argument and use the generated SIP directly when building a TAN-SIP WCS.

[TallJimbo]

FWIW, I think we should be probably using LinearTransform in all of those places, but I'll bow to local precedent here.

[r-owen]

Please add methods to return the coefficients as Eigen matrices

[r-owen]

Please add brief documentation of what this does. Why this name instead of, say, "makePolynomial"?

[TallJimbo]

afw changes from the last round of review comments are now ready too.

Do you anticipate any syntactic sugar for making a WCS from one of these (e.g. a new function)?

Not immediately; I'm not entirely certain how we'll want to use this class to construct a new WCS in detail - after initializing a SipApproximation, you might fit multiple times and change the grid size or order a few times before deciding you're happy with the results. We'll need to find some reasonable default logic for that, of course, but it's not something I have time to do on this ticket.