Coverage for python/lsst/meas/base/diaCalculationPlugins.py: 40%

1# This file is part of ap_association.

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.

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.sphgeom as sphgeom

39from astropy.timeseries import LombScargle

40from astropy.timeseries import LombScargleMultiband

41import math

42import statistics

44from .diaCalculation import (

45 DiaObjectCalculationPluginConfig,

46 DiaObjectCalculationPlugin)

47from .pluginRegistry import register

50__all__ = ("MeanDiaPositionConfig", "MeanDiaPosition",

51 "HTMIndexDiaPosition", "HTMIndexDiaPositionConfig",

52 "NumDiaSourcesDiaPlugin", "NumDiaSourcesDiaPluginConfig",

53 "SimpleSourceFlagDiaPlugin", "SimpleSourceFlagDiaPluginConfig",

54 "WeightedMeanDiaPsfFluxConfig", "WeightedMeanDiaPsfFlux",

55 "PercentileDiaPsfFlux", "PercentileDiaPsfFluxConfig",

56 "SigmaDiaPsfFlux", "SigmaDiaPsfFluxConfig",

57 "Chi2DiaPsfFlux", "Chi2DiaPsfFluxConfig",

58 "MadDiaPsfFlux", "MadDiaPsfFluxConfig",

59 "SkewDiaPsfFlux", "SkewDiaPsfFluxConfig",

60 "MinMaxDiaPsfFlux", "MinMaxDiaPsfFluxConfig",

61 "MaxSlopeDiaPsfFlux", "MaxSlopeDiaPsfFluxConfig",

62 "ErrMeanDiaPsfFlux", "ErrMeanDiaPsfFluxConfig",

63 "LinearFitDiaPsfFlux", "LinearFitDiaPsfFluxConfig",

64 "StetsonJDiaPsfFlux", "StetsonJDiaPsfFluxConfig",

65 "WeightedMeanDiaTotFlux", "WeightedMeanDiaTotFluxConfig",

66 "SigmaDiaTotFlux", "SigmaDiaTotFluxConfig",

67 "LombScarglePeriodogram", "LombScarglePeriodogramConfig",

68 "LombScarglePeriodogramMulti", "LombScarglePeriodogramMultiConfig")

71def catchWarnings(_func=None, *, warns=[]):

72 """Decorator for generically catching numpy warnings.

73 """

74 def decoratorCatchWarnings(func):

75 @functools.wraps(func)

76 def wrapperCatchWarnings(*args, **kwargs):

77 with warnings.catch_warnings():

78 for val in warns:

79 warnings.filterwarnings("ignore", val)

80 return func(*args, **kwargs)

81 return wrapperCatchWarnings

83 if _func is None: 83 ↛ 86line 83 didn't jump to line 86 because the condition on line 83 was always true

84 return decoratorCatchWarnings

85 else:

86 return decoratorCatchWarnings(_func)

89def compute_optimized_periodogram_grid(x0, oversampling_factor=5, nyquist_factor=100):

90 """

91 Computes an optimized periodogram frequency grid for a given time series.

93 Parameters

94 ----------

95 x0 : `array`

96 The input time axis.

97 oversampling_factor : `int`, optional

98 The oversampling factor for frequency grid.

99 nyquist_factor : `int`, optional

100 The Nyquist factor for frequency grid.

101

102 Returns

103 -------

104 frequencies : `array`

105 The computed optimized periodogram frequency grid.

106 """

107

108 num_points = len(x0)

109 baseline = np.max(x0) - np.min(x0)

110

111 # Calculate the frequency resolution based on oversampling factor and baseline

112 frequency_resolution = 1. / baseline / oversampling_factor

113

114 num_frequencies = int(

115 0.5 * oversampling_factor * nyquist_factor * num_points)

116 frequencies = frequency_resolution + \

117 frequency_resolution * np.arange(num_frequencies)

118

119 return frequencies

120

121

122class LombScarglePeriodogramConfig(DiaObjectCalculationPluginConfig):

123 pass

124

125

126@register("ap_lombScarglePeriodogram")

127class LombScarglePeriodogram(DiaObjectCalculationPlugin):

128 """Compute the single-band period of a DiaObject given a set of DiaSources.

129 """

130 ConfigClass = LombScarglePeriodogramConfig

131

132 plugType = "multi"

133 outputCols = ["period", "power"]

134 needsFilter = True

135

136 @classmethod

137 def getExecutionOrder(cls):

138 return cls.DEFAULT_CATALOGCALCULATION

139

140 @catchWarnings(warns=["All-NaN slice encountered"])

141 def calculate(self,

142 diaObjects,

143 diaSources,

144 filterDiaSources,

145 band):

146 """Compute the periodogram.

147

148 Parameters

149 ----------

150 diaObjects : `pandas.DataFrame`

151 Summary objects to store values in.

152 diaSources : `pandas.DataFrame` or `pandas.DataFrameGroupBy`

153 Catalog of DiaSources summarized by this DiaObject.

154 """

155

156 # Check and initialize output columns in diaObjects.

157 if (periodCol := f"{band}_period") not in diaObjects.columns:

158 diaObjects[periodCol] = np.nan

159 if (powerCol := f"{band}_power") not in diaObjects.columns:

160 diaObjects[powerCol] = np.nan

161

162 def _calculate_period(df, min_detections=5, nterms=1, oversampling_factor=5, nyquist_factor=100):

163 """Compute the Lomb-Scargle periodogram given a set of DiaSources.

164

165 Parameters

166 ----------

167 df : `pandas.DataFrame`

168 The input DataFrame.

169 min_detections : `int`, optional

170 The minimum number of detections.

171 nterms : `int`, optional

172 The number of terms in the Lomb-Scargle model.

173 oversampling_factor : `int`, optional

174 The oversampling factor for frequency grid.

175 nyquist_factor : `int`, optional

176 The Nyquist factor for frequency grid.

177

178 Returns

179 -------

180 pd_tab : `pandas.Series`

181 The output DataFrame with the Lomb-Scargle parameters.

182 """

183 tmpDf = df[~np.logical_or(np.isnan(df["psfFlux"]),

184 np.isnan(df["midpointMjdTai"]))]

185

186 if len(tmpDf) < min_detections:

187 return pd.Series({periodCol: np.nan, powerCol: np.nan})

188

189 time = tmpDf["midpointMjdTai"].to_numpy()

190 flux = tmpDf["psfFlux"].to_numpy()

191 flux_err = tmpDf["psfFluxErr"].to_numpy()

192

193 lsp = LombScargle(time, flux, dy=flux_err, nterms=nterms)

194 f_grid = compute_optimized_periodogram_grid(

195 time, oversampling_factor=oversampling_factor, nyquist_factor=nyquist_factor)

196 period = 1/f_grid

197 power = lsp.power(f_grid)

198

199 return pd.Series({periodCol: period[np.argmax(power)],

200 powerCol: np.max(power)})

201

