Coverage for tests/test_diaCalculationPlugins.py: 14%
404 statements
« prev ^ index » next coverage.py v7.10.6, created at 2025-10-14 01:52 -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/>.
22from astropy.stats import median_absolute_deviation
23import numpy as np
24import pandas as pd
25from scipy.stats import skew
26import unittest
28from lsst.meas.base import (
29 MeanDiaPosition, MeanDiaPositionConfig,
30 HTMIndexDiaPosition, HTMIndexDiaPositionConfig,
31 NumDiaSourcesDiaPlugin, NumDiaSourcesDiaPluginConfig,
32 SimpleSourceFlagDiaPlugin, SimpleSourceFlagDiaPluginConfig,
33 WeightedMeanDiaPsfFlux, WeightedMeanDiaPsfFluxConfig,
34 PercentileDiaPsfFlux, PercentileDiaPsfFluxConfig,
35 SigmaDiaPsfFlux, SigmaDiaPsfFluxConfig,
36 Chi2DiaPsfFlux, Chi2DiaPsfFluxConfig,
37 MadDiaPsfFlux, MadDiaPsfFluxConfig,
38 SkewDiaPsfFlux, SkewDiaPsfFluxConfig,
39 MinMaxDiaPsfFlux, MinMaxDiaPsfFluxConfig,
40 MaxSlopeDiaPsfFlux, MaxSlopeDiaPsfFluxConfig,
41 ErrMeanDiaPsfFlux, ErrMeanDiaPsfFluxConfig,
42 LinearFitDiaPsfFlux, LinearFitDiaPsfFluxConfig,
43 StetsonJDiaPsfFlux, StetsonJDiaPsfFluxConfig,
44 WeightedMeanDiaTotFlux, WeightedMeanDiaTotFluxConfig,
45 SigmaDiaTotFlux, SigmaDiaTotFluxConfig,
46 LombScarglePeriodogram, LombScarglePeriodogramConfig,
47 LombScarglePeriodogramMulti, LombScarglePeriodogramMultiConfig)
48import lsst.utils.tests
51def run_single_plugin(diaObjectCat,
52 diaObjectId,
53 diaSourceCat,
54 band,
55 plugin):
56 """Wrapper for running single plugins.
58 Reproduces some of the behavior of `lsst.ap.association.DiaCalcuation.run`
60 Parameters
61 ----------
62 diaObjectCat : `pandas.DataFrame`
63 Input object catalog to store data into and read from.
64 diaSourcesCat : `pandas.DataFrame`
65 DiaSource catalog to read data from and groupby on.
66 fitlerName : `str`
67 String name of the filter to process.
68 plugin : `lsst.ap.association.DiaCalculationPlugin`
69 Plugin to run.
70 """
71 diaObjectCat.set_index("diaObjectId", inplace=True, drop=False)
72 diaSourceCat.set_index(
73 ["diaObjectId", "band", "diaSourceId"],
74 inplace=True,
75 drop=False)
77 objDiaSources = diaSourceCat.loc[diaObjectId]
78 updatingFilterDiaSources = diaSourceCat.loc[
79 (diaObjectId, band), :
80 ]
82 plugin.calculate(diaObjects=diaObjectCat,
83 diaObjectId=diaObjectId,
84 diaSources=objDiaSources,
85 filterDiaSources=updatingFilterDiaSources,
86 band=band)
89def run_multi_plugin(diaObjectCat, diaSourceCat, band, plugin):
90 """Wrapper for running multi plugins.
92 Reproduces some of the behavior of `lsst.ap.association.DiaCalcuation.run`
94 Parameters
95 ----------
96 diaObjectCat : `pandas.DataFrame`
97 Input object catalog to store data into and read from.
98 diaSourcesCat : `pandas.DataFrame`
99 DiaSource catalog to read data from and groupby on.
100 filterName : `str`
101 String name of the filter to process.
102 plugin : `lsst.ap.association.DiaCalculationPlugin`
103 Plugin to run.
104 """
105 diaObjectCat.set_index("diaObjectId", inplace=True, drop=False)
106 diaSourceCat.set_index(
107 ["diaObjectId", "band", "diaSourceId"],
108 inplace=True,
109 drop=False)
111 updatingFilterDiaSources = diaSourceCat.loc[
112 (slice(None), band), :
113 ]
115 diaSourcesGB = diaSourceCat.groupby(level=0)
116 filterDiaSourcesGB = updatingFilterDiaSources.groupby(level=0)
118 plugin.calculate(diaObjects=diaObjectCat,
119 diaSources=diaSourcesGB,
120 filterDiaSources=filterDiaSourcesGB,
121 band=band)
124def run_multiband_plugin(diaObjectCat, diaSourceCat, plugin):
125 """Wrapper for running multi plugins.
127 Reproduces some of the behavior of `lsst.ap.association.DiaCalcuation.run`
129 Parameters
130 ----------
131 diaObjectCat : `pandas.DataFrame`
132 Input object catalog to store data into and read from.
133 diaSourcesCat : `pandas.DataFrame`
134 DiaSource catalog to read data from and groupby on.
135 plugin : `lsst.ap.association.DiaCalculationPlugin`
136 Plugin to run.
137 """
138 diaObjectCat.set_index("diaObjectId", inplace=True, drop=False)
139 diaSourceCat.set_index(
140 ["diaObjectId", "band", "diaSourceId"],
141 inplace=True,
142 drop=False)
144 diaSourcesGB = diaSourceCat.groupby(level=0)
146 plugin.calculate(diaObjects=diaObjectCat,
147 diaSources=diaSourcesGB,
148 )
151def make_diaObject_table(objId, plugin, default_value=None, band=None):
152 """Create a minimal diaObject table with columns required for the plugin
154 Parameters
155 ----------
156 objId : `int`
157 The diaObjectId
158 plugin : `lsst.ap.association.DiaCalculationPlugin`
159 The plugin that will be run.
160 default_value : None, optional
161 Value to set new columns to.
162 band : `str`, optional
163 Band designation to append to the plugin columns.
165 Returns
166 -------
167 diaObjects : `pandas.DataFrame`
168 Output catalog with the required columns for the plugin.
169 """
170 diaObjects = {"diaObjectId": [objId]}
171 for col in plugin.outputCols:
172 if band is not None:
173 diaObjects[f"{band}_{col}"] = default_value
174 else:
175 diaObjects[col] = default_value
176 return pd.DataFrame(diaObjects)
179class TestMeanPosition(unittest.TestCase):
181 def testCalculate(self):
182 """Test mean position calculation.
183 """
184 n_sources = 10
185 objId = 0
187 plug = MeanDiaPosition(MeanDiaPositionConfig(),
188 "ap_meanPosition",
189 None)
191 # Test expected means in RA.
192 diaObjects = pd.DataFrame({"diaObjectId": [objId]})
193 diaSources = pd.DataFrame(data={"ra": np.linspace(-1, 1, n_sources),
194 "dec": np.zeros(n_sources),
195 "midpointMjdTai": np.linspace(0, n_sources, n_sources),
196 "diaObjectId": n_sources * [objId],
197 "band": n_sources * ["g"],
198 "diaSourceId": np.arange(n_sources,
199 dtype=int)})
200 run_multi_plugin(diaObjects, diaSources, "g", plug)
202 self.assertAlmostEqual(diaObjects.loc[objId, "ra"], 0.0)
203 self.assertAlmostEqual(diaObjects.loc[objId, "dec"], 0.0)
205 # Test expected means in DEC.
206 diaObjects = pd.DataFrame({"diaObjectId": [objId]})
207 diaSources = pd.DataFrame(data={"ra": np.zeros(n_sources),
208 "dec": np.linspace(-1, 1, n_sources),
209 "midpointMjdTai": np.linspace(0, n_sources, n_sources),
210 "diaObjectId": n_sources * [objId],
211 "band": n_sources * ["g"],
212 "diaSourceId": np.arange(n_sources,
213 dtype=int)})
214 run_multi_plugin(diaObjects, diaSources, "g", plug)
216 self.assertAlmostEqual(diaObjects.loc[objId, "ra"], 0.0)
217 self.assertAlmostEqual(diaObjects.loc[objId, "dec"], 0.0)
219 # Test failure mode RA is nan.
220 diaObjects = pd.DataFrame({"diaObjectId": [objId]})
221 diaSources = pd.DataFrame(data={"ra": np.full(n_sources, np.nan),
222 "dec": np.zeros(n_sources),
223 "midpointMjdTai": np.linspace(0, n_sources, n_sources),
224 "diaObjectId": n_sources * [objId],
225 "band": n_sources * ["g"],
226 "diaSourceId": np.arange(n_sources,
227 dtype=int)})
228 run_multi_plugin(diaObjects, diaSources, "g", plug)
230 self.assertTrue(np.isnan(diaObjects.loc[objId, "ra"]))
231 self.assertTrue(np.isnan(diaObjects.loc[objId, "dec"]))
233 # Test failure mode DEC is nan.
234 diaObjects = pd.DataFrame({"diaObjectId": [objId]})
235 diaSources = pd.DataFrame(data={"ra": np.zeros(n_sources),
236 "dec": np.full(n_sources, np.nan),
237 "midpointMjdTai": np.linspace(0, n_sources, n_sources),
238 "diaObjectId": n_sources * [objId],
239 "band": n_sources * ["g"],
240 "diaSourceId": np.arange(n_sources,
241 dtype=int)})
242 run_multi_plugin(diaObjects, diaSources, "g", plug)
244 self.assertTrue(np.isnan(diaObjects.loc[objId, "ra"]))
245 self.assertTrue(np.isnan(diaObjects.loc[objId, "dec"]))
248class TestHTMIndexPosition(unittest.TestCase):
250 def testCalculate(self):
251 """Test HTMPixel assignment calculation.
252 """
253 # Test expected pixelId at RA, DEC = 0
254 objId = 0
255 n_sources = 10
256 diaObjects = pd.DataFrame({"diaObjectId": [objId]})
257 diaObjects.loc[objId, "ra"] = 0.
258 diaObjects.loc[objId, "dec"] = 0.
259 diaSources = pd.DataFrame(
260 data={"diaObjectId": n_sources * [objId],
261 "band": n_sources * ["g"],
262 "diaSourceId": np.arange(n_sources, dtype=int)})
263 plug = HTMIndexDiaPosition(HTMIndexDiaPositionConfig(),
264 "ap_HTMIndex",
265 None)
267 run_single_plugin(diaObjectCat=diaObjects,
268 diaObjectId=objId,
269 diaSourceCat=diaSources,
270 band="g",
271 plugin=plug)
272 self.assertEqual(diaObjects.at[objId, "pixelId"],
273 17042430230528)
275 # Test expected pixelId at some value of RA and DEC.
276 diaObjects = pd.DataFrame({"diaObjectId": [objId]})
277 diaObjects.loc[objId, "ra"] = 45.37
278 diaObjects.loc[objId, "dec"] = 13.67
279 diaSources = pd.DataFrame(
280 data={"diaObjectId": n_sources * [objId],
281 "band": n_sources * ["g"],
282 "diaSourceId": np.arange(n_sources, dtype=int)})
283 run_single_plugin(diaObjectCat=diaObjects,
284 diaObjectId=objId,
285 diaSourceCat=diaSources,
286 band="g",
287 plugin=plug)
288 self.assertEqual(diaObjects.at[objId, "pixelId"],
289 17450571968473)
292class TestNDiaSourcesDiaPlugin(unittest.TestCase):
294 def testCalculate(self):
295 """Test that the number of DiaSources is correct.
296 """
298 for n_sources in [1, 8, 10]:
299 # Test expected number of sources per object.
300 objId = 0
301 diaSources = pd.DataFrame(
302 data={"diaObjectId": n_sources * [objId],
303 "band": n_sources * ["g"],
304 "diaSourceId": np.arange(n_sources, dtype=int)})
305 plug = NumDiaSourcesDiaPlugin(NumDiaSourcesDiaPluginConfig(),
306 "ap_nDiaSources",
307 None)
308 diaObjects = make_diaObject_table(objId, plug, default_value=int(-1))
309 run_multi_plugin(diaObjects, diaSources, "g", plug)
311 self.assertEqual(n_sources, diaObjects.at[objId, "nDiaSources"])
314class TestSimpleSourceFlagDiaPlugin(unittest.TestCase):
316 def testCalculate(self):
317 """Test that DiaObject flags are set.
318 """
319 objId = 0
320 n_sources = 10
322 # Test expected flags, no flags set.
323 diaSources = pd.DataFrame(
324 data={"diaObjectId": n_sources * [objId],
325 "band": n_sources * ["g"],
326 "diaSourceId": np.arange(n_sources, dtype=int),
327 "flags": np.zeros(n_sources, dtype=np.uint64)})
328 plug = SimpleSourceFlagDiaPlugin(SimpleSourceFlagDiaPluginConfig(),
329 "ap_diaObjectFlag",
330 None)
332 diaObjects = make_diaObject_table(objId, plug, default_value=np.uint64)
333 run_multi_plugin(diaObjects, diaSources, "g", plug)
334 self.assertEqual(diaObjects.at[objId, "flags"], 0)
336 # Test expected flags, all flags set.
337 diaSources = pd.DataFrame(
338 data={"diaObjectId": n_sources * [objId],
339 "band": n_sources * ["g"],
340 "diaSourceId": np.arange(n_sources, dtype=int),
341 "flags": np.ones(n_sources, dtype=np.uint64)})
342 diaObjects = make_diaObject_table(objId, plug, default_value=np.uint64)
343 run_multi_plugin(diaObjects, diaSources, "g", plug)
344 self.assertEqual(diaObjects.at[objId, "flags"], 1)
346 # Test expected flags, random flags.
347 diaSources = pd.DataFrame(
348 data={"diaObjectId": n_sources * [objId],
349 "band": n_sources * ["g"],
350 "diaSourceId": np.arange(n_sources, dtype=int),
351 "flags": np.random.randint(0, 2 ** 16, size=n_sources)})
353 diaObjects = make_diaObject_table(objId, plug, default_value=np.uint64)
354 run_multi_plugin(diaObjects, diaSources, "g", plug)
355 self.assertEqual(diaObjects.at[objId, "flags"], 1)
357 # Test expected flags, one flag set.
358 flag_array = np.zeros(n_sources, dtype=np.uint64)
359 flag_array[4] = 256
360 diaSources = pd.DataFrame(
361 data={"diaObjectId": n_sources * [objId],
362 "band": n_sources * ["g"],
363 "diaSourceId": np.arange(n_sources, dtype=int),
364 "flags": flag_array})
365 diaObjects = make_diaObject_table(objId, plug, default_value=np.uint64)
366 run_multi_plugin(diaObjects, diaSources, "g", plug)
367 self.assertEqual(diaObjects.at[objId, "flags"], 1)
370class TestWeightedMeanDiaPsfFlux(unittest.TestCase):
372 def testCalculate(self):
373 """Test mean value calculation.
374 """
375 n_sources = 10
376 objId = 0
378 # Test expected mean.
379 diaObjects = pd.DataFrame({"diaObjectId": [objId]})
380 diaSources = pd.DataFrame(
381 data={"diaObjectId": n_sources * [objId],
382 "band": n_sources * ["u"],
383 "diaSourceId": np.arange(n_sources, dtype=int),
384 "psfFlux": np.linspace(-1, 1, n_sources),
385 "psfFluxErr": np.ones(n_sources)})
387 plug = WeightedMeanDiaPsfFlux(WeightedMeanDiaPsfFluxConfig(),
388 "ap_meanFlux",
389 None)
390 run_multi_plugin(diaObjects, diaSources, "u", plug)
392 self.assertAlmostEqual(diaObjects.loc[objId, "u_psfFluxMean"], 0.0)
393 self.assertAlmostEqual(diaObjects.loc[objId, "u_psfFluxMeanErr"],
394 np.sqrt(1 / n_sources))
395 self.assertEqual(diaObjects.loc[objId, "u_psfFluxNdata"], n_sources)
397 # Test expected mean with a nan value.
398 diaObjects = pd.DataFrame({"diaObjectId": [objId]})
399 fluxes = np.linspace(-1, 1, n_sources)
400 fluxes[4] = np.nan
401 diaSources = pd.DataFrame(
402 data={"diaObjectId": n_sources * [objId],
403 "band": n_sources * ["r"],
404 "diaSourceId": np.arange(n_sources, dtype=int),
405 "psfFlux": fluxes,
406 "psfFluxErr": np.ones(n_sources)})
407 run_multi_plugin(diaObjects, diaSources, "r", plug)
409 self.assertAlmostEqual(diaObjects.at[objId, "r_psfFluxMean"],
410 np.nanmean(fluxes))
411 self.assertAlmostEqual(diaObjects.at[objId, "r_psfFluxMeanErr"],
412 np.sqrt(1 / (n_sources - 1)))
413 self.assertEqual(diaObjects.loc[objId, "r_psfFluxNdata"], n_sources - 1)
416class TestPercentileDiaPsfFlux(unittest.TestCase):
418 def testCalculate(self):
419 """Test flux percentile calculation.
420 """
421 n_sources = 10
422 objId = 0
424 # Test expected percentile values.
425 fluxes = np.linspace(-1, 1, n_sources)
426 diaObjects = pd.DataFrame({"diaObjectId": [objId]})
427 diaSources = pd.DataFrame(
428 data={"diaObjectId": n_sources * [objId],
429 "band": n_sources * ["u"],
430 "diaSourceId": np.arange(n_sources, dtype=int),
431 "psfFlux": fluxes,
432 "psfFluxErr": np.ones(n_sources)})
434 plug = PercentileDiaPsfFlux(PercentileDiaPsfFluxConfig(),
435 "ap_percentileFlux",
436 None)
437 run_multi_plugin(diaObjects, diaSources, "u", plug)
438 for pTile, testVal in zip(plug.config.percentiles,
439 np.nanpercentile(
440 fluxes,
441 plug.config.percentiles)):
442 self.assertAlmostEqual(
443 diaObjects.at[objId, "u_psfFluxPercentile{:02d}".format(pTile)],
444 testVal)
446 # Test expected percentile values with a nan value.
447 diaObjects = pd.DataFrame({"diaObjectId": [objId]})
448 fluxes[4] = np.nan
449 diaSources = pd.DataFrame(
450 data={"diaObjectId": n_sources * [objId],
451 "band": n_sources * ["r"],
452 "diaSourceId": np.arange(n_sources, dtype=int),
453 "psfFlux": fluxes,
454 "psfFluxErr": np.ones(n_sources)})
455 run_multi_plugin(diaObjects, diaSources, "r", plug)
456 for pTile, testVal in zip(plug.config.percentiles,
457 np.nanpercentile(
458 fluxes,
459 plug.config.percentiles)):
460 self.assertAlmostEqual(
461 diaObjects.at[objId, "r_psfFluxPercentile{:02d}".format(pTile)],
462 testVal)
465class TestSigmaDiaPsfFlux(unittest.TestCase):
467 def testCalculate(self):
468 """Test flux scatter calculation.
469 """
470 n_sources = 10
471 objId = 0
473 # Test expected sigma scatter of fluxes.
474 fluxes = np.linspace(-1, 1, n_sources)
475 diaSources = pd.DataFrame(
476 data={"diaObjectId": n_sources * [objId],
477 "band": n_sources * ["u"],
478 "diaSourceId": np.arange(n_sources, dtype=int),
479 "psfFlux": fluxes,
480 "psfFluxErr": np.ones(n_sources)})
482 plug = SigmaDiaPsfFlux(SigmaDiaPsfFluxConfig(),
483 "ap_sigmaFlux",
484 None)
485 diaObjects = make_diaObject_table(objId, plug, band='u')
486 run_multi_plugin(diaObjects, diaSources, "u", plug)
487 self.assertAlmostEqual(diaObjects.at[objId, "u_psfFluxSigma"],
488 np.nanstd(fluxes, ddof=1))
490 # test one input, returns nan.
491 diaSources = pd.DataFrame(
492 data={"diaObjectId": 1 * [objId],
493 "band": 1 * ["g"],
494 "diaSourceId": [0],
495 "psfFlux": [fluxes[0]],
496 "psfFluxErr": [1.]})
498 diaObjects = make_diaObject_table(objId, plug, band='g')
499 run_multi_plugin(diaObjects, diaSources, "g", plug)
500 self.assertTrue(np.isnan(diaObjects.at[objId, "g_psfFluxSigma"]))
502 # Test expected sigma scatter of fluxes with a nan value.
503 fluxes[4] = np.nan
504 diaSources = pd.DataFrame(
505 data={"diaObjectId": n_sources * [objId],
506 "band": n_sources * ["r"],
507 "diaSourceId": np.arange(n_sources, dtype=int),
508 "psfFlux": fluxes,
509 "psfFluxErr": np.ones(n_sources)})
511 diaObjects = make_diaObject_table(objId, plug, band='r')
512 run_multi_plugin(diaObjects, diaSources, "r", plug)
513 self.assertAlmostEqual(diaObjects.at[objId, "r_psfFluxSigma"],
514 np.nanstd(fluxes, ddof=1))
517class TestChi2DiaPsfFlux(unittest.TestCase):
519 def testCalculate(self):
520 """Test flux chi2 calculation.
521 """
522 n_sources = 10
523 objId = 0
525 # Test expected chi^2 value.
526 fluxes = np.linspace(-1, 1, n_sources)
527 diaObjects = pd.DataFrame({"diaObjectId": [objId],
528 "u_psfFluxMean": [0.0]})
529 diaSources = pd.DataFrame(
530 data={"diaObjectId": n_sources * [objId],
531 "band": n_sources * ["u"],
532 "diaSourceId": np.arange(n_sources, dtype=int),
533 "psfFlux": fluxes,
534 "psfFluxErr": np.ones(n_sources)})
536 plug = Chi2DiaPsfFlux(Chi2DiaPsfFluxConfig(),
537 "ap_chi2Flux",
538 None)
539 run_multi_plugin(diaObjects, diaSources, "u", plug)
540 self.assertAlmostEqual(
541 diaObjects.loc[objId, "u_psfFluxChi2"],
542 np.nansum(((diaSources["psfFlux"]
543 - np.nanmean(diaSources["psfFlux"]))
544 / diaSources["psfFluxErr"]) ** 2))
546 # Test expected chi^2 value with a nan value set.
547 fluxes[4] = np.nan
548 diaObjects = pd.DataFrame({"diaObjectId": [objId],
549 "r_psfFluxMean": [np.nanmean(fluxes)]})
550 diaSources = pd.DataFrame(
551 data={"diaObjectId": n_sources * [objId],
552 "band": n_sources * ["r"],
553 "diaSourceId": np.arange(n_sources, dtype=int),
554 "psfFlux": fluxes,
555 "psfFluxErr": np.ones(n_sources)})
556 run_multi_plugin(diaObjects, diaSources, "r", plug)
557 self.assertAlmostEqual(
558 diaObjects.loc[objId, "r_psfFluxChi2"],
559 np.nansum(((diaSources["psfFlux"]
560 - np.nanmean(diaSources["psfFlux"]))
561 / diaSources["psfFluxErr"]) ** 2))
564class TestMadDiaPsfFlux(unittest.TestCase):
566 def testCalculate(self):
567 """Test flux median absolute deviation calculation.
568 """
569 n_sources = 10
570 objId = 0
572 # Test expected MAD value.
573 fluxes = np.linspace(-1, 1, n_sources)
574 diaObjects = pd.DataFrame({"diaObjectId": [objId]})
575 diaSources = pd.DataFrame(
576 data={"diaObjectId": n_sources * [objId],
577 "band": n_sources * ["u"],
578 "diaSourceId": np.arange(n_sources, dtype=int),
579 "psfFlux": fluxes,
580 "psfFluxErr": np.ones(n_sources)})
582 plug = MadDiaPsfFlux(MadDiaPsfFluxConfig(),
583 "ap_madFlux",
584 None)
585 run_multi_plugin(diaObjects, diaSources, "u", plug)
586 self.assertAlmostEqual(diaObjects.at[objId, "u_psfFluxMAD"],
587 median_absolute_deviation(fluxes,
588 ignore_nan=True))
590 # Test expected MAD value with a nan set.
591 fluxes[4] = np.nan
592 diaObjects = pd.DataFrame({"diaObjectId": [objId]})
593 diaSources = pd.DataFrame(
594 data={"diaObjectId": n_sources * [objId],
595 "band": n_sources * ["r"],
596 "diaSourceId": np.arange(n_sources, dtype=int),
597 "psfFlux": fluxes,
598 "psfFluxErr": np.ones(n_sources)})
599 run_multi_plugin(diaObjects, diaSources, "r", plug)
600 self.assertAlmostEqual(diaObjects.at[objId, "r_psfFluxMAD"],
601 median_absolute_deviation(fluxes,
602 ignore_nan=True))
605class TestSkewDiaPsfFlux(unittest.TestCase):
607 def testCalculate(self):
608 """Test flux skew calculation.
609 """
610 n_sources = 10
611 objId = 0
613 # Test expected skew value.
614 fluxes = np.linspace(-1, 1, n_sources)
615 diaObjects = pd.DataFrame({"diaObjectId": [objId]})
616 diaSources = pd.DataFrame(
617 data={"diaObjectId": n_sources * [objId],
618 "band": n_sources * ["u"],
619 "diaSourceId": np.arange(n_sources, dtype=int),
620 "psfFlux": fluxes,
621 "psfFluxErr": np.ones(n_sources)})
623 plug = SkewDiaPsfFlux(SkewDiaPsfFluxConfig(),
624 "ap_skewFlux",
625 None)
626 run_multi_plugin(diaObjects, diaSources, "u", plug)
627 self.assertAlmostEqual(
628 diaObjects.loc[objId, "u_psfFluxSkew"],
629 skew_wrapper(fluxes))
631 # Test expected skew value with a nan set.
632 fluxes[4] = np.nan
633 diaObjects = pd.DataFrame({"diaObjectId": [objId]})
634 diaSources = pd.DataFrame(
635 data={"diaObjectId": n_sources * [objId],
636 "band": n_sources * ["r"],
637 "diaSourceId": np.arange(n_sources, dtype=int),
638 "psfFlux": fluxes,
639 "psfFluxErr": np.ones(n_sources)})
640 run_multi_plugin(diaObjects, diaSources, "r", plug)
642 self.assertAlmostEqual(
643 diaObjects.at[objId, "r_psfFluxSkew"],
644 skew_wrapper(fluxes))
647class TestMinMaxDiaPsfFlux(unittest.TestCase):
649 def testCalculate(self):
650 """Test flux min/max calculation.
651 """
652 n_sources = 10
653 objId = 0
655 # Test expected MinMax fluxes.
656 fluxes = np.linspace(-1, 1, n_sources)
657 diaObjects = pd.DataFrame({"diaObjectId": [objId]})
658 diaSources = pd.DataFrame(
659 data={"diaObjectId": n_sources * [objId],
660 "band": n_sources * ["u"],
661 "diaSourceId": np.arange(n_sources, dtype=int),
662 "psfFlux": fluxes,
663 "psfFluxErr": np.ones(n_sources)})
665 plug = MinMaxDiaPsfFlux(MinMaxDiaPsfFluxConfig(),
666 "ap_minMaxFlux",
667 None)
668 run_multi_plugin(diaObjects, diaSources, "u", plug)
669 self.assertEqual(diaObjects.loc[objId, "u_psfFluxMin"], -1)
670 self.assertEqual(diaObjects.loc[objId, "u_psfFluxMax"], 1)
672 # Test expected MinMax fluxes with a nan set.
673 fluxes[4] = np.nan
674 diaObjects = pd.DataFrame({"diaObjectId": [objId]})
675 diaSources = pd.DataFrame(
676 data={"diaObjectId": n_sources * [objId],
677 "band": n_sources * ["r"],
678 "diaSourceId": np.arange(n_sources, dtype=int),
679 "psfFlux": fluxes,
680 "psfFluxErr": np.ones(n_sources)})
681 run_multi_plugin(diaObjects, diaSources, "r", plug)
682 self.assertEqual(diaObjects.loc[objId, "r_psfFluxMin"], -1)
683 self.assertEqual(diaObjects.loc[objId, "r_psfFluxMax"], 1)
686class TestMaxSlopeDiaPsfFlux(unittest.TestCase):
688 def testCalculate(self):
689 """Test flux maximum slope.
690 """
691 n_sources = 10
692 objId = 0
694 # Test max slope value.
695 fluxes = np.linspace(-1, 1, n_sources)
696 times = np.concatenate([np.linspace(0, 1, n_sources)[:-1], [1 - 1/90]])
697 diaSources = pd.DataFrame(
698 data={"diaObjectId": n_sources * [objId],
699 "band": n_sources * ["u"],
700 "diaSourceId": np.arange(n_sources, dtype=int),
701 "psfFlux": fluxes,
702 "psfFluxErr": np.ones(n_sources),
703 "midpointMjdTai": times})
705 plug = MaxSlopeDiaPsfFlux(MaxSlopeDiaPsfFluxConfig(),
706 "ap_maxSlopeFlux",
707 None)
708 diaObjects = make_diaObject_table(objId, plug, band='u')
709 run_multi_plugin(diaObjects, diaSources, "u", plug)
710 self.assertAlmostEqual(diaObjects.at[objId, "u_psfFluxMaxSlope"], 2 + 2/9)
712 # Test max slope value returns nan on 1 input.
713 diaSources = pd.DataFrame(
714 data={"diaObjectId": 1 * [objId],
715 "band": 1 * ["g"],
716 "diaSourceId": np.arange(1, dtype=int),
717 "psfFlux": fluxes[0],
718 "psfFluxErr": np.ones(1),
719 "midpointMjdTai": times[0]})
720 diaObjects = make_diaObject_table(objId, plug, band='g')
721 run_multi_plugin(diaObjects, diaSources, "g", plug)
722 self.assertTrue(np.isnan(diaObjects.at[objId, "g_psfFluxMaxSlope"]))
724 # Test max slope value inputing nan values.
725 fluxes[4] = np.nan
726 times[7] = np.nan
727 diaSources = pd.DataFrame(
728 data={"diaObjectId": n_sources * [objId],
729 "band": n_sources * ["r"],
730 "diaSourceId": np.arange(n_sources, dtype=int),
731 "psfFlux": fluxes,
732 "psfFluxErr": np.ones(n_sources),
733 "midpointMjdTai": times})
734 diaObjects = make_diaObject_table(objId, plug, band='r')
735 run_multi_plugin(diaObjects, diaSources, "r", plug)
736 self.assertAlmostEqual(diaObjects.at[objId, "r_psfFluxMaxSlope"], 2 + 2 / 9)
739class TestErrMeanDiaPsfFlux(unittest.TestCase):
741 def testCalculate(self):
742 """Test error mean calculation.
743 """
744 n_sources = 10
745 objId = 0
747 # Test mean of the errors.
748 fluxes = np.linspace(-1, 1, n_sources)
749 errors = np.linspace(1, 2, n_sources)
750 diaObjects = pd.DataFrame({"diaObjectId": [objId]})
751 diaSources = pd.DataFrame(
752 data={"diaObjectId": n_sources * [objId],
753 "band": n_sources * ["u"],
754 "diaSourceId": np.arange(n_sources, dtype=int),
755 "psfFlux": fluxes,
756 "psfFluxErr": errors})
758 plug = ErrMeanDiaPsfFlux(ErrMeanDiaPsfFluxConfig(),
759 "ap_errMeanFlux",
760 None)
761 run_multi_plugin(diaObjects, diaSources, "u", plug)
762 self.assertAlmostEqual(diaObjects.at[objId, "u_psfFluxErrMean"],
763 np.nanmean(errors))
765 # Test mean of the errors with input nan value.
766 errors[4] = np.nan
767 diaObjects = pd.DataFrame({"diaObjectId": [objId]})
768 diaSources = pd.DataFrame(
769 data={"diaObjectId": n_sources * [objId],
770 "band": n_sources * ["r"],
771 "diaSourceId": np.arange(n_sources, dtype=int),
772 "psfFlux": fluxes,
773 "psfFluxErr": errors})
774 run_multi_plugin(diaObjects, diaSources, "r", plug)
775 self.assertAlmostEqual(diaObjects.at[objId, "r_psfFluxErrMean"],
776 np.nanmean(errors))
779class TestLinearFitDiaPsfFlux(unittest.TestCase):
781 def testCalculate(self):
782 """Test a linear fit to flux vs time.
783 """
784 n_sources = 10
785 objId = 0
787 # Test best fit linear model.
788 fluxes = np.linspace(-1, 1, n_sources)
789 errors = np.linspace(1, 2, n_sources)
790 times = np.linspace(0, 1, n_sources)
791 diaObjects = pd.DataFrame({"diaObjectId": [objId]})
792 diaSources = pd.DataFrame(
793 data={"diaObjectId": n_sources * [objId],
794 "band": n_sources * ["u"],
795 "diaSourceId": np.arange(n_sources, dtype=int),
796 "psfFlux": fluxes,
797 "psfFluxErr": errors,
798 "midpointMjdTai": times})
800 plug = LinearFitDiaPsfFlux(LinearFitDiaPsfFluxConfig(),
801 "ap_LinearFit",
802 None)
803 run_multi_plugin(diaObjects, diaSources, "u", plug)
804 self.assertAlmostEqual(diaObjects.loc[objId, "u_psfFluxLinearSlope"],
805 2.)
806 self.assertAlmostEqual(diaObjects.loc[objId, "u_psfFluxLinearIntercept"],
807 -1.)
809 # Test best fit linear model with input nans.
810 fluxes[7] = np.nan
811 errors[4] = np.nan
812 times[2] = np.nan
813 diaObjects = pd.DataFrame({"diaObjectId": [objId]})
814 diaSources = pd.DataFrame(
815 data={"diaObjectId": n_sources * [objId],
816 "band": n_sources * ["r"],
817 "diaSourceId": np.arange(n_sources, dtype=int),
818 "psfFlux": fluxes,
819 "psfFluxErr": errors,
820 "midpointMjdTai": times})
821 run_multi_plugin(diaObjects, diaSources, "r", plug)
822 self.assertAlmostEqual(diaObjects.loc[objId, "r_psfFluxLinearSlope"], 2.)
823 self.assertAlmostEqual(diaObjects.loc[objId, "r_psfFluxLinearIntercept"],
824 -1.)
827class TestStetsonJDiaPsfFlux(unittest.TestCase):
829 def testCalculate(self):
830 """Test the stetsonJ statistic.
831 """
832 n_sources = 10
833 objId = 0
835 # Test stetsonJ calculation.
836 fluxes = np.linspace(-1, 1, n_sources)
837 errors = np.ones(n_sources)
838 diaObjects = pd.DataFrame({"diaObjectId": [objId],
839 "u_psfFluxMean": [np.nanmean(fluxes)]})
840 diaSources = pd.DataFrame(
841 data={"diaObjectId": n_sources * [objId],
842 "band": n_sources * ["u"],
843 "diaSourceId": np.arange(n_sources, dtype=int),
844 "psfFlux": fluxes,
845 "psfFluxErr": errors})
847 plug = StetsonJDiaPsfFlux(StetsonJDiaPsfFluxConfig(),
848 "ap_StetsonJ",
849 None)
850 run_multi_plugin(diaObjects, diaSources, "u", plug)
851 # Expected StetsonJ for the values created. Confirmed using Cesimum's
852 # implementation. http://github.com/cesium-ml/cesium
853 self.assertAlmostEqual(diaObjects.loc[objId, "u_psfFluxStetsonJ"],
854 -0.5958393936080928)
856 # Test stetsonJ calculation returns nan on single input.
857 diaObjects = pd.DataFrame({"diaObjectId": [objId],
858 "g_psfFluxMean": [np.nanmean(fluxes)]})
859 diaSources = pd.DataFrame(
860 data={"diaObjectId": 1 * [objId],
861 "band": 1 * ["g"],
862 "diaSourceId": np.arange(1, dtype=int),
863 "psfFlux": fluxes[0],
864 "psfFluxErr": errors[0]})
865 run_multi_plugin(diaObjects, diaSources, "g", plug)
866 self.assertTrue(np.isnan(diaObjects.at[objId, "g_psfFluxStetsonJ"]))
868 # Test stetsonJ calculation returns when nans are input.
869 fluxes[7] = np.nan
870 errors[4] = np.nan
871 nonNanMask = np.logical_and(~np.isnan(fluxes),
872 ~np.isnan(errors))
873 diaObjects = pd.DataFrame(
874 {"diaObjectId": [objId],
875 "r_psfFluxMean": [np.average(fluxes[nonNanMask],
876 weights=errors[nonNanMask])]})
877 diaSources = pd.DataFrame(
878 data={"diaObjectId": n_sources * [objId],
879 "band": n_sources * ["r"],
880 "diaSourceId": np.arange(n_sources, dtype=int),
881 "psfFlux": fluxes,
882 "psfFluxErr": errors})
883 run_multi_plugin(diaObjects, diaSources, "r", plug)
884 self.assertAlmostEqual(diaObjects.at[objId, "r_psfFluxStetsonJ"],
885 -0.5412797916187173)
888class TestWeightedMeanDiaTotFlux(unittest.TestCase):
890 def testCalculate(self):
891 """Test mean value calculation.
892 """
893 n_sources = 10
894 objId = 0
896 # Test test mean on scienceFlux.
897 fluxes = np.linspace(-1, 1, n_sources)
898 diaObjects = pd.DataFrame({"diaObjectId": [objId]})
899 diaSources = pd.DataFrame(
900 data={"diaObjectId": n_sources * [objId],
901 "band": n_sources * ["u"],
902 "diaSourceId": np.arange(n_sources, dtype=int),
903 "scienceFlux": fluxes,
904 "scienceFluxErr": np.ones(n_sources)})
906 plug = WeightedMeanDiaTotFlux(WeightedMeanDiaTotFluxConfig(),
907 "ap_meanTotFlux",
908 None)
909 run_multi_plugin(diaObjects, diaSources, "u", plug)
911 self.assertAlmostEqual(diaObjects.at[objId, "u_scienceFluxMean"], 0.0)
912 self.assertAlmostEqual(diaObjects.at[objId, "u_scienceFluxMeanErr"],
913 np.sqrt(1 / n_sources))
915 # Test test mean on scienceFlux with input nans
916 fluxes[4] = np.nan
917 diaObjects = pd.DataFrame({"diaObjectId": [objId]})
918 diaSources = pd.DataFrame(
919 data={"diaObjectId": n_sources * [objId],
920 "band": n_sources * ["r"],
921 "diaSourceId": np.arange(n_sources, dtype=int),
922 "scienceFlux": fluxes,
923 "scienceFluxErr": np.ones(n_sources)})
924 run_multi_plugin(diaObjects, diaSources, "r", plug)
926 self.assertAlmostEqual(diaObjects.at[objId, "r_scienceFluxMean"],
927 np.nanmean(fluxes))
928 self.assertAlmostEqual(diaObjects.at[objId, "r_scienceFluxMeanErr"],
929 np.sqrt(1 / (n_sources - 1)))
932def generatePeriodicData(n=10, period=10):
933 """Generate noisy, sinusoidally-varying periodic data for testing Lomb-
934 Scargle Periodogram.
936 The returned fluxes will have, within the errors, the passed-in period and
937 a power close to 1, because the fluxes are purely sinusoidal.
939 Parameters
940 ----------
941 n : int
942 Number of data points to generate.
943 period : float
944 Period of the periodic signal.
946 Returns
947 -------
948 t : np.ndarray
949 Time values.
950 y_obs : np.ndarray
951 Observed flux values.
952 """
953 np.random.seed(42)
955 t = np.linspace(-2*np.pi, 2*np.pi, n) + 100*np.random.random(n)
956 y = 10 + np.sin(2 * np.pi * t / period)
957 y_obs = np.random.normal(y, 0.001)
959 return t, y_obs
962class TestMultiLombScarglePeriodogram(lsst.utils.tests.TestCase):
964 def testCalculate(self):
965 """Test Mulitband Lomb Scargle Periodogram."""
966 n_sources = 10
967 objId = 0
969 # Create synthetic multi-band data
970 times, fluxes = generatePeriodicData(n_sources, period=10)
971 diaObjects = pd.DataFrame({"diaObjectId": [objId]})
972 diaSources = pd.DataFrame(
973 data={"diaObjectId": n_sources * [objId],
974 "band": n_sources//2 * ["u"] + n_sources//2 * ["g"],
975 "diaSourceId": np.arange(n_sources, dtype=int),
976 "midpointMjdTai": times,
977 "psfFlux": fluxes,
978 "psfFluxErr": 1e-3+np.zeros(n_sources)})
980 plugin = LombScarglePeriodogramMulti(LombScarglePeriodogramMultiConfig(),
981 "ap_lombScarglePeriodogramMulti",
982 None)
984 run_multiband_plugin(diaObjects, diaSources, plugin)
985 self.assertAlmostEqual(diaObjects.at[objId, "multiPeriod"], 10, delta=0.04)
986 self.assertAlmostEqual(diaObjects.at[objId, "multiPower"], 1, delta=1e-2)
987 # This implementation of LS returns a normalized power < 1.
988 self.assertLess(diaObjects.at[objId, "multiPower"], 1)
989 self.assertAlmostEqual(diaObjects.at[objId, "multiFap"], 0, delta=0.04)
990 # Note: The below values are empirical, but seem reasonable, and
991 # test that we get values for each band.
992 self.assertAlmostEqual(diaObjects.at[objId, "u_multiAmp"], 0.029, delta=0.01)
993 self.assertAlmostEqual(diaObjects.at[objId, "g_multiAmp"], 0.029, delta=0.01)
994 self.assertAlmostEqual(diaObjects.at[objId, "u_multiPhase"], -2.0, delta=0.2)
995 self.assertAlmostEqual(diaObjects.at[objId, "g_multiPhase"], 1.0, delta=0.1)
997 def testCalculateTwoSources(self):
998 """Test Mulitband Lomb Scargle Periodogram with 2 sources (minimum
999 detections = 5), which will result in NaN output."""
1000 objId = 0
1001 n_sources = 2
1002 times, fluxes = generatePeriodicData(n_sources, period=10)
1003 diaObjects = pd.DataFrame({"diaObjectId": [objId]})
1004 diaSources = pd.DataFrame(
1005 data={"diaObjectId": n_sources * [objId],
1006 "band": n_sources * ["u"],
1007 "diaSourceId": np.arange(n_sources, dtype=int),
1008 "midpointMjdTai": times,
1009 "psfFlux": fluxes,
1010 "psfFluxErr": 1e-3+np.zeros(n_sources)})
1012 plugin = LombScarglePeriodogramMulti(LombScarglePeriodogramMultiConfig(),
1013 "ap_lombScarglePeriodogramMulti",
1014 None)
1016 run_multi_plugin(diaObjects, diaSources, "u", plugin)
1017 self.assertTrue(np.isnan(diaObjects.at[objId, "multiPeriod"]))
1018 self.assertTrue(np.isnan(diaObjects.at[objId, "multiPower"]))
1019 self.assertTrue(np.isnan(diaObjects.at[objId, "multiFap"]))
1022class TestLombScarglePeriodogram(lsst.utils.tests.TestCase):
1024 def testCalculate(self):
1025 """Test Lomb Scargle Periodogram."""
1026 n_sources = 10
1027 objId = 0
1029 # Test period calculation.
1030 times, fluxes = generatePeriodicData(n_sources, period=10)
1031 diaObjects = pd.DataFrame({"diaObjectId": [objId]})
1032 diaSources = pd.DataFrame(
1033 data={"diaObjectId": n_sources * [objId],
1034 "band": n_sources * ["u"],
1035 "diaSourceId": np.arange(n_sources, dtype=int),
1036 "midpointMjdTai": times,
1037 "psfFlux": fluxes,
1038 "psfFluxErr": 1e-3+np.zeros(n_sources)})
1040 plugin = LombScarglePeriodogram(LombScarglePeriodogramConfig(),
1041 "ap_lombScarglePeriodogram",
1042 None)
1044 run_multi_plugin(diaObjects, diaSources, "u", plugin)
1045 self.assertAlmostEqual(diaObjects.at[objId, "u_period"], 10, delta=0.04)
1046 # This implementation of LS returns a normalized power < 1.
1047 self.assertAlmostEqual(diaObjects.at[objId, "u_power"], 1, delta=1e-2)
1048 self.assertLess(diaObjects.at[objId, "u_power"], 1)
1050 # Test that we get the same result with a NaN flux.
1051 diaSources.loc[4, "psfFlux"] = np.nan
1052 diaObjects = pd.DataFrame({"diaObjectId": [objId]})
1053 diaSources = pd.DataFrame(
1054 data={"diaObjectId": n_sources * [objId],
1055 "band": n_sources * ["r"],
1056 "diaSourceId": np.arange(n_sources, dtype=int),
1057 "midpointMjdTai": times,
1058 "psfFlux": fluxes,
1059 "psfFluxErr": np.ones(n_sources)})
1060 run_multi_plugin(diaObjects, diaSources, "r", plugin)
1061 self.assertAlmostEqual(diaObjects.at[objId, "r_period"], 10, delta=0.04)
1062 self.assertAlmostEqual(diaObjects.at[objId, "r_power"], 1, delta=1e-2)
1063 # This implementation of LS returns a normalized power < 1.
1064 self.assertLess(diaObjects.at[objId, "r_power"], 1)
1067class TestSigmaDiaTotFlux(unittest.TestCase):
1069 def testCalculate(self):
1070 """Test flux scatter calculation.
1071 """
1072 n_sources = 10
1073 objId = 0
1075 # Test test scatter on scienceFlux.
1076 fluxes = np.linspace(-1, 1, n_sources)
1077 diaSources = pd.DataFrame(
1078 data={"diaObjectId": n_sources * [objId],
1079 "band": n_sources * ["u"],
1080 "diaSourceId": np.arange(n_sources, dtype=int),
1081 "scienceFlux": fluxes,
1082 "scienceFluxErr": np.ones(n_sources)})
1084 plug = SigmaDiaTotFlux(SigmaDiaTotFluxConfig(),
1085 "ap_sigmaTotFlux",
1086 None)
1087 diaObjects = make_diaObject_table(objId, plug, band='u')
1088 run_multi_plugin(diaObjects, diaSources, "u", plug)
1089 self.assertAlmostEqual(diaObjects.at[objId, "u_scienceFluxSigma"],
1090 np.nanstd(fluxes, ddof=1))
1092 # Test test scatter on scienceFlux returns nan on 1 input.
1093 diaSources = pd.DataFrame(
1094 data={"diaObjectId": 1 * [objId],
1095 "band": 1 * ["g"],
1096 "diaSourceId": np.arange(1, dtype=int),
1097 "scienceFlux": fluxes[0],
1098 "scienceFluxErr": np.ones(1)})
1099 diaObjects = make_diaObject_table(objId, plug, band='g')
1100 run_multi_plugin(diaObjects, diaSources, "g", plug)
1101 self.assertTrue(np.isnan(diaObjects.at[objId, "g_scienceFluxSigma"]))
1103 # Test test scatter on scienceFlux takes input nans.
1104 fluxes[4] = np.nan
1105 diaSources = pd.DataFrame(
1106 data={"diaObjectId": n_sources * [objId],
1107 "band": n_sources * ["r"],
1108 "diaSourceId": np.arange(n_sources, dtype=int),
1109 "scienceFlux": fluxes,
1110 "scienceFluxErr": np.ones(n_sources)})
1111 diaObjects = make_diaObject_table(objId, plug, band='r')
1112 run_multi_plugin(diaObjects, diaSources, "r", plug)
1113 self.assertAlmostEqual(diaObjects.at[objId, "r_scienceFluxSigma"],
1114 np.nanstd(fluxes, ddof=1))
1117def skew_wrapper(values):
1118 """Compute scipy skew, omitting nans.
1120 This version works with both scipy<1.9 (where it erroneously returns a
1121 masked array) and scipy>=1.9 (where it correctly returns a float).
1123 Parameters
1124 ----------
1125 values : `np.ndarray`
1127 Returns
1128 -------
1129 skew_value : `float`
1130 """
1131 value = skew(values, bias=False, nan_policy="omit")
1132 if isinstance(value, np.ma.masked_array):
1133 return value.data
1134 else:
1135 return value
1138class MemoryTester(lsst.utils.tests.MemoryTestCase):
1139 pass
1142def setup_module(module):
1143 lsst.utils.tests.init()
1146if __name__ == "__main__": 1146 ↛ 1147line 1146 didn't jump to line 1147 because the condition on line 1146 was never true
1147 lsst.utils.tests.init()
1148 unittest.main()