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 # std() is inconsistent about what type it returns; force to double.

798 diaObjects.loc[:, column] = filterDiaSources.psfFlux.std().astype(np.float64)

799

800

801class Chi2DiaPsfFluxConfig(DiaObjectCalculationPluginConfig):

802 pass

803

804

805@register("ap_chi2Flux")

806class Chi2DiaPsfFlux(DiaObjectCalculationPlugin):

807 """Compute chi2 of diaSource fluxes.

808 """

809

810 ConfigClass = Chi2DiaPsfFluxConfig

811

812 # Required input Cols

813 inputCols = ["psfFluxMean"]

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

815 outputCols = ["psfFluxChi2"]

816 plugType = "multi"

817 needsFilter = True

818

819 @classmethod

820 def getExecutionOrder(cls):

821 return cls.FLUX_MOMENTS_CALCULATED

822

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

824 def calculate(self,

825 diaObjects,

826 diaSources,

827 filterDiaSources,

828 band,

829 **kwargs):

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

831

832 Parameters

833 ----------

834 diaObject : `dict`

835 Summary object to store values in.

836 diaSources : `pandas.DataFrame`

837 DataFrame representing all diaSources associated with this

838 diaObject.

839 filterDiaSources : `pandas.DataFrame`

840 DataFrame representing diaSources associated with this

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

842 band : `str`

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

844 **kwargs

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

846 """

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

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

849

850 def _chi2(df):

851 delta = (df["psfFlux"]

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

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

854

855 if column in diaObjects:

856 dtype = diaObjects[column].dtype

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

858 else:

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

860

861

862class MadDiaPsfFluxConfig(DiaObjectCalculationPluginConfig):

863 pass

864

865

866@register("ap_madFlux")

867class MadDiaPsfFlux(DiaObjectCalculationPlugin):

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

869 """

870

871 ConfigClass = MadDiaPsfFluxConfig

872

873 # Required input Cols

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

875 outputCols = ["psfFluxMAD"]

876 plugType = "multi"

877 needsFilter = True

878

879 @classmethod

880 def getExecutionOrder(cls):

881 return cls.DEFAULT_CATALOGCALCULATION

882

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

884 def calculate(self,

885 diaObjects,

886 diaSources,

887 filterDiaSources,

888 band,

889 **kwargs):

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

891

892 Parameters

893 ----------

894 diaObject : `dict`

895 Summary object to store values in.

896 diaSources : `pandas.DataFrame`

897 DataFrame representing all diaSources associated with this

898 diaObject.

899 filterDiaSources : `pandas.DataFrame`

900 DataFrame representing diaSources associated with this

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

902 band : `str`

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

904 **kwargs

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

906 """

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

908 if column in diaObjects:

909 dtype = diaObjects[column].dtype

910 diaObjects.loc[:, column] = \

911 filterDiaSources.psfFlux.apply(median_absolute_deviation,

912 ignore_nan=True).astype(dtype)

913 else:

914 diaObjects.loc[:, column] = \

915 filterDiaSources.psfFlux.apply(median_absolute_deviation,

916 ignore_nan=True)

917

918

919class SkewDiaPsfFluxConfig(DiaObjectCalculationPluginConfig):

920 pass

921

922

923@register("ap_skewFlux")

924class SkewDiaPsfFlux(DiaObjectCalculationPlugin):

925 """Compute the skew of diaSource fluxes.

926 """

927

928 ConfigClass = SkewDiaPsfFluxConfig

929

930 # Required input Cols

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

932 outputCols = ["psfFluxSkew"]

933 plugType = "multi"

934 needsFilter = True

935

936 @classmethod

937 def getExecutionOrder(cls):

938 return cls.DEFAULT_CATALOGCALCULATION

939

940 def calculate(self,

941 diaObjects,

942 diaSources,

943 filterDiaSources,

944 band,

945 **kwargs):

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

947

948 Parameters

949 ----------

950 diaObject : `dict`

951 Summary object to store values in.

952 diaSources : `pandas.DataFrame`

953 DataFrame representing all diaSources associated with this

954 diaObject.

955 filterDiaSources : `pandas.DataFrame`

956 DataFrame representing diaSources associated with this

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

958 band : `str`

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

960 **kwargs

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

962 """

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