202 diaObjects.loc[:, [periodCol, powerCol]

203 ] = filterDiaSources.apply(_calculate_period)

204

205

206class LombScarglePeriodogramMultiConfig(DiaObjectCalculationPluginConfig):

207 pass

208

209

210@register("ap_lombScarglePeriodogramMulti")

211class LombScarglePeriodogramMulti(DiaObjectCalculationPlugin):

212 """Compute the multi-band LombScargle periodogram of a DiaObject given a set of DiaSources.

213 """

214 ConfigClass = LombScarglePeriodogramMultiConfig

215

216 plugType = "multi"

217 outputCols = ["multiPeriod", "multiPower",

218 "multiFap", "multiAmp", "multiPhase"]

219 needsFilter = True

220

221 @classmethod

222 def getExecutionOrder(cls):

223 return cls.DEFAULT_CATALOGCALCULATION

224

225 @staticmethod

226 def calculate_baluev_fap(time, n, maxPeriod, zmax):

227 """Calculate the False-Alarm probability using the Baluev approximation.

228

229 Parameters

230 ----------

231 time : `array`

232 The input time axis.

233 n : `int`

234 The number of detections.

235 maxPeriod : `float`

236 The maximum period in the grid.

237 zmax : `float`

238 The maximum power in the grid.

239

240 Returns

241 -------

242 fap_estimate : `float`

243 The False-Alarm probability Baluev approximation.

244

245 Notes

246 ----------

247 .. [1] Baluev, R. V. 2008, MNRAS, 385, 1279

248 .. [2] Süveges, M., Guy, L.P., Eyer, L., et al. 2015, MNRAS, 450, 2052

249 """

250 if n <= 2:

251 return np.nan

252

253 gam_ratio = math.factorial(int((n - 1)/2)) / math.factorial(int((n - 2)/2))

254 fu = 1/maxPeriod

255 return gam_ratio * np.sqrt(

256 4*np.pi*statistics.variance(time)

257 ) * fu * (1-zmax)**((n-4)/2) * np.sqrt(zmax)

258

259 @staticmethod

260 def generate_lsp_params(lsp_model, fbest, bands):

261 """Generate the Lomb-Scargle parameters.

262 Parameters

263 ----------

264 lsp_model : `astropy.timeseries.LombScargleMultiband`

265 The Lomb-Scargle model.

266 fbest : `float`

267 The best period.

268 bands : `array`

269 The bands of the time series.

270

271 Returns

272 -------

273 Amp : `array`

274 The amplitude of the time series.

275 Ph : `array`

276 The phase of the time series.

277

278 Notes

279 ----------

280 .. [1] VanderPlas, J. T., & Ivezic, Z. 2015, ApJ, 812, 18

281 """

282 best_params = lsp_model.model_parameters(fbest, units=True)

283

284 name_params = [f"theta_base_{i}" for i in range(3)]

285 name_params += [f"theta_band_{band}_{i}" for band in np.unique(bands) for i in range(3)]

286

287 df_params = pd.DataFrame([best_params], columns=name_params)

288

289 unique_bands = np.unique(bands)

290

291 amplitude_band = [np.sqrt(df_params[f"theta_band_{band}_1"]**2

292 + df_params[f"theta_band_{band}_2"]**2)

293 for band in unique_bands]

294 phase_bands = [np.arctan2(df_params[f"theta_band_{band}_2"],

295 df_params[f"theta_band_{band}_1"]) for band in unique_bands]

296

297 amp = [a[0] for a in amplitude_band]

298 ph = [p[0] for p in phase_bands]

299

300 return amp, ph

301

302 @catchWarnings(warns=["All-NaN slice encountered"])

303 def calculate(self,

304 diaObjects,

305 diaSources,

306 **kwargs):

307 """Compute the multi-band LombScargle periodogram of a DiaObject given

308 a set of DiaSources.

309

310 Parameters

311 ----------

312 diaObjects : `pandas.DataFrame`

313 Summary objects to store values in.

314 diaSources : `pandas.DataFrame` or `pandas.DataFrameGroupBy`

315 Catalog of DiaSources summarized by this DiaObject.

316 **kwargs : `dict`

317 Unused kwargs that are always passed to a plugin.

318 """

319

320 bands_arr = diaSources['band'].unique().values

321 unique_bands = np.unique(np.concatenate(bands_arr))

322 # Check and initialize output columns in diaObjects.

323 if (periodCol := "multiPeriod") not in diaObjects.columns:

324 diaObjects[periodCol] = np.nan

325 if (powerCol := "multiPower") not in diaObjects.columns:

326 diaObjects[powerCol] = np.nan

327 if (fapCol := "multiFap") not in diaObjects.columns:

328 diaObjects[fapCol] = np.nan

329 ampCol = "multiAmp"

330 phaseCol = "multiPhase"

331 for i in range(len(unique_bands)):

332 ampCol_band = f"{unique_bands[i]}_{ampCol}"

333 if ampCol_band not in diaObjects.columns:

334 diaObjects[ampCol_band] = np.nan

335 phaseCol_band = f"{unique_bands[i]}_{phaseCol}"

336 if phaseCol_band not in diaObjects.columns:

337 diaObjects[phaseCol_band] = np.nan

338

339 def _calculate_period_multi(df, all_unique_bands,

340 min_detections=9, oversampling_factor=5, nyquist_factor=100):

341 """Calculate the multi-band Lomb-Scargle periodogram.

342

343 Parameters

344 ----------

345 df : `pandas.DataFrame`

346 The input DataFrame.

347 all_unique_bands : `list` of `str`

348 List of all bands present in the diaSource table that is being worked on.

349 min_detections : `int`, optional

350 The minimum number of detections, including all bands.

351 oversampling_factor : `int`, optional

352 The oversampling factor for frequency grid.

353 nyquist_factor : `int`, optional

354 The Nyquist factor for frequency grid.

355

356 Returns

357 -------

358 pd_tab : `pandas.Series`

359 The output DataFrame with the Lomb-Scargle parameters.

360 """

361 tmpDf = df[~np.logical_or(np.isnan(df["psfFlux"]),

362 np.isnan(df["midpointMjdTai"]))]

363

364 if (len(tmpDf)) < min_detections:

365 pd_tab_nodet = pd.Series({periodCol: np.nan,

366 powerCol: np.nan,

367 fapCol: np.nan})

368 for band in all_unique_bands:

369 pd_tab_nodet[f"{band}_{ampCol}"] = np.nan

370 pd_tab_nodet[f"{band}_{phaseCol}"] = np.nan

371

372 return pd_tab_nodet

373

374 time = tmpDf["midpointMjdTai"].to_numpy()

375 flux = tmpDf["psfFlux"].to_numpy()

376 flux_err = tmpDf["psfFluxErr"].to_numpy()

377 bands = tmpDf["band"].to_numpy()

378

379 lsp = LombScargleMultiband(time, flux, bands, dy=flux_err,

380 nterms_base=1, nterms_band=1)

