diff --git a/python/lsst/ip/diffim/dipoleFitTask.py b/python/lsst/ip/diffim/dipoleFitTask.py index 666ae0627..45444a606 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 @@ -78,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( @@ -859,43 +860,64 @@ 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 - centroid = ((fitParams['xcenPos'] + fitParams['xcenNeg']) / 2., - (fitParams['ycenPos'] + fitParams['ycenNeg']) / 2.) - dx, dy = fitParams['xcenPos'] - fitParams['xcenNeg'], fitParams['ycenPos'] - fitParams['ycenNeg'] + # TODO: We could include covariances, which could be derived from + # `fitResult.params[name].correl`, but those are correlations. + posCentroid = measBase.CentroidResult(fitResult.params['xcenPos'].value, + fitResult.params['ycenPos'].value, + fitResult.params['xcenPos'].stderr, + fitResult.params['ycenPos'].stderr) + negCentroid = measBase.CentroidResult(fitResult.params['xcenNeg'].value, + 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 + + 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`? 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 + # 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 = fitParams['fluxNeg'] - 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 = np.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(posCentroidX=fitParams['xcenPos'], posCentroidY=fitParams['ycenPos'], - negCentroidX=fitParams['xcenNeg'], negCentroidY=fitParams['ycenNeg'], - posFlux=fluxVal, negFlux=-fluxValNeg, posFluxErr=fluxErr, negFluxErr=fluxErrNeg, - centroidX=centroid[0], centroidY=centroid[1], orientation=angle, + out = Struct(posCentroid=posCentroid, negCentroid=negCentroid, centroid=centroid, + posFlux=posFlux, negFlux=negFlux, flux=flux, orientation=angle, signalToNoise=signalToNoise, chi2=fitResult.chisqr, redChi2=fitResult.redchi, nData=fitResult.ndata) @@ -918,36 +940,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() @@ -1010,11 +1051,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.", @@ -1129,25 +1170,19 @@ 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! - measRecord[self.posCentroidKey.getX()] = result.posCentroidX - measRecord[self.posCentroidKey.getY()] = result.posCentroidY + # TODO: add chi2, dipole classification + self.posFluxKey.set(measRecord, result.posFlux) + self.posCentroidKey.set(measRecord, result.posCentroid) + + self.negFluxKey.set(measRecord, result.negFlux) + self.negCentroidKey.set(measRecord, result.negCentroid) - measRecord[self.negFluxKey.getInstFlux()] = result.negFlux - measRecord[self.negFluxKey.getInstFluxErr()] = result.signalToNoise # to be changed to actual sigma! - measRecord[self.negCentroidKey.getX()] = result.negCentroidX - measRecord[self.negCentroidKey.getY()] = result.negCentroidY + self.fluxKey.set(measRecord, result.flux) + self.centroidKey.set(measRecord, result.centroid) - # 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] = math.sqrt((result.posCentroid.x - result.negCentroid.x)**2 + + (result.posCentroid.y - result.negCentroid.y)**2) measRecord[self.signalToNoiseKey] = result.signalToNoise measRecord[self.chi2dofKey] = result.redChi2 diff --git a/tests/test_dipoleFitter.py b/tests/test_dipoleFitter.py index 5da747eb4..09611b13e 100644 --- a/tests/test_dipoleFitter.py +++ b/tests/test_dipoleFitter.py @@ -113,12 +113,12 @@ 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.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,50 +156,88 @@ 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, 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) + + # 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) + # 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'], + .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'], + .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. # 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 - def testDipoleTask(self): """Test the dipole fitting singleFramePlugin. @@ -269,7 +307,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