964 if column in diaObjects:

965 dtype = diaObjects[column].dtype

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

967 else:

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

969

970

971class MinMaxDiaPsfFluxConfig(DiaObjectCalculationPluginConfig):

972 pass

973

974

975@register("ap_minMaxFlux")

976class MinMaxDiaPsfFlux(DiaObjectCalculationPlugin):

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

978 """

979

980 ConfigClass = MinMaxDiaPsfFluxConfig

981

982 # Required input Cols

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

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

985 plugType = "multi"

986 needsFilter = True

987

988 @classmethod

989 def getExecutionOrder(cls):

990 return cls.DEFAULT_CATALOGCALCULATION

991

992 def calculate(self,

993 diaObjects,

994 diaSources,

995 filterDiaSources,

996 band,

997 **kwargs):

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

999

1000 Parameters

1001 ----------

1002 diaObject : `dict`

1003 Summary object to store values in.

1004 diaSources : `pandas.DataFrame`

1005 DataFrame representing all diaSources associated with this

1006 diaObject.

1007 filterDiaSources : `pandas.DataFrame`

1008 DataFrame representing diaSources associated with this

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

1010 band : `str`

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

1012 **kwargs

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

1014 """

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

1016 if minName not in diaObjects.columns:

1017 diaObjects[minName] = np.nan

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

1019 if maxName not in diaObjects.columns:

1020 diaObjects[maxName] = np.nan

1021

1022 dtype = diaObjects[minName].dtype

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

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

1025

1026

1027class MaxSlopeDiaPsfFluxConfig(DiaObjectCalculationPluginConfig):

1028 pass

1029

1030

1031@register("ap_maxSlopeFlux")

1032class MaxSlopeDiaPsfFlux(DiaObjectCalculationPlugin):

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

1034 """

1035

1036 ConfigClass = MinMaxDiaPsfFluxConfig

1037

1038 # Required input Cols

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

1040 outputCols = ["psfFluxMaxSlope"]

1041 plugType = "multi"

1042 needsFilter = True

1043

1044 @classmethod

1045 def getExecutionOrder(cls):

1046 return cls.DEFAULT_CATALOGCALCULATION

1047

1048 def calculate(self,

1049 diaObjects,

1050 diaSources,

1051 filterDiaSources,

1052 band,

1053 **kwargs):

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

1055

1056 Parameters

1057 ----------

1058 diaObject : `dict`

1059 Summary object to store values in.

1060 diaSources : `pandas.DataFrame`

1061 DataFrame representing all diaSources associated with this

1062 diaObject.

1063 filterDiaSources : `pandas.DataFrame`

1064 DataFrame representing diaSources associated with this

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

1066 band : `str`

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

1068 **kwargs

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

1070 """

1071

1072 def _maxSlope(df):

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

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

1075 if len(tmpDf) < 2:

1076 return np.nan

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

1078 timeArgs = times.argsort()

1079 times = times[timeArgs]

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

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

1082

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

1084 if column in diaObjects:

1085 dtype = diaObjects[column].dtype

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

1087 else:

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

1089

1090

1091class ErrMeanDiaPsfFluxConfig(DiaObjectCalculationPluginConfig):

1092 pass

1093

1094

1095@register("ap_meanErrFlux")

1096class ErrMeanDiaPsfFlux(DiaObjectCalculationPlugin):

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

1098 """

1099

1100 ConfigClass = ErrMeanDiaPsfFluxConfig

1101

1102 # Required input Cols

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

1104 outputCols = ["psfFluxErrMean"]

1105 plugType = "multi"

1106 needsFilter = True

1107

1108 @classmethod

1109 def getExecutionOrder(cls):

1110 return cls.DEFAULT_CATALOGCALCULATION

1111

1112 def calculate(self,

1113 diaObjects,

1114 diaSources,

1115 filterDiaSources,

1116 band,

1117 **kwargs):

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

1119

1120 Parameters

1121 ----------

1122 diaObject : `dict`

1123 Summary object to store values in.

1124 diaSources : `pandas.DataFrame`

1125 DataFrame representing all diaSources associated with this

1126 diaObject.

1127 filterDiaSources : `pandas.DataFrame`

1128 DataFrame representing diaSources associated with this

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

1130 band : `str`

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

1132 **kwargs

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

1134 """

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