381

382 f_grid = compute_optimized_periodogram_grid(

383 time, oversampling_factor=oversampling_factor, nyquist_factor=nyquist_factor)

384 period = 1/f_grid

385 power = lsp.power(f_grid)

386

387 fap_estimate = self.calculate_baluev_fap(

388 time, len(time), period[np.argmax(power)], np.max(power))

389

390 params_table_new = self.generate_lsp_params(lsp, f_grid[np.argmax(power)], bands)

391

392 pd_tab = pd.Series({periodCol: period[np.argmax(power)],

393 powerCol: np.max(power),

394 fapCol: fap_estimate

395 })

396

397 # Initialize the per-band amplitude/phase columns as NaNs

398 for band in all_unique_bands:

399 pd_tab[f"{band}_{ampCol}"] = np.nan

400 pd_tab[f"{band}_{phaseCol}"] = np.nan

401

402 # Populate the values of only the bands that have data for this diaSource

403 unique_bands = np.unique(bands)

404 for i in range(len(unique_bands)):

405 pd_tab[f"{unique_bands[i]}_{ampCol}"] = params_table_new[0][i]

406 pd_tab[f"{unique_bands[i]}_{phaseCol}"] = params_table_new[1][i]

407

408 return pd_tab

409

410 columns_list = [periodCol, powerCol, fapCol]

411 for i in range(len(unique_bands)):

412 columns_list.append(f"{unique_bands[i]}_{ampCol}")

413 columns_list.append(f"{unique_bands[i]}_{phaseCol}")

414

415 diaObjects.loc[:, columns_list

416 ] = diaSources.apply(_calculate_period_multi, unique_bands)

417

418

419class MeanDiaPositionConfig(DiaObjectCalculationPluginConfig):

420 pass

421

422

423@register("ap_meanPosition")

424class MeanDiaPosition(DiaObjectCalculationPlugin):

425 """Compute the mean position of a DiaObject given a set of DiaSources.

426 """

427

428 ConfigClass = MeanDiaPositionConfig

429

430 plugType = 'multi'

431

432 outputCols = ["ra", "dec"]

433 needsFilter = False

434

435 @classmethod

436 def getExecutionOrder(cls):

437 return cls.DEFAULT_CATALOGCALCULATION

438

439 def calculate(self, diaObjects, diaSources, **kwargs):

440 """Compute the mean ra/dec position of the diaObject given the

441 diaSource locations.

442

443 Parameters

444 ----------

445 diaObjects : `pandas.DataFrame`

446 Summary objects to store values in.

447 diaSources : `pandas.DataFrame` or `pandas.DataFrameGroupBy`

448 Catalog of DiaSources summarized by this DiaObject.

449 **kwargs

450 Any additional keyword arguments that may be passed to the plugin.

451 """

452 for outCol in self.outputCols:

453 if outCol not in diaObjects.columns:

454 diaObjects[outCol] = np.nan

455

456 def _computeMeanPos(df):

457 aveCoord = geom.averageSpherePoint(

458 list(geom.SpherePoint(src["ra"], src["dec"], geom.degrees)

459 for idx, src in df.iterrows()))

460

461 return pd.Series({"ra": aveCoord.getRa().asDegrees(),

462 "dec": aveCoord.getDec().asDegrees()})

463

464 ans = diaSources.apply(_computeMeanPos)

465 diaObjects.loc[:, ["ra", "dec"]] = ans

466

467

468class HTMIndexDiaPositionConfig(DiaObjectCalculationPluginConfig):

469

470 htmLevel = pexConfig.Field(

471 dtype=int,

472 doc="Level of the HTM pixelization.",

473 default=20,

474 )

475

476

477@register("ap_HTMIndex")

478class HTMIndexDiaPosition(DiaObjectCalculationPlugin):

479 """Compute the mean position of a DiaObject given a set of DiaSources.

480

481 Notes

482 -----

483 This plugin was implemented to satisfy requirements of old APDB interface

484 which required ``pixelId`` column in DiaObject with HTM20 index. APDB

485 interface had migrated to not need that information, but we keep this

486 plugin in case it may be useful for something else.

487 """

488 ConfigClass = HTMIndexDiaPositionConfig

489

490 plugType = 'single'

491

492 inputCols = ["ra", "dec"]

493 outputCols = ["pixelId"]

494 needsFilter = False

495

496 def __init__(self, config, name, metadata):

497 DiaObjectCalculationPlugin.__init__(self, config, name, metadata)

498 self.pixelator = sphgeom.HtmPixelization(self.config.htmLevel)

499

500 @classmethod

501 def getExecutionOrder(cls):

502 return cls.FLUX_MOMENTS_CALCULATED

503

504 def calculate(self, diaObjects, diaObjectId, **kwargs):

505 """Compute the mean position of a DiaObject given a set of DiaSources

506

507 Parameters

508 ----------

509 diaObjects : `pandas.dataFrame`

510 Summary objects to store values in and read ra/dec from.

511 diaObjectId : `int`

512 Id of the diaObject to update.

513 **kwargs

514 Any additional keyword arguments that may be passed to the plugin.

515 """

516 sphPoint = geom.SpherePoint(

517 diaObjects.at[diaObjectId, "ra"] * geom.degrees,

518 diaObjects.at[diaObjectId, "dec"] * geom.degrees)

519 diaObjects.at[diaObjectId, "pixelId"] = self.pixelator.index(

520 sphPoint.getVector())

521

522

523class NumDiaSourcesDiaPluginConfig(DiaObjectCalculationPluginConfig):

524 pass

525

526

527@register("ap_nDiaSources")

528class NumDiaSourcesDiaPlugin(DiaObjectCalculationPlugin):

529 """Compute the total number of DiaSources associated with this DiaObject.

530 """

531

532 ConfigClass = NumDiaSourcesDiaPluginConfig

533 outputCols = ["nDiaSources"]

534 plugType = "multi"

535 needsFilter = False

536

537 @classmethod

538 def getExecutionOrder(cls):

539 return cls.DEFAULT_CATALOGCALCULATION

540

541 def calculate(self, diaObjects, diaSources, **kwargs):

542 """Compute the total number of DiaSources associated with this DiaObject.

543

544 Parameters

545 ----------

546 diaObject : `dict`

547 Summary object to store values in and read ra/dec from.

548 **kwargs

549 Any additional keyword arguments that may be passed to the plugin.

550 """

551 dtype = diaObjects["nDiaSources"].dtype

552 diaObjects.loc[:, "nDiaSources"] = diaSources.diaObjectId.count().astype(dtype)

553

554

555class SimpleSourceFlagDiaPluginConfig(DiaObjectCalculationPluginConfig):

556 pass

557

558

559@register("ap_diaObjectFlag")

560class SimpleSourceFlagDiaPlugin(DiaObjectCalculationPlugin):

561 """Find if any DiaSource is flagged.

562

563 Set the DiaObject flag if any DiaSource is flagged.

564 """

