-
Notifications
You must be signed in to change notification settings - Fork 158
/
pldcorrector.py
569 lines (507 loc) · 22.2 KB
/
pldcorrector.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
"""Defines a `PLDCorrector` class which provides a simple way to correct a
light curve by utilizing the pixel time series data contained within the
target's own Target Pixel File.
`PLDCorrector` builds upon `RegressionCorrector` by correlating the light curve
against a design matrix composed of the following elements:
* A background light curve to capture the dominant scattered light systematics.
* Background-corrected pixel time series to capture any residual systematics.
* Splines to capture the target's intrinsic variability.
"""
import logging
import warnings
from itertools import combinations_with_replacement as multichoose
import numpy as np
import matplotlib.pyplot as plt
from astropy.utils.decorators import deprecated, deprecated_renamed_argument
from .designmatrix import (
DesignMatrix,
DesignMatrixCollection,
SparseDesignMatrixCollection,
)
from .regressioncorrector import RegressionCorrector
from .designmatrix import create_spline_matrix, create_sparse_spline_matrix
from .. import MPLSTYLE
from ..utils import LightkurveDeprecationWarning
log = logging.getLogger(__name__)
__all__ = ["PLDCorrector", "TessPLDCorrector"]
class PLDCorrector(RegressionCorrector):
r"""Implements the Pixel Level Decorrelation (PLD) systematics removal method.
Special case of `.RegressionCorrector` where the `.DesignMatrix` is
composed of background-corrected pixel time series.
The design matrix also contains columns representing a spline in time
design to capture the intrinsic, long-term variability of the target.
Pixel Level Decorrelation (PLD) was developed by [1]_ to remove
systematic noise caused by spacecraft jitter for the Spitzer
Space Telescope. It was adapted to K2 data by [2]_ and [3]_
for the EVEREST pipeline [4]_.
For a detailed description and implementation of PLD, please refer to
these references. Lightkurve provides a reference implementation
of PLD that is less sophisticated than EVEREST, but is suitable
for quick-look analyses and detrending experiments.
Our simple implementation of PLD is performed by first calculating the
noise model for each cadence in time. This function goes up to arbitrary
order, and is represented by
.. math::
m_i = \sum_l a_l \frac{f_{il}}{\sum_k f_{ik}} + \sum_l \sum_m b_{lm} \frac{f_{il}f_{im}}{\left( \sum_k f_{ik} \right)^2} + ...
where
- :math:`m_i` is the noise model at time :math:`t_i`
- :math:`f_{il}` is the flux in the :math:`l^\text{th}` pixel at time :math:`t_i`
- :math:`a_l` is the first-order PLD coefficient on the linear term
- :math:`b_{lm}` is the second-order PLD coefficient on the :math:`l^\text{th}`,
:math:`m^\text{th}` pixel pair
We perform Principal Component Analysis (PCA) to reduce the number of
vectors in our final model to limit the set to best capture instrumental
noise. With a PCA-reduced set of vectors, we can construct a design matrix
containing fractional pixel fluxes.
To solve for the PLD model, we need to minimize the difference squared
.. math::
\chi^2 = \sum_i \frac{(y_i - m_i)^2}{\sigma_i^2},
where :math:`y_i` is the observed flux value at time :math:`t_i`, by solving
.. math::
\frac{\partial \chi^2}{\partial a_l} = 0.
The design matrix also contains columns representing a spline in time
design to capture the intrinsic, long-term variability of the target.
Examples
--------
Download the pixel data for GJ 9827 and obtain a PLD-corrected light curve:
>>> import lightkurve as lk
>>> tpf = lk.search_targetpixelfile("GJ9827").download() # doctest: +SKIP
>>> corrector = tpf.to_corrector('pld') # doctest: +SKIP
>>> lc = corrector.correct() # doctest: +SKIP
>>> lc.plot() # doctest: +SKIP
However, the above example will over-fit the small transits!
It is necessary to mask the transits using `corrector.correct(cadence_mask=...)`.
References
----------
.. [1] Deming et al. (2015), ads:2015ApJ...805..132D.
(arXiv:1411.7404)
.. [2] Luger et al. (2016), ads:2016AJ....152..100L
(arXiv:1607.00524)
.. [3] Luger et al. (2018), ads:2018AJ....156...99L
(arXiv:1702.05488)
.. [4] EVEREST pipeline webpage, https://rodluger.github.io/everest
"""
def __init__(self, tpf, aperture_mask=None):
if aperture_mask is None:
aperture_mask = tpf.create_threshold_mask(3)
self.aperture_mask = aperture_mask
lc = tpf.to_lightcurve(aperture_mask=aperture_mask)
# Remove cadences that have NaN flux (cf. #874). We don't simply call
# `lc.remove_nans()` here because we need to mask both lc & tpf.
nan_mask = np.isnan(lc.flux)
lc = lc[~nan_mask]
self.tpf = tpf[~nan_mask]
super().__init__(lc=lc)
def __repr__(self):
return "PLDCorrector (ID: {})".format(self.lc.label)
def create_design_matrix(
self,
pld_order=3,
pca_components=16,
pld_aperture_mask=None,
background_aperture_mask="background",
spline_n_knots=None,
spline_degree=3,
normalize_background_pixels=None,
sparse=False,
):
"""Returns a `.DesignMatrixCollection` containing a `DesignMatrix` object
for the background regressors, the PLD pixel component regressors, and
the spline regressors.
If the parameters `pld_order` and `pca_components` are None, their
value will be assigned based on the mission. K2 and TESS experience
different dominant sources of noise, and require different defaults.
For information about how the defaults were chosen, see Pull Request #746.
Parameters
----------
pld_order : int
The order of Pixel Level De-correlation to be performed. First order
(`n=1`) uses only the pixel fluxes to construct the design matrix.
Higher order populates the design matrix with columns constructed
from the products of pixel fluxes.
pca_components : int or tuple of int
Number of terms added to the design matrix for each order of PLD
pixel fluxes. Increasing this value may provide higher precision
at the expense of slower speed and/or overfitting.
If performing PLD with `pld_order > 1`, `pca_components` can be
a tuple containing the number of terms for each order of PLD.
If a single int is passed, the same number of terms will be used
for each order. If zero is passed, PCA will not be performed.
Defaults to 16 for K2 and 8 for TESS.
pld_aperture_mask : array-like, 'pipeline', 'all', 'threshold', or None
A boolean array describing the aperture such that `True` means
that the pixel will be used when selecting the PLD basis vectors.
If `None` or `all` are passed in, all pixels will be used.
If 'pipeline' is passed, the mask suggested by the official pipeline
will be returned.
If 'threshold' is passed, all pixels brighter than 3-sigma above
the median flux will be used.
spline_n_knots : int
Number of knots in spline.
spline_degree : int
Polynomial degree of spline.
sparse : bool
Whether to create `SparseDesignMatrix`.
Returns
-------
dm : `.DesignMatrixCollection`
`.DesignMatrixCollection` containing pixel, background, and spline
components.
"""
# Validate the inputs
pld_aperture_mask = self.tpf._parse_aperture_mask(pld_aperture_mask)
self.pld_aperture_mask = pld_aperture_mask
background_aperture_mask = self.tpf._parse_aperture_mask(
background_aperture_mask
)
self.background_aperture_mask = background_aperture_mask
if spline_n_knots is None:
# Default to a spline per 50 data points
spline_n_knots = int(len(self.lc) / 50)
if sparse:
DMC = SparseDesignMatrixCollection
spline = create_sparse_spline_matrix
else:
DMC = DesignMatrixCollection
spline = create_spline_matrix
# We set the width of all coefficient priors to 10 times the standard
# deviation to prevent the fit from going crazy.
prior_sigma = np.nanstd(self.lc.flux.value) * 10
# Flux normalize background components for K2 and not for TESS by default
bkg_pixels = self.tpf.flux[:, background_aperture_mask].reshape(
len(self.tpf.flux), -1
)
if normalize_background_pixels:
bkg_flux = np.nansum(self.tpf.flux[:, background_aperture_mask], -1)
bkg_pixels = np.array([r / f for r, f in zip(bkg_pixels, bkg_flux)])
else:
bkg_pixels = bkg_pixels.value
# Remove NaNs
bkg_pixels = np.array([r[np.isfinite(r)] for r in bkg_pixels])
# Create background design matrix
dm_bkg = DesignMatrix(bkg_pixels, name="background")
# Apply PCA
dm_bkg = dm_bkg.pca(pca_components)
# Set prior sigma to 10 * standard deviation
dm_bkg.prior_sigma = np.ones(dm_bkg.shape[1]) * prior_sigma
# Create a design matric containing splines plus a constant
dm_spline = spline(
self.lc.time.value, n_knots=spline_n_knots, degree=spline_degree
).append_constant()
# Set prior sigma to 10 * standard deviation
dm_spline.prior_sigma = np.ones(dm_spline.shape[1]) * prior_sigma
# Create a PLD matrix if there are pixels in the pld_aperture_mask
if np.sum(pld_aperture_mask) != 0:
# Flux normalize the PLD components
pld_pixels = self.tpf.flux[:, pld_aperture_mask].reshape(
len(self.tpf.flux), -1
)
pld_pixels = np.array(
[r / f for r, f in zip(pld_pixels, self.lc.flux.value)]
)
# Remove NaNs
pld_pixels = np.array([r[np.isfinite(r)] for r in pld_pixels])
# Use the DesignMatrix infrastructure to apply PCA to the regressors.
regressors_dm = DesignMatrix(pld_pixels)
if pca_components > 0:
regressors_dm = regressors_dm.pca(pca_components)
regressors_pld = regressors_dm.values
# Create a DesignMatrix for each PLD order
all_pld = []
for order in range(1, pld_order + 1):
reg_n = np.product(list(multichoose(regressors_pld.T, order)), axis=1).T
pld_n = DesignMatrix(
reg_n,
prior_sigma=np.ones(reg_n.shape[1]) * prior_sigma / reg_n.shape[1],
name=f"pld_order_{order}",
)
# Apply PCA.
if pca_components > 0:
pld_n = pld_n.pca(pca_components)
# Calling pca() resets the priors, so we set them again.
pld_n.prior_sigma = (
np.ones(pld_n.shape[1]) * prior_sigma / pca_components
)
all_pld.append(pld_n)
# Create the collection of DesignMatrix objects.
# DesignMatrix 1 contains the PLD pixel series
dm_pixels = DesignMatrixCollection(all_pld).to_designmatrix(
name="pixel_series"
)
with warnings.catch_warnings():
warnings.filterwarnings(
"ignore",
message=".*Not all matrices are `SparseDesignMatrix` objects..*",
)
dm_collection = DMC([dm_pixels, dm_bkg, dm_spline])
else:
with warnings.catch_warnings():
warnings.filterwarnings(
"ignore",
message=".*Not all matrices are `SparseDesignMatrix` objects..*",
)
dm_collection = DMC([dm_bkg, dm_spline])
return dm_collection
@deprecated_renamed_argument(
"n_pca_terms",
"pca_components",
"2.0",
warning_type=LightkurveDeprecationWarning,
)
@deprecated_renamed_argument(
"use_gp", None, "2.0", warning_type=LightkurveDeprecationWarning
)
@deprecated_renamed_argument(
"gp_timescale", None, "2.0", warning_type=LightkurveDeprecationWarning
)
@deprecated_renamed_argument(
"aperture_mask", None, "2.0", warning_type=LightkurveDeprecationWarning
)
def correct(
self,
pld_order=None,
pca_components=None,
pld_aperture_mask=None,
background_aperture_mask="background",
spline_n_knots=None,
spline_degree=5,
normalize_background_pixels=None,
restore_trend=True,
sparse=False,
cadence_mask=None,
sigma=5,
niters=5,
propagate_errors=False,
use_gp=None,
gp_timescale=None,
aperture_mask=None,
):
"""Returns a systematics-corrected light curve.
If the parameters `pld_order` and `pca_components` are None, their
value will be assigned based on the mission. K2 and TESS experience
different dominant sources of noise, and require different defaults.
For information about how the defaults were chosen, see PR #746 at
https://github.com/lightkurve/lightkurve/pull/746#issuecomment-658458270
Parameters
----------
pld_order : int
The order of Pixel Level De-correlation to be performed. First order
(`n=1`) uses only the pixel fluxes to construct the design matrix.
Higher order populates the design matrix with columns constructed
from the products of pixel fluxes. Default 3 for K2 and 1 for TESS.
pca_components : int
Number of terms added to the design matrix for each order of PLD
pixel fluxes. Increasing this value may provide higher precision
at the expense of slower speed and/or overfitting.
pld_aperture_mask : array-like, 'pipeline', 'all', 'threshold', or None
A boolean array describing the aperture such that `True` means
that the pixel will be used when selecting the PLD basis vectors.
If `None` or `all` are passed in, all pixels will be used.
If 'pipeline' is passed, the mask suggested by the official pipeline
will be returned.
If 'threshold' is passed, all pixels brighter than 3-sigma above
the median flux will be used.
spline_n_knots : int
Number of knots in spline.
spline_degree : int
Polynomial degree of spline.
restore_trend : bool
Whether to restore the long term spline trend to the light curve.
sparse : bool
Whether to create `SparseDesignMatrix`.
cadence_mask : np.ndarray of bools (optional)
Mask, where True indicates a cadence that should be used.
sigma : int (default 5)
Standard deviation at which to remove outliers from fitting
niters : int (default 5)
Number of iterations to fit and remove outliers
propagate_errors : bool (default False)
Whether to propagate the uncertainties from the regression. Default is False.
Setting to True will increase run time, but will sample from multivariate normal
distribution of weights.
use_gp, gp_timescale : DEPRECATED
As of Lightkurve v2.0 PLDCorrector uses splines instead of Gaussian Processes.
aperture_mask : DEPRECATED
As of Lightkurve v2.0 the `aperture_mask` parameter needs to be
passed to the class constructor.
Returns
-------
clc : `.LightCurve`
Noise-corrected `.LightCurve`.
"""
self.restore_trend = restore_trend
# Set mission-specific values for pld_order and pca_components
if pld_order is None:
if self.tpf.meta.get("MISSION") == "K2":
pld_order = 3
else:
pld_order = 1
if pca_components is None:
if self.tpf.meta.get("MISSION") == "K2":
pca_components = 16
else:
pca_components = 3
if pld_aperture_mask is None:
if self.tpf.meta.get("MISSION") == "K2":
# K2 noise is dominated by motion
pld_aperture_mask = "threshold"
else:
# TESS noise is dominated by background
pld_aperture_mask = "empty"
if normalize_background_pixels is None:
if self.tpf.meta.get("MISSION") == "K2":
normalize_background_pixels = True
else:
normalize_background_pixels = False
dm = self.create_design_matrix(
pld_aperture_mask=pld_aperture_mask,
background_aperture_mask=background_aperture_mask,
pld_order=pld_order,
pca_components=pca_components,
spline_n_knots=spline_n_knots,
spline_degree=spline_degree,
normalize_background_pixels=normalize_background_pixels,
sparse=sparse,
)
clc = super().correct(
dm,
cadence_mask=cadence_mask,
sigma=sigma,
niters=niters,
propagate_errors=propagate_errors,
)
if restore_trend:
clc += self.diagnostic_lightcurves["spline"] - np.median(
self.diagnostic_lightcurves["spline"].flux
)
return clc
def diagnose(self):
"""Returns diagnostic plots to assess the most recent call to `correct()`.
If `correct()` has not yet been called, a ``ValueError`` will be raised.
Returns
-------
`~matplotlib.axes.Axes`
The matplotlib axes object.
"""
if not getattr(self, "corrected_lc"):
raise ValueError(
"You need to call the `correct()` method "
"before you can call `diagnose()`."
)
names = self.diagnostic_lightcurves.keys()
# Plot the right version of corrected light curve
if self.restore_trend:
clc = (
self.corrected_lc
+ self.diagnostic_lightcurves["spline"]
- np.median(self.diagnostic_lightcurves["spline"].flux)
)
else:
clc = self.corrected_lc
uncorr_cdpp = self.lc.estimate_cdpp()
corr_cdpp = clc.estimate_cdpp()
# Get y-axis limits
ylim = [
min(min(self.lc.flux.value), min(clc.flux.value)),
max(max(self.lc.flux.value), max(clc.flux.value)),
]
# Use lightkurve plotting style
with plt.style.context(MPLSTYLE):
# Plot background model
_, axs = plt.subplots(3, figsize=(10, 9), sharex=True)
ax = axs[0]
self.lc.plot(
ax=ax,
normalize=False,
clip_outliers=True,
label=f"uncorrected ({uncorr_cdpp:.0f})",
)
ax.set_xlabel("")
ax.set_ylim(ylim) # use same ylim for all plots
# Plot pixel and spline components
ax = axs[1]
clc.plot(
ax=ax, normalize=False, alpha=0.4, label=f"corrected ({corr_cdpp:.0f})"
)
for key, color in zip(names, ["dodgerblue", "r", "C1"]):
if key in ["background", "spline", "pixel_series"]:
tmplc = (
self.diagnostic_lightcurves[key]
- np.median(self.diagnostic_lightcurves[key].flux)
+ np.median(self.lc.flux)
)
tmplc.plot(ax=ax, c=color)
ax.set_xlabel("")
ax.set_ylim(ylim)
# Plot final corrected light curve with outliers marked
ax = axs[2]
self.lc.plot(
ax=ax,
normalize=False,
alpha=0.2,
label=f"uncorrected ({uncorr_cdpp:.0f})",
)
clc[self.outlier_mask].scatter(
normalize=False, c="r", marker="x", s=10, label="outlier_mask", ax=ax
)
clc[~self.cadence_mask].scatter(
normalize=False,
c="dodgerblue",
marker="x",
s=10,
label="~cadence_mask",
ax=ax,
)
clc.plot(
normalize=False, ax=ax, c="k", label=f"corrected ({corr_cdpp:.0f})"
)
ax.set_ylim(ylim)
return axs
def diagnose_masks(self):
"""Show different aperture masks used by PLD in the most recent call to
`correct()`. If `correct()` has not yet been called, a ``ValueError``
will be raised.
Returns
-------
`~matplotlib.axes.Axes`
The matplotlib axes object.
"""
if not hasattr(self, "corrected_lc"):
raise ValueError(
"You need to call the `correct()` method "
"before you can call `diagnose()`."
)
# Use lightkurve plotting style
with plt.style.context(MPLSTYLE):
_, axs = plt.subplots(1, 3, figsize=(12, 4), sharey=True)
# Show light curve aperture mask
ax = axs[0]
self.tpf.plot(
ax=ax,
show_colorbar=False,
aperture_mask=self.aperture_mask,
title="aperture_mask",
)
# Show PLD pixel mask
ax = axs[1]
self.tpf.plot(
ax=ax,
show_colorbar=False,
aperture_mask=self.pld_aperture_mask,
title="pld_aperture_mask",
)
ax = axs[2]
self.tpf.plot(
ax=ax,
show_colorbar=False,
aperture_mask=self.background_aperture_mask,
title="background_aperture_mask",
)
return axs
# `TessPLDCorrector` was briefly introduced in Lightkurve v1.9
# but was removed in v2.0 in favor of a single generic `PLDCorrector`.
@deprecated(
"2.0", alternative="PLDCorrector", warning_type=LightkurveDeprecationWarning
)
class TessPLDCorrector(PLDCorrector):
pass