Skip to content

Commit

Permalink
Fix type instabilities in extirpolate
Browse files Browse the repository at this point in the history
Do not assert lengths of X and Y are the same, this is already done in
`_press_rybicki'; make N mandatory and do not change its value; make
`denominator' always an AbstractFloat.
  • Loading branch information
giordano committed Dec 15, 2016
1 parent 2908781 commit 2e07d53
Show file tree
Hide file tree
Showing 2 changed files with 6 additions and 10 deletions.
12 changes: 4 additions & 8 deletions src/extirpolation.jl
Original file line number Diff line number Diff line change
Expand Up @@ -30,13 +30,9 @@ end

function extirpolate{RE<:Real,NU<:Number}(X::AbstractVector{RE},
Y::AbstractVector{NU},
N::Integer=0, M::Integer=4)
@assert length(X) == length(Y)
x, y = collect(X), collect(Y)
if N <= 0
# Get the maximum of "X", `maximum' has a faster method for ranges.
N = trunc(Int, maximum(X) + 0.5*M + 1)
end
N::Integer, M::Integer=4)
x = collect(X)
y = copy(Y)
# 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)
Expand All @@ -50,7 +46,7 @@ function extirpolate{RE<:Real,NU<:Number}(X::AbstractVector{RE},
# 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)
denominator = float(factorial(M - 1))
@inbounds for j in 0:(M - 1)
if j > 0
denominator *= j/(j - M)
Expand Down
4 changes: 2 additions & 2 deletions test/runtests.jl
Original file line number Diff line number Diff line change
Expand Up @@ -129,8 +129,8 @@ LombScargle.add_at!(a, [3, 1, 3, 1, 2], 1:5)
# Test extirpolation function
x = linspace(0, 10)
y = sin(x)
@test_approx_eq LombScargle.extirpolate(x, y) [0.39537718210649553,3.979484140636793,4.833090108345013,0.506805556164743,-3.828112427525919,-4.748341359084166,-1.3022050566901917,3.3367666084342256,5.070478111668922,1.291245296032218,-0.8264466821981216,0.0,0.0]
@test_approx_eq LombScargle.extirpolate(x, y, 11) LombScargle.extirpolate(x, y)[1:11]
@test_approx_eq LombScargle.extirpolate(x, y, 13) [0.39537718210649553,3.979484140636793,4.833090108345013,0.506805556164743,-3.828112427525919,-4.748341359084166,-1.3022050566901917,3.3367666084342256,5.070478111668922,1.291245296032218,-0.8264466821981216,0.0,0.0]
@test_approx_eq LombScargle.extirpolate(x, y, 11) LombScargle.extirpolate(x, y, 13)[1:11]

# Test trig_sum
C, S = LombScargle.trig_sum(x, y, 1, 10)
Expand Down

0 comments on commit 2e07d53

Please sign in to comment.