565

566 ConfigClass = NumDiaSourcesDiaPluginConfig

567 outputCols = ["flags"]

568 plugType = "multi"

569 needsFilter = False

570

571 @classmethod

572 def getExecutionOrder(cls):

573 return cls.DEFAULT_CATALOGCALCULATION

574

575 def calculate(self, diaObjects, diaSources, **kwargs):

576 """Find if any DiaSource is flagged.

577

578 Set the DiaObject flag if any DiaSource is flagged.

579

580 Parameters

581 ----------

582 diaObject : `dict`

583 Summary object to store values in and read ra/dec from.

584 **kwargs

585 Any additional keyword arguments that may be passed to the plugin.

586 """

587 dtype = diaObjects["flags"].dtype

588 diaObjects.loc[:, "flags"] = diaSources.flags.any().astype(dtype)

589

590

591class WeightedMeanDiaPsfFluxConfig(DiaObjectCalculationPluginConfig):

592 pass

593

594

595@register("ap_meanFlux")

596class WeightedMeanDiaPsfFlux(DiaObjectCalculationPlugin):

597 """Compute the weighted mean and mean error on the point source fluxes

598 of the DiaSource measured on the difference image.

599

600 Additionally store number of usable data points.

601 """

602

603 ConfigClass = WeightedMeanDiaPsfFluxConfig

604 outputCols = ["psfFluxMean", "psfFluxMeanErr", "psfFluxNdata"]

605 plugType = "multi"

606 needsFilter = True

607

608 @classmethod

609 def getExecutionOrder(cls):

610 return cls.DEFAULT_CATALOGCALCULATION

611

612 @catchWarnings(warns=["invalid value encountered",

613 "divide by zero"])

614 def calculate(self,

615 diaObjects,

616 diaSources,

617 filterDiaSources,

618 band,

619 **kwargs):

620 """Compute the weighted mean and mean error of the point source flux.

621

622 Parameters

623 ----------

624 diaObject : `dict`

625 Summary object to store values in.

626 diaSources : `pandas.DataFrame`

627 DataFrame representing all diaSources associated with this

628 diaObject.

629 filterDiaSources : `pandas.DataFrame`

630 DataFrame representing diaSources associated with this

631 diaObject that are observed in the band pass ``band``.

632 band : `str`

633 Simple, string name of the filter for the flux being calculated.

634 **kwargs

635 Any additional keyword arguments that may be passed to the plugin.

636 """

637 meanName = "{}_psfFluxMean".format(band)

638 errName = "{}_psfFluxMeanErr".format(band)

639 nDataName = "{}_psfFluxNdata".format(band)

640 if meanName not in diaObjects.columns:

641 diaObjects[meanName] = np.nan

642 if errName not in diaObjects.columns:

643 diaObjects[errName] = np.nan

644 if nDataName not in diaObjects.columns:

645 diaObjects[nDataName] = 0

646

647 def _weightedMean(df):

648 tmpDf = df[~np.logical_or(np.isnan(df["psfFlux"]),

649 np.isnan(df["psfFluxErr"]))]

650 tot_weight = np.nansum(1 / tmpDf["psfFluxErr"] ** 2)

651 fluxMean = np.nansum(tmpDf["psfFlux"]

652 / tmpDf["psfFluxErr"] ** 2)

653 fluxMean /= tot_weight

654 if tot_weight > 0:

655 fluxMeanErr = np.sqrt(1 / tot_weight)

656 else:

657 fluxMeanErr = np.nan

658 nFluxData = len(tmpDf)

659

660 return pd.Series({meanName: fluxMean,

661 errName: fluxMeanErr,

662 nDataName: nFluxData},

663 dtype="object")

664 df = filterDiaSources.apply(_weightedMean).astype(diaObjects.dtypes[[meanName, errName, nDataName]])

665

666 diaObjects.loc[:, [meanName, errName, nDataName]] = df

667

668

669class PercentileDiaPsfFluxConfig(DiaObjectCalculationPluginConfig):

670 percentiles = pexConfig.ListField(

671 dtype=int,

672 default=[5, 25, 50, 75, 95],

673 doc="Percentiles to calculate to compute values for. Should be "

674 "integer values."

675 )

676

677

678@register("ap_percentileFlux")

679class PercentileDiaPsfFlux(DiaObjectCalculationPlugin):

680 """Compute percentiles of diaSource fluxes.

681 """

682

683 ConfigClass = PercentileDiaPsfFluxConfig

684 # Output columns are created upon instantiation of the class.

685 outputCols = []

686 plugType = "multi"

687 needsFilter = True

688

689 def __init__(self, config, name, metadata, **kwargs):

690 DiaObjectCalculationPlugin.__init__(self,

691 config,

692 name,

693 metadata,

694 **kwargs)

695 self.outputCols = ["psfFluxPercentile{:02d}".format(percent)

696 for percent in self.config.percentiles]

697

698 @classmethod

699 def getExecutionOrder(cls):

700 return cls.DEFAULT_CATALOGCALCULATION

701

702 @catchWarnings(warns=["All-NaN slice encountered"])

703 def calculate(self,

704 diaObjects,

705 diaSources,

706 filterDiaSources,

707 band,

708 **kwargs):

709 """Compute the percentile fluxes of the point source flux.

710

711 Parameters

712 ----------

713 diaObject : `dict`

714 Summary object to store values in.

715 diaSources : `pandas.DataFrame`

716 DataFrame representing all diaSources associated with this

717 diaObject.

718 filterDiaSources : `pandas.DataFrame`

719 DataFrame representing diaSources associated with this

720 diaObject that are observed in the band pass ``band``.

721 band : `str`

722 Simple, string name of the filter for the flux being calculated.

723 **kwargs

724 Any additional keyword arguments that may be passed to the plugin.

725 """

726 pTileNames = []

727 dtype = None

728 for tilePercent in self.config.percentiles:

729 pTileName = "{}_psfFluxPercentile{:02d}".format(band,

730 tilePercent)

731 pTileNames.append(pTileName)

732 if pTileName not in diaObjects.columns:

733 diaObjects[pTileName] = np.nan

734 elif dtype is None:

735 dtype = diaObjects[pTileName].dtype

736

737 def _fluxPercentiles(df):

738 pTiles = np.nanpercentile(df["psfFlux"], self.config.percentiles)

739 return pd.Series(

740 dict((tileName, pTile)

741 for tileName, pTile in zip(pTileNames, pTiles)))

742 if dtype is None:

743 diaObjects.loc[:, pTileNames] = filterDiaSources.apply(_fluxPercentiles)

744 else:

745 diaObjects.loc[:, pTileNames] = filterDiaSources.apply(_fluxPercentiles).astype(dtype)

746

747

748class SigmaDiaPsfFluxConfig(DiaObjectCalculationPluginConfig):

749 pass

750

751

752@register("ap_sigmaFlux")

