Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

meas_astrom/DM-3549: new SIP fitter #47

Merged
merged 8 commits into from Oct 26, 2016
Merged

meas_astrom/DM-3549: new SIP fitter #47

merged 8 commits into from Oct 26, 2016

Conversation

TallJimbo
Copy link
Member

This adds a new fitter for SIP distortions in the form of a number of C++ classes (several polynomial-based transform classes and a fitter) and a new Python task that is intended as a drop-in replacement for FitTanSipWcsTask.

Some caveats:

  • I have not yet enabled the new task, FitSipDistortionTask, as the default; I think it needs a bit more testing before we do that. I do however expect it to fully replace FitTanSipWcsTask eventually, so I have not tried to more strictly define the interface for WCS fitters or add a registry as we do with other pluggable tasks. I also have ignored some (minor) code duplication between the tasks (especially their config classes) that would be better to clean up with common base class if we expected them to coexist for a long time.
  • The new transform classes do not inherit from XYTransform or any other common base class, and that probably limits their usefulness to this particular application. I think we'll want to refactor that when we have our new WCS/transform library, but doing it before then is probably premature.

Copy link
Contributor

@r-owen r-owen left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I have gotten through the .h and .i files. I'll review more tomorrow.

So far my main request is to document some of the private code.


namespace lsst { namespace meas { namespace astrom { namespace detail {

inline int computePackedOffset(int order) {
Copy link
Contributor

@r-owen r-owen Oct 5, 2016

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please document this and all other functions and classes in this file. This one is not so obscure, but I have no idea what some of the things in this file do.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Done (locally; will commit and push after aggregating minor comments).

namespace lsst { namespace meas { namespace astrom { namespace detail {

inline int computePackedOffset(int order) {
return (order*(order+1))/2;
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Space around +?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Done.

* @param[in] bbox Bounding box for the grid in the input
* coordinates of the toInvert transform
* (or equivalently the output coordinates
* of the transform to be fit).
Copy link
Contributor

@r-owen r-owen Oct 5, 2016

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can you provide a concrete example of what the bbox will be for the typical case? Will it be the bbox of the exposure on which the sources were measured? The smallest bbox that encloses all sources after removing distortion? Or? Put another way: can you suggest how the user would usually compute this bbox?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In the way it's used in FitSipDistortionTask, the bounding box has to be in intermediate world coordinates, and given that we only start with the bounding box of the image and we're trying to find the mapping from pixel coordinates to intermediate world coordinates here, our only choice is to transform the image bounding box using (part of) the initial guess WCS.

I don't think it's appropriate to include that description in the docs for fromGrid, since it's more general than that, but I will add code comments to FitSipDistortionTask where I'm computing the bbox.


/**
* Perform a linear least-squares fit of the polynomial coefficients.
*
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Could you say a few words (perhaps in the class description) how you might call this to perform a fit, e.g. when to call fit, updateIntrinsicScatter and rejectOutliers, in what order, how often, etc.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Done. I've added a bit to each of these functions and a code block to the class docs describing the typical iteration pattern.

/**
* Return the current estimate of the intrinsic scatter.
*/
double getIntrinsicScatter() const { return _intrinsicScatter; }
Copy link
Contributor

@r-owen r-owen Oct 5, 2016

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Could you say a few words (preferably in the doc for this function or the previous) why you have separate updateIntrinsicScatter and getIntrinsicScatter methods? Naively I would expect you could manage with just one of them. Forgetting to call the former before calling the latter seems like an accident waiting to happen.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The intrinsic scatter is actually included in the uncertainty on the data points, so one might want to call this before updateIntrinsicScatter to obtain the scatter used in the last call to fit or rejectOutliers. Calling it after updateIntrinsicScatter tells you the scatter that will be used in the next call to fit or rejectOutliers. (I've added that to the docs for getIntrinsicScatter).

* Convert a PolynomialTransform to an equivalent SipReverseTransform.
*
* @param[in] poly PolynomialTransform to convert.
* @param[in] pixelOrigin CRPIX @f$(u_0,v_0)@f$
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Again, anywhere you mention CRPIX please be explicit about the zero convention (that's safer than putting it in one place and hoping folks stumble across it).

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Done.

};

/**
* Create a new TanWcs from a pair of SIP transforms and the sky origin.
Copy link
Contributor

@r-owen r-owen Oct 6, 2016

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I suggest you say "TAN SIP WCS" instead of a "TanWcs"; the latter is a poorly named class that should go away; the former is a more general concept that should remain after the WC rewrite.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Done.


%include "lsst/meas/astrom/matchOptimisticB.h"

%include "std_pair.i"

namespace std { // have to do this in namespace block so Swig can find size_t
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

weird. You can't just use std::size_t?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Correct. That was what I tried first. I think I could have also fixed it by not using std:: in the C++ header file, but I'm not sure that's valid C++ even though my compiler allows it.


%include "lsst/meas/astrom/matchOptimisticB.h"

%include "std_pair.i"
Copy link
Contributor

@r-owen r-owen Oct 6, 2016

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If this could safely go higher (e.g. right near %include "ndarray.i" I think it would make the code a bit easier to read; we usually put standard and 3rd party includes together near the top when we can get away with it.

Another option worth considering is moving the new code into a separate .i file. There is a fair amount of it, so it may enhance readability.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'll refactor into a new file.

%template(DoubleSizePair) pair<double,size_t>;
}

%declareNumPyConverters(ndarray::Array<double const,2,2>)
Copy link
Contributor

@r-owen r-owen Oct 6, 2016

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Here's another line I suggest that you move way up with the standard imports (and typically vector templates and such) if that's safe to do.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I was actually planning to put this in the custom .i file for the transform classes, since it's only needed by them.

Copy link
Contributor

@r-owen r-owen left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Python code reviewed. One debug "or True" needs to be removed, but otherwise just minor stuff, including requests for clarification.

default=1
)
maxScatterArcsec = lsst.pex.config.RangeField(
doc="maximum median scatter of a WCS fit beyond which the fit fails (arcsec); " +
Copy link
Contributor

@r-owen r-owen Oct 6, 2016

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

(picky) please caplitalize "Maximum" for consistency with the other doc strings.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Done (also capitalized "order" and "number" for the same reason).

)
refUncertainty = lsst.pex.config.Field(
doc="RMS uncertainty in reference catalog positions, in pixels. Will be added " +
"in quadrature with measured uncertainties in the fit.",
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please consider indenting continuation strings, for improved readability.

default=100,
)
gridBorder = lsst.pex.config.Field(
doc=("When setting the gird region, how much to extend the image "
Copy link
Contributor

@r-owen r-owen Oct 6, 2016

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please remove the parens as they are unnecessary and confusing (hinting that doc is being set to a tuple). I suggest using the same format for all multi-line strings, with a slight preference for the format you used above with "+" at the end of each line, since the "+", though optional, makes the intent clearer.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Done.


FitSipDistortionTask is a drop-in replacement for
:py:class:`lsst.meas.astrom.FitTanSinWcsTask`. It is built on fundamentally
stronger fitting algorithms, but has received significantly less testing.
Copy link
Contributor

@r-owen r-owen Oct 6, 2016

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If you have not already done so, please file a ticket to remove this caveat and make this the default fitter once it has been adequately tested. Otherwise this comment may live long past its relevant date.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.


- We set the CRVAL and CRPIX reference points to the mean position, while
holding the CD matrix fixed to its input value. This work is done by
the makeInitialWcs method.
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Summarizing what this means in English would be helpful. I think you mean that you adjust the mean position of the WCS. That said, I'm surprised that you choose to adjust both CRVAL and CRPIX instead of just one of them. I would have thought one would be treated as a given.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I thought this was a summary, though I see that I just said "mean position" when I should have said "mean position of all the matches". Does that clarify things?

I think we could add a config option to leave CRVAL and/or CRPIX fixed at the positions from the given initial WCS, but this matches the behavior of FitTanSipWcsTask and I didn't see anything wrong with that.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If it matches current behavior then I agree there seems no strong need to allow keeping them fixed.

# Fit the forward mapping to a grid of points created from the reverse
# transform. Because that grid needs to be defined in intermediate
# world coordinates, we first grow the bounding box and then apply the
# initial WCS's CRPIX and CD matrix.
Copy link
Contributor

@r-owen r-owen Oct 6, 2016

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't understand this explanation, starting from "Because that grid...". What is the connection between using intermediate world coordinates and needing to grow the bbox? Or is it that you intended to grow the bbox all along, but the sequencing of growing first, then applying... is driven by using intermediate world coordinates?

Also, I suggest replacing "needs to be" with "is" (or similar) because it is more of a choice (albeit an important choice driven by experiment) rather than a fundamental requirement.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm growing because I'm worried the transform I use to generate the box is only approximate, and I'd rather get a too-big box than a too-small one. I'll clarify.

# initial WCS's CRPIX and CD matrix.
gridBBoxPix = lsst.afw.geom.Box2D(bbox)
gridBBoxPix.grow(self.config.gridBorder)
gridBBox = lsst.afw.geom.Box2D()
Copy link
Contributor

@r-owen r-owen Oct 6, 2016

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please consider a longer name to better differentiate gridBBox from bbox, since they are in different coordinate systems; perhaps gridBBoxIwc.

If provided then coords are updated with the new WCS;
otherwise only the coords for sources in matches are updated.
Required fields are "slot_Centroid_x", "slot_Centroid_y",
"coord_ra", and "coord_dec".
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Don't you also need centroid uncertainty?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Good catch; fixed.

sErr = record.get(srcErrKey)
sEllipse = lsst.afw.geom.ellipses.Quadrupole(sErr[0, 0], sErr[1, 1], sErr[0, 1])
disp.dot(sEllipse, sx, sy, ctype=colors[1])
print("Dropping into debugger to allow inspection of display. Type 'continue' when done.")
Copy link
Contributor

@r-owen r-owen Oct 6, 2016

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why not just display each iteration in a new frame without stopping? It seems less disruptive.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If you think that's the more common pattern, I'm happy to change it. This is the first time I've actually written interactive display code in a pipeline task.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It's a slightly delicate issue. Dropping into the debugger is fairly rare, but by no means unheard of. However, the current display code handles frame numbers rather poorly because it does not support any kind of auto-increment. It's easy enough to keep track of a local frame number, but if later code wants to display something in a different frame...

Which reminds me: display is often an integer that states the starting frame number (though if we had auto-increment we'd rarely want that).

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I've added separate lsstDebug variables for the frame and whether to pause (I've also seen that as a way to specify the (starting) frame).

matches with no estimated Wcs, and the input Wcs is a wild guess.
We're using the best of each: positions from the matches, and scale
from the input Wcs.

Copy link
Contributor

@r-owen r-owen Oct 6, 2016

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This really surprises me. The initial WCS ought to include an estimate of optical distortion, and so ought to be pretty decent. Still...I assume you found you needed to do this instead of trusting the initial WCS. I just wish I knew why.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I actually just stole this method (including that doc text) from the old WCS fitter task. If it's no longer accurate, I'll change it.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I've softened the docs here to make it clear that we're being defensive rather than explicitly distrusting the input, and that it's the position we're least certain about.

std::swap(_pixelOrigin, other._pixelOrigin);
std::swap(_cdMatrix, other._cdMatrix);
_poly.swap(other._poly);
}
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Any particular reason you are not using ADL for swap?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I have to admit I'm not entirely clear when it's advantageous to do so; just I figured I didn't need it here. Are you advocating it just because you think it's a good pattern to follow, or is there a benefit I might not be aware of?

ndarray::EigenView<double,2,2> _xCoeffs;
ndarray::EigenView<double,2,2> _yCoeffs;
mutable Eigen::VectorXd _u; // workspace for operator() and linearize
mutable Eigen::VectorXd _v;
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thread safe?

Copy link
Member Author

@TallJimbo TallJimbo Oct 7, 2016

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hmm., I think I'll just document that instance of this class should be confined to a single thread, as I think that's probably how it ought to be used. Ultimately I expect it to be a moot point when we switch to using AST for these classes.

Edited: I was more careful than I remembered being: BinomialMatrix is already threadsafe.

Copy link
Contributor

@r-owen r-owen left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Finished the .cc files. Still have the unit tests.

PolynomialTransform PolynomialTransform::convert(SipForwardTransform const & other) {
PolynomialTransform poly = other.getPoly();
poly._xCoeffs(1, 0) += 1;
poly._yCoeffs(0, 1) += 1;
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Would you consider a brief comment as to why you are incrementing these coefficients?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Done. It's to account for the extra terms outside the sum in the SIP definition formula; you can fold those into the sum by adding one to these two coefficients.

throw LSST_EXCEPT(
pex::exceptions::LengthError,
"PolynomialTransform order must be >= 0"
);
Copy link
Contributor

@r-owen r-owen Oct 6, 2016

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If order < 1 then won't the ndarray::allocate(...) statements in the assignment list have caused problems by the time you get here? If so, then perhaps you should allocate here rather than in the assignment list.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Good catch. Done (need to move the Eigen initializations too, I think).

if (xCoeffs.getShape() != yCoeffs.getShape()) {
throw LSST_EXCEPT(
pex::exceptions::LengthError,
"X and Y coefficient matrices must have the same shape."
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please print the shapes

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Done.

if (_xCoeffs.cols() != _xCoeffs.rows()) {
throw LSST_EXCEPT(
pex::exceptions::LengthError,
"Coefficient matrices must be triangular, not trapezoidal."
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

trapezoidal? Aren't they rectangular? In any case, please print the offending values.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

For matrices, trapezoidal is to triangular as rectangular is to square. They're triangular because we don't want nonzero terms with p+q > order, though it's not worth out time to test that.

Offending values printed.

ScaledPolynomialTransform result(
sipForward.getPoly(),
afw::geom::AffineTransform(afw::geom::Point2D() - sipForward.getPixelOrigin()),
afw::geom::AffineTransform(sipForward.getCDMatrix())
Copy link
Contributor

@r-owen r-owen Oct 6, 2016

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I suggest Point2D(0, 0) (or 0.0, 0.0) so the reader doesn't have to look up what the default Point2D constructor does.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Fine; done.

#include "lsst/afw/table/Match.h"
#include "lsst/afw/math/LeastSquares.h"

#define LSST_ScaledPolynomialTransformFitter_TEST_IN_PLACE 1
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please document what this does. Do you really want it left true for production code?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Good catch. I do not want it left true for production code. I've added an extensive comment.

Eigen::ArrayXd fyy = 1.0/(syy - sxy.square()/sxx);
Eigen::ArrayXd fxy = -(sxy/sxx)*fyy;
#ifdef LSST_ScaledPolynomialTransformFitter_TEST_IN_PLACE
assert((sxx*fxx + sxy*fxy).isApproxToConstant(1.0));
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I suspect we will disagree about this, but this section is presently on by default and I'm very nervous about asserts being run in production code. I'd be much happier with raising an exception so the user gets a more useful message and other code has a chance to continue if the fitter fails.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hopefully this is resolved by the above change to disable this block entirely.

// Function that computes the -log likelihood of the current deltas with
// the variance modeled as described above.
auto logLikelihood = [this](double intrinsicVariance) {
// Uncertainties in the table right now include the old intrinsic scatter; need to
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is it important to use this syntax instead of defining a named function (perhaps in the anonymous namespace)? I think the latter would be easier to read.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

A lambda defined here has access to ScaledPolynomialTransformFitter's private members, which a free function in an anonymous namespace would not. A private member function would have access to those private members, but then wouldn't be usable as a functor as I could pass to brent_find_minima. So I'd have to have a lambda function anyway to call the private member function, or do a lot of friending.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I see. And C++ (even C++11) doesn't support defining a function inside a function.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think it does support that, actually, at least in some limited sense, but it's regarded as strongly deprecated in favor of lambdas.

OutlierRejectionControl const & ctrl
) {
if (!_keys.valid.isValid()) {
throw LSST_EXCEPT(
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I already admitted to not knowing what _keys.valid is (or _keys) but this test really surprises me. The very fact that the _keys.valid.isValid() is true tells you that you are trying to reject outliers? Why? A comment here might be helpful.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

_keys.valid is a Key for the field that's true if an object is outlier-rejected and false if it isn't. But if we construct with fromGrid, there's no uncertainty and hence no outlier rejection, and that field never gets added to the schema.

(I've added that as a comment).

if (!sipForward.getPixelOrigin().asEigen().isApprox(sipReverse.getPixelOrigin().asEigen())) {
throw LSST_EXCEPT(
pex::exceptions::InvalidParameterError,
"SIP forward and reverse transforms have inconsistent CRPIX"
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Print the values?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Done.

Copy link
Contributor

@r-owen r-owen left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Done. This is a really nice set of tests.

The main thing I ask is to make sure your python code is compatible with both python 2 and python 3 and that it flake8 reports it as clean. I believe meas_astrom has been updated to py3 compatibility, so you can test your new code with Jenkins py3 by building product lsst_py3 and skipping the demo.

#

import unittest
import numpy
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please change to import numpy as np

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Fine.

class TransformTestMixin(object):

def checkTransformEquivalence(self, a, b, rtol=1E-8):
for i in xrange(10):
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please use range and from builtins import range for python 3 compatibility. More generally, I suggest a futurize stage 2 pass on this code to make sure it all should be python 3 compatible.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Fixed. I'll run futurize if I can get it working quickly (haven't tried it before), and I'll certainly run the py3 Jenkins before merging to make sure I haven't missed anything.



class TransformTestMixin(object):

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

To improve the readability I suggest adding a makeRandom method that raises NotImlpementedError or using abc to make this an abstract base class, and documenting makeRandom (the API and the need to override it).

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Done.

self.assertClose(p.getXCoeffs(), xc, atol=0, rtol=0)
self.assertClose(p.getYCoeffs(), yc, atol=0, rtol=0)
# Test that the coefficients are not a view.
xc[0, 0] = 50.0
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It is not self-evident that the existing value cannot be 50. Please consider incrementing the existing value (after saving the initial value to use in the comparison test).

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Done.

composed1 = lsst.meas.astrom.compose(poly, affine)
composed2 = lsst.meas.astrom.compose(affine, poly)
self.checkTransformEquivalence(composed1, lambda p: poly(affine(p)))
self.checkTransformEquivalence(composed2, lambda p: affine(poly(p)))
Copy link
Contributor

@r-owen r-owen Oct 7, 2016

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The test doesn't quite seem to do what is advertised: check that the result of composing gives the expected result. I suggest you define a python transform that applies affine and then poly (or vice-versa) and check equivalence (similar to what ScaledPolynomialTransformTestCase .testConstruction does below)

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think you must have misread this - it's already doing exactly what you've asked; the "python transform that applies affine and then poly" is what the first of these lambdas is.

def t2(p):
sky = lsst.afw.coord.IcrsCoord(crval + lsst.afw.geom.Extent2D(p), lsst.afw.geom.degrees)
return wcs.skyToPixel(sky)
self.checkTransformEquivalence(t1, t2)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is it worth also checking fwd in the same way, with another pair of functions?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That's actually a lot harder - it requires actually inverting the polynomial part of rev (note that fwd isn't actually the inverse of rev right now). I do essentially that down in testFromGrid.

poly = makeRandomPolynomialTransform(order, sip=False)
return SipReverseTransform(origin, cd, poly)


Copy link
Contributor

@r-owen r-owen Oct 7, 2016

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is a very nice set of tests. However, I feel as if too much is contained in one file. I suggest you move the common code to a library module, such as lsst.meas.astrom.testUtils, including TransformTestMixin. Then divide the remaining subclasses of TransformTestMixin into individual test files with more descriptive filenames.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Also, this test has a number of minor flake8 linter errors. Please run flake8 on this and all the other new python code you added and touched.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

While it is a long file, I'd prefer to keep these together - I don't think there's anything intrinsically bad about having a few hundred lines of code in one file, and the test utility code is just limited to these tests rather than other stuff in meas_astrom. And, as we've discussed, I expect all of these classes to eventually be replaced by new AST-like things, and if these tests remain valuable then, I expect they'll be moved between packages and it will all be easier to deal with if it's all in one file.

I will run flake8.


class TransformTestMixin(object):

def checkTransformEquivalence(self, a, b, rtol=1E-8):
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please rename to assert..., perhaps assertTransformsNearlyEqual.

Please also consider moving this to afw.geom.testUtils, such that it becomes a standard part of lsst.utils.tests.TestCase with a name that reflects the fact that the transforms it tests are 2d (e.g. assert2DTransformsNearlyEqual) and with a user-specified domain.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I've renamed it as you've suggested.

I don't think I'm going to move it somewhere more general right now, though, because I think the domain/range questions are very important for testing transforms more generally, and I don't want to put in the work right now to make this utility function handle more general domains and ranges.


def testFromMatches(self):
# Setup artifical matches that correspond to a known (random) PolynomialTransform.
order = 3
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I suggest moving this code into a new method, e.g. makeMatchList, to make the test easier to read

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'd rather not; it'd only be used in one place, and I think the comments are sufficient to make it clear where in the code we're doing prep work rather than doing the tests. And I don't want to comb through this function to find which variables I set in the top half I'll need to pass out of the function unless there's a very good reason to split them up.

def t2(p):
sky = lsst.afw.coord.IcrsCoord(crval + lsst.afw.geom.Extent2D(p), lsst.afw.geom.degrees)
return wcs.skyToPixel(sky)
self.assertTransformsNearlyEqual(t1, t2)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I had asked if it was also worth copying this test for the reverse direction (using different words) and you replied "That's actually a lot harder - it requires actually inverting the polynomial part of rev (note that fwd isn't actually the inverse of rev right now). I do essentially that down in testFromGrid."

For the record, all I had in mind was literally testing rev the same way you test fwd, not testing them against each other (I realize they are not inverses of each other, since your documentation was clear on that point). I don't know if it would actually increase test coverage.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ah, now I understand. I think that's a good idea, and I started doing it - but to do it it looks like I'd need a Wcs method that doesn't exist: an inverse of skyTointermediateWorldCoord. (And now that I've said that, I can vaguely recall coming to the same conclusion back when I originally wrote this test).

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That makes sense. May I suggest you add a small comment to that effect?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Done.

* scatter estimation.
*
* In either case, the fitter creates affine transforms that map the input
* and output data points onto to [-1, 1]. In then fits a polynomial
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Typo: This should be "It then fits a polynomial...".

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Done.

@TallJimbo TallJimbo merged commit b143333 into master Oct 26, 2016
@ktlim ktlim deleted the tickets/DM-3549 branch August 25, 2018 06:15
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
4 participants