Coverage for python/lsst/meas/base/diaCalculationPlugins.py: 40%
537 statements
« prev ^ index » next coverage.py v7.11.0, created at 2025-11-07 02:09 -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 diaObjects.loc[:, column] = filterDiaSources.psfFlux.std()
800class Chi2DiaPsfFluxConfig(DiaObjectCalculationPluginConfig):
801 pass
804@register("ap_chi2Flux")
805class Chi2DiaPsfFlux(DiaObjectCalculationPlugin):
806 """Compute chi2 of diaSource fluxes.
807 """
809 ConfigClass = Chi2DiaPsfFluxConfig
811 # Required input Cols
812 inputCols = ["psfFluxMean"]
813 # Output columns are created upon instantiation of the class.
814 outputCols = ["psfFluxChi2"]
815 plugType = "multi"
816 needsFilter = True
818 @classmethod
819 def getExecutionOrder(cls):
820 return cls.FLUX_MOMENTS_CALCULATED
822 @catchWarnings(warns=["All-NaN slice encountered"])
823 def calculate(self,
824 diaObjects,
825 diaSources,
826 filterDiaSources,
827 band,
828 **kwargs):
829 """Compute the chi2 of the point source fluxes.
831 Parameters
832 ----------
833 diaObject : `dict`
834 Summary object to store values in.
835 diaSources : `pandas.DataFrame`
836 DataFrame representing all diaSources associated with this
837 diaObject.
838 filterDiaSources : `pandas.DataFrame`
839 DataFrame representing diaSources associated with this
840 diaObject that are observed in the band pass ``band``.
841 band : `str`
842 Simple, string name of the filter for the flux being calculated.
843 **kwargs
844 Any additional keyword arguments that may be passed to the plugin.
845 """
846 meanName = "{}_psfFluxMean".format(band)
847 column = "{}_psfFluxChi2".format(band)
849 def _chi2(df):
850 delta = (df["psfFlux"]
851 - diaObjects.at[df.diaObjectId.iat[0], meanName])
852 return np.nansum((delta / df["psfFluxErr"]) ** 2)
854 if column in diaObjects:
855 dtype = diaObjects[column].dtype
856 diaObjects.loc[:, column] = filterDiaSources.apply(_chi2).astype(dtype)
857 else:
858 diaObjects.loc[:, column] = filterDiaSources.apply(_chi2)
861class MadDiaPsfFluxConfig(DiaObjectCalculationPluginConfig):
862 pass
865@register("ap_madFlux")
866class MadDiaPsfFlux(DiaObjectCalculationPlugin):
867 """Compute median absolute deviation of diaSource fluxes.
868 """
870 ConfigClass = MadDiaPsfFluxConfig
872 # Required input Cols
873 # Output columns are created upon instantiation of the class.
874 outputCols = ["psfFluxMAD"]
875 plugType = "multi"
876 needsFilter = True
878 @classmethod
879 def getExecutionOrder(cls):
880 return cls.DEFAULT_CATALOGCALCULATION
882 @catchWarnings(warns=["All-NaN slice encountered"])
883 def calculate(self,
884 diaObjects,
885 diaSources,
886 filterDiaSources,
887 band,
888 **kwargs):
889 """Compute the median absolute deviation of the point source fluxes.
891 Parameters
892 ----------
893 diaObject : `dict`
894 Summary object to store values in.
895 diaSources : `pandas.DataFrame`
896 DataFrame representing all diaSources associated with this
897 diaObject.
898 filterDiaSources : `pandas.DataFrame`
899 DataFrame representing diaSources associated with this
900 diaObject that are observed in the band pass ``band``.
901 band : `str`
902 Simple, string name of the filter for the flux being calculated.
903 **kwargs
904 Any additional keyword arguments that may be passed to the plugin.
905 """
906 column = "{}_psfFluxMAD".format(band)
907 if column in diaObjects:
908 dtype = diaObjects[column].dtype
909 diaObjects.loc[:, column] = \
910 filterDiaSources.psfFlux.apply(median_absolute_deviation,
911 ignore_nan=True).astype(dtype)
912 else:
913 diaObjects.loc[:, column] = \
914 filterDiaSources.psfFlux.apply(median_absolute_deviation,
915 ignore_nan=True)
918class SkewDiaPsfFluxConfig(DiaObjectCalculationPluginConfig):
919 pass
922@register("ap_skewFlux")
923class SkewDiaPsfFlux(DiaObjectCalculationPlugin):
924 """Compute the skew of diaSource fluxes.
925 """
927 ConfigClass = SkewDiaPsfFluxConfig
929 # Required input Cols
930 # Output columns are created upon instantiation of the class.
931 outputCols = ["psfFluxSkew"]
932 plugType = "multi"
933 needsFilter = True
935 @classmethod
936 def getExecutionOrder(cls):
937 return cls.DEFAULT_CATALOGCALCULATION
939 def calculate(self,
940 diaObjects,
941 diaSources,
942 filterDiaSources,
943 band,
944 **kwargs):
945 """Compute the skew of the point source fluxes.
947 Parameters
948 ----------
949 diaObject : `dict`
950 Summary object to store values in.
951 diaSources : `pandas.DataFrame`
952 DataFrame representing all diaSources associated with this
953 diaObject.
954 filterDiaSources : `pandas.DataFrame`
955 DataFrame representing diaSources associated with this
956 diaObject that are observed in the band pass ``band``.
957 band : `str`
958 Simple, string name of the filter for the flux being calculated.
959 **kwargs
960 Any additional keyword arguments that may be passed to the plugin.
961 """
962 column = "{}_psfFluxSkew".format(band)
963 if column in diaObjects:
964 dtype = diaObjects[column].dtype
965 diaObjects.loc[:, column] = filterDiaSources.psfFlux.skew().astype(dtype)
966 else:
967 diaObjects.loc[:, column] = filterDiaSources.psfFlux.skew()
970class MinMaxDiaPsfFluxConfig(DiaObjectCalculationPluginConfig):
971 pass
974@register("ap_minMaxFlux")
975class MinMaxDiaPsfFlux(DiaObjectCalculationPlugin):
976 """Compute min/max of diaSource fluxes.
977 """
979 ConfigClass = MinMaxDiaPsfFluxConfig
981 # Required input Cols
982 # Output columns are created upon instantiation of the class.
983 outputCols = ["psfFluxMin", "psfFluxMax"]
984 plugType = "multi"
985 needsFilter = True
987 @classmethod
988 def getExecutionOrder(cls):
989 return cls.DEFAULT_CATALOGCALCULATION
991 def calculate(self,
992 diaObjects,
993 diaSources,
994 filterDiaSources,
995 band,
996 **kwargs):
997 """Compute min/max of the point source fluxes.
999 Parameters
1000 ----------
1001 diaObject : `dict`
1002 Summary object to store values in.
1003 diaSources : `pandas.DataFrame`
1004 DataFrame representing all diaSources associated with this
1005 diaObject.
1006 filterDiaSources : `pandas.DataFrame`
1007 DataFrame representing diaSources associated with this
1008 diaObject that are observed in the band pass ``band``.
1009 band : `str`
1010 Simple, string name of the filter for the flux being calculated.
1011 **kwargs
1012 Any additional keyword arguments that may be passed to the plugin.
1013 """
1014 minName = "{}_psfFluxMin".format(band)
1015 if minName not in diaObjects.columns:
1016 diaObjects[minName] = np.nan
1017 maxName = "{}_psfFluxMax".format(band)
1018 if maxName not in diaObjects.columns:
1019 diaObjects[maxName] = np.nan
1021 dtype = diaObjects[minName].dtype
1022 diaObjects.loc[:, minName] = filterDiaSources.psfFlux.min().astype(dtype)
1023 diaObjects.loc[:, maxName] = filterDiaSources.psfFlux.max().astype(dtype)
1026class MaxSlopeDiaPsfFluxConfig(DiaObjectCalculationPluginConfig):
1027 pass
1030@register("ap_maxSlopeFlux")
1031class MaxSlopeDiaPsfFlux(DiaObjectCalculationPlugin):
1032 """Compute the maximum ratio time ordered deltaFlux / deltaTime.
1033 """
1035 ConfigClass = MinMaxDiaPsfFluxConfig
1037 # Required input Cols
1038 # Output columns are created upon instantiation of the class.
1039 outputCols = ["psfFluxMaxSlope"]
1040 plugType = "multi"
1041 needsFilter = True
1043 @classmethod
1044 def getExecutionOrder(cls):
1045 return cls.DEFAULT_CATALOGCALCULATION
1047 def calculate(self,
1048 diaObjects,
1049 diaSources,
1050 filterDiaSources,
1051 band,
1052 **kwargs):
1053 """Compute the maximum ratio time ordered deltaFlux / deltaTime.
1055 Parameters
1056 ----------
1057 diaObject : `dict`
1058 Summary object to store values in.
1059 diaSources : `pandas.DataFrame`
1060 DataFrame representing all diaSources associated with this
1061 diaObject.
1062 filterDiaSources : `pandas.DataFrame`
1063 DataFrame representing diaSources associated with this
1064 diaObject that are observed in the band pass ``band``.
1065 band : `str`
1066 Simple, string name of the filter for the flux being calculated.
1067 **kwargs
1068 Any additional keyword arguments that may be passed to the plugin.
1069 """
1071 def _maxSlope(df):
1072 tmpDf = df[~np.logical_or(np.isnan(df["psfFlux"]),
1073 np.isnan(df["midpointMjdTai"]))]
1074 if len(tmpDf) < 2:
1075 return np.nan
1076 times = tmpDf["midpointMjdTai"].to_numpy()
1077 timeArgs = times.argsort()
1078 times = times[timeArgs]
1079 fluxes = tmpDf["psfFlux"].to_numpy()[timeArgs]
1080 return (np.diff(fluxes) / np.diff(times)).max()
1082 column = "{}_psfFluxMaxSlope".format(band)
1083 if column in diaObjects:
1084 dtype = diaObjects[column].dtype
1085 diaObjects.loc[:, column] = filterDiaSources.apply(_maxSlope).astype(dtype)
1086 else:
1087 diaObjects.loc[:, column] = filterDiaSources.apply(_maxSlope)
1090class ErrMeanDiaPsfFluxConfig(DiaObjectCalculationPluginConfig):
1091 pass
1094@register("ap_meanErrFlux")
1095class ErrMeanDiaPsfFlux(DiaObjectCalculationPlugin):
1096 """Compute the mean of the dia source errors.
1097 """
1099 ConfigClass = ErrMeanDiaPsfFluxConfig
1101 # Required input Cols
1102 # Output columns are created upon instantiation of the class.
1103 outputCols = ["psfFluxErrMean"]
1104 plugType = "multi"
1105 needsFilter = True
1107 @classmethod
1108 def getExecutionOrder(cls):
1109 return cls.DEFAULT_CATALOGCALCULATION
1111 def calculate(self,
1112 diaObjects,
1113 diaSources,
1114 filterDiaSources,
1115 band,
1116 **kwargs):
1117 """Compute the mean of the dia source errors.
1119 Parameters
1120 ----------
1121 diaObject : `dict`
1122 Summary object to store values in.
1123 diaSources : `pandas.DataFrame`
1124 DataFrame representing all diaSources associated with this
1125 diaObject.
1126 filterDiaSources : `pandas.DataFrame`
1127 DataFrame representing diaSources associated with this
1128 diaObject that are observed in the band pass ``band``.
1129 band : `str`
1130 Simple, string name of the filter for the flux being calculated.
1131 **kwargs
1132 Any additional keyword arguments that may be passed to the plugin.
1133 """
1134 column = "{}_psfFluxErrMean".format(band)
1135 if column in diaObjects:
1136 dtype = diaObjects[column].dtype
1137 diaObjects.loc[:, column] = filterDiaSources.psfFluxErr.mean().astype(dtype)
1138 else:
1139 diaObjects.loc[:, column] = filterDiaSources.psfFluxErr.mean()
1142class LinearFitDiaPsfFluxConfig(DiaObjectCalculationPluginConfig):
1143 pass
1146@register("ap_linearFit")
1147class LinearFitDiaPsfFlux(DiaObjectCalculationPlugin):
1148 """Compute fit a linear model to flux vs time.
1149 """
1151 ConfigClass = LinearFitDiaPsfFluxConfig
1153 # Required input Cols
1154 # Output columns are created upon instantiation of the class.
1155 outputCols = ["psfFluxLinearSlope", "psfFluxLinearIntercept"]
1156 plugType = "multi"
1157 needsFilter = True
1159 @classmethod
1160 def getExecutionOrder(cls):
1161 return cls.DEFAULT_CATALOGCALCULATION
1163 def calculate(self,
1164 diaObjects,
1165 diaSources,
1166 filterDiaSources,
1167 band,
1168 **kwargs):
1169 """Compute fit a linear model to flux vs time.
1171 Parameters
1172 ----------
1173 diaObject : `dict`
1174 Summary object to store values in.
1175 diaSources : `pandas.DataFrame`
1176 DataFrame representing all diaSources associated with this
1177 diaObject.
1178 filterDiaSources : `pandas.DataFrame`
1179 DataFrame representing diaSources associated with this
1180 diaObject that are observed in the band pass ``band``.
1181 band : `str`
1182 Simple, string name of the filter for the flux being calculated.
1183 **kwargs
1184 Any additional keyword arguments that may be passed to the plugin.
1185 """
1187 mName = "{}_psfFluxLinearSlope".format(band)
1188 if mName not in diaObjects.columns:
1189 diaObjects[mName] = np.nan
1190 bName = "{}_psfFluxLinearIntercept".format(band)
1191 if bName not in diaObjects.columns:
1192 diaObjects[bName] = np.nan
1193 dtype = diaObjects[mName].dtype
1195 def _linearFit(df):
1196 tmpDf = df[~np.logical_or(
1197 np.isnan(df["psfFlux"]),
1198 np.logical_or(np.isnan(df["psfFluxErr"]),
1199 np.isnan(df["midpointMjdTai"])))]
1200 if len(tmpDf) < 2:
1201 return pd.Series({mName: np.nan, bName: np.nan})
1202 fluxes = tmpDf["psfFlux"].to_numpy()
1203 errors = tmpDf["psfFluxErr"].to_numpy()
1204 times = tmpDf["midpointMjdTai"].to_numpy()
1205 A = np.array([times / errors, 1 / errors]).transpose()
1206 m, b = lsq_linear(A, fluxes / errors).x
1207 return pd.Series({mName: m, bName: b}, dtype=dtype)
1209 diaObjects.loc[:, [mName, bName]] = filterDiaSources.apply(_linearFit)
1212class StetsonJDiaPsfFluxConfig(DiaObjectCalculationPluginConfig):
1213 pass
1216@register("ap_stetsonJ")
1217class StetsonJDiaPsfFlux(DiaObjectCalculationPlugin):
1218 """Compute the StetsonJ statistic on the DIA point source fluxes.
1219 """
1221 ConfigClass = LinearFitDiaPsfFluxConfig
1223 # Required input Cols
1224 inputCols = ["psfFluxMean"]
1225 # Output columns are created upon instantiation of the class.
1226 outputCols = ["psfFluxStetsonJ"]
1227 plugType = "multi"
1228 needsFilter = True
1230 @classmethod
1231 def getExecutionOrder(cls):
1232 return cls.FLUX_MOMENTS_CALCULATED
1234 def calculate(self,
1235 diaObjects,
1236 diaSources,
1237 filterDiaSources,
1238 band,
1239 **kwargs):
1240 """Compute the StetsonJ statistic on the DIA point source fluxes.
1242 Parameters
1243 ----------
1244 diaObject : `dict`
1245 Summary object to store values in.
1246 diaSources : `pandas.DataFrame`
1247 DataFrame representing all diaSources associated with this
1248 diaObject.
1249 filterDiaSources : `pandas.DataFrame`
1250 DataFrame representing diaSources associated with this
1251 diaObject that are observed in the band pass ``band``.
1252 band : `str`
1253 Simple, string name of the filter for the flux being calculated.
1254 **kwargs
1255 Any additional keyword arguments that may be passed to the plugin.
1256 """
1257 meanName = "{}_psfFluxMean".format(band)
1259 def _stetsonJ(df):
1260 tmpDf = df[~np.logical_or(np.isnan(df["psfFlux"]),
1261 np.isnan(df["psfFluxErr"]))]
1262 if len(tmpDf) < 2:
1263 return np.nan
1264 fluxes = tmpDf["psfFlux"].to_numpy()
1265 errors = tmpDf["psfFluxErr"].to_numpy()
1267 return self._stetson_J(
1268 fluxes,
1269 errors,
1270 diaObjects.at[tmpDf.diaObjectId.iat[0], meanName])
1272 column = "{}_psfFluxStetsonJ".format(band)
1273 if column in diaObjects:
1274 dtype = diaObjects[column].dtype
1275 diaObjects.loc[:, column] = filterDiaSources.apply(_stetsonJ).astype(dtype)
1276 else:
1277 diaObjects.loc[:, column] = filterDiaSources.apply(_stetsonJ)
1279 def _stetson_J(self, fluxes, errors, mean=None):
1280 """Compute the single band stetsonJ statistic.
1282 Parameters
1283 ----------
1284 fluxes : `numpy.ndarray` (N,)
1285 Calibrated lightcurve flux values.
1286 errors : `numpy.ndarray` (N,)
1287 Errors on the calibrated lightcurve fluxes.
1288 mean : `float`
1289 Starting mean from previous plugin.
1291 Returns
1292 -------
1293 stetsonJ : `float`
1294 stetsonJ statistic for the input fluxes and errors.
1296 References
1297 ----------
1298 .. [1] Stetson, P. B., "On the Automatic Determination of Light-Curve
1299 Parameters for Cepheid Variables", PASP, 108, 851S, 1996
1300 """
1301 n_points = len(fluxes)
1302 flux_mean = self._stetson_mean(fluxes, errors, mean)
1303 delta_val = (
1304 np.sqrt(n_points / (n_points - 1)) * (fluxes - flux_mean) / errors)
1305 p_k = delta_val ** 2 - 1
1307 return np.mean(np.sign(p_k) * np.sqrt(np.fabs(p_k)))
1309 def _stetson_mean(self,
1310 values,
1311 errors,
1312 mean=None,
1313 alpha=2.,
1314 beta=2.,
1315 n_iter=20,
1316 tol=1e-6):
1317 """Compute the stetson mean of the fluxes which down-weights outliers.
1319 Weighted biased on an error weighted difference scaled by a constant
1320 (1/``a``) and raised to the power beta. Higher betas more harshly
1321 penalize outliers and ``a`` sets the number of sigma where a weighted
1322 difference of 1 occurs.
1324 Parameters
1325 ----------
1326 values : `numpy.dnarray`, (N,)
1327 Input values to compute the mean of.
1328 errors : `numpy.ndarray`, (N,)
1329 Errors on the input values.
1330 mean : `float`
1331 Starting mean value or None.
1332 alpha : `float`
1333 Scalar down-weighting of the fractional difference. lower->more
1334 clipping. (Default value is 2.)
1335 beta : `float`
1336 Power law slope of the used to down-weight outliers. higher->more
1337 clipping. (Default value is 2.)
1338 n_iter : `int`
1339 Number of iterations of clipping.
1340 tol : `float`
1341 Fractional and absolute tolerance goal on the change in the mean
1342 before exiting early. (Default value is 1e-6)
1344 Returns
1345 -------
1346 mean : `float`
1347 Weighted stetson mean result.
1349 References
1350 ----------
1351 .. [1] Stetson, P. B., "On the Automatic Determination of Light-Curve
1352 Parameters for Cepheid Variables", PASP, 108, 851S, 1996
1353 """
1354 n_points = len(values)
1355 n_factor = np.sqrt(n_points / (n_points - 1))
1356 inv_var = 1 / errors ** 2
1358 if mean is None:
1359 mean = np.average(values, weights=inv_var)
1360 for iter_idx in range(n_iter):
1361 chi = np.fabs(n_factor * (values - mean) / errors)
1362 tmp_mean = np.average(
1363 values,
1364 weights=inv_var / (1 + (chi / alpha) ** beta))
1365 diff = np.fabs(tmp_mean - mean)
1366 mean = tmp_mean
1367 if diff / mean < tol and diff < tol:
1368 break
1369 return mean
1372class WeightedMeanDiaTotFluxConfig(DiaObjectCalculationPluginConfig):
1373 pass
1376@register("ap_meanTotFlux")
1377class WeightedMeanDiaTotFlux(DiaObjectCalculationPlugin):
1378 """Compute the weighted mean and mean error on the point source fluxes
1379 forced photometered at the DiaSource location in the calibrated image.
1381 Additionally store number of usable data points.
1382 """
1384 ConfigClass = WeightedMeanDiaPsfFluxConfig
1385 outputCols = ["scienceFluxMean", "scienceFluxMeanErr"]
1386 plugType = "multi"
1387 needsFilter = True
1389 @classmethod
1390 def getExecutionOrder(cls):
1391 return cls.DEFAULT_CATALOGCALCULATION
1393 @catchWarnings(warns=["invalid value encountered",
1394 "divide by zero"])
1395 def calculate(self,
1396 diaObjects,
1397 diaSources,
1398 filterDiaSources,
1399 band,
1400 **kwargs):
1401 """Compute the weighted mean and mean error of the point source flux.
1403 Parameters
1404 ----------
1405 diaObject : `dict`
1406 Summary object to store values in.
1407 diaSources : `pandas.DataFrame`
1408 DataFrame representing all diaSources associated with this
1409 diaObject.
1410 filterDiaSources : `pandas.DataFrame`
1411 DataFrame representing diaSources associated with this
1412 diaObject that are observed in the band pass ``band``.
1413 band : `str`
1414 Simple, string name of the filter for the flux being calculated.
1415 **kwargs
1416 Any additional keyword arguments that may be passed to the plugin.
1417 """
1418 totMeanName = "{}_scienceFluxMean".format(band)
1419 if totMeanName not in diaObjects.columns:
1420 diaObjects[totMeanName] = np.nan
1421 totErrName = "{}_scienceFluxMeanErr".format(band)
1422 if totErrName not in diaObjects.columns:
1423 diaObjects[totErrName] = np.nan
1425 def _meanFlux(df):
1426 tmpDf = df[~np.logical_or(np.isnan(df["scienceFlux"]),
1427 np.isnan(df["scienceFluxErr"]))]
1428 tot_weight = np.nansum(1 / tmpDf["scienceFluxErr"] ** 2)
1429 fluxMean = np.nansum(tmpDf["scienceFlux"]
1430 / tmpDf["scienceFluxErr"] ** 2)
1431 fluxMean /= tot_weight
1432 fluxMeanErr = np.sqrt(1 / tot_weight)
1434 return pd.Series({totMeanName: fluxMean,
1435 totErrName: fluxMeanErr})
1437 df = filterDiaSources.apply(_meanFlux).astype(diaObjects.dtypes[[totMeanName, totErrName]])
1438 diaObjects.loc[:, [totMeanName, totErrName]] = df
1441class SigmaDiaTotFluxConfig(DiaObjectCalculationPluginConfig):
1442 pass
1445@register("ap_sigmaTotFlux")
1446class SigmaDiaTotFlux(DiaObjectCalculationPlugin):
1447 """Compute scatter of diaSource fluxes.
1448 """
1450 ConfigClass = SigmaDiaPsfFluxConfig
1451 # Output columns are created upon instantiation of the class.
1452 outputCols = ["scienceFluxSigma"]
1453 plugType = "multi"
1454 needsFilter = True
1456 @classmethod
1457 def getExecutionOrder(cls):
1458 return cls.DEFAULT_CATALOGCALCULATION
1460 def calculate(self,
1461 diaObjects,
1462 diaSources,
1463 filterDiaSources,
1464 band,
1465 **kwargs):
1466 """Compute the sigma fluxes of the point source flux measured on the
1467 calibrated image.
1469 Parameters
1470 ----------
1471 diaObject : `dict`
1472 Summary object to store values in.
1473 diaSources : `pandas.DataFrame`
1474 DataFrame representing all diaSources associated with this
1475 diaObject.
1476 filterDiaSources : `pandas.DataFrame`
1477 DataFrame representing diaSources associated with this
1478 diaObject that are observed in the band pass ``band``.
1479 band : `str`
1480 Simple, string name of the filter for the flux being calculated.
1481 **kwargs
1482 Any additional keyword arguments that may be passed to the plugin.
1483 """
1484 # Set "delta degrees of freedom (ddf)" to 1 to calculate the unbiased
1485 # estimator of scatter (i.e. 'N - 1' instead of 'N').
1486 column = "{}_scienceFluxSigma".format(band)
1487 if column in diaObjects:
1488 dtype = diaObjects[column].dtype
1489 diaObjects.loc[:, column] = filterDiaSources.scienceFlux.std().astype(dtype)
1490 else:
1492 diaObjects.loc[:, column] = filterDiaSources.scienceFlux.std()