Skip to content

Commit

Permalink
Use dot syntax for vectorial function whenever it is possible
Browse files Browse the repository at this point in the history
Now the fast method should eat even less memory than before (and thanks
to the use of in-place operators the situation had already much
improved).
  • Loading branch information
giordano committed Dec 24, 2016
1 parent fb850b9 commit c5c1b06
Show file tree
Hide file tree
Showing 2 changed files with 15 additions and 16 deletions.
16 changes: 8 additions & 8 deletions src/LombScargle.jl
Original file line number Diff line number Diff line change
Expand Up @@ -216,7 +216,7 @@ function _press_rybicki{R1<:Real,R2<:Real,R3<:Real,R4<:Real}(times::AbstractVect
if fit_mean
C, S = trig_sum!(grid, ifft_vec, plan, times, w, df, N, Nfft, f0,
1, Mfft)
tan_2ωτ = (S2 .- 2.0 * S .* C) ./ (C2 .- (C .* C .- S .* S))
tan_2ωτ = (S2 .- 2 .* S .* C) ./ (C2 .- (C .* C .- S .* S))
else
tan_2ωτ = S2 ./ C2
end
Expand All @@ -225,21 +225,21 @@ function _press_rybicki{R1<:Real,R2<:Real,R3<:Real,R4<:Real}(times::AbstractVect
# ωτ = 0.5 * atan(tan_2ωτ)
# S2w, C2w = sin(2 * ωτ), cos(2 * ωτ)
# Sw, Cw = sin(ωτ), cos(ωτ)
C2w = 1./(sqrt.(1.0 + tan_2ωτ .* tan_2ωτ))
C2w = 1 ./ (sqrt.(1 .+ 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)
Cw = sqrt.((1 .+ C2w) ./ 2)
Sw = sign(S2w) .* sqrt.((1 .- C2w) ./ 2)
#---------------------------------------------------------------------------
# 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)
CC = (1 .+ C2 .* C2w .+ S2 .* S2w) ./ 2
SS = (1 .- C2 .* C2w .- S2 .* S2w) ./ 2
if fit_mean
CC -= (C .* Cw .+ S .* Sw) .^ 2
SS -= (S .* Cw .- C .* Sw) .^ 2
CC .-= (C .* Cw .+ S .* Sw) .^ 2
SS .-= (S .* Cw .- C .* Sw) .^ 2
end
return (YC .* YC ./ CC .+ YS .* YS ./ SS)/YY
end
Expand Down
15 changes: 7 additions & 8 deletions src/extirpolation.jl
Original file line number Diff line number Diff line change
Expand Up @@ -39,13 +39,13 @@ function extirpolate!{RE<:Real,NU<:Number}(result,
fill!(result, zero(NU))
# first take care of the easy cases where x is an integer
integers = find(isinteger, x)
add_at!(result, mod(trunc(Int, x[integers]), N) + 1, y[integers])
add_at!(result, mod(trunc(Int, x[integers]), N) .+ 1, 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)
ilo = clamp(trunc(Int, x .- div(M, 2)), 0, N - M)
v = collect(0:(M - 1))
numerator = y .* [prod(x[j] - ilo[j] - v) for j in eachindex(x)]
denominator = float(factorial(M - 1))
Expand All @@ -54,7 +54,7 @@ function extirpolate!{RE<:Real,NU<:Number}(result,
denominator *= j/(j - M)
end
ind = ilo + (M - j)
add_at!(result, ind, numerator ./ (denominator * (x .- ind + 1)))
add_at!(result, ind, numerator ./ (denominator * (x .- ind .+ 1)))
end
return result
end
Expand All @@ -71,17 +71,16 @@ function trig_sum!{R1<:Real,R2<:Real}(grid, ifft_vec, ifft_plan,
@assert df > 0
t0 = minimum(t)
if f0 > 0
H = h .* exp.(2im * pi * f0 * (t - t0))
H = h .* exp.(2im .* pi .* f0 .* (t .- t0))
else
H = complex(h)
end
tnorm = mod(((t - t0) * Nfft * df), Nfft)
tnorm = mod(((t .- t0) .* Nfft .* df), Nfft)
extirpolate!(grid, tnorm, H, Nfft, Mfft)
A_mul_B!(ifft_vec, ifft_plan, grid)
fftgrid = Nfft * ifft_vec[1:N]
fftgrid = Nfft .* ifft_vec[1:N]
if t0 != 0
f = f0 + df * (0:N - 1)
fftgrid .*= exp.(2im * pi * t0 * f)
fftgrid .*= exp.(2im .* pi .* t0 .* (f0 .+ df .* (0:N - 1)))
end
C = real(fftgrid)
S = imag(fftgrid)
Expand Down

0 comments on commit c5c1b06

Please sign in to comment.