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-23786: Add bright star subtraction task #555

Merged
merged 3 commits into from Mar 31, 2023
Merged

Conversation

MorganSchmitz
Copy link
Contributor

@MorganSchmitz MorganSchmitz commented Jul 26, 2021

Copy link
Member

@arunkannawadi arunkannawadi left a comment

Choose a reason for hiding this comment

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

I've left a ton of comments, but almost all of them are trivial. My major complaint is that there is no test code to run this task. I noticed at least one place where running this task in accordance with the documentation would have still crashed it, so some unit testing is necessary before merging this.

dimensions=("band",),
)
skyCorr = cT.Input(
doc="Input Sky Correction to be subtracted from the calexp if doApplySkyCorr=True",
Copy link
Member

Choose a reason for hiding this comment

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

doApplySkyCorr should be within `` so it can link to the appropriate variable in the docs

Copy link
Contributor

Choose a reason for hiding this comment

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

Added - thanks!

)
inputExtendedPsf = cT.Input(
doc="Extended PSF model.",
name="extended_psf",
Copy link
Member

Choose a reason for hiding this comment

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

Has it been already defined elsewhere that the name of this input will be extended_psf? While it is within the rules, having some variables in camelCase and some in snake_case is going to get confusing.

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 this nomenclature has already been established in prior tickets. Future tickets should aim to standardize all the names associated with bright star subtraction.


def __init__(self, *, config=None):
super().__init__(config=config)
if not config.doApplySkyCorr:
Copy link
Member

Choose a reason for hiding this comment

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

I am not 100% sure on when the lines in Connections get triggered. Is there a test case that checks that both values of doApplySkyCorr do what we want it to do?

Copy link
Contributor

Choose a reason for hiding this comment

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

A future ticket which adds unit testing capability should run this with doApplySkyCorr set to both False and True (i.e., two separate tests).

"bilinear": "bilinear interpolation",
"lanczos3": "Lanczos kernel of order 3",
"lanczos4": "Lanczos kernel of order 4",
"lanczos5": "Lanczos kernel of order 5",
Copy link
Member

Choose a reason for hiding this comment

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

It might be preferable to allow up tolanczos7 as well. In general, the default value being one of the extremes can become inconvenient if we'd like to check the behavior around the default.

Copy link
Contributor

Choose a reason for hiding this comment

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

See comment above (added 6 and 7).

_DefaultName = "subtractBrightStars"

def __init__(self, initInputs=None, *args, **kwargs):
pipeBase.CmdLineTask.__init__(self, *args, **kwargs)
Copy link
Member

Choose a reason for hiding this comment

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

Can we not user super here?

Copy link
Contributor

Choose a reason for hiding this comment

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

Added super(), and updated the inheritance to also include pipeBase.PipelineTask to make this gen3 compatible.

model : `lsst.afw.image.MaskedImageF`
The extended PSF model, shifted (and potentially warped) to match
the bright star's positioning.
star : `lsst.meas.algorithms.brightStarStamps.BrightStarStamp`
Copy link
Member

Choose a reason for hiding this comment

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

Consider using ~lsst.meas.algorithms.brightStarStamps.BrightStarStamp (start with tilde within single backticks). I couldn't find the link in the style guide, but it keeps it short and still points to the right documentation.

Copy link
Contributor

Choose a reason for hiding this comment

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

Added - and on relevant line above - and throughout the entire file.

# "unnormalize" star
starIm *= star.annularFlux
# make sure both components use the star mask
XY = starIm.clone()
Copy link
Member

Choose a reason for hiding this comment

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

Variables cannot being with upper case letters, let alone being in all caps unless they are global variables: https://developer.lsst.io/python/style.html#naming-conventions. I also wouldn't rename this to xy because that still doesn't seem very indicative. I was expecting these to be multipled with spatial coordinates, probably to compute an inner product, but that is not what is going on in the code below. Perhaps something like cross and auto correlations?

Copy link
Contributor

Choose a reason for hiding this comment

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

Following discussion, agreed to change these to lower case (xx and xy, and their Sum equivalents), and added an in-line comment explaining these variables.

Copy link
Member

Choose a reason for hiding this comment

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

Or, use xx and xy variables but say the estimator of the scalingFactor (f) that minimizes (Y-fX)^2 is E[XY]/E[XX].

starIm = afwMath.rotateImageBy90(starIm, nb90Rots)
# "unnormalize" star
starIm *= star.annularFlux
# make sure both components use the star mask
Copy link
Member

Choose a reason for hiding this comment

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

More minor quibble - Comments must be full sentences, with starting letter capitalized and ending in a period. This applies to other comments in the file as well.

Copy link
Contributor

Choose a reason for hiding this comment

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

All comments in the file modified as suggested. This particular comment (# make sure both components use the star mask) removed, as it seemingly does not apply to this code block.

inputExtendedPsf : `lsst.pipe.tasks.extended_psf.ExtendedPsf`
Extended PSF model, produced by
`lsst.pipe.tasks.extended_psf.MeasureExtendedPsfTask`.
skyCorr : `lsst.afw.math.backgroundList.BackgroundList` or ``None``,
Copy link
Member

Choose a reason for hiding this comment

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

None should be within single backtics.

Copy link
Contributor

Choose a reason for hiding this comment

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

Modified.

to its corresponding input bright star.
"""
inputExpBBox = inputExposure.getBBox()
if self.config.doApplySkyCorr:
Copy link
Member

Choose a reason for hiding this comment

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

Do we have to check if skyCorr is not None here (default), or can applySkyCorr handle that case?

Copy link
Contributor

Choose a reason for hiding this comment

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

Modified run method doc-string, and added additional skyCorr check for None type.


class SubtractBrightStarsConnections(pipeBase.PipelineTaskConnections,
dimensions=("instrument", "visit", "detector"),
defaultTemplates={"outputExposureName": "subtracted",
Copy link
Contributor

Choose a reason for hiding this comment

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

Thinking long-term about this, is 'subtracted' a good descriptor for these outputs? I.e., many subtractions already take place across the pipeline (and perhaps future additional subtractive data products will also be defined), so perhaps this needs to be more descriptive? bright_star_subtracted_calexp? Although, now that I write it out, this seems too long... I'm not fully convinced either way on this, so I wanted to ping it back to you to hear your thoughts.

Copy link
Member

Choose a reason for hiding this comment

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

How about brightStar_subtracted_calexp?

Copy link
Contributor

Choose a reason for hiding this comment

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

Yes, that sounds good to me - adopted.

defaultTemplates={"outputExposureName": "subtracted",
"outputBackgroundName": "brightStars"}):
inputExposure = cT.Input(
doc="Input exposure from which to subtract bright star stamps",
Copy link
Contributor

Choose a reason for hiding this comment

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

(very very) minor comment - some of these doc-strings terminate with a period '.', whereas some don't.

Copy link
Contributor

Choose a reason for hiding this comment

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

Added.

doc="Extended PSF model.",
name="extended_psf",
storageClass="ExtendedPsf",
dimensions=("band",),
Copy link
Contributor

Choose a reason for hiding this comment

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

(another) (very very) minor comment - also, some tuples and/or connection blocks have terminating commas, and some don't.

Copy link
Contributor

Choose a reason for hiding this comment

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

Added trailing commas.

doc="Magnitude limit, in Gaia G; all stars brighter than this value will be subtracted",
default=18
)
warpingKernelName = pexConfig.ChoiceField(
Copy link
Contributor

Choose a reason for hiding this comment

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

Arun's comment from the other day was to also add lanczos7 here too - I'm not sure how difficult/different/desired this is in practice?

Copy link
Contributor

Choose a reason for hiding this comment

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

Added 6 and 7 as allowed options. Future unit testing capability should test that lanczos6 and lanczos7 work as expected.

doc="Mask planes that, if set, lead to associated pixels not being included in the computation of "
"the scaling factor. Ignored if scalingType is annularFlux, as the stamps are expected to "
"already be normalized.",
default=('BAD', 'CR', 'CROSSTALK', 'EDGE', 'NO_DATA', 'SAT', 'SUSPECT', 'UNMASKEDNAN')
Copy link
Contributor

Choose a reason for hiding this comment

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

Following our discussion yesterday, perhaps a note should be left here (either in the doc-string itself, or otherwise in a code comment block) reminding that secondary DETECTED sources (i.e., all sources other than the primary source) will be added to the BAD mask plane, and therefore, removing BAD from this list will also necessarily remove secondary source masking too.

Copy link
Contributor

Choose a reason for hiding this comment

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

Added text to doc string and and in-line comment with more detail.

Comment on lines 146 to 149
This task uses both an extended PSF model produced by
`lsst.pipe.tasks.extended_psf.MeasureExtendedPsfTask` and a set of bright
star stamps produced by
`lsst.pipe.tasks.processBrightStars.ProcessBrightStarsTask`.
Copy link
Contributor

Choose a reason for hiding this comment

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

These two might make more sense being swapped around, to indicate that ProcessBrightStarsTask must be run first, and then MeasureExtendedPsfTask subsequently run on the output.

Copy link
Contributor

Choose a reason for hiding this comment

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

Swapped.

Comment on lines 266 to 267
An Exposure containing every bright star's profile; its image can
then be subtracted from the input exposure.
Copy link
Contributor

Choose a reason for hiding this comment

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

Would this be more accurate reading something like An Exposure containing a scaled bright star model fit to every bright star profile; ...? (please re-word as necessary - I just wanted to communicate that the exposure doesn't contain the bright star profile itself, rather, some kind of model fit to said profile).

Copy link
Contributor

Choose a reason for hiding this comment

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

Modified as suggested.

# make a copy of the input model
model = inputExtendedPsf(dataId["detector"]).clone()
modelStampSize = model.getDimensions()
inv90Rots = 4 - inputBrightStarStamps.nb90Rots
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 sure this is obvious, but I can't see it right now - why is this 4 here?!

Copy link
Contributor

Choose a reason for hiding this comment

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

The 4 minus nb90Rots turns this into an inverse rotation.

Copy link
Member

Choose a reason for hiding this comment

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

I'm not sure if there's any check to ensure nb90Rots is an integer between 0-3. So this might be better off as

inv90Rots = 4 - inputBrightStarStamps.nb90Rots%4

inv90Rots = 4 - inputBrightStarStamps.nb90Rots
model = afwMath.rotateImageBy90(model, inv90Rots)
warpCont = afwMath.WarpingControl(self.config.warpingKernelName)
invImages = []
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 add a comment line above this explaining the purpose of the following for loop, e.g.,
# loop over bright stars, computing the inverse transformed and scaled postage stamp for each

Copy link
Contributor

Choose a reason for hiding this comment

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

Added.

Copy link
Member

@arunkannawadi arunkannawadi left a comment

Choose a reason for hiding this comment

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

I want to give this another pass, so not marking it as Approved for now. Happy to work on unit tests for this.

class SubtractBrightStarsConfig(PipelineTaskConfig, pipelineConnections=SubtractBrightStarsConnections):
"""Configuration parameters for SubtractBrightStarsTask"""

doWriteSubtractor = Field(
Copy link
Member

Choose a reason for hiding this comment

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

Either in this ticket, or in the refactoring ticket that should follow soon, it'd be good to make it as Field[bool] so type annotations can work.

Copy link
Contributor

Choose a reason for hiding this comment

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

Type annotations added to all Config options (including ChoiceField, which I think works now?)

default=("BAD", "CR", "CROSSTALK", "EDGE", "NO_DATA", "SAT", "SUSPECT", "UNMASKEDNAN"),
)
doApplySkyCorr = Field(
dtype=bool, doc="Apply full focal plane sky correction before extracting stars?", default=True
Copy link
Member

Choose a reason for hiding this comment

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

Is this a black thing? I'd like it personally if these were in their own individual lines. One way to make black wrap it would be to give a trailing comma, as is the case in previous Fields.

def applySkyCorr(self, calexp, skyCorr):
"""Apply correction to the sky background level.
Sky corrections can be generated with the 'skyCorrection.py'
executable in pipe_drivers. Because the sky model used by that
Copy link
Member

Choose a reason for hiding this comment

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

Is the above line still true?

Copy link
Contributor

Choose a reason for hiding this comment

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

Good catch! Modified.

xySum = makeStatistics(xy, self.statsFlag, self.statsControl).getValue()
xxSum = makeStatistics(xx, self.statsFlag, self.statsControl).getValue()
scalingFactor = xySum / xxSum if xxSum else 0
else:
Copy link
Member

Choose a reason for hiding this comment

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

Since scalingType is a ChoiceField, an invalid option is caught at config time and not in runtime. So this will never get triggered. You could keep it if you like, but this block is unnecessary IMO.

Copy link
Contributor

Choose a reason for hiding this comment

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

Removed. I'd like to replace this with a match during a refactor, as I think that's more readable, but leaving as-is for now sans the final else block.

# Compute the least squares scaling factor.
xySum = makeStatistics(xy, self.statsFlag, self.statsControl).getValue()
xxSum = makeStatistics(xx, self.statsFlag, self.statsControl).getValue()
scalingFactor = xySum / xxSum if xxSum else 0
Copy link
Member

Choose a reason for hiding this comment

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

Should this default to 1 if xxSum is None or NaN instead of 0? For the default value of inPlace=True, it looks like we might lose the model altogether if scaleModel was called and if makeStatistics behaved improperly.

Copy link
Contributor

Choose a reason for hiding this comment

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

Another good catch, changed to 1 👍

return scalingFactor

def runQuantum(self, butlerQC, inputRefs, outputRefs):
"""Run a single quantum of the task."""
Copy link
Member

Choose a reason for hiding this comment

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

You can remove this docstring and replace it with a comment "#Docstring inherited." If it didn't have a docstring, it'd inherit it automatically from the parent PipelineTask which is probably sufficient here.


Parameters
----------
inputExposure : `afwImage.exposure.exposure.ExposureF`
Copy link
Member

Choose a reason for hiding this comment

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

Suggested change
inputExposure : `afwImage.exposure.exposure.ExposureF`
inputExposure : `~lsst.afw.image.exposure.exposure.ExposureF`

Copy link
Member

Choose a reason for hiding this comment

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

Similarly in the Returns section as well.

# make a copy of the input model
model = inputExtendedPsf(dataId["detector"]).clone()
modelStampSize = model.getDimensions()
inv90Rots = 4 - inputBrightStarStamps.nb90Rots
Copy link
Member

Choose a reason for hiding this comment

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

I'm not sure if there's any check to ensure nb90Rots is an integer between 0-3. So this might be better off as

inv90Rots = 4 - inputBrightStarStamps.nb90Rots%4

@leeskelvin leeskelvin merged commit 3b994c1 into main Mar 31, 2023
1 check passed
@leeskelvin leeskelvin deleted the tickets/DM-23786 branch March 31, 2023 21:34
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

3 participants