Skip to content

Commit

Permalink
Separate method by Press & Rybicki from the generalised algorithm
Browse files Browse the repository at this point in the history
  • Loading branch information
@giordano
giordano committed Aug 17, 2016
1 parent 608fb3f commit 0d178d1
Show file tree
Hide file tree
Showing 2 changed files with 144 additions and 96 deletions.
223 changes: 127 additions & 96 deletions src/LombScargle.jl
View file
Expand Up @@ -247,113 +247,143 @@ function _generalised_lombscargle{R1<:Real,R2<:Real,R3<:Real,R4<:Real}(times::Ab
w::AbstractVector{R3},
freqs::AbstractVector{R4},
center_data::Bool,
fit_mean::Bool,
use_fft::Bool,
oversampling::Int,
Mfft::Int)
P_type = promote_type(float(R1), float(R2))
fit_mean::Bool)
P_type = promote_type(float(R1), float(R2), float(R3), float(R4))
P = Vector{P_type}(freqs)
nil = zero(P_type)
# If "center_data" keyword is true, subtract the mean from each point.
if center_data || fit_mean
if use_fft
y = signal - w⋅signal
else
y = signal - mean(signal)
end
if center_data
y = signal - mean(signal)
else
y = signal
end
YY = w⋅(y.^2)
if !use_fft && fit_mean
if fit_mean
Y = w⋅y
end

if use_fft
df = step(freqs)
N = length(freqs)
f0 = minimum(freqs)
Ch, Sh = trig_sum(times, w .* y, df, N, f0, 1, oversampling, Mfft)
C2, S2 = trig_sum(times, w, df, N, f0, 2, oversampling, Mfft)
@inbounds for n in eachindex(freqs)
ω = 2pi*freqs[n]

# Find τ for current angular frequency
C = S = CS = CC = nil
for i in eachindex(times)
ωt = ω*times[i]
W = w[i]
c = cos(ωt)
s = sin(ωt)
CS += W*c*s
CC += W*c*c
if fit_mean
C += W*c
S += W*s
end
end
if fit_mean
C, S = trig_sum(times, w, df, N, f0, 1, oversampling, Mfft)
tan_2ωτ = (S2 .- 2.0 * S .* C) ./ (C2 .- (C .* C .- S .* S))
CS -= C*S
SS = 1.0 - CC - S*S
CC -= C*C
else
tan_2ωτ = S2 ./ C2
SS = 1.0 - CC
end
τ = 0.5*atan2(2CS, CC - SS)/ω
# These quantities will be used below. See Press & Rybicki paper.
# NOTE: They can be computed in a faster way, for the time being we keep
# this simple and straightforward implementation.
cos_ωτ = cos(ω*τ)
sin_ωτ = sin(ω*τ)

# Now we can compute the power
YC_τ = YS_τ = CC_τ = nil
for i in eachindex(times)
ωt = ω*(times[i] - τ)
W = w[i]
c = cos(ωt)
s = sin(ωt)
YC_τ += W*y[i]*c
YS_τ += W*y[i]*s
CC_τ += W*c*c
end
C2w = 1./(sqrt(1.0 + tan_2ωτ .* tan_2ωτ))
S2w = tan_2ωτ .* C2w
Cw = sqrt(0.5 * (1.0 + C2w))
Sw = sqrt(0.5) .* sign(S2w) .* sqrt(1.0 - C2w)
YC = Ch .* Cw .+ Sh .* Sw
YS = Sh .* Cw .- Ch .* Sw
CC = 0.5 * (1.0 + C2 .* C2w .+ S2 .* S2w)
SS = 0.5 * (1.0 - C2 .* C2w .- S2 .* S2w)
if fit_mean
CC -= (C .* Cw .+ S .* Sw) .^ 2
SS -= (S .* Cw .- C .* Sw) .^ 2
# "C_τ" and "S_τ" are computed following equation (7) of Press &
# Rybicki, this formula simply comes from angle difference
# trigonometric identities.
C_τ = C*cos_ωτ + S*sin_ωτ
S_τ = S*cos_ωτ - C*sin_ωτ
YC_τ -= Y*C_τ
YS_τ -= Y*S_τ
SS_τ = 1.0 - CC_τ - S_τ*S_τ
CC_τ -= C_τ*C_τ
YY -= Y*Y
else
SS_τ = 1.0 - CC_τ
end
P = (YC .* YC ./ CC .+ YS .* YS ./ SS)/YY
P[n] = (abs2(YC_τ)/CC_τ + abs2(YS_τ)/SS_τ)/YY
end
return P, freqs
end

