Coverage for python/lsst/meas/base/diaCalculationPlugins.py: 40%
537 statements
« prev ^ index » next coverage.py v7.11.0, created at 2025-11-08 01:14 -0800
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/>.
22"""Plugins for use in DiaSource summary statistics.
24Output columns must be
25as defined in the schema of the Apdb both in name and units.
26"""
28import functools
29import warnings
31from astropy.stats import median_absolute_deviation
32import numpy as np
33import pandas as pd
34from scipy.optimize import lsq_linear
36import lsst.geom as geom
37import lsst.pex.config as pexConfig
38import lsst.sphgeom as sphgeom
39from astropy.timeseries import LombScargle
40from astropy.timeseries import LombScargleMultiband
41import math
42import statistics
44from .diaCalculation import (
45 DiaObjectCalculationPluginConfig,
46 DiaObjectCalculationPlugin)
47from .pluginRegistry import register
50__all__ = ("MeanDiaPositionConfig", "MeanDiaPosition",
51 "HTMIndexDiaPosition", "HTMIndexDiaPositionConfig",
52 "NumDiaSourcesDiaPlugin", "NumDiaSourcesDiaPluginConfig",
53 "SimpleSourceFlagDiaPlugin", "SimpleSourceFlagDiaPluginConfig",
54 "WeightedMeanDiaPsfFluxConfig", "WeightedMeanDiaPsfFlux",
55 "PercentileDiaPsfFlux", "PercentileDiaPsfFluxConfig",
56 "SigmaDiaPsfFlux", "SigmaDiaPsfFluxConfig",
57 "Chi2DiaPsfFlux", "Chi2DiaPsfFluxConfig",
58 "MadDiaPsfFlux", "MadDiaPsfFluxConfig",
59 "SkewDiaPsfFlux", "SkewDiaPsfFluxConfig",
60 "MinMaxDiaPsfFlux", "MinMaxDiaPsfFluxConfig",
61 "MaxSlopeDiaPsfFlux", "MaxSlopeDiaPsfFluxConfig",
62 "ErrMeanDiaPsfFlux", "ErrMeanDiaPsfFluxConfig",
63 "LinearFitDiaPsfFlux", "LinearFitDiaPsfFluxConfig",
64 "StetsonJDiaPsfFlux", "StetsonJDiaPsfFluxConfig",
65 "WeightedMeanDiaTotFlux", "WeightedMeanDiaTotFluxConfig",
66 "SigmaDiaTotFlux", "SigmaDiaTotFluxConfig",
67 "LombScarglePeriodogram", "LombScarglePeriodogramConfig",
68 "LombScarglePeriodogramMulti", "LombScarglePeriodogramMultiConfig")
71def catchWarnings(_func=None, *, warns=[]):
72 """Decorator for generically catching numpy warnings.
73 """
74 def decoratorCatchWarnings(func):
75 @functools.wraps(func)
76 def wrapperCatchWarnings(*args, **kwargs):
77 with warnings.catch_warnings():
78 for val in warns:
79 warnings.filterwarnings("ignore", val)
80 return func(*args, **kwargs)
81 return wrapperCatchWarnings
83 if _func is None: 83 ↛ 86line 83 didn't jump to line 86 because the condition on line 83 was always true
84 return decoratorCatchWarnings
85 else:
86 return decoratorCatchWarnings(_func)
89def compute_optimized_periodogram_grid(x0, oversampling_factor=5, nyquist_factor=100):
90 """
91 Computes an optimized periodogram frequency grid for a given time series.
93 Parameters
94 ----------
95 x0 : `array`
96 The input time axis.
97 oversampling_factor : `int`, optional
98 The oversampling factor for frequency grid.
99 nyquist_factor : `int`, optional
100 The Nyquist factor for frequency grid.
102 Returns
103 -------
104 frequencies : `array`
105 The computed optimized periodogram frequency grid.
106 """
108 num_points = len(x0)
109 baseline = np.max(x0) - np.min(x0)
111 # Calculate the frequency resolution based on oversampling factor and baseline
112 frequency_resolution = 1. / baseline / oversampling_factor
114 num_frequencies = int(
115 0.5 * oversampling_factor * nyquist_factor * num_points)
116 frequencies = frequency_resolution + \
117 frequency_resolution * np.arange(num_frequencies)
119 return frequencies
122class LombScarglePeriodogramConfig(DiaObjectCalculationPluginConfig):
123 pass
126@register("ap_lombScarglePeriodogram")
127class LombScarglePeriodogram(DiaObjectCalculationPlugin):
128 """Compute the single-band period of a DiaObject given a set of DiaSources.
129 """
130 ConfigClass = LombScarglePeriodogramConfig
132 plugType = "multi"
133 outputCols = ["period", "power"]
134 needsFilter = True
136 @classmethod
137 def getExecutionOrder(cls):
138 return cls.DEFAULT_CATALOGCALCULATION
140 @catchWarnings(warns=["All-NaN slice encountered"])
141 def calculate(self,
142 diaObjects,
143 diaSources,
144 filterDiaSources,
145 band):
146 """Compute the periodogram.
148 Parameters
149 ----------
150 diaObjects : `pandas.DataFrame`
151 Summary objects to store values in.
152 diaSources : `pandas.DataFrame` or `pandas.DataFrameGroupBy`
153 Catalog of DiaSources summarized by this DiaObject.
154 """
156 # Check and initialize output columns in diaObjects.
157 if (periodCol := f"{band}_period") not in diaObjects.columns:
158 diaObjects[periodCol] = np.nan
159 if (powerCol := f"{band}_power") not in diaObjects.columns:
160 diaObjects[powerCol] = np.nan
162 def _calculate_period(df, min_detections=5, nterms=1, oversampling_factor=5, nyquist_factor=100):
163 """Compute the Lomb-Scargle periodogram given a set of DiaSources.
165 Parameters
166 ----------
167 df : `pandas.DataFrame`
168 The input DataFrame.
169 min_detections : `int`, optional
170 The minimum number of detections.
171 nterms : `int`, optional
172 The number of terms in the Lomb-Scargle model.
173 oversampling_factor : `int`, optional
174 The oversampling factor for frequency grid.
175 nyquist_factor : `int`, optional
176 The Nyquist factor for frequency grid.
178 Returns
179 -------
180 pd_tab : `pandas.Series`
181 The output DataFrame with the Lomb-Scargle parameters.
182 """
183 tmpDf = df[~np.logical_or(np.isnan(df["psfFlux"]),
184 np.isnan(df["midpointMjdTai"]))]
186 if len(tmpDf) < min_detections:
187 return pd.Series({periodCol: np.nan, powerCol: np.nan})
189 time = tmpDf["midpointMjdTai"].to_numpy()
190 flux = tmpDf["psfFlux"].to_numpy()
191 flux_err = tmpDf["psfFluxErr"].to_numpy()
193 lsp = LombScargle(time, flux, dy=flux_err, nterms=nterms)
194 f_grid = compute_optimized_periodogram_grid(
195 time, oversampling_factor=oversampling_factor, nyquist_factor=nyquist_factor)
196 period = 1/f_grid
197 power = lsp.power(f_grid)
199 return pd.Series({periodCol: period[np.argmax(power)],
200 powerCol: np.max(power)})
202 diaObjects.loc[:, [periodCol, powerCol]
203 ] = filterDiaSources.apply(_calculate_period)
206class LombScarglePeriodogramMultiConfig(DiaObjectCalculationPluginConfig):
207 pass
210@register("ap_lombScarglePeriodogramMulti")
211class LombScarglePeriodogramMulti(DiaObjectCalculationPlugin):
212 """Compute the multi-band LombScargle periodogram of a DiaObject given a set of DiaSources.
213 """
214 ConfigClass = LombScarglePeriodogramMultiConfig
216 plugType = "multi"
217 outputCols = ["multiPeriod", "multiPower",
218 "multiFap", "multiAmp", "multiPhase"]
219 needsFilter = True
221 @classmethod
222 def getExecutionOrder(cls):
223 return cls.DEFAULT_CATALOGCALCULATION
225 @staticmethod
226 def calculate_baluev_fap(time, n, maxPeriod, zmax):
227 """Calculate the False-Alarm probability using the Baluev approximation.
229 Parameters
230 ----------
231 time : `array`
232 The input time axis.
233 n : `int`
234 The number of detections.
235 maxPeriod : `float`
236 The maximum period in the grid.
237 zmax : `float`
238 The maximum power in the grid.
240 Returns
241 -------
242 fap_estimate : `float`
243 The False-Alarm probability Baluev approximation.
245 Notes
246 ----------
247 .. [1] Baluev, R. V. 2008, MNRAS, 385, 1279
248 .. [2] Süveges, M., Guy, L.P., Eyer, L., et al. 2015, MNRAS, 450, 2052
249 """
250 if n <= 2:
251 return np.nan
253 gam_ratio = math.factorial(int((n - 1)/2)) / math.factorial(int((n - 2)/2))
254 fu = 1/maxPeriod
255 return gam_ratio * np.sqrt(
256 4*np.pi*statistics.variance(time)
257 ) * fu * (1-zmax)**((n-4)/2) * np.sqrt(zmax)
259 @staticmethod
260 def generate_lsp_params(lsp_model, fbest, bands):
261 """Generate the Lomb-Scargle parameters.
262 Parameters
263 ----------
264 lsp_model : `astropy.timeseries.LombScargleMultiband`
265 The Lomb-Scargle model.
266 fbest : `float`
267 The best period.
268 bands : `array`
269 The bands of the time series.
271 Returns
272 -------
273 Amp : `array`
274 The amplitude of the time series.
275 Ph : `array`
276 The phase of the time series.
278 Notes
279 ----------
280 .. [1] VanderPlas, J. T., & Ivezic, Z. 2015, ApJ, 812, 18
281 """
282 best_params = lsp_model.model_parameters(fbest, units=True)
284 name_params = [f"theta_base_{i}" for i in range(3)]
285 name_params += [f"theta_band_{band}_{i}" for band in np.unique(bands) for i in range(3)]
287 df_params = pd.DataFrame([best_params], columns=name_params)
289 unique_bands = np.unique(bands)
291 amplitude_band = [np.sqrt(df_params[f"theta_band_{band}_1"]**2
292 + df_params[f"theta_band_{band}_2"]**2)
293 for band in unique_bands]
294 phase_bands = [np.arctan2(df_params[f"theta_band_{band}_2"],
295 df_params[f"theta_band_{band}_1"]) for band in unique_bands]
297 amp = [a[0] for a in amplitude_band]
298 ph = [p[0] for p in phase_bands]
300 return amp, ph
302 @catchWarnings(warns=["All-NaN slice encountered"])
303 def calculate(self,
304 diaObjects,
305 diaSources,
306 **kwargs):
307 """Compute the multi-band LombScargle periodogram of a DiaObject given
308 a set of DiaSources.
310 Parameters
311 ----------
312 diaObjects : `pandas.DataFrame`
313 Summary objects to store values in.
314 diaSources : `pandas.DataFrame` or `pandas.DataFrameGroupBy`
315 Catalog of DiaSources summarized by this DiaObject.
316 **kwargs : `dict`
317 Unused kwargs that are always passed to a plugin.
318 """
320 bands_arr = diaSources['band'].unique().values
321 unique_bands = np.unique(np.concatenate(bands_arr))
322 # Check and initialize output columns in diaObjects.
323 if (periodCol := "multiPeriod") not in diaObjects.columns:
324 diaObjects[periodCol] = np.nan
325 if (powerCol := "multiPower") not in diaObjects.columns:
326 diaObjects[powerCol] = np.nan
327 if (fapCol := "multiFap") not in diaObjects.columns:
328 diaObjects[fapCol] = np.nan
329 ampCol = "multiAmp"
330 phaseCol = "multiPhase"
331 for i in range(len(unique_bands)):
332 ampCol_band = f"{unique_bands[i]}_{ampCol}"
333 if ampCol_band not in diaObjects.columns:
334 diaObjects[ampCol_band] = np.nan
335 phaseCol_band = f"{unique_bands[i]}_{phaseCol}"
336 if phaseCol_band not in diaObjects.columns:
337 diaObjects[phaseCol_band] = np.nan
339 def _calculate_period_multi(df, all_unique_bands,
340 min_detections=9, oversampling_factor=5, nyquist_factor=100):
341 """Calculate the multi-band Lomb-Scargle periodogram.
343 Parameters
344 ----------
345 df : `pandas.DataFrame`
346 The input DataFrame.
347 all_unique_bands : `list` of `str`
348 List of all bands present in the diaSource table that is being worked on.
349 min_detections : `int`, optional
350 The minimum number of detections, including all bands.
351 oversampling_factor : `int`, optional
352 The oversampling factor for frequency grid.
353 nyquist_factor : `int`, optional
354 The Nyquist factor for frequency grid.
356 Returns
357 -------
358 pd_tab : `pandas.Series`
359 The output DataFrame with the Lomb-Scargle parameters.
360 """
361 tmpDf = df[~np.logical_or(np.isnan(df["psfFlux"]),
362 np.isnan(df["midpointMjdTai"]))]
364 if (len(tmpDf)) < min_detections:
365 pd_tab_nodet = pd.Series({periodCol: np.nan,
366 powerCol: np.nan,
367 fapCol: np.nan})
368 for band in all_unique_bands:
369 pd_tab_nodet[f"{band}_{ampCol}"] = np.nan
370 pd_tab_nodet[f"{band}_{phaseCol}"] = np.nan
372 return pd_tab_nodet
374 time = tmpDf["midpointMjdTai"].to_numpy()
375 flux = tmpDf["psfFlux"].to_numpy()
376 flux_err = tmpDf["psfFluxErr"].to_numpy()
377 bands = tmpDf["band"].to_numpy()
379 lsp = LombScargleMultiband(time, flux, bands, dy=flux_err,
380 nterms_base=1, nterms_band=1)
382 f_grid = compute_optimized_periodogram_grid(
383 time, oversampling_factor=oversampling_factor, nyquist_factor=nyquist_factor)
384 period = 1/f_grid
385 power = lsp.power(f_grid)
387 fap_estimate = self.calculate_baluev_fap(
388 time, len(time), period[np.argmax(power)], np.max(power))
390 params_table_new = self.generate_lsp_params(lsp, f_grid[np.argmax(power)], bands)
392 pd_tab = pd.Series({periodCol: period[np.argmax(power)],
393 powerCol: np.max(power),
394 fapCol: fap_estimate
395 })
397 # Initialize the per-band amplitude/phase columns as NaNs
398 for band in all_unique_bands:
399 pd_tab[f"{band}_{ampCol}"] = np.nan
400 pd_tab[f"{band}_{phaseCol}"] = np.nan
402 # Populate the values of only the bands that have data for this diaSource
403 unique_bands = np.unique(bands)
404 for i in range(len(unique_bands)):
405 pd_tab[f"{unique_bands[i]}_{ampCol}"] = params_table_new[0][i]
406 pd_tab[f"{unique_bands[i]}_{phaseCol}"] = params_table_new[1][i]
408 return pd_tab
410 columns_list = [periodCol, powerCol, fapCol]
411 for i in range(len(unique_bands)):
412 columns_list.append(f"{unique_bands[i]}_{ampCol}")
413 columns_list.append(f"{unique_bands[i]}_{phaseCol}")
415 diaObjects.loc[:, columns_list
416 ] = diaSources.apply(_calculate_period_multi, unique_bands)
419class MeanDiaPositionConfig(DiaObjectCalculationPluginConfig):
420 pass
423@register("ap_meanPosition")
424class MeanDiaPosition(DiaObjectCalculationPlugin):
425 """Compute the mean position of a DiaObject given a set of DiaSources.
426 """
428 ConfigClass = MeanDiaPositionConfig
430 plugType = 'multi'
432 outputCols = ["ra", "dec"]
433 needsFilter = False
435 @classmethod
436 def getExecutionOrder(cls):
437 return cls.DEFAULT_CATALOGCALCULATION
439 def calculate(self, diaObjects, diaSources, **kwargs):
440 """Compute the mean ra/dec position of the diaObject given the
441 diaSource locations.
443 Parameters
444 ----------
445 diaObjects : `pandas.DataFrame`
446 Summary objects to store values in.
447 diaSources : `pandas.DataFrame` or `pandas.DataFrameGroupBy`
448 Catalog of DiaSources summarized by this DiaObject.
449 **kwargs
450 Any additional keyword arguments that may be passed to the plugin.
451 """
452 for outCol in self.outputCols:
453 if outCol not in diaObjects.columns:
454 diaObjects[outCol] = np.nan
456 def _computeMeanPos(df):
457 aveCoord = geom.averageSpherePoint(
458 list(geom.SpherePoint(src["ra"], src["dec"], geom.degrees)
459 for idx, src in df.iterrows()))
461 return pd.Series({"ra": aveCoord.getRa().asDegrees(),
462 "dec": aveCoord.getDec().asDegrees()})
464 ans = diaSources.apply(_computeMeanPos)
465 diaObjects.loc[:, ["ra", "dec"]] = ans
468class HTMIndexDiaPositionConfig(DiaObjectCalculationPluginConfig):
470 htmLevel = pexConfig.Field(
471 dtype=int,
472 doc="Level of the HTM pixelization.",
473 default=20,
474 )
477@register("ap_HTMIndex")
478class HTMIndexDiaPosition(DiaObjectCalculationPlugin):
479 """Compute the mean position of a DiaObject given a set of DiaSources.
481 Notes
482 -----
483 This plugin was implemented to satisfy requirements of old APDB interface
484 which required ``pixelId`` column in DiaObject with HTM20 index. APDB
485 interface had migrated to not need that information, but we keep this
486 plugin in case it may be useful for something else.
487 """
488 ConfigClass = HTMIndexDiaPositionConfig
490 plugType = 'single'
492 inputCols = ["ra", "dec"]
493 outputCols = ["pixelId"]
494 needsFilter = False
496 def __init__(self, config, name, metadata):
497 DiaObjectCalculationPlugin.__init__(self, config, name, metadata)
498 self.pixelator = sphgeom.HtmPixelization(self.config.htmLevel)
500 @classmethod
501 def getExecutionOrder(cls):
502 return cls.FLUX_MOMENTS_CALCULATED
504 def calculate(self, diaObjects, diaObjectId, **kwargs):
505 """Compute the mean position of a DiaObject given a set of DiaSources
507 Parameters
508 ----------
509 diaObjects : `pandas.dataFrame`
510 Summary objects to store values in and read ra/dec from.
511 diaObjectId : `int`
512 Id of the diaObject to update.
513 **kwargs
514 Any additional keyword arguments that may be passed to the plugin.
515 """
516 sphPoint = geom.SpherePoint(
517 diaObjects.at[diaObjectId, "ra"] * geom.degrees,
518 diaObjects.at[diaObjectId, "dec"] * geom.degrees)
519 diaObjects.at[diaObjectId, "pixelId"] = self.pixelator.index(
520 sphPoint.getVector())
523class NumDiaSourcesDiaPluginConfig(DiaObjectCalculationPluginConfig):
524 pass
527@register("ap_nDiaSources")
528class NumDiaSourcesDiaPlugin(DiaObjectCalculationPlugin):
529 """Compute the total number of DiaSources associated with this DiaObject.
530 """
532 ConfigClass = NumDiaSourcesDiaPluginConfig
533 outputCols = ["nDiaSources"]
534 plugType = "multi"
535 needsFilter = False
537 @classmethod
538 def getExecutionOrder(cls):
539 return cls.DEFAULT_CATALOGCALCULATION
541 def calculate(self, diaObjects, diaSources, **kwargs):
542 """Compute the total number of DiaSources associated with this DiaObject.
544 Parameters
545 ----------
546 diaObject : `dict`
547 Summary object to store values in and read ra/dec from.
548 **kwargs
549 Any additional keyword arguments that may be passed to the plugin.
550 """
551 dtype = diaObjects["nDiaSources"].dtype
552 diaObjects.loc[:, "nDiaSources"] = diaSources.diaObjectId.count().astype(dtype)
555class SimpleSourceFlagDiaPluginConfig(DiaObjectCalculationPluginConfig):
556 pass
559@register("ap_diaObjectFlag")
560class SimpleSourceFlagDiaPlugin(DiaObjectCalculationPlugin):
561 """Find if any DiaSource is flagged.
563 Set the DiaObject flag if any DiaSource is flagged.
564 """
566 ConfigClass = NumDiaSourcesDiaPluginConfig
567 outputCols = ["flags"]
568 plugType = "multi"
569 needsFilter = False
571 @classmethod
572 def getExecutionOrder(cls):
573 return cls.DEFAULT_CATALOGCALCULATION
575 def calculate(self, diaObjects, diaSources, **kwargs):
576 """Find if any DiaSource is flagged.
578 Set the DiaObject flag if any DiaSource is flagged.
580 Parameters
581 ----------
582 diaObject : `dict`
583 Summary object to store values in and read ra/dec from.
584 **kwargs
585 Any additional keyword arguments that may be passed to the plugin.
586 """
587 dtype = diaObjects["flags"].dtype
588 diaObjects.loc[:, "flags"] = diaSources.flags.any().astype(dtype)
591class WeightedMeanDiaPsfFluxConfig(DiaObjectCalculationPluginConfig):
592 pass
595@register("ap_meanFlux")
596class WeightedMeanDiaPsfFlux(DiaObjectCalculationPlugin):
597 """Compute the weighted mean and mean error on the point source fluxes
598 of the DiaSource measured on the difference image.
600 Additionally store number of usable data points.
601 """
603 ConfigClass = WeightedMeanDiaPsfFluxConfig
604 outputCols = ["psfFluxMean", "psfFluxMeanErr", "psfFluxNdata"]
605 plugType = "multi"
606 needsFilter = True
608 @classmethod
609 def getExecutionOrder(cls):
610 return cls.DEFAULT_CATALOGCALCULATION
612 @catchWarnings(warns=["invalid value encountered",
613 "divide by zero"])
614 def calculate(self,
615 diaObjects,
616 diaSources,
617 filterDiaSources,
618 band,
619 **kwargs):
620 """Compute the weighted mean and mean error of the point source flux.
622 Parameters
623 ----------
624 diaObject : `dict`
625 Summary object to store values in.
626 diaSources : `pandas.DataFrame`
627 DataFrame representing all diaSources associated with this
628 diaObject.
629 filterDiaSources : `pandas.DataFrame`
630 DataFrame representing diaSources associated with this
631 diaObject that are observed in the band pass ``band``.
632 band : `str`
633 Simple, string name of the filter for the flux being calculated.
634 **kwargs
635 Any additional keyword arguments that may be passed to the plugin.
636 """
637 meanName = "{}_psfFluxMean".format(band)
638 errName = "{}_psfFluxMeanErr".format(band)
639 nDataName = "{}_psfFluxNdata".format(band)
640 if meanName not in diaObjects.columns:
641 diaObjects[meanName] = np.nan
642 if errName not in diaObjects.columns:
643 diaObjects[errName] = np.nan
644 if nDataName not in diaObjects.columns:
645 diaObjects[nDataName] = 0
647 def _weightedMean(df):
648 tmpDf = df[~np.logical_or(np.isnan(df["psfFlux"]),
649 np.isnan(df["psfFluxErr"]))]
650 tot_weight = np.nansum(1 / tmpDf["psfFluxErr"] ** 2)
651 fluxMean = np.nansum(tmpDf["psfFlux"]
652 / tmpDf["psfFluxErr"] ** 2)
653 fluxMean /= tot_weight
654 if tot_weight > 0:
655 fluxMeanErr = np.sqrt(1 / tot_weight)
656 else:
657 fluxMeanErr = np.nan
658 nFluxData = len(tmpDf)
660 return pd.Series({meanName: fluxMean,
661 errName: fluxMeanErr,
662 nDataName: nFluxData},
663 dtype="object")
664 df = filterDiaSources.apply(_weightedMean).astype(diaObjects.dtypes[[meanName, errName, nDataName]])
666 diaObjects.loc[:, [meanName, errName, nDataName]] = df
669class PercentileDiaPsfFluxConfig(DiaObjectCalculationPluginConfig):
670 percentiles = pexConfig.ListField(
671 dtype=int,
672 default=[5, 25, 50, 75, 95],
673 doc="Percentiles to calculate to compute values for. Should be "
674 "integer values."
675 )
678@register("ap_percentileFlux")
679class PercentileDiaPsfFlux(DiaObjectCalculationPlugin):
680 """Compute percentiles of diaSource fluxes.
681 """
683 ConfigClass = PercentileDiaPsfFluxConfig
684 # Output columns are created upon instantiation of the class.
685 outputCols = []
686 plugType = "multi"
687 needsFilter = True
689 def __init__(self, config, name, metadata, **kwargs):
690 DiaObjectCalculationPlugin.__init__(self,
691 config,
692 name,
693 metadata,
694 **kwargs)
695 self.outputCols = ["psfFluxPercentile{:02d}".format(percent)
696 for percent in self.config.percentiles]
698 @classmethod
699 def getExecutionOrder(cls):
700 return cls.DEFAULT_CATALOGCALCULATION
702 @catchWarnings(warns=["All-NaN slice encountered"])
703 def calculate(self,
704 diaObjects,
705 diaSources,
706 filterDiaSources,
707 band,
708 **kwargs):
709 """Compute the percentile fluxes of the point source flux.
711 Parameters
712 ----------
713 diaObject : `dict`
714 Summary object to store values in.
715 diaSources : `pandas.DataFrame`
716 DataFrame representing all diaSources associated with this
717 diaObject.
718 filterDiaSources : `pandas.DataFrame`
719 DataFrame representing diaSources associated with this
720 diaObject that are observed in the band pass ``band``.
721 band : `str`
722 Simple, string name of the filter for the flux being calculated.
723 **kwargs
724 Any additional keyword arguments that may be passed to the plugin.
725 """
726 pTileNames = []
727 dtype = None
728 for tilePercent in self.config.percentiles:
729 pTileName = "{}_psfFluxPercentile{:02d}".format(band,
730 tilePercent)
731 pTileNames.append(pTileName)
732 if pTileName not in diaObjects.columns:
733 diaObjects[pTileName] = np.nan
734 elif dtype is None:
735 dtype = diaObjects[pTileName].dtype
737 def _fluxPercentiles(df):
738 pTiles = np.nanpercentile(df["psfFlux"], self.config.percentiles)
739 return pd.Series(
740 dict((tileName, pTile)
741 for tileName, pTile in zip(pTileNames, pTiles)))
742 if dtype is None:
743 diaObjects.loc[:, pTileNames] = filterDiaSources.apply(_fluxPercentiles)
744 else:
745 diaObjects.loc[:, pTileNames] = filterDiaSources.apply(_fluxPercentiles).astype(dtype)
748class SigmaDiaPsfFluxConfig(DiaObjectCalculationPluginConfig):
749 pass
752@register("ap_sigmaFlux")
753class SigmaDiaPsfFlux(DiaObjectCalculationPlugin):
754 """Compute scatter of diaSource fluxes.
755 """
757 ConfigClass = SigmaDiaPsfFluxConfig
758 # Output columns are created upon instantiation of the class.
759 outputCols = ["psfFluxSigma"]
760 plugType = "multi"
761 needsFilter = True
763 @classmethod
764 def getExecutionOrder(cls):
765 return cls.DEFAULT_CATALOGCALCULATION
767 def calculate(self,
768 diaObjects,
769 diaSources,
770 filterDiaSources,
771 band,
772 **kwargs):
773 """Compute the sigma fluxes of the point source flux.
775 Parameters
776 ----------
777 diaObject : `dict`
778 Summary object to store values in.
779 diaSources : `pandas.DataFrame`
780 DataFrame representing all diaSources associated with this
781 diaObject.
782 filterDiaSources : `pandas.DataFrame`
783 DataFrame representing diaSources associated with this
784 diaObject that are observed in the band pass ``band``.
785 band : `str`
786 Simple, string name of the filter for the flux being calculated.
787 **kwargs
788 Any additional keyword arguments that may be passed to the plugin.
789 """
790 # Set "delta degrees of freedom (ddf)" to 1 to calculate the unbiased
791 # estimator of scatter (i.e. 'N - 1' instead of 'N').
792 column = "{}_psfFluxSigma".format(band)
793 if column in diaObjects:
794 dtype = diaObjects[column].dtype
795 diaObjects.loc[:, column] = filterDiaSources.psfFlux.std().astype(dtype)
796 else:
797 # std() is inconsistent about what type it returns; force to double.
798 diaObjects.loc[:, column] = filterDiaSources.psfFlux.std().astype(np.float64)
801class Chi2DiaPsfFluxConfig(DiaObjectCalculationPluginConfig):
802 pass
805@register("ap_chi2Flux")
806class Chi2DiaPsfFlux(DiaObjectCalculationPlugin):
807 """Compute chi2 of diaSource fluxes.
808 """
810 ConfigClass = Chi2DiaPsfFluxConfig
812 # Required input Cols
813 inputCols = ["psfFluxMean"]
814 # Output columns are created upon instantiation of the class.
815 outputCols = ["psfFluxChi2"]
816 plugType = "multi"
817 needsFilter = True
819 @classmethod
820 def getExecutionOrder(cls):
821 return cls.FLUX_MOMENTS_CALCULATED
823 @catchWarnings(warns=["All-NaN slice encountered"])
824 def calculate(self,
825 diaObjects,
826 diaSources,
827 filterDiaSources,
828 band,
829 **kwargs):
830 """Compute the chi2 of the point source fluxes.
832 Parameters
833 ----------
834 diaObject : `dict`
835 Summary object to store values in.
836 diaSources : `pandas.DataFrame`
837 DataFrame representing all diaSources associated with this
838 diaObject.
839 filterDiaSources : `pandas.DataFrame`
840 DataFrame representing diaSources associated with this
841 diaObject that are observed in the band pass ``band``.
842 band : `str`
843 Simple, string name of the filter for the flux being calculated.
844 **kwargs
845 Any additional keyword arguments that may be passed to the plugin.
846 """
847 meanName = "{}_psfFluxMean".format(band)
848 column = "{}_psfFluxChi2".format(band)
850 def _chi2(df):
851 delta = (df["psfFlux"]
852 - diaObjects.at[df.diaObjectId.iat[0], meanName])
853 return np.nansum((delta / df["psfFluxErr"]) ** 2)
855 if column in diaObjects:
856 dtype = diaObjects[column].dtype
857 diaObjects.loc[:, column] = filterDiaSources.apply(_chi2).astype(dtype)
858 else:
859 diaObjects.loc[:, column] = filterDiaSources.apply(_chi2)
862class MadDiaPsfFluxConfig(DiaObjectCalculationPluginConfig):
863 pass
866@register("ap_madFlux")
867class MadDiaPsfFlux(DiaObjectCalculationPlugin):
868 """Compute median absolute deviation of diaSource fluxes.
869 """
871 ConfigClass = MadDiaPsfFluxConfig
873 # Required input Cols
874 # Output columns are created upon instantiation of the class.
875 outputCols = ["psfFluxMAD"]
876 plugType = "multi"
877 needsFilter = True
879 @classmethod
880 def getExecutionOrder(cls):
881 return cls.DEFAULT_CATALOGCALCULATION
883 @catchWarnings(warns=["All-NaN slice encountered"])
884 def calculate(self,
885 diaObjects,
886 diaSources,
887 filterDiaSources,
888 band,
889 **kwargs):
890 """Compute the median absolute deviation of the point source fluxes.
892 Parameters
893 ----------
894 diaObject : `dict`
895 Summary object to store values in.
896 diaSources : `pandas.DataFrame`
897 DataFrame representing all diaSources associated with this
898 diaObject.
899 filterDiaSources : `pandas.DataFrame`
900 DataFrame representing diaSources associated with this
901 diaObject that are observed in the band pass ``band``.
902 band : `str`
903 Simple, string name of the filter for the flux being calculated.
904 **kwargs
905 Any additional keyword arguments that may be passed to the plugin.
906 """
907 column = "{}_psfFluxMAD".format(band)
908 if column in diaObjects:
909 dtype = diaObjects[column].dtype
910 diaObjects.loc[:, column] = \
911 filterDiaSources.psfFlux.apply(median_absolute_deviation,
912 ignore_nan=True).astype(dtype)
913 else:
914 diaObjects.loc[:, column] = \
915 filterDiaSources.psfFlux.apply(median_absolute_deviation,
916 ignore_nan=True)
919class SkewDiaPsfFluxConfig(DiaObjectCalculationPluginConfig):
920 pass
923@register("ap_skewFlux")
924class SkewDiaPsfFlux(DiaObjectCalculationPlugin):
925 """Compute the skew of diaSource fluxes.
926 """
928 ConfigClass = SkewDiaPsfFluxConfig
930 # Required input Cols
931 # Output columns are created upon instantiation of the class.
932 outputCols = ["psfFluxSkew"]
933 plugType = "multi"
934 needsFilter = True
936 @classmethod
937 def getExecutionOrder(cls):
938 return cls.DEFAULT_CATALOGCALCULATION
940 def calculate(self,
941 diaObjects,
942 diaSources,
943 filterDiaSources,
944 band,
945 **kwargs):
946 """Compute the skew of the point source fluxes.
948 Parameters
949 ----------
950 diaObject : `dict`
951 Summary object to store values in.
952 diaSources : `pandas.DataFrame`
953 DataFrame representing all diaSources associated with this
954 diaObject.
955 filterDiaSources : `pandas.DataFrame`
956 DataFrame representing diaSources associated with this
957 diaObject that are observed in the band pass ``band``.
958 band : `str`
959 Simple, string name of the filter for the flux being calculated.
960 **kwargs
961 Any additional keyword arguments that may be passed to the plugin.
962 """
963 column = "{}_psfFluxSkew".format(band)
964 if column in diaObjects:
965 dtype = diaObjects[column].dtype
966 diaObjects.loc[:, column] = filterDiaSources.psfFlux.skew().astype(dtype)
967 else:
968 diaObjects.loc[:, column] = filterDiaSources.psfFlux.skew()
971class MinMaxDiaPsfFluxConfig(DiaObjectCalculationPluginConfig):
972 pass
975@register("ap_minMaxFlux")
976class MinMaxDiaPsfFlux(DiaObjectCalculationPlugin):
977 """Compute min/max of diaSource fluxes.
978 """
980 ConfigClass = MinMaxDiaPsfFluxConfig
982 # Required input Cols
983 # Output columns are created upon instantiation of the class.
984 outputCols = ["psfFluxMin", "psfFluxMax"]
985 plugType = "multi"
986 needsFilter = True
988 @classmethod
989 def getExecutionOrder(cls):
990 return cls.DEFAULT_CATALOGCALCULATION
992 def calculate(self,
993 diaObjects,
994 diaSources,
995 filterDiaSources,
996 band,
997 **kwargs):
998 """Compute min/max of the point source fluxes.
1000 Parameters
1001 ----------
1002 diaObject : `dict`
1003 Summary object to store values in.
1004 diaSources : `pandas.DataFrame`
1005 DataFrame representing all diaSources associated with this
1006 diaObject.
1007 filterDiaSources : `pandas.DataFrame`
1008 DataFrame representing diaSources associated with this
1009 diaObject that are observed in the band pass ``band``.
1010 band : `str`
1011 Simple, string name of the filter for the flux being calculated.
1012 **kwargs
1013 Any additional keyword arguments that may be passed to the plugin.
1014 """
1015 minName = "{}_psfFluxMin".format(band)
1016 if minName not in diaObjects.columns:
1017 diaObjects[minName] = np.nan
1018 maxName = "{}_psfFluxMax".format(band)
1019 if maxName not in diaObjects.columns:
1020 diaObjects[maxName] = np.nan
1022 dtype = diaObjects[minName].dtype
1023 diaObjects.loc[:, minName] = filterDiaSources.psfFlux.min().astype(dtype)
1024 diaObjects.loc[:, maxName] = filterDiaSources.psfFlux.max().astype(dtype)
1027class MaxSlopeDiaPsfFluxConfig(DiaObjectCalculationPluginConfig):
1028 pass
1031@register("ap_maxSlopeFlux")
1032class MaxSlopeDiaPsfFlux(DiaObjectCalculationPlugin):
1033 """Compute the maximum ratio time ordered deltaFlux / deltaTime.
1034 """
1036 ConfigClass = MinMaxDiaPsfFluxConfig
1038 # Required input Cols
1039 # Output columns are created upon instantiation of the class.
1040 outputCols = ["psfFluxMaxSlope"]
1041 plugType = "multi"
1042 needsFilter = True
1044 @classmethod
1045 def getExecutionOrder(cls):
1046 return cls.DEFAULT_CATALOGCALCULATION
1048 def calculate(self,
1049 diaObjects,
1050 diaSources,
1051 filterDiaSources,
1052 band,
1053 **kwargs):
1054 """Compute the maximum ratio time ordered deltaFlux / deltaTime.
1056 Parameters
1057 ----------
1058 diaObject : `dict`
1059 Summary object to store values in.
1060 diaSources : `pandas.DataFrame`
1061 DataFrame representing all diaSources associated with this
1062 diaObject.
1063 filterDiaSources : `pandas.DataFrame`
1064 DataFrame representing diaSources associated with this
1065 diaObject that are observed in the band pass ``band``.
1066 band : `str`
1067 Simple, string name of the filter for the flux being calculated.
1068 **kwargs
1069 Any additional keyword arguments that may be passed to the plugin.
1070 """
1072 def _maxSlope(df):
1073 tmpDf = df[~np.logical_or(np.isnan(df["psfFlux"]),
1074 np.isnan(df["midpointMjdTai"]))]
1075 if len(tmpDf) < 2:
1076 return np.nan
1077 times = tmpDf["midpointMjdTai"].to_numpy()
1078 timeArgs = times.argsort()
1079 times = times[timeArgs]
1080 fluxes = tmpDf["psfFlux"].to_numpy()[timeArgs]
1081 return (np.diff(fluxes) / np.diff(times)).max()
1083 column = "{}_psfFluxMaxSlope".format(band)
1084 if column in diaObjects:
1085 dtype = diaObjects[column].dtype
1086 diaObjects.loc[:, column] = filterDiaSources.apply(_maxSlope).astype(dtype)
1087 else:
1088 diaObjects.loc[:, column] = filterDiaSources.apply(_maxSlope)
1091class ErrMeanDiaPsfFluxConfig(DiaObjectCalculationPluginConfig):
1092 pass
1095@register("ap_meanErrFlux")
1096class ErrMeanDiaPsfFlux(DiaObjectCalculationPlugin):
1097 """Compute the mean of the dia source errors.
1098 """
1100 ConfigClass = ErrMeanDiaPsfFluxConfig
1102 # Required input Cols
1103 # Output columns are created upon instantiation of the class.
1104 outputCols = ["psfFluxErrMean"]
1105 plugType = "multi"
1106 needsFilter = True
1108 @classmethod
1109 def getExecutionOrder(cls):
1110 return cls.DEFAULT_CATALOGCALCULATION
1112 def calculate(self,
1113 diaObjects,
1114 diaSources,
1115 filterDiaSources,
1116 band,
1117 **kwargs):
1118 """Compute the mean of the dia source errors.
1120 Parameters
1121 ----------
1122 diaObject : `dict`
1123 Summary object to store values in.
1124 diaSources : `pandas.DataFrame`
1125 DataFrame representing all diaSources associated with this
1126 diaObject.
1127 filterDiaSources : `pandas.DataFrame`
1128 DataFrame representing diaSources associated with this
1129 diaObject that are observed in the band pass ``band``.
1130 band : `str`
1131 Simple, string name of the filter for the flux being calculated.
1132 **kwargs
1133 Any additional keyword arguments that may be passed to the plugin.
1134 """
1135 column = "{}_psfFluxErrMean".format(band)
1136 if column in diaObjects:
1137 dtype = diaObjects[column].dtype
1138 diaObjects.loc[:, column] = filterDiaSources.psfFluxErr.mean().astype(dtype)
1139 else:
1140 diaObjects.loc[:, column] = filterDiaSources.psfFluxErr.mean()
1143class LinearFitDiaPsfFluxConfig(DiaObjectCalculationPluginConfig):
1144 pass
1147@register("ap_linearFit")
1148class LinearFitDiaPsfFlux(DiaObjectCalculationPlugin):
1149 """Compute fit a linear model to flux vs time.
1150 """
1152 ConfigClass = LinearFitDiaPsfFluxConfig
1154 # Required input Cols
1155 # Output columns are created upon instantiation of the class.
1156 outputCols = ["psfFluxLinearSlope", "psfFluxLinearIntercept"]
1157 plugType = "multi"
1158 needsFilter = True
1160 @classmethod
1161 def getExecutionOrder(cls):
1162 return cls.DEFAULT_CATALOGCALCULATION
1164 def calculate(self,
1165 diaObjects,
1166 diaSources,
1167 filterDiaSources,
1168 band,
1169 **kwargs):
1170 """Compute fit a linear model to flux vs time.
1172 Parameters
1173 ----------
1174 diaObject : `dict`
1175 Summary object to store values in.
1176 diaSources : `pandas.DataFrame`
1177 DataFrame representing all diaSources associated with this
1178 diaObject.
1179 filterDiaSources : `pandas.DataFrame`
1180 DataFrame representing diaSources associated with this
1181 diaObject that are observed in the band pass ``band``.
1182 band : `str`
1183 Simple, string name of the filter for the flux being calculated.
1184 **kwargs
1185 Any additional keyword arguments that may be passed to the plugin.
1186 """
1188 mName = "{}_psfFluxLinearSlope".format(band)
1189 if mName not in diaObjects.columns:
1190 diaObjects[mName] = np.nan
1191 bName = "{}_psfFluxLinearIntercept".format(band)
1192 if bName not in diaObjects.columns:
1193 diaObjects[bName] = np.nan
1194 dtype = diaObjects[mName].dtype
1196 def _linearFit(df):
1197 tmpDf = df[~np.logical_or(
1198 np.isnan(df["psfFlux"]),
1199 np.logical_or(np.isnan(df["psfFluxErr"]),
1200 np.isnan(df["midpointMjdTai"])))]
1201 if len(tmpDf) < 2:
1202 return pd.Series({mName: np.nan, bName: np.nan})
1203 fluxes = tmpDf["psfFlux"].to_numpy()
1204 errors = tmpDf["psfFluxErr"].to_numpy()
1205 times = tmpDf["midpointMjdTai"].to_numpy()
1206 A = np.array([times / errors, 1 / errors]).transpose()
1207 m, b = lsq_linear(A, fluxes / errors).x
1208 return pd.Series({mName: m, bName: b}, dtype=dtype)
1210 diaObjects.loc[:, [mName, bName]] = filterDiaSources.apply(_linearFit)
1213class StetsonJDiaPsfFluxConfig(DiaObjectCalculationPluginConfig):
1214 pass
1217@register("ap_stetsonJ")
1218class StetsonJDiaPsfFlux(DiaObjectCalculationPlugin):
1219 """Compute the StetsonJ statistic on the DIA point source fluxes.
1220 """
1222 ConfigClass = LinearFitDiaPsfFluxConfig
1224 # Required input Cols
1225 inputCols = ["psfFluxMean"]
1226 # Output columns are created upon instantiation of the class.
1227 outputCols = ["psfFluxStetsonJ"]
1228 plugType = "multi"
1229 needsFilter = True
1231 @classmethod
1232 def getExecutionOrder(cls):
1233 return cls.FLUX_MOMENTS_CALCULATED
1235 def calculate(self,
1236 diaObjects,
1237 diaSources,
1238 filterDiaSources,
1239 band,
1240 **kwargs):
1241 """Compute the StetsonJ statistic on the DIA point source fluxes.
1243 Parameters
1244 ----------
1245 diaObject : `dict`
1246 Summary object to store values in.
1247 diaSources : `pandas.DataFrame`
1248 DataFrame representing all diaSources associated with this
1249 diaObject.
1250 filterDiaSources : `pandas.DataFrame`
1251 DataFrame representing diaSources associated with this
1252 diaObject that are observed in the band pass ``band``.
1253 band : `str`
1254 Simple, string name of the filter for the flux being calculated.
1255 **kwargs
1256 Any additional keyword arguments that may be passed to the plugin.
1257 """
1258 meanName = "{}_psfFluxMean".format(band)
1260 def _stetsonJ(df):
1261 tmpDf = df[~np.logical_or(np.isnan(df["psfFlux"]),
1262 np.isnan(df["psfFluxErr"]))]
1263 if len(tmpDf) < 2:
1264 return np.nan
1265 fluxes = tmpDf["psfFlux"].to_numpy()
1266 errors = tmpDf["psfFluxErr"].to_numpy()
1268 return self._stetson_J(
1269 fluxes,
1270 errors,
1271 diaObjects.at[tmpDf.diaObjectId.iat[0], meanName])
1273 column = "{}_psfFluxStetsonJ".format(band)
1274 if column in diaObjects:
1275 dtype = diaObjects[column].dtype
1276 diaObjects.loc[:, column] = filterDiaSources.apply(_stetsonJ).astype(dtype)
1277 else:
1278 diaObjects.loc[:, column] = filterDiaSources.apply(_stetsonJ)
1280 def _stetson_J(self, fluxes, errors, mean=None):
1281 """Compute the single band stetsonJ statistic.
1283 Parameters
1284 ----------
1285 fluxes : `numpy.ndarray` (N,)
1286 Calibrated lightcurve flux values.
1287 errors : `numpy.ndarray` (N,)
1288 Errors on the calibrated lightcurve fluxes.
1289 mean : `float`
1290 Starting mean from previous plugin.
1292 Returns
1293 -------
1294 stetsonJ : `float`
1295 stetsonJ statistic for the input fluxes and errors.
1297 References
1298 ----------
1299 .. [1] Stetson, P. B., "On the Automatic Determination of Light-Curve
1300 Parameters for Cepheid Variables", PASP, 108, 851S, 1996
1301 """
1302 n_points = len(fluxes)
1303 flux_mean = self._stetson_mean(fluxes, errors, mean)
1304 delta_val = (
1305 np.sqrt(n_points / (n_points - 1)) * (fluxes - flux_mean) / errors)
1306 p_k = delta_val ** 2 - 1
1308 return np.mean(np.sign(p_k) * np.sqrt(np.fabs(p_k)))
1310 def _stetson_mean(self,
1311 values,
1312 errors,
1313 mean=None,
1314 alpha=2.,
1315 beta=2.,
1316 n_iter=20,
1317 tol=1e-6):
1318 """Compute the stetson mean of the fluxes which down-weights outliers.
1320 Weighted biased on an error weighted difference scaled by a constant
1321 (1/``a``) and raised to the power beta. Higher betas more harshly
1322 penalize outliers and ``a`` sets the number of sigma where a weighted
1323 difference of 1 occurs.
1325 Parameters
1326 ----------
1327 values : `numpy.dnarray`, (N,)
1328 Input values to compute the mean of.
1329 errors : `numpy.ndarray`, (N,)
1330 Errors on the input values.
1331 mean : `float`
1332 Starting mean value or None.
1333 alpha : `float`
1334 Scalar down-weighting of the fractional difference. lower->more
1335 clipping. (Default value is 2.)
1336 beta : `float`
1337 Power law slope of the used to down-weight outliers. higher->more
1338 clipping. (Default value is 2.)
1339 n_iter : `int`
1340 Number of iterations of clipping.
1341 tol : `float`
1342 Fractional and absolute tolerance goal on the change in the mean
1343 before exiting early. (Default value is 1e-6)
1345 Returns
1346 -------
1347 mean : `float`
1348 Weighted stetson mean result.
1350 References
1351 ----------
1352 .. [1] Stetson, P. B., "On the Automatic Determination of Light-Curve
1353 Parameters for Cepheid Variables", PASP, 108, 851S, 1996
1354 """
1355 n_points = len(values)
1356 n_factor = np.sqrt(n_points / (n_points - 1))
1357 inv_var = 1 / errors ** 2
1359 if mean is None:
1360 mean = np.average(values, weights=inv_var)
1361 for iter_idx in range(n_iter):
1362 chi = np.fabs(n_factor * (values - mean) / errors)
1363 tmp_mean = np.average(
1364 values,
1365 weights=inv_var / (1 + (chi / alpha) ** beta))
1366 diff = np.fabs(tmp_mean - mean)
1367 mean = tmp_mean
1368 if diff / mean < tol and diff < tol:
1369 break
1370 return mean
1373class WeightedMeanDiaTotFluxConfig(DiaObjectCalculationPluginConfig):
1374 pass
1377@register("ap_meanTotFlux")
1378class WeightedMeanDiaTotFlux(DiaObjectCalculationPlugin):
1379 """Compute the weighted mean and mean error on the point source fluxes
1380 forced photometered at the DiaSource location in the calibrated image.
1382 Additionally store number of usable data points.
1383 """
1385 ConfigClass = WeightedMeanDiaPsfFluxConfig
1386 outputCols = ["scienceFluxMean", "scienceFluxMeanErr"]
1387 plugType = "multi"
1388 needsFilter = True
1390 @classmethod
1391 def getExecutionOrder(cls):
1392 return cls.DEFAULT_CATALOGCALCULATION
1394 @catchWarnings(warns=["invalid value encountered",
1395 "divide by zero"])
1396 def calculate(self,
1397 diaObjects,
1398 diaSources,
1399 filterDiaSources,
1400 band,
1401 **kwargs):
1402 """Compute the weighted mean and mean error of the point source flux.
1404 Parameters
1405 ----------
1406 diaObject : `dict`
1407 Summary object to store values in.
1408 diaSources : `pandas.DataFrame`
1409 DataFrame representing all diaSources associated with this
1410 diaObject.
1411 filterDiaSources : `pandas.DataFrame`
1412 DataFrame representing diaSources associated with this
1413 diaObject that are observed in the band pass ``band``.
1414 band : `str`
1415 Simple, string name of the filter for the flux being calculated.
1416 **kwargs
1417 Any additional keyword arguments that may be passed to the plugin.
1418 """
1419 totMeanName = "{}_scienceFluxMean".format(band)
1420 if totMeanName not in diaObjects.columns:
1421 diaObjects[totMeanName] = np.nan
1422 totErrName = "{}_scienceFluxMeanErr".format(band)
1423 if totErrName not in diaObjects.columns:
1424 diaObjects[totErrName] = np.nan
1426 def _meanFlux(df):
1427 tmpDf = df[~np.logical_or(np.isnan(df["scienceFlux"]),
1428 np.isnan(df["scienceFluxErr"]))]
1429 tot_weight = np.nansum(1 / tmpDf["scienceFluxErr"] ** 2)
1430 fluxMean = np.nansum(tmpDf["scienceFlux"]
1431 / tmpDf["scienceFluxErr"] ** 2)
1432 fluxMean /= tot_weight
1433 fluxMeanErr = np.sqrt(1 / tot_weight)
1435 return pd.Series({totMeanName: fluxMean,
1436 totErrName: fluxMeanErr})
1438 df = filterDiaSources.apply(_meanFlux).astype(diaObjects.dtypes[[totMeanName, totErrName]])
1439 diaObjects.loc[:, [totMeanName, totErrName]] = df
1442class SigmaDiaTotFluxConfig(DiaObjectCalculationPluginConfig):
1443 pass
1446@register("ap_sigmaTotFlux")
1447class SigmaDiaTotFlux(DiaObjectCalculationPlugin):
1448 """Compute scatter of diaSource fluxes.
1449 """
1451 ConfigClass = SigmaDiaPsfFluxConfig
1452 # Output columns are created upon instantiation of the class.
1453 outputCols = ["scienceFluxSigma"]
1454 plugType = "multi"
1455 needsFilter = True
1457 @classmethod
1458 def getExecutionOrder(cls):
1459 return cls.DEFAULT_CATALOGCALCULATION
1461 def calculate(self,
1462 diaObjects,
1463 diaSources,
1464 filterDiaSources,
1465 band,
1466 **kwargs):
1467 """Compute the sigma fluxes of the point source flux measured on the
1468 calibrated image.
1470 Parameters
1471 ----------
1472 diaObject : `dict`
1473 Summary object to store values in.
1474 diaSources : `pandas.DataFrame`
1475 DataFrame representing all diaSources associated with this
1476 diaObject.
1477 filterDiaSources : `pandas.DataFrame`
1478 DataFrame representing diaSources associated with this
1479 diaObject that are observed in the band pass ``band``.
1480 band : `str`
1481 Simple, string name of the filter for the flux being calculated.
1482 **kwargs
1483 Any additional keyword arguments that may be passed to the plugin.
1484 """
1485 # Set "delta degrees of freedom (ddf)" to 1 to calculate the unbiased
1486 # estimator of scatter (i.e. 'N - 1' instead of 'N').
1487 column = "{}_scienceFluxSigma".format(band)
1488 if column in diaObjects:
1489 dtype = diaObjects[column].dtype
1490 diaObjects.loc[:, column] = filterDiaSources.scienceFlux.std().astype(dtype)
1491 else:
1493 diaObjects.loc[:, column] = filterDiaSources.scienceFlux.std()