753class SigmaDiaPsfFlux(DiaObjectCalculationPlugin):

754 """Compute scatter of diaSource fluxes.

755 """

756

757 ConfigClass = SigmaDiaPsfFluxConfig

758 # Output columns are created upon instantiation of the class.

759 outputCols = ["psfFluxSigma"]

760 plugType = "multi"

761 needsFilter = True

762

763 @classmethod

764 def getExecutionOrder(cls):

765 return cls.DEFAULT_CATALOGCALCULATION

766

767 def calculate(self,

768 diaObjects,

769 diaSources,

770 filterDiaSources,

771 band,

772 **kwargs):

773 """Compute the sigma fluxes of the point source flux.

774

775 Parameters

776 ----------

777 diaObject : `dict`

778 Summary object to store values in.

779 diaSources : `pandas.DataFrame`

780 DataFrame representing all diaSources associated with this

781 diaObject.

782 filterDiaSources : `pandas.DataFrame`

783 DataFrame representing diaSources associated with this

784 diaObject that are observed in the band pass ``band``.

785 band : `str`

786 Simple, string name of the filter for the flux being calculated.

787 **kwargs

788 Any additional keyword arguments that may be passed to the plugin.

789 """

790 # Set "delta degrees of freedom (ddf)" to 1 to calculate the unbiased

791 # estimator of scatter (i.e. 'N - 1' instead of 'N').

792 column = "{}_psfFluxSigma".format(band)

793 if column in diaObjects:

794 dtype = diaObjects[column].dtype

795 diaObjects.loc[:, column] = filterDiaSources.psfFlux.std().astype(dtype)

796 else:

797 diaObjects.loc[:, column] = filterDiaSources.psfFlux.std()

798

799

800class Chi2DiaPsfFluxConfig(DiaObjectCalculationPluginConfig):

801 pass

802

803

804@register("ap_chi2Flux")

805class Chi2DiaPsfFlux(DiaObjectCalculationPlugin):

806 """Compute chi2 of diaSource fluxes.

807 """

808

809 ConfigClass = Chi2DiaPsfFluxConfig

810

811 # Required input Cols

812 inputCols = ["psfFluxMean"]

813 # Output columns are created upon instantiation of the class.

814 outputCols = ["psfFluxChi2"]

815 plugType = "multi"

816 needsFilter = True

817

818 @classmethod

819 def getExecutionOrder(cls):

820 return cls.FLUX_MOMENTS_CALCULATED

821

822 @catchWarnings(warns=["All-NaN slice encountered"])

823 def calculate(self,

824 diaObjects,

825 diaSources,

826 filterDiaSources,

827 band,

828 **kwargs):

829 """Compute the chi2 of the point source fluxes.

830

831 Parameters

832 ----------

833 diaObject : `dict`

834 Summary object to store values in.

835 diaSources : `pandas.DataFrame`

836 DataFrame representing all diaSources associated with this

837 diaObject.

838 filterDiaSources : `pandas.DataFrame`

839 DataFrame representing diaSources associated with this

840 diaObject that are observed in the band pass ``band``.

841 band : `str`

842 Simple, string name of the filter for the flux being calculated.

843 **kwargs

844 Any additional keyword arguments that may be passed to the plugin.

845 """

846 meanName = "{}_psfFluxMean".format(band)

847 column = "{}_psfFluxChi2".format(band)

848

849 def _chi2(df):

850 delta = (df["psfFlux"]

851 - diaObjects.at[df.diaObjectId.iat[0], meanName])

852 return np.nansum((delta / df["psfFluxErr"]) ** 2)

853

854 if column in diaObjects:

855 dtype = diaObjects[column].dtype

856 diaObjects.loc[:, column] = filterDiaSources.apply(_chi2).astype(dtype)

857 else:

858 diaObjects.loc[:, column] = filterDiaSources.apply(_chi2)

859

860

861class MadDiaPsfFluxConfig(DiaObjectCalculationPluginConfig):

862 pass

863

864

865@register("ap_madFlux")

866class MadDiaPsfFlux(DiaObjectCalculationPlugin):

867 """Compute median absolute deviation of diaSource fluxes.

868 """

869

870 ConfigClass = MadDiaPsfFluxConfig

871

872 # Required input Cols

873 # Output columns are created upon instantiation of the class.

874 outputCols = ["psfFluxMAD"]

875 plugType = "multi"

876 needsFilter = True

877

878 @classmethod

879 def getExecutionOrder(cls):

880 return cls.DEFAULT_CATALOGCALCULATION

881

882 @catchWarnings(warns=["All-NaN slice encountered"])

883 def calculate(self,

884 diaObjects,

885 diaSources,

886 filterDiaSources,

887 band,

888 **kwargs):

889 """Compute the median absolute deviation of the point source fluxes.

890

891 Parameters

892 ----------

893 diaObject : `dict`

894 Summary object to store values in.

895 diaSources : `pandas.DataFrame`

896 DataFrame representing all diaSources associated with this

897 diaObject.

898 filterDiaSources : `pandas.DataFrame`

899 DataFrame representing diaSources associated with this

900 diaObject that are observed in the band pass ``band``.

901 band : `str`

902 Simple, string name of the filter for the flux being calculated.

903 **kwargs

904 Any additional keyword arguments that may be passed to the plugin.

905 """

906 column = "{}_psfFluxMAD".format(band)

907 if column in diaObjects:

908 dtype = diaObjects[column].dtype

909 diaObjects.loc[:, column] = \

910 filterDiaSources.psfFlux.apply(median_absolute_deviation,

911 ignore_nan=True).astype(dtype)

912 else:

913 diaObjects.loc[:, column] = \

914 filterDiaSources.psfFlux.apply(median_absolute_deviation,

915 ignore_nan=True)

916

917

918class SkewDiaPsfFluxConfig(DiaObjectCalculationPluginConfig):

919 pass

920

921

922@register("ap_skewFlux")

923class SkewDiaPsfFlux(DiaObjectCalculationPlugin):

924 """Compute the skew of diaSource fluxes.

925 """

926

927 ConfigClass = SkewDiaPsfFluxConfig

928

929 # Required input Cols

930 # Output columns are created upon instantiation of the class.

931 outputCols = ["psfFluxSkew"]

932 plugType = "multi"

933 needsFilter = True

934

935 @classmethod

936 def getExecutionOrder(cls):

937 return cls.DEFAULT_CATALOGCALCULATION

938

939 def calculate(self,

940 diaObjects,

941 diaSources,

942 filterDiaSources,

943 band,

944 **kwargs):

945 """Compute the skew of the point source fluxes.

946

947 Parameters

948 ----------

949 diaObject : `dict`

950 Summary object to store values in.

951 diaSources : `pandas.DataFrame`

952 DataFrame representing all diaSources associated with this

953 diaObject.

954 filterDiaSources : `pandas.DataFrame`

955 DataFrame representing diaSources associated with this

956 diaObject that are observed in the band pass ``band``.

957 band : `str`

958 Simple, string name of the filter for the flux being calculated.

959 **kwargs

960 Any additional keyword arguments that may be passed to the plugin.

961 """

