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"]
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()))
461 return pd.Series({
"ra": aveCoord.getRa().asDegrees(),
462 "dec": aveCoord.getDec().asDegrees()})
464 ans = diaSources.apply(_computeMeanPos)
465 diaObjects.loc[:, [
"ra",
"dec"]] = ans
470 htmLevel = pexConfig.Field(
472 doc=
"Level of the HTM pixelization.",
477@register("ap_HTMIndex")
479 """Compute the mean position of a DiaObject given a set of DiaSources.
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.
488 ConfigClass = HTMIndexDiaPositionConfig
492 inputCols = [
"ra",
"dec"]
493 outputCols = [
"pixelId"]
497 DiaObjectCalculationPlugin.__init__(self, config, name, metadata)
505 """Compute the mean position of a DiaObject given a set of DiaSources
509 diaObjects : `pandas.dataFrame`
510 Summary objects to store values in and read ra/dec from.
512 Id of the diaObject to update.
514 Any additional keyword arguments that may be passed to the plugin.
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())
529 """Compute the total number of DiaSources associated with this DiaObject.
532 ConfigClass = NumDiaSourcesDiaPluginConfig
533 outputCols = [
"nDiaSources"]
542 """Compute the total number of DiaSources associated with this DiaObject.
547 Summary object to store values in and read ra/dec from.
549 Any additional keyword arguments that may be passed to the plugin.
551 dtype = diaObjects[
"nDiaSources"].dtype
552 diaObjects.loc[:,
"nDiaSources"] = diaSources.diaObjectId.count().astype(dtype)
561 """Find if any DiaSource is flagged.
563 Set the DiaObject flag if any DiaSource is flagged.
566 ConfigClass = NumDiaSourcesDiaPluginConfig
567 outputCols = [
"flags"]
576 """Find if any DiaSource is flagged.
578 Set the DiaObject flag if any DiaSource is flagged.
583 Summary object to store values in and read ra/dec from.
585 Any additional keyword arguments that may be passed to the plugin.
587 dtype = diaObjects[
"flags"].dtype
588 diaObjects.loc[:,
"flags"] = diaSources.flags.any().astype(dtype)
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.
603 ConfigClass = WeightedMeanDiaPsfFluxConfig
604 outputCols = [
"psfFluxMean",
"psfFluxMeanErr",
"psfFluxNdata"]
612 @catchWarnings(warns=[
"invalid value encountered",
620 """Compute the weighted mean and mean error of the point source flux.
625 Summary object to store values in.
626 diaSources : `pandas.DataFrame`
627 DataFrame representing all diaSources associated with this
629 filterDiaSources : `pandas.DataFrame`
630 DataFrame representing diaSources associated with this
631 diaObject that are observed in the band pass ``band``.
633 Simple, string name of the filter for the flux being calculated.
635 Any additional keyword arguments that may be passed to the plugin.
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
655 fluxMeanErr = np.sqrt(1 / tot_weight)
658 nFluxData = len(tmpDf)
660 return pd.Series({meanName: fluxMean,
661 errName: fluxMeanErr,
662 nDataName: nFluxData},
664 df = filterDiaSources.apply(_weightedMean).astype(diaObjects.dtypes[[meanName, errName, nDataName]])
666 diaObjects.loc[:, [meanName, errName, nDataName]] = df
670 percentiles = pexConfig.ListField(
672 default=[5, 25, 50, 75, 95],
673 doc=
"Percentiles to calculate to compute values for. Should be "
678@register("ap_percentileFlux")
680 """Compute percentiles of diaSource fluxes.
683 ConfigClass = PercentileDiaPsfFluxConfig
689 def __init__(self, config, name, metadata, **kwargs):
690 DiaObjectCalculationPlugin.__init__(self,
695 self.
outputCols = [
"psfFluxPercentile{:02d}".format(percent)
696 for percent
in self.
config.percentiles]
702 @catchWarnings(warns=["All-NaN slice encountered"])
709 """Compute the percentile fluxes of the point source flux.
714 Summary object to store values in.
715 diaSources : `pandas.DataFrame`
716 DataFrame representing all diaSources associated with this
718 filterDiaSources : `pandas.DataFrame`
719 DataFrame representing diaSources associated with this
720 diaObject that are observed in the band pass ``band``.
722 Simple, string name of the filter for the flux being calculated.
724 Any additional keyword arguments that may be passed to the plugin.
728 for tilePercent
in self.
config.percentiles:
729 pTileName =
"{}_psfFluxPercentile{:02d}".format(band,
731 pTileNames.append(pTileName)
732 if pTileName
not in diaObjects.columns:
733 diaObjects[pTileName] = np.nan
735 dtype = diaObjects[pTileName].dtype
737 def _fluxPercentiles(df):
738 pTiles = np.nanpercentile(df[
"psfFlux"], self.
config.percentiles)
740 dict((tileName, pTile)
741 for tileName, pTile
in zip(pTileNames, pTiles)))
743 diaObjects.loc[:, pTileNames] = filterDiaSources.apply(_fluxPercentiles)
745 diaObjects.loc[:, pTileNames] = filterDiaSources.apply(_fluxPercentiles).astype(dtype)
754 """Compute scatter of diaSource fluxes.
757 ConfigClass = SigmaDiaPsfFluxConfig
759 outputCols = [
"psfFluxSigma"]
773 """Compute the sigma fluxes of the point source flux.
778 Summary object to store values in.
779 diaSources : `pandas.DataFrame`
780 DataFrame representing all diaSources associated with this
782 filterDiaSources : `pandas.DataFrame`
783 DataFrame representing diaSources associated with this
784 diaObject that are observed in the band pass ``band``.
786 Simple, string name of the filter for the flux being calculated.
788 Any additional keyword arguments that may be passed to the plugin.
792 column =
"{}_psfFluxSigma".format(band)
793 if column
in diaObjects:
794 dtype = diaObjects[column].dtype
795 diaObjects.loc[:, column] = filterDiaSources.psfFlux.std().astype(dtype)
797 diaObjects.loc[:, column] = filterDiaSources.psfFlux.std()
806 """Compute chi2 of diaSource fluxes.
809 ConfigClass = Chi2DiaPsfFluxConfig
812 inputCols = [
"psfFluxMean"]
814 outputCols = [
"psfFluxChi2"]
822 @catchWarnings(warns=["All-NaN slice encountered"])
829 """Compute the chi2 of the point source fluxes.
834 Summary object to store values in.
835 diaSources : `pandas.DataFrame`
836 DataFrame representing all diaSources associated with this
838 filterDiaSources : `pandas.DataFrame`
839 DataFrame representing diaSources associated with this
840 diaObject that are observed in the band pass ``band``.
842 Simple, string name of the filter for the flux being calculated.
844 Any additional keyword arguments that may be passed to the plugin.
846 meanName =
"{}_psfFluxMean".format(band)
847 column =
"{}_psfFluxChi2".format(band)
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)
858 diaObjects.loc[:, column] = filterDiaSources.apply(_chi2)
867 """Compute median absolute deviation of diaSource fluxes.
870 ConfigClass = MadDiaPsfFluxConfig
874 outputCols = [
"psfFluxMAD"]
882 @catchWarnings(warns=["All-NaN slice encountered"])
889 """Compute the median absolute deviation of the point source fluxes.
894 Summary object to store values in.
895 diaSources : `pandas.DataFrame`
896 DataFrame representing all diaSources associated with this
898 filterDiaSources : `pandas.DataFrame`
899 DataFrame representing diaSources associated with this
900 diaObject that are observed in the band pass ``band``.
902 Simple, string name of the filter for the flux being calculated.
904 Any additional keyword arguments that may be passed to the plugin.
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)
913 diaObjects.loc[:, column] = \
914 filterDiaSources.psfFlux.apply(median_absolute_deviation,
924 """Compute the skew of diaSource fluxes.
927 ConfigClass = SkewDiaPsfFluxConfig
931 outputCols = [
"psfFluxSkew"]
945 """Compute the skew of the point source fluxes.
950 Summary object to store values in.
951 diaSources : `pandas.DataFrame`
952 DataFrame representing all diaSources associated with this
954 filterDiaSources : `pandas.DataFrame`
955 DataFrame representing diaSources associated with this
956 diaObject that are observed in the band pass ``band``.
958 Simple, string name of the filter for the flux being calculated.
960 Any additional keyword arguments that may be passed to the plugin.
962 column =
"{}_psfFluxSkew".format(band)
963 if column
in diaObjects:
964 dtype = diaObjects[column].dtype
965 diaObjects.loc[:, column] = filterDiaSources.psfFlux.skew().astype(dtype)
967 diaObjects.loc[:, column] = filterDiaSources.psfFlux.skew()
976 """Compute min/max of diaSource fluxes.
979 ConfigClass = MinMaxDiaPsfFluxConfig
983 outputCols = [
"psfFluxMin",
"psfFluxMax"]
997 """Compute min/max of the point source fluxes.
1002 Summary object to store values in.
1003 diaSources : `pandas.DataFrame`
1004 DataFrame representing all diaSources associated with this
1006 filterDiaSources : `pandas.DataFrame`
1007 DataFrame representing diaSources associated with this
1008 diaObject that are observed in the band pass ``band``.
1010 Simple, string name of the filter for the flux being calculated.
1012 Any additional keyword arguments that may be passed to the plugin.
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)
1032 """Compute the maximum ratio time ordered deltaFlux / deltaTime.
1035 ConfigClass = MinMaxDiaPsfFluxConfig
1039 outputCols = [
"psfFluxMaxSlope"]
1053 """Compute the maximum ratio time ordered deltaFlux / deltaTime.
1058 Summary object to store values in.
1059 diaSources : `pandas.DataFrame`
1060 DataFrame representing all diaSources associated with this
1062 filterDiaSources : `pandas.DataFrame`
1063 DataFrame representing diaSources associated with this
1064 diaObject that are observed in the band pass ``band``.
1066 Simple, string name of the filter for the flux being calculated.
1068 Any additional keyword arguments that may be passed to the plugin.
1072 tmpDf = df[~np.logical_or(np.isnan(df[
"psfFlux"]),
1073 np.isnan(df[
"midpointMjdTai"]))]
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)
1087 diaObjects.loc[:, column] = filterDiaSources.apply(_maxSlope)
1096 """Compute the mean of the dia source errors.
1099 ConfigClass = ErrMeanDiaPsfFluxConfig
1103 outputCols = [
"psfFluxErrMean"]
1117 """Compute the mean of the dia source errors.
1122 Summary object to store values in.
1123 diaSources : `pandas.DataFrame`
1124 DataFrame representing all diaSources associated with this
1126 filterDiaSources : `pandas.DataFrame`
1127 DataFrame representing diaSources associated with this
1128 diaObject that are observed in the band pass ``band``.
1130 Simple, string name of the filter for the flux being calculated.
1132 Any additional keyword arguments that may be passed to the plugin.
1134 column =
"{}_psfFluxErrMean".format(band)
1135 if column
in diaObjects:
1136 dtype = diaObjects[column].dtype
1137 diaObjects.loc[:, column] = filterDiaSources.psfFluxErr.mean().astype(dtype)
1139 diaObjects.loc[:, column] = filterDiaSources.psfFluxErr.mean()
1148 """Compute fit a linear model to flux vs time.
1151 ConfigClass = LinearFitDiaPsfFluxConfig
1155 outputCols = [
"psfFluxLinearSlope",
"psfFluxLinearIntercept"]
1169 """Compute fit a linear model to flux vs time.
1174 Summary object to store values in.
1175 diaSources : `pandas.DataFrame`
1176 DataFrame representing all diaSources associated with this
1178 filterDiaSources : `pandas.DataFrame`
1179 DataFrame representing diaSources associated with this
1180 diaObject that are observed in the band pass ``band``.
1182 Simple, string name of the filter for the flux being calculated.
1184 Any additional keyword arguments that may be passed to the plugin.
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
1196 tmpDf = df[~np.logical_or(
1197 np.isnan(df[
"psfFlux"]),
1198 np.logical_or(np.isnan(df[
"psfFluxErr"]),
1199 np.isnan(df[
"midpointMjdTai"])))]
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)
1218 """Compute the StetsonJ statistic on the DIA point source fluxes.
1221 ConfigClass = LinearFitDiaPsfFluxConfig
1224 inputCols = [
"psfFluxMean"]
1226 outputCols = [
"psfFluxStetsonJ"]
1240 """Compute the StetsonJ statistic on the DIA point source fluxes.
1245 Summary object to store values in.
1246 diaSources : `pandas.DataFrame`
1247 DataFrame representing all diaSources associated with this
1249 filterDiaSources : `pandas.DataFrame`
1250 DataFrame representing diaSources associated with this
1251 diaObject that are observed in the band pass ``band``.
1253 Simple, string name of the filter for the flux being calculated.
1255 Any additional keyword arguments that may be passed to the plugin.
1257 meanName =
"{}_psfFluxMean".format(band)
1260 tmpDf = df[~np.logical_or(np.isnan(df[
"psfFlux"]),
1261 np.isnan(df[
"psfFluxErr"]))]
1264 fluxes = tmpDf[
"psfFlux"].to_numpy()
1265 errors = tmpDf[
"psfFluxErr"].to_numpy()
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)
1277 diaObjects.loc[:, column] = filterDiaSources.apply(_stetsonJ)
1280 """Compute the single band stetsonJ statistic.
1284 fluxes : `numpy.ndarray` (N,)
1285 Calibrated lightcurve flux values.
1286 errors : `numpy.ndarray` (N,)
1287 Errors on the calibrated lightcurve fluxes.
1289 Starting mean from previous plugin.
1294 stetsonJ statistic for the input fluxes and errors.
1298 .. [1] Stetson, P. B., "On the Automatic Determination of Light-Curve
1299 Parameters for Cepheid Variables", PASP, 108, 851S, 1996
1301 n_points = len(fluxes)
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)))
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.
1326 values : `numpy.dnarray`, (N,)
1327 Input values to compute the mean of.
1328 errors : `numpy.ndarray`, (N,)
1329 Errors on the input values.
1331 Starting mean value or None.
1333 Scalar down-weighting of the fractional difference. lower->more
1334 clipping. (Default value is 2.)
1336 Power law slope of the used to down-weight outliers. higher->more
1337 clipping. (Default value is 2.)
1339 Number of iterations of clipping.
1341 Fractional and absolute tolerance goal on the change in the mean
1342 before exiting early. (Default value is 1e-6)
1347 Weighted stetson mean result.
1351 .. [1] Stetson, P. B., "On the Automatic Determination of Light-Curve
1352 Parameters for Cepheid Variables", PASP, 108, 851S, 1996
1354 n_points = len(values)
1355 n_factor = np.sqrt(n_points / (n_points - 1))
1356 inv_var = 1 / errors ** 2
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(
1364 weights=inv_var / (1 + (chi / alpha) ** beta))
1365 diff = np.fabs(tmp_mean - mean)
1367 if diff / mean < tol
and diff < tol:
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.
1384 ConfigClass = WeightedMeanDiaPsfFluxConfig
1385 outputCols = [
"scienceFluxMean",
"scienceFluxMeanErr"]
1393 @catchWarnings(warns=[
"invalid value encountered",
1401 """Compute the weighted mean and mean error of the point source flux.
1406 Summary object to store values in.
1407 diaSources : `pandas.DataFrame`
1408 DataFrame representing all diaSources associated with this
1410 filterDiaSources : `pandas.DataFrame`
1411 DataFrame representing diaSources associated with this
1412 diaObject that are observed in the band pass ``band``.
1414 Simple, string name of the filter for the flux being calculated.
1416 Any additional keyword arguments that may be passed to the plugin.
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
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
1447 """Compute scatter of diaSource fluxes.
1450 ConfigClass = SigmaDiaPsfFluxConfig
1452 outputCols = [
"scienceFluxSigma"]
1466 """Compute the sigma fluxes of the point source flux measured on the
1472 Summary object to store values in.
1473 diaSources : `pandas.DataFrame`
1474 DataFrame representing all diaSources associated with this
1476 filterDiaSources : `pandas.DataFrame`
1477 DataFrame representing diaSources associated with this
1478 diaObject that are observed in the band pass ``band``.
1480 Simple, string name of the filter for the flux being calculated.
1482 Any additional keyword arguments that may be passed to the plugin.
1486 column =
"{}_scienceFluxSigma".format(band)
1487 if column
in diaObjects:
1488 dtype = diaObjects[column].dtype
1489 diaObjects.loc[:, column] = filterDiaSources.scienceFlux.std().astype(dtype)
1492 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=())