LSSTApplications  19.0.0-14-gb0260a2+72efe9b372,20.0.0+7927753e06,20.0.0+8829bf0056,20.0.0+995114c5d2,20.0.0+b6f4b2abd1,20.0.0+bddc4f4cbe,20.0.0-1-g253301a+8829bf0056,20.0.0-1-g2b7511a+0d71a2d77f,20.0.0-1-g5b95a8c+7461dd0434,20.0.0-12-g321c96ea+23efe4bbff,20.0.0-16-gfab17e72e+fdf35455f6,20.0.0-2-g0070d88+ba3ffc8f0b,20.0.0-2-g4dae9ad+ee58a624b3,20.0.0-2-g61b8584+5d3db074ba,20.0.0-2-gb780d76+d529cf1a41,20.0.0-2-ged6426c+226a441f5f,20.0.0-2-gf072044+8829bf0056,20.0.0-2-gf1f7952+ee58a624b3,20.0.0-20-geae50cf+e37fec0aee,20.0.0-25-g3dcad98+544a109665,20.0.0-25-g5eafb0f+ee58a624b3,20.0.0-27-g64178ef+f1f297b00a,20.0.0-3-g4cc78c6+e0676b0dc8,20.0.0-3-g8f21e14+4fd2c12c9a,20.0.0-3-gbd60e8c+187b78b4b8,20.0.0-3-gbecbe05+48431fa087,20.0.0-38-ge4adf513+a12e1f8e37,20.0.0-4-g97dc21a+544a109665,20.0.0-4-gb4befbc+087873070b,20.0.0-4-gf910f65+5d3db074ba,20.0.0-5-gdfe0fee+199202a608,20.0.0-5-gfbfe500+d529cf1a41,20.0.0-6-g64f541c+d529cf1a41,20.0.0-6-g9a5b7a1+a1cd37312e,20.0.0-68-ga3f3dda+5fca18c6a4,20.0.0-9-g4aef684+e18322736b,w.2020.45
LSSTDataManagementBasePackage
jointcal.py
Go to the documentation of this file.
1 # This file is part of jointcal.
2 #
3 # Developed for the LSST Data Management System.
4 # This product includes software developed by the LSST Project
5 # (https://www.lsst.org).
6 # See the COPYRIGHT file at the top-level directory of this distribution
7 # for details of code ownership.
8 #
9 # This program is free software: you can redistribute it and/or modify
10 # it under the terms of the GNU General Public License as published by
11 # the Free Software Foundation, either version 3 of the License, or
12 # (at your option) any later version.
13 #
14 # This program is distributed in the hope that it will be useful,
15 # but WITHOUT ANY WARRANTY; without even the implied warranty of
16 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 # GNU General Public License for more details.
18 #
19 # You should have received a copy of the GNU General Public License
20 # along with this program. If not, see <https://www.gnu.org/licenses/>.
21 
22 import collections
23 import os
24 
25 import astropy.time
26 import numpy as np
27 import astropy.units as u
28 
29 import lsst.geom
30 import lsst.utils
31 import lsst.pex.config as pexConfig
32 import lsst.pipe.base as pipeBase
33 from lsst.afw.image import fluxErrFromABMagErr
34 import lsst.pex.exceptions as pexExceptions
35 import lsst.afw.table
36 import lsst.log
38 from lsst.pipe.tasks.colorterms import ColortermLibrary
39 from lsst.verify import Job, Measurement
40 
41 from lsst.meas.algorithms import LoadIndexedReferenceObjectsTask, ReferenceSourceSelectorTask
42 from lsst.meas.algorithms.sourceSelector import sourceSelectorRegistry
43 
44 from .dataIds import PerTractCcdDataIdContainer
45 
46 import lsst.jointcal
47 from lsst.jointcal import MinimizeResult
48 
49 __all__ = ["JointcalConfig", "JointcalRunner", "JointcalTask"]
50 
51 Photometry = collections.namedtuple('Photometry', ('fit', 'model'))
52 Astrometry = collections.namedtuple('Astrometry', ('fit', 'model', 'sky_to_tan_projection'))
53 
54 
55 # TODO: move this to MeasurementSet in lsst.verify per DM-12655.
56 def add_measurement(job, name, value):
57  meas = Measurement(job.metrics[name], value)
58  job.measurements.insert(meas)
59 
60 
61 class JointcalRunner(pipeBase.ButlerInitializedTaskRunner):
62  """Subclass of TaskRunner for jointcalTask
63 
64  jointcalTask.runDataRef() takes a number of arguments, one of which is a list of dataRefs
65  extracted from the command line (whereas most CmdLineTasks' runDataRef methods take
66  single dataRef, are are called repeatedly). This class transforms the processed
67  arguments generated by the ArgumentParser into the arguments expected by
68  Jointcal.runDataRef().
69 
70  See pipeBase.TaskRunner for more information.
71  """
72 
73  @staticmethod
74  def getTargetList(parsedCmd, **kwargs):
75  """
76  Return a list of tuples per tract, each containing (dataRefs, kwargs).
77 
78  Jointcal operates on lists of dataRefs simultaneously.
79  """
80  kwargs['profile_jointcal'] = parsedCmd.profile_jointcal
81  kwargs['butler'] = parsedCmd.butler
82 
83  # organize data IDs by tract
84  refListDict = {}
85  for ref in parsedCmd.id.refList:
86  refListDict.setdefault(ref.dataId["tract"], []).append(ref)
87  # we call runDataRef() once with each tract
88  result = [(refListDict[tract], kwargs) for tract in sorted(refListDict.keys())]
89  return result
90 
91  def __call__(self, args):
92  """
93  Parameters
94  ----------
95  args
96  Arguments for Task.runDataRef()
97 
98  Returns
99  -------
100  pipe.base.Struct
101  if self.doReturnResults is False:
102 
103  - ``exitStatus``: 0 if the task completed successfully, 1 otherwise.
104 
105  if self.doReturnResults is True:
106 
107  - ``result``: the result of calling jointcal.runDataRef()
108  - ``exitStatus``: 0 if the task completed successfully, 1 otherwise.
109  """
110  exitStatus = 0 # exit status for shell
111 
112  # NOTE: cannot call self.makeTask because that assumes args[0] is a single dataRef.
113  dataRefList, kwargs = args
114  butler = kwargs.pop('butler')
115  task = self.TaskClass(config=self.config, log=self.log, butler=butler)
116  result = None
117  try:
118  result = task.runDataRef(dataRefList, **kwargs)
119  exitStatus = result.exitStatus
120  job_path = butler.get('verify_job_filename')
121  result.job.write(job_path[0])
122  except Exception as e: # catch everything, sort it out later.
123  if self.doRaise:
124  raise e
125  else:
126  exitStatus = 1
127  eName = type(e).__name__
128  tract = dataRefList[0].dataId['tract']
129  task.log.fatal("Failed processing tract %s, %s: %s", tract, eName, e)
130 
131  # Put the butler back into kwargs for the other Tasks.
132  kwargs['butler'] = butler
133  if self.doReturnResults:
134  return pipeBase.Struct(result=result, exitStatus=exitStatus)
135  else:
136  return pipeBase.Struct(exitStatus=exitStatus)
137 
138 
139 class JointcalConfig(pexConfig.Config):
140  """Configuration for JointcalTask"""
141 
142  doAstrometry = pexConfig.Field(
143  doc="Fit astrometry and write the fitted result.",
144  dtype=bool,
145  default=True
146  )
147  doPhotometry = pexConfig.Field(
148  doc="Fit photometry and write the fitted result.",
149  dtype=bool,
150  default=True
151  )
152  coaddName = pexConfig.Field(
153  doc="Type of coadd, typically deep or goodSeeing",
154  dtype=str,
155  default="deep"
156  )
157  positionErrorPedestal = pexConfig.Field(
158  doc="Systematic term to apply to the measured position error (pixels)",
159  dtype=float,
160  default=0.02,
161  )
162  photometryErrorPedestal = pexConfig.Field(
163  doc="Systematic term to apply to the measured error on flux or magnitude as a "
164  "fraction of source flux or magnitude delta (e.g. 0.05 is 5% of flux or +50 millimag).",
165  dtype=float,
166  default=0.0,
167  )
168  # TODO: DM-6885 matchCut should be an geom.Angle
169  matchCut = pexConfig.Field(
170  doc="Matching radius between fitted and reference stars (arcseconds)",
171  dtype=float,
172  default=3.0,
173  )
174  minMeasurements = pexConfig.Field(
175  doc="Minimum number of associated measured stars for a fitted star to be included in the fit",
176  dtype=int,
177  default=2,
178  )
179  minMeasuredStarsPerCcd = pexConfig.Field(
180  doc="Minimum number of measuredStars per ccdImage before printing warnings",
181  dtype=int,
182  default=100,
183  )
184  minRefStarsPerCcd = pexConfig.Field(
185  doc="Minimum number of measuredStars per ccdImage before printing warnings",
186  dtype=int,
187  default=30,
188  )
189  allowLineSearch = pexConfig.Field(
190  doc="Allow a line search during minimization, if it is reasonable for the model"
191  " (models with a significant non-linear component, e.g. constrainedPhotometry).",
192  dtype=bool,
193  default=False
194  )
195  astrometrySimpleOrder = pexConfig.Field(
196  doc="Polynomial order for fitting the simple astrometry model.",
197  dtype=int,
198  default=3,
199  )
200  astrometryChipOrder = pexConfig.Field(
201  doc="Order of the per-chip transform for the constrained astrometry model.",
202  dtype=int,
203  default=1,
204  )
205  astrometryVisitOrder = pexConfig.Field(
206  doc="Order of the per-visit transform for the constrained astrometry model.",
207  dtype=int,
208  default=5,
209  )
210  useInputWcs = pexConfig.Field(
211  doc="Use the input calexp WCSs to initialize a SimpleAstrometryModel.",
212  dtype=bool,
213  default=True,
214  )
215  astrometryModel = pexConfig.ChoiceField(
216  doc="Type of model to fit to astrometry",
217  dtype=str,
218  default="constrained",
219  allowed={"simple": "One polynomial per ccd",
220  "constrained": "One polynomial per ccd, and one polynomial per visit"}
221  )
222  photometryModel = pexConfig.ChoiceField(
223  doc="Type of model to fit to photometry",
224  dtype=str,
225  default="constrainedMagnitude",
226  allowed={"simpleFlux": "One constant zeropoint per ccd and visit, fitting in flux space.",
227  "constrainedFlux": "Constrained zeropoint per ccd, and one polynomial per visit,"
228  " fitting in flux space.",
229  "simpleMagnitude": "One constant zeropoint per ccd and visit,"
230  " fitting in magnitude space.",
231  "constrainedMagnitude": "Constrained zeropoint per ccd, and one polynomial per visit,"
232  " fitting in magnitude space.",
233  }
234  )
235  applyColorTerms = pexConfig.Field(
236  doc="Apply photometric color terms to reference stars?"
237  "Requires that colorterms be set to a ColortermLibrary",
238  dtype=bool,
239  default=False
240  )
241  colorterms = pexConfig.ConfigField(
242  doc="Library of photometric reference catalog name to color term dict.",
243  dtype=ColortermLibrary,
244  )
245  photometryVisitOrder = pexConfig.Field(
246  doc="Order of the per-visit polynomial transform for the constrained photometry model.",
247  dtype=int,
248  default=7,
249  )
250  photometryDoRankUpdate = pexConfig.Field(
251  doc=("Do the rank update step during minimization. "
252  "Skipping this can help deal with models that are too non-linear."),
253  dtype=bool,
254  default=True,
255  )
256  astrometryDoRankUpdate = pexConfig.Field(
257  doc=("Do the rank update step during minimization (should not change the astrometry fit). "
258  "Skipping this can help deal with models that are too non-linear."),
259  dtype=bool,
260  default=True,
261  )
262  outlierRejectSigma = pexConfig.Field(
263  doc="How many sigma to reject outliers at during minimization.",
264  dtype=float,
265  default=5.0,
266  )
267  maxPhotometrySteps = pexConfig.Field(
268  doc="Maximum number of minimize iterations to take when fitting photometry.",
269  dtype=int,
270  default=20,
271  )
272  maxAstrometrySteps = pexConfig.Field(
273  doc="Maximum number of minimize iterations to take when fitting photometry.",
274  dtype=int,
275  default=20,
276  )
277  astrometryRefObjLoader = pexConfig.ConfigurableField(
278  target=LoadIndexedReferenceObjectsTask,
279  doc="Reference object loader for astrometric fit",
280  )
281  photometryRefObjLoader = pexConfig.ConfigurableField(
282  target=LoadIndexedReferenceObjectsTask,
283  doc="Reference object loader for photometric fit",
284  )
285  sourceSelector = sourceSelectorRegistry.makeField(
286  doc="How to select sources for cross-matching",
287  default="astrometry"
288  )
289  astrometryReferenceSelector = pexConfig.ConfigurableField(
290  target=ReferenceSourceSelectorTask,
291  doc="How to down-select the loaded astrometry reference catalog.",
292  )
293  photometryReferenceSelector = pexConfig.ConfigurableField(
294  target=ReferenceSourceSelectorTask,
295  doc="How to down-select the loaded photometry reference catalog.",
296  )
297  astrometryReferenceErr = pexConfig.Field(
298  doc=("Uncertainty on reference catalog coordinates [mas] to use in place of the `coord_*Err` fields. "
299  "If None, then raise an exception if the reference catalog is missing coordinate errors. "
300  "If specified, overrides any existing `coord_*Err` values."),
301  dtype=float,
302  default=None,
303  optional=True
304  )
305  writeInitMatrix = pexConfig.Field(
306  dtype=bool,
307  doc=("Write the pre/post-initialization Hessian and gradient to text files, for debugging. "
308  "The output files will be of the form 'astrometry_preinit-mat.txt', in the current directory. "
309  "Note that these files are the dense versions of the matrix, and so may be very large."),
310  default=False
311  )
312  writeChi2FilesInitialFinal = pexConfig.Field(
313  dtype=bool,
314  doc="Write .csv files containing the contributions to chi2 for the initialization and final fit.",
315  default=False
316  )
317  writeChi2FilesOuterLoop = pexConfig.Field(
318  dtype=bool,
319  doc="Write .csv files containing the contributions to chi2 for the outer fit loop.",
320  default=False
321  )
322  writeInitialModel = pexConfig.Field(
323  dtype=bool,
324  doc=("Write the pre-initialization model to text files, for debugging."
325  " Output is written to `initial[Astro|Photo]metryModel.txt` in the current working directory."),
326  default=False
327  )
328  debugOutputPath = pexConfig.Field(
329  dtype=str,
330  default=".",
331  doc=("Path to write debug output files to. Used by "
332  "`writeInitialModel`, `writeChi2Files*`, `writeInitMatrix`.")
333  )
334  sourceFluxType = pexConfig.Field(
335  dtype=str,
336  doc="Source flux field to use in source selection and to get fluxes from the catalog.",
337  default='Calib'
338  )
339 
340  def validate(self):
341  super().validate()
342  if self.doPhotometry and self.applyColorTerms and len(self.colorterms.data) == 0:
343  msg = "applyColorTerms=True requires the `colorterms` field be set to a ColortermLibrary."
344  raise pexConfig.FieldValidationError(JointcalConfig.colorterms, self, msg)
345  if self.doAstrometry and not self.doPhotometry and self.applyColorTerms:
346  msg = ("Only doing astrometry, but Colorterms are not applied for astrometry;"
347  "applyColorTerms=True will be ignored.")
348  lsst.log.warn(msg)
349 
350  def setDefaults(self):
351  # Use science source selector which can filter on extendedness, SNR, and whether blended
352  self.sourceSelector.name = 'science'
353  # Use only stars because aperture fluxes of galaxies are biased and depend on seeing
354  self.sourceSelector['science'].doUnresolved = True
355  # with dependable signal to noise ratio.
356  self.sourceSelector['science'].doSignalToNoise = True
357  # Min SNR must be > 0 because jointcal cannot handle negative fluxes,
358  # and S/N > 10 to use sources that are not too faint, and thus better measured.
359  self.sourceSelector['science'].signalToNoise.minimum = 10.
360  # Base SNR on CalibFlux because that is the flux jointcal that fits and must be positive
361  fluxField = f"slot_{self.sourceFluxType}Flux_instFlux"
362  self.sourceSelector['science'].signalToNoise.fluxField = fluxField
363  self.sourceSelector['science'].signalToNoise.errField = fluxField + "Err"
364  # Do not trust blended sources' aperture fluxes which also depend on seeing
365  self.sourceSelector['science'].doIsolated = True
366  # Do not trust either flux or centroid measurements with flags,
367  # chosen from the usual QA flags for stars)
368  self.sourceSelector['science'].doFlags = True
369  badFlags = ['base_PixelFlags_flag_edge', 'base_PixelFlags_flag_saturated',
370  'base_PixelFlags_flag_interpolatedCenter', 'base_SdssCentroid_flag',
371  'base_PsfFlux_flag', 'base_PixelFlags_flag_suspectCenter']
372  self.sourceSelector['science'].flags.bad = badFlags
373 
374  # Default to Gaia-DR2 (with proper motions) for astrometry and
375  # PS1-DR1 for photometry, with a reasonable initial filterMap.
376  self.astrometryRefObjLoader.ref_dataset_name = "gaia_dr2_20200414"
377  self.astrometryRefObjLoader.requireProperMotion = True
378  self.astrometryRefObjLoader.anyFilterMapsToThis = 'phot_g_mean'
379  self.photometryRefObjLoader.ref_dataset_name = "ps1_pv3_3pi_20170110"
380 
381 
382 def writeModel(model, filename, log):
383  """Write model to outfile."""
384  with open(filename, "w") as file:
385  file.write(repr(model))
386  log.info("Wrote %s to file: %s", model, filename)
387 
388 
389 class JointcalTask(pipeBase.CmdLineTask):
390  """Jointly astrometrically and photometrically calibrate a group of images."""
391 
392  ConfigClass = JointcalConfig
393  RunnerClass = JointcalRunner
394  _DefaultName = "jointcal"
395 
396  def __init__(self, butler=None, profile_jointcal=False, **kwargs):
397  """
398  Instantiate a JointcalTask.
399 
400  Parameters
401  ----------
402  butler : `lsst.daf.persistence.Butler`
403  The butler is passed to the refObjLoader constructor in case it is
404  needed. Ignored if the refObjLoader argument provides a loader directly.
405  Used to initialize the astrometry and photometry refObjLoaders.
406  profile_jointcal : `bool`
407  Set to True to profile different stages of this jointcal run.
408  """
409  pipeBase.CmdLineTask.__init__(self, **kwargs)
410  self.profile_jointcal = profile_jointcal
411  self.makeSubtask("sourceSelector")
412  if self.config.doAstrometry:
413  self.makeSubtask('astrometryRefObjLoader', butler=butler)
414  self.makeSubtask("astrometryReferenceSelector")
415  else:
417  if self.config.doPhotometry:
418  self.makeSubtask('photometryRefObjLoader', butler=butler)
419  self.makeSubtask("photometryReferenceSelector")
420  else:
422 
423  # To hold various computed metrics for use by tests
424  self.job = Job.load_metrics_package(subset='jointcal')
425 
426  # We don't currently need to persist the metadata.
427  # If we do in the future, we will have to add appropriate dataset templates
428  # to each obs package (the metadata template should look like `jointcal_wcs`).
429  def _getMetadataName(self):
430  return None
431 
432  @classmethod
433  def _makeArgumentParser(cls):
434  """Create an argument parser"""
435  parser = pipeBase.ArgumentParser(name=cls._DefaultName)
436  parser.add_argument("--profile_jointcal", default=False, action="store_true",
437  help="Profile steps of jointcal separately.")
438  parser.add_id_argument("--id", "calexp", help="data ID, e.g. --id visit=6789 ccd=0..9",
439  ContainerClass=PerTractCcdDataIdContainer)
440  return parser
441 
442  def _build_ccdImage(self, dataRef, associations, jointcalControl):
443  """
444  Extract the necessary things from this dataRef to add a new ccdImage.
445 
446  Parameters
447  ----------
448  dataRef : `lsst.daf.persistence.ButlerDataRef`
449  DataRef to extract info from.
450  associations : `lsst.jointcal.Associations`
451  Object to add the info to, to construct a new CcdImage
452  jointcalControl : `jointcal.JointcalControl`
453  Control object for associations management
454 
455  Returns
456  ------
457  namedtuple
458  ``wcs``
459  The TAN WCS of this image, read from the calexp
460  (`lsst.afw.geom.SkyWcs`).
461  ``key``
462  A key to identify this dataRef by its visit and ccd ids
463  (`namedtuple`).
464  ``filter``
465  This calexp's filter (`str`).
466  """
467  if "visit" in dataRef.dataId.keys():
468  visit = dataRef.dataId["visit"]
469  else:
470  visit = dataRef.getButler().queryMetadata("calexp", ("visit"), dataRef.dataId)[0]
471 
472  src = dataRef.get("src", flags=lsst.afw.table.SOURCE_IO_NO_FOOTPRINTS, immediate=True)
473 
474  visitInfo = dataRef.get('calexp_visitInfo')
475  detector = dataRef.get('calexp_detector')
476  ccdId = detector.getId()
477  photoCalib = dataRef.get('calexp_photoCalib')
478  tanWcs = dataRef.get('calexp_wcs')
479  bbox = dataRef.get('calexp_bbox')
480  filt = dataRef.get('calexp_filter')
481  filterName = filt.getName()
482 
483  goodSrc = self.sourceSelector.run(src)
484 
485  if len(goodSrc.sourceCat) == 0:
486  self.log.warn("No sources selected in visit %s ccd %s", visit, ccdId)
487  else:
488  self.log.info("%d sources selected in visit %d ccd %d", len(goodSrc.sourceCat), visit, ccdId)
489  associations.createCcdImage(goodSrc.sourceCat,
490  tanWcs,
491  visitInfo,
492  bbox,
493  filterName,
494  photoCalib,
495  detector,
496  visit,
497  ccdId,
498  jointcalControl)
499 
500  Result = collections.namedtuple('Result_from_build_CcdImage', ('wcs', 'key', 'filter'))
501  Key = collections.namedtuple('Key', ('visit', 'ccd'))
502  return Result(tanWcs, Key(visit, ccdId), filterName)
503 
504  def _getDebugPath(self, filename):
505  """Constructs a path to filename using the configured debug path.
506  """
507  return os.path.join(self.config.debugOutputPath, filename)
508 
509  @pipeBase.timeMethod
510  def runDataRef(self, dataRefs, profile_jointcal=False):
511  """
512  Jointly calibrate the astrometry and photometry across a set of images.
513 
514  Parameters
515  ----------
516  dataRefs : `list` of `lsst.daf.persistence.ButlerDataRef`
517  List of data references to the exposures to be fit.
518  profile_jointcal : `bool`
519  Profile the individual steps of jointcal.
520 
521  Returns
522  -------
523  result : `lsst.pipe.base.Struct`
524  Struct of metadata from the fit, containing:
525 
526  ``dataRefs``
527  The provided data references that were fit (with updated WCSs)
528  ``oldWcsList``
529  The original WCS from each dataRef
530  ``metrics``
531  Dictionary of internally-computed metrics for testing/validation.
532  """
533  if len(dataRefs) == 0:
534  raise ValueError('Need a non-empty list of data references!')
535 
536  exitStatus = 0 # exit status for shell
537 
538  sourceFluxField = "slot_%sFlux" % (self.config.sourceFluxType,)
539  jointcalControl = lsst.jointcal.JointcalControl(sourceFluxField)
540  associations = lsst.jointcal.Associations()
541 
542  visit_ccd_to_dataRef = {}
543  oldWcsList = []
544  filters = []
545  load_cat_prof_file = 'jointcal_build_ccdImage.prof' if profile_jointcal else ''
546  with pipeBase.cmdLineTask.profile(load_cat_prof_file):
547  # We need the bounding-box of the focal plane for photometry visit models.
548  # NOTE: we only need to read it once, because its the same for all exposures of a camera.
549  camera = dataRefs[0].get('camera', immediate=True)
550  self.focalPlaneBBox = camera.getFpBBox()
551  for ref in dataRefs:
552  result = self._build_ccdImage(ref, associations, jointcalControl)
553  oldWcsList.append(result.wcs)
554  visit_ccd_to_dataRef[result.key] = ref
555  filters.append(result.filter)
556  filters = collections.Counter(filters)
557 
558  associations.computeCommonTangentPoint()
559 
560  boundingCircle = associations.computeBoundingCircle()
561  center = lsst.geom.SpherePoint(boundingCircle.getCenter())
562  radius = lsst.geom.Angle(boundingCircle.getOpeningAngle().asRadians(), lsst.geom.radians)
563 
564  self.log.info(f"Data has center={center} with radius={radius.asDegrees()} degrees.")
565 
566  # Determine a default filter associated with the catalog. See DM-9093
567  defaultFilter = filters.most_common(1)[0][0]
568  self.log.debug("Using %s band for reference flux", defaultFilter)
569 
570  # TODO: need a better way to get the tract.
571  tract = dataRefs[0].dataId['tract']
572 
573  if self.config.doAstrometry:
574  astrometry = self._do_load_refcat_and_fit(associations, defaultFilter, center, radius,
575  name="astrometry",
576  refObjLoader=self.astrometryRefObjLoader,
577  referenceSelector=self.astrometryReferenceSelector,
578  fit_function=self._fit_astrometry,
579  profile_jointcal=profile_jointcal,
580  tract=tract)
581  self._write_astrometry_results(associations, astrometry.model, visit_ccd_to_dataRef)
582  else:
583  astrometry = Astrometry(None, None, None)
584 
585  if self.config.doPhotometry:
586  photometry = self._do_load_refcat_and_fit(associations, defaultFilter, center, radius,
587  name="photometry",
588  refObjLoader=self.photometryRefObjLoader,
589  referenceSelector=self.photometryReferenceSelector,
590  fit_function=self._fit_photometry,
591  profile_jointcal=profile_jointcal,
592  tract=tract,
593  filters=filters,
594  reject_bad_fluxes=True)
595  self._write_photometry_results(associations, photometry.model, visit_ccd_to_dataRef)
596  else:
597  photometry = Photometry(None, None)
598 
599  return pipeBase.Struct(dataRefs=dataRefs,
600  oldWcsList=oldWcsList,
601  job=self.job,
602  astrometryRefObjLoader=self.astrometryRefObjLoader,
603  photometryRefObjLoader=self.photometryRefObjLoader,
604  defaultFilter=defaultFilter,
605  exitStatus=exitStatus)
606 
607  def _get_refcat_coordinate_error_override(self, refCat, name):
608  """Check whether we should override the refcat coordinate errors, and
609  return the overridden error if necessary.
610 
611  Parameters
612  ----------
613  refCat : `lsst.afw.table.SimpleCatalog`
614  The reference catalog to check for a ``coord_raErr`` field.
615  name : `str`
616  Whether we are doing "astrometry" or "photometry".
617 
618  Returns
619  -------
620  refCoordErr : `float`
621  The refcat coordinate error to use, or NaN if we are not overriding
622  those fields.
623 
624  Raises
625  ------
626  lsst.pex.config.FieldValidationError
627  Raised if the refcat does not contain coordinate errors and
628  ``config.astrometryReferenceErr`` is not set.
629  """
630  # This value doesn't matter for photometry, so just set something to
631  # keep old refcats from causing problems.
632  if name.lower() == "photometry":
633  if 'coord_raErr' not in refCat.schema:
634  return 100
635  else:
636  return float('nan')
637 
638  if self.config.astrometryReferenceErr is None and 'coord_raErr' not in refCat.schema:
639  msg = ("Reference catalog does not contain coordinate errors, "
640  "and config.astrometryReferenceErr not supplied.")
641  raise pexConfig.FieldValidationError(JointcalConfig.astrometryReferenceErr,
642  self.config,
643  msg)
644 
645  if self.config.astrometryReferenceErr is not None and 'coord_raErr' in refCat.schema:
646  self.log.warn("Overriding reference catalog coordinate errors with %f/coordinate [mas]",
647  self.config.astrometryReferenceErr)
648 
649  if self.config.astrometryReferenceErr is None:
650  return float('nan')
651  else:
652  return self.config.astrometryReferenceErr
653 
654  def _compute_proper_motion_epoch(self, ccdImageList):
655  """Return the proper motion correction epoch of the provided images.
656 
657  Parameters
658  ----------
659  ccdImageList : `list` [`lsst.jointcal.CcdImage`]
660  The images to compute the appropriate epoch for.
661 
662  Returns
663  -------
664  epoch : `astropy.time.Time`
665  The date to use for proper motion corrections.
666  """
667  mjds = [ccdImage.getMjd() for ccdImage in ccdImageList]
668  return astropy.time.Time(np.mean(mjds), format='mjd', scale="tai")
669 
670  def _do_load_refcat_and_fit(self, associations, defaultFilter, center, radius,
671  filters=[],
672  tract="", profile_jointcal=False, match_cut=3.0,
673  reject_bad_fluxes=False, *,
674  name="", refObjLoader=None, referenceSelector=None,
675  fit_function=None):
676  """Load reference catalog, perform the fit, and return the result.
677 
678  Parameters
679  ----------
680  associations : `lsst.jointcal.Associations`
681  The star/reference star associations to fit.
682  defaultFilter : `str`
683  filter to load from reference catalog.
684  center : `lsst.geom.SpherePoint`
685  ICRS center of field to load from reference catalog.
686  radius : `lsst.geom.Angle`
687  On-sky radius to load from reference catalog.
688  name : `str`
689  Name of thing being fit: "astrometry" or "photometry".
690  refObjLoader : `lsst.meas.algorithms.LoadReferenceObjectsTask`
691  Reference object loader to use to load a reference catalog.
692  referenceSelector : `lsst.meas.algorithms.ReferenceSourceSelectorTask`
693  Selector to use to pick objects from the loaded reference catalog.
694  fit_function : callable
695  Function to call to perform fit (takes Associations object).
696  filters : `list` [`str`], optional
697  List of filters to load from the reference catalog.
698  tract : `str`, optional
699  Name of tract currently being fit.
700  profile_jointcal : `bool`, optional
701  Separately profile the fitting step.
702  match_cut : `float`, optional
703  Radius in arcseconds to find cross-catalog matches to during
704  associations.associateCatalogs.
705  reject_bad_fluxes : `bool`, optional
706  Reject refCat sources with NaN/inf flux or NaN/0 fluxErr.
707 
708  Returns
709  -------
710  result : `Photometry` or `Astrometry`
711  Result of `fit_function()`
712  """
713  self.log.info("====== Now processing %s...", name)
714  # TODO: this should not print "trying to invert a singular transformation:"
715  # if it does that, something's not right about the WCS...
716  associations.associateCatalogs(match_cut)
717  add_measurement(self.job, 'jointcal.associated_%s_fittedStars' % name,
718  associations.fittedStarListSize())
719 
720  applyColorterms = False if name.lower() == "astrometry" else self.config.applyColorTerms
721  epoch = self._compute_proper_motion_epoch(associations.getCcdImageList())
722  refCat, fluxField = self._load_reference_catalog(refObjLoader, referenceSelector,
723  center, radius, defaultFilter,
724  applyColorterms=applyColorterms,
725  epoch=epoch)
726  refCoordErr = self._get_refcat_coordinate_error_override(refCat, name)
727 
728  associations.collectRefStars(refCat,
729  self.config.matchCut*lsst.geom.arcseconds,
730  fluxField,
731  refCoordinateErr=refCoordErr,
732  rejectBadFluxes=reject_bad_fluxes)
733  add_measurement(self.job, 'jointcal.collected_%s_refStars' % name,
734  associations.refStarListSize())
735 
736  associations.prepareFittedStars(self.config.minMeasurements)
737 
738  self._check_star_lists(associations, name)
739  add_measurement(self.job, 'jointcal.selected_%s_refStars' % name,
740  associations.nFittedStarsWithAssociatedRefStar())
741  add_measurement(self.job, 'jointcal.selected_%s_fittedStars' % name,
742  associations.fittedStarListSize())
743  add_measurement(self.job, 'jointcal.selected_%s_ccdImages' % name,
744  associations.nCcdImagesValidForFit())
745 
746  load_cat_prof_file = 'jointcal_fit_%s.prof'%name if profile_jointcal else ''
747  dataName = "{}_{}".format(tract, defaultFilter)
748  with pipeBase.cmdLineTask.profile(load_cat_prof_file):
749  result = fit_function(associations, dataName)
750  # TODO DM-12446: turn this into a "butler save" somehow.
751  # Save reference and measurement chi2 contributions for this data
752  if self.config.writeChi2FilesInitialFinal:
753  baseName = self._getDebugPath(f"{name}_final_chi2-{dataName}")
754  result.fit.saveChi2Contributions(baseName+"{type}")
755  self.log.info("Wrote chi2 contributions files: %s", baseName)
756 
757  return result
758 
759  def _load_reference_catalog(self, refObjLoader, referenceSelector, center, radius, filterName,
760  applyColorterms=False, epoch=None):
761  """Load the necessary reference catalog sources, convert fluxes to
762  correct units, and apply color term corrections if requested.
763 
764  Parameters
765  ----------
766  refObjLoader : `lsst.meas.algorithms.LoadReferenceObjectsTask`
767  The reference catalog loader to use to get the data.
768  referenceSelector : `lsst.meas.algorithms.ReferenceSourceSelectorTask`
769  Source selector to apply to loaded reference catalog.
770  center : `lsst.geom.SpherePoint`
771  The center around which to load sources.
772  radius : `lsst.geom.Angle`
773  The radius around ``center`` to load sources in.
774  filterName : `str`
775  The name of the camera filter to load fluxes for.
776  applyColorterms : `bool`
777  Apply colorterm corrections to the refcat for ``filterName``?
778  epoch : `astropy.time.Time`, optional
779  Epoch to which to correct refcat proper motion and parallax,
780  or `None` to not apply such corrections.
781 
782  Returns
783  -------
784  refCat : `lsst.afw.table.SimpleCatalog`
785  The loaded reference catalog.
786  fluxField : `str`
787  The name of the reference catalog flux field appropriate for ``filterName``.
788  """
789  skyCircle = refObjLoader.loadSkyCircle(center,
790  radius,
791  filterName,
792  epoch=epoch)
793 
794  selected = referenceSelector.run(skyCircle.refCat)
795  # Need memory contiguity to get reference filters as a vector.
796  if not selected.sourceCat.isContiguous():
797  refCat = selected.sourceCat.copy(deep=True)
798  else:
799  refCat = selected.sourceCat
800 
801  if applyColorterms:
802  refCatName = refObjLoader.ref_dataset_name
803  self.log.info("Applying color terms for filterName=%r reference catalog=%s",
804  filterName, refCatName)
805  colorterm = self.config.colorterms.getColorterm(
806  filterName=filterName, photoCatName=refCatName, doRaise=True)
807 
808  refMag, refMagErr = colorterm.getCorrectedMagnitudes(refCat, filterName)
809  refCat[skyCircle.fluxField] = u.Magnitude(refMag, u.ABmag).to_value(u.nJy)
810  # TODO: I didn't want to use this, but I'll deal with it in DM-16903
811  refCat[skyCircle.fluxField+'Err'] = fluxErrFromABMagErr(refMagErr, refMag) * 1e9
812 
813  return refCat, skyCircle.fluxField
814 
815  def _check_star_lists(self, associations, name):
816  # TODO: these should be len(blah), but we need this properly wrapped first.
817  if associations.nCcdImagesValidForFit() == 0:
818  raise RuntimeError('No images in the ccdImageList!')
819  if associations.fittedStarListSize() == 0:
820  raise RuntimeError('No stars in the {} fittedStarList!'.format(name))
821  if associations.refStarListSize() == 0:
822  raise RuntimeError('No stars in the {} reference star list!'.format(name))
823 
824  def _logChi2AndValidate(self, associations, fit, model, chi2Label, writeChi2Name=None):
825  """Compute chi2, log it, validate the model, and return chi2.
826 
827  Parameters
828  ----------
829  associations : `lsst.jointcal.Associations`
830  The star/reference star associations to fit.
831  fit : `lsst.jointcal.FitterBase`
832  The fitter to use for minimization.
833  model : `lsst.jointcal.Model`
834  The model being fit.
835  chi2Label : `str`
836  Label to describe the chi2 (e.g. "Initialized", "Final").
837  writeChi2Name : `str`, optional
838  Filename prefix to write the chi2 contributions to.
839  Do not supply an extension: an appropriate one will be added.
840 
841  Returns
842  -------
843  chi2: `lsst.jointcal.Chi2Accumulator`
844  The chi2 object for the current fitter and model.
845 
846  Raises
847  ------
848  FloatingPointError
849  Raised if chi2 is infinite or NaN.
850  ValueError
851  Raised if the model is not valid.
852  """
853  if writeChi2Name is not None:
854  fullpath = self._getDebugPath(writeChi2Name)
855  fit.saveChi2Contributions(fullpath+"{type}")
856  self.log.info("Wrote chi2 contributions files: %s", fullpath)
857 
858  chi2 = fit.computeChi2()
859  self.log.info("%s %s", chi2Label, chi2)
860  self._check_stars(associations)
861  if not np.isfinite(chi2.chi2):
862  raise FloatingPointError(f'{chi2Label} chi2 is invalid: {chi2}')
863  if not model.validate(associations.getCcdImageList(), chi2.ndof):
864  raise ValueError("Model is not valid: check log messages for warnings.")
865  return chi2
866 
867  def _fit_photometry(self, associations, dataName=None):
868  """
869  Fit the photometric data.
870 
871  Parameters
872  ----------
873  associations : `lsst.jointcal.Associations`
874  The star/reference star associations to fit.
875  dataName : `str`
876  Name of the data being processed (e.g. "1234_HSC-Y"), for
877  identifying debugging files.
878 
879  Returns
880  -------
881  fit_result : `namedtuple`
882  fit : `lsst.jointcal.PhotometryFit`
883  The photometric fitter used to perform the fit.
884  model : `lsst.jointcal.PhotometryModel`
885  The photometric model that was fit.
886  """
887  self.log.info("=== Starting photometric fitting...")
888 
889  # TODO: should use pex.config.RegistryField here (see DM-9195)
890  if self.config.photometryModel == "constrainedFlux":
891  model = lsst.jointcal.ConstrainedFluxModel(associations.getCcdImageList(),
892  self.focalPlaneBBox,
893  visitOrder=self.config.photometryVisitOrder,
894  errorPedestal=self.config.photometryErrorPedestal)
895  # potentially nonlinear problem, so we may need a line search to converge.
896  doLineSearch = self.config.allowLineSearch
897  elif self.config.photometryModel == "constrainedMagnitude":
898  model = lsst.jointcal.ConstrainedMagnitudeModel(associations.getCcdImageList(),
899  self.focalPlaneBBox,
900  visitOrder=self.config.photometryVisitOrder,
901  errorPedestal=self.config.photometryErrorPedestal)
902  # potentially nonlinear problem, so we may need a line search to converge.
903  doLineSearch = self.config.allowLineSearch
904  elif self.config.photometryModel == "simpleFlux":
905  model = lsst.jointcal.SimpleFluxModel(associations.getCcdImageList(),
906  errorPedestal=self.config.photometryErrorPedestal)
907  doLineSearch = False # purely linear in model parameters, so no line search needed
908  elif self.config.photometryModel == "simpleMagnitude":
909  model = lsst.jointcal.SimpleMagnitudeModel(associations.getCcdImageList(),
910  errorPedestal=self.config.photometryErrorPedestal)
911  doLineSearch = False # purely linear in model parameters, so no line search needed
912 
913  fit = lsst.jointcal.PhotometryFit(associations, model)
914  # TODO DM-12446: turn this into a "butler save" somehow.
915  # Save reference and measurement chi2 contributions for this data
916  if self.config.writeChi2FilesInitialFinal:
917  baseName = f"photometry_initial_chi2-{dataName}"
918  else:
919  baseName = None
920  if self.config.writeInitialModel:
921  fullpath = self._getDebugPath("initialPhotometryModel.txt")
922  writeModel(model, fullpath, self.log)
923  self._logChi2AndValidate(associations, fit, model, "Initialized", writeChi2Name=baseName)
924 
925  def getChi2Name(whatToFit):
926  if self.config.writeChi2FilesOuterLoop:
927  return f"photometry_init-%s_chi2-{dataName}" % whatToFit
928  else:
929  return None
930 
931  # The constrained model needs the visit transform fit first; the chip
932  # transform is initialized from the singleFrame PhotoCalib, so it's close.
933  dumpMatrixFile = self._getDebugPath("photometry_preinit") if self.config.writeInitMatrix else ""
934  if self.config.photometryModel.startswith("constrained"):
935  # no line search: should be purely (or nearly) linear,
936  # and we want a large step size to initialize with.
937  fit.minimize("ModelVisit", dumpMatrixFile=dumpMatrixFile)
938  self._logChi2AndValidate(associations, fit, model, "Initialize ModelVisit",
939  writeChi2Name=getChi2Name("ModelVisit"))
940  dumpMatrixFile = "" # so we don't redo the output on the next step
941 
942  fit.minimize("Model", doLineSearch=doLineSearch, dumpMatrixFile=dumpMatrixFile)
943  self._logChi2AndValidate(associations, fit, model, "Initialize Model",
944  writeChi2Name=getChi2Name("Model"))
945 
946  fit.minimize("Fluxes") # no line search: always purely linear.
947  self._logChi2AndValidate(associations, fit, model, "Initialize Fluxes",
948  writeChi2Name=getChi2Name("Fluxes"))
949 
950  fit.minimize("Model Fluxes", doLineSearch=doLineSearch)
951  self._logChi2AndValidate(associations, fit, model, "Initialize ModelFluxes",
952  writeChi2Name=getChi2Name("ModelFluxes"))
953 
954  model.freezeErrorTransform()
955  self.log.debug("Photometry error scales are frozen.")
956 
957  chi2 = self._iterate_fit(associations,
958  fit,
959  self.config.maxPhotometrySteps,
960  "photometry",
961  "Model Fluxes",
962  doRankUpdate=self.config.photometryDoRankUpdate,
963  doLineSearch=doLineSearch,
964  dataName=dataName)
965 
966  add_measurement(self.job, 'jointcal.photometry_final_chi2', chi2.chi2)
967  add_measurement(self.job, 'jointcal.photometry_final_ndof', chi2.ndof)
968  return Photometry(fit, model)
969 
970  def _fit_astrometry(self, associations, dataName=None):
971  """
972  Fit the astrometric data.
973 
974  Parameters
975  ----------
976  associations : `lsst.jointcal.Associations`
977  The star/reference star associations to fit.
978  dataName : `str`
979  Name of the data being processed (e.g. "1234_HSC-Y"), for
980  identifying debugging files.
981 
982  Returns
983  -------
984  fit_result : `namedtuple`
985  fit : `lsst.jointcal.AstrometryFit`
986  The astrometric fitter used to perform the fit.
987  model : `lsst.jointcal.AstrometryModel`
988  The astrometric model that was fit.
989  sky_to_tan_projection : `lsst.jointcal.ProjectionHandler`
990  The model for the sky to tangent plane projection that was used in the fit.
991  """
992 
993  self.log.info("=== Starting astrometric fitting...")
994 
995  associations.deprojectFittedStars()
996 
997  # NOTE: need to return sky_to_tan_projection so that it doesn't get garbage collected.
998  # TODO: could we package sky_to_tan_projection and model together so we don't have to manage
999  # them so carefully?
1000  sky_to_tan_projection = lsst.jointcal.OneTPPerVisitHandler(associations.getCcdImageList())
1001 
1002  if self.config.astrometryModel == "constrained":
1003  model = lsst.jointcal.ConstrainedAstrometryModel(associations.getCcdImageList(),
1004  sky_to_tan_projection,
1005  chipOrder=self.config.astrometryChipOrder,
1006  visitOrder=self.config.astrometryVisitOrder)
1007  elif self.config.astrometryModel == "simple":
1008  model = lsst.jointcal.SimpleAstrometryModel(associations.getCcdImageList(),
1009  sky_to_tan_projection,
1010  self.config.useInputWcs,
1011  nNotFit=0,
1012  order=self.config.astrometrySimpleOrder)
1013 
1014  fit = lsst.jointcal.AstrometryFit(associations, model, self.config.positionErrorPedestal)
1015  # TODO DM-12446: turn this into a "butler save" somehow.
1016  # Save reference and measurement chi2 contributions for this data
1017  if self.config.writeChi2FilesInitialFinal:
1018  baseName = f"astrometry_initial_chi2-{dataName}"
1019  else:
1020  baseName = None
1021  if self.config.writeInitialModel:
1022  fullpath = self._getDebugPath("initialAstrometryModel.txt")
1023  writeModel(model, fullpath, self.log)
1024  self._logChi2AndValidate(associations, fit, model, "Initial", writeChi2Name=baseName)
1025 
1026  def getChi2Name(whatToFit):
1027  if self.config.writeChi2FilesOuterLoop:
1028  return f"astrometry_init-%s_chi2-{dataName}" % whatToFit
1029  else:
1030  return None
1031 
1032  dumpMatrixFile = self._getDebugPath("astrometry_preinit") if self.config.writeInitMatrix else ""
1033  # The constrained model needs the visit transform fit first; the chip
1034  # transform is initialized from the detector's cameraGeom, so it's close.
1035  if self.config.astrometryModel == "constrained":
1036  fit.minimize("DistortionsVisit", dumpMatrixFile=dumpMatrixFile)
1037  self._logChi2AndValidate(associations, fit, model, "Initialize DistortionsVisit",
1038  writeChi2Name=getChi2Name("DistortionsVisit"))
1039  dumpMatrixFile = "" # so we don't redo the output on the next step
1040 
1041  fit.minimize("Distortions", dumpMatrixFile=dumpMatrixFile)
1042  self._logChi2AndValidate(associations, fit, model, "Initialize Distortions",
1043  writeChi2Name=getChi2Name("Distortions"))
1044 
1045  fit.minimize("Positions")
1046  self._logChi2AndValidate(associations, fit, model, "Initialize Positions",
1047  writeChi2Name=getChi2Name("Positions"))
1048 
1049  fit.minimize("Distortions Positions")
1050  self._logChi2AndValidate(associations, fit, model, "Initialize DistortionsPositions",
1051  writeChi2Name=getChi2Name("DistortionsPositions"))
1052 
1053  chi2 = self._iterate_fit(associations,
1054  fit,
1055  self.config.maxAstrometrySteps,
1056  "astrometry",
1057  "Distortions Positions",
1058  doRankUpdate=self.config.astrometryDoRankUpdate,
1059  dataName=dataName)
1060 
1061  add_measurement(self.job, 'jointcal.astrometry_final_chi2', chi2.chi2)
1062  add_measurement(self.job, 'jointcal.astrometry_final_ndof', chi2.ndof)
1063 
1064  return Astrometry(fit, model, sky_to_tan_projection)
1065 
1066  def _check_stars(self, associations):
1067  """Count measured and reference stars per ccd and warn/log them."""
1068  for ccdImage in associations.getCcdImageList():
1069  nMeasuredStars, nRefStars = ccdImage.countStars()
1070  self.log.debug("ccdImage %s has %s measured and %s reference stars",
1071  ccdImage.getName(), nMeasuredStars, nRefStars)
1072  if nMeasuredStars < self.config.minMeasuredStarsPerCcd:
1073  self.log.warn("ccdImage %s has only %s measuredStars (desired %s)",
1074  ccdImage.getName(), nMeasuredStars, self.config.minMeasuredStarsPerCcd)
1075  if nRefStars < self.config.minRefStarsPerCcd:
1076  self.log.warn("ccdImage %s has only %s RefStars (desired %s)",
1077  ccdImage.getName(), nRefStars, self.config.minRefStarsPerCcd)
1078 
1079  def _iterate_fit(self, associations, fitter, max_steps, name, whatToFit,
1080  dataName="",
1081  doRankUpdate=True,
1082  doLineSearch=False):
1083  """Run fitter.minimize up to max_steps times, returning the final chi2.
1084 
1085  Parameters
1086  ----------
1087  associations : `lsst.jointcal.Associations`
1088  The star/reference star associations to fit.
1089  fitter : `lsst.jointcal.FitterBase`
1090  The fitter to use for minimization.
1091  max_steps : `int`
1092  Maximum number of steps to run outlier rejection before declaring
1093  convergence failure.
1094  name : {'photometry' or 'astrometry'}
1095  What type of data are we fitting (for logs and debugging files).
1096  whatToFit : `str`
1097  Passed to ``fitter.minimize()`` to define the parameters to fit.
1098  dataName : `str`, optional
1099  Descriptive name for this dataset (e.g. tract and filter),
1100  for debugging.
1101  doRankUpdate : `bool`, optional
1102  Do an Eigen rank update during minimization, or recompute the full
1103  matrix and gradient?
1104  doLineSearch : `bool`, optional
1105  Do a line search for the optimum step during minimization?
1106 
1107  Returns
1108  -------
1109  chi2: `lsst.jointcal.Chi2Statistic`
1110  The final chi2 after the fit converges, or is forced to end.
1111 
1112  Raises
1113  ------
1114  FloatingPointError
1115  Raised if the fitter fails with a non-finite value.
1116  RuntimeError
1117  Raised if the fitter fails for some other reason;
1118  log messages will provide further details.
1119  """
1120  dumpMatrixFile = self._getDebugPath(f"{name}_postinit") if self.config.writeInitMatrix else ""
1121  oldChi2 = lsst.jointcal.Chi2Statistic()
1122  oldChi2.chi2 = float("inf")
1123  for i in range(max_steps):
1124  if self.config.writeChi2FilesOuterLoop:
1125  writeChi2Name = f"{name}_iterate_{i}_chi2-{dataName}"
1126  else:
1127  writeChi2Name = None
1128  result = fitter.minimize(whatToFit,
1129  self.config.outlierRejectSigma,
1130  doRankUpdate=doRankUpdate,
1131  doLineSearch=doLineSearch,
1132  dumpMatrixFile=dumpMatrixFile)
1133  dumpMatrixFile = "" # clear it so we don't write the matrix again.
1134  chi2 = self._logChi2AndValidate(associations, fitter, fitter.getModel(),
1135  f"Fit iteration {i}", writeChi2Name=writeChi2Name)
1136 
1137  if result == MinimizeResult.Converged:
1138  if doRankUpdate:
1139  self.log.debug("fit has converged - no more outliers - redo minimization "
1140  "one more time in case we have lost accuracy in rank update.")
1141  # Redo minimization one more time in case we have lost accuracy in rank update
1142  result = fitter.minimize(whatToFit, self.config.outlierRejectSigma)
1143  chi2 = self._logChi2AndValidate(associations, fitter, fitter.getModel(), "Fit completed")
1144 
1145  # log a message for a large final chi2, TODO: DM-15247 for something better
1146  if chi2.chi2/chi2.ndof >= 4.0:
1147  self.log.error("Potentially bad fit: High chi-squared/ndof.")
1148 
1149  break
1150  elif result == MinimizeResult.Chi2Increased:
1151  self.log.warn("Still some outliers remaining but chi2 increased - retry")
1152  # Check whether the increase was large enough to cause trouble.
1153  chi2Ratio = chi2.chi2 / oldChi2.chi2
1154  if chi2Ratio > 1.5:
1155  self.log.warn('Significant chi2 increase by a factor of %.4g / %.4g = %.4g',
1156  chi2.chi2, oldChi2.chi2, chi2Ratio)
1157  # Based on a variety of HSC jointcal logs (see DM-25779), it
1158  # appears that chi2 increases more than a factor of ~2 always
1159  # result in the fit diverging rapidly and ending at chi2 > 1e10.
1160  # Using 10 as the "failure" threshold gives some room between
1161  # leaving a warning and bailing early.
1162  if chi2Ratio > 10:
1163  msg = ("Large chi2 increase between steps: fit likely cannot converge."
1164  " Try setting one or more of the `writeChi2*` config fields and looking"
1165  " at how individual star chi2-values evolve during the fit.")
1166  raise RuntimeError(msg)
1167  oldChi2 = chi2
1168  elif result == MinimizeResult.NonFinite:
1169  filename = self._getDebugPath("{}_failure-nonfinite_chi2-{}.csv".format(name, dataName))
1170  # TODO DM-12446: turn this into a "butler save" somehow.
1171  fitter.saveChi2Contributions(filename+"{type}")
1172  msg = "Nonfinite value in chi2 minimization, cannot complete fit. Dumped star tables to: {}"
1173  raise FloatingPointError(msg.format(filename))
1174  elif result == MinimizeResult.Failed:
1175  raise RuntimeError("Chi2 minimization failure, cannot complete fit.")
1176  else:
1177  raise RuntimeError("Unxepected return code from minimize().")
1178  else:
1179  self.log.error("%s failed to converge after %d steps"%(name, max_steps))
1180 
1181  return chi2
1182 
1183  def _write_astrometry_results(self, associations, model, visit_ccd_to_dataRef):
1184  """
1185  Write the fitted astrometric results to a new 'jointcal_wcs' dataRef.
1186 
1187  Parameters
1188  ----------
1189  associations : `lsst.jointcal.Associations`
1190  The star/reference star associations to fit.
1191  model : `lsst.jointcal.AstrometryModel`
1192  The astrometric model that was fit.
1193  visit_ccd_to_dataRef : `dict` of Key: `lsst.daf.persistence.ButlerDataRef`
1194  Dict of ccdImage identifiers to dataRefs that were fit.
1195  """
1196 
1197  ccdImageList = associations.getCcdImageList()
1198  for ccdImage in ccdImageList:
1199  # TODO: there must be a better way to identify this ccdImage than a visit,ccd pair?
1200  ccd = ccdImage.ccdId
1201  visit = ccdImage.visit
1202  dataRef = visit_ccd_to_dataRef[(visit, ccd)]
1203  self.log.info("Updating WCS for visit: %d, ccd: %d", visit, ccd)
1204  skyWcs = model.makeSkyWcs(ccdImage)
1205  try:
1206  dataRef.put(skyWcs, 'jointcal_wcs')
1207  except pexExceptions.Exception as e:
1208  self.log.fatal('Failed to write updated Wcs: %s', str(e))
1209  raise e
1210 
1211  def _write_photometry_results(self, associations, model, visit_ccd_to_dataRef):
1212  """
1213  Write the fitted photometric results to a new 'jointcal_photoCalib' dataRef.
1214 
1215  Parameters
1216  ----------
1217  associations : `lsst.jointcal.Associations`
1218  The star/reference star associations to fit.
1219  model : `lsst.jointcal.PhotometryModel`
1220  The photoometric model that was fit.
1221  visit_ccd_to_dataRef : `dict` of Key: `lsst.daf.persistence.ButlerDataRef`
1222  Dict of ccdImage identifiers to dataRefs that were fit.
1223  """
1224 
1225  ccdImageList = associations.getCcdImageList()
1226  for ccdImage in ccdImageList:
1227  # TODO: there must be a better way to identify this ccdImage than a visit,ccd pair?
1228  ccd = ccdImage.ccdId
1229  visit = ccdImage.visit
1230  dataRef = visit_ccd_to_dataRef[(visit, ccd)]
1231  self.log.info("Updating PhotoCalib for visit: %d, ccd: %d", visit, ccd)
1232  photoCalib = model.toPhotoCalib(ccdImage)
1233  try:
1234  dataRef.put(photoCalib, 'jointcal_photoCalib')
1235  except pexExceptions.Exception as e:
1236  self.log.fatal('Failed to write updated PhotoCalib: %s', str(e))
1237  raise e
lsst::afw::image
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects.
Definition: imageAlgorithm.dox:1
lsst::jointcal.jointcal.JointcalTask._write_photometry_results
def _write_photometry_results(self, associations, model, visit_ccd_to_dataRef)
Definition: jointcal.py:1211
lsst::log.log.logContinued.warn
def warn(fmt, *args)
Definition: logContinued.py:205
lsst::jointcal.jointcal.JointcalTask.__init__
def __init__(self, butler=None, profile_jointcal=False, **kwargs)
Definition: jointcal.py:396
lsst::log.log.logContinued.info
def info(fmt, *args)
Definition: logContinued.py:201
lsst::jointcal.jointcal.JointcalTask.astrometryRefObjLoader
astrometryRefObjLoader
Definition: jointcal.py:416
lsst::log.log.logContinued.error
def error(fmt, *args)
Definition: logContinued.py:213
lsst::jointcal.jointcal.JointcalTask._getDebugPath
def _getDebugPath(self, filename)
Definition: jointcal.py:504
lsst::jointcal.jointcal.Astrometry
Astrometry
Definition: jointcal.py:52
lsst::meas::algorithms.sourceSelector
Definition: sourceSelector.py:1
lsst::jointcal.jointcal.JointcalConfig
Definition: jointcal.py:139
lsst::jointcal.jointcal.JointcalTask._fit_photometry
def _fit_photometry(self, associations, dataName=None)
Definition: jointcal.py:867
lsst::jointcal.jointcal.add_measurement
def add_measurement(job, name, value)
Definition: jointcal.py:56
lsst::jointcal.jointcal.JointcalRunner.getTargetList
def getTargetList(parsedCmd, **kwargs)
Definition: jointcal.py:74
lsst.gdb.afw.printers.debug
bool debug
Definition: printers.py:9
lsst::jointcal.jointcal.JointcalTask._load_reference_catalog
def _load_reference_catalog(self, refObjLoader, referenceSelector, center, radius, filterName, applyColorterms=False, epoch=None)
Definition: jointcal.py:759
lsst::jointcal::SimpleMagnitudeModel
Definition: SimplePhotometryModel.h:123
lsst::jointcal.jointcal.JointcalTask._do_load_refcat_and_fit
def _do_load_refcat_and_fit(self, associations, defaultFilter, center, radius, filters=[], tract="", profile_jointcal=False, match_cut=3.0, reject_bad_fluxes=False, *name="", refObjLoader=None, referenceSelector=None, fit_function=None)
Definition: jointcal.py:670
ast::append
std::shared_ptr< FrameSet > append(FrameSet const &first, FrameSet const &second)
Construct a FrameSet that performs two transformations in series.
Definition: functional.cc:33
lsst.pex.config.history.format
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
Definition: history.py:174
lsst::jointcal.jointcal.JointcalTask.job
job
Definition: jointcal.py:424
lsst::jointcal.jointcal.Photometry
Photometry
Definition: jointcal.py:51
lsst::jointcal.jointcal.JointcalTask._iterate_fit
def _iterate_fit(self, associations, fitter, max_steps, name, whatToFit, dataName="", doRankUpdate=True, doLineSearch=False)
Definition: jointcal.py:1079
lsst::jointcal::ConstrainedMagnitudeModel
Definition: ConstrainedPhotometryModel.h:183
lsst::afw::image::fluxErrFromABMagErr
double fluxErrFromABMagErr(double magErr, double mag) noexcept
Compute flux error in Janskys from AB magnitude error and AB magnitude.
Definition: Calib.h:60
lsst::jointcal.jointcal.JointcalConfig.setDefaults
def setDefaults(self)
Definition: jointcal.py:350
lsst::jointcal::PhotometryFit
Class that handles the photometric least squares problem.
Definition: PhotometryFit.h:45
lsst::jointcal::ConstrainedFluxModel
Definition: ConstrainedPhotometryModel.h:137
lsst::jointcal.jointcal.JointcalTask._DefaultName
string _DefaultName
Definition: jointcal.py:394
lsst::jointcal::Chi2Statistic
Simple structure to accumulate chi2 and ndof.
Definition: Chi2.h:52
lsst::log.log.logContinued.fatal
def fatal(fmt, *args)
Definition: logContinued.py:217
lsst::jointcal.jointcal.JointcalRunner.__call__
def __call__(self, args)
Definition: jointcal.py:91
lsst::jointcal.jointcal.JointcalTask._write_astrometry_results
def _write_astrometry_results(self, associations, model, visit_ccd_to_dataRef)
Definition: jointcal.py:1183
lsst.pipe.tasks.assembleCoadd.run
def run(self, skyInfo, tempExpRefList, imageScalerList, weightList, altMaskList=None, mask=None, supplementaryData=None)
Definition: assembleCoadd.py:720
lsst::jointcal.jointcal.JointcalTask.profile_jointcal
profile_jointcal
Definition: jointcal.py:410
lsst::jointcal.jointcal.JointcalTask._check_star_lists
def _check_star_lists(self, associations, name)
Definition: jointcal.py:815
lsst::jointcal::OneTPPerVisitHandler
A projection handler in which all CCDs from the same visit have the same tangent point.
Definition: ProjectionHandler.h:80
lsst.pex.config
Definition: __init__.py:1
lsst::jointcal::SimpleAstrometryModel
A model where there is one independent transform per CcdImage.
Definition: SimpleAstrometryModel.h:62
lsst::jointcal.jointcal.JointcalTask
Definition: jointcal.py:389
lsst::jointcal.jointcal.JointcalTask._build_ccdImage
def _build_ccdImage(self, dataRef, associations, jointcalControl)
Definition: jointcal.py:442
lsst::log
Definition: Log.h:706
lsst::jointcal.jointcal.JointcalTask.focalPlaneBBox
focalPlaneBBox
Definition: jointcal.py:550
lsst::afw::table
Definition: table.dox:3
lsst::jointcal::SimpleFluxModel
Definition: SimplePhotometryModel.h:86
lsst::utils
Definition: Backtrace.h:29
lsst::jointcal
Definition: Associations.h:49
lsst::jointcal.jointcal.JointcalTask._check_stars
def _check_stars(self, associations)
Definition: jointcal.py:1066
lsst::jointcal.jointcal.JointcalTask._fit_astrometry
def _fit_astrometry(self, associations, dataName=None)
Definition: jointcal.py:970
lsst::jointcal.jointcal.JointcalConfig.colorterms
colorterms
Definition: jointcal.py:241
lsst::jointcal::Associations
The class that implements the relations between MeasuredStar and FittedStar.
Definition: Associations.h:54
lsst::geom
Definition: AffineTransform.h:36
lsst::jointcal.jointcal.JointcalTask.runDataRef
def runDataRef(self, dataRefs, profile_jointcal=False)
Definition: jointcal.py:510
type
table::Key< int > type
Definition: Detector.cc:163
lsst::jointcal::JointcalControl
Definition: JointcalControl.h:39
lsst::geom::Angle
A class representing an angle.
Definition: Angle.h:127
lsst::jointcal.jointcal.JointcalTask._compute_proper_motion_epoch
def _compute_proper_motion_epoch(self, ccdImageList)
Definition: jointcal.py:654
lsst::jointcal.jointcal.JointcalRunner
Definition: jointcal.py:61
lsst::jointcal.jointcal.JointcalTask.photometryRefObjLoader
photometryRefObjLoader
Definition: jointcal.py:421
lsst::jointcal::ConstrainedAstrometryModel
A multi-component model, fitting mappings for sensors and visits simultaneously.
Definition: ConstrainedAstrometryModel.h:62
lsst::jointcal::AstrometryFit
Class that handles the astrometric least squares problem.
Definition: AstrometryFit.h:78
lsst.pex::exceptions
Definition: Exception.h:37
lsst.pex::exceptions::Exception
Provides consistent interface for LSST exceptions.
Definition: Exception.h:107
lsst.pipe.tasks.colorterms
Definition: colorterms.py:1
lsst::geom::SpherePoint
Point in an unspecified spherical coordinate system.
Definition: SpherePoint.h:57
lsst::jointcal.jointcal.JointcalTask._logChi2AndValidate
def _logChi2AndValidate(self, associations, fit, model, chi2Label, writeChi2Name=None)
Definition: jointcal.py:824
lsst::jointcal.jointcal.JointcalConfig.astrometryRefObjLoader
astrometryRefObjLoader
Definition: jointcal.py:277
lsst::jointcal.jointcal.JointcalConfig.doPhotometry
doPhotometry
Definition: jointcal.py:147
lsst::jointcal.jointcal.JointcalConfig.validate
def validate(self)
Definition: jointcal.py:340
lsst::jointcal.jointcal.JointcalConfig.photometryRefObjLoader
photometryRefObjLoader
Definition: jointcal.py:281
lsst.pipe.base
Definition: __init__.py:1
lsst::jointcal.jointcal.writeModel
def writeModel(model, filename, log)
Definition: jointcal.py:382
lsst::meas::algorithms
Fit spatial kernel using approximate fluxes for candidates, and solving a linear system of equations.
Definition: CoaddBoundedField.h:34
lsst::jointcal.jointcal.JointcalConfig.sourceSelector
sourceSelector
Definition: jointcal.py:285
lsst::jointcal.jointcal.JointcalConfig.doAstrometry
doAstrometry
Definition: jointcal.py:142
lsst::jointcal.jointcal.JointcalConfig.applyColorTerms
applyColorTerms
Definition: jointcal.py:235
lsst::jointcal.jointcal.JointcalTask._get_refcat_coordinate_error_override
def _get_refcat_coordinate_error_override(self, refCat, name)
Definition: jointcal.py:607