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-9751: Verify the performance of new matchPessimisticB code on selected test fields #72

Merged
merged 2 commits into from Jun 6, 2017

Conversation

morriscb
Copy link
Contributor

Several stability improvements also made to matcher during this verification.

matchPessimisticB.py: Reordered cut calculations for clarity. Changed first matcher iteration to only use optimistic mode only after a match has been found on a previous match/fit iteration. Changed cut on match distance to be a 2 sigma clip instead of a hard cut.
pessimistic_pattern_matcher_b_3D.py: Added least squares fit to shift and rotation matrix in intermediate verify. removed distance cut from final verify. This cut is now done in matchPessimisticB.

Removed max_dist from _match_sources method.

Linter pass and editted sigma clip start.

@@ -377,6 +365,18 @@ def _doMatch(self, refCat, sourceCat, wcs, refFluxField, numUsableSources,
src_array[src_idx, :] = \
self._latlong_flux_to_xyz_mag(theta, phi, flux)

# Set configurable defaults when we encounter None type or set
# state based on pervious run of matchandfit.
Copy link
Contributor

Choose a reason for hiding this comment

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

What is "matchandfit"? Perhaps you mean a different method, or if it's a method on a different class, then please provide the class name (e.g. Foo.matchandfit)

pervious -> previous

maxShiftArcseconds = (self.config.maxOffsetPix *
wcs.pixelScale().asArcseconds())
else:
# We don't want to clamp down to hard on the allowed shift so
Copy link
Contributor

Choose a reason for hiding this comment

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

to hard -> too hard

if ref_array.shape[0] < self.config.numBrightStars + \
self.config.numPointsForShapeAttempt or \
src_array.shape[0] < self.config.numBrightStars + \
self.config.numPointsForShapeAttempt:
Copy link
Contributor

Choose a reason for hiding this comment

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

I find this rather hard to parse. What is the sum numBrightStars + numPointsForShapeAttempt? If you can find a short name that describes it, I suggest you use that as a variable name, assign the sum to it, and test against that (and similarly for the second sum). Would numRefStarsRequired suit?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Should I make this intermediate variable or would you recommend putting it in the config?

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 an intermediate variable since it is a simple sum of two config variables. The only reason to put it in the config would be if the user might want to use some other value, but unless that's likely it just adds clutter.

@@ -472,23 +476,82 @@ def _doMatch(self, refCat, sourceCat, wcs, refFluxField, numUsableSources,
match_tolerance.lastMatchedPattern = matcher_struct.pattern_idx
break

# The matcher returns all the nearst neighbors between the reference
# source catalogs. For the currrent astrometry solver we need to clip
# these down to give the best solution.
Copy link
Contributor

Choose a reason for hiding this comment

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

Clip these down to what? Could you expand on this a bit?

Also nearst -> nearest

# The matcher returns all the nearst neighbors between the reference
# source catalogs. For the currrent astrometry solver we need to clip
# these down to give the best solution.
distances = np.degrees(matcher_struct.distances) * 3600
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 np.rad2deg instead because np.degrees sounds too close to afwGeom's degrees, which produces an Angle when multiplied by a float.

Also I suggest distances_arcsec as the variable name. Normally we express angles as afwGeom Angle and don't worry about units, but when we stray from that it's safest to be specific about units. Similarly matcher_struct.distances should be distances_rad.

Parameters
----------
flattened_rot_matrix : float array
A flattened array representing a 3x3 rotation matrix.
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 you use a normal 2D array, instead. It makes the API easier to understand, and the matrix needs to be unflattened at some point -- might as well do it before calling this.

if you insist on the data being a 1-D array then you must explain the order in which the elements are stored.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

The matrix needs to be flattened for the fitter that is used later in the code (scipy.optimize.leastsq) which only works with flattened arrays as input (I checked). I added an extra line of explanation to the doc string.

"""Compute the squared differences for least squares fitting.

Given a flattened rotation matrix, a N point pattern from the source
catalog and the reference pattern the sources match to, compute return
Copy link
Contributor

Choose a reason for hiding this comment

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

"compute return" -> "compute"

float array
Array of differences between the points of the source pattern
rotated into the reference frame and the reference pattern rotated
into the source frame.
Copy link
Contributor

Choose a reason for hiding this comment

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

I'm afraid I don't know what this means. What is the "difference between the points..."?

Also, why do you rotate the source into the reference frame and also rotate the reference pattern into the source frame. Naively I would expect only one rotation would be required to compare patterns of points in the two frames.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Added a extra set of explanation to this. The reason the second rotation exist is to test that the fitted rotation is a valid rotation that is np.dot(matrix.transpose(), matrix) = identity. I could equivalently test the above holds instead of doing the inverse rotation on the reference pattern. This may be clearer. This would replace line 45 with diff_ident = np.dot(rot_matrix.transpose(), rot_matrix) - np.identity(3).

Copy link
Contributor

@r-owen r-owen Jun 1, 2017

Choose a reason for hiding this comment

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

If it would be clearer (and at least as robust) then it sounds like a good change

# Don't fill the array column wise if we are on the last as
# these pairs have already been filled by pervious
# these pairs have already been filled by pervious
Copy link
Contributor

Choose a reason for hiding this comment

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

pervious -> previous

@@ -157,7 +152,7 @@ def plot(catalog, symbol):

# Current skip refObj 15 as it matches to source 12.
# All other matches are correct.
Copy link
Contributor

Choose a reason for hiding this comment

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

It appears this comment is outdated, based on your change below

morriscb added a commit that referenced this pull request Jun 1, 2017
Incorportated comments made on pull request #72.
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.

A few minor requests, but basically good to go

match.second = sourceCat[match_id_pair[0]]
# We compute the true distance along and sphere instead
# and store it as an afw.geom.Angle object. The previous
# distances we used were aproximate.
Copy link
Contributor

Choose a reason for hiding this comment

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

Unfortunately you can't store the distance as an Angle because ReferenceMath.distance is a C++ double. The angle is almost certainly being silently cast to radians; please make the cast explicit (and update the code comment), e.g.:

match.distance = match.first.getCoord().angularSeparation(match.second.getCoord()).asRadians()

----------
flattened_rot_matrix : float array
A flattened array representing a 3x3 rotation matrix. The array is
flattened to comply with the API of scipy.optimize.leastsq
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 the oder of the data -- how do the 9 values in the array represent the 3x3 values in the matrix? Is it elements [0,0], [0,1], [0,2], [1,0], ... or the transpose or ...?

@@ -227,7 +271,7 @@ def match(self, source_array, n_check, n_match, n_agree,

# Initialize output struct.
output_match_struct = pipeBase.Struct(
matches=[],
match_ids=[],
Copy link
Contributor

Choose a reason for hiding this comment

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

Thanks!

bool
Return true if all of the points in our source pattern are within
max_dist_rad of their matched reference objects.
float array
Copy link
Contributor

Choose a reason for hiding this comment

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

float array or None

if ref_array.shape[0] < self.config.numBrightStars + \
self.config.numPointsForShapeAttempt or \
src_array.shape[0] < self.config.numBrightStars + \
self.config.numPointsForShapeAttempt:
Copy link
Contributor

Choose a reason for hiding this comment

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

Thanks for defining minRefAndSrcObjectsForConsensus but then please use it in the test. (flake8 should have warned you that the variable is unused?) Also, consider shortening the name to minObjectsForConsensus.

Verified matcher on visits from HSC, HITS, and CFHTLS.

Several stability improvements also made to matcher during this verification.

matchPessimisticB.py: Reordered cut calculations for clarity. Changed first matcher iteration to only use optimistic mode only after a match has been found on a previous match/fit iteration. Changed cut on match distance to be a 2 sigma clip instead of a hard cut.
pessimistic_pattern_matcher_b_3D.py: Added least squares fit to shift and rotation matrix in intermediate verify. removed distance cut from final verify. This cut is now done in matchPessimisticB.

Removed max_dist from _match_sources method.

Linter pass and editted sigma clip start.

Comment clarity and some changes to variable names.

Incorportated comments made on pull request #72.

Cleared up comments and linted code.

Fixed comment wording and added missing test to matchPessimistcB.

matchPessimisticB.py: Added maxMatchDistArcsec to final cut on output
matches.
@morriscb morriscb merged commit 8f2c18f into master Jun 6, 2017
@ktlim ktlim deleted the tickets/DM-9751 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
2 participants