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

DM-13272 write tests for astrometryModel Wcs output #69

Merged
merged 22 commits into from Apr 17, 2018

Conversation

parejkoj
Copy link
Collaborator

In addition to new tests of the AstrometryModel wcs generation, this also cleans up several file and class names and a variety of other minor code issues (e.g. pointers, types, error messages).

@parejkoj parejkoj force-pushed the tickets/DM-13272 branch 6 times, most recently from 783a8bf to 157d560 Compare February 15, 2018 21:25
Copy link
Contributor

@PaulPrice PaulPrice left a comment

Choose a reason for hiding this comment

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

Not yet done; got through "Make gtransfoCompose bare pointers into references".

@@ -603,14 +603,17 @@ def _fit_astrometry(self, associations, dataName=None):
sky_to_tan_projection = lsst.jointcal.OneTPPerVisitHandler(associations.getCcdImageList())

if self.config.astrometryModel == "constrainedPoly":
model = lsst.jointcal.ConstrainedPolyModel(associations.getCcdImageList(),
sky_to_tan_projection, self.config.useInputWcs, 0,
Copy link
Contributor

Choose a reason for hiding this comment

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

The 0 disappeared. Maybe on purpose?

Copy link
Collaborator Author

@parejkoj parejkoj Mar 18, 2018

Choose a reason for hiding this comment

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

Yes, that argument, nNotFit, no longer applies to this constructor (see the commit message).

<< " not found in constrainedPolyModel mapping list.");
LOGLS_ERROR(_log, "CcdImage with ccd/visit "
<< ccdImage.getCcdId() << "/" << ccdImage.getVisit()
<< " not found in constrainedAstrometryModel mapping list.");
Copy link
Contributor

Choose a reason for hiding this comment

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

Strange indentation.

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

I run clang-format automatically on save for my C++, so I'm not responsible for any formatting. Take it up with clang-format

@@ -48,7 +48,14 @@ SimpleAstrometryModel::SimpleAstrometryModel(CcdImageList const &ccdImageList,
}
Copy link
Contributor

Choose a reason for hiding this comment

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

Suspect the summary message for this commit is too long.