1136 if column in diaObjects:

1137 dtype = diaObjects[column].dtype

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

1139 else:

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

1141

1142

1143class LinearFitDiaPsfFluxConfig(DiaObjectCalculationPluginConfig):

1144 pass

1145

1146

1147@register("ap_linearFit")

1148class LinearFitDiaPsfFlux(DiaObjectCalculationPlugin):

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

1150 """

1151

1152 ConfigClass = LinearFitDiaPsfFluxConfig

1153

1154 # Required input Cols

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

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

1157 plugType = "multi"

1158 needsFilter = True

1159

1160 @classmethod

1161 def getExecutionOrder(cls):

1162 return cls.DEFAULT_CATALOGCALCULATION

1163

1164 def calculate(self,

1165 diaObjects,

1166 diaSources,

1167 filterDiaSources,

1168 band,

1169 **kwargs):

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

1171

1172 Parameters

1173 ----------

1174 diaObject : `dict`

1175 Summary object to store values in.

1176 diaSources : `pandas.DataFrame`

1177 DataFrame representing all diaSources associated with this

1178 diaObject.

1179 filterDiaSources : `pandas.DataFrame`

1180 DataFrame representing diaSources associated with this

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

1182 band : `str`

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

1184 **kwargs

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

1186 """

1187

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

1189 if mName not in diaObjects.columns:

1190 diaObjects[mName] = np.nan

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

1192 if bName not in diaObjects.columns:

1193 diaObjects[bName] = np.nan

1194 dtype = diaObjects[mName].dtype

1195

1196 def _linearFit(df):

