Skip to content

Commit

Permalink
Store fftgrid into the plan type for PR method
Browse files Browse the repository at this point in the history
  • Loading branch information
@giordano
giordano committed Nov 15, 2017
1 parent 2112130 commit 5443a27
Show file tree
Hide file tree
Showing 5 changed files with 42 additions and 38 deletions.
1 change: 0 additions & 1 deletion src/LombScargle.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,6 @@ __precompile__()
module LombScargle

using FFTW, Measurements
importall FFTW

export lombscargle

Expand Down
4 changes: 2 additions & 2 deletions src/extirpolation.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -56,7 +56,7 @@ function extirpolate!(result, x::AbstractVector{<:Real}, y::AbstractVector{NU},
return result
end

function trig_sum!(grid, bfft_vec, bfft_plan, t::AbstractVector{<:Real},
function trig_sum!(grid, fftgrid, bfft_vec, bfft_plan, t::AbstractVector{<:Real},
h::AbstractVector{<:Real}, df::Real, N::Integer,
Nfft::Integer, t0::Real, f0::Real=0.0,
freq_factor::Integer=1, Mfft::Integer=4)
Expand All @@ -70,7 +70,7 @@ function trig_sum!(grid, bfft_vec, bfft_plan, t::AbstractVector{<:Real},
tnorm = mod.(((t .- t0) .* Nfft .* df), Nfft)
extirpolate!(grid, tnorm, H, Nfft, Mfft)
A_mul_B!(bfft_vec, bfft_plan, grid)
fftgrid = @view(bfft_vec[1:N])
fftgrid .= @view(bfft_vec[1:N])
if t0 != 0
fftgrid .*= cis.(t0 .* (f0 .+ df .* (0:(N - 1))) .* 2 .* pi)
end
Expand Down
12 changes: 7 additions & 5 deletions src/planning.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -64,19 +64,21 @@ function _plan(times::AbstractVector{<:Real}, signal::AbstractVector{R1}, sumw::
y = signal
end
YY = w ⋅ (y .^ 2)
Nfft = nextpow2(length(frequencies) * oversampling)
N = length(frequencies)
Nfft = nextpow2(N * oversampling)
T = promote_type(float(R1), float(R2))
bfft_vect = Vector{Complex{T}}(Nfft)
bfft_grid = Vector{Complex{T}}(Nfft)
fftgrid = Vector{Complex{T}}(N)
bfft_plan = FFTW.plan_bfft(bfft_vect, flags = flags, timelimit = timelimit)
if fit_mean
return FastGLSPlan_fit_mean(times, signal, frequencies, sumw, w, y, YY,
bfft_vect, bfft_grid, bfft_plan, Mfft,
noise, normalization, Vector{T}(length(frequencies)))
fftgrid, bfft_vect, bfft_grid, bfft_plan, Mfft,
noise, normalization, Vector{T}(N))
else
return FastGLSPlan(times, signal, frequencies, sumw, w, y, YY,
bfft_vect, bfft_grid, bfft_plan, Mfft,
noise, normalization, Vector{T}(length(frequencies)))
fftgrid, bfft_vect, bfft_grid, bfft_plan, Mfft,
noise, normalization, Vector{T}(N))
end
else
return _plan_no_fast(times, signal, sumw, w, frequencies, with_errors,
Expand Down
60 changes: 31 additions & 29 deletions src/press-rybicki.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,38 +12,40 @@
#
### Code:

struct FastGLSPlan{T,A,B<:AbstractVector{T},C,D,E,F,G,H,I,J,K} <: PeriodogramPlan
struct FastGLSPlan{T,A,B<:AbstractVector{T},C,D,E,F,G,H,I,J,K,L} <: PeriodogramPlan
times::A
signal::B
freq::C
sumw::D
w::E
y::B
YY::T
bfft_vect::F
bfft_grid::G
bfft_plan::H
Mfft::I
noise::J
fftgrid::F
bfft_vect::G
bfft_grid::H
bfft_plan::I
Mfft::J
noise::K
norm::Symbol
P::K
P::L
end

struct FastGLSPlan_fit_mean{T,A,B<:AbstractVector{T},C,D,E,F,G,H,I,J,K} <: PeriodogramPlan
struct FastGLSPlan_fit_mean{T,A,B<:AbstractVector{T},C,D,E,F,G,H,I,J,K,L} <: PeriodogramPlan
times::A
signal::B
freq::C
sumw::D
w::E
y::B
YY::T
bfft_vect::F
bfft_grid::G
bfft_plan::H
Mfft::I
noise::J
fftgrid::F
bfft_vect::G
bfft_grid::H
bfft_plan::I
Mfft::J
noise::K
norm::Symbol
P::K
P::L
end

include("extirpolation.jl")
Expand All @@ -58,13 +60,13 @@ include("extirpolation.jl")

function _press_rybicki!(P, times::AbstractVector{<:Real}, y::AbstractVector{<:Real},
w::AbstractVector{<:Real}, t0, df, N , f0, YY::Real,
bfft_vec::AbstractVector{Complex{T}},
fftgrid, bfft_vec::AbstractVector{Complex{T}},
grid, plan, Nfft::Integer, Mfft::Integer) where {T<:AbstractFloat}
#---------------------------------------------------------------------------
# 1. compute functions of the time-shift tau at each frequency
Ch, Sh = trig_sum!(grid, bfft_vec, plan, times, w .* y, df, N, Nfft, t0, f0, 1, Mfft)
C2, S2 = trig_sum!(grid, bfft_vec, plan, times, w, df, N, Nfft, t0, f0, 2, Mfft)
#-----------------------------------------------------------------------
Ch, Sh = trig_sum!(grid, fftgrid, bfft_vec, plan, times, w .* y, df, N, Nfft, t0, f0, 1, Mfft)
C2, S2 = trig_sum!(grid, fftgrid, bfft_vec, plan, times, w, df, N, Nfft, t0, f0, 2, Mfft)
#---------------------------------------------------------------------------
# 2. Compute the periodogram, following Zechmeister & Kurster
# and using tricks from Press & Rybicki.
tan_2ωτ = S2 ./ C2
Expand All @@ -78,15 +80,15 @@ end

function _press_rybicki_fit_mean!(P, times::AbstractVector{<:Real},
y::AbstractVector{<:Real}, w::AbstractVector{<:Real},
t0, df, N , f0, YY::Real,
t0, df, N , f0, YY::Real, fftgrid,
bfft_vec::AbstractVector{Complex{T}}, grid, plan,
Nfft::Integer, Mfft::Integer) where {T<:AbstractFloat}
#---------------------------------------------------------------------------
# 1. compute functions of the time-shift tau at each frequency
Ch, Sh = trig_sum!(grid, bfft_vec, plan, times, w .* y, df, N, Nfft, t0, f0, 1, Mfft)
C2, S2 = trig_sum!(grid, bfft_vec, plan, times, w, df, N, Nfft, t0, f0, 2, Mfft)
C, S = trig_sum!(grid, bfft_vec, plan, times, w, df, N, Nfft, t0, f0, 1, Mfft)
#-----------------------------------------------------------------------
Ch, Sh = trig_sum!(grid, fftgrid, bfft_vec, plan, times, w .* y, df, N, Nfft, t0, f0, 1, Mfft)
C2, S2 = trig_sum!(grid, fftgrid, bfft_vec, plan, times, w, df, N, Nfft, t0, f0, 2, Mfft)
C, S = trig_sum!(grid, fftgrid, bfft_vec, plan, times, w, df, N, Nfft, t0, f0, 1, Mfft)
#---------------------------------------------------------------------------
# 2. Compute the periodogram, following Zechmeister & Kurster
# and using tricks from Press & Rybicki.
tan_2ωτ = (S2 .- 2 .* S .* C) ./ (C2 .- (C .* C .- S .* S))
Expand All @@ -101,18 +103,18 @@ function _press_rybicki_fit_mean!(P, times::AbstractVector{<:Real},
end

_periodogram!(p::FastGLSPlan) =
_press_rybicki!(p.P, p.times, p.y, p.w, minimum(p.times), step(p.freq),
length(p.freq), minimum(p.freq), p.YY, p.bfft_vect, p.bfft_grid,
p.bfft_plan, length(p.bfft_vect), p.Mfft)
_press_rybicki!(p.P, p.times, p.y, p.w, minimum(p.times), step(p.freq), length(p.freq),
minimum(p.freq), p.YY, p.fftgrid, p.bfft_vect, p.bfft_grid, p.bfft_plan,
length(p.bfft_vect), p.Mfft)
_periodogram!(times, p::FastGLSPlan) =
_press_rybicki!(p.P, times, p.y, p.w, minimum(p.times), step(p.freq),
length(p.freq), minimum(p.freq), p.YY, p.bfft_vect,
length(p.freq), minimum(p.freq), p.YY, p.fftgrid, p.bfft_vect,
p.bfft_grid, p.bfft_plan, length(p.bfft_vect), p.Mfft)
_periodogram!(p::FastGLSPlan_fit_mean) =
_press_rybicki_fit_mean!(p.P, p.times, p.y, p.w, minimum(p.times), step(p.freq),
length(p.freq), minimum(p.freq), p.YY, p.bfft_vect,
length(p.freq), minimum(p.freq), p.YY, p.fftgrid, p.bfft_vect,
p.bfft_grid, p.bfft_plan, length(p.bfft_vect), p.Mfft)
_periodogram!(times, p::FastGLSPlan_fit_mean) =
_press_rybicki_fit_mean!(p.P, times, p.y, p.w, minimum(p.times), step(p.freq),
length(p.freq), minimum(p.freq), p.YY, p.bfft_vect,
length(p.freq), minimum(p.freq), p.YY, p.fftgrid, p.bfft_vect,
p.bfft_grid, p.bfft_plan, length(p.bfft_vect), p.Mfft)
3 changes: 2 additions & 1 deletion test/runtests.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -170,10 +170,11 @@ end
@testset "LombScargle.trig_sum!" begin
N = 10
Nfft = nextpow2(5N)
fftgrid = Vector{Complex{Float64}}(N)
bfft_vec = Vector{Complex{Float64}}(Nfft)
p = plan_bfft(bfft_vec)
grid = similar(bfft_vec)
C, S = LombScargle.trig_sum!(grid, bfft_vec, p, x, y, 1, N, Nfft, minimum(x))
C, S = LombScargle.trig_sum!(grid, fftgrid, bfft_vec, p, x, y, 1, N, Nfft, minimum(x))
@test S ≈ [0.0,0.3753570125888358,0.08163980192703546,-0.10139634351774979,-0.4334223744905633,-2.7843373311769777,0.32405810159838055,0.05729663600471602,-0.13191736591325876,-0.5295781583202946]
@test C ≈ [8.708141477890015,-0.5402668064176129,-0.37460815054027985,-0.3793457539084364,-0.5972351546196192,14.612204307982497,-0.5020253140297526,-0.37724493022381034,-0.394096831764578,-0.6828241623474718]
end
Expand Down

0 comments on commit 5443a27

Please sign in to comment.