while (unsigned(pol.getNpar()) > 2 * nObj) {
LOGLS_WARN(_log, "Reducing polynomial degree from "
<< pol.getDegree() << ", due to too few sources (" << nObj
<< " vs. " << pol.getNpar() << " parameters)");
Copy link
Contributor

Choose a reason for hiding this comment

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

Can you write this so it only logs a warning once rather than every time around the loop?

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

This is a situation I hope never happens, and when it does I'd rather it be verbose about it.

while (unsigned(pol.getNpar()) > 2 * nObj) {
LOGLS_WARN(_log, "Reducing polynomial degree from "
<< pol.getDegree() << ", due to too few sources (" << nObj
<< " vs. " << pol.getNpar() << " parameters)");
Copy link
Contributor

Choose a reason for hiding this comment

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

Strange indentation.

degree=self.degree2)
self._prepModels()

def test_getNpar(self):
Copy link
Contributor

Choose a reason for hiding this comment

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

Underscore.


def checkParams(ccdImage, model, chipDegree, visitDegree):
result = model.getNpar(ccdImage)
failMsg = "ccdImage: %s, with chipDegree %s and visitDegree %s"%(ccdImage.getName(),
Copy link
Contributor

Choose a reason for hiding this comment

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

Space either side of % operator.


std::unique_ptr<Gtransfo> gtransfoCompose(const Gtransfo *left, const Gtransfo *right);
/**
* Returns a pointer to a composition of gtransfos, representing `left(right())`.
Copy link
Contributor

Choose a reason for hiding this comment

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

Allow me to again register my dissatisfaction with the name gtransfo...

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

I'm open to a new name (possibly AstrometryTransfo to match the newer PhotometryTransfo?) but that should happen on a different ticket, I think: it'd be a big change to manage.

Copy link
Contributor

Choose a reason for hiding this comment

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

Why would you want to call anything a "Transfo"????????????

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

As I said, I'm not renaming the object on this ticket.

@@ -181,8 +189,13 @@ class GtransfoIdentity : public Gtransfo {
// ClassDef(GtransfoIdentity,1)
};

//! Shorthand test to tell if a transfo belongs to the GtransfoIdentity class.
bool isIdentity(const Gtransfo *gtransfo);
Copy link
Contributor

Choose a reason for hiding this comment

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

Why remove this?

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

It wasn't used anywhere.

@@ -163,7 +171,7 @@ class GtransfoIdentity : public Gtransfo {
"GtransfoIdentity is the identity transformation: it cannot be fit to anything.");
}

std::unique_ptr<Gtransfo> reduceCompo(const Gtransfo *right) const { return right->clone(); }
std::unique_ptr<Gtransfo> reduceCompo(Gtransfo const &right) const { return right.clone(); }
Copy link
Contributor

Choose a reason for hiding this comment

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

So long as you're changing the API, can you change the name to something a bit more standard?

Copy link
Contributor

Choose a reason for hiding this comment

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

And this method needs to be documented.

@parejkoj
Copy link
Collaborator Author

Thanks for starting on the review. While rebasing to master (with the new SkyWcs), I discovered a failure in the test I've added. Russell and I have decided to try to implement a AstrometryModel.makeSkyWcs() that bypasses SIP entirely (which is what I've wanted to do for a while). If we can get that to work, it means we can remove a big block of code involved in this review. Russell and I are both on vacation next week, so won't be able to make progress on that until after the 25th.

Please hold off on reviewing further, and I'll ping you once things are in a better place. I'll also try to deal with your above comments as part of that.

@parejkoj parejkoj force-pushed the tickets/DM-13272 branch 4 times, most recently from 8d1c9d9 to ff1d7a0 Compare March 6, 2018 22:04
@parejkoj parejkoj force-pushed the tickets/DM-13272 branch 8 times, most recently from 7f9d024 to d2a5a5e Compare March 19, 2018 23:35
@parejkoj parejkoj force-pushed the tickets/DM-13272 branch 2 times, most recently from c601265 to 80b7551 Compare March 27, 2018 20:47
Copy link
Contributor

@PaulPrice PaulPrice left a comment

Choose a reason for hiding this comment

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

Got through "Cleanup AstrometryModel interface to use findMapping and getHashKey"

@@ -20,6 +20,9 @@ class Gtransfo;
the top of simplepolymodel.h */
class AstrometryModel {
public:
/// Return the number of parameters in the mapping of CcdImage
unsigned getNpar(CcdImage const &ccdImage) const { return findMapping(ccdImage)->getNpar(); }
Copy link
Contributor

Choose a reason for hiding this comment

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

Gtransfo::getNpar returns an int, not an unsigned.

mappingMapType _mappings;
typedef std::map<CcdIdType, std::unique_ptr<SimpleGtransfoMapping>> chipMapType;
chipMapType _chipMap;
typedef std::map<VisitIdType, std::unique_ptr<SimpleGtransfoMapping>> visitMapType;
visitMapType _visitMap;
const ProjectionHandler *_sky2TP;
bool _fittingChips, _fittingVisits;

AstrometryMapping *findMapping(CcdImage const &ccdImage) const;
Copy link
Contributor

Choose a reason for hiding this comment

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

Docstring.

@@ -80,14 +80,16 @@ class ConstrainedAstrometryModel : public AstrometryModel {
std::shared_ptr<TanSipPix2RaDec> produceSipWcs(CcdImage const &ccdImage) const;

private:
typedef std::map<const CcdImage *, std::unique_ptr<TwoTransfoMapping>> mappingMapType;
typedef std::map<CcdImageKey, std::unique_ptr<TwoTransfoMapping>> mappingMapType;
Copy link
Contributor

Choose a reason for hiding this comment

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

Shouldn't this variable name have a leading underscore?

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

I've sidestepped the whole thing by removing the typedefs: it's probably cleaner without them anyway (and they're not used elsewhere).

Copy link
Contributor

Choose a reason for hiding this comment

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

Oh, silly me --- it's not a variable name. I'm told that using X = Y is now preferred to typedef Y X.

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

Yes, using X=Y is the new standard. But these particular typedefs are all gone now.

@@ -37,6 +37,11 @@ void declareAstrometryMapping(py::module &mod) {
py::class_<AstrometryMapping, std::shared_ptr<AstrometryMapping>> cls(mod, "AstrometryMapping");

cls.def("getNpar", &AstrometryMapping::getNpar);
cls.def("transformPosAndErrors", [](AstrometryMapping const &self, jointcal::FatPoint &inPos) {
Copy link
Contributor

Choose a reason for hiding this comment

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

Why would you name it transformPosAndErrors? What else would you transform?

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

You could transform the positions or the errors separately (e.g. gtransfo.apply()). This does both, while handling the error derivatives correctly.

Copy link
Contributor

Choose a reason for hiding this comment

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

But the function signature distinguishes it from the existing apply method. So it sounds like this should just be an overload of apply, and you could have a separate method called applyErrors if you need it.

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

apply() doesn't "apply" to both positions and errors, just positions, and is not lifted into the Mappings and Models.

If you're suggesting to overload apply with FatPoint to take over this job, I'd be worried about having to downcast a FatPoint-derived thing (e.g. MeasuredStar) into a Point before passing it to apply when you want to just transform the position. I think it's clearer having them be distinct names, unless we were to totally rework the gtransfo API. Which is not this ticket.

"ccdImageList"_a, "projectionHandler"_a, "initFromWcs"_a, "chipDegree"_a = 3,
"visitDegree"_a = 2);

cls.def("getChipTransfo", &ConstrainedAstrometryModel::getChipTransfo,
Copy link
Contributor

Choose a reason for hiding this comment

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

The multiplication of Transfos is really annoying me. I think you should bring the codebase up to spec now. Otherwise you're just increasing the work you're going to have to do in the future.


AstrometryMapping *ConstrainedAstrometryModel::findMapping(CcdImage const &ccdImage) const {
auto i = _mappings.find(ccdImage.getHashKey());
if (i == _mappings.end())
Copy link
Contributor

Choose a reason for hiding this comment

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

i->first?

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

I'm confused by this comment. What are you suggesting be changed to i->first? This is the standard way to check that a thing isn't in a std::map.

Copy link
Contributor

Choose a reason for hiding this comment

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

No, I was the one that was confused. You're quite right.

@PaulPrice
Copy link
Contributor

The commit "Implement DM-13800: fix one chip in constrainedAstromModel" stands in stark contrast to a lot of the other commits: it bundles a bunch of changes, whereas most other changes are in their own commits.

@PaulPrice
Copy link
Contributor

Commit "Move applyTransfo(frame) into gtransfo" is another that should be orthogonalised.

Copy link
Contributor

@PaulPrice PaulPrice left a comment

Choose a reason for hiding this comment

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

This is the balance of the comments.

@@ -69,8 +71,6 @@ void declareBaseStar(py::module &mod) {

// cls.def("getFlux", &BaseStar::getFlux);
cls.def_property_readonly("flux", (double (BaseStar::*)() const) & BaseStar::getFlux);

utils::python::addOutputOp(cls, "__str__");
Copy link
Contributor

Choose a reason for hiding this comment

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

Why is this being removed, when it's being added everywhere else?

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

It got lifted into the parent class Point, so it's not needed here.

@@ -248,6 +248,17 @@ class GtransfoPoly : public Gtransfo {
//! Constructs a "polynomial image" from an existing transfo, over a specified domain
GtransfoPoly(const Gtransfo *gtransfo, const Frame &frame, unsigned degree, unsigned nPoint = 1000);

/**
* @brief Constructs an polynomial approximation to an afw::geom::Transform2To2.
Copy link
Contributor

Choose a reason for hiding this comment

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

"an polynomial" --> "a polynomial"

@@ -248,6 +248,17 @@ class GtransfoPoly : public Gtransfo {
//! Constructs a "polynomial image" from an existing transfo, over a specified domain
GtransfoPoly(const Gtransfo *gtransfo, const Frame &frame, unsigned degree, unsigned nPoint = 1000);

/**
* @brief Constructs an polynomial approximation to an afw::geom::Transform2To2.
Copy link
Contributor

Choose a reason for hiding this comment

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

Please spell out TransformPoint2ToPoint2 so it's more obvious and easier to find.

*
* @param[in] transform The transform to be approximated.
* @param[in] domain The valid domain of the transform.
* @param[in] degree The polynomial degree to use when approximating.
Copy link
Contributor

Choose a reason for hiding this comment

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

I think we use order in various other places instead of degree, and I think you should adopt the same.

* @param[in] transform The transform to be approximated.
* @param[in] domain The valid domain of the transform.
* @param[in] degree The polynomial degree to use when approximating.
* @param[in] nSteps The number of sample points per axis (nSteps^2 total points).
Copy link
Contributor

Choose a reason for hiding this comment

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

In general, I don't think it's a good idea to line things up over multiple lines, because when you want to add something you invariably have to touch multiple lines. I believe this is reflected in our standards, though probably not specifically in relation to docstrings, but I think it would be good to follow the general rule.

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

It looks like we only have one example of multiple @params in the dev docs, but it's not lined up. Might be good to be explicit about this here: https://developer.lsst.io/cpp/api-docs.html#function-method-parameters

I was using DoxyDoxygen to help autofill C++ docstrings, which preferred to align them. I've switched to DoxyDoc, which doesn't.

*
* @return The transformed frame.
*/
Frame applyTransfoToFrame(Frame const &inputframe, Gtransfo const &gtransfo, WhichTransformed const which);
Copy link
Contributor

Choose a reason for hiding this comment

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

Why isn't this an apply method on Gtransfo?

If WhichTransformed only has two possible values, why not make it a bool?


src = dataRef.get("src", flags=lsst.afw.table.SOURCE_IO_NO_FOOTPRINTS, immediate=True)
goodSrc = sourceSelector.selectSources(src)
# Need memory contiguity to do vector-like things on the sourceCat.
Copy link
Contributor

Choose a reason for hiding this comment

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

DM-12207 in afw. You reviewed it. (:

def testMakeSkyWcsModel2(self):
self._testMakeSkyWcs(self.model2, self.fitter2, self.inverseMaxDiff2)

def _testMakeSkyWcs(self, model, fitter, inverseMaxDiff):
Copy link
Contributor

Choose a reason for hiding this comment

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

Looks like this is a test that has been disabled by the leading underscore. It should be renamed.

"""Test converting this object to a FITS-SIP gtransfo."""

sipWcs = model.produceSipWcs(ccdImage)
sky2TP = model.getSky2TP(ccdImage)
Copy link
Contributor

Choose a reason for hiding this comment

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

You should fix here and now what you can, like your choice of variable names (e.g., sky2TP --> skyToTp).

points = []
for x in xx:
for y in yy:
points.append(lsst.afw.geom.Point2D(x, y))
Copy link
Contributor

Choose a reason for hiding this comment

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

points = [lsst.afw.geom.Point2D(*xy) for xy in itertools.product(xx, yy)]

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

Oh, I totally forgot about itertools.product. Nice catch!

Now that there are Simple and Constrained PhotometryModels, it's good to have
the same naming conventions for AstrometryModels too.

Remove unused and undocumented `nNotFit` arg from ConstrainedAstro constructor.
@parejkoj
Copy link
Collaborator Author

parejkoj commented Apr 11, 2018

I believe I've responded to all your comments (either inline, or by just fixing the thing). The only remaining points of contention I see are:

  1. degree vs. order for the maximum monomial of the polynomial.
  2. Renaming/refactoring transformPosAndErrors.
  3. Refactoring the initialization minimize() calls in jointcal.py.

For 1., my informal poll on slack#dm-tea-time had a very slight advantage for degree. @r-owen said he may have been responsible for most of the existing uses of order in the stack. Per my argument above, I find order to be confusing because it can also imply the number of parameters in the model. I'll expand on my poll in tea-time and see if that leads anywhere.

For 2. I'd rather work that change into DM-4044, since that ticket may entail refactoring the inheritance of the various Point->Star relationships into composition, which would mean rethinking both apply and transformPosAndErrors.

For 3. when I tried to refactor and consolidate the fitter initialization between astrometry and photometry, I ended up down a rabbit hole because of the different choices of whatToFit between those. I marked it in DM-8046 as something to explore as part of that ticket, and I'd rather not try to rework it here.

@parejkoj
Copy link
Collaborator Author

Oh, and as to the comment about the two large commits:

  1. "Move applyTransfo": this one had some leftover statements in it that became no longer relevant after a rebase, so I've reworded it.
  2. "Implement DM-13800": all of the changes in this ticket were necessary to make it work. It might be possible to split it up into a couple of commits, but it would be difficult because they're intertwined or not meaningful when not together.

@TallJimbo
Copy link
Member

On order vs degree, I wouldn't have an ab initio preference between them for standard polynomials, but we do have a pretty strong precedent for using order in the stack, and I think that moves the needle in its favor (and does so more than going for external consistency with scipy/numpy by using degree). However, I don't think fixing that should be a high priority if it's in any way difficult, and I could be persuaded that changing the rest of the stack to use degree is a better choice (but also at very low priority).

@r-owen
Copy link
Contributor

r-owen commented Apr 16, 2018

I agree with everything @TallJimbo said, but have a stronger bias to stick with the term order in our stack because I fear that changing the stack to use degree instead of order will be much more trouble than it is worth. The term is widely used, including in many configs. Here are some examples I found: background.approxOrderX/Y, wcsFittter.order, and psfDeterminer.spatialOrder, measureApCorr.fitConfig.orderX/Y, models['DoubleShapelet'].wings.order.

@parejkoj
Copy link
Collaborator Author

I believe I've successfully changed all of jointcal's use of degree to order in this context. I'm waiting for a scons run to finish before I push.

parejkoj and others added 20 commits April 16, 2018 22:43
Make ConstrainedAstrometryModel use an unordered_map: almost all accesses are
random, so we don't need the orderedness of `map`. This makes it the same as
ConstrainedPhotometryModel.

Expand pybind11 interface of Point, FatPoint, AstrometryMapping, gtransfo,
CcdImage, AstrometryModels.
Rename reduceCompo -> composeAndReduce, and document.
makeSkyWcs() bypasses creation of a SIP and just stuffs all the gtransfos in a
model into an AST FrameDict.

BUGFIX! in inversePolyTransfo: the grid was not properly sampling the domain,
but was assuming it started at 0,0 and was recentering so it didn't include
the bounds of the box.

For inversePolyTransfo():
* More useful messages when there is a problem with the fit.
* Allow bailing when fitting if chi2 increases.
* Add TRACE logging for inversion steps.

Improve error reporting in getChipTransfo/getVisitTransfo

Expand gtransfo pybind11, adding astshim includes for toAstMap().
Support reading/writing polynomials to strings via python.
Increased stream output precision.
Fit the pixel->focal cameraGeom to initialize the chip transfos.
Add GtransfoPoly constructor to approximate afw.Transform2To2.
Remove initFromWCS option to constrainedAstromModel, since it does nothing now.

Better default polynomial degrees for constrained model.

First fitting step of constrained models is now `DistortionsVisit`, since we've
already got the chips in a reasonable state during initialization. This makes
the first step much better (gets closer to minimum-chi2) than it was.

Add toFit getter/setter for simpleAstroMap and use it in constrained model

Use shared instead of unique ptr for shared transfos; avoids the dangerous raw
pointers that I got to segfault.

Tweak chi2 and tighten RMS thresholds for cfht and decam constrained tests
due to new model, defaults, and DistortionsVisit step.
Restrict DECam to good ccds, correct expected test values to match. The DECam
testdata now produce quite reasonable fits: I should have restricted to those
ccds from the start.
Move AstroUtils::applyTransfo()->Gtransfo::frame() and cleanup cruft.
Make it, and the necessary enums visible to python.
Use it in ConstrainedAstrometryModel::toSkyWcs().
After a straw poll on slack and subsequent discussion, I've made jointcal
consistent with the rest of the stack by changing `degree` to `order` as the
term for the maximum monomial.
@parejkoj parejkoj merged commit ab4c137 into master Apr 17, 2018
@ktlim ktlim deleted the tickets/DM-13272 branch August 25, 2018 06:44
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
5 participants