962 column = "{}_psfFluxSkew".format(band)

963 if column in diaObjects:

964 dtype = diaObjects[column].dtype

965 diaObjects.loc[:, column] = filterDiaSources.psfFlux.skew().astype(dtype)

966 else:

967 diaObjects.loc[:, column] = filterDiaSources.psfFlux.skew()

968

969

970class MinMaxDiaPsfFluxConfig(DiaObjectCalculationPluginConfig):

971 pass

972

973

974@register("ap_minMaxFlux")

975class MinMaxDiaPsfFlux(DiaObjectCalculationPlugin):

976 """Compute min/max of diaSource fluxes.

977 """

978

979 ConfigClass = MinMaxDiaPsfFluxConfig

980

981 # Required input Cols

982 # Output columns are created upon instantiation of the class.

983 outputCols = ["psfFluxMin", "psfFluxMax"]

984 plugType = "multi"

985 needsFilter = True

986

987 @classmethod

988 def getExecutionOrder(cls):

989 return cls.DEFAULT_CATALOGCALCULATION

990

991 def calculate(self,

992 diaObjects,

993 diaSources,

994 filterDiaSources,

995 band,

996 **kwargs):

997 """Compute min/max of the point source fluxes.

998

999 Parameters

1000 ----------

1001 diaObject : `dict`

1002 Summary object to store values in.

1003 diaSources : `pandas.DataFrame`

1004 DataFrame representing all diaSources associated with this

1005 diaObject.

1006 filterDiaSources : `pandas.DataFrame`

1007 DataFrame representing diaSources associated with this

1008 diaObject that are observed in the band pass ``band``.

1009 band : `str`

1010 Simple, string name of the filter for the flux being calculated.

1011 **kwargs

1012 Any additional keyword arguments that may be passed to the plugin.

1013 """

1014 minName = "{}_psfFluxMin".format(band)

1015 if minName not in diaObjects.columns:

1016 diaObjects[minName] = np.nan

1017 maxName = "{}_psfFluxMax".format(band)

1018 if maxName not in diaObjects.columns:

1019 diaObjects[maxName] = np.nan

1020

1021 dtype = diaObjects[minName].dtype

1022 diaObjects.loc[:, minName] = filterDiaSources.psfFlux.min().astype(dtype)

1023 diaObjects.loc[:, maxName] = filterDiaSources.psfFlux.max().astype(dtype)

1024

1025

1026class MaxSlopeDiaPsfFluxConfig(DiaObjectCalculationPluginConfig):

1027 pass

1028

1029

1030@register("ap_maxSlopeFlux")

1031class MaxSlopeDiaPsfFlux(DiaObjectCalculationPlugin):

1032 """Compute the maximum ratio time ordered deltaFlux / deltaTime.

1033 """

1034

1035 ConfigClass = MinMaxDiaPsfFluxConfig

1036

1037 # Required input Cols

1038 # Output columns are created upon instantiation of the class.

1039 outputCols = ["psfFluxMaxSlope"]

1040 plugType = "multi"

1041 needsFilter = True

1042

1043 @classmethod

1044 def getExecutionOrder(cls):

1045 return cls.DEFAULT_CATALOGCALCULATION

1046

1047 def calculate(self,

1048 diaObjects,

1049 diaSources,

1050 filterDiaSources,

1051 band,

1052 **kwargs):

1053 """Compute the maximum ratio time ordered deltaFlux / deltaTime.

1054

1055 Parameters

1056 ----------

1057 diaObject : `dict`

1058 Summary object to store values in.

1059 diaSources : `pandas.DataFrame`

1060 DataFrame representing all diaSources associated with this

1061 diaObject.

1062 filterDiaSources : `pandas.DataFrame`

1063 DataFrame representing diaSources associated with this

1064 diaObject that are observed in the band pass ``band``.

1065 band : `str`

1066 Simple, string name of the filter for the flux being calculated.

1067 **kwargs

1068 Any additional keyword arguments that may be passed to the plugin.

1069 """

1070

1071 def _maxSlope(df):

1072 tmpDf = df[~np.logical_or(np.isnan(df["psfFlux"]),

1073 np.isnan(df["midpointMjdTai"]))]

1074 if len(tmpDf) < 2:

1075 return np.nan

1076 times = tmpDf["midpointMjdTai"].to_numpy()

1077 timeArgs = times.argsort()

1078 times = times[timeArgs]

1079 fluxes = tmpDf["psfFlux"].to_numpy()[timeArgs]

1080 return (np.diff(fluxes) / np.diff(times)).max()

1081

1082 column = "{}_psfFluxMaxSlope".format(band)

1083 if column in diaObjects:

1084 dtype = diaObjects[column].dtype

1085 diaObjects.loc[:, column] = filterDiaSources.apply(_maxSlope).astype(dtype)

1086 else:

1087 diaObjects.loc[:, column] = filterDiaSources.apply(_maxSlope)

1088

1089

1090class ErrMeanDiaPsfFluxConfig(DiaObjectCalculationPluginConfig):

1091 pass

1092

1093

1094@register("ap_meanErrFlux")

1095class ErrMeanDiaPsfFlux(DiaObjectCalculationPlugin):

1096 """Compute the mean of the dia source errors.

1097 """

1098

1099 ConfigClass = ErrMeanDiaPsfFluxConfig

1100

1101 # Required input Cols

1102 # Output columns are created upon instantiation of the class.

1103 outputCols = ["psfFluxErrMean"]

1104 plugType = "multi"

1105 needsFilter = True

1106

1107 @classmethod

1108 def getExecutionOrder(cls):

1109 return cls.DEFAULT_CATALOGCALCULATION

1110

1111 def calculate(self,

1112 diaObjects,

1113 diaSources,

1114 filterDiaSources,

1115 band,

1116 **kwargs):

1117 """Compute the mean of the dia source errors.

1118

1119 Parameters

1120 ----------

1121 diaObject : `dict`

1122 Summary object to store values in.

1123 diaSources : `pandas.DataFrame`

1124 DataFrame representing all diaSources associated with this

1125 diaObject.

1126 filterDiaSources : `pandas.DataFrame`

1127 DataFrame representing diaSources associated with this

1128 diaObject that are observed in the band pass ``band``.

1129 band : `str`

1130 Simple, string name of the filter for the flux being calculated.

1131 **kwargs

1132 Any additional keyword arguments that may be passed to the plugin.

1133 """

1134 column = "{}_psfFluxErrMean".format(band)

1135 if column in diaObjects:

1136 dtype = diaObjects[column].dtype

1137 diaObjects.loc[:, column] = filterDiaSources.psfFluxErr.mean().astype(dtype)

1138 else:

1139 diaObjects.loc[:, column] = filterDiaSources.psfFluxErr.mean()

