LSSTApplications  10.0-2-g4f67435,11.0.rc2+1,11.0.rc2+12,11.0.rc2+3,11.0.rc2+4,11.0.rc2+5,11.0.rc2+6,11.0.rc2+7,11.0.rc2+8
LSSTDataManagementBasePackage
kernelCandidateQa.py
Go to the documentation of this file.
1 #
2 # LSST Data Management System
3 # Copyright 2008, 2009, 2010 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 """Quality Assessment class for Kernel Candidates"""
24 import numpy as np
25 import numpy.ma as ma
26 import lsst.afw.geom as afwGeom
27 import lsst.afw.image as afwImage
28 import lsst.afw.table as afwTable
29 import lsst.afw.math as afwMath
30 from . import diffimLib
31 from .utils import calcCentroid, calcWidth
32 
33 class KernelCandidateQa(object):
34 
35  def __init__(self, nKernelSpatial):
36  """Class to undertake QA of KernelCandidates after modeling of
37  the Psf-matching kernel. Both directly--fitted diffim (LOCAL)
38  and spatially--interpolated kernel diffim (SPATIAL) metrics
39  are calculated, based on the distribution of residuals in the
40  KernelCandidates stamp.
41 
42  @param nKernelSpatial : Number of terms in the spatial model; needed to initialize per-basis QA arrays
43  """
44  self.fields = []
45  self.fields.append(afwTable.Field["PointD"]("RegisterRefPosition",
46  "Position of reference object for registration (radians)."))
47  #TODO check units of the following angles
48  self.fields.append(afwTable.Field["Angle"]("RegisterResidualBearing",
49  "Angle of residual wrt declination parallel in radians"))
50 
51  self.fields.append(afwTable.Field["Angle"]("RegisterResidualDistance",
52  "Offset of residual in radians"))
53  metricMap = self.makeMetricMap()
54 
55  for kType in ("LOCAL", "SPATIAL"):
56  for k in metricMap:
57  commentAndUnit = metricMap[k]['comment']
58  self.fields.append(afwTable.Field[metricMap[k]['type']](k%(kType), *commentAndUnit))
59 
60  self.fields.append(afwTable.Field["I"]("KCKernelStatus_LOCAL",
61  "Status of the KernelCandidate"))
62 
63  self.fields.append(afwTable.Field["ArrayD"]("KernelCoeffValues_LOCAL",
64  "Original basis coefficients",
65  nKernelSpatial))
66 
67  self.fields.append(afwTable.Field["F"]("BackgroundValue_LOCAL",
68  "Evaluation of background model at this point"))
69 
70  self.fields.append(afwTable.Field["F"]("KCDiffimMseKernel_SPATIAL",
71  "Mean squared error of spatial kernel estimate"))
72 
73  def makeMetricMap(self):
74  nameList = ['KCDiffimMean_%s', 'KCDiffimMedian_%s', 'KCDiffimIQR_%s', 'KCDiffimStDev_%s',
75  'KCDiffimKSD_%s', 'KCDiffimKSProb_%s', 'KCDiffimADA2_%s', 'KCDiffimADCrit_%s',
76  'KCDiffimADSig_%s', 'KCDiffimChiSq_%s', 'KCDiffimMseResids_%s', 'KCKernelCentX_%s',
77  'KCKernelCentY_%s', 'KCKernelStdX_%s', 'KCKernelStdY_%s', 'KernelCandidateId_%s']
78  typeList = ['F', 'F', 'F', 'F', 'F', 'F', 'F', 'ArrayD', 'ArrayD', 'F', 'F', 'F',
79  'F', 'F', 'F', 'I']
80  commentList = [("Mean of KernelCandidate diffim", "sigma"),
81  ("Median of KernelCandidate diffim", "sigma"),
82  ("Inner quartile range of KernelCandidate diffim", "sigma"),
83  ("Standard deviation of KernelCandidate diffim","sigma"),
84  ("D from K-S test of diffim pixels relative to Normal", ),
85  ("Prob from K-S test of diffim pixels relative to Normal", "likelihood"),
86  ("Anderson-Darling test statistic of diffim pixels relative to Normal", ),
87  ("Critical values for the significance levels in KCDiffimADSig. If A2 is greater "+\
88  "than this number, hypothesis that the distributions are similar can be rejected.", 5),
89  ("Anderson-Darling significance levels for the Normal distribution", 5),
90  ("Reduced chi^2 of the residual.", "likelihood"),
91  ("Mean squared error in diffim : Variance + Bias**2",),
92  ("Centroid in X for this Kernel", "pixels"),
93  ("Centroid in Y for this Kernel", "pixels"),
94  ("Standard deviation in X for this Kernel", "pixels"),
95  ("Standard deviation in Y for this Kernel", "pixels"),
96  ("Id for this KernelCandidate",)]
97  metricMap = {}
98  for name, mtype, comment in zip(nameList, typeList, commentList):
99  metricMap[name] = {'type':mtype, 'comment':comment}
100 
101  return metricMap
102 
103 
104  def addToSchema(self, inSourceCatalog):
105  """Add the to-be-generated QA keys to the Source schema"""
106  schema = inSourceCatalog.getSchema()
107  inKeys = []
108  psfDef = inSourceCatalog.getPsfFluxDefinition()
109  centroidDef = inSourceCatalog.getCentroidDefinition()
110  shapeDef = inSourceCatalog.getShapeDefinition()
111  for n in schema.getNames():
112  inKeys.append(schema[n].asKey())
113 
114  for field in self.fields:
115  schema.addField(field)
116  outSourceCatalog = afwTable.SourceCatalog(schema)
117  for source in inSourceCatalog:
118  rec = outSourceCatalog.addNew()
119  for k in inKeys:
120  if k.getTypeString() == 'Coord':
121  rec.setCoord(source.getCoord())
122  else:
123  setter = getattr(rec, "set"+k.getTypeString())
124  getter = getattr(source, "get"+k.getTypeString())
125  setter(k, getter(k))
126  outSourceCatalog.definePsfFlux(psfDef)
127  outSourceCatalog.defineCentroid(centroidDef)
128  outSourceCatalog.defineShape(shapeDef)
129  return outSourceCatalog
130 
131  def _calculateStats(self, di, dof=0.):
132  """Calculate the core QA statistics on a difference image"""
133  mask = di.getMask()
134  maskArr = di.getMask().getArray()
135 
136  # Create a mask using BAD, SAT, NO_DATA, EDGE bits. Keep detections
137  maskArr &= mask.getPlaneBitMask(["BAD", "SAT", "NO_DATA", "EDGE"])
138 
139  # Mask out values based on maskArr
140  diArr = ma.array(di.getImage().getArray(), mask=maskArr)
141  varArr = ma.array(di.getVariance().getArray(), mask=maskArr)
142 
143  # Normalize by sqrt variance, units are in sigma
144  diArr /= np.sqrt(varArr)
145  mean = diArr.mean()
146 
147  # This is the maximum-likelihood extimate of the variance stdev**2
148  stdev = diArr.std()
149  median = ma.extras.median(diArr)
150 
151  # Compute IQR of just un-masked data
152  data = ma.getdata(diArr[~diArr.mask])
153  iqr = np.percentile(data, 75.) - np.percentile(data, 25.)
154 
155  #Calculte chisquare of the residual
156  chisq=np.sum(np.power(data, 2.))
157 
158  # Mean squared error: variance + bias**2
159  # Bias = |data - model| = mean of diffim
160  # Variance = |(data - model)**2| = mean of diffim**2
161  bias = mean
162  variance = np.power(data, 2.).mean()
163  mseResids = bias**2 + variance
164 
165  # If scipy is not set up, return zero for the stats
166  try:
167  #In try block because of risk of divide by zero
168  rchisq=chisq/(len(data)-1-dof)
169  # K-S test on the diffim to a Normal distribution
170  import scipy.stats
171  D, prob = scipy.stats.kstest(data, 'norm')
172 
173  A2, crit, sig = scipy.stats.anderson(data, 'norm')
174  # Anderson Darling statistic cand be inf for really non-Gaussian distributions.
175  if np.isinf(A2) or np.isnan(A2):
176  A2 = 9999.
177  except ZeroDivisionError:
178  D = 0.
179  prob = 0.
180  A2 = 0.
181  crit = np.zeros(5)
182  sig = np.zeros(5)
183  rchisq = 0
184 
185  return {"mean": mean, "stdev": stdev, "median": median, "iqr": iqr,
186  "D": D, "prob": prob, "A2": A2, "crit": crit, "sig": sig,
187  "rchisq": rchisq, "mseResids": mseResids}
188 
189 
190  def apply(self, candidateList, spatialKernel, spatialBackground, dof=0):
191  """Evaluate the QA metrics for all KernelCandidates in the
192  candidateList; set the values of the metrics in their
193  associated Sources"""
194  for kernelCandidate in candidateList:
195  source = kernelCandidate.getSource()
196  schema = source.schema
197 
198  # Calculate ORIG stats (original basis fit)
199  if kernelCandidate.getStatus() != afwMath.SpatialCellCandidate.UNKNOWN:
200  kType = getattr(diffimLib.KernelCandidateF, "ORIG")
201  di = kernelCandidate.getDifferenceImage(kType)
202  kernelValues = kernelCandidate.getKernel(kType).getKernelParameters()
203  kernelValues = np.asarray(kernelValues)
204 
205  lkim = kernelCandidate.getKernelImage(kType)
206  centx, centy = calcCentroid(lkim.getArray())
207  stdx, stdy = calcWidth(lkim.getArray(), centx, centy)
208  # NOTE
209  # What is the difference between kernelValues and solution?
210 
211  localResults = self._calculateStats(di, dof=dof)
212 
213  metrics = {"KCDiffimMean_LOCAL":localResults["mean"],
214  "KCDiffimMedian_LOCAL":localResults["median"],
215  "KCDiffimIQR_LOCAL":localResults["iqr"],
216  "KCDiffimStDev_LOCAL":localResults["stdev"],
217  "KCDiffimKSD_LOCAL":localResults["D"],
218  "KCDiffimKSProb_LOCAL":localResults["prob"],
219  "KCDiffimADA2_LOCAL":localResults["A2"],
220  "KCDiffimADCrit_LOCAL":localResults["crit"],
221  "KCDiffimADSig_LOCAL":localResults["sig"],
222  "KCDiffimChiSq_LOCAL":localResults["rchisq"],
223  "KCDiffimMseResids_LOCAL":localResults["mseResids"],
224  "KCKernelCentX_LOCAL":centx,
225  "KCKernelCentY_LOCAL":centy,
226  "KCKernelStdX_LOCAL":stdx,
227  "KCKernelStdY_LOCAL":stdy,
228  "KernelCandidateId_LOCAL":kernelCandidate.getId(),
229  "KernelCoeffValues_LOCAL":kernelValues}
230  for k in metrics.keys():
231  key = schema[k].asKey()
232  setter = getattr(source, "set"+key.getTypeString())
233  setter(key, metrics[k])
234  else:
235  try:
236  kType = getattr(diffimLib.KernelCandidateF, "ORIG")
237  lkim = kernelCandidate.getKernelImage(kType)
238  except Exception:
239  lkim = None
240 
241  # Calculate spatial model evaluated at each position, for
242  # all candidates
243  skim = afwImage.ImageD(spatialKernel.getDimensions())
244  spatialKernel.computeImage(skim, False, kernelCandidate.getXCenter(),
245  kernelCandidate.getYCenter())
246  centx, centy = calcCentroid(skim.getArray())
247  stdx, stdy = calcWidth(skim.getArray(), centx, centy)
248 
249  sk = afwMath.FixedKernel(skim)
250  sbg = spatialBackground(kernelCandidate.getXCenter(), kernelCandidate.getYCenter())
251  di = kernelCandidate.getDifferenceImage(sk, sbg)
252  spatialResults = self._calculateStats(di, dof=dof)
253 
254  # Kernel mse
255  if lkim is not None:
256  skim -= lkim
257  bias = np.mean(skim.getArray())
258  variance = np.mean(np.power(skim.getArray(), 2.))
259  mseKernel = bias**2 + variance
260  else:
261  mseKernel = -99.999
262 
263  metrics = {"KCDiffimMean_SPATIAL":spatialResults["mean"],
264  "KCDiffimMedian_SPATIAL":spatialResults["median"],
265  "KCDiffimIQR_SPATIAL":spatialResults["iqr"],
266  "KCDiffimStDev_SPATIAL":spatialResults["stdev"],
267  "KCDiffimKSD_SPATIAL":spatialResults["D"],
268  "KCDiffimKSProb_SPATIAL":spatialResults["prob"],
269  "KCDiffimADA2_SPATIAL":spatialResults["A2"],
270  "KCDiffimADCrit_SPATIAL":spatialResults["crit"],
271  "KCDiffimADSig_SPATIAL":spatialResults["sig"],
272  "KCDiffimChiSq_SPATIAL":spatialResults["rchisq"],
273  "KCDiffimMseResids_SPATIAL":spatialResults["mseResids"],
274  "KCDiffimMseKernel_SPATIAL":mseKernel,
275  "KCKernelCentX_SPATIAL":centx,
276  "KCKernelCentY_SPATIAL":centy,
277  "KCKernelStdX_SPATIAL":stdx,
278  "KCKernelStdY_SPATIAL":stdy,
279  "KernelCandidateId_SPATIAL":kernelCandidate.getId()}
280  for k in metrics.keys():
281  key = schema[k].asKey()
282  setter = getattr(source, "set"+key.getTypeString())
283  setter(key, metrics[k])
284 
285  def aggregate(self, sourceCatalog, metadata, wcsresids, diaSources=None):
286  """Generate aggregate metrics (e.g. total numbers of false
287  positives) from all the Sources in the sourceCatalog"""
288  for source in sourceCatalog:
289  sourceId = source.getId()
290  if sourceId in wcsresids.keys():
291  #Note that the residuals are not delta RA, delta Dec
292  #From the source code "bearing (angle wrt a declination parallel) and distance
293  coord, resids = wcsresids[sourceId]
294  key = source.schema["RegisterResidualBearing"].asKey()
295  setter = getattr(source, "set"+key.getTypeString())
296  setter(key, resids[0])
297  key = source.schema["RegisterResidualDistance"].asKey()
298  setter = getattr(source, "set"+key.getTypeString())
299  setter(key, resids[1])
300  key = source.schema["RegisterRefPosition"].asKey()
301  setter = getattr(source, "set"+key.getTypeString())
302  setter(key, afwGeom.Point2D(coord.getRa().asRadians(),
303  coord.getDec().asRadians()))
304  if diaSources:
305  metadata.add("NFalsePositivesTotal", len(diaSources))
306  nRefMatch = 0
307  nSrcMatch = 0
308  nunmatched = 0
309  for source in diaSources:
310  refId = source.get("refMatchId")
311  srcId = source.get("srcMatchId")
312  if refId > 0:
313  nRefMatch += 1
314  if srcId > 0:
315  nSrcMatch += 1
316  if refId == 0 and srcId == 0:
317  nunmatched += 1
318  metadata.add("NFalsePositivesRefAssociated", nRefMatch)
319  metadata.add("NFalsePositivesSrcAssociated", nSrcMatch)
320  metadata.add("NFalsePositivesUnassociated", nunmatched)
321  for kType in ("LOCAL", "SPATIAL"):
322  for sName in ("KCDiffimMean", "KCDiffimMedian", "KCDiffimIQR", "KCDiffimStDev",
323  "KCDiffimKSProb", "KCDiffimADSig", "KCDiffimChiSq",
324  "KCDiffimMseResids", "KCDiffimMseKernel"):
325  if sName == "KCDiffimMseKernel" and kType == "LOCAL":
326  continue
327  kName = "%s_%s" % (sName, kType)
328  vals = np.array([s.get(kName) for s in sourceCatalog])
329  idx = np.isfinite(vals)
330  metadata.add("%s_MEAN" % (kName), np.mean(vals[idx]))
331  metadata.add("%s_MEDIAN" % (kName), np.median(vals[idx]))
332  metadata.add("%s_STDEV" % (kName), np.std(vals[idx]))
333 
Custom catalog class for record/table subclasses that are guaranteed to have an ID, and should generally be sorted by that ID.
Definition: fwd.h:55
A description of a field in a table.
Definition: Field.h:22
A kernel created from an Image.
Definition: Kernel.h:551