1197 tmpDf = df[~np.logical_or(

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

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

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

1201 if len(tmpDf) < 2:

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

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

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

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

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

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

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

1209

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

1211

1212

1213class StetsonJDiaPsfFluxConfig(DiaObjectCalculationPluginConfig):

1214 pass

1215

1216

1217@register("ap_stetsonJ")

1218class StetsonJDiaPsfFlux(DiaObjectCalculationPlugin):

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

1220 """

1221

1222 ConfigClass = LinearFitDiaPsfFluxConfig

1223

1224 # Required input Cols

1225 inputCols = ["psfFluxMean"]

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

1227 outputCols = ["psfFluxStetsonJ"]

1228 plugType = "multi"

1229 needsFilter = True

1230

1231 @classmethod

1232 def getExecutionOrder(cls):

1233 return cls.FLUX_MOMENTS_CALCULATED

1234

1235 def calculate(self,

1236 diaObjects,

1237 diaSources,

1238 filterDiaSources,

1239 band,

1240 **kwargs):

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

1242

1243 Parameters

1244 ----------

1245 diaObject : `dict`

1246 Summary object to store values in.

1247 diaSources : `pandas.DataFrame`

1248 DataFrame representing all diaSources associated with this

1249 diaObject.

1250 filterDiaSources : `pandas.DataFrame`

1251 DataFrame representing diaSources associated with this

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

1253 band : `str`

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

1255 **kwargs

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

1257 """

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

1259

1260 def _stetsonJ(df):

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

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

1263 if len(tmpDf) < 2:

1264 return np.nan

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

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

1267

1268 return self._stetson_J(

1269 fluxes,

1270 errors,

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

1272

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

1274 if column in diaObjects:

1275 dtype = diaObjects[column].dtype

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

1277 else:

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

1279

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

1281 """Compute the single band stetsonJ statistic.

1282

1283 Parameters

1284 ----------

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

1286 Calibrated lightcurve flux values.

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

1288 Errors on the calibrated lightcurve fluxes.

1289 mean : `float`

1290 Starting mean from previous plugin.

1291

1292 Returns

1293 -------

1294 stetsonJ : `float`

1295 stetsonJ statistic for the input fluxes and errors.

1296

1297 References

1298 ----------

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

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

1301 """

1302 n_points = len(fluxes)

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

1304 delta_val = (

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

1306 p_k = delta_val ** 2 - 1

1307

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

1309

1310 def _stetson_mean(self,

1311 values,

1312 errors,

1313 mean=None,

1314 alpha=2.,

1315 beta=2.,

1316 n_iter=20,

1317 tol=1e-6):

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

1319

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

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

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

1323 difference of 1 occurs.

1324

1325 Parameters

1326 ----------

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

1328 Input values to compute the mean of.

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

1330 Errors on the input values.

1331 mean : `float`

1332 Starting mean value or None.

1333 alpha : `float`

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

1335 clipping. (Default value is 2.)

1336 beta : `float`

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

1338 clipping. (Default value is 2.)

1339 n_iter : `int`

1340 Number of iterations of clipping.

1341 tol : `float`

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

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

1344

1345 Returns

1346 -------

1347 mean : `float`

1348 Weighted stetson mean result.

1349

1350 References

1351 ----------

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

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

1354 """

1355 n_points = len(values)

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

1357 inv_var = 1 / errors ** 2

1358

1359 if mean is None:

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

1361 for iter_idx in range(n_iter):

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

1363 tmp_mean = np.average(

1364 values,

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

1366 diff = np.fabs(tmp_mean - mean)

1367 mean = tmp_mean

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

1369 break

1370 return mean

1371

1372

1373class WeightedMeanDiaTotFluxConfig(DiaObjectCalculationPluginConfig):

1374 pass

1375

1376

1377@register("ap_meanTotFlux")

1378class WeightedMeanDiaTotFlux(DiaObjectCalculationPlugin):

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

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

1381

1382 Additionally store number of usable data points.

1383 """

1384

1385 ConfigClass = WeightedMeanDiaPsfFluxConfig

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

1387 plugType = "multi"

1388 needsFilter = True

1389

1390 @classmethod

1391 def getExecutionOrder(cls):

1392 return cls.DEFAULT_CATALOGCALCULATION

1393

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

1395 "divide by zero"])

1396 def calculate(self,

1397 diaObjects,

1398 diaSources,

1399 filterDiaSources,

1400 band,

1401 **kwargs):

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

1403

1404 Parameters

1405 ----------

1406 diaObject : `dict`

1407 Summary object to store values in.

1408 diaSources : `pandas.DataFrame`

1409 DataFrame representing all diaSources associated with this

1410 diaObject.

1411 filterDiaSources : `pandas.DataFrame`

1412 DataFrame representing diaSources associated with this

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

1414 band : `str`

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

1416 **kwargs

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

1418 """

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

1420 if totMeanName not in diaObjects.columns:

1421 diaObjects[totMeanName] = np.nan

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

1423 if totErrName not in diaObjects.columns:

1424 diaObjects[totErrName] = np.nan

1425

1426 def _meanFlux(df):

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

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

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

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

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

1432 fluxMean /= tot_weight

1433 fluxMeanErr = np.sqrt(1 / tot_weight)

1434

1435 return pd.Series({totMeanName: fluxMean,

1436 totErrName: fluxMeanErr})

1437

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

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

1440

1441

1442class SigmaDiaTotFluxConfig(DiaObjectCalculationPluginConfig):

1443 pass

1444

1445

1446@register("ap_sigmaTotFlux")

1447class SigmaDiaTotFlux(DiaObjectCalculationPlugin):

1448 """Compute scatter of diaSource fluxes.

1449 """

1450

1451 ConfigClass = SigmaDiaPsfFluxConfig

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

1453 outputCols = ["scienceFluxSigma"]

1454 plugType = "multi"

1455 needsFilter = True

1456

1457 @classmethod

1458 def getExecutionOrder(cls):

1459 return cls.DEFAULT_CATALOGCALCULATION

1460

1461 def calculate(self,

1462 diaObjects,

1463 diaSources,

1464 filterDiaSources,

1465 band,

1466 **kwargs):

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

1468 calibrated image.

1469

1470 Parameters

1471 ----------

1472 diaObject : `dict`

1473 Summary object to store values in.

1474 diaSources : `pandas.DataFrame`

1475 DataFrame representing all diaSources associated with this

1476 diaObject.

1477 filterDiaSources : `pandas.DataFrame`

1478 DataFrame representing diaSources associated with this

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

1480 band : `str`

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

1482 **kwargs

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

1484 """

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

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

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

1488 if column in diaObjects:

1489 dtype = diaObjects[column].dtype

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

1491 else:

1492

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