1140

1141

1142class LinearFitDiaPsfFluxConfig(DiaObjectCalculationPluginConfig):

1143 pass

1144

1145

1146@register("ap_linearFit")

1147class LinearFitDiaPsfFlux(DiaObjectCalculationPlugin):

1148 """Compute fit a linear model to flux vs time.

1149 """

1150

1151 ConfigClass = LinearFitDiaPsfFluxConfig

1152

1153 # Required input Cols

1154 # Output columns are created upon instantiation of the class.

1155 outputCols = ["psfFluxLinearSlope", "psfFluxLinearIntercept"]

1156 plugType = "multi"

1157 needsFilter = True

1158

1159 @classmethod

1160 def getExecutionOrder(cls):

1161 return cls.DEFAULT_CATALOGCALCULATION

1162

1163 def calculate(self,

1164 diaObjects,

1165 diaSources,

1166 filterDiaSources,

1167 band,

1168 **kwargs):

1169 """Compute fit a linear model to flux vs time.

1170

1171 Parameters

1172 ----------

1173 diaObject : `dict`

1174 Summary object to store values in.

1175 diaSources : `pandas.DataFrame`

1176 DataFrame representing all diaSources associated with this

1177 diaObject.

1178 filterDiaSources : `pandas.DataFrame`

1179 DataFrame representing diaSources associated with this

1180 diaObject that are observed in the band pass ``band``.

1181 band : `str`

1182 Simple, string name of the filter for the flux being calculated.

1183 **kwargs

1184 Any additional keyword arguments that may be passed to the plugin.

1185 """

1186

1187 mName = "{}_psfFluxLinearSlope".format(band)

1188 if mName not in diaObjects.columns:

1189 diaObjects[mName] = np.nan

1190 bName = "{}_psfFluxLinearIntercept".format(band)

1191 if bName not in diaObjects.columns:

1192 diaObjects[bName] = np.nan

1193 dtype = diaObjects[mName].dtype

1194

1195 def _linearFit(df):

1196 tmpDf = df[~np.logical_or(

1197 np.isnan(df["psfFlux"]),

1198 np.logical_or(np.isnan(df["psfFluxErr"]),

1199 np.isnan(df["midpointMjdTai"])))]

1200 if len(tmpDf) < 2:

1201 return pd.Series({mName: np.nan, bName: np.nan})

1202 fluxes = tmpDf["psfFlux"].to_numpy()

1203 errors = tmpDf["psfFluxErr"].to_numpy()

1204 times = tmpDf["midpointMjdTai"].to_numpy()

1205 A = np.array([times / errors, 1 / errors]).transpose()

1206 m, b = lsq_linear(A, fluxes / errors).x

1207 return pd.Series({mName: m, bName: b}, dtype=dtype)

1208

1209 diaObjects.loc[:, [mName, bName]] = filterDiaSources.apply(_linearFit)

1210

1211

1212class StetsonJDiaPsfFluxConfig(DiaObjectCalculationPluginConfig):

1213 pass

1214

1215

1216@register("ap_stetsonJ")

1217class StetsonJDiaPsfFlux(DiaObjectCalculationPlugin):

1218 """Compute the StetsonJ statistic on the DIA point source fluxes.

1219 """

1220

1221 ConfigClass = LinearFitDiaPsfFluxConfig

1222

1223 # Required input Cols

1224 inputCols = ["psfFluxMean"]

1225 # Output columns are created upon instantiation of the class.

1226 outputCols = ["psfFluxStetsonJ"]

1227 plugType = "multi"

1228 needsFilter = True

1229

1230 @classmethod

1231 def getExecutionOrder(cls):

1232 return cls.FLUX_MOMENTS_CALCULATED

1233

1234 def calculate(self,

1235 diaObjects,

1236 diaSources,

1237 filterDiaSources,

1238 band,

1239 **kwargs):

1240 """Compute the StetsonJ statistic on the DIA point source fluxes.

1241

1242 Parameters

1243 ----------

1244 diaObject : `dict`

1245 Summary object to store values in.

1246 diaSources : `pandas.DataFrame`

1247 DataFrame representing all diaSources associated with this

1248 diaObject.

1249 filterDiaSources : `pandas.DataFrame`

1250 DataFrame representing diaSources associated with this

1251 diaObject that are observed in the band pass ``band``.

1252 band : `str`

1253 Simple, string name of the filter for the flux being calculated.

1254 **kwargs

1255 Any additional keyword arguments that may be passed to the plugin.

1256 """

1257 meanName = "{}_psfFluxMean".format(band)

1258

1259 def _stetsonJ(df):

1260 tmpDf = df[~np.logical_or(np.isnan(df["psfFlux"]),

1261 np.isnan(df["psfFluxErr"]))]

1262 if len(tmpDf) < 2:

1263 return np.nan

1264 fluxes = tmpDf["psfFlux"].to_numpy()

1265 errors = tmpDf["psfFluxErr"].to_numpy()

1266

1267 return self._stetson_J(

1268 fluxes,

1269 errors,

1270 diaObjects.at[tmpDf.diaObjectId.iat[0], meanName])

1271

1272 column = "{}_psfFluxStetsonJ".format(band)

1273 if column in diaObjects:

1274 dtype = diaObjects[column].dtype

1275 diaObjects.loc[:, column] = filterDiaSources.apply(_stetsonJ).astype(dtype)

1276 else:

1277 diaObjects.loc[:, column] = filterDiaSources.apply(_stetsonJ)

1278

1279 def _stetson_J(self, fluxes, errors, mean=None):

1280 """Compute the single band stetsonJ statistic.

1281

1282 Parameters

1283 ----------

1284 fluxes : `numpy.ndarray` (N,)

1285 Calibrated lightcurve flux values.

1286 errors : `numpy.ndarray` (N,)

1287 Errors on the calibrated lightcurve fluxes.

1288 mean : `float`

1289 Starting mean from previous plugin.

1290

1291 Returns

1292 -------

1293 stetsonJ : `float`

1294 stetsonJ statistic for the input fluxes and errors.

1295

1296 References

1297 ----------

1298 .. [1] Stetson, P. B., "On the Automatic Determination of Light-Curve

1299 Parameters for Cepheid Variables", PASP, 108, 851S, 1996

1300 """

1301 n_points = len(fluxes)

1302 flux_mean = self._stetson_mean(fluxes, errors, mean)

1303 delta_val = (

1304 np.sqrt(n_points / (n_points - 1)) * (fluxes - flux_mean) / errors)

1305 p_k = delta_val ** 2 - 1

1306

1307 return np.mean(np.sign(p_k) * np.sqrt(np.fabs(p_k)))

1308

1309 def _stetson_mean(self,

1310 values,

1311 errors,

1312 mean=None,

1313 alpha=2.,

1314 beta=2.,

1315 n_iter=20,

1316 tol=1e-6):