# Fast, but approximate, method to compute the Lomb-Scargle periodogram for
# evenly spaced data. See
# * Press, W. H., Rybicki, G. B. 1989, ApJ, 338, 277 (URL:
# http://dx.doi.org/10.1086/167197,
# Bibcode: http://adsabs.harvard.edu/abs/1989ApJ...338..277P)
# This is adapted from Astropy implementation of the method. See
# `lombscargle_fast' function.
function _press_rybicki{R1<:Real,R2<:Real,R3<:Real,R4<:Real}(times::AbstractVector{R1},
signal::AbstractVector{R2},
w::AbstractVector{R3},
freqs::AbstractVector{R4},
center_data::Bool,
fit_mean::Bool,
oversampling::Int,
Mfft::Int)
P_type = promote_type(float(R1), float(R2))
P = Vector{P_type}(freqs)
nil = zero(P_type)
# If "center_data" keyword is true, subtract the mean from each point.
if center_data || fit_mean
y = signal - w⋅signal
else
@inbounds for n in eachindex(freqs)
ω = 2pi*freqs[n]

# Find τ for current angular frequency
C = S = CS = CC = nil
for i in eachindex(times)
ωt = ω*times[i]
W = w[i]
c = cos(ωt)
s = sin(ωt)
CS += W*c*s
CC += W*c*c
if fit_mean
C += W*c
S += W*s
end
end
if fit_mean
CS -= C*S
SS = 1.0 - CC - S*S
CC -= C*C
else
SS = 1.0 - CC
end
τ = 0.5*atan2(2CS, CC - SS)/ω
# These quantities will be used below. See Press & Rybicki paper.
# NOTE: They can be computed in a faster way, for the time being we keep
# this simple and straightforward implementation.
cos_ωτ = cos(ω*τ)
sin_ωτ = sin(ω*τ)

# Now we can compute the power
YC_τ = YS_τ = CC_τ = nil
for i in eachindex(times)
ωt = ω*(times[i] - τ)
W = w[i]
c = cos(ωt)
s = sin(ωt)
YC_τ += W*y[i]*c
YS_τ += W*y[i]*s
CC_τ += W*c*c
end
if fit_mean
# "C_τ" and "S_τ" are computed following equation (7) of Press &
# Rybicki, this formula simply comes from angle difference
# trigonometric identities.
C_τ = C*cos_ωτ + S*sin_ωτ
S_τ = S*cos_ωτ - C*sin_ωτ
YC_τ -= Y*C_τ
YS_τ -= Y*S_τ
SS_τ = 1.0 - CC_τ - S_τ*S_τ
CC_τ -= C_τ*C_τ
YY -= Y*Y
else
SS_τ = 1.0 - CC_τ
end
P[n] = (abs2(YC_τ)/CC_τ + abs2(YS_τ)/SS_τ)/YY
end
y = signal
end

