Merge pull request #1 from Deduction42/extrapolations

Extrapolations

.gitignore

@@ -22,3 +22,4 @@ docs/site/
 # committed for packages, but should be committed for applications that require a static
 # environment.
 Manifest.toml
+.CondaPkg/

.vscode/settings.json

@@ -1,2 +1,3 @@
 {
+    "julia.environmentPath": "C:\\Users\\ruben\\.julia\\environments\\v1.10"
 }

ToDo.txt

@@ -5,4 +5,8 @@
     - Reconstruct local cumulative integrals to interpolate
         - Use the existing time_integral(ts::AbstractTimeSeries{T}, Δt::TimeInterval) implementation which is efficient
         - Consider adding an "indhint" argument to make it more efficient (with an indhint output)
-    - Modify the interpolated results by subtracting the first value
+    - Modify the interpolated results by subtracting the first value
+
+(3) Allow "merge" to make use of a splatable constructor function
+    - No constrcutor function should return a "Tuple"
+    - Tuple interpolation will need to broadcast math operators

src/TimeRecords.jl

@@ -7,6 +7,7 @@ module TimeRecords
         TimeInterval, 
         TimeSeries, 
         interpolate,
+        strictinterp,
         time_averages,
         time_integrals,
         cumulative_integral,

src/_TimeRecord.jl

@@ -28,20 +28,27 @@ Base.:+(Δt::TimeRecord, x::Real) = TimeRecord(Δt.t, Δt.v+x)
 Base.:-(Δt::TimeRecord, x::Real) = TimeRecord(Δt.t, Δt.v-x)
 Base.:*(Δt::TimeRecord, x::Real) = TimeRecord(Δt.t, Δt.v*x)
 Base.:/(Δt::TimeRecord, x::Real) = TimeRecord(Δt.t, Δt.v/x)
+Base.promote_rule(T1::Type{TimeRecord{R1}}, T2::Type{TimeRecord{R2}}) where {R1,R2} = TimeRecord{promote_rule(R1,R2)}
+Base.convert(::Type{TimeRecord{T}}, x::TimeRecord) where T = TimeRecord{T}(timestamp(x), convert(T,value(x)))
+Base.promote_typejoin(T1::Type{TimeRecord{R1}}, T2::Type{TimeRecord{R2}}) where {R1,R2} = TimeRecord{Base.promote_typejoin(R1,R2)}
+Base.typejoin(T1::Type{TimeRecord{R1}}, T2::Type{TimeRecord{R2}}) where {R1,R2} = TimeRecord{Base.typejoin(R1,R2)}
 
+"""
+Merge multiple time record with the same timestamp and apply the function "f" to the results
+"""
+function Base.merge(f::Union{Type,Function}, vtr::TimeRecord...)
+    mtr = merge(vtr...)
+    return TimeRecord(vtr[1].t, f(value(mtr)...))
+end 
 
 """
-Merge multiple time record with the same timestamp into a single static vector
+Merge multiple time records with the same timestamp as a tuple
 """
 function Base.merge(vtr::TimeRecord...)
     if !mapreduce(timestamp, isequal, vtr)
         error("Cannot merge time records for different timestamps")
     end
-
-    T = promote_type(recordtype.(vtr)...)
-    N = length(vtr)
-
-    return TimeRecord(vtr[1].t, SVector{N,T}(value.(vtr)...))
+    return TimeRecord(vtr[begin].t, value.(vtr))
 end
 
 

src/_TimeSeries.jl

@@ -14,6 +14,8 @@ function Base.getindex(ts::T, ind) where T <: AbstractTimeSeries
     return T(getindex(ts.records, ind), issorted=issorted(ind))
 end
 
+
+
 Base.setindex!(ts::AbstractTimeSeries, x, ind) = setindex!(ts.records, x, ind)
 Base.size(ts::AbstractTimeSeries)           = (length(ts.records),)
 Base.firstindex(ts::AbstractTimeSeries)     = 1
@@ -60,7 +62,47 @@ Constructs a time series from two vectors (unix timestamps, values)
 TimeSeries(v::AbstractVector{TimeRecord{T}}; issorted=false) where T = TimeSeries{T}(v, issorted=issorted)
 TimeSeries(t::AbstractVector{<:Real}, v::AbstractVector{T}; issorted=false) where T = TimeSeries{T}(TimeRecord{T}.(t, v), issorted=issorted)
 
