From 7e142cffefedb610cef091edf182e54475039696 Mon Sep 17 00:00:00 2001 From: Abhimat Gautam Date: Sun, 2 Jun 2019 19:59:27 -0700 Subject: [PATCH 1/6] Implemented multiband support for trended Lomb-Scargle --- gatspy/periodic/__init__.py | 4 +- .../trended_lomb_scargle_multiband.py | 106 ++++++++++++++++++ 2 files changed, 109 insertions(+), 1 deletion(-) create mode 100644 gatspy/periodic/trended_lomb_scargle_multiband.py diff --git a/gatspy/periodic/__init__.py b/gatspy/periodic/__init__.py index 7d20402..9cb28cd 100644 --- a/gatspy/periodic/__init__.py +++ b/gatspy/periodic/__init__.py @@ -6,7 +6,8 @@ __all__ = ['LombScargle', 'LombScargleFast', 'LombScargleAstroML', 'LombScargleMultiband', 'LombScargleMultibandFast', - 'TrendedLombScargle', 'SuperSmoother', 'SuperSmootherMultiband', + 'TrendedLombScargle', 'TrendedLombScargleMultiband', + 'SuperSmoother', 'SuperSmootherMultiband', 'RRLyraeTemplateModeler', 'RRLyraeTemplateModelerMultiband', 'NaiveMultiband'] @@ -14,6 +15,7 @@ from .lomb_scargle_fast import * from .lomb_scargle_multiband import * from .trended_lomb_scargle import * +from .trended_lomb_scargle_multiband import * from .supersmoother import * from .template_modeler import * from .naive_multiband import * diff --git a/gatspy/periodic/trended_lomb_scargle_multiband.py b/gatspy/periodic/trended_lomb_scargle_multiband.py new file mode 100644 index 0000000..fdc4dbc --- /dev/null +++ b/gatspy/periodic/trended_lomb_scargle_multiband.py @@ -0,0 +1,106 @@ +from __future__ import division, print_function, absolute_import + +__all__ = ['TrendedLombScargleMultiband'] + +import warnings + +import numpy as np + +from .modeler import PeriodicModeler +from .lomb_scargle import LombScargle +from .lomb_scargle_multiband import LombScargleMultiband + + +class TrendedLombScargleMultiband(LombScargleMultiband): + """Trended Lomb-Scargle Periodogram Implementation + + This is a generalized periodogram implementation using the matrix formalism + outlined in VanderPlas & Ivezic 2015. It fits both a floating mean + and a trend parameter (as opposed to the `LombScargle` class, which + fits only the mean). + + Parameters + ---------- + optimizer : PeriodicOptimizer instance + Optimizer to use to find the best period. If not specified, the + LinearScanOptimizer will be used. + center_data : boolean (default = True) + If True, then compute the weighted mean of the input data and subtract + before fitting the model. + fit_offset : boolean (default = True) + If True, then fit a floating-mean sinusoid model. + Nterms : int (default = 1) + Number of Fourier frequencies to fit in the model + regularization : float, vector or None (default = None) + If specified, then add this regularization penalty to the + least squares fit. + regularize_by_trace : boolean (default = True) + If True, multiply regularization by the trace of the matrix + fit_period : bool (optional) + If True, then fit for the best period when fit() method is called. + optimizer_kwds : dict (optional) + Dictionary of keyword arguments for constructing the optimizer. For + example, silence optimizer output with `optimizer_kwds={"quiet": True}`. + + Examples + -------- + >>> rng = np.random.RandomState(0) + >>> t = 100 * rng.rand(100) + >>> dy = 0.1 + >>> omega = 10 + >>> slope = 2. + >>> y = np.sin(omega * t) + slope * t + dy * rng.randn(100) + >>> ls = TrendedLombScargle().fit(t, y, dy) + >>> ls.optimizer.period_range = (0.2, 1.2) + >>> ls.best_period + Finding optimal frequency: + - Estimated peak width = 0.0639 + - Using 5 steps per peak; omega_step = 0.0128 + - User-specified period range: 0.2 to 1.2 + - Computing periods at 2051 steps + Zooming-in on 5 candidate peaks: + - Computing periods at 1000 steps + 0.62827068275990694 + >>> ls.predict([0, 0.5]) + array([-0.01144474, 0.07567192]) + + See Also + -------- + LombScargle + LombScargleAstroML + LombScargleMultiband + LombScargleMultibandFast + """ + + def _construct_X(self, omega, weighted=True, **kwargs): + t = kwargs.get('t', self.t) + dy = kwargs.get('dy', self.dy) + filts = kwargs.get('filts', self.filts) + + # X is a huge-ass matrix that has lots of zeros depending on + # which filters are present... + # + # huge-ass, quantitatively speaking, is of shape + # [len(t), (1 + 1 + 2 * Nterms_base + nfilts * (1 + 1 + 2 * Nterms_band))] + + # TODO: the building of the matrix could be more efficient + cols = [np.ones(len(t))] + cols.append(t) # coefficients for trend parameter + cols = sum(([np.sin((i + 1) * omega * t), + np.cos((i + 1) * omega * t)] + for i in range(self.Nterms_base)), cols) + + for filt in self.unique_filts_: + cols.append(np.ones(len(t))) + cols.append(t) # coefficients for trend parameter (in filt) + cols = sum(([np.sin((i + 1) * omega * t), + np.cos((i + 1) * omega * t)] + for i in range(self.Nterms_band)), cols) + mask = (filts == filt) + for i in range(-1 - 2 * self.Nterms_band, 0): + cols[i][~mask] = 0 + + if weighted: + return np.transpose(np.vstack(cols) / dy) + else: + return np.transpose(np.vstack(cols)) From 3a3e1769b79d09c437c74476b75f6672160f628c Mon Sep 17 00:00:00 2001 From: Abhimat Gautam Date: Mon, 3 Jun 2019 11:12:24 -0700 Subject: [PATCH 2/6] Fixed multiband and design matrix for trended multiband --- .../trended_lomb_scargle_multiband.py | 23 +++++++++++++++---- 1 file changed, 18 insertions(+), 5 deletions(-) diff --git a/gatspy/periodic/trended_lomb_scargle_multiband.py b/gatspy/periodic/trended_lomb_scargle_multiband.py index fdc4dbc..657c3ec 100644 --- a/gatspy/periodic/trended_lomb_scargle_multiband.py +++ b/gatspy/periodic/trended_lomb_scargle_multiband.py @@ -72,6 +72,19 @@ class TrendedLombScargleMultiband(LombScargleMultiband): LombScargleMultibandFast """ + def _construct_regularization(self): + if self.reg_base is None and self.reg_band is None: + reg = 0 + else: + Nbase = 2 + 2 * self.Nterms_base + Nband = 2 + 2 * self.Nterms_band + reg = np.zeros(Nbase + len(self.unique_filts_) * Nband) + if self.reg_base is not None: + reg[:Nbase] = self.reg_base + if self.reg_band is not None: + reg[Nbase:] = self.reg_band + return reg + def _construct_X(self, omega, weighted=True, **kwargs): t = kwargs.get('t', self.t) dy = kwargs.get('dy', self.dy) @@ -85,21 +98,21 @@ def _construct_X(self, omega, weighted=True, **kwargs): # TODO: the building of the matrix could be more efficient cols = [np.ones(len(t))] - cols.append(t) # coefficients for trend parameter + cols.append(np.copy(t)) # coefficients for trend parameter cols = sum(([np.sin((i + 1) * omega * t), np.cos((i + 1) * omega * t)] for i in range(self.Nterms_base)), cols) - + for filt in self.unique_filts_: cols.append(np.ones(len(t))) - cols.append(t) # coefficients for trend parameter (in filt) + cols.append(np.copy(t)) # coefficients for trend parameter (in filt) cols = sum(([np.sin((i + 1) * omega * t), np.cos((i + 1) * omega * t)] for i in range(self.Nterms_band)), cols) mask = (filts == filt) - for i in range(-1 - 2 * self.Nterms_band, 0): + for i in range(-2 - (2 * self.Nterms_band), 0): cols[i][~mask] = 0 - + if weighted: return np.transpose(np.vstack(cols) / dy) else: From d8ac05bd4780d209f89a01243f20a6b7995130cb Mon Sep 17 00:00:00 2001 From: Abhimat Gautam Date: Wed, 5 Jun 2019 14:26:12 -0700 Subject: [PATCH 3/6] Cleaned up documentation for trended_lomb_scargle_multiband --- .../trended_lomb_scargle_multiband.py | 62 ++++++------------- 1 file changed, 19 insertions(+), 43 deletions(-) diff --git a/gatspy/periodic/trended_lomb_scargle_multiband.py b/gatspy/periodic/trended_lomb_scargle_multiband.py index 657c3ec..fa7ec1e 100644 --- a/gatspy/periodic/trended_lomb_scargle_multiband.py +++ b/gatspy/periodic/trended_lomb_scargle_multiband.py @@ -12,64 +12,40 @@ class TrendedLombScargleMultiband(LombScargleMultiband): - """Trended Lomb-Scargle Periodogram Implementation + """Trended Multiband Lomb-Scargle Periodogram Implementation - This is a generalized periodogram implementation using the matrix formalism - outlined in VanderPlas & Ivezic 2015. It fits both a floating mean - and a trend parameter (as opposed to the `LombScargle` class, which - fits only the mean). + This is a generalized multiband periodogram implementation using the matrix + formalism outlined in VanderPlas & Ivezic 2015. It fits both a floating mean + and a trend parameter in each band. Parameters ---------- optimizer : PeriodicOptimizer instance Optimizer to use to find the best period. If not specified, the LinearScanOptimizer will be used. + Nterms_base : integer (default = 1) + number of frequency terms to use for the base model common to all bands + Nterms_band : integer (default = 1) + number of frequency terms to use for the residuals between the base + model and each individual band + reg_base : float or None (default = None) + amount of regularization to use on the base model parameters + reg_band : float or None (default = 1E-6) + amount of regularization to use on the band model parameters + regularize_by_trace : bool (default = True) + if True, then regularization is expressed in units of the trace of + the normal matrix center_data : boolean (default = True) - If True, then compute the weighted mean of the input data and subtract - before fitting the model. - fit_offset : boolean (default = True) - If True, then fit a floating-mean sinusoid model. - Nterms : int (default = 1) - Number of Fourier frequencies to fit in the model - regularization : float, vector or None (default = None) - If specified, then add this regularization penalty to the - least squares fit. - regularize_by_trace : boolean (default = True) - If True, multiply regularization by the trace of the matrix - fit_period : bool (optional) - If True, then fit for the best period when fit() method is called. + if True, then center the y data prior to the fit optimizer_kwds : dict (optional) Dictionary of keyword arguments for constructing the optimizer. For example, silence optimizer output with `optimizer_kwds={"quiet": True}`. - Examples - -------- - >>> rng = np.random.RandomState(0) - >>> t = 100 * rng.rand(100) - >>> dy = 0.1 - >>> omega = 10 - >>> slope = 2. - >>> y = np.sin(omega * t) + slope * t + dy * rng.randn(100) - >>> ls = TrendedLombScargle().fit(t, y, dy) - >>> ls.optimizer.period_range = (0.2, 1.2) - >>> ls.best_period - Finding optimal frequency: - - Estimated peak width = 0.0639 - - Using 5 steps per peak; omega_step = 0.0128 - - User-specified period range: 0.2 to 1.2 - - Computing periods at 2051 steps - Zooming-in on 5 candidate peaks: - - Computing periods at 1000 steps - 0.62827068275990694 - >>> ls.predict([0, 0.5]) - array([-0.01144474, 0.07567192]) - See Also -------- LombScargle - LombScargleAstroML LombScargleMultiband - LombScargleMultibandFast + TrendedLombScargle """ def _construct_regularization(self): @@ -94,7 +70,7 @@ def _construct_X(self, omega, weighted=True, **kwargs): # which filters are present... # # huge-ass, quantitatively speaking, is of shape - # [len(t), (1 + 1 + 2 * Nterms_base + nfilts * (1 + 1 + 2 * Nterms_band))] + # [len(t), (2 + 2 * Nterms_base + nfilts * (2 + 2 * Nterms_band))] # TODO: the building of the matrix could be more efficient cols = [np.ones(len(t))] From 5aba05a839347eef1552cd108b8d3301d3ce63e0 Mon Sep 17 00:00:00 2001 From: Abhimat Gautam Date: Thu, 6 Jun 2019 22:53:40 -0700 Subject: [PATCH 4/6] Basic test for multiband trended Lomb-Scargle --- .../test_trended_lomb_scargle_multiband.py | 41 +++++++++++++++++++ 1 file changed, 41 insertions(+) create mode 100755 gatspy/periodic/tests/test_trended_lomb_scargle_multiband.py diff --git a/gatspy/periodic/tests/test_trended_lomb_scargle_multiband.py b/gatspy/periodic/tests/test_trended_lomb_scargle_multiband.py new file mode 100755 index 0000000..5a7242b --- /dev/null +++ b/gatspy/periodic/tests/test_trended_lomb_scargle_multiband.py @@ -0,0 +1,41 @@ +from __future__ import division + +import numpy as np +from numpy.testing import assert_allclose, assert_, assert_raises +from nose import SkipTest + +from .. import (LombScargle, TrendedLombScargleMultiband) + + +def _generate_data(N=100, period=1, theta=[10, 2, 3], dy=1, rseed=0): + """Generate some data for testing""" + rng = np.random.RandomState(rseed) + t = 20 * period * rng.rand(N) + omega = 2 * np.pi / period + y = theta[0] + theta[1] * np.sin(omega * t) + theta[2] * np.cos(omega * t) + dy = dy * (0.5 + rng.rand(N)) + y += dy * rng.randn(N) + return t, y, dy + + +def test_multiband_lomb_scargle_linear_trend(N=100, period=1, slope=5): + """Test whether the generalized lomb-scargle properly fits data with a + linear trend component (single filter) + """ + t, y, dy = _generate_data(N, period) + omega = 2 * np.pi / period + + # Model and optimizer for standard Lomb-Scargle + model = LombScargle(center_data=False, fit_offset=True) + model.fit(t, y, dy) + model.optimizer.period_range = (period - 0.5, period + 0.5) + y_hat = model.predict(t) + + # Model and optimizer for trended multiband Lomb-Scargle + y_mb_t = y + slope * t + model_mb_t = TrendedLombScargleMultiband(center_data=False) + model_mb_t.fit(t, y, dy, filts=1) + model_mb_t.optimizer.period_range = (period - 0.5, period + 0.5) + y_hat_mb_t = model_mb_t.predict(t) + + assert_allclose(y_hat, y_hat_mb_t - slope * t, rtol=5E-2) From bca9c65fb105a0f9b038766f10f2836983f40336 Mon Sep 17 00:00:00 2001 From: Abhimat Gautam Date: Thu, 12 Jan 2023 13:46:27 -0800 Subject: [PATCH 5/6] New class for polynomial trends --- gatspy/periodic/__init__.py | 1 + .../trended_lomb_scargle_multiband.py | 112 +++++++++++++++++- 2 files changed, 112 insertions(+), 1 deletion(-) diff --git a/gatspy/periodic/__init__.py b/gatspy/periodic/__init__.py index 9cb28cd..454c612 100644 --- a/gatspy/periodic/__init__.py +++ b/gatspy/periodic/__init__.py @@ -7,6 +7,7 @@ __all__ = ['LombScargle', 'LombScargleFast', 'LombScargleAstroML', 'LombScargleMultiband', 'LombScargleMultibandFast', 'TrendedLombScargle', 'TrendedLombScargleMultiband', + 'PolyTrend_LombScargleMultiband', 'SuperSmoother', 'SuperSmootherMultiband', 'RRLyraeTemplateModeler', 'RRLyraeTemplateModelerMultiband', 'NaiveMultiband'] diff --git a/gatspy/periodic/trended_lomb_scargle_multiband.py b/gatspy/periodic/trended_lomb_scargle_multiband.py index fa7ec1e..dd35b1e 100644 --- a/gatspy/periodic/trended_lomb_scargle_multiband.py +++ b/gatspy/periodic/trended_lomb_scargle_multiband.py @@ -1,6 +1,6 @@ from __future__ import division, print_function, absolute_import -__all__ = ['TrendedLombScargleMultiband'] +__all__ = ['TrendedLombScargleMultiband', 'PolyTrend_LombScargleMultiband'] import warnings @@ -93,3 +93,113 @@ def _construct_X(self, omega, weighted=True, **kwargs): return np.transpose(np.vstack(cols) / dy) else: return np.transpose(np.vstack(cols)) + + +class PolyTrend_LombScargleMultiband(LombScargleMultiband): + """Trended Multiband Lomb-Scargle Periodogram Implementation + + This is a generalized multiband periodogram implementation using the matrix + formalism outlined in VanderPlas & Ivezic 2015. It fits both a floating mean + and a trend parameter in each band. + + Parameters + ---------- + optimizer : PeriodicOptimizer instance + Optimizer to use to find the best period. If not specified, the + LinearScanOptimizer will be used. + Nterms_base : integer (default = 1) + number of frequency terms to use for the base model common to all bands + Nterms_band : integer (default = 1) + number of frequency terms to use for the residuals between the base + model and each individual band + reg_base : float or None (default = None) + amount of regularization to use on the base model parameters + reg_band : float or None (default = 1E-6) + amount of regularization to use on the band model parameters + regularize_by_trace : bool (default = True) + if True, then regularization is expressed in units of the trace of + the normal matrix + center_data : boolean (default = True) + if True, then center the y data prior to the fit + optimizer_kwds : dict (optional) + Dictionary of keyword arguments for constructing the optimizer. For + example, silence optimizer output with `optimizer_kwds={"quiet": True}`. + + See Also + -------- + LombScargle + LombScargleMultiband + TrendedLombScargle + """ + + def __init__(self, optimizer=None, + poly_trend_order=2, + Nterms_base=1, Nterms_band=1, + reg_base=None, reg_band=1E-6, regularize_by_trace=True, + center_data=True, fit_period=False, optimizer_kwds=None): + self.poly_trend_order = poly_trend_order + + LombScargleMultiband.__init__( + self, optimizer=optimizer, + Nterms_base=Nterms_base, Nterms_band=Nterms_band, + reg_base=reg_base, reg_band=reg_band, + regularize_by_trace=regularize_by_trace, + center_data=center_data, fit_period=fit_period, + optimizer_kwds=optimizer_kwds, + ) + + def _construct_regularization(self): + if self.reg_base is None and self.reg_band is None: + reg = 0 + else: + Nbase = (self.poly_trend_order + 1) + 2 * self.Nterms_base + Nband = (self.poly_trend_order + 1) + 2 * self.Nterms_band + reg = np.zeros(Nbase + len(self.unique_filts_) * Nband) + if self.reg_base is not None: + reg[:Nbase] = self.reg_base + if self.reg_band is not None: + reg[Nbase:] = self.reg_band + return reg + + def _construct_X(self, omega, weighted=True, **kwargs): + t = kwargs.get('t', self.t) + dy = kwargs.get('dy', self.dy) + filts = kwargs.get('filts', self.filts) + + # X is a huge-ass matrix that has lots of zeros depending on + # which filters are present... + # + # huge-ass, quantitatively speaking, is of shape + # [len(t), ((poly_trend_order + 1) + 2 * Nterms_base + + # nfilts * ((poly_trend_order + 1) + 2 * Nterms_band))] + + # TODO: the building of the matrix could be more efficient + cols = [np.ones(len(t))] # Coefficients for constant param + + for cur_poly_trend_order in range(1, self.poly_trend_order + 1): + # Coefficients for current polynomial trend order + cols.append(np.copy(t**cur_poly_trend_order)) + + cols = sum(([np.sin((i + 1) * omega * t), + np.cos((i + 1) * omega * t)] + for i in range(self.Nterms_base)), cols) + + for filt in self.unique_filts_: + cols.append(np.ones(len(t))) + + for cur_poly_trend_order in range(1, self.poly_trend_order + 1): + # Coefficients for current polynomial trend order + cols.append(np.copy(t**cur_poly_trend_order)) + + cols = sum(([np.sin((i + 1) * omega * t), + np.cos((i + 1) * omega * t)] + for i in range(self.Nterms_band)), cols) + mask = (filts == filt) + for i in range((-1 * (self.poly_trend_order + 1)) - + (2 * self.Nterms_band), 0): + cols[i][~mask] = 0 + + if weighted: + return np.transpose(np.vstack(cols) / dy) + else: + return np.transpose(np.vstack(cols)) From 51613eff9f3d5bc58dd73aabf5882e4aafd15e65 Mon Sep 17 00:00:00 2001 From: Abhimat Gautam Date: Sun, 15 Jan 2023 15:53:57 -0800 Subject: [PATCH 6/6] Added support for separate orders for base and band poly trends --- .../trended_lomb_scargle_multiband.py | 19 ++++++++++--------- 1 file changed, 10 insertions(+), 9 deletions(-) diff --git a/gatspy/periodic/trended_lomb_scargle_multiband.py b/gatspy/periodic/trended_lomb_scargle_multiband.py index dd35b1e..f732c8d 100644 --- a/gatspy/periodic/trended_lomb_scargle_multiband.py +++ b/gatspy/periodic/trended_lomb_scargle_multiband.py @@ -133,11 +133,12 @@ class PolyTrend_LombScargleMultiband(LombScargleMultiband): """ def __init__(self, optimizer=None, - poly_trend_order=2, + poly_trend_order_base=2, poly_trend_order_band=1, Nterms_base=1, Nterms_band=1, reg_base=None, reg_band=1E-6, regularize_by_trace=True, center_data=True, fit_period=False, optimizer_kwds=None): - self.poly_trend_order = poly_trend_order + self.poly_trend_order_base = poly_trend_order_base + self.poly_trend_order_band = poly_trend_order_band LombScargleMultiband.__init__( self, optimizer=optimizer, @@ -152,8 +153,8 @@ def _construct_regularization(self): if self.reg_base is None and self.reg_band is None: reg = 0 else: - Nbase = (self.poly_trend_order + 1) + 2 * self.Nterms_base - Nband = (self.poly_trend_order + 1) + 2 * self.Nterms_band + Nbase = (self.poly_trend_order_base + 1) + 2 * self.Nterms_base + Nband = (self.poly_trend_order_band + 1) + 2 * self.Nterms_band reg = np.zeros(Nbase + len(self.unique_filts_) * Nband) if self.reg_base is not None: reg[:Nbase] = self.reg_base @@ -170,13 +171,13 @@ def _construct_X(self, omega, weighted=True, **kwargs): # which filters are present... # # huge-ass, quantitatively speaking, is of shape - # [len(t), ((poly_trend_order + 1) + 2 * Nterms_base + - # nfilts * ((poly_trend_order + 1) + 2 * Nterms_band))] + # [len(t), ((poly_trend_order_base + 1) + 2 * Nterms_base + + # nfilts * ((poly_trend_order_band + 1) + 2 * Nterms_band))] # TODO: the building of the matrix could be more efficient cols = [np.ones(len(t))] # Coefficients for constant param - for cur_poly_trend_order in range(1, self.poly_trend_order + 1): + for cur_poly_trend_order in range(1, self.poly_trend_order_base + 1): # Coefficients for current polynomial trend order cols.append(np.copy(t**cur_poly_trend_order)) @@ -187,7 +188,7 @@ def _construct_X(self, omega, weighted=True, **kwargs): for filt in self.unique_filts_: cols.append(np.ones(len(t))) - for cur_poly_trend_order in range(1, self.poly_trend_order + 1): + for cur_poly_trend_order in range(1, self.poly_trend_order_band + 1): # Coefficients for current polynomial trend order cols.append(np.copy(t**cur_poly_trend_order)) @@ -195,7 +196,7 @@ def _construct_X(self, omega, weighted=True, **kwargs): np.cos((i + 1) * omega * t)] for i in range(self.Nterms_band)), cols) mask = (filts == filt) - for i in range((-1 * (self.poly_trend_order + 1)) - + for i in range((-1 * (self.poly_trend_order_band + 1)) - (2 * self.Nterms_band), 0): cols[i][~mask] = 0