LSST Applications  21.0.0-172-gfb10e10a+18fedfabac,22.0.0+297cba6710,22.0.0+80564b0ff1,22.0.0+8d77f4f51a,22.0.0+a28f4c53b1,22.0.0+dcf3732eb2,22.0.1-1-g7d6de66+2a20fdde0d,22.0.1-1-g8e32f31+297cba6710,22.0.1-1-geca5380+7fa3b7d9b6,22.0.1-12-g44dc1dc+2a20fdde0d,22.0.1-15-g6a90155+515f58c32b,22.0.1-16-g9282f48+790f5f2caa,22.0.1-2-g92698f7+dcf3732eb2,22.0.1-2-ga9b0f51+7fa3b7d9b6,22.0.1-2-gd1925c9+bf4f0e694f,22.0.1-24-g1ad7a390+a9625a72a8,22.0.1-25-g5bf6245+3ad8ecd50b,22.0.1-25-gb120d7b+8b5510f75f,22.0.1-27-g97737f7+2a20fdde0d,22.0.1-32-gf62ce7b1+aa4237961e,22.0.1-4-g0b3f228+2a20fdde0d,22.0.1-4-g243d05b+871c1b8305,22.0.1-4-g3a563be+32dcf1063f,22.0.1-4-g44f2e3d+9e4ab0f4fa,22.0.1-42-gca6935d93+ba5e5ca3eb,22.0.1-5-g15c806e+85460ae5f3,22.0.1-5-g58711c4+611d128589,22.0.1-5-g75bb458+99c117b92f,22.0.1-6-g1c63a23+7fa3b7d9b6,22.0.1-6-g50866e6+84ff5a128b,22.0.1-6-g8d3140d+720564cf76,22.0.1-6-gd805d02+cc5644f571,22.0.1-8-ge5750ce+85460ae5f3,master-g6e05de7fdc+babf819c66,master-g99da0e417a+8d77f4f51a,w.2021.48
LSST Data Management Base Package
diaCalculationPlugins.py
Go to the documentation of this file.
1 # This file is part of ap_association.
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 """Plugins for use in DiaSource summary statistics.
23 
24 Output columns must be
25 as defined in the schema of the Apdb both in name and units.
26 """
27 
28 import functools
29 import warnings
30 
31 from astropy.stats import median_absolute_deviation
32 import numpy as np
33 import pandas as pd
34 from scipy.optimize import lsq_linear
35 
36 import lsst.geom as geom
37 import lsst.pex.config as pexConfig
38 import lsst.sphgeom as sphgeom
39 
40 from .diaCalculation import (
41  DiaObjectCalculationPluginConfig,
42  DiaObjectCalculationPlugin)
43 from .pluginRegistry import register
44 
45 __all__ = ("MeanDiaPositionConfig", "MeanDiaPosition",
46  "HTMIndexDiaPosition", "HTMIndexDiaPositionConfig",
47  "NumDiaSourcesDiaPlugin", "NumDiaSourcesDiaPluginConfig",
48  "SimpleSourceFlagDiaPlugin", "SimpleSourceFlagDiaPluginConfig",
49  "WeightedMeanDiaPsFluxConfig", "WeightedMeanDiaPsFlux",
50  "PercentileDiaPsFlux", "PercentileDiaPsFluxConfig",
51  "SigmaDiaPsFlux", "SigmaDiaPsFluxConfig",
52  "Chi2DiaPsFlux", "Chi2DiaPsFluxConfig",
53  "MadDiaPsFlux", "MadDiaPsFluxConfig",
54  "SkewDiaPsFlux", "SkewDiaPsFluxConfig",
55  "MinMaxDiaPsFlux", "MinMaxDiaPsFluxConfig",
56  "MaxSlopeDiaPsFlux", "MaxSlopeDiaPsFluxConfig",
57  "ErrMeanDiaPsFlux", "ErrMeanDiaPsFluxConfig",
58  "LinearFitDiaPsFlux", "LinearFitDiaPsFluxConfig",
59  "StetsonJDiaPsFlux", "StetsonJDiaPsFluxConfig",
60  "WeightedMeanDiaTotFlux", "WeightedMeanDiaTotFluxConfig",
61  "SigmaDiaTotFlux", "SigmaDiaTotFluxConfig")
62 
63 
64 def catchWarnings(_func=None, *, warns=[]):
65  """Decorator for generically catching numpy warnings.
66  """
67  def decoratorCatchWarnings(func):
68  @functools.wraps(func)
69  def wrapperCatchWarnings(*args, **kwargs):
70  with warnings.catch_warnings():
71  for val in warns:
72  warnings.filterwarnings("ignore", val)
73  return func(*args, **kwargs)
74  return wrapperCatchWarnings
75 
76  if _func is None:
77  return decoratorCatchWarnings
78  else:
79  return decoratorCatchWarnings(_func)
80 
81 
83  pass
84 
85 
86 @register("ap_meanPosition")
88  """Compute the mean position of a DiaObject given a set of DiaSources.
89  """
90 
91  ConfigClass = MeanDiaPositionConfig
92 
93  plugType = 'multi'
94 
95  outputCols = ["ra", "decl", "radecTai"]
96  needsFilter = False
97 
98  @classmethod
100  return cls.DEFAULT_CATALOGCALCULATIONDEFAULT_CATALOGCALCULATION
101 
102  def calculate(self, diaObjects, diaSources, **kwargs):
103  """Compute the mean ra/dec position of the diaObject given the
104  diaSource locations.
105 
106  Parameters
107  ----------
108  diaObjects : `pandas.DataFrame`
109  Summary objects to store values in.
110  diaSources : `pandas.DataFrame` or `pandas.DataFrameGroupBy`
111  Catalog of DiaSources summarized by this DiaObject.
112  **kwargs
113  Any additional keyword arguments that may be passed to the plugin.
114  """
115  for outCol in self.outputColsoutputColsoutputCols:
116  if outCol not in diaObjects.columns:
117  diaObjects[outCol] = np.nan
118 
119  def _computeMeanPos(df):
120  aveCoord = geom.averageSpherePoint(
121  list(geom.SpherePoint(src["ra"], src["decl"], geom.degrees)
122  for idx, src in df.iterrows()))
123  ra = aveCoord.getRa().asDegrees()
124  decl = aveCoord.getDec().asDegrees()
125  if np.isnan(ra) or np.isnan(decl):
126  radecTai = np.nan
127  else:
128  radecTai = df["midPointTai"].max()
129 
130  return pd.Series({"ra": aveCoord.getRa().asDegrees(),
131  "decl": aveCoord.getDec().asDegrees(),
132  "radecTai": radecTai})
133 
134  ans = diaSources.apply(_computeMeanPos)
135  diaObjects.loc[:, ["ra", "decl", "radecTai"]] = ans
136 
137 
139 
140  htmLevel = pexConfig.Field(
141  dtype=int,
142  doc="Level of the HTM pixelization.",
143  default=20,
144  )
145 
146 
147 @register("ap_HTMIndex")
149  """Compute the mean position of a DiaObject given a set of DiaSources.
150 
151  Notes
152  -----
153  This plugin was implemented to satisfy requirements of old APDB interface
154  which required ``pixelId`` column in DiaObject with HTM20 index. APDB
155  interface had migrated to not need that information, but we keep this
156  plugin in case it may be useful for something else.
157  """
158  ConfigClass = HTMIndexDiaPositionConfig
159 
160  plugType = 'single'
161 
162  inputCols = ["ra", "decl"]
163  outputCols = ["pixelId"]
164  needsFilter = False
165 
166  def __init__(self, config, name, metadata):
167  DiaObjectCalculationPlugin.__init__(self, config, name, metadata)
168  self.pixelatorpixelator = sphgeom.HtmPixelization(self.configconfig.htmLevel)
169 
170  @classmethod
172  return cls.FLUX_MOMENTS_CALCULATEDFLUX_MOMENTS_CALCULATED
173 
174  def calculate(self, diaObjects, diaObjectId, **kwargs):
175  """Compute the mean position of a DiaObject given a set of DiaSources
176 
177  Parameters
178  ----------
179  diaObjects : `pandas.dataFrame`
180  Summary objects to store values in and read ra/decl from.
181  diaObjectId : `int`
182  Id of the diaObject to update.
183  **kwargs
184  Any additional keyword arguments that may be passed to the plugin.
185  """
186  sphPoint = geom.SpherePoint(
187  diaObjects.at[diaObjectId, "ra"] * geom.degrees,
188  diaObjects.at[diaObjectId, "decl"] * geom.degrees)
189  diaObjects.at[diaObjectId, "pixelId"] = self.pixelatorpixelator.index(
190  sphPoint.getVector())
191 
192 
194  pass
195 
196 
197 @register("ap_nDiaSources")
199  """Compute the total number of DiaSources associated with this DiaObject.
200  """
201 
202  ConfigClass = NumDiaSourcesDiaPluginConfig
203  outputCols = ["nDiaSources"]
204  plugType = "multi"
205  needsFilter = False
206 
207  @classmethod
209  return cls.DEFAULT_CATALOGCALCULATIONDEFAULT_CATALOGCALCULATION
210 
211  def calculate(self, diaObjects, diaSources, **kwargs):
212  """Compute the total number of DiaSources associated with this DiaObject.
213 
214  Parameters
215  ----------
216  diaObject : `dict`
217  Summary object to store values in and read ra/decl from.
218  **kwargs
219  Any additional keyword arguments that may be passed to the plugin.
220  """
221  diaObjects.loc[:, "nDiaSources"] = diaSources.diaObjectId.count()
222 
223 
225  pass
226 
227 
228 @register("ap_diaObjectFlag")
230  """Find if any DiaSource is flagged.
231 
232  Set the DiaObject flag if any DiaSource is flagged.
233  """
234 
235  ConfigClass = NumDiaSourcesDiaPluginConfig
236  outputCols = ["flags"]
237  plugType = "multi"
238  needsFilter = False
239 
240  @classmethod
242  return cls.DEFAULT_CATALOGCALCULATIONDEFAULT_CATALOGCALCULATION
243 
244  def calculate(self, diaObjects, diaSources, **kwargs):
245  """Find if any DiaSource is flagged.
246 
247  Set the DiaObject flag if any DiaSource is flagged.
248 
249  Parameters
250  ----------
251  diaObject : `dict`
252  Summary object to store values in and read ra/decl from.
253  **kwargs
254  Any additional keyword arguments that may be passed to the plugin.
255  """
256  diaObjects.loc[:, "flags"] = diaSources.flags.any().astype(np.uint64)
257 
258 
260  pass
261 
262 
263 @register("ap_meanFlux")
265  """Compute the weighted mean and mean error on the point source fluxes
266  of the DiaSource measured on the difference image.
267 
268  Additionally store number of usable data points.
269  """
270 
271  ConfigClass = WeightedMeanDiaPsFluxConfig
272  outputCols = ["PSFluxMean", "PSFluxMeanErr", "PSFluxNdata"]
273  plugType = "multi"
274  needsFilter = True
275 
276  @classmethod
278  return cls.DEFAULT_CATALOGCALCULATIONDEFAULT_CATALOGCALCULATION
279 
280  @catchWarnings(warns=["invalid value encountered", "divide by zero"])
281  def calculate(self,
282  diaObjects,
283  diaSources,
284  filterDiaSources,
285  filterName,
286  **kwargs):
287  """Compute the weighted mean and mean error of the point source flux.
288 
289  Parameters
290  ----------
291  diaObject : `dict`
292  Summary object to store values in.
293  diaSources : `pandas.DataFrame`
294  DataFrame representing all diaSources associated with this
295  diaObject.
296  filterDiaSources : `pandas.DataFrame`
297  DataFrame representing diaSources associated with this
298  diaObject that are observed in the band pass ``filterName``.
299  filterName : `str`
300  Simple, string name of the filter for the flux being calculated.
301  **kwargs
302  Any additional keyword arguments that may be passed to the plugin.
303  """
304  meanName = "{}PSFluxMean".format(filterName)
305  errName = "{}PSFluxMeanErr".format(filterName)
306  nDataName = "{}PSFluxNdata".format(filterName)
307  if meanName not in diaObjects.columns:
308  diaObjects[meanName] = np.nan
309  if errName not in diaObjects.columns:
310  diaObjects[errName] = np.nan
311  if nDataName not in diaObjects.columns:
312  diaObjects[nDataName] = 0
313 
314  def _weightedMean(df):
315  tmpDf = df[~np.logical_or(np.isnan(df["psFlux"]),
316  np.isnan(df["psFluxErr"]))]
317  tot_weight = np.nansum(1 / tmpDf["psFluxErr"] ** 2)
318  fluxMean = np.nansum(tmpDf["psFlux"]
319  / tmpDf["psFluxErr"] ** 2)
320  fluxMean /= tot_weight
321  if tot_weight > 0:
322  fluxMeanErr = np.sqrt(1 / tot_weight)
323  else:
324  fluxMeanErr = np.nan
325  nFluxData = len(tmpDf)
326 
327  return pd.Series({meanName: fluxMean,
328  errName: fluxMeanErr,
329  nDataName: nFluxData},
330  dtype="object")
331 
332  diaObjects.loc[:, [meanName, errName, nDataName]] = \
333  filterDiaSources.apply(_weightedMean)
334 
335 
337  percentiles = pexConfig.ListField(
338  dtype=int,
339  default=[5, 25, 50, 75, 95],
340  doc="Percentiles to calculate to compute values for. Should be "
341  "integer values."
342  )
343 
344 
345 @register("ap_percentileFlux")
347  """Compute percentiles of diaSource fluxes.
348  """
349 
350  ConfigClass = PercentileDiaPsFluxConfig
351  # Output columns are created upon instantiation of the class.
352  outputCols = []
353  plugType = "multi"
354  needsFilter = True
355 
356  def __init__(self, config, name, metadata, **kwargs):
357  DiaObjectCalculationPlugin.__init__(self,
358  config,
359  name,
360  metadata,
361  **kwargs)
362  self.outputColsoutputColsoutputColsoutputCols = ["PSFluxPercentile{:02d}".format(percent)
363  for percent in self.configconfig.percentiles]
364 
365  @classmethod
366  def getExecutionOrder(cls):
367  return cls.DEFAULT_CATALOGCALCULATIONDEFAULT_CATALOGCALCULATION
368 
369  @catchWarnings(warns=["All-NaN slice encountered"])
370  def calculate(self,
371  diaObjects,
372  diaSources,
373  filterDiaSources,
374  filterName,
375  **kwargs):
376  """Compute the percentile fluxes of the point source flux.
377 
378  Parameters
379  ----------
380  diaObject : `dict`
381  Summary object to store values in.
382  diaSources : `pandas.DataFrame`
383  DataFrame representing all diaSources associated with this
384  diaObject.
385  filterDiaSources : `pandas.DataFrame`
386  DataFrame representing diaSources associated with this
387  diaObject that are observed in the band pass ``filterName``.
388  filterName : `str`
389  Simple, string name of the filter for the flux being calculated.
390  **kwargs
391  Any additional keyword arguments that may be passed to the plugin.
392  """
393  pTileNames = []
394  for tilePercent in self.configconfig.percentiles:
395  pTileName = "{}PSFluxPercentile{:02d}".format(filterName,
396  tilePercent)
397  pTileNames.append(pTileName)
398  if pTileName not in diaObjects.columns:
399  diaObjects[pTileName] = np.nan
400 
401  def _fluxPercentiles(df):
402  pTiles = np.nanpercentile(df["psFlux"], self.configconfig.percentiles)
403  return pd.Series(
404  dict((tileName, pTile)
405  for tileName, pTile in zip(pTileNames, pTiles)))
406 
407  diaObjects.loc[:, pTileNames] = filterDiaSources.apply(_fluxPercentiles)
408 
409 
411  pass
412 
413 
414 @register("ap_sigmaFlux")
416  """Compute scatter of diaSource fluxes.
417  """
418 
419  ConfigClass = SigmaDiaPsFluxConfig
420  # Output columns are created upon instantiation of the class.
421  outputCols = ["PSFluxSigma"]
422  plugType = "multi"
423  needsFilter = True
424 
425  @classmethod
426  def getExecutionOrder(cls):
427  return cls.DEFAULT_CATALOGCALCULATIONDEFAULT_CATALOGCALCULATION
428 
429  def calculate(self,
430  diaObjects,
431  diaSources,
432  filterDiaSources,
433  filterName,
434  **kwargs):
435  """Compute the sigma fluxes of the point source flux.
436 
437  Parameters
438  ----------
439  diaObject : `dict`
440  Summary object to store values in.
441  diaSources : `pandas.DataFrame`
442  DataFrame representing all diaSources associated with this
443  diaObject.
444  filterDiaSources : `pandas.DataFrame`
445  DataFrame representing diaSources associated with this
446  diaObject that are observed in the band pass ``filterName``.
447  filterName : `str`
448  Simple, string name of the filter for the flux being calculated.
449  **kwargs
450  Any additional keyword arguments that may be passed to the plugin.
451  """
452  # Set "delta degrees of freedom (ddf)" to 1 to calculate the unbiased
453  # estimator of scatter (i.e. 'N - 1' instead of 'N').
454  diaObjects.loc[:, "{}PSFluxSigma".format(filterName)] = \
455  filterDiaSources.psFlux.std()
456 
457 
459  pass
460 
461 
462 @register("ap_chi2Flux")
464  """Compute chi2 of diaSource fluxes.
465  """
466 
467  ConfigClass = Chi2DiaPsFluxConfig
468 
469  # Required input Cols
470  inputCols = ["PSFluxMean"]
471  # Output columns are created upon instantiation of the class.
472  outputCols = ["PSFluxChi2"]
473  plugType = "multi"
474  needsFilter = True
475 
476  @classmethod
477  def getExecutionOrder(cls):
478  return cls.FLUX_MOMENTS_CALCULATEDFLUX_MOMENTS_CALCULATED
479 
480  @catchWarnings(warns=["All-NaN slice encountered"])
481  def calculate(self,
482  diaObjects,
483  diaSources,
484  filterDiaSources,
485  filterName,
486  **kwargs):
487  """Compute the chi2 of the point source fluxes.
488 
489  Parameters
490  ----------
491  diaObject : `dict`
492  Summary object to store values in.
493  diaSources : `pandas.DataFrame`
494  DataFrame representing all diaSources associated with this
495  diaObject.
496  filterDiaSources : `pandas.DataFrame`
497  DataFrame representing diaSources associated with this
498  diaObject that are observed in the band pass ``filterName``.
499  filterName : `str`
500  Simple, string name of the filter for the flux being calculated.
501  **kwargs
502  Any additional keyword arguments that may be passed to the plugin.
503  """
504  meanName = "{}PSFluxMean".format(filterName)
505 
506  def _chi2(df):
507  delta = (df["psFlux"]
508  - diaObjects.at[df.diaObjectId.iat[0], meanName])
509  return np.nansum((delta / df["psFluxErr"]) ** 2)
510 
511  diaObjects.loc[:, "{}PSFluxChi2".format(filterName)] = \
512  filterDiaSources.apply(_chi2)
513 
514 
516  pass
517 
518 
519 @register("ap_madFlux")
521  """Compute median absolute deviation of diaSource fluxes.
522  """
523 
524  ConfigClass = MadDiaPsFluxConfig
525 
526  # Required input Cols
527  # Output columns are created upon instantiation of the class.
528  outputCols = ["PSFluxMAD"]
529  plugType = "multi"
530  needsFilter = True
531 
532  @classmethod
533  def getExecutionOrder(cls):
534  return cls.DEFAULT_CATALOGCALCULATIONDEFAULT_CATALOGCALCULATION
535 
536  @catchWarnings(warns=["All-NaN slice encountered"])
537  def calculate(self,
538  diaObjects,
539  diaSources,
540  filterDiaSources,
541  filterName,
542  **kwargs):
543  """Compute the median absolute deviation of the point source fluxes.
544 
545  Parameters
546  ----------
547  diaObject : `dict`
548  Summary object to store values in.
549  diaSources : `pandas.DataFrame`
550  DataFrame representing all diaSources associated with this
551  diaObject.
552  filterDiaSources : `pandas.DataFrame`
553  DataFrame representing diaSources associated with this
554  diaObject that are observed in the band pass ``filterName``.
555  filterName : `str`
556  Simple, string name of the filter for the flux being calculated.
557  **kwargs
558  Any additional keyword arguments that may be passed to the plugin.
559  """
560  diaObjects.loc[:, "{}PSFluxMAD".format(filterName)] = \
561  filterDiaSources.psFlux.apply(median_absolute_deviation,
562  ignore_nan=True)
563 
564 
566  pass
567 
568 
569 @register("ap_skewFlux")
571  """Compute the skew of diaSource fluxes.
572  """
573 
574  ConfigClass = SkewDiaPsFluxConfig
575 
576  # Required input Cols
577  # Output columns are created upon instantiation of the class.
578  outputCols = ["PSFluxSkew"]
579  plugType = "multi"
580  needsFilter = True
581 
582  @classmethod
583  def getExecutionOrder(cls):
584  return cls.DEFAULT_CATALOGCALCULATIONDEFAULT_CATALOGCALCULATION
585 
586  def calculate(self,
587  diaObjects,
588  diaSources,
589  filterDiaSources,
590  filterName,
591  **kwargs):
592  """Compute the skew of the point source fluxes.
593 
594  Parameters
595  ----------
596  diaObject : `dict`
597  Summary object to store values in.
598  diaSources : `pandas.DataFrame`
599  DataFrame representing all diaSources associated with this
600  diaObject.
601  filterDiaSources : `pandas.DataFrame`
602  DataFrame representing diaSources associated with this
603  diaObject that are observed in the band pass ``filterName``.
604  filterName : `str`
605  Simple, string name of the filter for the flux being calculated.
606  **kwargs
607  Any additional keyword arguments that may be passed to the plugin.
608  """
609  diaObjects.loc[:, "{}PSFluxSkew".format(filterName)] = \
610  filterDiaSources.psFlux.skew()
611 
612 
614  pass
615 
616 
617 @register("ap_minMaxFlux")
619  """Compute min/max of diaSource fluxes.
620  """
621 
622  ConfigClass = MinMaxDiaPsFluxConfig
623 
624  # Required input Cols
625  # Output columns are created upon instantiation of the class.
626  outputCols = ["PSFluxMin", "PSFluxMax"]
627  plugType = "multi"
628  needsFilter = True
629 
630  @classmethod
631  def getExecutionOrder(cls):
632  return cls.DEFAULT_CATALOGCALCULATIONDEFAULT_CATALOGCALCULATION
633 
634  def calculate(self,
635  diaObjects,
636  diaSources,
637  filterDiaSources,
638  filterName,
639  **kwargs):
640  """Compute min/max of the point source fluxes.
641 
642  Parameters
643  ----------
644  diaObject : `dict`
645  Summary object to store values in.
646  diaSources : `pandas.DataFrame`
647  DataFrame representing all diaSources associated with this
648  diaObject.
649  filterDiaSources : `pandas.DataFrame`
650  DataFrame representing diaSources associated with this
651  diaObject that are observed in the band pass ``filterName``.
652  filterName : `str`
653  Simple, string name of the filter for the flux being calculated.
654  **kwargs
655  Any additional keyword arguments that may be passed to the plugin.
656  """
657  minName = "{}PSFluxMin".format(filterName)
658  if minName not in diaObjects.columns:
659  diaObjects[minName] = np.nan
660  maxName = "{}PSFluxMax".format(filterName)
661  if maxName not in diaObjects.columns:
662  diaObjects[maxName] = np.nan
663 
664  diaObjects.loc[:, minName] = filterDiaSources.psFlux.min()
665  diaObjects.loc[:, maxName] = filterDiaSources.psFlux.max()
666 
667 
669  pass
670 
671 
672 @register("ap_maxSlopeFlux")
674  """Compute the maximum ratio time ordered deltaFlux / deltaTime.
675  """
676 
677  ConfigClass = MinMaxDiaPsFluxConfig
678 
679  # Required input Cols
680  # Output columns are created upon instantiation of the class.
681  outputCols = ["PSFluxMaxSlope"]
682  plugType = "multi"
683  needsFilter = True
684 
685  @classmethod
686  def getExecutionOrder(cls):
687  return cls.DEFAULT_CATALOGCALCULATIONDEFAULT_CATALOGCALCULATION
688 
689  def calculate(self,
690  diaObjects,
691  diaSources,
692  filterDiaSources,
693  filterName,
694  **kwargs):
695  """Compute the maximum ratio time ordered deltaFlux / deltaTime.
696 
697  Parameters
698  ----------
699  diaObject : `dict`
700  Summary object to store values in.
701  diaSources : `pandas.DataFrame`
702  DataFrame representing all diaSources associated with this
703  diaObject.
704  filterDiaSources : `pandas.DataFrame`
705  DataFrame representing diaSources associated with this
706  diaObject that are observed in the band pass ``filterName``.
707  filterName : `str`
708  Simple, string name of the filter for the flux being calculated.
709  **kwargs
710  Any additional keyword arguments that may be passed to the plugin.
711  """
712 
713  def _maxSlope(df):
714  tmpDf = df[~np.logical_or(np.isnan(df["psFlux"]),
715  np.isnan(df["midPointTai"]))]
716  if len(tmpDf) < 2:
717  return np.nan
718  times = tmpDf["midPointTai"].to_numpy()
719  timeArgs = times.argsort()
720  times = times[timeArgs]
721  fluxes = tmpDf["psFlux"].to_numpy()[timeArgs]
722  return (np.diff(fluxes) / np.diff(times)).max()
723 
724  diaObjects.loc[:, "{}PSFluxMaxSlope".format(filterName)] = \
725  filterDiaSources.apply(_maxSlope)
726 
727 
729  pass
730 
731 
732 @register("ap_meanErrFlux")
734  """Compute the mean of the dia source errors.
735  """
736 
737  ConfigClass = ErrMeanDiaPsFluxConfig
738 
739  # Required input Cols
740  # Output columns are created upon instantiation of the class.
741  outputCols = ["PSFluxErrMean"]
742  plugType = "multi"
743  needsFilter = True
744 
745  @classmethod
746  def getExecutionOrder(cls):
747  return cls.DEFAULT_CATALOGCALCULATIONDEFAULT_CATALOGCALCULATION
748 
749  def calculate(self,
750  diaObjects,
751  diaSources,
752  filterDiaSources,
753  filterName,
754  **kwargs):
755  """Compute the mean of the dia source errors.
756 
757  Parameters
758  ----------
759  diaObject : `dict`
760  Summary object to store values in.
761  diaSources : `pandas.DataFrame`
762  DataFrame representing all diaSources associated with this
763  diaObject.
764  filterDiaSources : `pandas.DataFrame`
765  DataFrame representing diaSources associated with this
766  diaObject that are observed in the band pass ``filterName``.
767  filterName : `str`
768  Simple, string name of the filter for the flux being calculated.
769  **kwargs
770  Any additional keyword arguments that may be passed to the plugin.
771  """
772  diaObjects.loc[:, "{}PSFluxErrMean".format(filterName)] = \
773  filterDiaSources.psFluxErr.mean()
774 
775 
777  pass
778 
779 
780 @register("ap_linearFit")
782  """Compute fit a linear model to flux vs time.
783  """
784 
785  ConfigClass = LinearFitDiaPsFluxConfig
786 
787  # Required input Cols
788  # Output columns are created upon instantiation of the class.
789  outputCols = ["PSFluxLinearSlope", "PSFluxLinearIntercept"]
790  plugType = "multi"
791  needsFilter = True
792 
793  @classmethod
794  def getExecutionOrder(cls):
795  return cls.DEFAULT_CATALOGCALCULATIONDEFAULT_CATALOGCALCULATION
796 
797  def calculate(self,
798  diaObjects,
799  diaSources,
800  filterDiaSources,
801  filterName,
802  **kwargs):
803  """Compute fit a linear model to flux vs time.
804 
805  Parameters
806  ----------
807  diaObject : `dict`
808  Summary object to store values in.
809  diaSources : `pandas.DataFrame`
810  DataFrame representing all diaSources associated with this
811  diaObject.
812  filterDiaSources : `pandas.DataFrame`
813  DataFrame representing diaSources associated with this
814  diaObject that are observed in the band pass ``filterName``.
815  filterName : `str`
816  Simple, string name of the filter for the flux being calculated.
817  **kwargs
818  Any additional keyword arguments that may be passed to the plugin.
819  """
820 
821  mName = "{}PSFluxLinearSlope".format(filterName)
822  if mName not in diaObjects.columns:
823  diaObjects[mName] = np.nan
824  bName = "{}PSFluxLinearIntercept".format(filterName)
825  if bName not in diaObjects.columns:
826  diaObjects[bName] = np.nan
827 
828  def _linearFit(df):
829  tmpDf = df[~np.logical_or(
830  np.isnan(df["psFlux"]),
831  np.logical_or(np.isnan(df["psFluxErr"]),
832  np.isnan(df["midPointTai"])))]
833  if len(tmpDf) < 2:
834  return pd.Series({mName: np.nan, bName: np.nan})
835  fluxes = tmpDf["psFlux"].to_numpy()
836  errors = tmpDf["psFluxErr"].to_numpy()
837  times = tmpDf["midPointTai"].to_numpy()
838  A = np.array([times / errors, 1 / errors]).transpose()
839  m, b = lsq_linear(A, fluxes / errors).x
840  return pd.Series({mName: m, bName: b})
841 
842  diaObjects.loc[:, [mName, bName]] = filterDiaSources.apply(_linearFit)
843 
844 
846  pass
847 
848 
849 @register("ap_stetsonJ")
851  """Compute the StetsonJ statistic on the DIA point source fluxes.
852  """
853 
854  ConfigClass = LinearFitDiaPsFluxConfig
855 
856  # Required input Cols
857  inputCols = ["PSFluxMean"]
858  # Output columns are created upon instantiation of the class.
859  outputCols = ["PSFluxStetsonJ"]
860  plugType = "multi"
861  needsFilter = True
862 
863  @classmethod
864  def getExecutionOrder(cls):
865  return cls.FLUX_MOMENTS_CALCULATEDFLUX_MOMENTS_CALCULATED
866 
867  def calculate(self,
868  diaObjects,
869  diaSources,
870  filterDiaSources,
871  filterName,
872  **kwargs):
873  """Compute the StetsonJ statistic on the DIA point source fluxes.
874 
875  Parameters
876  ----------
877  diaObject : `dict`
878  Summary object to store values in.
879  diaSources : `pandas.DataFrame`
880  DataFrame representing all diaSources associated with this
881  diaObject.
882  filterDiaSources : `pandas.DataFrame`
883  DataFrame representing diaSources associated with this
884  diaObject that are observed in the band pass ``filterName``.
885  filterName : `str`
886  Simple, string name of the filter for the flux being calculated.
887  **kwargs
888  Any additional keyword arguments that may be passed to the plugin.
889  """
890  meanName = "{}PSFluxMean".format(filterName)
891 
892  def _stetsonJ(df):
893  tmpDf = df[~np.logical_or(np.isnan(df["psFlux"]),
894  np.isnan(df["psFluxErr"]))]
895  if len(tmpDf) < 2:
896  return np.nan
897  fluxes = tmpDf["psFlux"].to_numpy()
898  errors = tmpDf["psFluxErr"].to_numpy()
899 
900  return self._stetson_J_stetson_J(
901  fluxes,
902  errors,
903  diaObjects.at[tmpDf.diaObjectId.iat[0], meanName])
904 
905  diaObjects.loc[:, "{}PSFluxStetsonJ".format(filterName)] = \
906  filterDiaSources.apply(_stetsonJ)
907 
908  def _stetson_J(self, fluxes, errors, mean=None):
909  """Compute the single band stetsonJ statistic.
910 
911  Parameters
912  ----------
913  fluxes : `numpy.ndarray` (N,)
914  Calibrated lightcurve flux values.
915  errors : `numpy.ndarray` (N,)
916  Errors on the calibrated lightcurve fluxes.
917  mean : `float`
918  Starting mean from previous plugin.
919 
920  Returns
921  -------
922  stetsonJ : `float`
923  stetsonJ statistic for the input fluxes and errors.
924 
925  References
926  ----------
927  .. [1] Stetson, P. B., "On the Automatic Determination of Light-Curve
928  Parameters for Cepheid Variables", PASP, 108, 851S, 1996
929  """
930  n_points = len(fluxes)
931  flux_mean = self._stetson_mean_stetson_mean(fluxes, errors, mean)
932  delta_val = (
933  np.sqrt(n_points / (n_points - 1)) * (fluxes - flux_mean) / errors)
934  p_k = delta_val ** 2 - 1
935 
936  return np.mean(np.sign(p_k) * np.sqrt(np.fabs(p_k)))
937 
938  def _stetson_mean(self,
939  values,
940  errors,
941  mean=None,
942  alpha=2.,
943  beta=2.,
944  n_iter=20,
945  tol=1e-6):
946  """Compute the stetson mean of the fluxes which down-weights outliers.
947 
948  Weighted biased on an error weighted difference scaled by a constant
949  (1/``a``) and raised to the power beta. Higher betas more harshly
950  penalize outliers and ``a`` sets the number of sigma where a weighted
951  difference of 1 occurs.
952 
953  Parameters
954  ----------
955  values : `numpy.dnarray`, (N,)
956  Input values to compute the mean of.
957  errors : `numpy.ndarray`, (N,)
958  Errors on the input values.
959  mean : `float`
960  Starting mean value or None.
961  alpha : `float`
962  Scalar down-weighting of the fractional difference. lower->more
963  clipping. (Default value is 2.)
964  beta : `float`
965  Power law slope of the used to down-weight outliers. higher->more
966  clipping. (Default value is 2.)
967  n_iter : `int`
968  Number of iterations of clipping.
969  tol : `float`
970  Fractional and absolute tolerance goal on the change in the mean
971  before exiting early. (Default value is 1e-6)
972 
973  Returns
974  -------
975  mean : `float`
976  Weighted stetson mean result.
977 
978  References
979  ----------
980  .. [1] Stetson, P. B., "On the Automatic Determination of Light-Curve
981  Parameters for Cepheid Variables", PASP, 108, 851S, 1996
982  """
983  n_points = len(values)
984  n_factor = np.sqrt(n_points / (n_points - 1))
985  inv_var = 1 / errors ** 2
986 
987  if mean is None:
988  mean = np.average(values, weights=inv_var)
989  for iter_idx in range(n_iter):
990  chi = np.fabs(n_factor * (values - mean) / errors)
991  tmp_mean = np.average(
992  values,
993  weights=inv_var / (1 + (chi / alpha) ** beta))
994  diff = np.fabs(tmp_mean - mean)
995  mean = tmp_mean
996  if diff / mean < tol and diff < tol:
997  break
998  return mean
999 
1000 
1002  pass
1003 
1004 
1005 @register("ap_meanTotFlux")
1007  """Compute the weighted mean and mean error on the point source fluxes
1008  forced photometered at the DiaSource location in the calibrated image.
1009 
1010  Additionally store number of usable data points.
1011  """
1012 
1013  ConfigClass = WeightedMeanDiaPsFluxConfig
1014  outputCols = ["TOTFluxMean", "TOTFluxMeanErr"]
1015  plugType = "multi"
1016  needsFilter = True
1018  @classmethod
1019  def getExecutionOrder(cls):
1020  return cls.DEFAULT_CATALOGCALCULATIONDEFAULT_CATALOGCALCULATION
1021 
1022  @catchWarnings(warns=["invalid value encountered", "divide by zero"])
1023  def calculate(self,
1024  diaObjects,
1025  diaSources,
1026  filterDiaSources,
1027  filterName,
1028  **kwargs):
1029  """Compute the weighted mean and mean error of the point source flux.
1030 
1031  Parameters
1032  ----------
1033  diaObject : `dict`
1034  Summary object to store values in.
1035  diaSources : `pandas.DataFrame`
1036  DataFrame representing all diaSources associated with this
1037  diaObject.
1038  filterDiaSources : `pandas.DataFrame`
1039  DataFrame representing diaSources associated with this
1040  diaObject that are observed in the band pass ``filterName``.
1041  filterName : `str`
1042  Simple, string name of the filter for the flux being calculated.
1043  **kwargs
1044  Any additional keyword arguments that may be passed to the plugin.
1045  """
1046  totMeanName = "{}TOTFluxMean".format(filterName)
1047  if totMeanName not in diaObjects.columns:
1048  diaObjects[totMeanName] = np.nan
1049  totErrName = "{}TOTFluxMeanErr".format(filterName)
1050  if totErrName not in diaObjects.columns:
1051  diaObjects[totErrName] = np.nan
1052 
1053  def _meanFlux(df):
1054  tmpDf = df[~np.logical_or(np.isnan(df["totFlux"]),
1055  np.isnan(df["totFluxErr"]))]
1056  tot_weight = np.nansum(1 / tmpDf["totFluxErr"] ** 2)
1057  fluxMean = np.nansum(tmpDf["totFlux"]
1058  / tmpDf["totFluxErr"] ** 2)
1059  fluxMean /= tot_weight
1060  fluxMeanErr = np.sqrt(1 / tot_weight)
1061 
1062  return pd.Series({totMeanName: fluxMean,
1063  totErrName: fluxMeanErr})
1064 
1065  diaObjects.loc[:, [totMeanName, totErrName]] = \
1066  filterDiaSources.apply(_meanFlux)
1067 
1068 
1070  pass
1072 
1073 @register("ap_sigmaTotFlux")
1075  """Compute scatter of diaSource fluxes.
1076  """
1077 
1078  ConfigClass = SigmaDiaPsFluxConfig
1079  # Output columns are created upon instantiation of the class.
1080  outputCols = ["TOTFluxSigma"]
1081  plugType = "multi"
1082  needsFilter = True
1084  @classmethod
1085  def getExecutionOrder(cls):
1086  return cls.DEFAULT_CATALOGCALCULATIONDEFAULT_CATALOGCALCULATION
1088  def calculate(self,
1089  diaObjects,
1090  diaSources,
1091  filterDiaSources,
1092  filterName,
1093  **kwargs):
1094  """Compute the sigma fluxes of the point source flux measured on the
1095  calibrated image.
1096 
1097  Parameters
1098  ----------
1099  diaObject : `dict`
1100  Summary object to store values in.
1101  diaSources : `pandas.DataFrame`
1102  DataFrame representing all diaSources associated with this
1103  diaObject.
1104  filterDiaSources : `pandas.DataFrame`
1105  DataFrame representing diaSources associated with this
1106  diaObject that are observed in the band pass ``filterName``.
1107  filterName : `str`
1108  Simple, string name of the filter for the flux being calculated.
1109  **kwargs
1110  Any additional keyword arguments that may be passed to the plugin.
1111  """
1112  # Set "delta degrees of freedom (ddf)" to 1 to calculate the unbiased
1113  # estimator of scatter (i.e. 'N - 1' instead of 'N').
1114  diaObjects.loc[:, "{}TOTFluxSigma".format(filterName)] = \
1115  filterDiaSources.totFlux.std()
1116 
int max
Point in an unspecified spherical coordinate system.
Definition: SpherePoint.h:57
def calculate(self, diaObjects, diaSources, filterDiaSources, filterName, **kwargs)
def calculate(self, diaObjects, diaSources, filterDiaSources, filterName, **kwargs)
def calculate(self, diaObjects, diaObjectId, **kwargs)
def calculate(self, diaObjects, diaSources, filterDiaSources, filterName, **kwargs)
def calculate(self, diaObjects, diaSources, filterDiaSources, filterName, **kwargs)
def calculate(self, diaObjects, diaSources, filterDiaSources, filterName, **kwargs)
def calculate(self, diaObjects, diaSources, **kwargs)
def calculate(self, diaObjects, diaSources, filterDiaSources, filterName, **kwargs)
def calculate(self, diaObjects, diaSources, filterDiaSources, filterName, **kwargs)
def __init__(self, config, name, metadata, **kwargs)
def calculate(self, diaObjects, diaSources, filterDiaSources, filterName, **kwargs)
def calculate(self, diaObjects, diaSources, filterDiaSources, filterName, **kwargs)
def calculate(self, diaObjects, diaSources, filterDiaSources, filterName, **kwargs)
def _stetson_mean(self, values, errors, mean=None, alpha=2., beta=2., n_iter=20, tol=1e-6)
def calculate(self, diaObjects, diaSources, filterDiaSources, filterName, **kwargs)
def calculate(self, diaObjects, diaSources, filterDiaSources, filterName, **kwargs)
def calculate(self, diaObjects, diaSources, filterDiaSources, filterName, **kwargs)
HtmPixelization provides HTM indexing of points and regions.
daf::base::PropertyList * list
Definition: fits.cc:913
SpherePoint averageSpherePoint(std::vector< SpherePoint > const &coords)
Return the average of a list of coordinates.
Definition: SpherePoint.cc:235
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
Definition: history.py:174