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-39776: Add global astrometric fit task #30

Merged
merged 1 commit into from Feb 14, 2024
Merged

Conversation

cmsaunders
Copy link
Collaborator

No description provided.

@arunkannawadi
Copy link
Member

Would you be opposed to blackening the file? I'm thinking of adding black and isort formatting on DM-39243.

@cmsaunders
Copy link
Collaborator Author

Hmm, I don't have very strong feelings about using black, but would it mean the code has to be black compliant every time you make a commit? From my experience with analysis_tools, I do find that frustrating during development.

@timj
Copy link
Member

timj commented Oct 20, 2023

Once a repo is black-ified it has to stay that way and GitHub actions have to enforce it. This is a per-repo decision that should generally include discussions within the science pipelines team.

@arunkannawadi
Copy link
Member

I'll ask more broadly on Monday, but this is the only other file that has work in progress and one that I'm not working on.

@arunkannawadi
Copy link
Member

And no, I don't mean pre-commit. So you can still push code without it being black compliant, but you'll have to make it compliant before merging to main.

@@ -619,6 +625,8 @@ def _get_exposure_info(
data comes from the same instrument.
refEpoch : `float`
Epoch of the reference objects in MJD.
fieldRegions : `dict` of `lsst.sphgeom.ConvexPolygon`
Copy link
Contributor

Choose a reason for hiding this comment

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

`dict` [`KEY?`, `lsst.sphgeom.ConvexPolygon`]

https://developer.lsst.io/python/numpydoc.html#dictionary-types

@@ -619,6 +625,8 @@ def _get_exposure_info(
data comes from the same instrument.
refEpoch : `float`
Epoch of the reference objects in MJD.
fieldRegions : `dict` of `lsst.sphgeom.ConvexPolygon`
Dictionary of regions encompassing each group of input visits.
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 the key?

@@ -496,9 +500,10 @@ def run(
)

# Add the reference catalog to the associator
medianEpoch = astropy.time.Time(exposureInfo.medianEpoch, format="decimalyear").mjd
medianEpoch = astropy.time.Time(exposureInfo.medianEpoch, format='decimalyear').mjd
Copy link
Contributor

Choose a reason for hiding this comment

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

Why?

(decErr, "dec", 1),
(raPMErr, "pmra", 3),
(decPMErr, "pmdec", 4),
(parallaxErr, "parallax", 2)
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 the trailing comma?


__all__ = ["GbdesAstrometricFitConnections", "GbdesAstrometricFitConfig", "GbdesAstrometricFitTask"]
__all__ = ['GbdesAstrometricFitConnections', 'GbdesAstrometricFitConfig', 'GbdesAstrometricFitTask',
'GbdesGlobalAstrometricFitTask']
Copy link
Contributor

Choose a reason for hiding this comment

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

Why " -> '?

@@ -111,11 +115,11 @@ def _make_ref_covariance_matrix(
for i, pi in enumerate(positionParameters):
for j, pj in enumerate(positionParameters):
if i == j:
cov[:, k] = (refCat[f"{pi}Err"] ** 2 * inputUnit**2).to_value(units[j] * units[j])
cov[:, k] = ((refCat[f'{pi}Err'].value)**2 * inputUnit**2).to(units[j] * units[j]).value
Copy link
Contributor

Choose a reason for hiding this comment

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

Again, why " -> ' (here and below)?

And why does this one need the to_value -> .value change but the two below don't?

@cmsaunders cmsaunders force-pushed the tickets/DM-39776 branch 2 times, most recently from ae5176c to d14f708 Compare December 15, 2023 17:21
@@ -750,11 +751,7 @@ def _get_exposure_info(
ras.append(raDec.getRa().asRadians())
decs.append(raDec.getDec().asRadians())
if fieldRegions:
Copy link
Contributor

Choose a reason for hiding this comment

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

You check if fieldRegions: here, but if fieldRegions is None: below. Is there a reason for the inconsistency?

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

No, I just can never remember which one is preferred and random-walk between them.

extensionVisits.append(-1)
extensionDetectors.append(-1)
extensionType.append("REFERENCE")
# not used. There needs to be a separate catalog for each field (?)
Copy link
Contributor

Choose a reason for hiding this comment

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

Is this still a ??

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

No, removed it.

exposureInfo = pipeBase.Struct(
visits=visits, detectors=detectors, ras=ras, decs=decs, medianEpoch=medianEpoch
)

Copy link
Contributor

Choose a reason for hiding this comment

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

This seems to directly duplicate the lines above...

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

Weird--maybe a rebasing failure?

center : `lsst.geom.SpherePoint`, optional
Center of the circle in which to load reference objects.
radius : `lsst.sphgeom._sphgeom.Angle`, optional
Radius of the circle in which to load reference objects.
Copy link
Contributor

Choose a reason for hiding this comment

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

Note that center must also be set and is ignored if region is not None?

associations : `wcsfit.FoFClass`, optional
Object to which to add the catalog of reference objects.
center : `lsst.geom.SpherePoint`, optional
Center of the circle in which to load reference objects.
Copy link
Contributor

Choose a reason for hiding this comment

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

Note that radius must also be set? Also that this will be ignored if region is set??

WCS-fitting object.
sourceDict : `dict`
Dictionary containing the source centroids for each visit.
extensionInfo : `lsst.pipe.base.Struct`
Copy link
Contributor

Choose a reason for hiding this comment

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

----------
wcsf : `wcsfit.WCSFit`
WCS-fitting object.
sourceDict : `dict`
Copy link
Contributor

Choose a reason for hiding this comment

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

dict [?, ?]?

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 did

sourceDict : `dict` [`int`, [`int`, [`str`, `list` [`float`]]]]

I'm not sure if getting this in-depth is the best practice, but I think it's accurate.


self.log.info("Fit the WCSs")
# Set up a YAML-type string using the config variables and a sample
# visit
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 this for?

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

Added "which will be used to tell gbdes what form of model to fit."

self.log.info("Fit the WCSs")
# Set up a YAML-type string using the config variables and a sample
# visit
inputYAML = self.make_yaml(inputVisitSummaries[0])
Copy link
Contributor

Choose a reason for hiding this comment

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

camelCase?

)
outputCatalog = pipeBase.connectionTypes.Output(
doc=(
"Source table with stars used in fit, along with residuals in pixel coordinates and tangent "
Copy link
Contributor

Choose a reason for hiding this comment

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

"stars" -> "sources"?

record.set(raErrKey, 0.000001)
record.set(decErrKey, 0.000001)
record.set(pmraErrKey, 1e-10)
record.set(pmdecErrKey, 1e-10)
Copy link
Contributor

Choose a reason for hiding this comment

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

Why are these changing?

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 think I was having trouble with the fit being accurate enough. In the mean time, I have improved another aspect for the fitting though, so I should test this again to see if this is still necessary.

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

Maybe because of other changes since this PR started, this change doesn't seem to be necessary. I have to run jenkins to check macos too though.

starRAs = np.delete(starRAs, neighbors)
starDecs = np.delete(starDecs, neighbors)
cls.nStars = len(starRAs)
starIds = np.arange(cls.nStars)
Copy link
Contributor

Choose a reason for hiding this comment

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

This block is duplicated above. Worth making it it's own function?


for region in self.fieldRegions.values():
refCat, refCov = self.task._load_refcat(
self.refObjectLoader, self.extensionInfo, region=region, epoch=2015
Copy link
Contributor

Choose a reason for hiding this comment

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

Isn't epoch usually expressed in MJD (or do you really mean 1864? 😉)

# Set random seed
np.random.seed(1234)

# Fraction of simulated stars in the reference catalog and science
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 the meaning/purpose of "simulated" stars here?

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 moved this block later so it should make more sense.


@classmethod
def _make_isolatedStars(cls, allStarIds, allStarRAs, allStarDecs, trueWCSs, inScienceFraction):
# Take a subset of the simulated data
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 having trouble understanding what is meant by "simulated data" for these tests and how these are used vs. the real data you are also using (so it's hard to get a feel for how much of the task is getting tested). A few words of explanation would go a long way (for me at least!)

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 think these were notes to myself that I missed converting into docs. I added full documentation here, plus a short explanation in the class documentation for the parent class:

    This class tests `GbdesAstrometricFit` using real `visitSummaryTable`s
    from HSC RC2 processing, with simulated sources for those visits.

The only real data in the tests is the visitSummaries. The sources and true object positions are all simulated from some "true WCS", which is known to be realistic given the input information in the visitSummaries. The overarching test of the task, which happens in test_run(), checks that, given the inputSummaries and simulated source catalogs, you can recover the "trueWCS" to a certain accuracy.

@@ -721,6 +750,16 @@ def _get_exposure_info(
raDec = visitInfo.getBoresightRaDec()
ras.append(raDec.getRa().asRadians())
decs.append(raDec.getDec().asRadians())
if fieldRegions:
inField = [r for r, region in fieldRegions.items() if region.contains(raDec.getVector())]
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 not sure why breaking as soon as len(inField) > 1 doesn't count as the same sanity check, but ok.

@@ -1407,3 +1477,496 @@ def _compute_model_params(self, wcsf):
modelParams[key] = np.array(value)

return modelParams


class GbdesGlobalAstrometricFitConnections(
Copy link
Contributor

Choose a reason for hiding this comment

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

Cool...thanks for checking with the masses!


refConfig = LoadReferenceObjectsConfig()
refConfig.anyFilterMapsToThis = "phot_g_mean"
refConfig.requireProperMotion = True
Copy link
Contributor

Choose a reason for hiding this comment

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

But what if another survey wanted/needed to use something else? One use case I'm thinking of is LSSTCam-imSim...

Parameters
----------
inputVisitSummaries : `list` of `lsst.afw.table.ExposureCatalog`
List of catalogs with per-detector summary information.
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 it may have changed in the past...this is "current" truth 😉


for i in range(clusters.n_clusters_):
fieldInd = clusters.labels_ == i
fieldDetectors = sum([allDetectorCorners[f] for f, fInd in enumerate(fieldInd) if fInd], [])
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 the sum part that is confusing me...but if it's doing what you want (and now describe), ok!

@cmsaunders cmsaunders merged commit 6ef0a89 into main Feb 14, 2024
3 checks passed
@cmsaunders cmsaunders deleted the tickets/DM-39776 branch February 14, 2024 15:13
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

4 participants