Skip to content

Commit

Permalink
🔧 Use Julian Day as internal date representation (#3)
Browse files Browse the repository at this point in the history
This commit added compatibility with automatic differentiation packages.
  • Loading branch information
jmurphy6895 committed Aug 28, 2023
1 parent e266cfa commit d2e767d
Show file tree
Hide file tree
Showing 10 changed files with 244 additions and 191 deletions.
6 changes: 5 additions & 1 deletion Project.toml
Original file line number Diff line number Diff line change
Expand Up @@ -21,6 +21,10 @@ julia = "1.6"
[extras]
Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267"
Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"


[targets]
test = ["Test", "Logging"]
test = ["Test", "Logging", "ForwardDiff", "ReverseDiff", "Zygote"]
1 change: 1 addition & 0 deletions docs/src/man/api.md
Original file line number Diff line number Diff line change
Expand Up @@ -79,6 +79,7 @@ signature:

```julia
SpaceIndices.space_index(::Val{:index}, instant::DateTime; kwargs...) -> Number
SpaceIndices.space_index(::Val{:index}, jd::Number; kwargs...) -> Number
```

where the space `index` for the `instant` will be returned.
Expand Down
1 change: 1 addition & 0 deletions src/API.md
Original file line number Diff line number Diff line change
Expand Up @@ -73,6 +73,7 @@ signature:

```julia
SpaceIndices.space_index(::Val{:index}, instant::DateTime; kwargs...) -> Number
SpaceIndices.space_index(::Val{:index}, jd::Number; kwargs...) -> Number
```

where the space `index` for the `instant` will be returned.
Expand Down
6 changes: 3 additions & 3 deletions src/api.jl
Original file line number Diff line number Diff line change
Expand Up @@ -101,16 +101,16 @@ urls

"""
space_index(::Val{:index}, jd::Number; kwargs...) -> Number
space_index(::Val{:index}, instant::DateTime; kwargs...) -> Number
space_index(::Val{:index}, date::DateTime; kwargs...) -> Number
Get the space `index` for the Julian day `jd` or the `instant`. The latter must be an object
of type `DateTime`. `kwargs...` can be used to pass additional configuration for the space
index.
"""
space_index

function space_index(index::Val, jd::Number)
return space_index(index, julian2datetime(jd))
function space_index(index::Val, date::DateTime)
return space_index(index, datetime2julian(date))
end

"""
Expand Down
40 changes: 22 additions & 18 deletions src/interpolations.jl
Original file line number Diff line number Diff line change
Expand Up @@ -11,21 +11,25 @@
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

"""
constant_interpolation(knots::Vector{Date}, values::AbstractVector, x::Date) -> eltype(values)
constant_interpolation(knots::AbstractVector, values::AbstractVector, x) -> eltype(values)
Perform a constant interpolation at `x` of `values` evaluated at `knots`. The interpolation
returns `value(knots[k-1])` in which `knots[k-1] <= x < knots[k]`.
"""
function constant_interpolation(knots::AbstractVector{Tk}, values::AbstractVector, x::Tk) where Tk
function constant_interpolation(knots::AbstractVector, values::AbstractVector, x)
# First, we need to verify if `x` is inside the domain.
num_knots = length(knots)
knots_beg = first(knots)
knots_end = last(knots)

if !(knots_beg <= x <= knots_end)

x_dt = julian2datetime(x)
knots_beg_dt = julian2datetime(knots_beg)
knots_end_dt = julian2datetime(knots_end)

throw(ArgumentError("""
There is no available data for x = $(x)!
The available interval is: x ∈ [$knots_beg, $knots_end]."""
There is no available data for x = $(x_dt)!
The available interval is: x ∈ [$knots_beg_dt, $knots_end_dt]."""
))
end

Expand All @@ -37,46 +41,46 @@ function constant_interpolation(knots::AbstractVector{Tk}, values::AbstractVecto
end

"""
linear_interpolation(knots::AbstractVector{Tk}, values::AbstractVector{Tv}, x::Tk) where {Tk, Tv}
linear_interpolation(knots::AbstractVector, values::AbstractVector, x)
Perform a linear interpolation at `x` of `values` evaluated at `knots`.
"""
function linear_interpolation(knots::AbstractVector{Tk}, values::AbstractVector{Tv}, x::Tk) where {Tk, Tv}
function linear_interpolation(knots::AbstractVector, values::AbstractVector, x)
# First, we need to verify if `x` is inside the domain.
num_knots = length(knots)
knots_beg = first(knots)
knots_end = last(knots)

if !(knots_beg <= x <= knots_end)
x_dt = julian2datetime(x)
knots_beg_dt = julian2datetime(knots_beg)
knots_end_dt = julian2datetime(knots_end)

throw(ArgumentError("""
There is no available data for x = $(x)!
The available interval is: x ∈ [$knots_beg -- $knots_end]."""
There is no available data for x = $(x_dt)!
The available interval is: x ∈ [$knots_beg_dt -- $knots_end_dt]."""
))
end

# Type of the returned value.
T = float(Tv)

# Find the vector index related to the request interval using binary search. We can
# apply this algorithm because we assume that `knots` are unique and increasing.
id = _binary_search(knots, x)

# If we are at the knot precisely, just return it.
if x == knots[id]
return Tv(values[id])
return values[id]

else
# Here, we need to perform the interpolation using the adjacent knots.
x₀ = knots[id]
x₁ = knots[id + 1]

y₀ = T(values[id])
y₁ = T(values[id + 1])
y₀ = values[id]
y₁ = values[id + 1]

Δy = y₁ - y₀
Δx = x₁ - x₀

y = y₀ + Δy * T((x - x₀) / Δx)
y = y₀ + Δy * (x - x₀) / Δx

return y
end
Expand All @@ -88,7 +92,7 @@ end

# Perform a interval binary search of `x` in `v`. It means that this function returns `k`
# such that `v[k] <= x < v[k + 1]`.
function _binary_search(v::AbstractVector{T}, x::T) where T
function _binary_search(v::AbstractVector, x)
num_elements = length(v)
low = 1
high = num_elements
Expand Down
Loading

0 comments on commit d2e767d

Please sign in to comment.