Coverage for python/lsst/meas/base/diaCalculationPlugins.py: 92%
624 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-30 01:57 -0700
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.pipe.base as pipeBase
39import lsst.sphgeom as sphgeom
40from astropy.timeseries import LombScargle
41from astropy.timeseries import LombScargleMultiband
42import math
43import statistics
45from .diaCalculation import (
46 DiaObjectCalculationPluginConfig,
47 DiaObjectCalculationPlugin)
48from .pluginRegistry import register
51__all__ = ("MeanDiaPositionConfig", "MeanDiaPosition",
52 "FirstLastDiaSourceTimeConfig", "FirstLastDiaSourceTime",
53 "HTMIndexDiaPosition", "HTMIndexDiaPositionConfig",
54 "NumDiaSourcesDiaPlugin", "NumDiaSourcesDiaPluginConfig",
55 "SimpleSourceFlagDiaPlugin", "SimpleSourceFlagDiaPluginConfig",
56 "WeightedMeanDiaPsfFluxConfig", "WeightedMeanDiaPsfFlux",
57 "PercentileDiaPsfFlux", "PercentileDiaPsfFluxConfig",
58 "SigmaDiaPsfFlux", "SigmaDiaPsfFluxConfig",
59 "Chi2DiaPsfFlux", "Chi2DiaPsfFluxConfig",
60 "MadDiaPsfFlux", "MadDiaPsfFluxConfig",
61 "SkewDiaPsfFlux", "SkewDiaPsfFluxConfig",
62 "MinMaxDiaPsfFlux", "MinMaxDiaPsfFluxConfig",
63 "MaxSlopeDiaPsfFlux", "MaxSlopeDiaPsfFluxConfig",
64 "ErrMeanDiaPsfFlux", "ErrMeanDiaPsfFluxConfig",
65 "LinearFitDiaPsfFlux", "LinearFitDiaPsfFluxConfig",
66 "StetsonJDiaPsfFlux", "StetsonJDiaPsfFluxConfig",
67 "WeightedMeanDiaTotFlux", "WeightedMeanDiaTotFluxConfig",
68 "SigmaDiaTotFlux", "SigmaDiaTotFluxConfig",
69 "LombScarglePeriodogram", "LombScarglePeriodogramConfig",
70 "LombScarglePeriodogramMulti", "LombScarglePeriodogramMultiConfig",
71 "UnphysicalDiaSourceSeparation")
74def catchWarnings(_func=None, *, warns=[]):
75 """Decorator for generically catching numpy warnings.
76 """
77 def decoratorCatchWarnings(func):
78 @functools.wraps(func)
79 def wrapperCatchWarnings(*args, **kwargs):
80 with warnings.catch_warnings():
81 for val in warns:
82 warnings.filterwarnings("ignore", val)
83 return func(*args, **kwargs)
84 return wrapperCatchWarnings
86 if _func is None: 86 ↛ 89line 86 didn't jump to line 89 because the condition on line 86 was always true
87 return decoratorCatchWarnings
88 else:
89 return decoratorCatchWarnings(_func)
92def typeSafePandasAssignment(
93 target,
94 source,
95 columns,
96 default_dtype=np.float64,
97 int_fill_value=0,
98 # TODO DM-53254: Remove the force_int_to_float hack.
99 force_int_to_float=False,
100):
101 """
102 Assign from a source dataframe to a target dataframe in a type safe way.
104 Parameters
105 ----------
106 target : `pd.DataFrame`
107 Target pandas dataframe.
108 source : `pd.DataFrame` or `pd.Series`
109 Grouped source dataframe.
110 columns : `list` [`str`]
111 List of columns to transfer.
112 default_dtype : `np.dtype`, optional
113 Default datatype (if not in target).
114 int_fill_value : `int`, optional
115 Fill value for integer columns to avoid pandas insisting
116 that everything should be float-ified as nans.
117 force_int_to_float : `bool`, optional
118 Force integer columns to float columns? Use this option
119 for backwards compatibility for old pandas misfeatures which
120 are expected by some downstream processes.
121 """
122 is_series = isinstance(source, pd.Series)
123 for col in columns:
124 if is_series:
125 source_col = source
126 else:
127 source_col = source[col]
129 matched_length = False
130 if col in target.columns:
131 target_dtype = target[col].dtype
132 matched_length = len(target) == len(source)
133 else:
134 target_dtype = default_dtype
136 if (matched_length or pd.api.types.is_float_dtype(target_dtype)) and not force_int_to_float:
137 # If we have a matched length or float, we can do a
138 # straight assignment here.
139 target.loc[:, col] = source_col.astype(target_dtype)
140 else:
141 # With mis-matched integers, we must do this dance to preserve types.
142 # Note that this may lose precision with very large numbers.
144 # Convert to float
145 target[col] = target[col].astype(np.float64)
146 # Set the column, casting to float.
147 target.loc[:, col] = source_col.astype(np.float64)
148 if not force_int_to_float:
149 # Convert back to integer
150 target[col] = target[col].fillna(int_fill_value).astype(target_dtype)
153def compute_optimized_periodogram_grid(x0, oversampling_factor=5, nyquist_factor=100):
154 """
155 Computes an optimized periodogram frequency grid for a given time series.
157 Parameters
158 ----------
159 x0 : `array`
160 The input time axis.
161 oversampling_factor : `int`, optional
162 The oversampling factor for frequency grid.
163 nyquist_factor : `int`, optional
164 The Nyquist factor for frequency grid.
166 Returns
167 -------
168 frequencies : `array`
169 The computed optimized periodogram frequency grid.
170 """
172 num_points = len(x0)
173 baseline = np.max(x0) - np.min(x0)
175 # Calculate the frequency resolution based on oversampling factor and baseline
176 frequency_resolution = 1. / baseline / oversampling_factor
178 num_frequencies = int(
179 0.5 * oversampling_factor * nyquist_factor * num_points)
180 frequencies = frequency_resolution + \
181 frequency_resolution * np.arange(num_frequencies)
183 return frequencies
186class LombScarglePeriodogramConfig(DiaObjectCalculationPluginConfig):
187 pass
190@register("ap_lombScarglePeriodogram")
191class LombScarglePeriodogram(DiaObjectCalculationPlugin):
192 """Compute the single-band period of a DiaObject given a set of DiaSources.
193 """
194 ConfigClass = LombScarglePeriodogramConfig
196 plugType = "multi"
197 outputCols = ["period", "power"]
198 needsFilter = True
200 @classmethod
201 def getExecutionOrder(cls):
202 return cls.DEFAULT_CATALOGCALCULATION
204 @catchWarnings(warns=["All-NaN slice encountered"])
205 def calculate(self,
206 diaObjects,
207 diaSources,
208 filterDiaSources,
209 band):
210 """Compute the periodogram.
212 Parameters
213 ----------
214 diaObjects : `pandas.DataFrame`
215 Summary objects to store values in.
216 diaSources : `pandas.DataFrame` or `pandas.DataFrameGroupBy`
217 Catalog of DiaSources summarized by this DiaObject.
218 """
220 # Check and initialize output columns in diaObjects.
221 if (periodCol := f"{band}_period") not in diaObjects.columns: 221 ↛ 223line 221 didn't jump to line 223 because the condition on line 221 was always true
222 diaObjects[periodCol] = np.nan
223 if (powerCol := f"{band}_power") not in diaObjects.columns: 223 ↛ 226line 223 didn't jump to line 226 because the condition on line 223 was always true
224 diaObjects[powerCol] = np.nan
226 def _calculate_period(df, min_detections=5, nterms=1, oversampling_factor=5, nyquist_factor=100):
227 """Compute the Lomb-Scargle periodogram given a set of DiaSources.
229 Parameters
230 ----------
231 df : `pandas.DataFrame`
232 The input DataFrame.
233 min_detections : `int`, optional
234 The minimum number of detections.
235 nterms : `int`, optional
236 The number of terms in the Lomb-Scargle model.
237 oversampling_factor : `int`, optional
238 The oversampling factor for frequency grid.
239 nyquist_factor : `int`, optional
240 The Nyquist factor for frequency grid.
242 Returns
243 -------
244 pd_tab : `pandas.Series`
245 The output DataFrame with the Lomb-Scargle parameters.
246 """
247 tmpDf = df[~np.logical_or(np.isnan(df["psfFlux"]),
248 np.isnan(df["midpointMjdTai"]))]
250 if len(tmpDf) < min_detections: 250 ↛ 251line 250 didn't jump to line 251 because the condition on line 250 was never true
251 return pd.Series({periodCol: np.nan, powerCol: np.nan})
253 time = tmpDf["midpointMjdTai"].to_numpy()
254 flux = tmpDf["psfFlux"].to_numpy()
255 flux_err = tmpDf["psfFluxErr"].to_numpy()
257 lsp = LombScargle(time, flux, dy=flux_err, nterms=nterms)
258 f_grid = compute_optimized_periodogram_grid(
259 time, oversampling_factor=oversampling_factor, nyquist_factor=nyquist_factor)
260 period = 1/f_grid
261 power = lsp.power(f_grid)
263 return pd.Series({periodCol: period[np.argmax(power)],
264 powerCol: np.max(power)})
266 diaObjects.loc[:, [periodCol, powerCol]
267 ] = filterDiaSources.apply(_calculate_period)
270class LombScarglePeriodogramMultiConfig(DiaObjectCalculationPluginConfig):
271 pass
274@register("ap_lombScarglePeriodogramMulti")
275class LombScarglePeriodogramMulti(DiaObjectCalculationPlugin):
276 """Compute the multi-band LombScargle periodogram of a DiaObject given a set of DiaSources.
277 """
278 ConfigClass = LombScarglePeriodogramMultiConfig
280 plugType = "multi"
281 outputCols = ["multiPeriod", "multiPower",
282 "multiFap", "multiAmp", "multiPhase"]
283 needsFilter = True
285 @classmethod
286 def getExecutionOrder(cls):
287 return cls.DEFAULT_CATALOGCALCULATION
289 @staticmethod
290 def calculate_baluev_fap(time, n, maxPeriod, zmax):
291 """Calculate the False-Alarm probability using the Baluev approximation.
293 Parameters
294 ----------
295 time : `array`
296 The input time axis.
297 n : `int`
298 The number of detections.
299 maxPeriod : `float`
300 The maximum period in the grid.
301 zmax : `float`
302 The maximum power in the grid.
304 Returns
305 -------
306 fap_estimate : `float`
307 The False-Alarm probability Baluev approximation.
309 Notes
310 ----------
311 .. [1] Baluev, R. V. 2008, MNRAS, 385, 1279
312 .. [2] Süveges, M., Guy, L.P., Eyer, L., et al. 2015, MNRAS, 450, 2052
313 """
314 if n <= 2: 314 ↛ 315line 314 didn't jump to line 315 because the condition on line 314 was never true
315 return np.nan
317 gam_ratio = math.factorial(int((n - 1)/2)) / math.factorial(int((n - 2)/2))
318 fu = 1/maxPeriod
319 return gam_ratio * np.sqrt(
320 4*np.pi*statistics.variance(time)
321 ) * fu * (1-zmax)**((n-4)/2) * np.sqrt(zmax)
323 @staticmethod
324 def generate_lsp_params(lsp_model, fbest, bands):
325 """Generate the Lomb-Scargle parameters.
326 Parameters
327 ----------
328 lsp_model : `astropy.timeseries.LombScargleMultiband`
329 The Lomb-Scargle model.
330 fbest : `float`
331 The best period.
332 bands : `array`
333 The bands of the time series.
335 Returns
336 -------
337 Amp : `array`
338 The amplitude of the time series.
339 Ph : `array`
340 The phase of the time series.
342 Notes
343 ----------
344 .. [1] VanderPlas, J. T., & Ivezic, Z. 2015, ApJ, 812, 18
345 """
346 best_params = lsp_model.model_parameters(fbest, units=True)
348 name_params = [f"theta_base_{i}" for i in range(3)]
349 name_params += [f"theta_band_{band}_{i}" for band in np.unique(bands) for i in range(3)]
351 df_params = pd.DataFrame([best_params], columns=name_params)
353 unique_bands = np.unique(bands)
355 amplitude_band = [np.sqrt(df_params[f"theta_band_{band}_1"]**2
356 + df_params[f"theta_band_{band}_2"]**2)
357 for band in unique_bands]
358 phase_bands = [np.arctan2(df_params[f"theta_band_{band}_2"],
359 df_params[f"theta_band_{band}_1"]) for band in unique_bands]
361 amp = [a[0] for a in amplitude_band]
362 ph = [p[0] for p in phase_bands]
364 return amp, ph
366 @catchWarnings(warns=["All-NaN slice encountered"])
367 def calculate(self,
368 diaObjects,
369 diaSources,
370 **kwargs):
371 """Compute the multi-band LombScargle periodogram of a DiaObject given
372 a set of DiaSources.
374 Parameters
375 ----------
376 diaObjects : `pandas.DataFrame`
377 Summary objects to store values in.
378 diaSources : `pandas.DataFrame` or `pandas.DataFrameGroupBy`
379 Catalog of DiaSources summarized by this DiaObject.
380 **kwargs : `dict`
381 Unused kwargs that are always passed to a plugin.
382 """
384 bands_arr = diaSources['band'].unique().values
385 unique_bands = np.unique(np.concatenate(bands_arr))
386 # Check and initialize output columns in diaObjects.
387 if (periodCol := "multiPeriod") not in diaObjects.columns: 387 ↛ 389line 387 didn't jump to line 389 because the condition on line 387 was always true
388 diaObjects[periodCol] = np.nan
389 if (powerCol := "multiPower") not in diaObjects.columns: 389 ↛ 391line 389 didn't jump to line 391 because the condition on line 389 was always true
390 diaObjects[powerCol] = np.nan
391 if (fapCol := "multiFap") not in diaObjects.columns: 391 ↛ 393line 391 didn't jump to line 393 because the condition on line 391 was always true
392 diaObjects[fapCol] = np.nan
393 ampCol = "multiAmp"
394 phaseCol = "multiPhase"
395 for i in range(len(unique_bands)):
396 ampCol_band = f"{unique_bands[i]}_{ampCol}"
397 if ampCol_band not in diaObjects.columns: 397 ↛ 399line 397 didn't jump to line 399 because the condition on line 397 was always true
398 diaObjects[ampCol_band] = np.nan
399 phaseCol_band = f"{unique_bands[i]}_{phaseCol}"
400 if phaseCol_band not in diaObjects.columns: 400 ↛ 395line 400 didn't jump to line 395 because the condition on line 400 was always true
401 diaObjects[phaseCol_band] = np.nan
403 def _calculate_period_multi(df, all_unique_bands,
404 min_detections=9, oversampling_factor=5, nyquist_factor=100):
405 """Calculate the multi-band Lomb-Scargle periodogram.
407 Parameters
408 ----------
409 df : `pandas.DataFrame`
410 The input DataFrame.
411 all_unique_bands : `list` of `str`
412 List of all bands present in the diaSource table that is being worked on.
413 min_detections : `int`, optional
414 The minimum number of detections, including all bands.
415 oversampling_factor : `int`, optional
416 The oversampling factor for frequency grid.
417 nyquist_factor : `int`, optional
418 The Nyquist factor for frequency grid.
420 Returns
421 -------
422 pd_tab : `pandas.Series`
423 The output DataFrame with the Lomb-Scargle parameters.
424 """
425 tmpDf = df[~np.logical_or(np.isnan(df["psfFlux"]),
426 np.isnan(df["midpointMjdTai"]))]
428 if (len(tmpDf)) < min_detections:
429 pd_tab_nodet = pd.Series({periodCol: np.nan,
430 powerCol: np.nan,
431 fapCol: np.nan})
432 for band in all_unique_bands:
433 pd_tab_nodet[f"{band}_{ampCol}"] = np.nan
434 pd_tab_nodet[f"{band}_{phaseCol}"] = np.nan
436 return pd_tab_nodet
438 time = tmpDf["midpointMjdTai"].to_numpy()
439 flux = tmpDf["psfFlux"].to_numpy()
440 flux_err = tmpDf["psfFluxErr"].to_numpy()
441 bands = tmpDf["band"].to_numpy()
443 lsp = LombScargleMultiband(time, flux, bands, dy=flux_err,
444 nterms_base=1, nterms_band=1)
446 f_grid = compute_optimized_periodogram_grid(
447 time, oversampling_factor=oversampling_factor, nyquist_factor=nyquist_factor)
448 period = 1/f_grid
449 power = lsp.power(f_grid)
451 fap_estimate = self.calculate_baluev_fap(
452 time, len(time), period[np.argmax(power)], np.max(power))
454 params_table_new = self.generate_lsp_params(lsp, f_grid[np.argmax(power)], bands)
456 pd_tab = pd.Series({periodCol: period[np.argmax(power)],
457 powerCol: np.max(power),
458 fapCol: fap_estimate
459 })
461 # Initialize the per-band amplitude/phase columns as NaNs
462 for band in all_unique_bands:
463 pd_tab[f"{band}_{ampCol}"] = np.nan
464 pd_tab[f"{band}_{phaseCol}"] = np.nan
466 # Populate the values of only the bands that have data for this diaSource
467 unique_bands = np.unique(bands)
468 for i in range(len(unique_bands)):
469 pd_tab[f"{unique_bands[i]}_{ampCol}"] = params_table_new[0][i]
470 pd_tab[f"{unique_bands[i]}_{phaseCol}"] = params_table_new[1][i]
472 return pd_tab
474 columns_list = [periodCol, powerCol, fapCol]
475 for i in range(len(unique_bands)):
476 columns_list.append(f"{unique_bands[i]}_{ampCol}")
477 columns_list.append(f"{unique_bands[i]}_{phaseCol}")
479 diaObjects.loc[:, columns_list
480 ] = diaSources.apply(_calculate_period_multi, unique_bands)
483class UnphysicalDiaSourceSeparation(pipeBase.AlgorithmError):
484 """Raised if associated DiaSources are unphysically separated.
486 Parameters
487 ----------
488 separation : `float`
489 Observed separation in arseconds.
490 max_allowed_separation : `float`
491 Configured maximum separation in arcseconds.
492 """
494 def __init__(self, separation, max_allowed_separation) -> None:
495 self._separation = separation
496 self._max_allowed_separation = max_allowed_separation
497 super().__init__(f"Observed DiaSource separation {separation} exceeds allowed value of "
498 f"{max_allowed_separation}")
500 @property
501 def metadata(self) -> dict:
502 return {
503 "separation": self._separation,
504 "max_allowed_separation": self._max_allowed_separation,
505 }
508class MeanDiaPositionConfig(DiaObjectCalculationPluginConfig):
509 MaxAllowedDiaSourceSeparation = pexConfig.Field(
510 dtype=float,
511 default=3.0,
512 doc="Max allowed separation of associated DiaSources in arcsec. "
513 "Raises if unphysical separation is found. "
514 )
517@register("ap_meanPosition")
518class MeanDiaPosition(DiaObjectCalculationPlugin):
519 """Compute the mean position and position uncertainty of a DiaObject
520 from its associated DiaSources.
522 The reported (ra, dec) is the inverse-variance weighted spherical
523 mean of the per-source positions, using only the DiaSources whose
524 per-source coordinate errors are finite. The reported covariance
525 is the covariance of that weighted-mean estimator,
527 C_formal = (sum_i inv(C_i))^-1,
529 where ``C_i`` is the per-source 2x2 covariance, multiplied by a
530 *chi-squared scale factor* that inflates the result when the
531 empirical scatter is larger than the per-source errors predict:
533 chi2 = sum_i r_i^T inv(C_i) r_i (r_i = tangent-plane
534 residual from the mean)
535 dof = 2 * (N - 1)
536 C_obj = max(1, chi2 / dof) * C_formal
538 When the data agree with the per-source errors (chi2 ~ dof) the
539 scale factor is unity and ``C_obj == C_formal``. When the residuals
540 exceed what the errors predict, the covariance is inflated.
541 """
543 ConfigClass = MeanDiaPositionConfig
545 plugType = 'multi'
547 outputCols = ["ra", "dec", "raErr", "decErr", "ra_dec_Cov"]
548 needsFilter = False
550 @classmethod
551 def getExecutionOrder(cls):
552 return cls.DEFAULT_CATALOGCALCULATION
554 def calculate(self, diaObjects, diaSources, **kwargs):
555 """Compute the mean ra/dec position and its uncertainty for each
556 DiaObject.
558 Parameters
559 ----------
560 diaObjects : `pandas.DataFrame`
561 Summary objects to store values in.
562 diaSources : `pandas.DataFrame` or `pandas.DataFrameGroupBy`
563 Catalog of DiaSources summarized by this DiaObject.
564 **kwargs
565 Any additional keyword arguments that may be passed to the plugin.
566 """
567 for outCol in self.outputCols:
568 if outCol not in diaObjects.columns: 568 ↛ 567line 568 didn't jump to line 567 because the condition on line 568 was always true
569 diaObjects[outCol] = np.nan
571 maxAllowedSep = self.config.MaxAllowedDiaSourceSeparation
573 def _column(df, name):
574 """Read a column as float64, or return all-NaN if absent."""
575 if name in df.columns:
576 return df[name].to_numpy(dtype=np.float64)
577 return np.full(len(df), np.nan)
579 def _computeMeanPos(df):
580 coords = list(geom.SpherePoint(src["ra"], src["dec"], geom.degrees)
581 for idx, src in df.iterrows())
582 # Quick check of the unweighted mean to catch unphysical objects
583 # and as a fallback.
584 unweightedAvg = geom.averageSpherePoint(coords)
585 maxSep = max(unweightedAvg.separation(coord).asArcseconds() for coord in coords)
586 if maxSep > maxAllowedSep:
587 raise UnphysicalDiaSourceSeparation(maxSep, maxAllowedSep)
589 nSrc = len(coords)
590 raErr = _column(df, "raErr")
591 decErr = _column(df, "decErr")
592 raDecCov = _column(df, "ra_dec_Cov")
593 finiteDiag = (np.isfinite(raErr) & (raErr > 0)) & (np.isfinite(decErr) & (decErr > 0))
594 nDiag = int(finiteDiag.sum())
596 # Tangent-plane offsets from the unweighted mean reference,
597 # in degrees (arc-length east, north).
598 offsets = np.array(
599 [[off.asDegrees() for off in unweightedAvg.getTangentPlaneOffset(coord)] for coord in coords]
600 ) if nSrc >= 1 else np.zeros((0, 2))
602 aveCoord = unweightedAvg
603 varRa = np.nan
604 varDec = np.nan
605 raDecCovObj = np.nan
607 if nDiag >= 1:
608 # Weighted-mean path: include only sources with finite
609 # per-source coordinate errors.
610 idx = np.where(finiteDiag)[0]
611 sigmaRaSq = raErr[idx]**2
612 sigmaDecSq = decErr[idx]**2
613 # Use full 2x2 weights only when every included source
614 # has a finite ra_dec_Cov; otherwise use diagonal weights
615 # for all of them and emit NaN for the output covariance.
616 haveFullCov = bool(np.isfinite(raDecCov[idx]).all())
618 C = np.zeros((nDiag, 2, 2), dtype=np.float64)
619 C[:, 0, 0] = sigmaRaSq
620 C[:, 1, 1] = sigmaDecSq
621 if haveFullCov:
622 C[:, 0, 1] = raDecCov[idx]
623 C[:, 1, 0] = raDecCov[idx]
624 try:
625 W = np.linalg.inv(C)
626 except np.linalg.LinAlgError:
627 # Singular per-source covariance: drop off-diagonal
628 # and retry with diagonal weights only.
629 haveFullCov = False
630 C[:, 0, 1] = 0.0
631 C[:, 1, 0] = 0.0
632 W = np.linalg.inv(C)
634 W_sum = W.sum(axis=0)
635 C_formal = np.linalg.inv(W_sum)
637 # Weighted-mean tangent-plane offset from the unweighted
638 # reference: mu_offset = (sum W_i)^-1 (sum W_i r_i). Then
639 # apply that offset to the unweighted mean to get the
640 # weighted spherical mean. SpherePoint.offset takes a
641 # bearing measured counter-clockwise from east (bearing 0
642 # is due east, bearing pi/2 is due north).
643 mu_offset = C_formal @ np.einsum('nij,nj->i', W, offsets[idx])
644 east, north = float(mu_offset[0]), float(mu_offset[1])
645 sep_deg = float(np.hypot(east, north))
646 if sep_deg > 0.0:
647 bearing = geom.Angle(float(np.arctan2(north, east)), geom.radians)
648 aveCoord = unweightedAvg.offset(bearing, sep_deg*geom.degrees)
650 varRaFormal = float(C_formal[0, 0])
651 varDecFormal = float(C_formal[1, 1])
652 covFormal = float(C_formal[0, 1]) if haveFullCov else np.nan
654 if nDiag == 1:
655 # Single included source: pass its own 2x2 covariance through.
656 varRa = varRaFormal
657 varDec = varDecFormal
658 raDecCovObj = covFormal
659 else:
660 # Inflate the variance if the reduced chi**2 fit of the
661 # diaSource coordinates is greater than 1, to account for
662 # scatter.
663 residuals = offsets[idx] - mu_offset
664 if haveFullCov:
665 # Vectorized Σₙ rₙᵀ Wₙ rₙ (r = residuals).
666 chi2 = float(np.einsum('ni,nij,nj->', residuals, W, residuals))
667 else:
668 chi2 = float(np.sum(residuals[:, 0]**2/sigmaRaSq + residuals[:, 1]**2/sigmaDecSq))
669 dof = 2*(nDiag - 1)
670 scale = max(1.0, chi2/dof)
671 varRa = varRaFormal*scale
672 varDec = varDecFormal*scale
673 raDecCovObj = (covFormal*scale if np.isfinite(covFormal) else np.nan)
674 else:
675 # No source has finite per-source coordinate errors. Fall
676 # back to the unweighted mean position plus the empirical
677 # SEM (when N >= 2), and warn that the per-source errors
678 # were not usable.
679 warnings.warn(
680 "No DiaSources with finite coordinate errors; falling "
681 "back to the unweighted mean position with scatter-only "
682 "uncertainty.",
683 stacklevel=2,
684 )
685 if nSrc >= 2:
686 sampleCov = np.cov(offsets, rowvar=False, ddof=1)
687 semCov = sampleCov/nSrc
688 varRa = float(semCov[0, 0])
689 varDec = float(semCov[1, 1])
690 raDecCovObj = float(semCov[0, 1])
691 # else: nSrc == 1 with no error -- outputs stay NaN.
693 raErrObj = float(np.sqrt(varRa)) if np.isfinite(varRa) else np.nan
694 decErrObj = float(np.sqrt(varDec)) if np.isfinite(varDec) else np.nan
696 return pd.Series({"ra": aveCoord.getRa().asDegrees(),
697 "dec": aveCoord.getDec().asDegrees(),
698 "raErr": raErrObj,
699 "decErr": decErrObj,
700 "ra_dec_Cov": raDecCovObj})
702 ans = diaSources.apply(_computeMeanPos)
703 typeSafePandasAssignment(diaObjects, ans, ["ra", "dec", "raErr", "decErr", "ra_dec_Cov"])
706class FirstLastDiaSourceTimeConfig(DiaObjectCalculationPluginConfig):
707 pass
710@register("ap_firstLastDiaSourceTime")
711class FirstLastDiaSourceTime(DiaObjectCalculationPlugin):
712 """Determine the min and max time DiaSources have been observed.
713 """
715 ConfigClass = FirstLastDiaSourceTimeConfig
717 plugType = 'multi'
719 outputCols = ["firstDiaSourceMjdTai", "lastDiaSourceMjdTai"]
720 needsFilter = False
722 @classmethod
723 def getExecutionOrder(cls):
724 return cls.DEFAULT_CATALOGCALCULATION
726 def calculate(self, diaObjects, diaSources, **kwargs):
727 """Store the first and last DiaSource times per DiaObject.
729 Stores the smaller of the stored value of ``firstDiaSourceMjdTai``
730 and the minimum ``midpointMjdTai`` in the provided DiaSources, and
731 the larger of the stored value of ``lastDiaSourceMjdTai`` and the
732 maximum ``midpointMjdTai`` in the provided DiaSources. This allows
733 the stored values to persist even when only a shorter history of
734 DiaSources is loaded from the APDB. DiaSources with NaN
735 ``midpointMjdTai`` are ignored.
737 Parameters
738 ----------
739 diaObjects : `pandas.DataFrame`
740 Summary objects to store values in, which may include prior
741 values of ``firstDiaSourceMjdTai`` and ``lastDiaSourceMjdTai``
742 for some objects.
743 diaSources : `pandas.DataFrame` or `pandas.DataFrameGroupBy`
744 Catalog of DiaSources summarized by these DiaObjects.
745 **kwargs
746 Any additional keyword arguments that may be passed to the plugin.
747 """
748 for outCol in self.outputCols:
749 if outCol not in diaObjects.columns:
750 diaObjects[outCol] = np.nan
752 times = diaSources["midpointMjdTai"]
753 # groupby min/max skips NaN values by default.
754 firstTime = times.min()
755 lastTime = times.max()
757 # Combine with any previously stored values; fmin/fmax ignore NaN
758 # unless both inputs are NaN.
759 storedFirst = diaObjects.loc[firstTime.index, "firstDiaSourceMjdTai"].to_numpy(dtype=np.float64)
760 storedLast = diaObjects.loc[lastTime.index, "lastDiaSourceMjdTai"].to_numpy(dtype=np.float64)
761 # Update only objects with DiaSources, and keep the stored values of
762 # all others.
763 diaObjects.loc[firstTime.index, "firstDiaSourceMjdTai"] = np.fmin(
764 storedFirst, firstTime.to_numpy(dtype=np.float64))
765 diaObjects.loc[lastTime.index, "lastDiaSourceMjdTai"] = np.fmax(
766 storedLast, lastTime.to_numpy(dtype=np.float64))
769class HTMIndexDiaPositionConfig(DiaObjectCalculationPluginConfig):
771 htmLevel = pexConfig.Field(
772 dtype=int,
773 doc="Level of the HTM pixelization.",
774 default=20,
775 )
778@register("ap_HTMIndex")
779class HTMIndexDiaPosition(DiaObjectCalculationPlugin):
780 """Compute the mean position of a DiaObject given a set of DiaSources.
782 Notes
783 -----
784 This plugin was implemented to satisfy requirements of old APDB interface
785 which required ``pixelId`` column in DiaObject with HTM20 index. APDB
786 interface had migrated to not need that information, but we keep this
787 plugin in case it may be useful for something else.
788 """
789 ConfigClass = HTMIndexDiaPositionConfig
791 plugType = 'single'
793 inputCols = ["ra", "dec"]
794 outputCols = ["pixelId"]
795 needsFilter = False
797 def __init__(self, config, name, metadata):
798 DiaObjectCalculationPlugin.__init__(self, config, name, metadata)
799 self.pixelator = sphgeom.HtmPixelization(self.config.htmLevel)
801 @classmethod
802 def getExecutionOrder(cls):
803 return cls.FLUX_MOMENTS_CALCULATED
805 def calculate(self, diaObjects, diaObjectId, **kwargs):
806 """Compute the mean position of a DiaObject given a set of DiaSources
808 Parameters
809 ----------
810 diaObjects : `pandas.dataFrame`
811 Summary objects to store values in and read ra/dec from.
812 diaObjectId : `int`
813 Id of the diaObject to update.
814 **kwargs
815 Any additional keyword arguments that may be passed to the plugin.
816 """
817 sphPoint = geom.SpherePoint(
818 diaObjects.at[diaObjectId, "ra"] * geom.degrees,
819 diaObjects.at[diaObjectId, "dec"] * geom.degrees)
820 diaObjects.at[diaObjectId, "pixelId"] = self.pixelator.index(
821 sphPoint.getVector())
824class NumDiaSourcesDiaPluginConfig(DiaObjectCalculationPluginConfig):
825 pass
828@register("ap_nDiaSources")
829class NumDiaSourcesDiaPlugin(DiaObjectCalculationPlugin):
830 """Compute the total number of DiaSources associated with this DiaObject.
831 """
833 ConfigClass = NumDiaSourcesDiaPluginConfig
834 outputCols = ["nDiaSources"]
835 plugType = "multi"
836 needsFilter = False
838 @classmethod
839 def getExecutionOrder(cls):
840 return cls.DEFAULT_CATALOGCALCULATION
842 def calculate(self, diaObjects, diaSources, **kwargs):
843 """Compute the total number of DiaSources associated with this DiaObject.
845 Parameters
846 ----------
847 diaObject : `dict`
848 Summary object to store values in and read ra/dec from.
849 **kwargs
850 Any additional keyword arguments that may be passed to the plugin.
851 """
852 typeSafePandasAssignment(diaObjects, diaSources.diaObjectId.count(), ["nDiaSources"])
855class SimpleSourceFlagDiaPluginConfig(DiaObjectCalculationPluginConfig):
856 pass
859@register("ap_diaObjectFlag")
860class SimpleSourceFlagDiaPlugin(DiaObjectCalculationPlugin):
861 """Find if any DiaSource is flagged.
863 Set the DiaObject flag if any DiaSource is flagged.
864 """
866 ConfigClass = NumDiaSourcesDiaPluginConfig
867 outputCols = ["flags"]
868 plugType = "multi"
869 needsFilter = False
871 @classmethod
872 def getExecutionOrder(cls):
873 return cls.DEFAULT_CATALOGCALCULATION
875 def calculate(self, diaObjects, diaSources, **kwargs):
876 """Find if any DiaSource is flagged.
878 Set the DiaObject flag if any DiaSource is flagged.
880 Parameters
881 ----------
882 diaObject : `dict`
883 Summary object to store values in and read ra/dec from.
884 **kwargs
885 Any additional keyword arguments that may be passed to the plugin.
886 """
887 typeSafePandasAssignment(diaObjects, diaSources.flags.any(), ["flags"], default_dtype=np.uint64)
890class WeightedMeanDiaPsfFluxConfig(DiaObjectCalculationPluginConfig):
891 pass
894@register("ap_meanFlux")
895class WeightedMeanDiaPsfFlux(DiaObjectCalculationPlugin):
896 """Compute the weighted mean and mean error on the point source fluxes
897 of the DiaSource measured on the difference image.
899 Additionally store number of usable data points.
900 """
902 ConfigClass = WeightedMeanDiaPsfFluxConfig
903 outputCols = ["psfFluxMean", "psfFluxMeanErr", "psfFluxNdata"]
904 plugType = "multi"
905 needsFilter = True
907 @classmethod
908 def getExecutionOrder(cls):
909 return cls.DEFAULT_CATALOGCALCULATION
911 @catchWarnings(warns=["invalid value encountered",
912 "divide by zero"])
913 def calculate(self,
914 diaObjects,
915 diaSources,
916 filterDiaSources,
917 band,
918 **kwargs):
919 """Compute the weighted mean and mean error of the point source flux.
921 Parameters
922 ----------
923 diaObject : `dict`
924 Summary object to store values in.
925 diaSources : `pandas.DataFrame`
926 DataFrame representing all diaSources associated with this
927 diaObject.
928 filterDiaSources : `pandas.DataFrame`
929 DataFrame representing diaSources associated with this
930 diaObject that are observed in the band pass ``band``.
931 band : `str`
932 Simple, string name of the filter for the flux being calculated.
933 **kwargs
934 Any additional keyword arguments that may be passed to the plugin.
935 """
936 meanName = "{}_psfFluxMean".format(band)
937 errName = "{}_psfFluxMeanErr".format(band)
938 nDataName = "{}_psfFluxNdata".format(band)
939 if meanName not in diaObjects.columns: 939 ↛ 941line 939 didn't jump to line 941 because the condition on line 939 was always true
940 diaObjects[meanName] = np.nan
941 if errName not in diaObjects.columns: 941 ↛ 943line 941 didn't jump to line 943 because the condition on line 941 was always true
942 diaObjects[errName] = np.nan
943 if nDataName not in diaObjects.columns: 943 ↛ 946line 943 didn't jump to line 946 because the condition on line 943 was always true
944 diaObjects[nDataName] = 0
946 def _weightedMean(df):
947 tmpDf = df[~np.logical_or(np.isnan(df["psfFlux"]),
948 np.isnan(df["psfFluxErr"]))]
949 tot_weight = np.nansum(1 / tmpDf["psfFluxErr"] ** 2)
950 fluxMean = np.nansum(tmpDf["psfFlux"]
951 / tmpDf["psfFluxErr"] ** 2)
952 fluxMean /= tot_weight
953 if tot_weight > 0: 953 ↛ 956line 953 didn't jump to line 956 because the condition on line 953 was always true
954 fluxMeanErr = np.sqrt(1 / tot_weight)
955 else:
956 fluxMeanErr = np.nan
957 nFluxData = len(tmpDf)
959 return pd.Series({meanName: fluxMean,
960 errName: fluxMeanErr,
961 nDataName: nFluxData},
962 dtype="object")
963 df = filterDiaSources.apply(_weightedMean).astype(diaObjects.dtypes[[meanName, errName, nDataName]])
964 # TODO DM-53254: Remove the force_int_to_float hack.
965 typeSafePandasAssignment(diaObjects, df, [meanName, errName, nDataName], force_int_to_float=True)
968class PercentileDiaPsfFluxConfig(DiaObjectCalculationPluginConfig):
969 percentiles = pexConfig.ListField(
970 dtype=int,
971 default=[5, 25, 50, 75, 95],
972 doc="Percentiles to calculate to compute values for. Should be "
973 "integer values."
974 )
977@register("ap_percentileFlux")
978class PercentileDiaPsfFlux(DiaObjectCalculationPlugin):
979 """Compute percentiles of diaSource fluxes.
980 """
982 ConfigClass = PercentileDiaPsfFluxConfig
983 # Output columns are created upon instantiation of the class.
984 outputCols = []
985 plugType = "multi"
986 needsFilter = True
988 def __init__(self, config, name, metadata, **kwargs):
989 DiaObjectCalculationPlugin.__init__(self,
990 config,
991 name,
992 metadata,
993 **kwargs)
994 self.outputCols = ["psfFluxPercentile{:02d}".format(percent)
995 for percent in self.config.percentiles]
997 @classmethod
998 def getExecutionOrder(cls):
999 return cls.DEFAULT_CATALOGCALCULATION
1001 @catchWarnings(warns=["All-NaN slice encountered"])
1002 def calculate(self,
1003 diaObjects,
1004 diaSources,
1005 filterDiaSources,
1006 band,
1007 **kwargs):
1008 """Compute the percentile fluxes of the point source flux.
1010 Parameters
1011 ----------
1012 diaObject : `dict`
1013 Summary object to store values in.
1014 diaSources : `pandas.DataFrame`
1015 DataFrame representing all diaSources associated with this
1016 diaObject.
1017 filterDiaSources : `pandas.DataFrame`
1018 DataFrame representing diaSources associated with this
1019 diaObject that are observed in the band pass ``band``.
1020 band : `str`
1021 Simple, string name of the filter for the flux being calculated.
1022 **kwargs
1023 Any additional keyword arguments that may be passed to the plugin.
1024 """
1025 pTileNames = []
1026 for tilePercent in self.config.percentiles:
1027 pTileName = "{}_psfFluxPercentile{:02d}".format(band,
1028 tilePercent)
1029 pTileNames.append(pTileName)
1030 if pTileName not in diaObjects.columns: 1030 ↛ 1026line 1030 didn't jump to line 1026 because the condition on line 1030 was always true
1031 diaObjects[pTileName] = np.nan
1033 def _fluxPercentiles(df):
1034 pTiles = np.nanpercentile(df["psfFlux"], self.config.percentiles)
1035 return pd.Series(
1036 dict((tileName, pTile)
1037 for tileName, pTile in zip(pTileNames, pTiles)))
1039 typeSafePandasAssignment(
1040 diaObjects,
1041 filterDiaSources.apply(_fluxPercentiles),
1042 pTileNames,
1043 )
1046class SigmaDiaPsfFluxConfig(DiaObjectCalculationPluginConfig):
1047 pass
1050@register("ap_sigmaFlux")
1051class SigmaDiaPsfFlux(DiaObjectCalculationPlugin):
1052 """Compute scatter of diaSource fluxes.
1053 """
1055 ConfigClass = SigmaDiaPsfFluxConfig
1056 # Output columns are created upon instantiation of the class.
1057 outputCols = ["psfFluxSigma"]
1058 plugType = "multi"
1059 needsFilter = True
1061 @classmethod
1062 def getExecutionOrder(cls):
1063 return cls.DEFAULT_CATALOGCALCULATION
1065 def calculate(self,
1066 diaObjects,
1067 diaSources,
1068 filterDiaSources,
1069 band,
1070 **kwargs):
1071 """Compute the sigma fluxes of the point source flux.
1073 Parameters
1074 ----------
1075 diaObject : `dict`
1076 Summary object to store values in.
1077 diaSources : `pandas.DataFrame`
1078 DataFrame representing all diaSources associated with this
1079 diaObject.
1080 filterDiaSources : `pandas.DataFrame`
1081 DataFrame representing diaSources associated with this
1082 diaObject that are observed in the band pass ``band``.
1083 band : `str`
1084 Simple, string name of the filter for the flux being calculated.
1085 **kwargs
1086 Any additional keyword arguments that may be passed to the plugin.
1087 """
1088 # Set "delta degrees of freedom (ddf)" to 1 to calculate the unbiased
1089 # estimator of scatter (i.e. 'N - 1' instead of 'N').
1090 column = "{}_psfFluxSigma".format(band)
1092 typeSafePandasAssignment(
1093 diaObjects,
1094 filterDiaSources.psfFlux.std(),
1095 [column],
1096 )
1099class Chi2DiaPsfFluxConfig(DiaObjectCalculationPluginConfig):
1100 pass
1103@register("ap_chi2Flux")
1104class Chi2DiaPsfFlux(DiaObjectCalculationPlugin):
1105 """Compute chi2 of diaSource fluxes.
1106 """
1108 ConfigClass = Chi2DiaPsfFluxConfig
1110 # Required input Cols
1111 inputCols = ["psfFluxMean"]
1112 # Output columns are created upon instantiation of the class.
1113 outputCols = ["psfFluxChi2"]
1114 plugType = "multi"
1115 needsFilter = True
1117 @classmethod
1118 def getExecutionOrder(cls):
1119 return cls.FLUX_MOMENTS_CALCULATED
1121 @catchWarnings(warns=["All-NaN slice encountered"])
1122 def calculate(self,
1123 diaObjects,
1124 diaSources,
1125 filterDiaSources,
1126 band,
1127 **kwargs):
1128 """Compute the chi2 of the point source fluxes.
1130 Parameters
1131 ----------
1132 diaObject : `dict`
1133 Summary object to store values in.
1134 diaSources : `pandas.DataFrame`
1135 DataFrame representing all diaSources associated with this
1136 diaObject.
1137 filterDiaSources : `pandas.DataFrame`
1138 DataFrame representing diaSources associated with this
1139 diaObject that are observed in the band pass ``band``.
1140 band : `str`
1141 Simple, string name of the filter for the flux being calculated.
1142 **kwargs
1143 Any additional keyword arguments that may be passed to the plugin.
1144 """
1145 meanName = "{}_psfFluxMean".format(band)
1146 column = "{}_psfFluxChi2".format(band)
1148 def _chi2(df):
1149 delta = (df["psfFlux"]
1150 - diaObjects.at[df.diaObjectId.iat[0], meanName])
1151 return np.nansum((delta / df["psfFluxErr"]) ** 2)
1153 typeSafePandasAssignment(
1154 diaObjects,
1155 filterDiaSources.apply(_chi2),
1156 [column],
1157 )
1160class MadDiaPsfFluxConfig(DiaObjectCalculationPluginConfig):
1161 pass
1164@register("ap_madFlux")
1165class MadDiaPsfFlux(DiaObjectCalculationPlugin):
1166 """Compute median absolute deviation of diaSource fluxes.
1167 """
1169 ConfigClass = MadDiaPsfFluxConfig
1171 # Required input Cols
1172 # Output columns are created upon instantiation of the class.
1173 outputCols = ["psfFluxMAD"]
1174 plugType = "multi"
1175 needsFilter = True
1177 @classmethod
1178 def getExecutionOrder(cls):
1179 return cls.DEFAULT_CATALOGCALCULATION
1181 @catchWarnings(warns=["All-NaN slice encountered"])
1182 def calculate(self,
1183 diaObjects,
1184 diaSources,
1185 filterDiaSources,
1186 band,
1187 **kwargs):
1188 """Compute the median absolute deviation of the point source fluxes.
1190 Parameters
1191 ----------
1192 diaObject : `dict`
1193 Summary object to store values in.
1194 diaSources : `pandas.DataFrame`
1195 DataFrame representing all diaSources associated with this
1196 diaObject.
1197 filterDiaSources : `pandas.DataFrame`
1198 DataFrame representing diaSources associated with this
1199 diaObject that are observed in the band pass ``band``.
1200 band : `str`
1201 Simple, string name of the filter for the flux being calculated.
1202 **kwargs
1203 Any additional keyword arguments that may be passed to the plugin.
1204 """
1205 column = "{}_psfFluxMAD".format(band)
1207 typeSafePandasAssignment(
1208 diaObjects,
1209 filterDiaSources.psfFlux.apply(median_absolute_deviation, ignore_nan=True),
1210 [column],
1211 )
1214class SkewDiaPsfFluxConfig(DiaObjectCalculationPluginConfig):
1215 pass
1218@register("ap_skewFlux")
1219class SkewDiaPsfFlux(DiaObjectCalculationPlugin):
1220 """Compute the skew of diaSource fluxes.
1221 """
1223 ConfigClass = SkewDiaPsfFluxConfig
1225 # Required input Cols
1226 # Output columns are created upon instantiation of the class.
1227 outputCols = ["psfFluxSkew"]
1228 plugType = "multi"
1229 needsFilter = True
1231 @classmethod
1232 def getExecutionOrder(cls):
1233 return cls.DEFAULT_CATALOGCALCULATION
1235 def calculate(self,
1236 diaObjects,
1237 diaSources,
1238 filterDiaSources,
1239 band,
1240 **kwargs):
1241 """Compute the skew of the point source fluxes.
1243 Parameters
1244 ----------
1245 diaObject : `dict`
1246 Summary object to store values in.
1247 diaSources : `pandas.DataFrame`
1248 DataFrame representing all diaSources associated with this
1249 diaObject.
1250 filterDiaSources : `pandas.DataFrame`
1251 DataFrame representing diaSources associated with this
1252 diaObject that are observed in the band pass ``band``.
1253 band : `str`
1254 Simple, string name of the filter for the flux being calculated.
1255 **kwargs
1256 Any additional keyword arguments that may be passed to the plugin.
1257 """
1258 column = "{}_psfFluxSkew".format(band)
1260 typeSafePandasAssignment(
1261 diaObjects,
1262 filterDiaSources.psfFlux.skew(),
1263 [column],
1264 )
1267class MinMaxDiaPsfFluxConfig(DiaObjectCalculationPluginConfig):
1268 pass
1271@register("ap_minMaxFlux")
1272class MinMaxDiaPsfFlux(DiaObjectCalculationPlugin):
1273 """Compute min/max of diaSource fluxes.
1274 """
1276 ConfigClass = MinMaxDiaPsfFluxConfig
1278 # Required input Cols
1279 # Output columns are created upon instantiation of the class.
1280 outputCols = ["psfFluxMin", "psfFluxMax"]
1281 plugType = "multi"
1282 needsFilter = True
1284 @classmethod
1285 def getExecutionOrder(cls):
1286 return cls.DEFAULT_CATALOGCALCULATION
1288 def calculate(self,
1289 diaObjects,
1290 diaSources,
1291 filterDiaSources,
1292 band,
1293 **kwargs):
1294 """Compute min/max of the point source fluxes.
1296 Parameters
1297 ----------
1298 diaObject : `dict`
1299 Summary object to store values in.
1300 diaSources : `pandas.DataFrame`
1301 DataFrame representing all diaSources associated with this
1302 diaObject.
1303 filterDiaSources : `pandas.DataFrame`
1304 DataFrame representing diaSources associated with this
1305 diaObject that are observed in the band pass ``band``.
1306 band : `str`
1307 Simple, string name of the filter for the flux being calculated.
1308 **kwargs
1309 Any additional keyword arguments that may be passed to the plugin.
1310 """
1311 minName = "{}_psfFluxMin".format(band)
1312 if minName not in diaObjects.columns: 1312 ↛ 1314line 1312 didn't jump to line 1314 because the condition on line 1312 was always true
1313 diaObjects[minName] = np.nan
1314 maxName = "{}_psfFluxMax".format(band)
1315 if maxName not in diaObjects.columns: 1315 ↛ 1318line 1315 didn't jump to line 1318 because the condition on line 1315 was always true
1316 diaObjects[maxName] = np.nan
1318 typeSafePandasAssignment(
1319 diaObjects,
1320 filterDiaSources.psfFlux.min(),
1321 [minName],
1322 )
1323 typeSafePandasAssignment(
1324 diaObjects,
1325 filterDiaSources.psfFlux.max(),
1326 [maxName],
1327 )
1330class MaxSlopeDiaPsfFluxConfig(DiaObjectCalculationPluginConfig):
1331 pass
1334@register("ap_maxSlopeFlux")
1335class MaxSlopeDiaPsfFlux(DiaObjectCalculationPlugin):
1336 """Compute the maximum ratio time ordered deltaFlux / deltaTime.
1337 """
1339 ConfigClass = MinMaxDiaPsfFluxConfig
1341 # Required input Cols
1342 # Output columns are created upon instantiation of the class.
1343 outputCols = ["psfFluxMaxSlope"]
1344 plugType = "multi"
1345 needsFilter = True
1347 @classmethod
1348 def getExecutionOrder(cls):
1349 return cls.DEFAULT_CATALOGCALCULATION
1351 def calculate(self,
1352 diaObjects,
1353 diaSources,
1354 filterDiaSources,
1355 band,
1356 **kwargs):
1357 """Compute the maximum ratio time ordered deltaFlux / deltaTime.
1359 Parameters
1360 ----------
1361 diaObject : `dict`
1362 Summary object to store values in.
1363 diaSources : `pandas.DataFrame`
1364 DataFrame representing all diaSources associated with this
1365 diaObject.
1366 filterDiaSources : `pandas.DataFrame`
1367 DataFrame representing diaSources associated with this
1368 diaObject that are observed in the band pass ``band``.
1369 band : `str`
1370 Simple, string name of the filter for the flux being calculated.
1371 **kwargs
1372 Any additional keyword arguments that may be passed to the plugin.
1373 """
1375 def _maxSlope(df):
1376 tmpDf = df[~np.logical_or(np.isnan(df["psfFlux"]),
1377 np.isnan(df["midpointMjdTai"]))]
1378 if len(tmpDf) < 2:
1379 return np.nan
1380 times = tmpDf["midpointMjdTai"].to_numpy()
1381 timeArgs = times.argsort()
1382 times = times[timeArgs]
1383 fluxes = tmpDf["psfFlux"].to_numpy()[timeArgs]
1384 return (np.diff(fluxes) / np.diff(times)).max()
1386 column = "{}_psfFluxMaxSlope".format(band)
1388 typeSafePandasAssignment(
1389 diaObjects,
1390 filterDiaSources.apply(_maxSlope),
1391 [column],
1392 )
1395class ErrMeanDiaPsfFluxConfig(DiaObjectCalculationPluginConfig):
1396 pass
1399@register("ap_meanErrFlux")
1400class ErrMeanDiaPsfFlux(DiaObjectCalculationPlugin):
1401 """Compute the mean of the dia source errors.
1402 """
1404 ConfigClass = ErrMeanDiaPsfFluxConfig
1406 # Required input Cols
1407 # Output columns are created upon instantiation of the class.
1408 outputCols = ["psfFluxErrMean"]
1409 plugType = "multi"
1410 needsFilter = True
1412 @classmethod
1413 def getExecutionOrder(cls):
1414 return cls.DEFAULT_CATALOGCALCULATION
1416 def calculate(self,
1417 diaObjects,
1418 diaSources,
1419 filterDiaSources,
1420 band,
1421 **kwargs):
1422 """Compute the mean of the dia source errors.
1424 Parameters
1425 ----------
1426 diaObject : `dict`
1427 Summary object to store values in.
1428 diaSources : `pandas.DataFrame`
1429 DataFrame representing all diaSources associated with this
1430 diaObject.
1431 filterDiaSources : `pandas.DataFrame`
1432 DataFrame representing diaSources associated with this
1433 diaObject that are observed in the band pass ``band``.
1434 band : `str`
1435 Simple, string name of the filter for the flux being calculated.
1436 **kwargs
1437 Any additional keyword arguments that may be passed to the plugin.
1438 """
1439 column = "{}_psfFluxErrMean".format(band)
1441 typeSafePandasAssignment(
1442 diaObjects,
1443 filterDiaSources.psfFluxErr.mean(),
1444 [column],
1445 # Note that the schemas expect this to be single-precision.
1446 default_dtype=np.float32,
1447 )
1450class LinearFitDiaPsfFluxConfig(DiaObjectCalculationPluginConfig):
1451 pass
1454@register("ap_linearFit")
1455class LinearFitDiaPsfFlux(DiaObjectCalculationPlugin):
1456 """Compute fit a linear model to flux vs time.
1457 """
1459 ConfigClass = LinearFitDiaPsfFluxConfig
1461 # Required input Cols
1462 # Output columns are created upon instantiation of the class.
1463 outputCols = ["psfFluxLinearSlope", "psfFluxLinearIntercept"]
1464 plugType = "multi"
1465 needsFilter = True
1467 @classmethod
1468 def getExecutionOrder(cls):
1469 return cls.DEFAULT_CATALOGCALCULATION
1471 def calculate(self,
1472 diaObjects,
1473 diaSources,
1474 filterDiaSources,
1475 band,
1476 **kwargs):
1477 """Compute fit a linear model to flux vs time.
1479 Parameters
1480 ----------
1481 diaObject : `dict`
1482 Summary object to store values in.
1483 diaSources : `pandas.DataFrame`
1484 DataFrame representing all diaSources associated with this
1485 diaObject.
1486 filterDiaSources : `pandas.DataFrame`
1487 DataFrame representing diaSources associated with this
1488 diaObject that are observed in the band pass ``band``.
1489 band : `str`
1490 Simple, string name of the filter for the flux being calculated.
1491 **kwargs
1492 Any additional keyword arguments that may be passed to the plugin.
1493 """
1495 mName = "{}_psfFluxLinearSlope".format(band)
1496 if mName not in diaObjects.columns: 1496 ↛ 1498line 1496 didn't jump to line 1498 because the condition on line 1496 was always true
1497 diaObjects[mName] = np.nan
1498 bName = "{}_psfFluxLinearIntercept".format(band)
1499 if bName not in diaObjects.columns: 1499 ↛ 1501line 1499 didn't jump to line 1501 because the condition on line 1499 was always true
1500 diaObjects[bName] = np.nan
1501 dtype = diaObjects[mName].dtype
1503 def _linearFit(df):
1504 tmpDf = df[~np.logical_or(
1505 np.isnan(df["psfFlux"]),
1506 np.logical_or(np.isnan(df["psfFluxErr"]),
1507 np.isnan(df["midpointMjdTai"])))]
1508 if len(tmpDf) < 2: 1508 ↛ 1509line 1508 didn't jump to line 1509 because the condition on line 1508 was never true
1509 return pd.Series({mName: np.nan, bName: np.nan})
1510 fluxes = tmpDf["psfFlux"].to_numpy()
1511 errors = tmpDf["psfFluxErr"].to_numpy()
1512 times = tmpDf["midpointMjdTai"].to_numpy()
1513 A = np.array([times / errors, 1 / errors]).transpose()
1514 m, b = lsq_linear(A, fluxes / errors).x
1515 return pd.Series({mName: m, bName: b}, dtype=dtype)
1517 typeSafePandasAssignment(
1518 diaObjects,
1519 filterDiaSources.apply(_linearFit),
1520 [mName, bName],
1521 )
1524class StetsonJDiaPsfFluxConfig(DiaObjectCalculationPluginConfig):
1525 pass
1528@register("ap_stetsonJ")
1529class StetsonJDiaPsfFlux(DiaObjectCalculationPlugin):
1530 """Compute the StetsonJ statistic on the DIA point source fluxes.
1531 """
1533 ConfigClass = LinearFitDiaPsfFluxConfig
1535 # Required input Cols
1536 inputCols = ["psfFluxMean"]
1537 # Output columns are created upon instantiation of the class.
1538 outputCols = ["psfFluxStetsonJ"]
1539 plugType = "multi"
1540 needsFilter = True
1542 @classmethod
1543 def getExecutionOrder(cls):
1544 return cls.FLUX_MOMENTS_CALCULATED
1546 def calculate(self,
1547 diaObjects,
1548 diaSources,
1549 filterDiaSources,
1550 band,
1551 **kwargs):
1552 """Compute the StetsonJ statistic on the DIA point source fluxes.
1554 Parameters
1555 ----------
1556 diaObject : `dict`
1557 Summary object to store values in.
1558 diaSources : `pandas.DataFrame`
1559 DataFrame representing all diaSources associated with this
1560 diaObject.
1561 filterDiaSources : `pandas.DataFrame`
1562 DataFrame representing diaSources associated with this
1563 diaObject that are observed in the band pass ``band``.
1564 band : `str`
1565 Simple, string name of the filter for the flux being calculated.
1566 **kwargs
1567 Any additional keyword arguments that may be passed to the plugin.
1568 """
1569 meanName = "{}_psfFluxMean".format(band)
1571 def _stetsonJ(df):
1572 tmpDf = df[~np.logical_or(np.isnan(df["psfFlux"]),
1573 np.isnan(df["psfFluxErr"]))]
1574 if len(tmpDf) < 2:
1575 return np.nan
1576 fluxes = tmpDf["psfFlux"].to_numpy()
1577 errors = tmpDf["psfFluxErr"].to_numpy()
1579 return self._stetson_J(
1580 fluxes,
1581 errors,
1582 diaObjects.at[tmpDf.diaObjectId.iat[0], meanName])
1584 column = "{}_psfFluxStetsonJ".format(band)
1585 typeSafePandasAssignment(
1586 diaObjects,
1587 filterDiaSources.apply(_stetsonJ),
1588 [column],
1589 )
1591 def _stetson_J(self, fluxes, errors, mean=None):
1592 """Compute the single band stetsonJ statistic.
1594 Parameters
1595 ----------
1596 fluxes : `numpy.ndarray` (N,)
1597 Calibrated lightcurve flux values.
1598 errors : `numpy.ndarray` (N,)
1599 Errors on the calibrated lightcurve fluxes.
1600 mean : `float`
1601 Starting mean from previous plugin.
1603 Returns
1604 -------
1605 stetsonJ : `float`
1606 stetsonJ statistic for the input fluxes and errors.
1608 References
1609 ----------
1610 .. [1] Stetson, P. B., "On the Automatic Determination of Light-Curve
1611 Parameters for Cepheid Variables", PASP, 108, 851S, 1996
1612 """
1613 n_points = len(fluxes)
1614 flux_mean = self._stetson_mean(fluxes, errors, mean)
1615 delta_val = (
1616 np.sqrt(n_points / (n_points - 1)) * (fluxes - flux_mean) / errors)
1617 p_k = delta_val ** 2 - 1
1619 return np.mean(np.sign(p_k) * np.sqrt(np.fabs(p_k)))
1621 def _stetson_mean(self,
1622 values,
1623 errors,
1624 mean=None,
1625 alpha=2.,
1626 beta=2.,
1627 n_iter=20,
1628 tol=1e-6):
1629 """Compute the stetson mean of the fluxes which down-weights outliers.
1631 Weighted biased on an error weighted difference scaled by a constant
1632 (1/``a``) and raised to the power beta. Higher betas more harshly
1633 penalize outliers and ``a`` sets the number of sigma where a weighted
1634 difference of 1 occurs.
1636 Parameters
1637 ----------
1638 values : `numpy.dnarray`, (N,)
1639 Input values to compute the mean of.
1640 errors : `numpy.ndarray`, (N,)
1641 Errors on the input values.
1642 mean : `float`
1643 Starting mean value or None.
1644 alpha : `float`
1645 Scalar down-weighting of the fractional difference. lower->more
1646 clipping. (Default value is 2.)
1647 beta : `float`
1648 Power law slope of the used to down-weight outliers. higher->more
1649 clipping. (Default value is 2.)
1650 n_iter : `int`
1651 Number of iterations of clipping.
1652 tol : `float`
1653 Fractional and absolute tolerance goal on the change in the mean
1654 before exiting early. (Default value is 1e-6)
1656 Returns
1657 -------
1658 mean : `float`
1659 Weighted stetson mean result.
1661 References
1662 ----------
1663 .. [1] Stetson, P. B., "On the Automatic Determination of Light-Curve
1664 Parameters for Cepheid Variables", PASP, 108, 851S, 1996
1665 """
1666 n_points = len(values)
1667 n_factor = np.sqrt(n_points / (n_points - 1))
1668 inv_var = 1 / errors ** 2
1670 if mean is None: 1670 ↛ 1671line 1670 didn't jump to line 1671 because the condition on line 1670 was never true
1671 mean = np.average(values, weights=inv_var)
1672 for iter_idx in range(n_iter): 1672 ↛ 1681line 1672 didn't jump to line 1681 because the loop on line 1672 didn't complete
1673 chi = np.fabs(n_factor * (values - mean) / errors)
1674 tmp_mean = np.average(
1675 values,
1676 weights=inv_var / (1 + (chi / alpha) ** beta))
1677 diff = np.fabs(tmp_mean - mean)
1678 mean = tmp_mean
1679 if diff / mean < tol and diff < tol:
1680 break
1681 return mean
1684class WeightedMeanDiaTotFluxConfig(DiaObjectCalculationPluginConfig):
1685 pass
1688@register("ap_meanTotFlux")
1689class WeightedMeanDiaTotFlux(DiaObjectCalculationPlugin):
1690 """Compute the weighted mean and mean error on the point source fluxes
1691 forced photometered at the DiaSource location in the calibrated image.
1693 Additionally store number of usable data points.
1694 """
1696 ConfigClass = WeightedMeanDiaPsfFluxConfig
1697 outputCols = ["scienceFluxMean", "scienceFluxMeanErr"]
1698 plugType = "multi"
1699 needsFilter = True
1701 @classmethod
1702 def getExecutionOrder(cls):
1703 return cls.DEFAULT_CATALOGCALCULATION
1705 @catchWarnings(warns=["invalid value encountered",
1706 "divide by zero"])
1707 def calculate(self,
1708 diaObjects,
1709 diaSources,
1710 filterDiaSources,
1711 band,
1712 **kwargs):
1713 """Compute the weighted mean and mean error of the point source flux.
1715 Parameters
1716 ----------
1717 diaObject : `dict`
1718 Summary object to store values in.
1719 diaSources : `pandas.DataFrame`
1720 DataFrame representing all diaSources associated with this
1721 diaObject.
1722 filterDiaSources : `pandas.DataFrame`
1723 DataFrame representing diaSources associated with this
1724 diaObject that are observed in the band pass ``band``.
1725 band : `str`
1726 Simple, string name of the filter for the flux being calculated.
1727 **kwargs
1728 Any additional keyword arguments that may be passed to the plugin.
1729 """
1730 totMeanName = "{}_scienceFluxMean".format(band)
1731 if totMeanName not in diaObjects.columns: 1731 ↛ 1733line 1731 didn't jump to line 1733 because the condition on line 1731 was always true
1732 diaObjects[totMeanName] = np.nan
1733 totErrName = "{}_scienceFluxMeanErr".format(band)
1734 if totErrName not in diaObjects.columns: 1734 ↛ 1737line 1734 didn't jump to line 1737 because the condition on line 1734 was always true
1735 diaObjects[totErrName] = np.nan
1737 def _meanFlux(df):
1738 tmpDf = df[~np.logical_or(np.isnan(df["scienceFlux"]),
1739 np.isnan(df["scienceFluxErr"]))]
1740 tot_weight = np.nansum(1 / tmpDf["scienceFluxErr"] ** 2)
1741 fluxMean = np.nansum(tmpDf["scienceFlux"]
1742 / tmpDf["scienceFluxErr"] ** 2)
1743 fluxMean /= tot_weight
1744 fluxMeanErr = np.sqrt(1 / tot_weight)
1746 return pd.Series({totMeanName: fluxMean,
1747 totErrName: fluxMeanErr})
1749 df = filterDiaSources.apply(_meanFlux).astype(diaObjects.dtypes[[totMeanName, totErrName]])
1750 typeSafePandasAssignment(diaObjects, df, [totMeanName, totErrName])
1753class SigmaDiaTotFluxConfig(DiaObjectCalculationPluginConfig):
1754 pass
1757@register("ap_sigmaTotFlux")
1758class SigmaDiaTotFlux(DiaObjectCalculationPlugin):
1759 """Compute scatter of diaSource fluxes.
1760 """
1762 ConfigClass = SigmaDiaPsfFluxConfig
1763 # Output columns are created upon instantiation of the class.
1764 outputCols = ["scienceFluxSigma"]
1765 plugType = "multi"
1766 needsFilter = True
1768 @classmethod
1769 def getExecutionOrder(cls):
1770 return cls.DEFAULT_CATALOGCALCULATION
1772 def calculate(self,
1773 diaObjects,
1774 diaSources,
1775 filterDiaSources,
1776 band,
1777 **kwargs):
1778 """Compute the sigma fluxes of the point source flux measured on the
1779 calibrated image.
1781 Parameters
1782 ----------
1783 diaObject : `dict`
1784 Summary object to store values in.
1785 diaSources : `pandas.DataFrame`
1786 DataFrame representing all diaSources associated with this
1787 diaObject.
1788 filterDiaSources : `pandas.DataFrame`
1789 DataFrame representing diaSources associated with this
1790 diaObject that are observed in the band pass ``band``.
1791 band : `str`
1792 Simple, string name of the filter for the flux being calculated.
1793 **kwargs
1794 Any additional keyword arguments that may be passed to the plugin.
1795 """
1796 # Set "delta degrees of freedom (ddf)" to 1 to calculate the unbiased
1797 # estimator of scatter (i.e. 'N - 1' instead of 'N').
1798 column = "{}_scienceFluxSigma".format(band)
1799 typeSafePandasAssignment(diaObjects, filterDiaSources.scienceFlux.std(), [column])