From c74e0f50b4f98486980febd0cb7d9b9bd16ac2f5 Mon Sep 17 00:00:00 2001 From: John Parejko Date: Mon, 12 May 2025 11:23:46 -0700 Subject: [PATCH 1/7] Cleanup SourceRecord use No need to extract fields, just use the record itself; makes for easier test debugging. --- tests/test_dipoleFitter.py | 45 +++++++++++++++++++------------------- 1 file changed, 22 insertions(+), 23 deletions(-) diff --git a/tests/test_dipoleFitter.py b/tests/test_dipoleFitter.py index 5da747eb4..fca56e855 100644 --- a/tests/test_dipoleFitter.py +++ b/tests/test_dipoleFitter.py @@ -161,44 +161,43 @@ def _checkTaskOutput(self, dipoleTestImage, sources, rtol=None): if rtol is None: rtol = dipoleTestImage.rtol offsets = dipoleTestImage.offsets - for i, r1 in enumerate(sources): - result = r1.extract("ip_diffim_DipoleFit*") - self.assertFloatsAlmostEqual((result['ip_diffim_DipoleFit_pos_instFlux'] - + abs(result['ip_diffim_DipoleFit_neg_instFlux']))/2., + for i, record in enumerate(sources): + self.assertFloatsAlmostEqual((record['ip_diffim_DipoleFit_pos_instFlux'] + + abs(record['ip_diffim_DipoleFit_neg_instFlux']))/2., dipoleTestImage.flux[i], rtol=rtol) - self.assertFloatsAlmostEqual(result['ip_diffim_DipoleFit_pos_x'], + self.assertFloatsAlmostEqual(record['ip_diffim_DipoleFit_pos_x'], dipoleTestImage.xc[i] + offsets[i], rtol=rtol) - self.assertFloatsAlmostEqual(result['ip_diffim_DipoleFit_pos_y'], + self.assertFloatsAlmostEqual(record['ip_diffim_DipoleFit_pos_y'], dipoleTestImage.yc[i] + offsets[i], rtol=rtol) - self.assertFloatsAlmostEqual(result['ip_diffim_DipoleFit_neg_x'], + self.assertFloatsAlmostEqual(record['ip_diffim_DipoleFit_neg_x'], dipoleTestImage.xc[i] - offsets[i], rtol=rtol) - self.assertFloatsAlmostEqual(result['ip_diffim_DipoleFit_neg_y'], + self.assertFloatsAlmostEqual(record['ip_diffim_DipoleFit_neg_y'], dipoleTestImage.yc[i] - offsets[i], rtol=rtol) # Note this is dependent on the noise (variance) being realistic in the image. # otherwise it throws off the chi2 estimate, which is used for classification: - self.assertTrue(result['ip_diffim_DipoleFit_classification']) + self.assertTrue(record['ip_diffim_DipoleFit_classification']) # compare to the original ip_diffim_PsfDipoleFlux measurements - result2 = r1.extract("ip_diffim_PsfDipoleFlux*") - self.assertFloatsAlmostEqual((result['ip_diffim_DipoleFit_pos_instFlux'] - + abs(result['ip_diffim_DipoleFit_neg_instFlux']))/2., - (result2['ip_diffim_PsfDipoleFlux_pos_instFlux'] - + abs(result2['ip_diffim_PsfDipoleFlux_neg_instFlux']))/2., + # result2 = record.extract("ip_diffim_PsfDipoleFlux*") + self.assertFloatsAlmostEqual((record['ip_diffim_DipoleFit_pos_instFlux'] + + abs(record['ip_diffim_DipoleFit_neg_instFlux']))/2., + (record['ip_diffim_PsfDipoleFlux_pos_instFlux'] + + abs(record['ip_diffim_PsfDipoleFlux_neg_instFlux']))/2., rtol=rtol) - self.assertFloatsAlmostEqual(result['ip_diffim_DipoleFit_pos_x'], - result2['ip_diffim_PsfDipoleFlux_pos_centroid_x'], + self.assertFloatsAlmostEqual(record['ip_diffim_DipoleFit_pos_x'], + record['ip_diffim_PsfDipoleFlux_pos_centroid_x'], rtol=rtol) - self.assertFloatsAlmostEqual(result['ip_diffim_DipoleFit_pos_y'], - result2['ip_diffim_PsfDipoleFlux_pos_centroid_y'], + self.assertFloatsAlmostEqual(record['ip_diffim_DipoleFit_pos_y'], + record['ip_diffim_PsfDipoleFlux_pos_centroid_y'], rtol=rtol) - self.assertFloatsAlmostEqual(result['ip_diffim_DipoleFit_neg_x'], - result2['ip_diffim_PsfDipoleFlux_neg_centroid_x'], + self.assertFloatsAlmostEqual(record['ip_diffim_DipoleFit_neg_x'], + record['ip_diffim_PsfDipoleFlux_neg_centroid_x'], rtol=rtol) - self.assertFloatsAlmostEqual(result['ip_diffim_DipoleFit_neg_y'], - result2['ip_diffim_PsfDipoleFlux_neg_centroid_y'], + self.assertFloatsAlmostEqual(record['ip_diffim_DipoleFit_neg_y'], + record['ip_diffim_PsfDipoleFlux_neg_centroid_y'], rtol=rtol) - return result + return record.extract("ip_diffim_DipoleFit*") def testDipoleTask(self): """Test the dipole fitting singleFramePlugin. From 3f0d1d82c2ae1542a4e1c3e99941ffc19092d433 Mon Sep 17 00:00:00 2001 From: John Parejko Date: Mon, 12 May 2025 11:34:26 -0700 Subject: [PATCH 2/7] Add dipole centroid uncertainties, with basic test * Cleanup the return struct from `fitDipole` to use CentroidResult, so that we can just put that in the source record directly. * I'm not keen on how I'm using the SDSS Centroid uncertainties as the test expected values, but I don't know a better approach. * Not including the correlations/covariances, as that math is trickier. --- python/lsst/ip/diffim/dipoleFitTask.py | 36 ++++++++++------ tests/test_dipoleFitter.py | 59 +++++++++++++++++++++----- 2 files changed, 71 insertions(+), 24 deletions(-) diff --git a/python/lsst/ip/diffim/dipoleFitTask.py b/python/lsst/ip/diffim/dipoleFitTask.py index 666ae0627..8d59c7ec8 100644 --- a/python/lsst/ip/diffim/dipoleFitTask.py +++ b/python/lsst/ip/diffim/dipoleFitTask.py @@ -20,6 +20,7 @@ # see . # +import math import logging import numpy as np import warnings @@ -863,8 +864,18 @@ def fitDipole(self, source, tol=1e-7, rel_weight=0.1, if fitParams['flux'] <= 1.: # usually around 0.1 -- the minimum flux allowed -- i.e. bad fit. return None, fitResult - centroid = ((fitParams['xcenPos'] + fitParams['xcenNeg']) / 2., - (fitParams['ycenPos'] + fitParams['ycenNeg']) / 2.) + # TODO: We could include covariances, which could be derived from + # `fitResult.params[name].correl`, but those are correlations. + posCentroid = measBase.CentroidResult(fitParams['xcenPos'], fitParams['ycenPos'], + fitResult.params['xcenPos'].stderr, + fitResult.params['ycenPos'].stderr) + negCentroid = measBase.CentroidResult(fitParams['xcenNeg'], fitParams['ycenNeg'], + fitResult.params['xcenNeg'].stderr, + fitResult.params['ycenNeg'].stderr) + centroid = measBase.CentroidResult((fitParams['xcenPos'] + fitParams['xcenNeg']) / 2, + (fitParams['ycenPos'] + fitParams['ycenNeg']) / 2., + math.sqrt(posCentroid.xErr**2 + negCentroid.xErr**2), + math.sqrt(posCentroid.yErr**2 + negCentroid.yErr**2)) dx, dy = fitParams['xcenPos'] - fitParams['xcenNeg'], fitParams['ycenPos'] - fitParams['ycenNeg'] angle = np.arctan2(dy, dx) @@ -892,10 +903,9 @@ def computeSumVariance(exposure, footprint): except ZeroDivisionError: # catch divide by zero - should never happen. signalToNoise = np.nan - out = Struct(posCentroidX=fitParams['xcenPos'], posCentroidY=fitParams['ycenPos'], - negCentroidX=fitParams['xcenNeg'], negCentroidY=fitParams['ycenNeg'], + out = Struct(posCentroid=posCentroid, negCentroid=negCentroid, centroid=centroid, posFlux=fluxVal, negFlux=-fluxValNeg, posFluxErr=fluxErr, negFluxErr=fluxErrNeg, - centroidX=centroid[0], centroidY=centroid[1], orientation=angle, + orientation=angle, signalToNoise=signalToNoise, chi2=fitResult.chisqr, redChi2=fitResult.redchi, nData=fitResult.ndata) @@ -1010,11 +1020,11 @@ def _setupSchema(self, config, name, schema, metadata): self.posCentroidKey = measBase.CentroidResultKey.addFields(schema, schema.join(name, "pos"), "Dipole positive lobe centroid position.", - measBase.UncertaintyEnum.NO_UNCERTAINTY) + measBase.UncertaintyEnum.SIGMA_ONLY) self.negCentroidKey = measBase.CentroidResultKey.addFields(schema, schema.join(name, "neg"), "Dipole negative lobe centroid position.", - measBase.UncertaintyEnum.NO_UNCERTAINTY) + measBase.UncertaintyEnum.SIGMA_ONLY) self.centroidKey = measBase.CentroidResultKey.addFields(schema, name, "Dipole centroid position.", @@ -1133,21 +1143,21 @@ def measureDipoles(self, measRecord, exposure, posExp=None, negExp=None): # Add the relevant values to the measRecord measRecord[self.posFluxKey.getInstFlux()] = result.posFlux measRecord[self.posFluxKey.getInstFluxErr()] = result.signalToNoise # to be changed to actual sigma! - measRecord[self.posCentroidKey.getX()] = result.posCentroidX + self.posCentroidKey.set(measRecord, result.posCentroid) measRecord[self.posCentroidKey.getY()] = result.posCentroidY measRecord[self.negFluxKey.getInstFlux()] = result.negFlux measRecord[self.negFluxKey.getInstFluxErr()] = result.signalToNoise # to be changed to actual sigma! - measRecord[self.negCentroidKey.getX()] = result.negCentroidX + self.negCentroidKey.set(measRecord, result.negCentroid) measRecord[self.negCentroidKey.getY()] = result.negCentroidY # Dia source flux: average of pos+neg measRecord[self.fluxKey.getInstFlux()] = (abs(result.posFlux) + abs(result.negFlux))/2. measRecord[self.orientationKey] = result.orientation - measRecord[self.separationKey] = np.sqrt((result.posCentroidX - result.negCentroidX)**2. - + (result.posCentroidY - result.negCentroidY)**2.) - measRecord[self.centroidKey.getX()] = result.centroidX - measRecord[self.centroidKey.getY()] = result.centroidY + measRecord[self.separationKey] = np.sqrt((result.posCentroid.x - result.negCentroid.x)**2 + + (result.posCentroid.y - result.negCentroid.y)**2) + + self.centroidKey.set(measRecord, result.centroid) measRecord[self.signalToNoiseKey] = result.signalToNoise measRecord[self.chi2dofKey] = result.redChi2 diff --git a/tests/test_dipoleFitter.py b/tests/test_dipoleFitter.py index fca56e855..a2baf27a8 100644 --- a/tests/test_dipoleFitter.py +++ b/tests/test_dipoleFitter.py @@ -115,10 +115,10 @@ def testDipoleAlgorithm(self): self.assertFloatsAlmostEqual((result.posFlux + abs(result.negFlux))/2., dipoleTestImage.flux[i], rtol=rtol) - self.assertFloatsAlmostEqual(result.posCentroidX, dipoleTestImage.xc[i] + offsets[i], rtol=rtol) - self.assertFloatsAlmostEqual(result.posCentroidY, dipoleTestImage.yc[i] + offsets[i], rtol=rtol) - self.assertFloatsAlmostEqual(result.negCentroidX, dipoleTestImage.xc[i] - offsets[i], rtol=rtol) - self.assertFloatsAlmostEqual(result.negCentroidY, dipoleTestImage.yc[i] - offsets[i], rtol=rtol) + self.assertFloatsAlmostEqual(result.posCentroid.x, dipoleTestImage.xc[i] + offsets[i], rtol=rtol) + self.assertFloatsAlmostEqual(result.posCentroid.y, dipoleTestImage.yc[i] + offsets[i], rtol=rtol) + self.assertFloatsAlmostEqual(result.negCentroid.x, dipoleTestImage.xc[i] - offsets[i], rtol=rtol) + self.assertFloatsAlmostEqual(result.negCentroid.y, dipoleTestImage.yc[i] - offsets[i], rtol=rtol) def _runDetection(self, dipoleTestImage, maxFootprintArea=None): """Run 'diaSource' detection on the diffim, including merging of @@ -147,7 +147,7 @@ def _runDetection(self, dipoleTestImage, maxFootprintArea=None): measureTask.run(sources, testImage.diffim, testImage.posImage, testImage.negImage) return sources - def _checkTaskOutput(self, dipoleTestImage, sources, rtol=None): + def _checkTaskOutput(self, dipoleTestImage, sources, noPreImages=False): """Compare the fluxes/centroids in `sources` are entered into the correct slots of the catalog, and have values that are very close to the input values for both dipoles in the @@ -156,10 +156,18 @@ def _checkTaskOutput(self, dipoleTestImage, sources, rtol=None): Also test that the resulting fluxes are close to those generated by the existing ip_diffim_DipoleMeasurement task (PsfDipoleFit). - """ - if rtol is None: - rtol = dipoleTestImage.rtol + Parameters + ---------- + dipoleTestImage : `DipoleTestImage` + Image that was generated for the test. + sources : `lsst.afw.table.SourceCatalog` + Source catalog measured on the test image. + noPreImages : bool, optional + Set to True if there were no positive/negative images generated + for the test, to use much looser uncertainty tolerances. + """ + rtol = dipoleTestImage.rtol offsets = dipoleTestImage.offsets for i, record in enumerate(sources): self.assertFloatsAlmostEqual((record['ip_diffim_DipoleFit_pos_instFlux'] @@ -173,6 +181,37 @@ def _checkTaskOutput(self, dipoleTestImage, sources, rtol=None): dipoleTestImage.xc[i] - offsets[i], rtol=rtol) self.assertFloatsAlmostEqual(record['ip_diffim_DipoleFit_neg_y'], dipoleTestImage.yc[i] - offsets[i], rtol=rtol) + + # check uncertainties + # TODO: how to compute the correct uncertainty to test here? + # self.assertFloatsAlmostEqual(record['ip_diffim_DipoleFit_instFluxErr'], + # ???, rtol=rtol) + # emperical uncertainty values: not clear how to compute proper expected values. + with self.subTest(i=i, type="pos_xErr"): + self.assertFloatsAlmostEqual(record['ip_diffim_DipoleFit_pos_xErr'], + .53*record["base_SdssCentroid_xErr"], + rtol=0.1 if not noPreImages else 0.5) + with self.subTest(i=i, type="pos_yErr"): + self.assertFloatsAlmostEqual(record['ip_diffim_DipoleFit_pos_yErr'], + .53*record["base_SdssCentroid_yErr"], + rtol=0.1 if not noPreImages else 0.5) + with self.subTest(i=i, type="neg_xErr"): + self.assertFloatsAlmostEqual(record['ip_diffim_DipoleFit_neg_xErr'], + .53*record["base_SdssCentroid_xErr"], + rtol=0.1 if not noPreImages else 0.5) + with self.subTest(i=i, type="neg_yErr"): + self.assertFloatsAlmostEqual(record['ip_diffim_DipoleFit_neg_yErr'], + .53*record["base_SdssCentroid_yErr"], + rtol=0.1 if not noPreImages else 0.5) + with self.subTest(i=i, type="xErr"): + self.assertFloatsAlmostEqual(record['ip_diffim_DipoleFit_xErr'], + .74*record["base_SdssCentroid_xErr"], + rtol=0.06 if not noPreImages else 0.5) + with self.subTest(i=i, type="yErr"): + self.assertFloatsAlmostEqual(record['ip_diffim_DipoleFit_yErr'], + .75*record["base_SdssCentroid_yErr"], + rtol=0.06 if not noPreImages else 0.5) + # Note this is dependent on the noise (variance) being realistic in the image. # otherwise it throws off the chi2 estimate, which is used for classification: self.assertTrue(record['ip_diffim_DipoleFit_classification']) @@ -197,8 +236,6 @@ def _checkTaskOutput(self, dipoleTestImage, sources, rtol=None): record['ip_diffim_PsfDipoleFlux_neg_centroid_y'], rtol=rtol) - return record.extract("ip_diffim_DipoleFit*") - def testDipoleTask(self): """Test the dipole fitting singleFramePlugin. @@ -268,7 +305,7 @@ def testDipoleTaskNoPreSubImages(self): dipoleTestImage = DipoleTestImage() dipoleTestImage.testImage.posImage = dipoleTestImage.testImage.negImage = None sources = self._runDetection(dipoleTestImage) - self._checkTaskOutput(dipoleTestImage, sources) + self._checkTaskOutput(dipoleTestImage, sources, noPreImages=True) def testDipoleEdge(self): """Test the too-close-to-image-edge scenario for dipole fitting From 241d72961ae6820ba69d0daee24dbc385308cad0 Mon Sep 17 00:00:00 2001 From: John Parejko Date: Mon, 12 May 2025 12:02:35 -0700 Subject: [PATCH 3/7] Refactor displayFitResults Add centroids and uncertainties. Use plt.subplots to cleanup display2dArray with better labels. --- python/lsst/ip/diffim/dipoleFitTask.py | 57 +++++++++++++++++--------- 1 file changed, 38 insertions(+), 19 deletions(-) diff --git a/python/lsst/ip/diffim/dipoleFitTask.py b/python/lsst/ip/diffim/dipoleFitTask.py index 8d59c7ec8..09f6b133a 100644 --- a/python/lsst/ip/diffim/dipoleFitTask.py +++ b/python/lsst/ip/diffim/dipoleFitTask.py @@ -928,36 +928,55 @@ def displayFitResults(self, footprint, result): self.log.warning('Unable to import matplotlib: %s', err) raise err - def display2dArray(arr, title='Data', extent=None): + def display2dArray(ax, arr, x, y, xErr, yErr, title, extent=None): """Use `matplotlib.pyplot.imshow` to display a 2-D array with a given coordinate range. """ - fig = plt.imshow(arr, origin='lower', interpolation='none', cmap='gray', extent=extent) - plt.title(title) - plt.colorbar(fig, cmap='gray') + fig = ax.imshow(arr, origin='lower', interpolation='none', cmap='gray', extent=extent) + ax.set_title(title) + ax.errorbar(x["total"], y["total"], xErr["total"], yErr["total"], c="cyan") + ax.errorbar(x["Pos"], y["Pos"], xErr["Pos"], yErr["Pos"], c="green") + ax.errorbar(x["Neg"], y["Neg"], xErr["Neg"], yErr["Neg"], c="red") return fig z = result.data fit = result.best_fit bbox = footprint.getBBox() extent = (bbox.getBeginX(), bbox.getEndX(), bbox.getBeginY(), bbox.getEndY()) + if z.shape[0] == 3: - plt.figure(figsize=(8, 8)) - for i in range(3): - plt.subplot(3, 3, i*3+1) - display2dArray(z[i, :], 'Data', extent=extent) - plt.subplot(3, 3, i*3+2) - display2dArray(fit[i, :], 'Model', extent=extent) - plt.subplot(3, 3, i*3+3) - display2dArray(z[i, :] - fit[i, :], 'Residual', extent=extent) + x, y, xErr, yErr = {}, {}, {}, {} + for name in ("Pos", "Neg"): + x[name] = result.best_values[f"xcen{name}"] + y[name] = result.best_values[f"ycen{name}"] + xErr[name] = result.params[f"xcen{name}"].stderr + yErr[name] = result.params[f"ycen{name}"].stderr + x["total"] = (result.best_values["xcenPos"] + result.best_values["xcenNeg"])/2 + y["total"] = (result.best_values["ycenPos"] + result.best_values["ycenNeg"])/2 + xErr["total"] = math.sqrt(xErr["Pos"]**2 + xErr["Neg"]**2) + yErr["total"] = math.sqrt(yErr["Pos"]**2 + yErr["Neg"]**2) + + fig, axes = plt.subplots(nrows=3, ncols=3, sharex=True, sharey=True, figsize=(8, 8)) + for i, label in enumerate(("total", "Pos", "Neg")): + display2dArray(axes[i][0], z[i, :], x, y, xErr, yErr, + f'Data {label}', extent=extent) + display2dArray(axes[i][1], fit[i, :], x, y, xErr, yErr, + f'Model {label}', extent=extent) + display2dArray(axes[i][2], z[i, :] - fit[i, :], x, y, xErr, yErr, + f'Residual {label}', extent=extent) + + plt.setp(axes[i][1].get_yticklabels(), visible=False) + plt.setp(axes[i][2].get_yticklabels(), visible=False) + if i != 2: # remove top two row x-axis labels + plt.setp(axes[i][0].get_xticklabels(), visible=False) + plt.setp(axes[i][1].get_xticklabels(), visible=False) + plt.setp(axes[i][2].get_xticklabels(), visible=False) else: - plt.figure(figsize=(8, 2.5)) - plt.subplot(1, 3, 1) - display2dArray(z, 'Data', extent=extent) - plt.subplot(1, 3, 2) - display2dArray(fit, 'Model', extent=extent) - plt.subplot(1, 3, 3) - display2dArray(z - fit, 'Residual', extent=extent) + fig, axes = plt.subplots(nrows=3, ncols=3, sharex=True, sharey=True, figsize=(10, 2.5)) + display2dArray(axes[0], z, 'Data', extent=extent) + display2dArray(axes[1], 'Model', extent=extent) + display2dArray(axes[2], z - fit, 'Residual', extent=extent) + fig.tight_layout(pad=0, w_pad=0, h_pad=0) plt.show() From 1c642678e501adf0194b66b9d2986d1888c58302 Mon Sep 17 00:00:00 2001 From: John Parejko Date: Mon, 12 May 2025 12:06:17 -0700 Subject: [PATCH 4/7] Use math.sqrt for single values: it's much faster --- python/lsst/ip/diffim/dipoleFitTask.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/python/lsst/ip/diffim/dipoleFitTask.py b/python/lsst/ip/diffim/dipoleFitTask.py index 09f6b133a..cd54efb54 100644 --- a/python/lsst/ip/diffim/dipoleFitTask.py +++ b/python/lsst/ip/diffim/dipoleFitTask.py @@ -79,7 +79,7 @@ class DipoleFitPluginConfig(measBase.SingleFramePluginConfig): # Config params for classification of detected diaSources as dipole or not minSn = pexConfig.Field( - dtype=float, default=np.sqrt(2) * 5.0, + dtype=float, default=math.sqrt(2) * 5.0, doc="Minimum quadrature sum of positive+negative lobe S/N to be considered a dipole") maxFluxRatio = pexConfig.Field( @@ -882,7 +882,7 @@ def fitDipole(self, source, tol=1e-7, rel_weight=0.1, # Exctract flux value, compute signalToNoise from flux/variance_within_footprint # Also extract the stderr of flux estimate. def computeSumVariance(exposure, footprint): - return np.sqrt(np.nansum(exposure[footprint.getBBox(), afwImage.PARENT].variance.array)) + return math.sqrt(np.nansum(exposure[footprint.getBBox(), afwImage.PARENT].variance.array)) fluxVal = fluxVar = fitParams['flux'] fluxErr = fluxErrNeg = fitResult.params['flux'].stderr @@ -899,7 +899,7 @@ def computeSumVariance(exposure, footprint): fluxVarNeg = computeSumVariance(self.negImage, source.getFootprint()) try: - signalToNoise = np.sqrt((fluxVal/fluxVar)**2 + (fluxValNeg/fluxVarNeg)**2) + signalToNoise = math.sqrt((fluxVal/fluxVar)**2 + (fluxValNeg/fluxVarNeg)**2) except ZeroDivisionError: # catch divide by zero - should never happen. signalToNoise = np.nan @@ -1173,8 +1173,8 @@ def measureDipoles(self, measRecord, exposure, posExp=None, negExp=None): # Dia source flux: average of pos+neg measRecord[self.fluxKey.getInstFlux()] = (abs(result.posFlux) + abs(result.negFlux))/2. measRecord[self.orientationKey] = result.orientation - measRecord[self.separationKey] = np.sqrt((result.posCentroid.x - result.negCentroid.x)**2 - + (result.posCentroid.y - result.negCentroid.y)**2) + measRecord[self.separationKey] = math.sqrt((result.posCentroid.x - result.negCentroid.x)**2 + + (result.posCentroid.y - result.negCentroid.y)**2) self.centroidKey.set(measRecord, result.centroid) From 70a4d7aa13f1b4e7f68a677eeca474897a0583ec Mon Sep 17 00:00:00 2001 From: John Parejko Date: Mon, 12 May 2025 12:37:23 -0700 Subject: [PATCH 5/7] Refactor to remove deprecated `best_values` lmfit now prefers `params` in the MinimizerResult, which include both the best fit `value` and the `stderr` (and other things). Switching to this makes it less confusing, since I needed to use `fitResult.params` in addition to the existing `fitParams`, which were not at all the same thing! --- python/lsst/ip/diffim/dipoleFitTask.py | 20 ++++++++++++-------- 1 file changed, 12 insertions(+), 8 deletions(-) diff --git a/python/lsst/ip/diffim/dipoleFitTask.py b/python/lsst/ip/diffim/dipoleFitTask.py index cd54efb54..41217f8b4 100644 --- a/python/lsst/ip/diffim/dipoleFitTask.py +++ b/python/lsst/ip/diffim/dipoleFitTask.py @@ -860,23 +860,26 @@ def fitDipole(self, source, tol=1e-7, rel_weight=0.1, fp = source.getFootprint() self.displayFitResults(fp, fitResult) - fitParams = fitResult.best_values - if fitParams['flux'] <= 1.: # usually around 0.1 -- the minimum flux allowed -- i.e. bad fit. + # usually around 0.1 -- the minimum flux allowed -- i.e. bad fit. + if fitResult.params['flux'].value <= 1.: return None, fitResult # TODO: We could include covariances, which could be derived from # `fitResult.params[name].correl`, but those are correlations. - posCentroid = measBase.CentroidResult(fitParams['xcenPos'], fitParams['ycenPos'], + posCentroid = measBase.CentroidResult(fitResult.params['xcenPos'].value, + fitResult.params['ycenPos'].value, fitResult.params['xcenPos'].stderr, fitResult.params['ycenPos'].stderr) - negCentroid = measBase.CentroidResult(fitParams['xcenNeg'], fitParams['ycenNeg'], + negCentroid = measBase.CentroidResult(fitResult.params['xcenNeg'].value, + fitResult.params['ycenNeg'].value, fitResult.params['xcenNeg'].stderr, fitResult.params['ycenNeg'].stderr) - centroid = measBase.CentroidResult((fitParams['xcenPos'] + fitParams['xcenNeg']) / 2, - (fitParams['ycenPos'] + fitParams['ycenNeg']) / 2., + centroid = measBase.CentroidResult((fitResult.params['xcenPos'] + fitResult.params['xcenNeg']) / 2, + (fitResult.params['ycenPos'] + fitResult.params['ycenNeg']) / 2., math.sqrt(posCentroid.xErr**2 + negCentroid.xErr**2), math.sqrt(posCentroid.yErr**2 + negCentroid.yErr**2)) - dx, dy = fitParams['xcenPos'] - fitParams['xcenNeg'], fitParams['ycenPos'] - fitParams['ycenNeg'] + dx = fitResult.params['xcenPos'].value - fitResult.params['xcenNeg'].value + dy = fitResult.params['ycenPos'].value - fitResult.params['ycenNeg'].value angle = np.arctan2(dy, dx) # Exctract flux value, compute signalToNoise from flux/variance_within_footprint @@ -885,6 +888,7 @@ def computeSumVariance(exposure, footprint): return math.sqrt(np.nansum(exposure[footprint.getBBox(), afwImage.PARENT].variance.array)) fluxVal = fluxVar = fitParams['flux'] + fluxVal = fluxVar = fitResult.params['flux'].value fluxErr = fluxErrNeg = fitResult.params['flux'].stderr if self.posImage is not None: fluxVar = computeSumVariance(self.posImage, source.getFootprint()) @@ -893,7 +897,7 @@ def computeSumVariance(exposure, footprint): fluxValNeg, fluxVarNeg = fluxVal, fluxVar if separateNegParams: - fluxValNeg = fitParams['fluxNeg'] + fluxValNeg = fitResult.params['fluxNeg'].value fluxErrNeg = fitResult.params['fluxNeg'].stderr if self.negImage is not None: fluxVarNeg = computeSumVariance(self.negImage, source.getFootprint()) From 35aadda1df8744d239bf203fb2203b536c57f773 Mon Sep 17 00:00:00 2001 From: John Parejko Date: Mon, 12 May 2025 13:05:44 -0700 Subject: [PATCH 6/7] Refactor flux handling to use FluxResult This makes it more explicit that the positive/negative/total flux values are going to be the same unless separateNegParams is set, which it is not by default from the config. --- python/lsst/ip/diffim/dipoleFitTask.py | 38 ++++++++++++-------------- tests/test_dipoleFitter.py | 2 +- 2 files changed, 18 insertions(+), 22 deletions(-) diff --git a/python/lsst/ip/diffim/dipoleFitTask.py b/python/lsst/ip/diffim/dipoleFitTask.py index 41217f8b4..de5762324 100644 --- a/python/lsst/ip/diffim/dipoleFitTask.py +++ b/python/lsst/ip/diffim/dipoleFitTask.py @@ -884,32 +884,34 @@ def fitDipole(self, source, tol=1e-7, rel_weight=0.1, # Exctract flux value, compute signalToNoise from flux/variance_within_footprint # Also extract the stderr of flux estimate. + # TODO: should this instead use the lmfit-computed uncertainty from + # `lmfitResult.result.uvars['flux'].std_dev`? def computeSumVariance(exposure, footprint): return math.sqrt(np.nansum(exposure[footprint.getBBox(), afwImage.PARENT].variance.array)) - fluxVal = fluxVar = fitParams['flux'] - fluxVal = fluxVar = fitResult.params['flux'].value - fluxErr = fluxErrNeg = fitResult.params['flux'].stderr + # NOTE: These will all be the same unless separateNegParams=True! + flux = measBase.FluxResult(fitResult.params["flux"].value, fitResult.params["flux"].stderr) + posFlux = measBase.FluxResult(fitResult.params["flux"].value, fitResult.params["flux"].stderr) + negFlux = measBase.FluxResult(fitResult.params["flux"].value, fitResult.params["flux"].stderr) if self.posImage is not None: fluxVar = computeSumVariance(self.posImage, source.getFootprint()) else: fluxVar = computeSumVariance(self.diffim, source.getFootprint()) + fluxVarNeg = fluxVar - fluxValNeg, fluxVarNeg = fluxVal, fluxVar if separateNegParams: - fluxValNeg = fitResult.params['fluxNeg'].value - fluxErrNeg = fitResult.params['fluxNeg'].stderr + negFlux.instFlux = fitResult.params['fluxNeg'].value + negFlux.instFluxErr = fitResult.params['fluxNeg'].stderr if self.negImage is not None: fluxVarNeg = computeSumVariance(self.negImage, source.getFootprint()) try: - signalToNoise = math.sqrt((fluxVal/fluxVar)**2 + (fluxValNeg/fluxVarNeg)**2) + signalToNoise = math.sqrt((posFlux.instFlux/fluxVar)**2 + (negFlux.instFlux/fluxVarNeg)**2) except ZeroDivisionError: # catch divide by zero - should never happen. signalToNoise = np.nan out = Struct(posCentroid=posCentroid, negCentroid=negCentroid, centroid=centroid, - posFlux=fluxVal, negFlux=-fluxValNeg, posFluxErr=fluxErr, negFluxErr=fluxErrNeg, - orientation=angle, + posFlux=posFlux, negFlux=negFlux, flux=flux, orientation=angle, signalToNoise=signalToNoise, chi2=fitResult.chisqr, redChi2=fitResult.redchi, nData=fitResult.ndata) @@ -1162,26 +1164,20 @@ def measureDipoles(self, measRecord, exposure, posExp=None, negExp=None): self.fail(measRecord, measBase.MeasurementError("bad dipole fit", self.FAILURE_FIT)) return - # add chi2, coord/flux uncertainties (TBD), dipole classification - # Add the relevant values to the measRecord - measRecord[self.posFluxKey.getInstFlux()] = result.posFlux - measRecord[self.posFluxKey.getInstFluxErr()] = result.signalToNoise # to be changed to actual sigma! + # TODO: add chi2, dipole classification + self.posFluxKey.set(measRecord, result.posFlux) self.posCentroidKey.set(measRecord, result.posCentroid) - measRecord[self.posCentroidKey.getY()] = result.posCentroidY - measRecord[self.negFluxKey.getInstFlux()] = result.negFlux - measRecord[self.negFluxKey.getInstFluxErr()] = result.signalToNoise # to be changed to actual sigma! + self.negFluxKey.set(measRecord, result.negFlux) self.negCentroidKey.set(measRecord, result.negCentroid) - measRecord[self.negCentroidKey.getY()] = result.negCentroidY - # Dia source flux: average of pos+neg - measRecord[self.fluxKey.getInstFlux()] = (abs(result.posFlux) + abs(result.negFlux))/2. + self.fluxKey.set(measRecord, result.flux) + self.centroidKey.set(measRecord, result.centroid) + measRecord[self.orientationKey] = result.orientation measRecord[self.separationKey] = math.sqrt((result.posCentroid.x - result.negCentroid.x)**2 + (result.posCentroid.y - result.negCentroid.y)**2) - self.centroidKey.set(measRecord, result.centroid) - measRecord[self.signalToNoiseKey] = result.signalToNoise measRecord[self.chi2dofKey] = result.redChi2 diff --git a/tests/test_dipoleFitter.py b/tests/test_dipoleFitter.py index a2baf27a8..138be2cd0 100644 --- a/tests/test_dipoleFitter.py +++ b/tests/test_dipoleFitter.py @@ -113,7 +113,7 @@ def testDipoleAlgorithm(self): s, rel_weight=0.5, separateNegParams=False, verbose=verbose, display=display) - self.assertFloatsAlmostEqual((result.posFlux + abs(result.negFlux))/2., + self.assertFloatsAlmostEqual((result.posFlux.instFlux + abs(result.negFlux.instFlux))/2., dipoleTestImage.flux[i], rtol=rtol) self.assertFloatsAlmostEqual(result.posCentroid.x, dipoleTestImage.xc[i] + offsets[i], rtol=rtol) self.assertFloatsAlmostEqual(result.posCentroid.y, dipoleTestImage.yc[i] + offsets[i], rtol=rtol) From 6d86b0c8c16c8d778b1040e98d1659e0ad898e2c Mon Sep 17 00:00:00 2001 From: John Parejko Date: Tue, 20 May 2025 14:38:07 -0700 Subject: [PATCH 7/7] Include positive/negative covariance in total centroid error --- python/lsst/ip/diffim/dipoleFitTask.py | 12 +++++++++--- tests/test_dipoleFitter.py | 6 ++++-- 2 files changed, 13 insertions(+), 5 deletions(-) diff --git a/python/lsst/ip/diffim/dipoleFitTask.py b/python/lsst/ip/diffim/dipoleFitTask.py index de5762324..45444a606 100644 --- a/python/lsst/ip/diffim/dipoleFitTask.py +++ b/python/lsst/ip/diffim/dipoleFitTask.py @@ -874,15 +874,21 @@ def fitDipole(self, source, tol=1e-7, rel_weight=0.1, fitResult.params['ycenNeg'].value, fitResult.params['xcenNeg'].stderr, fitResult.params['ycenNeg'].stderr) + xposIdx = fitResult.var_names.index("xcenPos") + yposIdx = fitResult.var_names.index("ycenPos") + xnegIdx = fitResult.var_names.index("xcenNeg") + ynegIdx = fitResult.var_names.index("ycenNeg") centroid = measBase.CentroidResult((fitResult.params['xcenPos'] + fitResult.params['xcenNeg']) / 2, (fitResult.params['ycenPos'] + fitResult.params['ycenNeg']) / 2., - math.sqrt(posCentroid.xErr**2 + negCentroid.xErr**2), - math.sqrt(posCentroid.yErr**2 + negCentroid.yErr**2)) + math.sqrt(posCentroid.xErr**2 + negCentroid.xErr**2 + + 2*fitResult.covar[xposIdx, xnegIdx]) / 2., + math.sqrt(posCentroid.yErr**2 + negCentroid.yErr**2 + + 2*fitResult.covar[yposIdx, ynegIdx]) / 2.) dx = fitResult.params['xcenPos'].value - fitResult.params['xcenNeg'].value dy = fitResult.params['ycenPos'].value - fitResult.params['ycenNeg'].value angle = np.arctan2(dy, dx) - # Exctract flux value, compute signalToNoise from flux/variance_within_footprint + # Extract flux value, compute signalToNoise from flux/variance_within_footprint # Also extract the stderr of flux estimate. # TODO: should this instead use the lmfit-computed uncertainty from # `lmfitResult.result.uvars['flux'].std_dev`? diff --git a/tests/test_dipoleFitter.py b/tests/test_dipoleFitter.py index 138be2cd0..09611b13e 100644 --- a/tests/test_dipoleFitter.py +++ b/tests/test_dipoleFitter.py @@ -203,13 +203,15 @@ def _checkTaskOutput(self, dipoleTestImage, sources, noPreImages=False): self.assertFloatsAlmostEqual(record['ip_diffim_DipoleFit_neg_yErr'], .53*record["base_SdssCentroid_yErr"], rtol=0.1 if not noPreImages else 0.5) + # NOTE: these are smaller than the positive/negative uncertainties, + # because the of those covariance is negative! with self.subTest(i=i, type="xErr"): self.assertFloatsAlmostEqual(record['ip_diffim_DipoleFit_xErr'], - .74*record["base_SdssCentroid_xErr"], + .34*record["base_SdssCentroid_xErr"], rtol=0.06 if not noPreImages else 0.5) with self.subTest(i=i, type="yErr"): self.assertFloatsAlmostEqual(record['ip_diffim_DipoleFit_yErr'], - .75*record["base_SdssCentroid_yErr"], + .35*record["base_SdssCentroid_yErr"], rtol=0.06 if not noPreImages else 0.5) # Note this is dependent on the noise (variance) being realistic in the image.