Fix TT2000 leap-second read; promote DateTime/Timestamp to the CDF epoch types

CHANGELOG.md

@@ -2,6 +2,21 @@
 
 ## [Unreleased]
 
+### Changed
+
+- **Breaking**: `Date`, `DateTime` and `Timestamp` promote to the CDF epoch type (`TT2000`, `Epoch16`, `Epoch`) instead of the epoch promoting to `DateTime`, so mixed comparisons and differences are exact: to the nanosecond, and across leap seconds for `TT2000` ([#74]).
+- **Breaking**: `DateTime(::Epoch16)` floors to the millisecond instead of rounding ([#74]).
+
+### Added
+
+- `Timestamp{P}(x)` conversions (Durations.jl) for all CDF epochs
+- String constructors keep every fractional digit: nanoseconds for `TT2000`, picoseconds for `Epoch16` ([#74]).
+- Comparison between different CDF epoch types, `Epoch16` ordering and `Epoch16 ± Period`, `Date`/`Time` and day-based accessors ([#74]).
+
+### Fixed
+
+- `FILLVAL` display: `TT2000` uses the ISTP fill value (`typemin(Int64)`) rather than `9999`, and `Epoch16` fill no longer throws ([#74]).
+
 ## [0.2.0] - 2026-06-11
 
 ### Changed
@@ -11,4 +26,5 @@
 [Unreleased]: https://github.com/JuliaSpacePhysics/CommonDataFormat.jl/compare/v0.2.0...HEAD
 [0.2.0]: https://github.com/JuliaSpacePhysics/CommonDataFormat.jl/releases/tag/v0.2.0
 
+[#74]: https://github.com/JuliaSpacePhysics/CommonDataFormat.jl/pull/74
 [#46]: https://github.com/JuliaSpacePhysics/CommonDataFormat.jl/pull/46

Project.toml

@@ -9,6 +9,7 @@ projects = ["benchmark", "docs", "test"]
 [deps]
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 DiskArrays = "3c3547ce-8d99-4f5e-a174-61eb10b00ae3"
+Durations = "41ac09f5-0a95-4810-aa18-c88bf48ee130"
 LibDeflate = "9255714d-24a7-4b30-8ea3-d46a97f7e13b"
 Mmap = "a63ad114-7e13-5084-954f-fe012c677804"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
@@ -23,6 +24,7 @@ CommonDataFormatCommonDataModelExt = ["CommonDataModel"]
 CommonDataModel = "0.4"
 Dates = "1"
 DiskArrays = "0.4"
+Durations = "1.4"
 LibDeflate = "0.4.3"
 Mmap = "1"
 PrecompileTools = "1"

benchmark/Project.toml

@@ -1,6 +1,7 @@
 [deps]
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 CommonDataFormat = "a9737db6-c05c-4e48-868b-6bc41491d9d9"
+Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
 
 [compat]

benchmark/benchmarks.jl

@@ -1,5 +1,6 @@
 using BenchmarkTools
 using CommonDataFormat
+using Dates
 using Downloads
 
 const SUITE = BenchmarkGroup()
@@ -34,3 +35,17 @@ let ds = CDFDataset(MMS_FILE), var = ds["mms1_scm_acb_gse_scsrvy_srvy_l2"]
     g["sum_var_access"] = @benchmarkable sum($var[:, 100:100000])
     g["full_load"] = @benchmarkable full_load(MMS_FILE)
 end
+
+# Inputs stay within what every epoch type could always represent (millisecond strings, modern
+# dates) so the suite also runs on older revisions for comparison.
+let tt = TT2000(DateTime(2020)) .+ Millisecond.(1:10_000), dt = DateTime.(tt),
+        e16 = Epoch16.(dt), strs = string.(dt)
+    g = SUITE["epochs"] = BenchmarkGroup()
+    g["TT2000_to_DateTime"] = @benchmarkable DateTime.($tt)
+    g["DateTime_to_TT2000"] = @benchmarkable TT2000.($dt)
+    g["TT2000_vs_DateTime"] = @benchmarkable count(<($(dt[end ÷ 2])), $tt)
+    g["Epoch16_vs_DateTime"] = @benchmarkable count(==($(dt[end ÷ 2])), $e16)
+    g["TT2000_parse"] = @benchmarkable TT2000.($strs)
+    g["Epoch16_parse"] = @benchmarkable Epoch16.($strs)
+    g["TT2000_string"] = @benchmarkable foreach(x -> print(devnull, x), $tt)
+end

src/epochs.jl

@@ -7,6 +7,7 @@
 
 import Base: promote_rule, -, +
 using Dates: value, toms, tons
+using Durations: Timestamp
 
 include("leap_second.jl")
 
@@ -63,52 +64,70 @@ end
 
 TT2000(instant::Int64) = TT2000(Nanosecond(instant))
 
-TT2000(dt::TimeType) = convert(TT2000, dt)
-Epoch(dt::TimeType) = convert(Epoch, dt)
-
-fillvalue(::Epoch) = -1.0e31
-fillvalue(::Epoch16) = -1.0e31
-fillvalue(::TT2000) = 9999
+# ISTP FILLVAL
+_isfill(x::Epoch) = x.instant == -1.0e31
+_isfill(x::Epoch16) = x.seconds == -1.0e31
+_isfill(x::TT2000) = x.instant.value == typemin(Int64)
 
 (-)(epoch::Epoch, other::Epoch) = Millisecond(round(Int64, epoch.instant - other.instant))
 (-)(epoch::Epoch16, other::Epoch16) = Picosecond((epoch.seconds - other.seconds) * 1.0e12 + epoch.picoseconds - other.picoseconds)
 (+)(tt2000::TT2000, other::Period) = TT2000(value(tt2000) + tons(other))
 (-)(tt2000::TT2000, other::Period) = TT2000(value(tt2000) - tons(other))
 (+)(epoch::Epoch, other::Period) = Epoch(value(epoch) + toms(other))
 (-)(epoch::Epoch, other::Period) = Epoch(value(epoch) - toms(other))
-
-# Conversion to DateTime
-function Dates.DateTime(epoch::Epoch)
-    return DateTime(0) + Millisecond(round(Int64, epoch.instant))
+# Whole seconds are split off first: a Float64 picosecond count is only exact within ~2.5 hours.
+function (+)(epoch::Epoch16, other::Period)
+    s, ns = fldmod(tons(other), 1_000_000_000)
+    carry, ps = fldmod(epoch.picoseconds + 1.0e3 * ns, 1.0e12)
+    return Epoch16(epoch.seconds + s + carry, ps)
 end
+(-)(epoch::Epoch16, other::Period) = epoch + (-other)
 
-function Dates.DateTime(epoch::Epoch16)
-    s_since_unix = epoch.seconds - EPOCH_OFFSET_SECONDS
-    total_ms = s_since_unix * 1.0e3 + epoch.picoseconds / 1.0e9
-    return DateTime(1970) + Millisecond(round(Int64, total_ms))
-end
+Base.isless(x::Epoch16, y::Epoch16) = isless((x.seconds, x.picoseconds), (y.seconds, y.picoseconds))
+
+# UTC nanoseconds since the Unix epoch. Int128: Epoch16 spans years 0-9999 and TT2000 reaches
+# 2292, both beyond an Int64 nanosecond count.
+_unix_ns(x::TT2000) = utc_from_tt(Int128(x.instant.value) + TT2000_OFFSET)
+_unix_ns(x::Epoch16) = (Int128(x.seconds) - Int128(EPOCH_OFFSET_SECONDS)) * 1_000_000_000 + fld(Int128(x.picoseconds), 1000)
+_unix_ns(dt::Timestamp{P}) where {P} = Int128(value(dt)) * tons(P(1))
+_unix_ns(dt::TimeType) = (Int128(value(DateTime(dt))) - Dates.UNIXEPOCH) * 1_000_000
+
+# Int128 division is a slow library call, and nearly every instant fits an Int64.
+_fld(ns::Int128, d::Int64) = typemin(Int64) <= ns <= typemax(Int64) ? fld(ns % Int64, d) : Int64(fld(ns, d))
+
+# Floored to the millisecond, so calendar fields never run ahead of the instant.
+Dates.DateTime(x::Union{TT2000, Epoch16}) = DateTime(Dates.UTM(_fld(_unix_ns(x), 1_000_000) + Dates.UNIXEPOCH))
+Dates.DateTime(epoch::Epoch) = DateTime(0) + Millisecond(round(Int64, epoch.instant))
+
+Timestamp{P}(x::CDFDateTime) where {P} = convert(Timestamp{P}, x)
+Base.convert(::Type{Timestamp}, x::CDFDateTime) = convert(Timestamp{Nanosecond}, x)
+Base.convert(::Type{Timestamp{P}}, x::Union{TT2000, Epoch16}) where {P} = convert(Timestamp{P}, Nanosecond(Int64(_unix_ns(x))))
+Base.convert(::Type{Timestamp{P}}, x::Epoch) where {P} = convert(Timestamp{P}, DateTime(x))
+Base.convert(::Type{DateTime}, x::CDFDateTime) = DateTime(x)
+
+# Calendar fields and printing. Epoch and Epoch16 span years 0-9999, far beyond
+# Timestamp{Nanosecond} (1677-2262), so only TT2000 gets nanosecond fields.
+_calendar(x::TT2000) = Timestamp{Nanosecond}(x)
+_calendar(x::CDFDateTime) = DateTime(x)
+
+TT2000(dt::TimeType) = convert(TT2000, dt)
+Epoch(dt::TimeType) = convert(Epoch, dt)
+Epoch16(dt::TimeType) = convert(Epoch16, dt)
 
-function Dates.DateTime(epoch::TT2000)
-    # TT2000 to Unix time with leap second correction
-    ns_from_1970 = epoch.instant.value + TT2000_OFFSET
-    leap_seconds_ns = leap_second(ns_from_1970)
-    return DateTime(1970) + Nanosecond(ns_from_1970 - leap_seconds_ns)
+function Base.convert(::Type{TT2000}, dt::TimeType)
+    ns_since_unix = _unix_ns(dt)
+    return TT2000(Int64(ns_since_unix - TT2000_OFFSET + leap_second(ns_since_unix)))
 end
 
-# Conversion from TimeType
-function Epoch16(dt::DateTime)
-    ms_since_unix = value(dt - DateTime(1970, 1, 1))
-    s_since_unix = div(ms_since_unix, 1000)
-    s_total = s_since_unix + EPOCH_OFFSET_SECONDS
-    ps_component = rem(ms_since_unix, 1000) * 1000000000  # Convert nanoseconds remainder to picoseconds
-    return Epoch16(s_total, ps_component)
+function Base.convert(::Type{Epoch16}, dt::TimeType)
+    s, ns = fldmod(_unix_ns(dt), 1_000_000_000)
+    return Epoch16(Float64(s + Int128(EPOCH_OFFSET_SECONDS)), Float64(ns * 1000))
 end
 
-function Base.convert(::Type{TT2000}, dt::TimeType)
-    ns_since_unix = (DateTime(dt) - DateTime(1970, 1, 1)).value * 1_000_000
-    leap_seconds_ns = leap_second(ns_since_unix)
-    tt2000_value = ns_since_unix - TT2000_OFFSET + leap_seconds_ns
-    return TT2000(tt2000_value)
+# Int64 path for the common case; the Int128 division above is several times slower.
+function Base.convert(::Type{Epoch16}, dt::DateTime)
+    s, ms = fldmod(value(dt) - Dates.UNIXEPOCH, 1000)
+    return Epoch16(s + EPOCH_OFFSET_SECONDS, ms * 1.0e9)
 end
 
 function Base.convert(::Type{Epoch}, dt::TimeType)
@@ -120,38 +139,70 @@ for t in (:Epoch, :Epoch16, :TT2000)
     @eval Base.convert(::Type{$t}, dt::$t) = dt
 end
 
-for f in (:year, :month, :day, :hour, :minute, :second, :millisecond)
-    @eval Dates.$f(epoch::CDFDateTime) = Dates.$f(DateTime(epoch))
+for f in (:Date, :Time, :days, :year, :month, :day, :hour, :minute, :second, :millisecond)
+    @eval Dates.$f(epoch::CDFDateTime) = Dates.$f(_calendar(epoch))
+end
+for f in (:microsecond, :nanosecond)
+    @eval Dates.$f(epoch::TT2000) = Dates.$f(_calendar(epoch))
+end
+
+# Callers guarantee `r` is in bounds.
+@inline function _parse_int(b, r)
+    x = 0
+    for i in r
+        d = @inbounds(b[i]) - 0x30
+        d <= 0x09 || throw(ArgumentError("invalid digit in date string"))
+        x = 10x + d
+    end
+    return x
+end
+
+# Returns the whole seconds as a DateTime (full year range, unlike Timestamp{Nanosecond}) and the
+# fraction in picoseconds. "yyyy-mm-ddTHH:MM:SS[.f...]" is read straight from the bytes so no
+# fractional digit is dropped; any other shape goes to the Dates parser.
+function _parse_epoch(s::AbstractString)
+    b = codeunits(s)
+    n = length(b)
+    iso = @inbounds (n == 19 || 21 <= n <= 32 && b[20] == UInt8('.')) && b[5] == b[8] == UInt8('-') &&
+        b[11] == UInt8('T') && b[14] == b[17] == UInt8(':')
+    iso || return DateTime(s), 0
+    dt = DateTime(
+        _parse_int(b, 1:4), _parse_int(b, 6:7), _parse_int(b, 9:10),
+        _parse_int(b, 12:13), _parse_int(b, 15:16), _parse_int(b, 18:19)
+    )
+    return dt, _parse_int(b, 21:n) * 10^(32 - n)
 end
 
 Epoch(s::AbstractString) = Epoch(DateTime(s))
-Epoch16(s::AbstractString) = Epoch16(DateTime(s))
-TT2000(s::AbstractString) = TT2000(DateTime(s))
+function Epoch16(s::AbstractString)
+    dt, ps = _parse_epoch(s)
+    epoch = Epoch16(dt)
+    return Epoch16(epoch.seconds, epoch.picoseconds + ps)
+end
+function TT2000(s::AbstractString)
+    dt, ps = _parse_epoch(s)
+    ns, r = divrem(ps, 1000)
+    iszero(r) || throw(InexactError(:TT2000, TT2000, s))
+    return TT2000(dt) + Nanosecond(ns)
+end
 
 Dates.value(epoch::Epoch) = epoch.instant
 Dates.value(epoch::Epoch16) = ComplexF64(epoch.seconds, epoch.picoseconds)
 Dates.value(epoch::TT2000) = epoch.instant.value
 
 function Base.floor(x::T, p::Union{DatePeriod, TimePeriod}) where {T <: CDFDateTime}
-    return convert(T, floor(convert(DateTime, x), p))
+    return convert(T, floor(_calendar(x), p))
 end
 
-function Base.show(io::IO, epoch::CDFDateTime)
-    fillval = fillvalue(epoch)
-    return if fillval == Dates.value(epoch)
-        print(io, "FILLVAL")
-    else
-        print(io, DateTime(epoch))
-    end
-end
-function Base.show(io::IO, epoch::Epoch16)
-    return print(io, DateTime(epoch))
-end
+Base.show(io::IO, epoch::CDFDateTime) = _isfill(epoch) ? print(io, "FILLVAL") : print(io, _calendar(epoch))
 
-Base.promote_rule(::Type{<:CDFDateTime}, ::Type{Dates.DateTime}) = Dates.DateTime
-Base.promote_rule(::Type{T}, ::Type{Dates.Date}) where {T <: CDFDateTime} = T
-# Comment out because of invalidation
-Base.convert(::Type{Dates.DateTime}, x::CDFDateTime) = Dates.DateTime(x)
+# Promote towards the CDF type: UTC -> TT2000 is injective, whereas TT2000 -> UTC folds each
+# inserted leap second onto the following second, so comparing in UTC would equate distinct instants.
+Base.promote_rule(::Type{T}, ::Type{<:Union{Date, DateTime}}) where {T <: CDFDateTime} = T
+Base.promote_rule(::Type{T}, ::Type{<:Timestamp}) where {T <: Union{TT2000, Epoch16}} = T
+Base.promote_rule(::Type{Epoch}, ::Type{T}) where {T <: Timestamp} = promote_type(DateTime, T)
+Base.promote_rule(::Type{TT2000}, ::Type{<:Union{Epoch, Epoch16}}) = TT2000
+Base.promote_rule(::Type{Epoch16}, ::Type{Epoch}) = Epoch16
 Base.bswap(x::Epoch) = Epoch(Base.bswap(x.instant))
 Base.bswap(x::Epoch16) = Epoch16(Base.bswap(x.seconds), Base.bswap(x.picoseconds))
 Base.bswap(x::TT2000) = TT2000(Base.bswap(x.instant.value))

src/leap_second.jl

@@ -5,7 +5,6 @@
 # TT2000 offset constant from C++ implementation
 const TT2000_OFFSET = Int64(946727967816000000)  # nanoseconds from 1970 to J2000 with corrections
 const UNIX_EPOCH = DateTime(1970, 1, 1)
-const PRE1972_CUTOFF = DateTime(1972, 1, 1)
 const MJD_BASE = 2_400_000.5
 const NS_IN_SECOND = Int64(1_000_000_000)
 
@@ -75,13 +74,21 @@ function leap_seconds_pre1972(date)
     return entry[2] + (mjd - entry[3]) * entry[4]
 end
 
-function leap_second(ns_from_1970::Int64)
-    dt = UNIX_EPOCH + Nanosecond(ns_from_1970)
-    if dt >= PRE1972_CUTOFF
+# Branch and floor on the nanosecond count: `DateTime + Nanosecond` rounds to the millisecond,
+# which would push instants just below a day boundary into the next day's offset.
+function leap_second(ns_from_1970::Integer)
+    if ns_from_1970 >= LEAP_SECONDS_TT2000[1][1]
         idx = findlast(x -> x[1] <= ns_from_1970, LEAP_SECONDS_TT2000)
         return LEAP_SECONDS_TT2000[idx][2]
     else
-        seconds = leap_seconds_pre1972(dt)
+        seconds = leap_seconds_pre1972(UNIX_EPOCH + Millisecond(Int64(fld(ns_from_1970, 1_000_000))))
         return floor(Int64, seconds * NS_IN_SECOND)
     end
 end
+
+# Inverse of `utc + leap_second(utc)`. The second lookup corrects for the TT scale being ahead by the
+# offset; a TT instant inside an inserted leap second maps onto the following UTC second (as cdflib does).
+function utc_from_tt(tt_ns_from_1970::Integer)
+    utc = tt_ns_from_1970 - leap_second(tt_ns_from_1970)
+    return tt_ns_from_1970 - leap_second(utc)
+end

src/precompile.jl

@@ -19,7 +19,7 @@ PrecompileTools.@setup_workload begin
                 Array(var)
                 var[1]
                 var[ntuple(_ -> Colon(), ndims(var) - 1)..., 1:1]
-                eltype(var) <: CDFDateTime && DateTime.(var[1:1])
+                eltype(var) <: CDFDateTime && (DateTime.(var[1:1]); Timestamp.(var[1:1]); string(var[1]))
             end
         end
     end

test/Project.toml

@@ -3,6 +3,7 @@ Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 CommonDataFormat = "a9737db6-c05c-4e48-868b-6bc41491d9d9"
 CommonDataModel = "1fbeeb36-5f17-413c-809b-666fb144f157"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
+Durations = "41ac09f5-0a95-4810-aa18-c88bf48ee130"
 Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
 JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b"
 JuliaC = "acedd4c2-ced6-4a15-accc-2607eb759ba2"