LSST Applications  21.0.0+04719a4bac,21.0.0-1-ga51b5d4+f5e6047307,21.0.0-11-g2b59f77+a9c1acf22d,21.0.0-11-ga42c5b2+86977b0b17,21.0.0-12-gf4ce030+76814010d2,21.0.0-13-g1721dae+760e7a6536,21.0.0-13-g3a573fe+768d78a30a,21.0.0-15-g5a7caf0+f21cbc5713,21.0.0-16-g0fb55c1+b60e2d390c,21.0.0-19-g4cded4ca+71a93a33c0,21.0.0-2-g103fe59+bb20972958,21.0.0-2-g45278ab+04719a4bac,21.0.0-2-g5242d73+3ad5d60fb1,21.0.0-2-g7f82c8f+8babb168e8,21.0.0-2-g8f08a60+06509c8b61,21.0.0-2-g8faa9b5+616205b9df,21.0.0-2-ga326454+8babb168e8,21.0.0-2-gde069b7+5e4aea9c2f,21.0.0-2-gecfae73+1d3a86e577,21.0.0-2-gfc62afb+3ad5d60fb1,21.0.0-25-g1d57be3cd+e73869a214,21.0.0-3-g357aad2+ed88757d29,21.0.0-3-g4a4ce7f+3ad5d60fb1,21.0.0-3-g4be5c26+3ad5d60fb1,21.0.0-3-g65f322c+e0b24896a3,21.0.0-3-g7d9da8d+616205b9df,21.0.0-3-ge02ed75+a9c1acf22d,21.0.0-4-g591bb35+a9c1acf22d,21.0.0-4-g65b4814+b60e2d390c,21.0.0-4-gccdca77+0de219a2bc,21.0.0-4-ge8a399c+6c55c39e83,21.0.0-5-gd00fb1e+05fce91b99,21.0.0-6-gc675373+3ad5d60fb1,21.0.0-64-g1122c245+4fb2b8f86e,21.0.0-7-g04766d7+cd19d05db2,21.0.0-7-gdf92d54+04719a4bac,21.0.0-8-g5674e7b+d1bd76f71f,master-gac4afde19b+a9c1acf22d,w.2021.13
LSST Data Management Base Package
kernelCandidateQa.py
Go to the documentation of this file.
1 #
2 # LSST Data Management System
3 # Copyright 2008-2016 LSST Corporation.
4 #
5 # This product includes software developed by the
6 # LSST Project (http://www.lsst.org/).
7 #
8 # This program is free software: you can redistribute it and/or modify
9 # it under the terms of the GNU General Public License as published by
10 # the Free Software Foundation, either version 3 of the License, or
11 # (at your option) any later version.
12 #
13 # This program is distributed in the hope that it will be useful,
14 # but WITHOUT ANY WARRANTY; without even the implied warranty of
15 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 # GNU General Public License for more details.
17 #
18 # You should have received a copy of the LSST License Statement and
19 # the GNU General Public License along with this program. If not,
20 # see <http://www.lsstcorp.org/LegalNotices/>.
21 #
22 
23 __all__ = ["KernelCandidateQa"]
24 
25 import numpy as np
26 import numpy.ma as ma
27 
28 import lsst.afw.image as afwImage
29 import lsst.afw.table as afwTable
30 import lsst.afw.math as afwMath
31 import lsst.geom as geom
32 from . import diffimLib
33 from .utils import calcCentroid, calcWidth
34 
35 
37  """Quality Assessment class for Kernel Candidates"""
38 
39  def __init__(self, nKernelSpatial):
40  """Class to undertake QA of KernelCandidates after modeling of
41  the Psf-matching kernel. Both directly--fitted diffim (LOCAL)
42  and spatially--interpolated kernel diffim (SPATIAL) metrics
43  are calculated, based on the distribution of residuals in the
44  KernelCandidates stamp.
45 
46  Parameters
47  ----------
48  nKernelSpatial : `int`
49  Number of terms in the spatial model; needed to initialize per-basis QA arrays
50  """
51  self.fieldsfields = []
52  self.fieldsfields.append(afwTable.Field["PointD"](
53  "RegisterRefPosition",
54  "Position of reference object for registration (radians)."))
55  # TODO check units of the following angles
56  self.fieldsfields.append(afwTable.Field["Angle"]("RegisterResidualBearing",
57  "Angle of residual wrt declination parallel in radians"))
58 
59  self.fieldsfields.append(afwTable.Field["Angle"]("RegisterResidualDistance",
60  "Offset of residual in radians"))
61  metricMap = self.makeMetricMapmakeMetricMap()
62 
63  for kType in ("LOCAL", "SPATIAL"):
64  for k in metricMap:
65  commentAndUnit = metricMap[k]['comment']
66  self.fieldsfields.append(afwTable.Field[metricMap[k]['type']](k%(kType), *commentAndUnit))
67 
68  self.fieldsfields.append(afwTable.Field["I"]("KCKernelStatus_LOCAL",
69  "Status of the KernelCandidate"))
70 
71  self.fieldsfields.append(afwTable.Field["ArrayD"]("KernelCoeffValues_LOCAL",
72  "Original basis coefficients",
73  nKernelSpatial))
74 
75  self.fieldsfields.append(afwTable.Field["F"]("BackgroundValue_LOCAL",
76  "Evaluation of background model at this point"))
77 
78  self.fieldsfields.append(afwTable.Field["F"]("KCDiffimMseKernel_SPATIAL",
79  "Mean squared error of spatial kernel estimate"))
80 
81  def makeMetricMap(self):
82  nameList = ['KCDiffimMean_%s', 'KCDiffimMedian_%s', 'KCDiffimIQR_%s', 'KCDiffimStDev_%s',
83  'KCDiffimKSD_%s', 'KCDiffimKSProb_%s', 'KCDiffimADA2_%s', 'KCDiffimADCrit_%s',
84  'KCDiffimADSig_%s', 'KCDiffimChiSq_%s', 'KCDiffimMseResids_%s', 'KCKernelCentX_%s',
85  'KCKernelCentY_%s', 'KCKernelStdX_%s', 'KCKernelStdY_%s', 'KernelCandidateId_%s']
86  typeList = ['F', 'F', 'F', 'F', 'F', 'F', 'F', 'ArrayD', 'ArrayD', 'F', 'F', 'F',
87  'F', 'F', 'F', 'I']
88  commentList = [
89  ("Mean of KernelCandidate diffim", "sigma"),
90  ("Median of KernelCandidate diffim", "sigma"),
91  ("Inner quartile range of KernelCandidate diffim", "sigma"),
92  ("Standard deviation of KernelCandidate diffim", "sigma"),
93  ("D from K-S test of diffim pixels relative to Normal", ),
94  ("Prob from K-S test of diffim pixels relative to Normal", "likelihood"),
95  ("Anderson-Darling test statistic of diffim pixels relative to Normal", ),
96  ("Critical values for the significance levels in KCDiffimADSig. If A2 is greater "
97  "than this number, hypothesis that the distributions are similar can be rejected.", 5),
98  ("Anderson-Darling significance levels for the Normal distribution", 5),
99  ("Reduced chi^2 of the residual.", "likelihood"),
100  ("Mean squared error in diffim : Variance + Bias**2",),
101  ("Centroid in X for this Kernel", "pixel"),
102  ("Centroid in Y for this Kernel", "pixel"),
103  ("Standard deviation in X for this Kernel", "pixel"),
104  ("Standard deviation in Y for this Kernel", "pixel"),
105  ("Id for this KernelCandidate",),
106  ]
107  metricMap = {}
108  for name, mtype, comment in zip(nameList, typeList, commentList):
109  metricMap[name] = {'type': mtype, 'comment': comment}
110 
111  return metricMap
112 
113  def addToSchema(self, inSourceCatalog):
114  """Add the to-be-generated QA keys to the Source schema"""
115  schema = inSourceCatalog.getSchema()
116  inKeys = []
117  psfDef = inSourceCatalog.schema.getAliasMap().get("slot_PsfFlux")
118  centroidDef = inSourceCatalog.getCentroidDefinition()
119  shapeDef = inSourceCatalog.getShapeDefinition()
120  for n in schema.getNames():
121  inKeys.append(schema[n].asKey())
122 
123  for field in self.fieldsfields:
124  schema.addField(field)
125  outSourceCatalog = afwTable.SourceCatalog(schema)
126  for source in inSourceCatalog:
127  rec = outSourceCatalog.addNew()
128  for k in inKeys:
129  if k.getTypeString() == 'Coord':
130  rec.setCoord(source.getCoord())
131  else:
132  setter = getattr(rec, "set"+k.getTypeString())
133  getter = getattr(source, "get"+k.getTypeString())
134  setter(k, getter(k))
135  outSourceCatalog.definePsfFlux(psfDef)
136  outSourceCatalog.defineCentroid(centroidDef)
137  outSourceCatalog.defineShape(shapeDef)
138  return outSourceCatalog
139 
140  @staticmethod
141  def _calculateStats(di, dof=0.):
142  """Calculate the core QA statistics on a difference image"""
143  mask = di.getMask()
144  maskArr = di.getMask().getArray()
145 
146  # Create a mask using BAD, SAT, NO_DATA, EDGE bits. Keep detections
147  maskArr &= mask.getPlaneBitMask(["BAD", "SAT", "NO_DATA", "EDGE"])
148 
149  # Mask out values based on maskArr
150  diArr = ma.array(di.getImage().getArray(), mask=maskArr)
151  varArr = ma.array(di.getVariance().getArray(), mask=maskArr)
152 
153  # Normalize by sqrt variance, units are in sigma
154  diArr /= np.sqrt(varArr)
155  mean = diArr.mean()
156 
157  # This is the maximum-likelihood extimate of the variance stdev**2
158  stdev = diArr.std()
159  median = ma.extras.median(diArr)
160 
161  # Compute IQR of just un-masked data
162  data = ma.getdata(diArr[~diArr.mask])
163  iqr = np.percentile(data, 75.) - np.percentile(data, 25.)
164 
165  # Calculte chisquare of the residual
166  chisq = np.sum(np.power(data, 2.))
167 
168  # Mean squared error: variance + bias**2
169  # Bias = |data - model| = mean of diffim
170  # Variance = |(data - model)**2| = mean of diffim**2
171  bias = mean
172  variance = np.power(data, 2.).mean()
173  mseResids = bias**2 + variance
174 
175  # If scipy is not set up, return zero for the stats
176  try:
177  # In try block because of risk of divide by zero
178  rchisq = chisq/(len(data) - 1 - dof)
179  # K-S test on the diffim to a Normal distribution
180  import scipy.stats
181  D, prob = scipy.stats.kstest(data, 'norm')
182 
183  A2, crit, sig = scipy.stats.anderson(data, 'norm')
184  # Anderson Darling statistic cand be inf for really non-Gaussian distributions.
185  if np.isinf(A2) or np.isnan(A2):
186  A2 = 9999.
187  except ZeroDivisionError:
188  D = 0.
189  prob = 0.
190  A2 = 0.
191  crit = np.zeros(5)
192  sig = np.zeros(5)
193  rchisq = 0
194 
195  return {"mean": mean, "stdev": stdev, "median": median, "iqr": iqr,
196  "D": D, "prob": prob, "A2": A2, "crit": crit, "sig": sig,
197  "rchisq": rchisq, "mseResids": mseResids}
198 
199  @classmethod
200  def apply(cls, candidateList, spatialKernel, spatialBackground, dof=0):
201  """Evaluate the QA metrics for all KernelCandidates in the
202  candidateList; set the values of the metrics in their
203  associated Sources"""
204  for kernelCandidate in candidateList:
205  source = kernelCandidate.getSource()
206  schema = source.schema
207 
208  # Calculate ORIG stats (original basis fit)
209  if kernelCandidate.getStatus() != afwMath.SpatialCellCandidate.UNKNOWN:
210  kType = getattr(diffimLib.KernelCandidateF, "ORIG")
211  di = kernelCandidate.getDifferenceImage(kType)
212  kernelValues = kernelCandidate.getKernel(kType).getKernelParameters()
213  kernelValues = np.asarray(kernelValues)
214 
215  lkim = kernelCandidate.getKernelImage(kType)
216  centx, centy = calcCentroid(lkim.getArray())
217  stdx, stdy = calcWidth(lkim.getArray(), centx, centy)
218  # NOTE
219  # What is the difference between kernelValues and solution?
220 
221  localResults = cls._calculateStats_calculateStats(di, dof=dof)
222 
223  metrics = {"KCDiffimMean_LOCAL": localResults["mean"],
224  "KCDiffimMedian_LOCAL": localResults["median"],
225  "KCDiffimIQR_LOCAL": localResults["iqr"],
226  "KCDiffimStDev_LOCAL": localResults["stdev"],
227  "KCDiffimKSD_LOCAL": localResults["D"],
228  "KCDiffimKSProb_LOCAL": localResults["prob"],
229  "KCDiffimADA2_LOCAL": localResults["A2"],
230  "KCDiffimADCrit_LOCAL": localResults["crit"],
231  "KCDiffimADSig_LOCAL": localResults["sig"],
232  "KCDiffimChiSq_LOCAL": localResults["rchisq"],
233  "KCDiffimMseResids_LOCAL": localResults["mseResids"],
234  "KCKernelCentX_LOCAL": centx,
235  "KCKernelCentY_LOCAL": centy,
236  "KCKernelStdX_LOCAL": stdx,
237  "KCKernelStdY_LOCAL": stdy,
238  "KernelCandidateId_LOCAL": kernelCandidate.getId(),
239  "KernelCoeffValues_LOCAL": kernelValues}
240  for k in metrics:
241  key = schema[k].asKey()
242  setter = getattr(source, "set"+key.getTypeString())
243  setter(key, metrics[k])
244  else:
245  try:
246  kType = getattr(diffimLib.KernelCandidateF, "ORIG")
247  lkim = kernelCandidate.getKernelImage(kType)
248  except Exception:
249  lkim = None
250 
251  # Calculate spatial model evaluated at each position, for
252  # all candidates
253  skim = afwImage.ImageD(spatialKernel.getDimensions())
254  spatialKernel.computeImage(skim, False, kernelCandidate.getXCenter(),
255  kernelCandidate.getYCenter())
256  centx, centy = calcCentroid(skim.getArray())
257  stdx, stdy = calcWidth(skim.getArray(), centx, centy)
258 
259  sk = afwMath.FixedKernel(skim)
260  sbg = spatialBackground(kernelCandidate.getXCenter(), kernelCandidate.getYCenter())
261  di = kernelCandidate.getDifferenceImage(sk, sbg)
262  spatialResults = cls._calculateStats_calculateStats(di, dof=dof)
263 
264  # Kernel mse
265  if lkim is not None:
266  skim -= lkim
267  bias = np.mean(skim.getArray())
268  variance = np.mean(np.power(skim.getArray(), 2.))
269  mseKernel = bias**2 + variance
270  else:
271  mseKernel = -99.999
272 
273  metrics = {"KCDiffimMean_SPATIAL": spatialResults["mean"],
274  "KCDiffimMedian_SPATIAL": spatialResults["median"],
275  "KCDiffimIQR_SPATIAL": spatialResults["iqr"],
276  "KCDiffimStDev_SPATIAL": spatialResults["stdev"],
277  "KCDiffimKSD_SPATIAL": spatialResults["D"],
278  "KCDiffimKSProb_SPATIAL": spatialResults["prob"],
279  "KCDiffimADA2_SPATIAL": spatialResults["A2"],
280  "KCDiffimADCrit_SPATIAL": spatialResults["crit"],
281  "KCDiffimADSig_SPATIAL": spatialResults["sig"],
282  "KCDiffimChiSq_SPATIAL": spatialResults["rchisq"],
283  "KCDiffimMseResids_SPATIAL": spatialResults["mseResids"],
284  "KCDiffimMseKernel_SPATIAL": mseKernel,
285  "KCKernelCentX_SPATIAL": centx,
286  "KCKernelCentY_SPATIAL": centy,
287  "KCKernelStdX_SPATIAL": stdx,
288  "KCKernelStdY_SPATIAL": stdy,
289  "KernelCandidateId_SPATIAL": kernelCandidate.getId()}
290  for k in metrics:
291  key = schema[k].asKey()
292  setter = getattr(source, "set"+key.getTypeString())
293  setter(key, metrics[k])
294 
295  @staticmethod
296  def aggregate(sourceCatalog, metadata, wcsresids, diaSources=None):
297  """Generate aggregate metrics (e.g. total numbers of false
298  positives) from all the Sources in the sourceCatalog"""
299  for source in sourceCatalog:
300  sourceId = source.getId()
301  if sourceId in wcsresids:
302  # Note that the residuals are not delta RA, delta Dec
303  # From the source code "bearing (angle wrt a declination parallel) and distance
304  coord, resids = wcsresids[sourceId]
305  key = source.schema["RegisterResidualBearing"].asKey()
306  setter = getattr(source, "set"+key.getTypeString())
307  setter(key, resids[0])
308  key = source.schema["RegisterResidualDistance"].asKey()
309  setter = getattr(source, "set"+key.getTypeString())
310  setter(key, resids[1])
311  key = source.schema["RegisterRefPosition"].asKey()
312  setter = getattr(source, "set"+key.getTypeString())
313  setter(key, geom.Point2D(coord.getRa().asRadians(),
314  coord.getDec().asRadians()))
315  if diaSources:
316  metadata.add("NFalsePositivesTotal", len(diaSources))
317  nRefMatch = 0
318  nSrcMatch = 0
319  nunmatched = 0
320  for source in diaSources:
321  refId = source.get("refMatchId")
322  srcId = source.get("srcMatchId")
323  if refId > 0:
324  nRefMatch += 1
325  if srcId > 0:
326  nSrcMatch += 1
327  if refId == 0 and srcId == 0:
328  nunmatched += 1
329  metadata.add("NFalsePositivesRefAssociated", nRefMatch)
330  metadata.add("NFalsePositivesSrcAssociated", nSrcMatch)
331  metadata.add("NFalsePositivesUnassociated", nunmatched)
332  for kType in ("LOCAL", "SPATIAL"):
333  for sName in ("KCDiffimMean", "KCDiffimMedian", "KCDiffimIQR", "KCDiffimStDev",
334  "KCDiffimKSProb", "KCDiffimADSig", "KCDiffimChiSq",
335  "KCDiffimMseResids", "KCDiffimMseKernel"):
336  if sName == "KCDiffimMseKernel" and kType == "LOCAL":
337  continue
338  kName = "%s_%s" % (sName, kType)
339  vals = np.array([s.get(kName) for s in sourceCatalog])
340  idx = np.isfinite(vals)
341  metadata.add("%s_MEAN" % (kName), np.mean(vals[idx]))
342  metadata.add("%s_MEDIAN" % (kName), np.median(vals[idx]))
343  metadata.add("%s_STDEV" % (kName), np.std(vals[idx]))
A kernel created from an Image.
Definition: Kernel.h:472
def apply(cls, candidateList, spatialKernel, spatialBackground, dof=0)
def aggregate(sourceCatalog, metadata, wcsresids, diaSources=None)
std::shared_ptr< FrameSet > append(FrameSet const &first, FrameSet const &second)
Construct a FrameSet that performs two transformations in series.
Definition: functional.cc:33
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects.
def calcWidth(arr, centx, centy)
Definition: utils.py:829
def calcCentroid(arr)
Definition: utils.py:814
A description of a field in a table.
Definition: Field.h:24