1317 """Compute the stetson mean of the fluxes which down-weights outliers.

1318

1319 Weighted biased on an error weighted difference scaled by a constant

1320 (1/``a``) and raised to the power beta. Higher betas more harshly

1321 penalize outliers and ``a`` sets the number of sigma where a weighted

1322 difference of 1 occurs.

1323

1324 Parameters

1325 ----------

1326 values : `numpy.dnarray`, (N,)

1327 Input values to compute the mean of.

1328 errors : `numpy.ndarray`, (N,)

1329 Errors on the input values.

1330 mean : `float`

1331 Starting mean value or None.

1332 alpha : `float`

1333 Scalar down-weighting of the fractional difference. lower->more

1334 clipping. (Default value is 2.)

1335 beta : `float`

1336 Power law slope of the used to down-weight outliers. higher->more

1337 clipping. (Default value is 2.)

1338 n_iter : `int`

1339 Number of iterations of clipping.

1340 tol : `float`

1341 Fractional and absolute tolerance goal on the change in the mean

1342 before exiting early. (Default value is 1e-6)

1343

1344 Returns

1345 -------

1346 mean : `float`

1347 Weighted stetson mean result.

1348

1349 References

1350 ----------

1351 .. [1] Stetson, P. B., "On the Automatic Determination of Light-Curve

1352 Parameters for Cepheid Variables", PASP, 108, 851S, 1996

1353 """

1354 n_points = len(values)

1355 n_factor = np.sqrt(n_points / (n_points - 1))

1356 inv_var = 1 / errors ** 2

1357

1358 if mean is None:

1359 mean = np.average(values, weights=inv_var)

1360 for iter_idx in range(n_iter):

1361 chi = np.fabs(n_factor * (values - mean) / errors)

1362 tmp_mean = np.average(

1363 values,

1364 weights=inv_var / (1 + (chi / alpha) ** beta))

1365 diff = np.fabs(tmp_mean - mean)

1366 mean = tmp_mean

1367 if diff / mean < tol and diff < tol:

1368 break

1369 return mean

1370

1371

1372class WeightedMeanDiaTotFluxConfig(DiaObjectCalculationPluginConfig):

1373 pass

1374

1375

1376@register("ap_meanTotFlux")

1377class WeightedMeanDiaTotFlux(DiaObjectCalculationPlugin):

1378 """Compute the weighted mean and mean error on the point source fluxes

1379 forced photometered at the DiaSource location in the calibrated image.

1380

1381 Additionally store number of usable data points.

1382 """

1383

1384 ConfigClass = WeightedMeanDiaPsfFluxConfig

1385 outputCols = ["scienceFluxMean", "scienceFluxMeanErr"]

1386 plugType = "multi"

1387 needsFilter = True

1388

1389 @classmethod

1390 def getExecutionOrder(cls):

1391 return cls.DEFAULT_CATALOGCALCULATION

1392

1393 @catchWarnings(warns=["invalid value encountered",

1394 "divide by zero"])

1395 def calculate(self,

1396 diaObjects,

1397 diaSources,

1398 filterDiaSources,

1399 band,

1400 **kwargs):

1401 """Compute the weighted mean and mean error of the point source flux.

1402

1403 Parameters

1404 ----------

1405 diaObject : `dict`

1406 Summary object to store values in.

1407 diaSources : `pandas.DataFrame`

1408 DataFrame representing all diaSources associated with this

1409 diaObject.

1410 filterDiaSources : `pandas.DataFrame`

1411 DataFrame representing diaSources associated with this

1412 diaObject that are observed in the band pass ``band``.

1413 band : `str`

1414 Simple, string name of the filter for the flux being calculated.

1415 **kwargs

1416 Any additional keyword arguments that may be passed to the plugin.

1417 """

1418 totMeanName = "{}_scienceFluxMean".format(band)

1419 if totMeanName not in diaObjects.columns:

1420 diaObjects[totMeanName] = np.nan

1421 totErrName = "{}_scienceFluxMeanErr".format(band)

1422 if totErrName not in diaObjects.columns:

1423 diaObjects[totErrName] = np.nan

1424

1425 def _meanFlux(df):

1426 tmpDf = df[~np.logical_or(np.isnan(df["scienceFlux"]),

1427 np.isnan(df["scienceFluxErr"]))]

1428 tot_weight = np.nansum(1 / tmpDf["scienceFluxErr"] ** 2)

1429 fluxMean = np.nansum(tmpDf["scienceFlux"]

1430 / tmpDf["scienceFluxErr"] ** 2)

1431 fluxMean /= tot_weight

1432 fluxMeanErr = np.sqrt(1 / tot_weight)

1433

1434 return pd.Series({totMeanName: fluxMean,

1435 totErrName: fluxMeanErr})

1436

1437 df = filterDiaSources.apply(_meanFlux).astype(diaObjects.dtypes[[totMeanName, totErrName]])

1438 diaObjects.loc[:, [totMeanName, totErrName]] = df

1439

1440

1441class SigmaDiaTotFluxConfig(DiaObjectCalculationPluginConfig):

1442 pass

1443

1444

1445@register("ap_sigmaTotFlux")

1446class SigmaDiaTotFlux(DiaObjectCalculationPlugin):

1447 """Compute scatter of diaSource fluxes.

1448 """

1449

1450 ConfigClass = SigmaDiaPsfFluxConfig

1451 # Output columns are created upon instantiation of the class.

1452 outputCols = ["scienceFluxSigma"]

1453 plugType = "multi"

1454 needsFilter = True

1455

1456 @classmethod

1457 def getExecutionOrder(cls):

1458 return cls.DEFAULT_CATALOGCALCULATION

1459

1460 def calculate(self,

1461 diaObjects,

1462 diaSources,

1463 filterDiaSources,

1464 band,

1465 **kwargs):

1466 """Compute the sigma fluxes of the point source flux measured on the

1467 calibrated image.

1468

1469 Parameters

1470 ----------

1471 diaObject : `dict`

1472 Summary object to store values in.

1473 diaSources : `pandas.DataFrame`

1474 DataFrame representing all diaSources associated with this

1475 diaObject.

1476 filterDiaSources : `pandas.DataFrame`

1477 DataFrame representing diaSources associated with this

1478 diaObject that are observed in the band pass ``band``.

1479 band : `str`

1480 Simple, string name of the filter for the flux being calculated.

1481 **kwargs

1482 Any additional keyword arguments that may be passed to the plugin.

1483 """

1484 # Set "delta degrees of freedom (ddf)" to 1 to calculate the unbiased

1485 # estimator of scatter (i.e. 'N - 1' instead of 'N').

1486 column = "{}_scienceFluxSigma".format(band)

1487 if column in diaObjects:

1488 dtype = diaObjects[column].dtype

1489 diaObjects.loc[:, column] = filterDiaSources.scienceFlux.std().astype(dtype)

1490 else:

1491

1492 diaObjects.loc[:, column] = filterDiaSources.scienceFlux.std()