Comparing Lomb-Scargle Implementations

We want to make sure here that the astroML versions, which are based on optimized algorithms, match the results of the multiband_LS versions, which are based on less-optimal algorithms.

(Note that this is also checked in the multiband_LS unit tests, but it's nice to have visual confirmation here)



In [1]:

    
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt

import seaborn
seaborn.set()



In [2]:

    
# import path is up one directory
import os, sys; sys.path.append(os.path.abspath('..'))

from multiband_LS import LombScargle, LombScargleAstroML
from astroML.time_series import multiterm_periodogram



In [3]:

    
from astroML.datasets import fetch_LINEAR_sample
data = fetch_LINEAR_sample()









    



[10003298 10004892 10013411 ...,  9984569  9987252   999528]



In [4]:

    
lcid = 10022663
t, y, dy = data[lcid].T



In [5]:

    
plt.errorbar(t, y, dy, fmt='o');

Standard Lomb-Scargle



In [6]:

    
period = np.linspace(0.61, 0.62, 5000)

P1 = LombScargleAstroML(fit_offset=False).fit(t, y, dy).periodogram(period)
P2 = LombScargle(fit_offset=False).fit(t, y, dy).periodogram(period)
print("methods match:", np.allclose(P1, P2))

plt.plot(period, P1, label='astroML')
plt.plot(period, P2, label='multiterm')
plt.legend();









    



methods match: True

Timings



In [7]:

    
%timeit LombScargleAstroML(fit_offset=False).fit(t, y, dy).periodogram(period)
%timeit LombScargle(fit_offset=False).fit(t, y, dy).periodogram(period)









    



10 loops, best of 3: 129 ms per loop
1 loops, best of 3: 428 ms per loop

Generalized Lomb-Scargle



In [8]:

    
P1 = LombScargleAstroML(fit_offset=True).fit(t, y, dy).periodogram(period)
P2 = LombScargle(fit_offset=True).fit(t, y, dy).periodogram(period)
print("methods match:", np.allclose(P1, P2))

plt.plot(period, P1, label='astroML')
plt.plot(period, P2, label='multiterm')
plt.legend();









    



methods match: True

Timings



In [9]:

    
%timeit LombScargleAstroML(fit_offset=True).fit(t, y, dy).periodogram(period)
%timeit LombScargle(fit_offset=True).fit(t, y, dy).periodogram(period)









    



10 loops, best of 3: 139 ms per loop
1 loops, best of 3: 526 ms per loop

Multiterm Fits



In [10]:

    
P1 = multiterm_periodogram(t, y, dy, 2 * np.pi / period, n_terms=3)
P2 = LombScargle(fit_offset=True, Nterms=3).fit(t, y, dy).periodogram(period)
print("methods match:", np.allclose(P1, P2))

plt.plot(period, P1, label='astroML')
plt.plot(period, P2, label='multiterm')
plt.legend();









    



methods match: True

Timings



In [11]:

    
%timeit multiterm_periodogram(t, y, dy, 2 * np.pi / period, n_terms=3)
%timeit LombScargle(fit_offset=True, Nterms=3).fit(t, y, dy).periodogram(period)









    



1 loops, best of 3: 988 ms per loop
1 loops, best of 3: 898 ms per loop

Poorly-behaved data

Here we'll check out the difference between standard and generalized Lomb-Scargle...



In [12]:

    
np.random.seed(0)

lcid = 10022663
t, y, dy = data[lcid].T

# Find best period
P_best = LombScargleAstroML().fit(t, y, dy).best_period

# Fit a sinusoid
y = 5 * np.mean(dy) * np.sin(2 * np.pi * t / P_best) + np.random.normal(0, dy)
ymean = y.mean()









    



Finding optimal frequency:
 - Using omega_step = 0.00064
 - Computing periods at 40956 steps from 0.20 to 1.20
Zooming-in on 5 candidate peaks:
 - Computing periods at 1010 steps



In [13]:

    
plt.errorbar(t % P_best, y, dy, fmt='o');



In [14]:

    
# Screw-up lomb-scargle by using only the top half of the data

mask = (y > ymean)
t, y, dy = t[mask], y[mask], dy[mask]
plt.errorbar(t % P_best, y, dy, fmt='o');

Standard



In [15]:

    
period = np.linspace(0.61, 0.62, 5000)

P1 = LombScargleAstroML(fit_offset=False).fit(t, y, dy).periodogram(period)
P2 = LombScargle(fit_offset=False).fit(t, y, dy).periodogram(period)
print("methods match:", np.allclose(P1, P2))

print(period[np.argmax(P1)])

plt.plot(period, P1, label='astroML')
plt.plot(period, P2, label='multiterm')
plt.legend();









    



methods match: True
0.611622324465

Generalized



In [16]:

    
P1 = LombScargleAstroML(fit_offset=True).fit(t, y, dy).periodogram(period)
P2 = LombScargle(fit_offset=True).fit(t, y, dy).periodogram(period)
print("methods match:", np.allclose(P1, P2))

plt.plot(period, P1, label='astroML')
plt.plot(period, P2, label='multiterm')
plt.legend();









    



methods match: True