+# =======================================================================================
+# Timeseries views
+# =======================================================================================
+"""
+View of a timeseries
+"""
+struct TimeSeriesView{T, P, I, LinIndex} <: AbstractTimeSeries{T}
+    records :: SubArray{TimeRecord{T}, 1, P, I, LinIndex}
+end
+
+Base.view(ts::AbstractTimeSeries, ind::Any) = error("View of AbstractTimeSeries can only be indexed by a UnitRange or AbstractVector{Bool}")
+Base.view(ts::AbstractTimeSeries, ind::UnitRange) = TimeSeriesView(view(ts.records, ind))
+Base.view(ts::AbstractTimeSeries, ind::AbstractVector{Bool}) = TimeSeriesView(view(ts.records, ind))
+
+
+# =======================================================================================
+# Merging functionality through extrapolation
+# =======================================================================================
+"""
+Merges a set of timeseries to a common set of timestamps through extrapolation
+Produces a StaticVector for each timestamp
+"""
+function Base.merge(t::AbstractVector{<:Real}, vts::AbstractTimeSeries...; order=1)
+    ts_extrap = map(ts->interpolate(ts,t, order=order), vts)
+    return _merge_records(ts_extrap...)
+end
+
+function Base.merge(f::Union{Function,Type}, t::AbstractVector{<:Real}, vts::AbstractTimeSeries...; order=1)
+    ts_extrap = map(ts->interpolate(ts,t, order=order), vts)
+    return _merge_records(f, ts_extrap...)
+end
 
+function _merge_records(f::Union{Function,Type}, uts::AbstractTimeSeries...)
+    return TimeSeries([merge(f, r...) for r in zip(uts...)])
+end
+
+function _merge_records(uts::AbstractTimeSeries...)
+    return TimeSeries([merge(r...) for r in zip(uts...)])
+end
+
+#=
 # =======================================================================================
 # Stateful timeseries
 # =======================================================================================
@@ -93,5 +135,5 @@ function take_next!(ts::StatefulTimeSeries)
     return current_value(ts)
 end
 
-
+=#
 

src/__assembly.jl

@@ -3,6 +3,7 @@ using Dates
 
 include(joinpath(@__DIR__, "_TimeRecord.jl"))
 include(joinpath(@__DIR__, "_TimeSeries.jl"))
+include(joinpath(@__DIR__, "find.jl"))
 include(joinpath(@__DIR__, "interpolate.jl"))
 include(joinpath(@__DIR__, "time_integral.jl"))
 include(joinpath(@__DIR__, "tabular_form.jl"))

src/find.jl

