LSST Applications g063fba187b+fee0456c91,g0f08755f38+ea96e5a5a3,g1653933729+a8ce1bb630,g168dd56ebc+a8ce1bb630,g1a2382251a+90257ff92a,g20f6ffc8e0+ea96e5a5a3,g217e2c1bcf+937a289c59,g28da252d5a+daa7da44eb,g2bbee38e9b+253935c60e,g2bc492864f+253935c60e,g3156d2b45e+6e55a43351,g32e5bea42b+31359a2a7a,g347aa1857d+253935c60e,g35bb328faa+a8ce1bb630,g3a166c0a6a+253935c60e,g3b1af351f3+a8ce1bb630,g3e281a1b8c+c5dd892a6c,g414038480c+416496e02f,g41af890bb2+afe91b1188,g599934f4f4+0db33f7991,g7af13505b9+e36de7bce6,g80478fca09+da231ba887,g82479be7b0+a4516e59e3,g858d7b2824+ea96e5a5a3,g89c8672015+f4add4ffd5,g9125e01d80+a8ce1bb630,ga5288a1d22+bc6ab8dfbd,gb58c049af0+d64f4d3760,gc28159a63d+253935c60e,gcab2d0539d+3f2b72788c,gcf0d15dbbd+4ea9c45075,gda6a2b7d83+4ea9c45075,gdaeeff99f8+1711a396fd,ge79ae78c31+253935c60e,gef2f8181fd+3031e3cf99,gf0baf85859+c1f95f4921,gfa517265be+ea96e5a5a3,gfa999e8aa5+17cd334064,w.2024.50
LSST Data Management Base Package
Loading...
Searching...
No Matches
diaCalculationPlugins.py
Go to the documentation of this file.
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/>.
21
22"""Plugins for use in DiaSource summary statistics.
23
24Output columns must be
25as defined in the schema of the Apdb both in name and units.
26"""
27
28import functools
29import warnings
30
31from astropy.stats import median_absolute_deviation
32import numpy as np
33import pandas as pd
34from scipy.optimize import lsq_linear
35
36import lsst.geom as geom
37import lsst.pex.config as pexConfig
38import lsst.sphgeom as sphgeom
39
40from .diaCalculation import (
41 DiaObjectCalculationPluginConfig,
42 DiaObjectCalculationPlugin)
43from .pluginRegistry import register
44
45__all__ = ("MeanDiaPositionConfig", "MeanDiaPosition",
46 "HTMIndexDiaPosition", "HTMIndexDiaPositionConfig",
47 "NumDiaSourcesDiaPlugin", "NumDiaSourcesDiaPluginConfig",
48 "SimpleSourceFlagDiaPlugin", "SimpleSourceFlagDiaPluginConfig",
49 "WeightedMeanDiaPsfFluxConfig", "WeightedMeanDiaPsfFlux",
50 "PercentileDiaPsfFlux", "PercentileDiaPsfFluxConfig",
51 "SigmaDiaPsfFlux", "SigmaDiaPsfFluxConfig",
52 "Chi2DiaPsfFlux", "Chi2DiaPsfFluxConfig",
53 "MadDiaPsfFlux", "MadDiaPsfFluxConfig",
54 "SkewDiaPsfFlux", "SkewDiaPsfFluxConfig",
55 "MinMaxDiaPsfFlux", "MinMaxDiaPsfFluxConfig",
56 "MaxSlopeDiaPsfFlux", "MaxSlopeDiaPsfFluxConfig",
57 "ErrMeanDiaPsfFlux", "ErrMeanDiaPsfFluxConfig",
58 "LinearFitDiaPsfFlux", "LinearFitDiaPsfFluxConfig",
59 "StetsonJDiaPsfFlux", "StetsonJDiaPsfFluxConfig",
60 "WeightedMeanDiaTotFlux", "WeightedMeanDiaTotFluxConfig",
61 "SigmaDiaTotFlux", "SigmaDiaTotFluxConfig")
62
63
64def catchWarnings(_func=None, *, warns=[]):
65 """Decorator for generically catching numpy warnings.
66 """
67 def decoratorCatchWarnings(func):
68 @functools.wraps(func)
69 def wrapperCatchWarnings(*args, **kwargs):
70 with warnings.catch_warnings():
71 for val in warns:
72 warnings.filterwarnings("ignore", val)
73 return func(*args, **kwargs)
74 return wrapperCatchWarnings
75
76 if _func is None:
77 return decoratorCatchWarnings
78 else:
79 return decoratorCatchWarnings(_func)
80
81
83 pass
84
85
86@register("ap_meanPosition")
88 """Compute the mean position of a DiaObject given a set of DiaSources.
89 """
90
91 ConfigClass = MeanDiaPositionConfig
92
93 plugType = 'multi'
94
95 outputCols = ["ra", "dec", "radecMjdTai"]
96 needsFilter = False
97
98 @classmethod
101
102 def calculate(self, diaObjects, diaSources, **kwargs):
103 """Compute the mean ra/dec position of the diaObject given the
104 diaSource locations.
105
106 Parameters
107 ----------
108 diaObjects : `pandas.DataFrame`
109 Summary objects to store values in.
110 diaSources : `pandas.DataFrame` or `pandas.DataFrameGroupBy`
111 Catalog of DiaSources summarized by this DiaObject.
112 **kwargs
113 Any additional keyword arguments that may be passed to the plugin.
114 """
115 for outCol in self.outputColsoutputCols:
116 if outCol not in diaObjects.columns:
117 diaObjects[outCol] = np.nan
118
119 def _computeMeanPos(df):
120 aveCoord = geom.averageSpherePoint(
121 list(geom.SpherePoint(src["ra"], src["dec"], geom.degrees)
122 for idx, src in df.iterrows()))
123 ra = aveCoord.getRa().asDegrees()
124 dec = aveCoord.getDec().asDegrees()
125 if np.isnan(ra) or np.isnan(dec):
126 radecMjdTai = np.nan
127 else:
128 radecMjdTai = df["midpointMjdTai"].max()
129
130 return pd.Series({"ra": aveCoord.getRa().asDegrees(),
131 "dec": aveCoord.getDec().asDegrees(),
132 "radecMjdTai": radecMjdTai})
133
134 ans = diaSources.apply(_computeMeanPos)
135 diaObjects.loc[:, ["ra", "dec", "radecMjdTai"]] = ans
136
137
139
140 htmLevel = pexConfig.Field(
141 dtype=int,
142 doc="Level of the HTM pixelization.",
143 default=20,
144 )
145
146
147@register("ap_HTMIndex")
149 """Compute the mean position of a DiaObject given a set of DiaSources.
150
151 Notes
152 -----
153 This plugin was implemented to satisfy requirements of old APDB interface
154 which required ``pixelId`` column in DiaObject with HTM20 index. APDB
155 interface had migrated to not need that information, but we keep this
156 plugin in case it may be useful for something else.
157 """
158 ConfigClass = HTMIndexDiaPositionConfig
159
160 plugType = 'single'
161
162 inputCols = ["ra", "dec"]
163 outputCols = ["pixelId"]
164 needsFilter = False
165
166 def __init__(self, config, name, metadata):
167 DiaObjectCalculationPlugin.__init__(self, config, name, metadata)
169
170 @classmethod
172 return cls.FLUX_MOMENTS_CALCULATED
173
174 def calculate(self, diaObjects, diaObjectId, **kwargs):
175 """Compute the mean position of a DiaObject given a set of DiaSources
176
177 Parameters
178 ----------
179 diaObjects : `pandas.dataFrame`
180 Summary objects to store values in and read ra/dec from.
181 diaObjectId : `int`
182 Id of the diaObject to update.
183 **kwargs
184 Any additional keyword arguments that may be passed to the plugin.
185 """
186 sphPoint = geom.SpherePoint(
187 diaObjects.at[diaObjectId, "ra"] * geom.degrees,
188 diaObjects.at[diaObjectId, "dec"] * geom.degrees)
189 diaObjects.at[diaObjectId, "pixelId"] = self.pixelator.index(
190 sphPoint.getVector())
191
192
194 pass
195
196
197@register("ap_nDiaSources")
199 """Compute the total number of DiaSources associated with this DiaObject.
200 """
201
202 ConfigClass = NumDiaSourcesDiaPluginConfig
203 outputCols = ["nDiaSources"]
204 plugType = "multi"
205 needsFilter = False
206
207 @classmethod
210
211 def calculate(self, diaObjects, diaSources, **kwargs):
212 """Compute the total number of DiaSources associated with this DiaObject.
213
214 Parameters
215 ----------
216 diaObject : `dict`
217 Summary object to store values in and read ra/dec from.
218 **kwargs
219 Any additional keyword arguments that may be passed to the plugin.
220 """
221 dtype = diaObjects["nDiaSources"].dtype
222 diaObjects.loc[:, "nDiaSources"] = diaSources.diaObjectId.count().astype(dtype)
223
224
226 pass
227
228
229@register("ap_diaObjectFlag")
231 """Find if any DiaSource is flagged.
232
233 Set the DiaObject flag if any DiaSource is flagged.
234 """
235
236 ConfigClass = NumDiaSourcesDiaPluginConfig
237 outputCols = ["flags"]
238 plugType = "multi"
239 needsFilter = False
240
241 @classmethod
244
245 def calculate(self, diaObjects, diaSources, **kwargs):
246 """Find if any DiaSource is flagged.
247
248 Set the DiaObject flag if any DiaSource is flagged.
249
250 Parameters
251 ----------
252 diaObject : `dict`
253 Summary object to store values in and read ra/dec from.
254 **kwargs
255 Any additional keyword arguments that may be passed to the plugin.
256 """
257 dtype = diaObjects["flags"].dtype
258 diaObjects.loc[:, "flags"] = diaSources.flags.any().astype(dtype)
259
260
262 pass
263
264
265@register("ap_meanFlux")
267 """Compute the weighted mean and mean error on the point source fluxes
268 of the DiaSource measured on the difference image.
269
270 Additionally store number of usable data points.
271 """
272
273 ConfigClass = WeightedMeanDiaPsfFluxConfig
274 outputCols = ["psfFluxMean", "psfFluxMeanErr", "psfFluxNdata"]
275 plugType = "multi"
276 needsFilter = True
277
278 @classmethod
281
282 @catchWarnings(warns=["invalid value encountered",
283 "divide by zero"])
284 def calculate(self,
285 diaObjects,
286 diaSources,
287 filterDiaSources,
288 band,
289 **kwargs):
290 """Compute the weighted mean and mean error of the point source flux.
291
292 Parameters
293 ----------
294 diaObject : `dict`
295 Summary object to store values in.
296 diaSources : `pandas.DataFrame`
297 DataFrame representing all diaSources associated with this
298 diaObject.
299 filterDiaSources : `pandas.DataFrame`
300 DataFrame representing diaSources associated with this
301 diaObject that are observed in the band pass ``band``.
302 band : `str`
303 Simple, string name of the filter for the flux being calculated.
304 **kwargs
305 Any additional keyword arguments that may be passed to the plugin.
306 """
307 meanName = "{}_psfFluxMean".format(band)
308 errName = "{}_psfFluxMeanErr".format(band)
309 nDataName = "{}_psfFluxNdata".format(band)
310 if meanName not in diaObjects.columns:
311 diaObjects[meanName] = np.nan
312 if errName not in diaObjects.columns:
313 diaObjects[errName] = np.nan
314 if nDataName not in diaObjects.columns:
315 diaObjects[nDataName] = 0
316
317 def _weightedMean(df):
318 tmpDf = df[~np.logical_or(np.isnan(df["psfFlux"]),
319 np.isnan(df["psfFluxErr"]))]
320 tot_weight = np.nansum(1 / tmpDf["psfFluxErr"] ** 2)
321 fluxMean = np.nansum(tmpDf["psfFlux"]
322 / tmpDf["psfFluxErr"] ** 2)
323 fluxMean /= tot_weight
324 if tot_weight > 0:
325 fluxMeanErr = np.sqrt(1 / tot_weight)
326 else:
327 fluxMeanErr = np.nan
328 nFluxData = len(tmpDf)
329
330 return pd.Series({meanName: fluxMean,
331 errName: fluxMeanErr,
332 nDataName: nFluxData},
333 dtype="object")
334 df = filterDiaSources.apply(_weightedMean).astype(diaObjects.dtypes[[meanName, errName, nDataName]])
335
336 diaObjects.loc[:, [meanName, errName, nDataName]] = df
337
338
340 percentiles = pexConfig.ListField(
341 dtype=int,
342 default=[5, 25, 50, 75, 95],
343 doc="Percentiles to calculate to compute values for. Should be "
344 "integer values."
345 )
346
347
348@register("ap_percentileFlux")
350 """Compute percentiles of diaSource fluxes.
351 """
352
353 ConfigClass = PercentileDiaPsfFluxConfig
354 # Output columns are created upon instantiation of the class.
355 outputCols = []
356 plugType = "multi"
357 needsFilter = True
358
359 def __init__(self, config, name, metadata, **kwargs):
360 DiaObjectCalculationPlugin.__init__(self,
361 config,
362 name,
363 metadata,
364 **kwargs)
365 self.outputColsoutputColsoutputCols = ["psfFluxPercentile{:02d}".format(percent)
366 for percent in self.config.percentiles]
367
368 @classmethod
371
372 @catchWarnings(warns=["All-NaN slice encountered"])
373 def calculate(self,
374 diaObjects,
375 diaSources,
376 filterDiaSources,
377 band,
378 **kwargs):
379 """Compute the percentile fluxes of the point source flux.
380
381 Parameters
382 ----------
383 diaObject : `dict`
384 Summary object to store values in.
385 diaSources : `pandas.DataFrame`
386 DataFrame representing all diaSources associated with this
387 diaObject.
388 filterDiaSources : `pandas.DataFrame`
389 DataFrame representing diaSources associated with this
390 diaObject that are observed in the band pass ``band``.
391 band : `str`
392 Simple, string name of the filter for the flux being calculated.
393 **kwargs
394 Any additional keyword arguments that may be passed to the plugin.
395 """
396 pTileNames = []
397 dtype = None
398 for tilePercent in self.config.percentiles:
399 pTileName = "{}_psfFluxPercentile{:02d}".format(band,
400 tilePercent)
401 pTileNames.append(pTileName)
402 if pTileName not in diaObjects.columns:
403 diaObjects[pTileName] = np.nan
404 elif dtype is None:
405 dtype = diaObjects[pTileName].dtype
406
407 def _fluxPercentiles(df):
408 pTiles = np.nanpercentile(df["psfFlux"], self.config.percentiles)
409 return pd.Series(
410 dict((tileName, pTile)
411 for tileName, pTile in zip(pTileNames, pTiles)))
412 if dtype is None:
413 diaObjects.loc[:, pTileNames] = filterDiaSources.apply(_fluxPercentiles)
414 else:
415 diaObjects.loc[:, pTileNames] = filterDiaSources.apply(_fluxPercentiles).astype(dtype)
416
417
419 pass
420
421
422@register("ap_sigmaFlux")
424 """Compute scatter of diaSource fluxes.
425 """
426
427 ConfigClass = SigmaDiaPsfFluxConfig
428 # Output columns are created upon instantiation of the class.
429 outputCols = ["psfFluxSigma"]
430 plugType = "multi"
431 needsFilter = True
432
433 @classmethod
436
437 def calculate(self,
438 diaObjects,
439 diaSources,
440 filterDiaSources,
441 band,
442 **kwargs):
443 """Compute the sigma fluxes of the point source flux.
444
445 Parameters
446 ----------
447 diaObject : `dict`
448 Summary object to store values in.
449 diaSources : `pandas.DataFrame`
450 DataFrame representing all diaSources associated with this
451 diaObject.
452 filterDiaSources : `pandas.DataFrame`
453 DataFrame representing diaSources associated with this
454 diaObject that are observed in the band pass ``band``.
455 band : `str`
456 Simple, string name of the filter for the flux being calculated.
457 **kwargs
458 Any additional keyword arguments that may be passed to the plugin.
459 """
460 # Set "delta degrees of freedom (ddf)" to 1 to calculate the unbiased
461 # estimator of scatter (i.e. 'N - 1' instead of 'N').
462 column = "{}_psfFluxSigma".format(band)
463 if column in diaObjects:
464 dtype = diaObjects[column].dtype
465 diaObjects.loc[:, column] = filterDiaSources.psfFlux.std().astype(dtype)
466 else:
467 diaObjects.loc[:, column] = filterDiaSources.psfFlux.std()
468
469
471 pass
472
473
474@register("ap_chi2Flux")
476 """Compute chi2 of diaSource fluxes.
477 """
478
479 ConfigClass = Chi2DiaPsfFluxConfig
480
481 # Required input Cols
482 inputCols = ["psfFluxMean"]
483 # Output columns are created upon instantiation of the class.
484 outputCols = ["psfFluxChi2"]
485 plugType = "multi"
486 needsFilter = True
487
488 @classmethod
490 return cls.FLUX_MOMENTS_CALCULATED
491
492 @catchWarnings(warns=["All-NaN slice encountered"])
493 def calculate(self,
494 diaObjects,
495 diaSources,
496 filterDiaSources,
497 band,
498 **kwargs):
499 """Compute the chi2 of the point source fluxes.
500
501 Parameters
502 ----------
503 diaObject : `dict`
504 Summary object to store values in.
505 diaSources : `pandas.DataFrame`
506 DataFrame representing all diaSources associated with this
507 diaObject.
508 filterDiaSources : `pandas.DataFrame`
509 DataFrame representing diaSources associated with this
510 diaObject that are observed in the band pass ``band``.
511 band : `str`
512 Simple, string name of the filter for the flux being calculated.
513 **kwargs
514 Any additional keyword arguments that may be passed to the plugin.
515 """
516 meanName = "{}_psfFluxMean".format(band)
517 column = "{}_psfFluxChi2".format(band)
518
519 def _chi2(df):
520 delta = (df["psfFlux"]
521 - diaObjects.at[df.diaObjectId.iat[0], meanName])
522 return np.nansum((delta / df["psfFluxErr"]) ** 2)
523
524 if column in diaObjects:
525 dtype = diaObjects[column].dtype
526 diaObjects.loc[:, column] = filterDiaSources.apply(_chi2).astype(dtype)
527 else:
528 diaObjects.loc[:, column] = filterDiaSources.apply(_chi2)
529
530
532 pass
533
534
535@register("ap_madFlux")
537 """Compute median absolute deviation of diaSource fluxes.
538 """
539
540 ConfigClass = MadDiaPsfFluxConfig
541
542 # Required input Cols
543 # Output columns are created upon instantiation of the class.
544 outputCols = ["psfFluxMAD"]
545 plugType = "multi"
546 needsFilter = True
547
548 @classmethod
551
552 @catchWarnings(warns=["All-NaN slice encountered"])
553 def calculate(self,
554 diaObjects,
555 diaSources,
556 filterDiaSources,
557 band,
558 **kwargs):
559 """Compute the median absolute deviation of the point source fluxes.
560
561 Parameters
562 ----------
563 diaObject : `dict`
564 Summary object to store values in.
565 diaSources : `pandas.DataFrame`
566 DataFrame representing all diaSources associated with this
567 diaObject.
568 filterDiaSources : `pandas.DataFrame`
569 DataFrame representing diaSources associated with this
570 diaObject that are observed in the band pass ``band``.
571 band : `str`
572 Simple, string name of the filter for the flux being calculated.
573 **kwargs
574 Any additional keyword arguments that may be passed to the plugin.
575 """
576 column = "{}_psfFluxMAD".format(band)
577 if column in diaObjects:
578 dtype = diaObjects[column].dtype
579 diaObjects.loc[:, column] = \
580 filterDiaSources.psfFlux.apply(median_absolute_deviation,
581 ignore_nan=True).astype(dtype)
582 else:
583 diaObjects.loc[:, column] = \
584 filterDiaSources.psfFlux.apply(median_absolute_deviation,
585 ignore_nan=True)
586
587
589 pass
590
591
592@register("ap_skewFlux")
594 """Compute the skew of diaSource fluxes.
595 """
596
597 ConfigClass = SkewDiaPsfFluxConfig
598
599 # Required input Cols
600 # Output columns are created upon instantiation of the class.
601 outputCols = ["psfFluxSkew"]
602 plugType = "multi"
603 needsFilter = True
604
605 @classmethod
608
609 def calculate(self,
610 diaObjects,
611 diaSources,
612 filterDiaSources,
613 band,
614 **kwargs):
615 """Compute the skew of the point source fluxes.
616
617 Parameters
618 ----------
619 diaObject : `dict`
620 Summary object to store values in.
621 diaSources : `pandas.DataFrame`
622 DataFrame representing all diaSources associated with this
623 diaObject.
624 filterDiaSources : `pandas.DataFrame`
625 DataFrame representing diaSources associated with this
626 diaObject that are observed in the band pass ``band``.
627 band : `str`
628 Simple, string name of the filter for the flux being calculated.
629 **kwargs
630 Any additional keyword arguments that may be passed to the plugin.
631 """
632 column = "{}_psfFluxSkew".format(band)
633 if column in diaObjects:
634 dtype = diaObjects[column].dtype
635 diaObjects.loc[:, column] = filterDiaSources.psfFlux.skew().astype(dtype)
636 else:
637 diaObjects.loc[:, column] = filterDiaSources.psfFlux.skew()
638
639
641 pass
642
643
644@register("ap_minMaxFlux")
646 """Compute min/max of diaSource fluxes.
647 """
648
649 ConfigClass = MinMaxDiaPsfFluxConfig
650
651 # Required input Cols
652 # Output columns are created upon instantiation of the class.
653 outputCols = ["psfFluxMin", "psfFluxMax"]
654 plugType = "multi"
655 needsFilter = True
656
657 @classmethod
660
661 def calculate(self,
662 diaObjects,
663 diaSources,
664 filterDiaSources,
665 band,
666 **kwargs):
667 """Compute min/max of the point source fluxes.
668
669 Parameters
670 ----------
671 diaObject : `dict`
672 Summary object to store values in.
673 diaSources : `pandas.DataFrame`
674 DataFrame representing all diaSources associated with this
675 diaObject.
676 filterDiaSources : `pandas.DataFrame`
677 DataFrame representing diaSources associated with this
678 diaObject that are observed in the band pass ``band``.
679 band : `str`
680 Simple, string name of the filter for the flux being calculated.
681 **kwargs
682 Any additional keyword arguments that may be passed to the plugin.
683 """
684 minName = "{}_psfFluxMin".format(band)
685 if minName not in diaObjects.columns:
686 diaObjects[minName] = np.nan
687 maxName = "{}_psfFluxMax".format(band)
688 if maxName not in diaObjects.columns:
689 diaObjects[maxName] = np.nan
690
691 dtype = diaObjects[minName].dtype
692 diaObjects.loc[:, minName] = filterDiaSources.psfFlux.min().astype(dtype)
693 diaObjects.loc[:, maxName] = filterDiaSources.psfFlux.max().astype(dtype)
694
695
697 pass
698
699
700@register("ap_maxSlopeFlux")
702 """Compute the maximum ratio time ordered deltaFlux / deltaTime.
703 """
704
705 ConfigClass = MinMaxDiaPsfFluxConfig
706
707 # Required input Cols
708 # Output columns are created upon instantiation of the class.
709 outputCols = ["psfFluxMaxSlope"]
710 plugType = "multi"
711 needsFilter = True
712
713 @classmethod
716
717 def calculate(self,
718 diaObjects,
719 diaSources,
720 filterDiaSources,
721 band,
722 **kwargs):
723 """Compute the maximum ratio time ordered deltaFlux / deltaTime.
724
725 Parameters
726 ----------
727 diaObject : `dict`
728 Summary object to store values in.
729 diaSources : `pandas.DataFrame`
730 DataFrame representing all diaSources associated with this
731 diaObject.
732 filterDiaSources : `pandas.DataFrame`
733 DataFrame representing diaSources associated with this
734 diaObject that are observed in the band pass ``band``.
735 band : `str`
736 Simple, string name of the filter for the flux being calculated.
737 **kwargs
738 Any additional keyword arguments that may be passed to the plugin.
739 """
740
741 def _maxSlope(df):
742 tmpDf = df[~np.logical_or(np.isnan(df["psfFlux"]),
743 np.isnan(df["midpointMjdTai"]))]
744 if len(tmpDf) < 2:
745 return np.nan
746 times = tmpDf["midpointMjdTai"].to_numpy()
747 timeArgs = times.argsort()
748 times = times[timeArgs]
749 fluxes = tmpDf["psfFlux"].to_numpy()[timeArgs]
750 return (np.diff(fluxes) / np.diff(times)).max()
751
752 column = "{}_psfFluxMaxSlope".format(band)
753 if column in diaObjects:
754 dtype = diaObjects[column].dtype
755 diaObjects.loc[:, column] = filterDiaSources.apply(_maxSlope).astype(dtype)
756 else:
757 diaObjects.loc[:, column] = filterDiaSources.apply(_maxSlope)
758
759
761 pass
762
763
764@register("ap_meanErrFlux")
766 """Compute the mean of the dia source errors.
767 """
768
769 ConfigClass = ErrMeanDiaPsfFluxConfig
770
771 # Required input Cols
772 # Output columns are created upon instantiation of the class.
773 outputCols = ["psfFluxErrMean"]
774 plugType = "multi"
775 needsFilter = True
776
777 @classmethod
780
781 def calculate(self,
782 diaObjects,
783 diaSources,
784 filterDiaSources,
785 band,
786 **kwargs):
787 """Compute the mean of the dia source errors.
788
789 Parameters
790 ----------
791 diaObject : `dict`
792 Summary object to store values in.
793 diaSources : `pandas.DataFrame`
794 DataFrame representing all diaSources associated with this
795 diaObject.
796 filterDiaSources : `pandas.DataFrame`
797 DataFrame representing diaSources associated with this
798 diaObject that are observed in the band pass ``band``.
799 band : `str`
800 Simple, string name of the filter for the flux being calculated.
801 **kwargs
802 Any additional keyword arguments that may be passed to the plugin.
803 """
804 column = "{}_psfFluxErrMean".format(band)
805 if column in diaObjects:
806 dtype = diaObjects[column].dtype
807 diaObjects.loc[:, column] = filterDiaSources.psfFluxErr.mean().astype(dtype)
808 else:
809 diaObjects.loc[:, column] = filterDiaSources.psfFluxErr.mean()
810
811
813 pass
814
815
816@register("ap_linearFit")
818 """Compute fit a linear model to flux vs time.
819 """
820
821 ConfigClass = LinearFitDiaPsfFluxConfig
822
823 # Required input Cols
824 # Output columns are created upon instantiation of the class.
825 outputCols = ["psfFluxLinearSlope", "psfFluxLinearIntercept"]
826 plugType = "multi"
827 needsFilter = True
828
829 @classmethod
832
833 def calculate(self,
834 diaObjects,
835 diaSources,
836 filterDiaSources,
837 band,
838 **kwargs):
839 """Compute fit a linear model to flux vs time.
840
841 Parameters
842 ----------
843 diaObject : `dict`
844 Summary object to store values in.
845 diaSources : `pandas.DataFrame`
846 DataFrame representing all diaSources associated with this
847 diaObject.
848 filterDiaSources : `pandas.DataFrame`
849 DataFrame representing diaSources associated with this
850 diaObject that are observed in the band pass ``band``.
851 band : `str`
852 Simple, string name of the filter for the flux being calculated.
853 **kwargs
854 Any additional keyword arguments that may be passed to the plugin.
855 """
856
857 mName = "{}_psfFluxLinearSlope".format(band)
858 if mName not in diaObjects.columns:
859 diaObjects[mName] = np.nan
860 bName = "{}_psfFluxLinearIntercept".format(band)
861 if bName not in diaObjects.columns:
862 diaObjects[bName] = np.nan
863 dtype = diaObjects[mName].dtype
864
865 def _linearFit(df):
866 tmpDf = df[~np.logical_or(
867 np.isnan(df["psfFlux"]),
868 np.logical_or(np.isnan(df["psfFluxErr"]),
869 np.isnan(df["midpointMjdTai"])))]
870 if len(tmpDf) < 2:
871 return pd.Series({mName: np.nan, bName: np.nan})
872 fluxes = tmpDf["psfFlux"].to_numpy()
873 errors = tmpDf["psfFluxErr"].to_numpy()
874 times = tmpDf["midpointMjdTai"].to_numpy()
875 A = np.array([times / errors, 1 / errors]).transpose()
876 m, b = lsq_linear(A, fluxes / errors).x
877 return pd.Series({mName: m, bName: b}, dtype=dtype)
878
879 diaObjects.loc[:, [mName, bName]] = filterDiaSources.apply(_linearFit)
880
881
883 pass
884
885
886@register("ap_stetsonJ")
888 """Compute the StetsonJ statistic on the DIA point source fluxes.
889 """
890
891 ConfigClass = LinearFitDiaPsfFluxConfig
892
893 # Required input Cols
894 inputCols = ["psfFluxMean"]
895 # Output columns are created upon instantiation of the class.
896 outputCols = ["psfFluxStetsonJ"]
897 plugType = "multi"
898 needsFilter = True
899
900 @classmethod
902 return cls.FLUX_MOMENTS_CALCULATED
903
904 def calculate(self,
905 diaObjects,
906 diaSources,
907 filterDiaSources,
908 band,
909 **kwargs):
910 """Compute the StetsonJ statistic on the DIA point source fluxes.
911
912 Parameters
913 ----------
914 diaObject : `dict`
915 Summary object to store values in.
916 diaSources : `pandas.DataFrame`
917 DataFrame representing all diaSources associated with this
918 diaObject.
919 filterDiaSources : `pandas.DataFrame`
920 DataFrame representing diaSources associated with this
921 diaObject that are observed in the band pass ``band``.
922 band : `str`
923 Simple, string name of the filter for the flux being calculated.
924 **kwargs
925 Any additional keyword arguments that may be passed to the plugin.
926 """
927 meanName = "{}_psfFluxMean".format(band)
928
929 def _stetsonJ(df):
930 tmpDf = df[~np.logical_or(np.isnan(df["psfFlux"]),
931 np.isnan(df["psfFluxErr"]))]
932 if len(tmpDf) < 2:
933 return np.nan
934 fluxes = tmpDf["psfFlux"].to_numpy()
935 errors = tmpDf["psfFluxErr"].to_numpy()
936
937 return self._stetson_J(
938 fluxes,
939 errors,
940 diaObjects.at[tmpDf.diaObjectId.iat[0], meanName])
941
942 column = "{}_psfFluxStetsonJ".format(band)
943 if column in diaObjects:
944 dtype = diaObjects[column].dtype
945 diaObjects.loc[:, column] = filterDiaSources.apply(_stetsonJ).astype(dtype)
946 else:
947 diaObjects.loc[:, column] = filterDiaSources.apply(_stetsonJ)
948
949 def _stetson_J(self, fluxes, errors, mean=None):
950 """Compute the single band stetsonJ statistic.
951
952 Parameters
953 ----------
954 fluxes : `numpy.ndarray` (N,)
955 Calibrated lightcurve flux values.
956 errors : `numpy.ndarray` (N,)
957 Errors on the calibrated lightcurve fluxes.
958 mean : `float`
959 Starting mean from previous plugin.
960
961 Returns
962 -------
963 stetsonJ : `float`
964 stetsonJ statistic for the input fluxes and errors.
965
966 References
967 ----------
968 .. [1] Stetson, P. B., "On the Automatic Determination of Light-Curve
969 Parameters for Cepheid Variables", PASP, 108, 851S, 1996
970 """
971 n_points = len(fluxes)
972 flux_mean = self._stetson_mean(fluxes, errors, mean)
973 delta_val = (
974 np.sqrt(n_points / (n_points - 1)) * (fluxes - flux_mean) / errors)
975 p_k = delta_val ** 2 - 1
976
977 return np.mean(np.sign(p_k) * np.sqrt(np.fabs(p_k)))
978
980 values,
981 errors,
982 mean=None,
983 alpha=2.,
984 beta=2.,
985 n_iter=20,
986 tol=1e-6):
987 """Compute the stetson mean of the fluxes which down-weights outliers.
988
989 Weighted biased on an error weighted difference scaled by a constant
990 (1/``a``) and raised to the power beta. Higher betas more harshly
991 penalize outliers and ``a`` sets the number of sigma where a weighted
992 difference of 1 occurs.
993
994 Parameters
995 ----------
996 values : `numpy.dnarray`, (N,)
997 Input values to compute the mean of.
998 errors : `numpy.ndarray`, (N,)
999 Errors on the input values.
1000 mean : `float`
1001 Starting mean value or None.
1002 alpha : `float`
1003 Scalar down-weighting of the fractional difference. lower->more
1004 clipping. (Default value is 2.)
1005 beta : `float`
1006 Power law slope of the used to down-weight outliers. higher->more
1007 clipping. (Default value is 2.)
1008 n_iter : `int`
1009 Number of iterations of clipping.
1010 tol : `float`
1011 Fractional and absolute tolerance goal on the change in the mean
1012 before exiting early. (Default value is 1e-6)
1013
1014 Returns
1015 -------
1016 mean : `float`
1017 Weighted stetson mean result.
1018
1019 References
1020 ----------
1021 .. [1] Stetson, P. B., "On the Automatic Determination of Light-Curve
1022 Parameters for Cepheid Variables", PASP, 108, 851S, 1996
1023 """
1024 n_points = len(values)
1025 n_factor = np.sqrt(n_points / (n_points - 1))
1026 inv_var = 1 / errors ** 2
1027
1028 if mean is None:
1029 mean = np.average(values, weights=inv_var)
1030 for iter_idx in range(n_iter):
1031 chi = np.fabs(n_factor * (values - mean) / errors)
1032 tmp_mean = np.average(
1033 values,
1034 weights=inv_var / (1 + (chi / alpha) ** beta))
1035 diff = np.fabs(tmp_mean - mean)
1036 mean = tmp_mean
1037 if diff / mean < tol and diff < tol:
1038 break
1039 return mean
1040
1041
1043 pass
1044
1045
1046@register("ap_meanTotFlux")
1048 """Compute the weighted mean and mean error on the point source fluxes
1049 forced photometered at the DiaSource location in the calibrated image.
1050
1051 Additionally store number of usable data points.
1052 """
1053
1054 ConfigClass = WeightedMeanDiaPsfFluxConfig
1055 outputCols = ["scienceFluxMean", "scienceFluxMeanErr"]
1056 plugType = "multi"
1057 needsFilter = True
1058
1059 @classmethod
1062
1063 @catchWarnings(warns=["invalid value encountered",
1064 "divide by zero"])
1065 def calculate(self,
1066 diaObjects,
1067 diaSources,
1068 filterDiaSources,
1069 band,
1070 **kwargs):
1071 """Compute the weighted mean and mean error of the point source flux.
1072
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 totMeanName = "{}_scienceFluxMean".format(band)
1089 if totMeanName not in diaObjects.columns:
1090 diaObjects[totMeanName] = np.nan
1091 totErrName = "{}_scienceFluxMeanErr".format(band)
1092 if totErrName not in diaObjects.columns:
1093 diaObjects[totErrName] = np.nan
1094
1095 def _meanFlux(df):
1096 tmpDf = df[~np.logical_or(np.isnan(df["scienceFlux"]),
1097 np.isnan(df["scienceFluxErr"]))]
1098 tot_weight = np.nansum(1 / tmpDf["scienceFluxErr"] ** 2)
1099 fluxMean = np.nansum(tmpDf["scienceFlux"]
1100 / tmpDf["scienceFluxErr"] ** 2)
1101 fluxMean /= tot_weight
1102 fluxMeanErr = np.sqrt(1 / tot_weight)
1103
1104 return pd.Series({totMeanName: fluxMean,
1105 totErrName: fluxMeanErr})
1106
1107 df = filterDiaSources.apply(_meanFlux).astype(diaObjects.dtypes[[totMeanName, totErrName]])
1108 diaObjects.loc[:, [totMeanName, totErrName]] = df
1109
1110
1112 pass
1113
1114
1115@register("ap_sigmaTotFlux")
1117 """Compute scatter of diaSource fluxes.
1118 """
1119
1120 ConfigClass = SigmaDiaPsfFluxConfig
1121 # Output columns are created upon instantiation of the class.
1122 outputCols = ["scienceFluxSigma"]
1123 plugType = "multi"
1124 needsFilter = True
1125
1126 @classmethod
1129
1130 def calculate(self,
1131 diaObjects,
1132 diaSources,
1133 filterDiaSources,
1134 band,
1135 **kwargs):
1136 """Compute the sigma fluxes of the point source flux measured on the
1137 calibrated image.
1138
1139 Parameters
1140 ----------
1141 diaObject : `dict`
1142 Summary object to store values in.
1143 diaSources : `pandas.DataFrame`
1144 DataFrame representing all diaSources associated with this
1145 diaObject.
1146 filterDiaSources : `pandas.DataFrame`
1147 DataFrame representing diaSources associated with this
1148 diaObject that are observed in the band pass ``band``.
1149 band : `str`
1150 Simple, string name of the filter for the flux being calculated.
1151 **kwargs
1152 Any additional keyword arguments that may be passed to the plugin.
1153 """
1154 # Set "delta degrees of freedom (ddf)" to 1 to calculate the unbiased
1155 # estimator of scatter (i.e. 'N - 1' instead of 'N').
1156 column = "{}_scienceFluxSigma".format(band)
1157 if column in diaObjects:
1158 dtype = diaObjects[column].dtype
1159 diaObjects.loc[:, column] = filterDiaSources.scienceFlux.std().astype(dtype)
1160 else:
1161
1162 diaObjects.loc[:, column] = filterDiaSources.scienceFlux.std()
int max
Point in an unspecified spherical coordinate system.
Definition SpherePoint.h:57
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, filterDiaSources, band, **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, filterDiaSources, band, **kwargs)
calculate(self, diaObjects, diaSources, filterDiaSources, band, **kwargs)
calculate(self, diaObjects, diaSources, filterDiaSources, band, **kwargs)
calculate(self, diaObjects, diaSources, filterDiaSources, band, **kwargs)
_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)
HtmPixelization provides HTM indexing of points and regions.
SpherePoint averageSpherePoint(std::vector< SpherePoint > const &coords)
Return the average of a list of coordinates.