progress

Author: akoumjian

src/adam_core/time/tests/test_time.py

@@ -7,7 +7,7 @@
 import pytest
 import quivr as qv
 
-from ..time import Timestamp
+from ..time import NANOS_PER_DAY, Timestamp
 
 
 class Wrapper(qv.Table):
@@ -699,14 +699,19 @@ def test_Timestamp_rescale(scale1, scale2):
     Test rescaling by using round trip calculations
     """
     times = Timestamp.from_kwargs(
-        days=[-57023, 0, 51544, 103088, 164178, 68000, 68000, 68010, 68020],  # Spans from ~1702 to ~2308
-        nanos=[100_000_000, 200_000_000, 300_000_000, 400_000_000, 500_000_000, 1, 2, 3, 4],
+        days=[-57032, -36525, -2, -1, 0, 51544, 103088, 164178, 68000, 68000, 68010, 68020],  # Spans from ~1702 to ~2308
+        nanos=[50_000, 0, 123, 100_000_000, 200_000_000, 300_000_000, 400_000_000, 500_000_000, 1, 2, 3, 4],
         scale=scale1,
     )
     rescaled = times.rescale(scale2)
     round_tripped = rescaled.rescale(scale1)
     assert rescaled.scale == scale2
     assert round_tripped.scale == scale1
-    print(times.difference(round_tripped))
-    assert pc.all(times.equals(round_tripped, precision="us")).as_py()
+    print(scale1, scale2)
+    mjds = list(zip(times.mjd().to_pylist(), rescaled.mjd().to_pylist(), round_tripped.mjd().to_pylist()))
+    days_diff, nanos_diff = times.difference(round_tripped)
+    diff = list(zip(days_diff.to_pylist(), nanos_diff.to_pylist()))
+    for i in range(len(mjds)):
+        print(mjds[i], diff[i])
+    assert pc.all(times.equals(round_tripped, precision="ns")).as_py()
 

src/adam_core/time/time.py

@@ -21,6 +21,10 @@
 MJD_EPOCH_IN_TAI = hifitime.Epoch("1858-11-17T00:00:00 TAI")
 MJD_EPOCH_IN_UTC = hifitime.Epoch("1858-11-17T00:00:00 UTC")
 
+DAYS_PER_CENTURY = 36_525
+NANOS_PER_DAY = 86_400_000_000_000
+NANOS_PER_CENTURY = DAYS_PER_CENTURY * NANOS_PER_DAY
+
 class Timestamp(qv.Table):
     # Scale, the rate at which time passes:
     scale = qv.StringAttribute(default="tai")
@@ -503,30 +507,43 @@ def rescale(self, new_scale: str) -> Timestamp:
 
         # Get only the values needed based on the input scale
         if self.scale == "tdb":
-            day_values = self.jd().to_numpy()
-            init_epoch_fn = lambda val: hifitime.Epoch.init_from_jde_tdb(float(val))
+            def init_epoch_fn(days, nanos):
+                total_nanos = int(days) * NANOS_PER_DAY + int(nanos)
+                duration_since_mjd_epoch = hifitime.Duration.from_total_nanoseconds(total_nanos)
+                epoch = MJD_EPOCH_IN_TDB + duration_since_mjd_epoch
+                return epoch
         else:  # tai or utc
-            day_values = self.mjd().to_numpy()
             if self.scale == "tai":
-                init_epoch_fn = lambda val: hifitime.Epoch.init_from_mjd_tai(float(val))
+                def init_epoch_fn(days, nanos):
+                    total_nanos = int(days) * NANOS_PER_DAY + int(nanos)
+                    duration_since_mjd_epoch = hifitime.Duration.from_total_nanoseconds(total_nanos)
+                    epoch = MJD_EPOCH_IN_TAI + duration_since_mjd_epoch
+                    return epoch
             else:  # utc
-                init_epoch_fn = lambda val: hifitime.Epoch.init_from_mjd_utc(float(val))
+                def init_epoch_fn(days, nanos):
+                    total_nanos = int(days) * NANOS_PER_DAY + int(nanos)
+                    duration_since_mjd_epoch = hifitime.Duration.from_total_nanoseconds(total_nanos)
+                    epoch = MJD_EPOCH_IN_UTC + duration_since_mjd_epoch
+                    return epoch
 
         # Determine the result calculation based on the new scale
         if new_scale == "tai":
-            calculate_result = lambda epoch: epoch.to_mjd_tai_days()
+            calculate_result = _calculate_tai_result
         elif new_scale == "utc":
-            calculate_result = lambda epoch: epoch.to_mjd_utc_days()
+            calculate_result = _calculate_utc_result
         elif new_scale == "tdb":
-            calculate_result = lambda epoch: epoch.timedelta(MJD_EPOCH_IN_TDB).to_unit(hifitime.Unit.Day)
+            calculate_result = _calculate_tdb_result
         else:
             raise ValueError(f"Unknown scale: {new_scale}")
 
         # Process each timestamp individually
-        result_values = [calculate_result(init_epoch_fn(val)) for val in day_values]
-
+        inputs = zip(self.days.to_numpy(), self.nanos.to_numpy())
+        init_epochs = [init_epoch_fn(*val) for val in inputs]
+        result_values = [calculate_result(init_epoch, self.scale) for init_epoch in init_epochs]
+        result_days = [val[0] for val in result_values]
+        result_nanos = [val[1] for val in result_values]
         # Convert result list back to a Timestamp
-        return self.from_mjd(pa.array(result_values), scale=new_scale)
+        return self.from_kwargs(days=result_days, nanos=result_nanos, scale=new_scale)
 
     def link(
         self, other: Timestamp, precision: str = "ns"
@@ -554,8 +571,6 @@ def link(
         rounded = self.rounded(precision)
         other_rounded = other.rounded(precision)
 
-        import pdb; pdb.set_trace()
-
         left_keys = {"days": rounded.days, "nanos": rounded.nanos}
         right_keys = {"days": other_rounded.days, "nanos": other_rounded.nanos}
         return qv.MultiKeyLinkage(self, other, left_keys, right_keys)
@@ -582,3 +597,59 @@ def _duration_arrays_within_tolerance(
         pc.greater_equal(pc.abs(delta_nanos), 86400 * 1e9 - max_nanos_deviation),
     )
     return pc.or_(cond1, cond2)
+
+def _calculate_tai_result(epoch, from_scale: str):
+    epoch_new_time_scale = epoch.to_time_scale(hifitime.TimeScale.TAI)
+    tai_duration_since_mjd_epoch = epoch_new_time_scale.timedelta(MJD_EPOCH_IN_TAI)
+    
+    total_nanos = tai_duration_since_mjd_epoch.total_nanoseconds()
+    
+    # Integer division rounds towards negative infinity for negative numbers
+    days = total_nanos // NANOS_PER_DAY
+    
+    if days < -36525:
+        total_nanos += NANOS_PER_CENTURY
+        days = total_nanos // NANOS_PER_DAY
+
+    # Modulo will give us a positive remainder
+    nanos = total_nanos % NANOS_PER_DAY
+    
+    return int(days), int(nanos)
+
+def _calculate_utc_result(epoch, from_scale: str):
+    epoch_new_time_scale = epoch.to_time_scale(hifitime.TimeScale.UTC)
+    utc_duration_since_mjd_epoch = epoch_new_time_scale.timedelta(MJD_EPOCH_IN_UTC)
+    
+    total_nanos = utc_duration_since_mjd_epoch.total_nanoseconds()
+    
+    # Integer division rounds towards negative infinity for negative numbers
+    days = total_nanos // NANOS_PER_DAY
+
+    # Special case for the Julian century boundary
+    if days < -36525:
+        total_nanos += NANOS_PER_CENTURY
+        days = total_nanos // NANOS_PER_DAY
+
+    # Modulo will give us a positive remainder
+    nanos = total_nanos % NANOS_PER_DAY
+    
+    return int(days), int(nanos)
+
+def _calculate_tdb_result(epoch, from_scale: str):
+    epoch_new_time_scale = epoch.to_time_scale(hifitime.TimeScale.TDB)
+    tdb_duration_since_mjd_epoch = epoch_new_time_scale.timedelta(MJD_EPOCH_IN_TDB)
+    
+    total_nanos = tdb_duration_since_mjd_epoch.total_nanoseconds()
+    
+    # Integer division rounds towards negative infinity for negative numbers
+    days = total_nanos // NANOS_PER_DAY
+    
+    # Special case for the Julian century boundary
+    if days < -36525:
+        total_nanos += NANOS_PER_CENTURY
+        days = total_nanos // NANOS_PER_DAY
+    
+    # Modulo will give us a positive remainder
+    nanos = total_nanos % NANOS_PER_DAY
+    
+    return int(days), int(nanos)