@@ -0,0 +1,145 @@
+"""
+findouter(ts::AbstractTimeSeries, Δt::TimeInterval, indhint::Integer)
+
+Finds the indices of the time series (ts) surrounding the time interval (Δt) or nearest (if Δt extends beyond (ts)) 
+"""
+function findouter(ts::AbstractTimeSeries, Δt::TimeInterval, indhint::Integer)
+
+    after_begin_excl(x::TimeRecord)  = Δt[begin] < timestamp(x)
+    after_end_incl(x::TimeRecord)    = Δt[end] <= timestamp(x)
+    before_end_excl(x::TimeRecord)   = timestamp(x) < Δt[end]
+    before_begin_incl(x::TimeRecord) = timestamp(x) <= Δt[begin]
+
+    ind0 = firstindex(ts)
+    indN = lastindex(ts)
+    indhint = clamp(indhint, ind0, indN)
+
+    if timestamp(ts[indhint]) < Δt[begin] #Hint occurs before interval, walk forward to find it
+        indL = something(findnext(after_begin_excl, ts, indhint), indN+1) - 1
+        indH = something(findnext(after_end_incl, ts, indL), indN)
+        return indL:indH
+
+    elseif Δt[end] < timestamp(ts[indhint]) #Hint occurs after interval, walk backward to find it
+        indH = something(findprev(before_end_excl, ts, indhint), ind0-1) + 1
+        indL = something(findprev(before_begin_incl, ts, indH), ind0)
+        return indL:indH
+
+    else #Hint occurs inside interval walk backward to find lower, and walk forward to find uppper
+        indL = something(findprev(before_begin_incl, ts, indhint), ind0)
+        indH = something(findnext(after_end_incl, ts, indhint), indN)
+        return indL:indH
+    end
+end
+
+findouter(ts::AbstractTimeSeries, Δt::TimeInterval) = findouter(ts, Δt, firstindex(ts))
+
+"""
+findouter(ts::AbstractTimeSeries, Δt::TimeInterval, indhint::Base.RefValue{<:Integer})
+
+Finds the indices of the time series (ts) surrounding the time interval (Δt) or nearest (if Δt extends beyond (ts)) 
+Stores upper limit of result inside indhint for future reference
+"""
+function findouter(ts::AbstractTimeSeries, Δt::TimeInterval, indhint::Base.RefValue{<:Integer})
+    inds = findouter(ts, Δt, indhint[])
+    if length(inds) > 0
+        indhint[] = inds[end]
+    end
+    return inds
+end
+
+"""
+findinner(ts::AbstractTimeSeries, Δt::TimeInterval, indhint)
+
+Finds the indices of the TimeSeries inside the time interval Δt::TimeInterval
+Stores upper limit of result inside indhint for future reference if indhint is a RefValue
+"""
+findinner(ts::AbstractTimeSeries, Δt::TimeInterval, indhint) = _findinner(findouter(ts, Δt, indhint))
+findinner(ts::AbstractTimeSeries, Δt::TimeInterval) = findinner(ts, Δt, firstindex(ts))
+
+_findinner(outer::UnitRange) = (outer[begin]+1):(outer[end]-1)
+
+
+
+"""
+findbounds(ts::AbstractTimeSeries, t::Real, indhint::Integer)
+
+Finds the index of the TimeRecord before and after t::Real; indhint is the first index searched for
+"""
+function findbounds(ts::AbstractTimeSeries, t::Real, indhint::Integer)
+    earlier_than_t(x::TimeRecord) = timestamp(x) <= t
+    later_than_t(x::TimeRecord)   = timestamp(x) >= t
+
+    ind0 = firstindex(ts)
+    indN = lastindex(ts)
+    indhint = clamp(indhint, ind0, indN)
+
+    if earlier_than_t(ts[indhint]) #Walk forward in time if indhint record occurs earlier than t
+        indH = findnext(later_than_t, ts, indhint)
+        indL = something(indH, indN+1) - 1
+        indL = max(indL, ind0)
+        return (indL, indH)
+
+    else #Walk backwards in time if indhint reccord occurs later than t
+        indL = findprev(earlier_than_t, ts, indhint)
+        indH = something(indL, firstindex(ts)-1) + 1
+        indH = min(indH, indN)
+        return (indL, indH)
+    end
+end
+
+"""
+findbounds(ts::AbstractTimeSeries, t::Real, indhint::RefValue{<:Integer})
+
+Finds the index of the TimeRecord before and after t::Real; indhint is the first index searched for
+The upper limit of the boundary is saved in indhint (unless it's nothing, then the lower boundary is saved)
+"""
+function findbounds(ts::AbstractTimeSeries, t::Real, indhint::Base.RefValue{<:Integer})
+    (lb, ub) = findbounds(ts, t, indhint[])
+    indhint[] = something(ub, lb)
+    return (lb, ub)
+end
+
+
+"""
+findbounds(ts::AbstractTimeSeries, t::Real)
+
+Finds bounding indices for timeseries (ts) at time (t) using a bisection method
+"""
+function findbounds(ts::AbstractTimeSeries, t::Real)
+    (lb, ub) = (firstindex(ts), lastindex(ts))
+
+    if t < timestamp(ts[lb])
+        return (nothing, lb)
+    elseif timestamp(ts[ub]) < t
+        return (ub, nothing)
+    end
+
+    while (ub-lb) > 1
+        mb = ceil(Int64, 0.5*(lb+ub))
+        if timestamp(ts[mb]) < t
+            (lb, ub) = (mb, ub)
+        else
+            (lb, ub) = (lb, mb)
+        end
+    end
+
+    return (lb, ub)
+end
+
+findbounds(ts::AbstractTimeSeries, t::Real, indhint::Nothing) = findbounds(ts, t)
+
+"""
+findnearest(ts::AbstractTimeSeries, t::Real, indhint::RefValue{<:Integer})
+
+Behaves like findbounds except that it always returns integer boundaries within the timeseries bounds
+"""
+findnearest(ts::AbstractTimeSeries, t::Real, indhint) = clampbounds(findbounds(ts, t, indhint))
+
+"""
+findnearest(Tuple{Union{Integer,Nothing}, Union{Integer,Nothing}})
+
+Clamps the result of findbounds so that Nothing is not inside the results
+"""
+clampbounds(t::Tuple{Nothing, <:Integer})   = (t[2],t[2])
+clampbounds(t::Tuple{<:Integer, Nothing})   = (t[1],t[1])
+clampbounds(t::Tuple{<:Integer, <:Integer}) = t