22"""Plugins for use in DiaSource summary statistics.
25as defined in the schema of the Apdb both in name and units.
31from astropy.stats
import median_absolute_deviation
34from scipy.optimize
import lsq_linear
39from astropy.timeseries
import LombScargle
40from astropy.timeseries
import LombScargleMultiband
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")
72 """Decorator for generically catching numpy warnings.
74 def decoratorCatchWarnings(func):
75 @functools.wraps(func)
76 def wrapperCatchWarnings(*args, **kwargs):
77 with warnings.catch_warnings():
79 warnings.filterwarnings(
"ignore", val)
80 return func(*args, **kwargs)
81 return wrapperCatchWarnings
84 return decoratorCatchWarnings
86 return decoratorCatchWarnings(_func)
91 Computes an optimized periodogram frequency grid for a given time series.
97 oversampling_factor : `int`, optional
98 The oversampling factor for frequency grid.
99 nyquist_factor : `int`, optional
100 The Nyquist factor for frequency grid.
104 frequencies : `array`
105 The computed optimized periodogram frequency grid.
109 baseline = np.max(x0) - np.min(x0)
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)
126@
register(
"ap_lombScarglePeriodogram")
128 """Compute the single-band period of a DiaObject given a set of DiaSources.
130 ConfigClass = LombScarglePeriodogramConfig
133 outputCols = [
"period",
"power"]
140 @catchWarnings(warns=["All-NaN slice encountered"])
146 """Compute the periodogram.
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.
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.
167 df : `pandas.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.
180 pd_tab : `pandas.Series`
181 The output DataFrame with the Lomb-Scargle parameters.
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)
195 time, oversampling_factor=oversampling_factor, nyquist_factor=nyquist_factor)
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)
210@
register(
"ap_lombScarglePeriodogramMulti")
212 """Compute the multi-band LombScargle periodogram of a DiaObject given a set of DiaSources.
214 ConfigClass = LombScarglePeriodogramMultiConfig
217 outputCols = [
"multiPeriod",
"multiPower",
218 "multiFap",
"multiAmp",
"multiPhase"]
227 """Calculate the False-Alarm probability using the Baluev approximation.
234 The number of detections.
236 The maximum period in the grid.
238 The maximum power in the grid.
242 fap_estimate : `float`
243 The False-Alarm probability Baluev approximation.
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
253 gam_ratio = math.factorial(int((n - 1)/2)) / math.factorial(int((n - 2)/2))
255 return gam_ratio * np.sqrt(
256 4*np.pi*statistics.variance(time)
257 ) * fu * (1-zmax)**((n-4)/2) * np.sqrt(zmax)
261 """Generate the Lomb-Scargle parameters.
264 lsp_model : `astropy.timeseries.LombScargleMultiband`
265 The Lomb-Scargle model.
269 The bands of the time series.
274 The amplitude of the time series.
276 The phase of the time series.
280 .. [1] VanderPlas, J. T., & Ivezic, Z. 2015, ApJ, 812, 18
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]
302 @catchWarnings(warns=["All-NaN slice encountered"])
307 """Compute the multi-band LombScargle periodogram of a DiaObject given
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.
317 Unused kwargs that are always passed to a plugin.
320 bands_arr = diaSources[
'band'].unique().values
321 unique_bands = np.unique(np.concatenate(bands_arr))
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
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.
345 df : `pandas.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.
358 pd_tab : `pandas.Series`
359 The output DataFrame with the Lomb-Scargle parameters.
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,
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
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)
383 time, oversampling_factor=oversampling_factor, nyquist_factor=nyquist_factor)
385 power = lsp.power(f_grid)
388 time, len(time), period[np.argmax(power)], np.max(power))
392 pd_tab = pd.Series({periodCol: period[np.argmax(power)],
393 powerCol: np.max(power),
398 for band
in all_unique_bands:
399 pd_tab[f
"{band}_{ampCol}"] = np.nan
400 pd_tab[f
"{band}_{phaseCol}"] = np.nan
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]
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)
425 """Compute the mean position of a DiaObject given a set of DiaSources.
428 ConfigClass = MeanDiaPositionConfig
432 outputCols = [
"ra",
"dec",
"radecMjdTai"]
440 """Compute the mean ra/dec position of the diaObject given the
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.
450 Any additional keyword arguments that may be passed to the plugin.
453 if outCol
not in diaObjects.columns:
454 diaObjects[outCol] = np.nan
456 def _computeMeanPos(df):
459 for idx, src
in df.iterrows()))
460 ra = aveCoord.getRa().asDegrees()
461 dec = aveCoord.getDec().asDegrees()
462 if np.isnan(ra)
or np.isnan(dec):
465 radecMjdTai = df[
"midpointMjdTai"].max()
467 return pd.Series({
"ra": aveCoord.getRa().asDegrees(),
468 "dec": aveCoord.getDec().asDegrees(),
469 "radecMjdTai": radecMjdTai})
471 ans = diaSources.apply(_computeMeanPos)
472 diaObjects.loc[:, [
"ra",
"dec",
"radecMjdTai"]] = ans
477 htmLevel = pexConfig.Field(
479 doc=
"Level of the HTM pixelization.",
484@register("ap_HTMIndex")
486 """Compute the mean position of a DiaObject given a set of DiaSources.
490 This plugin was implemented to satisfy requirements of old APDB interface
491 which required ``pixelId`` column in DiaObject with HTM20 index. APDB
492 interface had migrated to not need that information, but we keep this
493 plugin in case it may be useful for something else.
495 ConfigClass = HTMIndexDiaPositionConfig
499 inputCols = [
"ra",
"dec"]
500 outputCols = [
"pixelId"]
504 DiaObjectCalculationPlugin.__init__(self, config, name, metadata)
512 """Compute the mean position of a DiaObject given a set of DiaSources
516 diaObjects : `pandas.dataFrame`
517 Summary objects to store values in and read ra/dec from.
519 Id of the diaObject to update.
521 Any additional keyword arguments that may be passed to the plugin.
524 diaObjects.at[diaObjectId,
"ra"] * geom.degrees,
525 diaObjects.at[diaObjectId,
"dec"] * geom.degrees)
526 diaObjects.at[diaObjectId,
"pixelId"] = self.
pixelator.index(
527 sphPoint.getVector())
536 """Compute the total number of DiaSources associated with this DiaObject.
539 ConfigClass = NumDiaSourcesDiaPluginConfig
540 outputCols = [
"nDiaSources"]
549 """Compute the total number of DiaSources associated with this DiaObject.
554 Summary object to store values in and read ra/dec from.
556 Any additional keyword arguments that may be passed to the plugin.
558 dtype = diaObjects[
"nDiaSources"].dtype
559 diaObjects.loc[:,
"nDiaSources"] = diaSources.diaObjectId.count().astype(dtype)
568 """Find if any DiaSource is flagged.
570 Set the DiaObject flag if any DiaSource is flagged.
573 ConfigClass = NumDiaSourcesDiaPluginConfig
574 outputCols = [
"flags"]
583 """Find if any DiaSource is flagged.
585 Set the DiaObject flag if any DiaSource is flagged.
590 Summary object to store values in and read ra/dec from.
592 Any additional keyword arguments that may be passed to the plugin.
594 dtype = diaObjects[
"flags"].dtype
595 diaObjects.loc[:,
"flags"] = diaSources.flags.any().astype(dtype)
604 """Compute the weighted mean and mean error on the point source fluxes
605 of the DiaSource measured on the difference image.
607 Additionally store number of usable data points.
610 ConfigClass = WeightedMeanDiaPsfFluxConfig
611 outputCols = [
"psfFluxMean",
"psfFluxMeanErr",
"psfFluxNdata"]
619 @catchWarnings(warns=[
"invalid value encountered",
627 """Compute the weighted mean and mean error of the point source flux.
632 Summary object to store values in.
633 diaSources : `pandas.DataFrame`
634 DataFrame representing all diaSources associated with this
636 filterDiaSources : `pandas.DataFrame`
637 DataFrame representing diaSources associated with this
638 diaObject that are observed in the band pass ``band``.
640 Simple, string name of the filter for the flux being calculated.
642 Any additional keyword arguments that may be passed to the plugin.
644 meanName =
"{}_psfFluxMean".format(band)
645 errName =
"{}_psfFluxMeanErr".format(band)
646 nDataName =
"{}_psfFluxNdata".format(band)
647 if meanName
not in diaObjects.columns:
648 diaObjects[meanName] = np.nan
649 if errName
not in diaObjects.columns:
650 diaObjects[errName] = np.nan
651 if nDataName
not in diaObjects.columns:
652 diaObjects[nDataName] = 0
654 def _weightedMean(df):
655 tmpDf = df[~np.logical_or(np.isnan(df[
"psfFlux"]),
656 np.isnan(df[
"psfFluxErr"]))]
657 tot_weight = np.nansum(1 / tmpDf[
"psfFluxErr"] ** 2)
658 fluxMean = np.nansum(tmpDf[
"psfFlux"]
659 / tmpDf[
"psfFluxErr"] ** 2)
660 fluxMean /= tot_weight
662 fluxMeanErr = np.sqrt(1 / tot_weight)
665 nFluxData = len(tmpDf)
667 return pd.Series({meanName: fluxMean,
668 errName: fluxMeanErr,
669 nDataName: nFluxData},
671 df = filterDiaSources.apply(_weightedMean).astype(diaObjects.dtypes[[meanName, errName, nDataName]])
673 diaObjects.loc[:, [meanName, errName, nDataName]] = df
677 percentiles = pexConfig.ListField(
679 default=[5, 25, 50, 75, 95],
680 doc=
"Percentiles to calculate to compute values for. Should be "
685@register("ap_percentileFlux")
687 """Compute percentiles of diaSource fluxes.
690 ConfigClass = PercentileDiaPsfFluxConfig
696 def __init__(self, config, name, metadata, **kwargs):
697 DiaObjectCalculationPlugin.__init__(self,
702 self.
outputCols = [
"psfFluxPercentile{:02d}".format(percent)
703 for percent
in self.
config.percentiles]
709 @catchWarnings(warns=["All-NaN slice encountered"])
716 """Compute the percentile fluxes of the point source flux.
721 Summary object to store values in.
722 diaSources : `pandas.DataFrame`
723 DataFrame representing all diaSources associated with this
725 filterDiaSources : `pandas.DataFrame`
726 DataFrame representing diaSources associated with this
727 diaObject that are observed in the band pass ``band``.
729 Simple, string name of the filter for the flux being calculated.
731 Any additional keyword arguments that may be passed to the plugin.
735 for tilePercent
in self.
config.percentiles:
736 pTileName =
"{}_psfFluxPercentile{:02d}".format(band,
738 pTileNames.append(pTileName)
739 if pTileName
not in diaObjects.columns:
740 diaObjects[pTileName] = np.nan
742 dtype = diaObjects[pTileName].dtype
744 def _fluxPercentiles(df):
745 pTiles = np.nanpercentile(df[
"psfFlux"], self.
config.percentiles)
747 dict((tileName, pTile)
748 for tileName, pTile
in zip(pTileNames, pTiles)))
750 diaObjects.loc[:, pTileNames] = filterDiaSources.apply(_fluxPercentiles)
752 diaObjects.loc[:, pTileNames] = filterDiaSources.apply(_fluxPercentiles).astype(dtype)
761 """Compute scatter of diaSource fluxes.
764 ConfigClass = SigmaDiaPsfFluxConfig
766 outputCols = [
"psfFluxSigma"]
780 """Compute the sigma fluxes of the point source flux.
785 Summary object to store values in.
786 diaSources : `pandas.DataFrame`
787 DataFrame representing all diaSources associated with this
789 filterDiaSources : `pandas.DataFrame`
790 DataFrame representing diaSources associated with this
791 diaObject that are observed in the band pass ``band``.
793 Simple, string name of the filter for the flux being calculated.
795 Any additional keyword arguments that may be passed to the plugin.
799 column =
"{}_psfFluxSigma".format(band)
800 if column
in diaObjects:
801 dtype = diaObjects[column].dtype
802 diaObjects.loc[:, column] = filterDiaSources.psfFlux.std().astype(dtype)
804 diaObjects.loc[:, column] = filterDiaSources.psfFlux.std()
813 """Compute chi2 of diaSource fluxes.
816 ConfigClass = Chi2DiaPsfFluxConfig
819 inputCols = [
"psfFluxMean"]
821 outputCols = [
"psfFluxChi2"]
829 @catchWarnings(warns=["All-NaN slice encountered"])
836 """Compute the chi2 of the point source fluxes.
841 Summary object to store values in.
842 diaSources : `pandas.DataFrame`
843 DataFrame representing all diaSources associated with this
845 filterDiaSources : `pandas.DataFrame`
846 DataFrame representing diaSources associated with this
847 diaObject that are observed in the band pass ``band``.
849 Simple, string name of the filter for the flux being calculated.
851 Any additional keyword arguments that may be passed to the plugin.
853 meanName =
"{}_psfFluxMean".format(band)
854 column =
"{}_psfFluxChi2".format(band)
857 delta = (df[
"psfFlux"]
858 - diaObjects.at[df.diaObjectId.iat[0], meanName])
859 return np.nansum((delta / df[
"psfFluxErr"]) ** 2)
861 if column
in diaObjects:
862 dtype = diaObjects[column].dtype
863 diaObjects.loc[:, column] = filterDiaSources.apply(_chi2).astype(dtype)
865 diaObjects.loc[:, column] = filterDiaSources.apply(_chi2)
874 """Compute median absolute deviation of diaSource fluxes.
877 ConfigClass = MadDiaPsfFluxConfig
881 outputCols = [
"psfFluxMAD"]
889 @catchWarnings(warns=["All-NaN slice encountered"])
896 """Compute the median absolute deviation of the point source fluxes.
901 Summary object to store values in.
902 diaSources : `pandas.DataFrame`
903 DataFrame representing all diaSources associated with this
905 filterDiaSources : `pandas.DataFrame`
906 DataFrame representing diaSources associated with this
907 diaObject that are observed in the band pass ``band``.
909 Simple, string name of the filter for the flux being calculated.
911 Any additional keyword arguments that may be passed to the plugin.
913 column =
"{}_psfFluxMAD".format(band)
914 if column
in diaObjects:
915 dtype = diaObjects[column].dtype
916 diaObjects.loc[:, column] = \
917 filterDiaSources.psfFlux.apply(median_absolute_deviation,
918 ignore_nan=
True).astype(dtype)
920 diaObjects.loc[:, column] = \
921 filterDiaSources.psfFlux.apply(median_absolute_deviation,
931 """Compute the skew of diaSource fluxes.
934 ConfigClass = SkewDiaPsfFluxConfig
938 outputCols = [
"psfFluxSkew"]
952 """Compute the skew of the point source fluxes.
957 Summary object to store values in.
958 diaSources : `pandas.DataFrame`
959 DataFrame representing all diaSources associated with this
961 filterDiaSources : `pandas.DataFrame`
962 DataFrame representing diaSources associated with this
963 diaObject that are observed in the band pass ``band``.
965 Simple, string name of the filter for the flux being calculated.
967 Any additional keyword arguments that may be passed to the plugin.
969 column =
"{}_psfFluxSkew".format(band)
970 if column
in diaObjects:
971 dtype = diaObjects[column].dtype
972 diaObjects.loc[:, column] = filterDiaSources.psfFlux.skew().astype(dtype)
974 diaObjects.loc[:, column] = filterDiaSources.psfFlux.skew()
983 """Compute min/max of diaSource fluxes.
986 ConfigClass = MinMaxDiaPsfFluxConfig
990 outputCols = [
"psfFluxMin",
"psfFluxMax"]
1004 """Compute min/max of the point source fluxes.
1009 Summary object to store values in.
1010 diaSources : `pandas.DataFrame`
1011 DataFrame representing all diaSources associated with this
1013 filterDiaSources : `pandas.DataFrame`
1014 DataFrame representing diaSources associated with this
1015 diaObject that are observed in the band pass ``band``.
1017 Simple, string name of the filter for the flux being calculated.
1019 Any additional keyword arguments that may be passed to the plugin.
1021 minName =
"{}_psfFluxMin".format(band)
1022 if minName
not in diaObjects.columns:
1023 diaObjects[minName] = np.nan
1024 maxName =
"{}_psfFluxMax".format(band)
1025 if maxName
not in diaObjects.columns:
1026 diaObjects[maxName] = np.nan
1028 dtype = diaObjects[minName].dtype
1029 diaObjects.loc[:, minName] = filterDiaSources.psfFlux.min().astype(dtype)
1030 diaObjects.loc[:, maxName] = filterDiaSources.psfFlux.max().astype(dtype)
1039 """Compute the maximum ratio time ordered deltaFlux / deltaTime.
1042 ConfigClass = MinMaxDiaPsfFluxConfig
1046 outputCols = [
"psfFluxMaxSlope"]
1060 """Compute the maximum ratio time ordered deltaFlux / deltaTime.
1065 Summary object to store values in.
1066 diaSources : `pandas.DataFrame`
1067 DataFrame representing all diaSources associated with this
1069 filterDiaSources : `pandas.DataFrame`
1070 DataFrame representing diaSources associated with this
1071 diaObject that are observed in the band pass ``band``.
1073 Simple, string name of the filter for the flux being calculated.
1075 Any additional keyword arguments that may be passed to the plugin.
1079 tmpDf = df[~np.logical_or(np.isnan(df[
"psfFlux"]),
1080 np.isnan(df[
"midpointMjdTai"]))]
1083 times = tmpDf[
"midpointMjdTai"].to_numpy()
1084 timeArgs = times.argsort()
1085 times = times[timeArgs]
1086 fluxes = tmpDf[
"psfFlux"].to_numpy()[timeArgs]
1087 return (np.diff(fluxes) / np.diff(times)).max()
1089 column =
"{}_psfFluxMaxSlope".format(band)
1090 if column
in diaObjects:
1091 dtype = diaObjects[column].dtype
1092 diaObjects.loc[:, column] = filterDiaSources.apply(_maxSlope).astype(dtype)
1094 diaObjects.loc[:, column] = filterDiaSources.apply(_maxSlope)
1103 """Compute the mean of the dia source errors.
1106 ConfigClass = ErrMeanDiaPsfFluxConfig
1110 outputCols = [
"psfFluxErrMean"]
1124 """Compute the mean of the dia source errors.
1129 Summary object to store values in.
1130 diaSources : `pandas.DataFrame`
1131 DataFrame representing all diaSources associated with this
1133 filterDiaSources : `pandas.DataFrame`
1134 DataFrame representing diaSources associated with this
1135 diaObject that are observed in the band pass ``band``.
1137 Simple, string name of the filter for the flux being calculated.
1139 Any additional keyword arguments that may be passed to the plugin.
1141 column =
"{}_psfFluxErrMean".format(band)
1142 if column
in diaObjects:
1143 dtype = diaObjects[column].dtype
1144 diaObjects.loc[:, column] = filterDiaSources.psfFluxErr.mean().astype(dtype)
1146 diaObjects.loc[:, column] = filterDiaSources.psfFluxErr.mean()
1155 """Compute fit a linear model to flux vs time.
1158 ConfigClass = LinearFitDiaPsfFluxConfig
1162 outputCols = [
"psfFluxLinearSlope",
"psfFluxLinearIntercept"]
1176 """Compute fit a linear model to flux vs time.
1181 Summary object to store values in.
1182 diaSources : `pandas.DataFrame`
1183 DataFrame representing all diaSources associated with this
1185 filterDiaSources : `pandas.DataFrame`
1186 DataFrame representing diaSources associated with this
1187 diaObject that are observed in the band pass ``band``.
1189 Simple, string name of the filter for the flux being calculated.
1191 Any additional keyword arguments that may be passed to the plugin.
1194 mName =
"{}_psfFluxLinearSlope".format(band)
1195 if mName
not in diaObjects.columns:
1196 diaObjects[mName] = np.nan
1197 bName =
"{}_psfFluxLinearIntercept".format(band)
1198 if bName
not in diaObjects.columns:
1199 diaObjects[bName] = np.nan
1200 dtype = diaObjects[mName].dtype
1203 tmpDf = df[~np.logical_or(
1204 np.isnan(df[
"psfFlux"]),
1205 np.logical_or(np.isnan(df[
"psfFluxErr"]),
1206 np.isnan(df[
"midpointMjdTai"])))]
1208 return pd.Series({mName: np.nan, bName: np.nan})
1209 fluxes = tmpDf[
"psfFlux"].to_numpy()
1210 errors = tmpDf[
"psfFluxErr"].to_numpy()
1211 times = tmpDf[
"midpointMjdTai"].to_numpy()
1212 A = np.array([times / errors, 1 / errors]).transpose()
1213 m, b = lsq_linear(A, fluxes / errors).x
1214 return pd.Series({mName: m, bName: b}, dtype=dtype)
1216 diaObjects.loc[:, [mName, bName]] = filterDiaSources.apply(_linearFit)
1225 """Compute the StetsonJ statistic on the DIA point source fluxes.
1228 ConfigClass = LinearFitDiaPsfFluxConfig
1231 inputCols = [
"psfFluxMean"]
1233 outputCols = [
"psfFluxStetsonJ"]
1247 """Compute the StetsonJ statistic on the DIA point source fluxes.
1252 Summary object to store values in.
1253 diaSources : `pandas.DataFrame`
1254 DataFrame representing all diaSources associated with this
1256 filterDiaSources : `pandas.DataFrame`
1257 DataFrame representing diaSources associated with this
1258 diaObject that are observed in the band pass ``band``.
1260 Simple, string name of the filter for the flux being calculated.
1262 Any additional keyword arguments that may be passed to the plugin.
1264 meanName =
"{}_psfFluxMean".format(band)
1267 tmpDf = df[~np.logical_or(np.isnan(df[
"psfFlux"]),
1268 np.isnan(df[
"psfFluxErr"]))]
1271 fluxes = tmpDf[
"psfFlux"].to_numpy()
1272 errors = tmpDf[
"psfFluxErr"].to_numpy()
1277 diaObjects.at[tmpDf.diaObjectId.iat[0], meanName])
1279 column =
"{}_psfFluxStetsonJ".format(band)
1280 if column
in diaObjects:
1281 dtype = diaObjects[column].dtype
1282 diaObjects.loc[:, column] = filterDiaSources.apply(_stetsonJ).astype(dtype)
1284 diaObjects.loc[:, column] = filterDiaSources.apply(_stetsonJ)
1287 """Compute the single band stetsonJ statistic.
1291 fluxes : `numpy.ndarray` (N,)
1292 Calibrated lightcurve flux values.
1293 errors : `numpy.ndarray` (N,)
1294 Errors on the calibrated lightcurve fluxes.
1296 Starting mean from previous plugin.
1301 stetsonJ statistic for the input fluxes and errors.
1305 .. [1] Stetson, P. B., "On the Automatic Determination of Light-Curve
1306 Parameters for Cepheid Variables", PASP, 108, 851S, 1996
1308 n_points = len(fluxes)
1311 np.sqrt(n_points / (n_points - 1)) * (fluxes - flux_mean) / errors)
1312 p_k = delta_val ** 2 - 1
1314 return np.mean(np.sign(p_k) * np.sqrt(np.fabs(p_k)))
1324 """Compute the stetson mean of the fluxes which down-weights outliers.
1326 Weighted biased on an error weighted difference scaled by a constant
1327 (1/``a``) and raised to the power beta. Higher betas more harshly
1328 penalize outliers and ``a`` sets the number of sigma where a weighted
1329 difference of 1 occurs.
1333 values : `numpy.dnarray`, (N,)
1334 Input values to compute the mean of.
1335 errors : `numpy.ndarray`, (N,)
1336 Errors on the input values.
1338 Starting mean value or None.
1340 Scalar down-weighting of the fractional difference. lower->more
1341 clipping. (Default value is 2.)
1343 Power law slope of the used to down-weight outliers. higher->more
1344 clipping. (Default value is 2.)
1346 Number of iterations of clipping.
1348 Fractional and absolute tolerance goal on the change in the mean
1349 before exiting early. (Default value is 1e-6)
1354 Weighted stetson mean result.
1358 .. [1] Stetson, P. B., "On the Automatic Determination of Light-Curve
1359 Parameters for Cepheid Variables", PASP, 108, 851S, 1996
1361 n_points = len(values)
1362 n_factor = np.sqrt(n_points / (n_points - 1))
1363 inv_var = 1 / errors ** 2
1366 mean = np.average(values, weights=inv_var)
1367 for iter_idx
in range(n_iter):
1368 chi = np.fabs(n_factor * (values - mean) / errors)
1369 tmp_mean = np.average(
1371 weights=inv_var / (1 + (chi / alpha) ** beta))
1372 diff = np.fabs(tmp_mean - mean)
1374 if diff / mean < tol
and diff < tol:
1385 """Compute the weighted mean and mean error on the point source fluxes
1386 forced photometered at the DiaSource location in the calibrated image.
1388 Additionally store number of usable data points.
1391 ConfigClass = WeightedMeanDiaPsfFluxConfig
1392 outputCols = [
"scienceFluxMean",
"scienceFluxMeanErr"]
1400 @catchWarnings(warns=[
"invalid value encountered",
1408 """Compute the weighted mean and mean error of the point source flux.
1413 Summary object to store values in.
1414 diaSources : `pandas.DataFrame`
1415 DataFrame representing all diaSources associated with this
1417 filterDiaSources : `pandas.DataFrame`
1418 DataFrame representing diaSources associated with this
1419 diaObject that are observed in the band pass ``band``.
1421 Simple, string name of the filter for the flux being calculated.
1423 Any additional keyword arguments that may be passed to the plugin.
1425 totMeanName =
"{}_scienceFluxMean".format(band)
1426 if totMeanName
not in diaObjects.columns:
1427 diaObjects[totMeanName] = np.nan
1428 totErrName =
"{}_scienceFluxMeanErr".format(band)
1429 if totErrName
not in diaObjects.columns:
1430 diaObjects[totErrName] = np.nan
1433 tmpDf = df[~np.logical_or(np.isnan(df[
"scienceFlux"]),
1434 np.isnan(df[
"scienceFluxErr"]))]
1435 tot_weight = np.nansum(1 / tmpDf[
"scienceFluxErr"] ** 2)
1436 fluxMean = np.nansum(tmpDf[
"scienceFlux"]
1437 / tmpDf[
"scienceFluxErr"] ** 2)
1438 fluxMean /= tot_weight
1439 fluxMeanErr = np.sqrt(1 / tot_weight)
1441 return pd.Series({totMeanName: fluxMean,
1442 totErrName: fluxMeanErr})
1444 df = filterDiaSources.apply(_meanFlux).astype(diaObjects.dtypes[[totMeanName, totErrName]])
1445 diaObjects.loc[:, [totMeanName, totErrName]] = df
1454 """Compute scatter of diaSource fluxes.
1457 ConfigClass = SigmaDiaPsfFluxConfig
1459 outputCols = [
"scienceFluxSigma"]
1473 """Compute the sigma fluxes of the point source flux measured on the
1479 Summary object to store values in.
1480 diaSources : `pandas.DataFrame`
1481 DataFrame representing all diaSources associated with this
1483 filterDiaSources : `pandas.DataFrame`
1484 DataFrame representing diaSources associated with this
1485 diaObject that are observed in the band pass ``band``.
1487 Simple, string name of the filter for the flux being calculated.
1489 Any additional keyword arguments that may be passed to the plugin.
1493 column =
"{}_scienceFluxSigma".format(band)
1494 if column
in diaObjects:
1495 dtype = diaObjects[column].dtype
1496 diaObjects.loc[:, column] = filterDiaSources.scienceFlux.std().astype(dtype)
1499 diaObjects.loc[:, column] = filterDiaSources.scienceFlux.std()
Point in an unspecified spherical coordinate system.
float FLUX_MOMENTS_CALCULATED
calculate(self, diaObjects, diaSources, filterDiaSources, band, **kwargs)
calculate(self, diaObjects, diaSources, filterDiaSources, band, **kwargs)
calculate(self, diaObjects, diaObjectId, **kwargs)
__init__(self, config, name, metadata)
calculate(self, diaObjects, diaSources, filterDiaSources, band, **kwargs)
calculate(self, diaObjects, diaSources, filterDiaSources, band)
calculate_baluev_fap(time, n, maxPeriod, zmax)
generate_lsp_params(lsp_model, fbest, bands)
calculate(self, diaObjects, diaSources, **kwargs)
calculate(self, diaObjects, diaSources, filterDiaSources, band, **kwargs)
calculate(self, diaObjects, diaSources, filterDiaSources, band, **kwargs)
calculate(self, diaObjects, diaSources, **kwargs)
calculate(self, diaObjects, diaSources, filterDiaSources, band, **kwargs)
calculate(self, diaObjects, diaSources, **kwargs)
__init__(self, config, name, metadata, **kwargs)
calculate(self, diaObjects, diaSources, filterDiaSources, band, **kwargs)
calculate(self, diaObjects, diaSources, filterDiaSources, band, **kwargs)
calculate(self, diaObjects, diaSources, filterDiaSources, band, **kwargs)
calculate(self, diaObjects, diaSources, **kwargs)
calculate(self, diaObjects, diaSources, filterDiaSources, band, **kwargs)
calculate(self, diaObjects, diaSources, filterDiaSources, band, **kwargs)
_stetson_J(self, fluxes, errors, mean=None)
_stetson_mean(self, values, errors, mean=None, alpha=2., beta=2., n_iter=20, tol=1e-6)
calculate(self, diaObjects, diaSources, filterDiaSources, band, **kwargs)
calculate(self, diaObjects, diaSources, filterDiaSources, band, **kwargs)
float DEFAULT_CATALOGCALCULATION
HtmPixelization provides HTM indexing of points and regions.
SpherePoint averageSpherePoint(std::vector< SpherePoint > const &coords)
Return the average of a list of coordinates.
compute_optimized_periodogram_grid(x0, oversampling_factor=5, nyquist_factor=100)
catchWarnings(_func=None, *, warns=[])
register(name, shouldApCorr=False, apCorrList=())