df = step(freqs)
N = length(freqs)
f0 = minimum(freqs)
#---------------------------------------------------------------------------
# 1. compute functions of the time-shift tau at each frequency
Ch, Sh = trig_sum(times, w .* y, df, N, f0, 1, oversampling, Mfft)
C2, S2 = trig_sum(times, w, df, N, f0, 2, oversampling, Mfft)
if fit_mean
C, S = trig_sum(times, w, df, N, f0, 1, oversampling, Mfft)
tan_2ωτ = (S2 .- 2.0 * S .* C) ./ (C2 .- (C .* C .- S .* S))
else
tan_2ωτ = S2 ./ C2
end
# This is what we're computing below; the straightforward way is slower and
# less stable, so we use trig identities instead
#
# ωτ = 0.5 * atan(tan_2ωτ)
# S2w, C2w = sin(2 * ωτ), cos(2 * ωτ)
# Sw, Cw = sin(ωτ), cos(ωτ)
C2w = 1./(sqrt(1.0 + tan_2ωτ .* tan_2ωτ))
S2w = tan_2ωτ .* C2w
Cw = sqrt(0.5 * (1.0 + C2w))
Sw = sqrt(0.5) .* sign(S2w) .* sqrt(1.0 - C2w)
#---------------------------------------------------------------------------
# 2. Compute the periodogram, following Zechmeister & Kurster
# and using tricks from Press & Rybicki.
YY = w⋅(y.^2)
YC = Ch .* Cw .+ Sh .* Sw
YS = Sh .* Cw .- Ch .* Sw
CC = 0.5 * (1.0 + C2 .* C2w .+ S2 .* S2w)
SS = 0.5 * (1.0 - C2 .* C2w .- S2 .* S2w)
if fit_mean
CC -= (C .* Cw .+ S .* Sw) .^ 2
SS -= (S .* Cw .- C .* Sw) .^ 2
end
P = (YC .* YC ./ CC .+ YS .* YS ./ SS)/YY
return P, freqs
end

Expand Down Expand Up @@ -409,12 +439,13 @@ function _lombscargle{R1<:Real,R2<:Real,R3<:Real,R4<:Real}(times::AbstractVector
maximum_frequency=maximum_frequency))
@assert length(times) == length(signal) == length(w)
if fit_mean || with_errors
P, f = _generalised_lombscargle(times, signal, w, frequencies,
center_data, fit_mean,
choose_use_fft(length(frequencies),
floatrange,
use_fft),
oversampling, Mfft)
if choose_use_fft(length(frequencies), floatrange, use_fft)
P, f = _press_rybicki(times, signal, w, frequencies, center_data,
fit_mean, oversampling, Mfft)
else
P, f = _generalised_lombscargle(times, signal, w, frequencies,
center_data, fit_mean)
end
else
P, f = _lombscargle_orig(times, signal, frequencies, center_data)
end
Expand Down
17 changes: 17 additions & 0 deletions src/extirpolation.jl
View file
Expand Up @@ -9,6 +9,15 @@
#
# License is MIT "Expat".
#
### Commentary:
#
# Utility functions used for implementing the fast algorithm proposed here:
# * Press, W. H., Rybicki, G. B. 1989, ApJ, 338, 277 (URL:
# http://dx.doi.org/10.1086/167197,
# Bibcode: http://adsabs.harvard.edu/abs/1989ApJ...338..277P)
# The actual implementation is adapted from Astropy project. See
# stats/lombscargle/implementations/utils.py in Astropy source code.
#
### Code:

function add_at!{N1<:Number,N2<:Number,N3<:Number}(arr::AbstractVector{N1},
Expand All @@ -28,11 +37,17 @@ function extirpolate{RE<:Real,NU<:Number}(X::AbstractVector{RE},
# Get the maximum of "X", `maximum' has a faster method for ranges.
N = trunc(Int, maximum(X) + 0.5*M + 1)
end
# Now use legendre polynomial weights to populate the results array; This is
# an efficient recursive implementation (See Press et al. 1989)
result = zeros(NU, N)
# first take care of the easy cases where x is an integer
integers = find(isinteger, x)
add_at!(result, trunc(Int, x[integers]), y[integers])
deleteat!(x, integers)
deleteat!(y, integers)
# For each remaining x, find the index describing the extirpolation range.
# i.e. ilo[i] < x[i] < ilo[i] + M with x[i] in the center, adjusted so that
# the limits are within the range 0...N
ilo = clamp(trunc(Int, x - div(M, 2)), 0, N - M)
numerator = y .* [prod(x[j] - ilo[j] - (0:(M - 1))) for j in eachindex(x)]
denominator = factorial(M - 1)
Expand All @@ -55,6 +70,8 @@ function trig_sum{R1<:Real,R2<:Real}(t::AbstractVector{R1},
df *= freq_factor
f0 *= freq_factor
@assert df > 0
# `ifft' can take arrays of any length in input, but it's faster when the
# length is exactly a power of 2.
Nfft = nextpow2(N * oversampling)
t0 = minimum(t)
if f0 > 0
Expand Down

0 comments on commit 0d178d1

Please sign in to comment.