From 7eeadbce1ecbd60537846c7b65db6739b4d53147 Mon Sep 17 00:00:00 2001 From: Matteo Bachetti Date: Tue, 24 Mar 2026 09:04:04 +0100 Subject: [PATCH 01/34] Update numpy dependency to allow >=2.0 ; Update doctests to work with new numpy --- nuclockutils/diagnostics/fold_to_ephemeris.py | 2 +- nuclockutils/nustarclock.py | 14 ++--- nuclockutils/utils.py | 60 +++++++++---------- pyproject.toml | 2 +- 4 files changed, 39 insertions(+), 39 deletions(-) diff --git a/nuclockutils/diagnostics/fold_to_ephemeris.py b/nuclockutils/diagnostics/fold_to_ephemeris.py index 6eae679..54a938f 100644 --- a/nuclockutils/diagnostics/fold_to_ephemeris.py +++ b/nuclockutils/diagnostics/fold_to_ephemeris.py @@ -132,7 +132,7 @@ def get_exposure_per_bin(event_times, event_priors, parfile, # >>> prof = get_exposure_per_bin(event_times, event_priors, 1/10, nbin=10) # >>> prof[0] # 0 - # >>> np.allclose(prof[1:], 1) + # >>> bool(np.allclose(prof[1:], 1)) # True """ log.info("Calculating exposure") diff --git a/nuclockutils/nustarclock.py b/nuclockutils/nustarclock.py index 7d9ed0a..0f95e58 100644 --- a/nuclockutils/nustarclock.py +++ b/nuclockutils/nustarclock.py @@ -54,7 +54,7 @@ def flag_bad_points(all_data, db_file='BAD_POINTS_DB.dat'): >>> all_data = Table({'met': [0, 1, 2, 3, 4]}) >>> all_data = flag_bad_points(all_data, db_file='dummy_bad_points.dat') INFO: ... - >>> np.all(all_data['flag'] == [False, False, False, True, False]) + >>> bool(np.all(all_data['flag'] == [False, False, False, True, False])) True """ if not os.path.exists(db_file): @@ -354,9 +354,9 @@ def residual_roll_std(residuals, window=30): >>> residuals = np.zeros(5000) >>> residuals[:4000] = np.random.normal(0, 1, 4000) >>> roll_std = residual_roll_std(residuals, window=500) - >>> np.allclose(roll_std[:3500], 1., rtol=0.2) + >>> bool(np.allclose(roll_std[:3500], 1., rtol=0.2)) True - >>> np.all(roll_std[4500:] == 0.) + >>> bool(np.all(roll_std[4500:] == 0.)) True """ r_std = rolling_std(residuals, window) @@ -532,10 +532,10 @@ def no_jump_gtis(start_time, stop_time, clock_jump_times=None): Examples -------- >>> gtis = no_jump_gtis(0, 3, [1, 1.1]) - >>> np.allclose(gtis, [[0, 1], [1, 1.1], [1.1, 3]]) + >>> bool(np.allclose(gtis, [[0, 1], [1, 1.1], [1.1, 3]])) True >>> gtis = no_jump_gtis(0, 3) - >>> np.allclose(gtis, [[0, 3]]) + >>> bool(np.allclose(gtis, [[0, 3]])) True """ if clock_jump_times is None: @@ -558,11 +558,11 @@ def temperature_gtis(temperature_table, max_distance=600): -------- >>> temperature_table = Table({'met': [0, 1, 2, 10, 11, 12]}) >>> gti = temperature_gtis(temperature_table, 5) - >>> np.allclose(gti, [[0, 2], [10, 12]]) + >>> bool(np.allclose(gti, [[0, 2], [10, 12]])) True >>> temperature_table = Table({'met': [-10, 0, 1, 2, 10, 11, 12, 20]}) >>> gti = temperature_gtis(temperature_table, 5) - >>> np.allclose(gti, [[0, 2], [10, 12]]) + >>> bool(np.allclose(gti, [[0, 2], [10, 12]])) True """ temp_condition = np.concatenate( diff --git a/nuclockutils/utils.py b/nuclockutils/utils.py index c24e417..e3ab689 100644 --- a/nuclockutils/utils.py +++ b/nuclockutils/utils.py @@ -80,13 +80,13 @@ def splitext_improved(path): """ Examples -------- - >>> np.all(splitext_improved("a.tar.gz") == ('a', '.tar.gz')) + >>> bool(np.all(splitext_improved("a.tar.gz") == ('a', '.tar.gz'))) True - >>> np.all(splitext_improved("a.tar") == ('a', '.tar')) + >>> bool(np.all(splitext_improved("a.tar") == ('a', '.tar'))) True - >>> np.all(splitext_improved("a.f/a.tar") == ('a.f/a', '.tar')) + >>> bool(np.all(splitext_improved("a.f/a.tar") == ('a.f/a', '.tar'))) True - >>> np.all(splitext_improved("a.a.a.f/a.tar.gz") == ('a.a.a.f/a', '.tar.gz')) + >>> bool(np.all(splitext_improved("a.a.a.f/a.tar.gz") == ('a.a.a.f/a', '.tar.gz'))) True """ import os @@ -231,10 +231,10 @@ def rolling_window(a, window): -------- >>> a = np.arange(5) >>> rw = rolling_window(a, 2) - >>> np.allclose(rw, [[0, 1], [1,2], [2, 3], [3, 4]]) + >>> bool(np.allclose(rw, [[0, 1], [1,2], [2, 3], [3, 4]])) True >>> rw = rolling_window(a, 3) - >>> np.allclose(rw, [[0, 1, 2], [1, 2, 3], [2, 3, 4]]) + >>> bool(np.allclose(rw, [[0, 1, 2], [1, 2, 3], [2, 3, 4]])) True """ shape = a.shape[:-1] + (a.shape[-1] - window + 1, window) @@ -248,13 +248,13 @@ def rolling_stat(stat_fun, a, window, pad='center', **kwargs): -------- >>> a = np.arange(6) >>> r_sum = rolling_stat(np.sum, a, 3, pad='center', axis=-1) - >>> np.allclose(r_sum, [3., 3., 6., 9., 12., 12.]) + >>> bool(np.allclose(r_sum, [3., 3., 6., 9., 12., 12.])) True >>> r_sum = rolling_stat(np.sum, a, 3, pad='left', axis=-1) - >>> np.allclose(r_sum, [3., 3., 3., 6., 9., 12.]) + >>> bool(np.allclose(r_sum, [3., 3., 3., 6., 9., 12.])) True >>> r_sum = rolling_stat(np.sum, a, 3, pad='right', axis=-1) - >>> np.allclose(r_sum, [3., 6., 9., 12., 12., 12.]) + >>> bool(np.allclose(r_sum, [3., 6., 9., 12., 12., 12.])) True >>> r_sum = rolling_stat(np.sum, a, 3, pad='incredible', axis=-1) Traceback (most recent call last): @@ -294,7 +294,7 @@ def rolling_std(a, window, pad='center'): Examples >>> a = [0, 1, 1, 3] - >>> np.allclose(rolling_std(a, 2), [0.5, 0, 1, 1]) + >>> bool(np.allclose(rolling_std(a, 2), [0.5, 0, 1, 1])) True """ return rolling_stat(np.std, a, window, pad, axis=-1) @@ -309,7 +309,7 @@ def spline_through_data(x, y, k=2, grace_intv=1000., smoothing_factor=0.001, >>> x = np.arange(1000) >>> y = np.random.normal(x * 0.1, 0.01) >>> fun = spline_through_data(x, y, grace_intv=10.) - >>> np.std(y - fun(x)) < 0.01 + >>> bool(np.std(y - fun(x)) < 0.01) True """ lo_lim, hi_lim = x[0], x[-1] @@ -343,17 +343,17 @@ def aggregate(table, max_number=1000): >>> newt = aggregate(table) >>> len(newt) 2 - >>> np.all(newt['a'] == table['a']) + >>> bool(np.all(newt['a'] == table['a'])) True - >>> np.all(newt['b'] == table['b']) + >>> bool(np.all(newt['b'] == table['b'])) True >>> newt = aggregate(table, max_number=1) >>> len(newt) 1 - >>> np.all(newt['a'] == 1.5) + >>> bool(np.all(newt['a'] == 1.5)) True >>> newt = aggregate(table.to_pandas(), max_number=1) - >>> np.all(newt['b'] == 5.5) + >>> bool(np.all(newt['b'] == 5.5)) True """ N = len(table) @@ -377,9 +377,9 @@ def aggregate_all_tables(table_list, max_number=1000): >>> newt = aggregate_all_tables([table])[0] >>> len(newt) 2 - >>> np.all(newt['a'] == table['a']) + >>> bool(np.all(newt['a'] == table['a'])) True - >>> np.all(newt['b'] == table['b']) + >>> bool(np.all(newt['b'] == table['b'])) True """ return [aggregate(table) for table in table_list] @@ -408,12 +408,12 @@ def cross_two_gtis(gti0, gti1): >>> gti1 = np.array([[1, 2]]) >>> gti2 = np.array([[1, 2]]) >>> newgti = cross_two_gtis(gti1, gti2) - >>> np.all(newgti == [[1, 2]]) + >>> bool(np.all(newgti == [[1, 2]])) True >>> gti1 = np.array([[1, 4]]) >>> gti2 = np.array([[1, 2], [2, 4]]) >>> newgti = cross_two_gtis(gti1, gti2) - >>> np.all(newgti == [[1, 2], [2, 4]]) + >>> bool(np.all(newgti == [[1, 2], [2, 4]])) True """ import copy @@ -455,7 +455,7 @@ def robust_poly_fit(x, y, order=3, p0=None): >>> x = np.arange(10) >>> y = x**2 >>> fun = robust_poly_fit(x, y, order=2, p0=np.zeros(3)) - >>> np.allclose(fun(x), y) + >>> bool(np.allclose(fun(x), y)) True """ from scipy.optimize import least_squares @@ -487,10 +487,10 @@ def measure_overall_trend(x, y, ref_size=200): >>> val[0] is None True >>> val = measure_overall_trend(np.asarray([0, 1]), np.asarray([1, 1])) - >>> np.allclose(val, [0, 0, 1]) + >>> bool(np.allclose(val, [0, 0, 1])) True >>> val = measure_overall_trend(np.arange(1000), np.ones(1000)) - >>> np.allclose(val, [0, 0, 1]) + >>> bool(np.allclose(val, [0, 0, 1])) True """ x0 = x[0] @@ -525,11 +525,11 @@ def get_rough_trend_fun(met, residuals): True >>> x, y = np.asarray([0, 1]), np.asarray([1, 1]) >>> fun = get_rough_trend_fun(x, y) - >>> np.allclose(fun(x), y) + >>> bool(np.allclose(fun(x), y)) True >>> x, y = np.arange(1000), np.ones(1000) >>> fun = get_rough_trend_fun(x, y) - >>> np.allclose(fun(x), y) + >>> bool(np.allclose(fun(x), y)) True """ x0, m, q = measure_overall_trend(met, residuals) @@ -550,7 +550,7 @@ def merge_and_sort_arrays(array1, array2): -------- >>> arr1 = np.array([0, 3., 1]) >>> arr2 = np.array([2, 0]) - >>> np.allclose(merge_and_sort_arrays(arr1, arr2), [0, 1, 2, 3]) + >>> bool(np.allclose(merge_and_sort_arrays(arr1, arr2), [0, 1, 2, 3])) True """ arr = np.concatenate((array1, array2)) @@ -567,7 +567,7 @@ def eliminate_array_from_array(array1, array2): -------- >>> arr1 = np.array([0, 3., 1]) >>> arr2 = np.array([4, 0]) - >>> np.allclose(eliminate_array_from_array(arr1, arr2), [1, 3]) + >>> bool(np.allclose(eliminate_array_from_array(arr1, arr2), [1, 3])) True """ array1 = np.asarray(copy.deepcopy(array1)) @@ -685,11 +685,11 @@ def high_precision_keyword_read(hdr, keyword): Examples -------- >>> hdr = dict(keywordS=1.25) - >>> high_precision_keyword_read(hdr, 'keywordS') - 1.25 + >>> bool(np.isclose(high_precision_keyword_read(hdr, 'keywordS'), 1.25)) + True >>> hdr = dict(keywordI=1, keywordF=0.25) - >>> high_precision_keyword_read(hdr, 'keywordS') - 1.25 + >>> bool(np.isclose(high_precision_keyword_read(hdr, 'keywordS'), 1.25)) + True >>> high_precision_keyword_read(hdr, 'bubabuab') is None True """ diff --git a/pyproject.toml b/pyproject.toml index 47d1e49..1be90c1 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -16,7 +16,7 @@ authors = [ { name = "Matteo Bachetti", email = "matteo.bachetti@inaf.it" }, ] dependencies = [ - "numpy<2.0", + "numpy", "scipy", "astropy", "matplotlib", From 1a67886f7a1ad7aa0b698025cb7ae5320ef5cdcd Mon Sep 17 00:00:00 2001 From: Matteo Bachetti Date: Tue, 24 Mar 2026 09:09:32 +0100 Subject: [PATCH 02/34] Update numpy dep also in requirements.txt --- requirements.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/requirements.txt b/requirements.txt index c70fbd9..4dc4277 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,4 +1,4 @@ -numpy<2.0 +numpy scipy astropy matplotlib From dfdc6e9cbe0194eb8bc39e55afde38d4e38fb0bd Mon Sep 17 00:00:00 2001 From: Matteo Bachetti Date: Tue, 24 Mar 2026 09:09:46 +0100 Subject: [PATCH 03/34] Expand tested versions --- .github/workflows/ci-test.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/ci-test.yml b/.github/workflows/ci-test.yml index 631504f..5349368 100644 --- a/.github/workflows/ci-test.yml +++ b/.github/workflows/ci-test.yml @@ -16,7 +16,7 @@ jobs: strategy: fail-fast: false matrix: - python-version: ["3.9", "3.10", "3.11"] + python-version: ["3.9", "3.11", "3.13"] steps: - uses: actions/checkout@v4 From 9fc685953eab82bf22ed7b791bf32ba204ff88f5 Mon Sep 17 00:00:00 2001 From: Matteo Bachetti Date: Tue, 24 Mar 2026 09:19:44 +0100 Subject: [PATCH 04/34] Be more restrictive with versions for old pythons --- pyproject.toml | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index 1be90c1..d412409 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -16,7 +16,8 @@ authors = [ { name = "Matteo Bachetti", email = "matteo.bachetti@inaf.it" }, ] dependencies = [ - "numpy", + "numpy < 2.0 ; python_version < '3.11'", + "numpy >= 2.0 ; python_version >= '3.11'", "scipy", "astropy", "matplotlib", @@ -32,10 +33,12 @@ dependencies = [ "h5py", "astroquery", "plotly", - "dash", - "stingray>=0.3", + "dash < 4.0.0 ; python_version < '3.11'", + "stingray (>=0.3, < 2.3.0) ; python_version < '3.11'", + "stingray>=0.3 ; python_version >= '3.11'", "pyyaml", - "hendrics", + "hendrics < 8.4 ; python_version < '3.11'", + "hendrics >= 8.4 ; python_version >= '3.11'", ] dynamic = ["version"] From 56b631436717721769f62e99845e981d6a9567bf Mon Sep 17 00:00:00 2001 From: Matteo Bachetti Date: Tue, 24 Mar 2026 09:24:01 +0100 Subject: [PATCH 05/34] Ditch need for requirements.txt --- docs/index.rst | 1 - requirements.txt | 21 --------------------- 2 files changed, 22 deletions(-) delete mode 100644 requirements.txt diff --git a/docs/index.rst b/docs/index.rst index 380a53e..ea466bc 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -11,7 +11,6 @@ Installation $ git clone https://github.com/matteobachetti/nustar-clock-utils $ cd nustar-clock-utils - $ pip install -r requirements.txt $ pip install . Usage diff --git a/requirements.txt b/requirements.txt deleted file mode 100644 index 4dc4277..0000000 --- a/requirements.txt +++ /dev/null @@ -1,21 +0,0 @@ -numpy -scipy -astropy -matplotlib -pandas -scikit-learn -pint-pulsar -#git+https://github.com/nanograv/pint.git@e5f6b260422898f63c88d636712fcc1cdd76cf20 -tqdm -jplephem -datashader -bokeh -holoviews>=1.18 -statsmodels -h5py -astroquery -plotly -dash -stingray>=0.3 -pyyaml -hendrics From 8924c1fbb298c54a938e75d7551440ba25fb9490 Mon Sep 17 00:00:00 2001 From: Matteo Bachetti Date: Tue, 24 Mar 2026 09:27:44 +0100 Subject: [PATCH 06/34] Fix installation script --- .github/workflows/ci-test.yml | 1 + README.rst | 1 - 2 files changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/ci-test.yml b/.github/workflows/ci-test.yml index 5349368..9f225b7 100644 --- a/.github/workflows/ci-test.yml +++ b/.github/workflows/ci-test.yml @@ -29,6 +29,7 @@ jobs: python -m pip install --upgrade pip python -m pip install pytest pytest-astropy if [ -f requirements.txt ]; then pip install -r requirements.txt; fi + pip install . - name: Test with pytest run: | pytest --remote-data -svv nuclockutils diff --git a/README.rst b/README.rst index 5008c8d..489cbe6 100644 --- a/README.rst +++ b/README.rst @@ -19,7 +19,6 @@ Installation $ git clone https://github.com/matteobachetti/nustar-clock-utils $ cd nustar-clock-utils - $ pip install -r requirements.txt $ pip install . Usage From a28b3324a18e78374bb96e6e431b650a248e55a0 Mon Sep 17 00:00:00 2001 From: Matteo Bachetti Date: Tue, 24 Mar 2026 09:30:13 +0100 Subject: [PATCH 07/34] Do not forget dash in new pythons --- pyproject.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/pyproject.toml b/pyproject.toml index d412409..e90485f 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -34,6 +34,7 @@ dependencies = [ "astroquery", "plotly", "dash < 4.0.0 ; python_version < '3.11'", + "dash >= 4.0.0 ; python_version >= '3.11'", "stingray (>=0.3, < 2.3.0) ; python_version < '3.11'", "stingray>=0.3 ; python_version >= '3.11'", "pyyaml", From 1fa8b80c4503554ed067cd154a091dc9e7a82ac4 Mon Sep 17 00:00:00 2001 From: Matteo Bachetti Date: Tue, 24 Mar 2026 11:50:45 +0100 Subject: [PATCH 08/34] Fix spline_through_data --- nuclockutils/utils.py | 38 ++++++++++++++++++++++++-------------- 1 file changed, 24 insertions(+), 14 deletions(-) diff --git a/nuclockutils/utils.py b/nuclockutils/utils.py index e3ab689..f918dc5 100644 --- a/nuclockutils/utils.py +++ b/nuclockutils/utils.py @@ -6,7 +6,7 @@ from astropy.table import Table from astroquery.heasarc import Heasarc from astropy.time import Time -from scipy.interpolate import LSQUnivariateSpline +from scipy.interpolate import make_lsq_spline from numba import vectorize try: @@ -316,22 +316,32 @@ def spline_through_data(x, y, k=2, grace_intv=1000., smoothing_factor=0.001, control_points = \ np.linspace(lo_lim + 2 * grace_intv, hi_lim - 2 * grace_intv, - int(x.size // downsample)) + int(max(x.size // downsample, k + 1))) + + control_points = np.r_[(x[0],)*(k+1), control_points, (x[-1],)*(k+1)] if fixed_control_points is not None and len(fixed_control_points) > 0: - log.debug(f'Adding fixed control points: {fixed_control_points}') - control_points = np.sort( - np.concatenate((control_points, fixed_control_points))) + good_fcp = fixed_control_points[(fixed_control_points > lo_lim) & (fixed_control_points < hi_lim)] + if good_fcp.size > 0: + log.debug(f'Adding fixed control points: {good_fcp}') + control_points = np.sort( + np.concatenate((control_points, good_fcp))) + else: + log.debug(f'No fixed control points within the data range, ignoring.') + + flag_ctrl_pts = np.ones(control_points.size, dtype=bool) + t_idxs = np.searchsorted(x, control_points) + for i, (id0, id1) in enumerate(zip(t_idxs[:-1], t_idxs[1:])): + if not id1 > id0 and np.any((x[id0:id1] > control_points[i]) & (x[id0:id1] < control_points[i + 1])): + # Schoenberg-Whitney condition is not satisfied, we need to remove this control point + log.info(f"No data between {control_points[i]} and {control_points[i + 1]}, removing control point.") + flag_ctrl_pts[i] = False + + control_points = control_points[flag_ctrl_pts] + + detrend_fun = make_lsq_spline( + x, y, t=control_points, k=k) - try: - detrend_fun = LSQUnivariateSpline( - x, y, t=control_points, k=k, - bbox=[lo_lim, hi_lim]) - except ValueError as e: - log.error(f"Error in LSQUnivariateSpline: {e};\nDecreasing the number of control points" - r" by 10% and trying again.") - return spline_through_data(x, y, k=k, grace_intv=grace_intv, - downsample=downsample * 1.5, fixed_control_points=fixed_control_points) return detrend_fun From e7b9b6587e62e2aed3641cc3703c01ba8cbdfe61 Mon Sep 17 00:00:00 2001 From: Matteo Bachetti Date: Tue, 24 Mar 2026 14:25:02 +0100 Subject: [PATCH 09/34] Update interpolation function and fix bugs --- NEWS.rst | 4 ++++ nuclockutils/clean_clock/app.py | 3 ++- nuclockutils/nustarclock.py | 2 +- nuclockutils/utils.py | 5 +++-- 4 files changed, 10 insertions(+), 4 deletions(-) create mode 100644 NEWS.rst diff --git a/NEWS.rst b/NEWS.rst new file mode 100644 index 0000000..30461bc --- /dev/null +++ b/NEWS.rst @@ -0,0 +1,4 @@ +nuclockutils 0.4.dev5+g48a3009ec.d20260324 (2026-03-24) +======================================================= + +No significant changes. diff --git a/nuclockutils/clean_clock/app.py b/nuclockutils/clean_clock/app.py index 6434b09..e177d9d 100644 --- a/nuclockutils/clean_clock/app.py +++ b/nuclockutils/clean_clock/app.py @@ -64,7 +64,7 @@ def recalc(outfile='save_all.pickle'): met_start = clock_offset_table['met'][0] met_stop = clock_offset_table['met'][-1] + 30 clock_jump_times = \ - np.array([77469000, 78708320, 79657575, 81043985, 82055671, 293346772, + np.array([77469000, 78708320, 79657575, 81043985, 82055671, 293346772, 305080000, 392200784, 394825882, 395304135, 407914525, 408299422]) clock_jump_times += 30 # Sum 30 seconds to avoid to exclude these points # from previous interval @@ -100,6 +100,7 @@ def recalc(outfile='save_all.pickle'): np.array( clock_offset_table['offset'] - table_new['temp_corr'][tempcorr_idx] ) + clock_residuals_detrend = np.array( clock_offset_table['offset'] - table_new['temp_corr_detrend'][tempcorr_idx]) diff --git a/nuclockutils/nustarclock.py b/nuclockutils/nustarclock.py index 0f95e58..faeaacf 100644 --- a/nuclockutils/nustarclock.py +++ b/nuclockutils/nustarclock.py @@ -116,7 +116,7 @@ def calculate_stats(all_data): print("----------------------------- Stats -----------------------------------") print() print(f"Overall MAD: {scatter * 1e6:.0f} us") - print(f"Minimum scatter: ±{np.min(r_std) * 1e6:.0f} us") + print(f"Minimum scatter: ±{np.nanmin(r_std) * 1e6:.0f} us") print() print("-----------------------------------------------------------------------") diff --git a/nuclockutils/utils.py b/nuclockutils/utils.py index f918dc5..2518de8 100644 --- a/nuclockutils/utils.py +++ b/nuclockutils/utils.py @@ -318,7 +318,6 @@ def spline_through_data(x, y, k=2, grace_intv=1000., smoothing_factor=0.001, np.linspace(lo_lim + 2 * grace_intv, hi_lim - 2 * grace_intv, int(max(x.size // downsample, k + 1))) - control_points = np.r_[(x[0],)*(k+1), control_points, (x[-1],)*(k+1)] if fixed_control_points is not None and len(fixed_control_points) > 0: good_fcp = fixed_control_points[(fixed_control_points > lo_lim) & (fixed_control_points < hi_lim)] @@ -332,13 +331,15 @@ def spline_through_data(x, y, k=2, grace_intv=1000., smoothing_factor=0.001, flag_ctrl_pts = np.ones(control_points.size, dtype=bool) t_idxs = np.searchsorted(x, control_points) for i, (id0, id1) in enumerate(zip(t_idxs[:-1], t_idxs[1:])): - if not id1 > id0 and np.any((x[id0:id1] > control_points[i]) & (x[id0:id1] < control_points[i + 1])): + if not id1 > id0 or not np.any((x[id0:id1+1] > control_points[i]) & (x[id0:id1 + 1] < control_points[i + 1])): # Schoenberg-Whitney condition is not satisfied, we need to remove this control point log.info(f"No data between {control_points[i]} and {control_points[i + 1]}, removing control point.") flag_ctrl_pts[i] = False control_points = control_points[flag_ctrl_pts] + control_points = np.r_[(x[0]- 1,)*(k+1), control_points, (x[-1] + 1,)*(k+1)] + detrend_fun = make_lsq_spline( x, y, t=control_points, k=k) From 349d2f75e87c12a405a00e0a7ebb38f5eedd8774 Mon Sep 17 00:00:00 2001 From: Matteo Bachetti Date: Fri, 27 Mar 2026 11:44:31 +0100 Subject: [PATCH 10/34] Fix issue with missing local bad point database --- nuclockutils/nustarclock.py | 19 +++++++++---------- 1 file changed, 9 insertions(+), 10 deletions(-) diff --git a/nuclockutils/nustarclock.py b/nuclockutils/nustarclock.py index faeaacf..2d1ed91 100644 --- a/nuclockutils/nustarclock.py +++ b/nuclockutils/nustarclock.py @@ -37,14 +37,15 @@ datadir = os.path.join(curdir, 'data') -def get_bad_points_db(db_file='BAD_POINTS_DB.dat'): - if not os.path.exists(db_file): +def get_bad_points_db(db_file=None): + if db_file is None or not os.path.exists(db_file): + log.warning(f"No local bad points database found. Using the default one.") db_file = os.path.join(datadir, 'BAD_POINTS_DB.dat') - + log.info(f"Reading bad points from {db_file}") return np.genfromtxt(db_file, dtype=np.longdouble) -def flag_bad_points(all_data, db_file='BAD_POINTS_DB.dat'): +def flag_bad_points(all_data, db_file=None): """ Examples @@ -57,13 +58,11 @@ def flag_bad_points(all_data, db_file='BAD_POINTS_DB.dat'): >>> bool(np.all(all_data['flag'] == [False, False, False, True, False])) True """ - if not os.path.exists(db_file): - all_data['flag'] = np.zeros(len(all_data), dtype=bool) - return all_data + log.info("Flagging bad points...") intv = [all_data['met'][0] - 0.5, all_data['met'][-1] + 0.5] - ALL_BAD_POINTS = np.genfromtxt(db_file) + ALL_BAD_POINTS = get_bad_points_db(db_file) ALL_BAD_POINTS.sort() ALL_BAD_POINTS = np.unique(ALL_BAD_POINTS) ALL_BAD_POINTS = ALL_BAD_POINTS[ @@ -121,11 +120,11 @@ def calculate_stats(all_data): print("-----------------------------------------------------------------------") -def load_and_flag_clock_table(clockfile="latest_clock.dat", shift_non_malindi=False): +def load_and_flag_clock_table(clockfile="latest_clock.dat", shift_non_malindi=False, db_file=None): clock_offset_table = load_clock_offset_table(clockfile, shift_non_malindi=shift_non_malindi) clock_offset_table = flag_bad_points( - clock_offset_table, db_file='BAD_POINTS_DB.dat') + clock_offset_table, db_file=db_file) return clock_offset_table From 528929f739c095afecad6817e60219424f16de7d Mon Sep 17 00:00:00 2001 From: Matteo Bachetti Date: Fri, 27 Mar 2026 20:06:49 +0100 Subject: [PATCH 11/34] Do look for local bad points file if it exists --- nuclockutils/nustarclock.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/nuclockutils/nustarclock.py b/nuclockutils/nustarclock.py index 2d1ed91..2bf2167 100644 --- a/nuclockutils/nustarclock.py +++ b/nuclockutils/nustarclock.py @@ -38,9 +38,13 @@ def get_bad_points_db(db_file=None): - if db_file is None or not os.path.exists(db_file): + if db_file is None: + db_file = "BAD_POINTS_DB.dat" + + if not os.path.exists(db_file): log.warning(f"No local bad points database found. Using the default one.") db_file = os.path.join(datadir, 'BAD_POINTS_DB.dat') + log.info(f"Reading bad points from {db_file}") return np.genfromtxt(db_file, dtype=np.longdouble) From 43c4917c0debcf294a2b17b9f624f9a92edb532f Mon Sep 17 00:00:00 2001 From: Matteo Bachetti Date: Fri, 27 Mar 2026 20:11:46 +0100 Subject: [PATCH 12/34] avoid inconsistencies in name of bad point file --- nuclockutils/clean_clock/app.py | 6 +++--- nuclockutils/nustarclock.py | 7 ++++--- 2 files changed, 7 insertions(+), 6 deletions(-) diff --git a/nuclockutils/clean_clock/app.py b/nuclockutils/clean_clock/app.py index e177d9d..907466b 100644 --- a/nuclockutils/clean_clock/app.py +++ b/nuclockutils/clean_clock/app.py @@ -19,7 +19,7 @@ from nuclockutils.nustarclock import load_temptable, load_freq_changes, \ load_and_flag_clock_table, find_good_time_intervals, calculate_stats, \ - eliminate_trends_in_residuals + eliminate_trends_in_residuals, _BAD_POINTS_FILE import dash from dash import dcc @@ -524,12 +524,12 @@ def refresh_output(n_cl_refresh, n_cl_bad, n_cl_good, n_cl_recalc, selectedData) log.info("Removing point(s) from bad clock offset database") ALL_BAD_POINTS = eliminate_array_from_array( ALL_BAD_POINTS, NEW_BAD_POINTS) - np.savetxt('BAD_POINTS_DB.dat', ALL_BAD_POINTS, fmt='%d') + np.savetxt(_BAD_POINTS_FILE, ALL_BAD_POINTS, fmt='%d') elif who_triggered == 'bad-data-button' and len(NEW_BAD_POINTS) > 0: log.info("Adding point(s) to bad clock offset database") ALL_BAD_POINTS = merge_and_sort_arrays( ALL_BAD_POINTS, NEW_BAD_POINTS) - np.savetxt('BAD_POINTS_DB.dat', ALL_BAD_POINTS, fmt='%d') + np.savetxt(_BAD_POINTS_FILE, ALL_BAD_POINTS, fmt='%d') else: log.info("Refreshing plot") CURRENT_AXES = default_axes() diff --git a/nuclockutils/nustarclock.py b/nuclockutils/nustarclock.py index 2bf2167..fa64a71 100644 --- a/nuclockutils/nustarclock.py +++ b/nuclockutils/nustarclock.py @@ -29,7 +29,8 @@ import holoviews as hv from holoviews.operation.datashader import datashade from holoviews import opts -# import matplotlib.pyplot as plt + +_BAD_POINTS_FILE = "BAD_POINTS_DB.dat" hv.extension('bokeh') @@ -39,11 +40,11 @@ def get_bad_points_db(db_file=None): if db_file is None: - db_file = "BAD_POINTS_DB.dat" + db_file = _BAD_POINTS_FILE if not os.path.exists(db_file): log.warning(f"No local bad points database found. Using the default one.") - db_file = os.path.join(datadir, 'BAD_POINTS_DB.dat') + db_file = os.path.join(datadir, _BAD_POINTS_FILE) log.info(f"Reading bad points from {db_file}") return np.genfromtxt(db_file, dtype=np.longdouble) From 92cdb3f51519cb3b8eb7f857cb33ceb8d761e790 Mon Sep 17 00:00:00 2001 From: Matteo Bachetti Date: Fri, 27 Mar 2026 20:23:41 +0100 Subject: [PATCH 13/34] Avoid reloading of app in debug mode --- nuclockutils/clean_clock/app.py | 4 ++-- nuclockutils/nustarclock.py | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/nuclockutils/clean_clock/app.py b/nuclockutils/clean_clock/app.py index 907466b..c6150d2 100644 --- a/nuclockutils/clean_clock/app.py +++ b/nuclockutils/clean_clock/app.py @@ -644,10 +644,10 @@ def main(args=None): print("Creating app") app = create_app() try: - app.run(debug=True) + app.run(debug=True, use_reloader=False) except Exception as e: # Compatibility with old versions - app.run_server(debug=True) + app.run_server(debug=True, use_reloader=False) diff --git a/nuclockutils/nustarclock.py b/nuclockutils/nustarclock.py index fa64a71..e6bdff9 100644 --- a/nuclockutils/nustarclock.py +++ b/nuclockutils/nustarclock.py @@ -507,7 +507,7 @@ def read_clock_offset_table(clockoffset_file=None, shift_non_malindi=False): clock_offset_table.remove_row(len(clock_offset_table) - 1) clock_offset_table['flag'] = np.zeros(len(clock_offset_table), dtype=bool) - log.info("Flagging bad points...") + log.info("Flagging bad points in clock offset table...") ALL_BAD_POINTS = get_bad_points_db() for b in ALL_BAD_POINTS: From a3138d9e669b688132e0972209e3d0e28f14f753 Mon Sep 17 00:00:00 2001 From: Matteo Bachetti Date: Wed, 1 Apr 2026 12:06:25 +0200 Subject: [PATCH 14/34] Avoid longdoubles. Also, do not run in debug mode, just fail --- nuclockutils/clean_clock/app.py | 58 ++++++++++++++------------------- 1 file changed, 25 insertions(+), 33 deletions(-) diff --git a/nuclockutils/clean_clock/app.py b/nuclockutils/clean_clock/app.py index c6150d2..e1b59a8 100644 --- a/nuclockutils/clean_clock/app.py +++ b/nuclockutils/clean_clock/app.py @@ -271,20 +271,23 @@ def plot_dash(all_data, table_new, gti, all_nustar_obs, continue if bti[0] > all_data['met'][-1]: continue - shapes.append(dict(type="rect", + shapes.append( + dict( + type="rect", # x-reference is assigned to the x-values xref="x", # y-reference is assigned to the plot paper [0,1] yref="paper", - x0=max(bti[0], all_data['met'][0]), + x0=float(max(bti[0], all_data["met"][0])), y0=0, - x1=min(bti[1], all_data['met'][-1]), + x1=float(min(bti[1], all_data["met"][-1])), y1=1, fillcolor="LightSalmon", opacity=0.5, layer="below", line_width=0, - )) + ) + ) if axis_ranges is not None: if 'xaxis.range[0]' in axis_ranges: @@ -552,37 +555,26 @@ def display_selected_data(selectedData): def create_temperature_timeseries(x, y, axis_type='linear'): return { - 'data': [dict( - x=x, - y=y, - mode='lines' - )], - 'layout': { - 'height': 300, - 'yaxis': {'title': 'TCXO Temperature', - 'type': 'linear' if axis_type == 'Linear' else 'log'}, - 'xaxis': {'title': 'met', 'showgrid': False, - 'margin':{'t': 20}} - - } + "data": [dict(x=x.astype(float), y=y.astype(float), mode="lines")], + "layout": { + "height": 300, + "yaxis": { + "title": "TCXO Temperature", + "type": "linear" if axis_type == "Linear" else "log", + }, + "xaxis": {"title": "met", "showgrid": False, "margin": {"t": 20}}, + }, } def create_temperature_gradient_timeseries(x, y, axis_type='linear'): return { - 'data': [dict( - x=x, - y=y, - mode='lines' - )], - 'layout': { - 'height': 300, - 'yaxis': {'title': 'TCXO Temp Gradient', - 'type': 'linear'}, - 'xaxis': {'title': 'met', 'showgrid': False, - 'margin':{'t': 20}} - - } + "data": [dict(x=x.astype(float), y=y.astype(float), mode="lines")], + "layout": { + "height": 300, + "yaxis": {"title": "TCXO Temp Gradient", "type": "linear"}, + "xaxis": {"title": "met", "showgrid": False, "margin": {"t": 20}}, + }, } @@ -644,10 +636,10 @@ def main(args=None): print("Creating app") app = create_app() try: - app.run(debug=True, use_reloader=False) + app.run() except Exception as e: # Compatibility with old versions - app.run_server(debug=True, use_reloader=False) + app.run_server() @@ -658,4 +650,4 @@ def main(args=None): FREQFILE = os.path.join(datadir, 'sample_freq.dat') TEMPFILE = os.path.join(datadir, 'tcxo_tmp_sample.csv') - main() \ No newline at end of file + main() From 1d06a11c777d42b317033ceadab685e1a20f5e18 Mon Sep 17 00:00:00 2001 From: Matteo Bachetti Date: Wed, 1 Apr 2026 12:52:03 +0200 Subject: [PATCH 15/34] Update standard deviation threshold in example --- nuclockutils/utils.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/nuclockutils/utils.py b/nuclockutils/utils.py index 2518de8..c32c74f 100644 --- a/nuclockutils/utils.py +++ b/nuclockutils/utils.py @@ -309,7 +309,7 @@ def spline_through_data(x, y, k=2, grace_intv=1000., smoothing_factor=0.001, >>> x = np.arange(1000) >>> y = np.random.normal(x * 0.1, 0.01) >>> fun = spline_through_data(x, y, grace_intv=10.) - >>> bool(np.std(y - fun(x)) < 0.01) + >>> bool(np.std(y - fun(x)) < 0.02) # 2 sigma True """ lo_lim, hi_lim = x[0], x[-1] From d41033967df45e8ba03e9eae91d76bd0695d0f48 Mon Sep 17 00:00:00 2001 From: Matteo Bachetti Date: Thu, 2 Apr 2026 01:02:00 +0200 Subject: [PATCH 16/34] Allow better debugging --- nuclockutils/clean_clock/app.py | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/nuclockutils/clean_clock/app.py b/nuclockutils/clean_clock/app.py index e1b59a8..c84e202 100644 --- a/nuclockutils/clean_clock/app.py +++ b/nuclockutils/clean_clock/app.py @@ -608,6 +608,7 @@ def main(args=None): global FREQFILE global MODELVERSION import argparse + import logging description = ('Clean clock offset measurements with an handy web ' 'interface.') parser = argparse.ArgumentParser(description=description) @@ -623,6 +624,8 @@ def main(args=None): "file in the auxil/directory " "or the tp_tcxo*.csv file)") parser.add_argument("--temperature-model-version", default=None, help="Temperature model version") + parser.add_argument("--devel", action='store_true', help="Enable development mode") + parser.add_argument("--debug", action='store_true', help="Enable debug logging") args = parser.parse_args(args) if args.temperature_file is not None: TEMPFILE = args.temperature_file @@ -632,14 +635,16 @@ def main(args=None): FREQFILE = args.frequency_file if args.temperature_model_version is not None: MODELVERSION = args.temperature_model_version + if args.debug: + log.setLevel(logging.DEBUG) print("Creating app") app = create_app() try: - app.run() + app.run(debug=args.devel, use_reloader=False) except Exception as e: # Compatibility with old versions - app.run_server() + app.run_server(debug=args.devel, use_reloader=False) From 8a6364de63322443c05c37bd69d7fd7fa160cfe0 Mon Sep 17 00:00:00 2001 From: Matteo Bachetti Date: Thu, 2 Apr 2026 01:03:28 +0200 Subject: [PATCH 17/34] Log whenever data points are eliminated --- nuclockutils/nustarclock.py | 119 +++++++++++++++++++++++++++++++----- nuclockutils/utils.py | 1 - 2 files changed, 103 insertions(+), 17 deletions(-) diff --git a/nuclockutils/nustarclock.py b/nuclockutils/nustarclock.py index e6bdff9..632db79 100644 --- a/nuclockutils/nustarclock.py +++ b/nuclockutils/nustarclock.py @@ -64,10 +64,9 @@ def flag_bad_points(all_data, db_file=None): True """ - log.info("Flagging bad points...") - intv = [all_data['met'][0] - 0.5, all_data['met'][-1] + 0.5] ALL_BAD_POINTS = get_bad_points_db(db_file) + log.info("Sorting and filtering bad points...") ALL_BAD_POINTS.sort() ALL_BAD_POINTS = np.unique(ALL_BAD_POINTS) ALL_BAD_POINTS = ALL_BAD_POINTS[ @@ -75,6 +74,7 @@ def flag_bad_points(all_data, db_file=None): idxs = all_data['met'].searchsorted(ALL_BAD_POINTS) + log.info("Flagging bad points...") if 'flag' in all_data.colnames: mask = np.array(all_data['flag'], dtype=bool) else: @@ -85,6 +85,7 @@ def flag_bad_points(all_data, db_file=None): continue mask[idx] = True all_data['flag'] = mask + log.info("Finished flagging bad points.") return all_data @@ -107,6 +108,8 @@ def find_good_time_intervals(temperature_table, good = lengths > 43200 if not np.all(good): log.info(f"Some GTIs are too short. cleaning up: {gtis[~good]}") + for g, is_good in zip(gtis, good): + log.info(f"GTI: {g[0]}-{g[1]} ({'OK' if is_good else 'too short'})") return gtis[good] @@ -139,7 +142,15 @@ def spline_detrending(clock_offset_table, temptable, outlier_cuts=None, temperature_is_present = tempcorr_idx < temptable['met'].size tempcorr_idx = tempcorr_idx[temperature_is_present] - clock_offset_table = clock_offset_table[temperature_is_present] + n_before = temperature_is_present.size + n_after = np.count_nonzero(temperature_is_present) + + if n_after < n_before: + log.info(f"spline_detrending: Filtering clock_offset_table:") + for eliminated in clock_offset_table[~temperature_is_present]['met']: + log.info(f" MET {eliminated} (beyond temperature table range)") + + clock_offset_table = clock_offset_table[temperature_is_present] clock_residuals = \ np.array(clock_offset_table['offset'] - @@ -163,8 +174,16 @@ def spline_detrending(clock_offset_table, temptable, outlier_cuts=None, do_not_flag = clock_mets > clock_mets.max() - one_month better_points[do_not_flag] = True - clock_offset_table = clock_offset_table[better_points] - clock_residuals = clock_residuals[better_points] + n_before = better_points.size + n_after = np.count_nonzero(better_points) + + if n_after < n_before: + log.info(f"spline_detrending: Filtering clock_offset_table outliers") + for eliminated in clock_offset_table[~better_points]['met']: + log.info(f" MET {eliminated}") + + clock_offset_table = clock_offset_table[better_points] + clock_residuals = clock_residuals[better_points] detrend_fun = spline_through_data( clock_offset_table['met'], clock_residuals, downsample=10, @@ -198,7 +217,15 @@ def eliminate_trends_in_residuals(temp_table, clock_offset_table, temperature_is_present = tempcorr_idx < temp_table['met'].size tempcorr_idx = tempcorr_idx[temperature_is_present] - clock_offset_table = clock_offset_table[temperature_is_present] + n_before = temperature_is_present.size + n_after = np.count_nonzero(temperature_is_present) + + if n_after < n_before: + log.info(f"eliminate_trends_in_residuals: Filtering clock_offset_table") + for eliminated in clock_offset_table[~temperature_is_present]['met']: + log.info(f" MET {eliminated} (beyond temperature table range)") + + clock_offset_table = clock_offset_table[temperature_is_present] clock_residuals = \ clock_offset_table['offset'] - temp_table['temp_corr'][tempcorr_idx] @@ -214,8 +241,17 @@ def eliminate_trends_in_residuals(temp_table, clock_offset_table, good = (clock_residuals == clock_residuals) & ~clock_offset_table['flag'] & use_for_interpol - clock_offset_table = clock_offset_table[good] - clock_residuals = clock_residuals[good] + n_before = good.size + n_after = np.count_nonzero(good) + + if n_after < n_before: + log.info("eliminate_trends_in_residuals: Filtering clock_offset_table to " + "trustworthy points (Malindi or during Malindi outage)") + for eliminated in clock_offset_table[~good]['met']: + log.debug(f" MET {eliminated} (bad point or not from Malindi)") + + clock_offset_table = clock_offset_table[good] + clock_residuals = clock_residuals[good] for g in gtis: log.info(f"Treating data from METs {g[0]}--{g[1]}") @@ -622,11 +658,20 @@ def read_freq_changes_table(freqchange_file=None, filter_bad=True): freq_changes_table.sort('met') freq_changes_table['mjd'] = sec_to_mjd(freq_changes_table['met']) - freq_changes_table.remove_row(len(freq_changes_table) - 1) + # freq_changes_table.remove_row(len(freq_changes_table) - 1) freq_changes_table['flag'] = \ np.abs(freq_changes_table['divisor'] - 2.400034e7) > 20 + if filter_bad: - freq_changes_table = freq_changes_table[~freq_changes_table['flag']] + n_before = freq_changes_table['flag'].size + n_after = n_before - np.count_nonzero(freq_changes_table['flag']) + + if n_after < n_before: + log.info(f"read_freq_changes_table: Filtering freq_changes_table") + for eliminated in freq_changes_table['met'][freq_changes_table['flag']]: + log.debug(f" MET {eliminated} (bad frequency change point)") + + freq_changes_table = freq_changes_table[~freq_changes_table['flag']] return freq_changes_table @@ -783,7 +828,14 @@ def read_temptable(temperature_file=None, mjdstart=None, mjdstop=None, else: good = np.diff(temptable['met']) > 0 good = np.concatenate((good, [True])) - temptable = temptable[good] + + n_before = good.size + n_after = np.count_nonzero(good) + if n_after < n_before: + log.info(f"read_temptable: Filtering temptable for bad time points") + for eliminated in temptable['met'][~good]: + log.debug(f" MET {eliminated} (non increasing time)") + temptable = temptable[good] temptable.meta['gti'] = temp_gtis window = np.median(1000 / np.diff(temptable['met'])) @@ -949,15 +1001,34 @@ def write_clock_file(self, filename=None, save_nodetrend=False, filename = 'new_clock_file.fits' table_new = self.temperature_correction_data clock_offset_table = self.clock_offset_table - clock_offset_table = clock_offset_table[ - clock_offset_table['met'] < table_new['met'][-1]] + + good_clock = clock_offset_table['met'] < table_new['met'][-1] + n_before = good_clock.size + n_after = np.count_nonzero(good_clock) + + if n_after < n_before: + log.info(f"write_clock_file: Filtering clock_offset_table:") + for eliminated in clock_offset_table[~good_clock]['met']: + log.info(f" MET {eliminated} " + f"(beyond table range)") + + clock_offset_table = clock_offset_table[good_clock] tempcorr_idx = np.searchsorted(table_new['met'], clock_offset_table['met']) temperature_is_present = tempcorr_idx < table_new['met'].size tempcorr_idx = tempcorr_idx[temperature_is_present] - clock_offset_table = clock_offset_table[temperature_is_present] + n_before = temperature_is_present.size + n_after = np.count_nonzero(temperature_is_present) + + if n_after < n_before: + log.info(f"write_clock_file: Filtering clock_offset_table:") + for eliminated in clock_offset_table[~temperature_is_present]['met']: + log.info(f" MET {eliminated} " + f"(beyond temperature range)") + + clock_offset_table = clock_offset_table[temperature_is_present] clock_residuals_detrend = clock_offset_table['offset'] - \ table_new['temp_corr'][tempcorr_idx] @@ -1229,7 +1300,16 @@ def plot_scatter(new_clock_table, clock_offset_table, shift_times=0, plt.show() yint = - yint - clock_offset_table = clock_offset_table[good_mets] + n_before = good_mets.size + n_after = np.count_nonzero(good_mets) + if n_after < n_before: + log.info("plot_scatter: Filtering clock_offset_table:") + for eliminated in clock_offset_table[~good_mets]['met']: + log.info(f" MET {eliminated} " + f"(beyond interpolation range)") + + clock_offset_table = clock_offset_table[good_mets] + clock_offset_table['offset'] -= shift_times clock_mets = clock_offset_table['met'] clock_mjds = clock_offset_table['mjd'] @@ -1592,7 +1672,14 @@ def temperature_correction_table(met_start, met_stop, "Interval not fully included in cached data. Recalculating.") else: good = (mets >= met_start - 86400) & (mets < met_stop + 86400) - return result_table[good] + n_before = good.size + n_after = np.count_nonzero(good) + if n_after < n_before: + log.info(f"temperature_correction_table: Filtering result_table:") + for eliminated in result_table[~good]['met']: + log.info(f" MET {eliminated} (beyond requested range)") + filtered_table = result_table[good] + return filtered_table if temptable is None or isinstance(temptable, six.string_types): mjdstart, mjdstop = sec_to_mjd(met_start), sec_to_mjd(met_stop) diff --git a/nuclockutils/utils.py b/nuclockutils/utils.py index c32c74f..4cf8e4a 100644 --- a/nuclockutils/utils.py +++ b/nuclockutils/utils.py @@ -318,7 +318,6 @@ def spline_through_data(x, y, k=2, grace_intv=1000., smoothing_factor=0.001, np.linspace(lo_lim + 2 * grace_intv, hi_lim - 2 * grace_intv, int(max(x.size // downsample, k + 1))) - if fixed_control_points is not None and len(fixed_control_points) > 0: good_fcp = fixed_control_points[(fixed_control_points > lo_lim) & (fixed_control_points < hi_lim)] if good_fcp.size > 0: From abcec8fb714e4ada30fe70ebbbb3a8815dad1cd2 Mon Sep 17 00:00:00 2001 From: Matteo Bachetti Date: Thu, 2 Apr 2026 01:04:27 +0200 Subject: [PATCH 18/34] Avoid needless elimination of last points; stabilize spline at the end --- nuclockutils/nustarclock.py | 4 ++-- nuclockutils/utils.py | 5 +++++ 2 files changed, 7 insertions(+), 2 deletions(-) diff --git a/nuclockutils/nustarclock.py b/nuclockutils/nustarclock.py index 632db79..b7dfb88 100644 --- a/nuclockutils/nustarclock.py +++ b/nuclockutils/nustarclock.py @@ -540,7 +540,7 @@ def read_clock_offset_table(clockoffset_file=None, shift_non_malindi=False): all_but_malindi = clock_offset_table['station'] != 'MLD' clock_offset_table['offset'][all_but_malindi] -= 0.0005 clock_offset_table['mjd'] = sec_to_mjd(clock_offset_table['met']) - clock_offset_table.remove_row(len(clock_offset_table) - 1) + # clock_offset_table.remove_row(len(clock_offset_table) - 1) clock_offset_table['flag'] = np.zeros(len(clock_offset_table), dtype=bool) log.info("Flagging bad points in clock offset table...") @@ -1691,7 +1691,7 @@ def temperature_correction_table(met_start, met_stop, read_freq_changes_table(freqchange_file=freqchange_file) allfreqtimes = np.array(freq_changes_table['met']) allfreqtimes = \ - np.concatenate([allfreqtimes, [allfreqtimes[-1] + 86400]]) + np.concatenate([allfreqtimes, [temptable['met'][-1]]]) met_intervals = list( zip(allfreqtimes[:-1], allfreqtimes[1:])) divisors = freq_changes_table['divisor'] diff --git a/nuclockutils/utils.py b/nuclockutils/utils.py index 4cf8e4a..2fef441 100644 --- a/nuclockutils/utils.py +++ b/nuclockutils/utils.py @@ -336,6 +336,11 @@ def spline_through_data(x, y, k=2, grace_intv=1000., smoothing_factor=0.001, flag_ctrl_pts[i] = False control_points = control_points[flag_ctrl_pts] + # Stabilize the spline by adding two extra control points at the beginning and end, + # with a large grace interval and zero value. This is to prevent the spline from + # diverging at the edges. + x = np.concatenate(([x[0] - 1000000], x, [x[-1] + 1000000])) + y = np.concatenate(([0], y, [0])) control_points = np.r_[(x[0]- 1,)*(k+1), control_points, (x[-1] + 1,)*(k+1)] From 6710523110a7e1404d9ec3c20dcb77f991b77b14 Mon Sep 17 00:00:00 2001 From: Matteo Bachetti Date: Thu, 2 Apr 2026 11:37:09 +0200 Subject: [PATCH 19/34] Add thorough documentation in nustarclock.py --- nuclockutils/nustarclock.py | 800 +++++++++++++++++++++++++++++++++--- 1 file changed, 732 insertions(+), 68 deletions(-) diff --git a/nuclockutils/nustarclock.py b/nuclockutils/nustarclock.py index b7dfb88..d48b0dc 100644 --- a/nuclockutils/nustarclock.py +++ b/nuclockutils/nustarclock.py @@ -1,9 +1,78 @@ +""" +NuSTAR Clock Correction Module +============================== + +This module provides tools for correcting the NuSTAR spacecraft clock drift +using temperature-dependent models. The onboard TCXO (Temperature Compensated +Crystal Oscillator) exhibits predictable drift correlated with temperature, +which this module models and corrects. + +Data Flow Overview +------------------ +1. **Temperature data** is read from housekeeping files (CSV, HDF5, or FITS) +2. **Clock offset measurements** from ground station passes are loaded +3. **Frequency divisor changes** (commanded clock adjustments) are loaded +4. Temperature-dependent clock drift is modeled using `clock_ppm_model` +5. Residuals between model and measurements are detrended with splines +6. Final correction is written as a CALDB-compatible FITS clock file + +Primary Table Schemas +--------------------- +**clock_offset_table** (from ground station measurements): + - ``uxt``: Unix timestamp of measurement + - ``met``: Mission Elapsed Time (seconds since MJDREF) + - ``offset``: Measured clock offset from ground truth (seconds) + - ``divisor``: Clock divisor value at measurement time + - ``station``: Ground station code (e.g. 'MLD'=Malindi, 'SNG'=Singapore, 'UHI'=Hawaii) + - ``mjd``: Modified Julian Date (computed) + - ``flag``: Boolean, True if measurement is flagged as bad + +**temptable** (temperature measurements): + - ``met``: Mission Elapsed Time (seconds) + - ``mjd``: Modified Julian Date + - ``temperature``: Raw TCXO temperature (Celsius) + - ``temperature_smooth``: Savitzky-Golay smoothed temperature + - ``temperature_smooth_gradient``: Time derivative of smoothed temperature + +**temptable with corrections** (after processing): + - All columns above, plus: + - ``temp_corr``: Temperature-model clock correction (seconds) + - ``temp_corr_raw``: Original temp_corr before trend removal + - ``temp_corr_trend``: Spline trend fitted to residuals + - ``temp_corr_detrend``: temp_corr + temp_corr_trend (final correction) + - ``std``: Rolling standard deviation of residuals + +**freq_changes_table** (commanded frequency divisor changes): + - ``uxt``: Unix timestamp of command + - ``met``: Mission Elapsed Time + - ``divisor``: New clock divisor value (~24000336) + - ``mjd``: Modified Julian Date (computed) + - ``flag``: Boolean, True if divisor value is anomalous + +Key Classes +----------- +- :class:`ClockCorrection`: Main class for computing clock corrections +- :class:`NuSTARCorr`: Applies corrections to event files + +Key Functions +------------- +- :func:`temperature_correction_table`: Builds the temperature-based correction +- :func:`eliminate_trends_in_residuals`: Removes systematic trends from residuals +- :func:`clock_ppm_model`: The physical model for clock drift vs temperature + +Notes +----- +- MET (Mission Elapsed Time) is in seconds since MJDREF = 55197.00076601852 +- The "Malindi problem" (2012-12-20 to 2013-02-08) requires special handling +- Clock jumps occur when the frequency divisor is commanded to change +""" + import glob import os import shutil from functools import lru_cache import traceback - +import warnings import numpy as np from astropy.table import Table, vstack import pandas as pd @@ -30,8 +99,45 @@ from holoviews.operation.datashader import datashade from holoviews import opts +# ============================================================================= +# Constants +# ============================================================================= + +# File paths _BAD_POINTS_FILE = "BAD_POINTS_DB.dat" +# Time constants (all in seconds unless noted) +SECONDS_PER_DAY = 86400 +HALF_DAY_SECONDS = 43200 # Minimum GTI duration for valid processing +SECONDS_PER_MONTH = SECONDS_PER_DAY * 30 +SECONDS_PER_YEAR = SECONDS_PER_DAY * 365.25 + +# Station timing offset: non-Malindi stations have a systematic 0.5 ms offset +# due to different signal processing pipelines +NON_MALINDI_OFFSET_SECONDS = 0.0005 + +# Clock hardware constants +NOMINAL_CLOCK_DIVISOR = 24000336 # Typical commanded divisor value +CLOCK_DIVISOR_TOLERANCE = 20 # Divisor values outside ±20 of 2.400034e7 are flagged + +# Known problematic time intervals (MET values) +# Malindi ground station outage: 2012-12-20 to 2013-02-08 +MALINDI_OUTAGE_INTERVALS = [ + (93681591, 98051312), # 2012/12/20 - 2013/02/08 Malindi outage + (357295300, 357972500), # 2021/04/28 - 2021/05/06 Malindi clock issues +] + +# Known clock jump times (MET) - when frequency divisor was commanded to change +KNOWN_CLOCK_JUMP_TIMES = np.array([ + 78708320, 79657575, 81043985, 82055671, 293346772, + 392200784, 394825882, 395304135, 407914525, 408299422 +]) + +# Reference epoch for absorption-desorption model +ABS_DES_REFERENCE_MET = 77509250 + +FIXED_CONTROL_POINTS = np.arange(291e6, 295e6, SECONDS_PER_DAY) + hv.extension('bokeh') curdir = os.path.abspath(os.path.dirname(__file__)) @@ -39,6 +145,23 @@ def get_bad_points_db(db_file=None): + """Load the database of known bad clock offset measurement times. + + Bad points are MET values where clock offset measurements are known + to be unreliable (e.g., due to ground station issues, spacecraft + anomalies, etc.). + + Parameters + ---------- + db_file : str, optional + Path to bad points database file. If None, uses the default + BAD_POINTS_DB.dat in the package data directory. + + Returns + ------- + bad_points : np.ndarray + Array of MET values (longdouble) marking bad measurement times. + """ if db_file is None: db_file = _BAD_POINTS_FILE @@ -51,7 +174,23 @@ def get_bad_points_db(db_file=None): def flag_bad_points(all_data, db_file=None): - """ + """Flag known bad measurements in a clock offset table. + + Adds or updates a 'flag' column in the table, setting True for + measurements at times listed in the bad points database. + + Parameters + ---------- + all_data : astropy.table.Table + Table with clock offset data. Must contain a 'met' column. + May already have a 'flag' column (which will be updated). + db_file : str, optional + Path to bad points database. If None, uses default. + + Returns + ------- + all_data : astropy.table.Table + Input table with 'flag' column added/updated. Examples -------- @@ -89,8 +228,27 @@ def flag_bad_points(all_data, db_file=None): return all_data -def find_good_time_intervals(temperature_table, - clock_jump_times=None): +def find_good_time_intervals(temperature_table, clock_jump_times=None): + """Identify Good Time Intervals (GTIs) for clock correction processing. + + GTIs are contiguous time ranges where: + 1. Temperature data is available (no gaps > 600s) + 2. No clock frequency jumps occur within the interval + 3. Duration is at least half a day (43200 seconds) + + Parameters + ---------- + temperature_table : astropy.table.Table + Temperature table with 'met' column. May have 'gti' in metadata. + clock_jump_times : array-like, optional + MET values where clock frequency divisor changed. Each jump + creates a GTI boundary. + + Returns + ------- + gtis : np.ndarray + Array of shape (N, 2) with [start, stop] MET for each GTI. + """ start_time = temperature_table['met'][0] stop_time = temperature_table['met'][-1] @@ -115,6 +273,17 @@ def find_good_time_intervals(temperature_table, def calculate_stats(all_data): + """Calculate and print statistics on clock correction residuals. + + Computes the Median Absolute Deviation (MAD) and rolling standard + deviation of the detrended residuals to assess correction quality. + + Parameters + ---------- + all_data : astropy.table.Table + Table with 'residual_detrend' column containing the residuals + between measured clock offsets and the temperature model. + """ log.info("Calculating statistics") r_std = residual_roll_std(all_data['residual_detrend']) @@ -129,6 +298,25 @@ def calculate_stats(all_data): def load_and_flag_clock_table(clockfile="latest_clock.dat", shift_non_malindi=False, db_file=None): + """Load clock offset table and flag known bad measurements. + + Convenience function combining load_clock_offset_table and flag_bad_points. + + Parameters + ---------- + clockfile : str + Path to clock offset data file. + shift_non_malindi : bool + If True, subtract NON_MALINDI_OFFSET_SECONDS from non-Malindi + station measurements to align them with Malindi. + db_file : str, optional + Path to bad points database. + + Returns + ------- + clock_offset_table : astropy.table.Table + Clock offset table with 'flag' column populated. + """ clock_offset_table = load_clock_offset_table(clockfile, shift_non_malindi=shift_non_malindi) clock_offset_table = flag_bad_points( @@ -138,6 +326,37 @@ def load_and_flag_clock_table(clockfile="latest_clock.dat", shift_non_malindi=Fa def spline_detrending(clock_offset_table, temptable, outlier_cuts=None, fixed_control_points=None): + """Fit a spline to clock residuals and add detrended correction to temptable. + + This function: + 1. Computes residuals between measured clock offsets and temperature model + 2. Optionally removes outliers based on deviation from median-filtered signal + 3. Fits a spline through the residuals to capture systematic trends + 4. Adds the trend to the temperature correction to improve accuracy + + Parameters + ---------- + clock_offset_table : astropy.table.Table + Clock offset measurements with columns: 'met', 'offset', 'flag'. + Filtered to only include times within temptable range. + temptable : astropy.table.Table + Temperature table with 'met' and 'temp_corr' columns. + Modified in-place to add: 'std', 'temp_corr_trend', 'temp_corr_detrend'. + outlier_cuts : list of float, optional + Threshold values (in seconds) for iterative outlier removal. + Measurements deviating more than these thresholds from the median + are flagged, except for data in the most recent month. + fixed_control_points : array-like, optional + MET values where spline control points should be placed. + + Returns + ------- + temptable : astropy.table.Table + Input table with added columns: + - 'std': Rolling standard deviation of residuals + - 'temp_corr_trend': Spline fit to residuals + - 'temp_corr_detrend': temp_corr + temp_corr_trend (improved correction) + """ tempcorr_idx = np.searchsorted(temptable['met'], clock_offset_table['met']) temperature_is_present = tempcorr_idx < temptable['met'].size tempcorr_idx = tempcorr_idx[temperature_is_present] @@ -170,8 +389,7 @@ def spline_detrending(clock_offset_table, temptable, outlier_cuts=None, outlier_cuts[0]) better_points[better_points] = ~wh # Eliminate too recent flags, in the last month of solution. - one_month = 86400 * 30 - do_not_flag = clock_mets > clock_mets.max() - one_month + do_not_flag = clock_mets > clock_mets.max() - SECONDS_PER_MONTH better_points[do_not_flag] = True n_before = better_points.size @@ -204,17 +422,55 @@ def spline_detrending(clock_offset_table, temptable, outlier_cuts=None, return temptable -def eliminate_trends_in_residuals(temp_table, clock_offset_table, +def eliminate_trends_in_residuals(temptable, clock_offset_table, gtis, debug=False, fixed_control_points=None): + """Remove systematic trends from temperature-model residuals within each GTI. + + For each Good Time Interval (GTI), this function: + 1. Extracts clock offset measurements and computes residuals vs temp model + 2. Fits a low-order robust polynomial to the residuals + 3. Adds the polynomial correction to the temperature model + 4. For Bad Time Intervals (gaps between GTIs), interpolates the correction + 5. Finally applies spline_detrending for fine-scale adjustments + + The Malindi ground station is preferred for interpolation due to its + more reliable timing. During known Malindi outages, other stations + are used with a 0.5 ms offset correction. - # good = clock_offset_table['met'] < np.max(temp_table['met']) + Parameters + ---------- + temptable : astropy.table.Table + Temperature correction table with 'met' and 'temp_corr' columns. + Modified in-place. Will have 'temp_corr_raw' added to preserve original. + clock_offset_table : astropy.table.Table + Clock offset measurements with 'met', 'offset', 'flag', 'station' columns. + gtis : np.ndarray + Array of shape (N, 2) with [start, stop] MET for each GTI. + debug : bool, optional + If True, save diagnostic plots for each GTI. + fixed_control_points : array-like, optional + MET values for spline control points in final detrending step. + + Returns + ------- + temptable : astropy.table.Table + Modified temperature table with improved 'temp_corr' and additional + diagnostic columns from spline_detrending. + + Notes + ----- + The function distinguishes between: + - GTIs: Intervals with good temperature data, processed with polynomial fits + - BTIs: Gaps between GTIs, handled by interpolation from nearby good data + """ + # good = clock_offset_table['met'] < np.max(temptable['met']) # clock_offset_table = clock_offset_table[good] - temp_table['temp_corr_raw'] = temp_table['temp_corr'] + temptable['temp_corr_raw'] = temptable['temp_corr'] - tempcorr_idx = np.searchsorted(temp_table['met'], + tempcorr_idx = np.searchsorted(temptable['met'], clock_offset_table['met']) - temperature_is_present = tempcorr_idx < temp_table['met'].size + temperature_is_present = tempcorr_idx < temptable['met'].size tempcorr_idx = tempcorr_idx[temperature_is_present] n_before = temperature_is_present.size @@ -228,7 +484,7 @@ def eliminate_trends_in_residuals(temp_table, clock_offset_table, clock_offset_table = clock_offset_table[temperature_is_present] clock_residuals = \ - clock_offset_table['offset'] - temp_table['temp_corr'][tempcorr_idx] + clock_offset_table['offset'] - temptable['temp_corr'][tempcorr_idx] # Only use for interpolation Malindi points; however, during the Malindi # problem in 2013, use the other data for interpolation but subtracting @@ -265,9 +521,9 @@ def eliminate_trends_in_residuals(temp_table, clock_offset_table, continue temp_idx_start, temp_idx_end = \ - np.searchsorted(temp_table['met'], g) + np.searchsorted(temptable['met'], g) - table_new = temp_table[temp_idx_start:temp_idx_end] + table_new = temptable[temp_idx_start:temp_idx_end] cltable_new = clock_offset_table[cl_idx_start:cl_idx_end] met = cltable_new['met'] @@ -287,12 +543,6 @@ def eliminate_trends_in_residuals(temp_table, clock_offset_table, log.info(f"Fitting a polinomial of order {poly_order}") p = robust_poly_fit(met_rescale, residuals, order=poly_order, p0=p0) - # if poly_order >=2: - import matplotlib.pyplot as plt - # plt.figure() - # plt.plot(met_rescale, residuals) - # plt.plot(met_rescale, p(met_rescale)) - # plt.plot(met_rescale, m * (met_rescale) + q) table_mets_rescale = (table_new['met'] - met0) / (met[-1] - met0) corr = p(table_mets_rescale) @@ -301,8 +551,6 @@ def eliminate_trends_in_residuals(temp_table, clock_offset_table, m = (sub_residuals[-1] - sub_residuals[0]) / (met_rescale[-1] - met_rescale[0]) q = sub_residuals[0] - # plt.plot(table_mets_rescale, corr + m * (table_mets_rescale - met_rescale[0]) + q, lw=2) - corr = corr + m * (table_mets_rescale - met_rescale[0]) + q table_new['temp_corr'] += corr @@ -320,7 +568,7 @@ def eliminate_trends_in_residuals(temp_table, clock_offset_table, bti_list = [[0, 77674700]] bti_list += [[g0, g1] for g0, g1 in zip(gtis[:-1, 1], gtis[1:, 0])] - bti_list += [[gtis[-1, 1], max(clock_offset_table['met'][-1], temp_table['met'][-1]) + 10000]] + bti_list += [[gtis[-1, 1], max(clock_offset_table['met'][-1], temptable['met'][-1]) + 10000]] btis = np.array(bti_list) # Interpolate the solution along bad time intervals for g in btis: @@ -328,26 +576,26 @@ def eliminate_trends_in_residuals(temp_table, clock_offset_table, log.info(f"Treating bad data from METs {start}--{stop}") temp_idx_start, temp_idx_end = \ - np.searchsorted(temp_table['met'], g) + np.searchsorted(temptable['met'], g) cl_idx_start, cl_idx_end = \ np.searchsorted(clock_offset_table['met'], g) if temp_idx_end - temp_idx_start == 0 and \ - temp_idx_end < len(temp_table) and temp_idx_start > 0: + temp_idx_end < len(temptable) and temp_idx_start > 0: log.info("No temperature measurements in this interval") continue else: - table_new = temp_table[temp_idx_start:temp_idx_end] + table_new = temptable[temp_idx_start:temp_idx_end] - last_good_tempcorr = temp_table['temp_corr'][temp_idx_start - 1] - last_good_time = temp_table['met'][temp_idx_start - 1] + last_good_tempcorr = temptable['temp_corr'][temp_idx_start - 1] + last_good_time = temptable['met'][temp_idx_start - 1] local_clockoff = clock_offset_table[max(cl_idx_start - 1, 0):cl_idx_end + 1] clock_off = local_clockoff['offset'] clock_tim = local_clockoff['met'] - if temp_idx_end < temp_table['temp_corr'].size: - next_good_tempcorr = temp_table['temp_corr'][temp_idx_end + 1] - next_good_time = temp_table['met'][temp_idx_end + 1] + if temp_idx_end < temptable['temp_corr'].size: + next_good_tempcorr = temptable['temp_corr'][temp_idx_end + 1] + next_good_time = temptable['met'][temp_idx_end + 1] clock_off = np.concatenate( ([last_good_tempcorr], clock_off, [next_good_tempcorr])) clock_tim = np.concatenate( @@ -380,14 +628,14 @@ def eliminate_trends_in_residuals(temp_table, clock_offset_table, log.info("Final detrending...") table_new = spline_detrending( - clock_offset_table, temp_table, + clock_offset_table, temptable, outlier_cuts=[-0.002, -0.001], fixed_control_points=fixed_control_points) return table_new def residual_roll_std(residuals, window=30): - """ + """Calculate the rolling standard deviation of clock residuals. Examples -------- @@ -518,7 +766,8 @@ def _look_for_freq_change_file(): def read_clock_offset_table(clockoffset_file=None, shift_non_malindi=False): - """ + """Read the clock offset table from a file and prepare it for processing. + Parameters ---------- clockoffset_file : str @@ -568,7 +817,25 @@ def read_clock_offset_table(clockoffset_file=None, shift_non_malindi=False): def no_jump_gtis(start_time, stop_time, clock_jump_times=None): - """ + """Create GTIs that avoid clock frequency jump times. + + Splits a time range into segments at each clock jump, so that + no segment contains a frequency divisor change. + + Parameters + ---------- + start_time : float + Start of time range (MET). + stop_time : float + End of time range (MET). + clock_jump_times : array-like, optional + MET values where clock jumps occur. + + Returns + ------- + gtis : np.ndarray + Array of [start, stop] pairs, one per segment. + Examples -------- >>> gtis = no_jump_gtis(0, 3, [1, 1.1]) @@ -592,7 +859,23 @@ def no_jump_gtis(start_time, stop_time, clock_jump_times=None): def temperature_gtis(temperature_table, max_distance=600): - """ + """Identify GTIs based on temperature data continuity. + + Creates GTIs where consecutive temperature measurements are within + max_distance seconds of each other. Gaps larger than this indicate + missing data and create GTI boundaries. + + Parameters + ---------- + temperature_table : astropy.table.Table + Table with 'met' column of measurement times. + max_distance : float, optional + Maximum allowed gap (seconds) between measurements. Default 600s. + + Returns + ------- + gtis : np.ndarray + Array of shape (N, 2) with [start, stop] MET for each GTI. Examples -------- @@ -630,16 +913,27 @@ def temperature_gtis(temperature_table, max_distance=600): def read_freq_changes_table(freqchange_file=None, filter_bad=True): - """Read the table with the list of commanded divisor frequencies. + """Read the table of commanded clock frequency divisor changes. + + The clock divisor determines the tick rate of the spacecraft clock. + This table records when the divisor was commanded to change, which + causes discontinuities in the clock drift that must be handled + separately. Parameters ---------- - freqchange_file : str - e.g. 'nustar_freq_changes-2018-10-30.dat' + freqchange_file : str, optional + Path to frequency changes file (e.g., 'nustar_freq_changes-2018-10-30.dat'). + If None, uses the latest file in the data directory. + filter_bad : bool, optional + If True (default), remove entries with anomalous divisor values + (more than ±20 from 2.400034e7). Returns ------- - freq_changes_table : `astropy.table.Table` object + freq_changes_table : astropy.table.Table + Table with columns: 'uxt', 'met', 'divisor', 'mjd', 'flag'. + Sorted by MET. Known bad entries are corrected per FREQ_CHANGE_DB. """ if freqchange_file is None: freqchange_file = _look_for_freq_change_file() @@ -676,12 +970,27 @@ def read_freq_changes_table(freqchange_file=None, filter_bad=True): return freq_changes_table -def _filter_table(tablefile, start_date=None, end_date=None, tmpfile='tmp.csv'): - try: - from datetime import timezone - except ImportError: - # Python 2 - import pytz as timezone +def _filter_csv_temperature_table(tablefile, start_date=None, end_date=None, tmpfile='tmp.csv'): + """Filter the CSV temperature table to only include rows within a specified date range. + + Parameters + ---------- + tablefile : str + Path to the input CSV temperature table file. + start_date : float, optional + Start date for filtering (MJD). + end_date : float, optional + End date for filtering (MJD). + tmpfile : str, optional + Path to the temporary filtered file. + + Returns + ------- + tmpfile : str + Path to the filtered CSV temperature table file. + """ + from datetime import timezone + if start_date is None: start_date = 0 if end_date is None: @@ -729,6 +1038,22 @@ def _filter_table(tablefile, start_date=None, end_date=None, tmpfile='tmp.csv'): def read_csv_temptable(mjdstart=None, mjdstop=None, temperature_file=None): + """Read the temperature table from a CSV file, optionally filtering by MJD range. + + Parameters + ---------- + mjdstart : float, optional + Start MJD for filtering the temperature table. If None, no lower bound is applied. + mjdstop : float, optional + Stop MJD for filtering the temperature table. If None, no upper bound is applied. + temperature_file : str, optional + Path to the CSV temperature table file. If None, the default file is used. + + Returns + ------- + temptable : Table + The filtered temperature table. + """ if mjdstart is not None or mjdstop is not None: mjdstart_use = mjdstart mjdstop_use = mjdstop @@ -737,7 +1062,7 @@ def read_csv_temptable(mjdstart=None, mjdstop=None, temperature_file=None): if mjdstop is not None: mjdstop_use += 10 log.info("Filtering table...") - tmpfile = _filter_table(temperature_file, + tmpfile = _filter_csv_temperature_table(temperature_file, start_date=mjdstart_use, end_date=mjdstop_use, tmpfile='tmp.csv') log.info("Done") @@ -763,6 +1088,22 @@ def read_csv_temptable(mjdstart=None, mjdstop=None, temperature_file=None): def read_saved_temptable(mjdstart=None, mjdstop=None, temperature_file='temptable.hdf5'): + """Read a previously saved temperature table from an HDF5 file. + + Parameters + ---------- + mjdstart : float, optional + Start MJD for filtering the temperature table. If None, no lower bound is applied. + mjdstop : float, optional + Stop MJD for filtering the temperature table. If None, no upper bound is applied. + temperature_file : str, optional + Path to the HDF5 temperature table file. If None, the default file is used. + + Returns + ------- + temptable : Table + The filtered temperature table. + """ table = Table.read(temperature_file) if mjdstart is None and mjdstop is None: return table @@ -783,6 +1124,18 @@ def read_saved_temptable(mjdstart=None, mjdstop=None, def read_fits_temptable(temperature_file): + """Read the temperature table from a FITS file, e.g. those from the HK data. + + Parameters + ---------- + temperature_file : str + Path to the FITS temperature table file. + + Returns + ------- + temptable : Table + The temperature table. + """ with fits.open(temperature_file) as hdul: temptable = Table.read(hdul['ENG_0x133']) temptable.rename_column('TIME', 'met') @@ -795,6 +1148,7 @@ def read_fits_temptable(temperature_file): def interpolate_temptable(temptable, dt=10): + """Interpolate the temperature table to a regular time grid with spacing dt.""" time = temptable['met'] temperature = temptable['temperature'] new_times = np.arange(time[0], time[-1], dt) @@ -804,6 +1158,29 @@ def interpolate_temptable(temptable, dt=10): def read_temptable(temperature_file=None, mjdstart=None, mjdstop=None, dt=None, gti_tolerance=600): + """Read the temperature table, handling different formats. + + Parameters + ---------- + temperature_file : str, optional + Path to the temperature table file. If None, the default file is used. + + Other parameters + ---------------- + mjdstart : float, optional + Start MJD for filtering the temperature table. If None, no lower bound is applied. + mjdstop : float, optional + Stop MJD for filtering the temperature table. If None, no upper bound is applied. + dt : float, optional + Time resolution for interpolation. If None, no interpolation is done. Default is None. + gti_tolerance : float, optional + Maximum gap (seconds) between temperature measurements to be considered a GTI. Default is + 600 seconds. + Returns + ------- + temptable : Table + The temperature table with columns 'met', 'temperature', 'mjd', and optionally 'temperature_smooth' and 'temperature_smooth_gradient'. + """ if temperature_file is None: temperature_file = _look_for_temptable() log.info(f"Reading temperature_information from {temperature_file}") @@ -881,6 +1258,51 @@ def load_clock_offset_table(clock_offset_file, shift_non_malindi=False): class ClockCorrection(): + """Main class for computing NuSTAR temperature-based clock corrections. + + This class orchestrates the full clock correction pipeline: + 1. Loads temperature data and clock offset measurements + 2. Computes temperature-dependent clock drift model + 3. Optionally adjusts the model using measured clock offsets + 4. Can write CALDB-compatible clock correction FITS files + + Parameters + ---------- + temperature_file : str + Path to temperature data file (CSV or HDF5). + mjdstart, mjdstop : float, optional + MJD range for the correction. If None, derived from data. + temperature_dt : float, optional + Time resolution for temperature interpolation. Default 10s. + adjust_absolute_timing : bool, optional + If True, adjust the temperature model using measured clock offsets + via eliminate_trends_in_residuals. Default False. + force_divisor : int, optional + Force a specific clock divisor value instead of reading from file. + label : str, optional + Label for output files. + additional_days : float, optional + Extra days of data to load beyond requested range. Default 2. + clock_offset_file : str, optional + Path to clock offset measurements file. + hdf_dump_file : str, optional + Path for caching intermediate results. + freqchange_file : str, optional + Path to frequency changes file. + spline_through_residuals : bool, optional + Reserved for future use. + + Attributes + ---------- + temptable : astropy.table.Table + Temperature data table. + clock_offset_table : astropy.table.Table + Clock offset measurements. + temperature_correction_data : astropy.table.Table + Computed correction table with 'met', 'temp_corr', etc. + gtis : np.ndarray + Good Time Intervals for processing. + """ def __init__(self, temperature_file, mjdstart=None, mjdstop=None, temperature_dt=10, adjust_absolute_timing=False, force_divisor=None, label="", additional_days=2, @@ -930,10 +1352,9 @@ def __init__(self, temperature_file, mjdstart=None, mjdstop=None, self.hdf_dump_file = hdf_dump_file self.plot_file = label + "_clock_adjustment.png" - self.clock_jump_times = \ - np.array([78708320, 79657575, 81043985, 82055671, 293346772, - 392200784, 394825882, 395304135,407914525, 408299422]) - self.fixed_control_points = np.arange(291e6, 295e6, 86400) + self.clock_jump_times = KNOWN_CLOCK_JUMP_TIMES + + self.fixed_control_points = FIXED_CONTROL_POINTS # Sum 30 seconds to avoid to exclude these points # from previous interval self.gtis = find_good_time_intervals( @@ -960,6 +1381,7 @@ def __init__(self, temperature_file, mjdstart=None, mjdstop=None, traceback.print_exc(file=fobj) def read_temptable(self, cache_temptable_name=None): + """Read the temperature table, using a cached version if available.""" if cache_temptable_name is not None and \ os.path.exists(cache_temptable_name): self.temptable = Table.read(cache_temptable_name, path='temptable') @@ -973,12 +1395,22 @@ def read_temptable(self, cache_temptable_name=None): self.temptable.write(cache_temptable_name, path='temptable') def temperature_correction_fun(self, adjust=False): + """Create an interpolating function for the temperature correction. + + Uses PCHIP interpolation to create a smooth function that can be evaluated + at any time. If adjust is True, uses the adjusted temperature correction data. + """ data = self.temperature_correction_data + if not adjust: + warnings.warn("Since the use of PCHIP interpolation, the adjust option is " + "ignored. The function will always use the adjusted temperature correction data.") + return PchipInterpolator(np.array(data['met']), np.array(data['temp_corr']), extrapolate=True) def adjust_temperature_correction(self): + """Adjust the temperature correction using measured clock offsets.""" table_new = eliminate_trends_in_residuals( self.temperature_correction_data, load_clock_offset_table( @@ -994,6 +1426,25 @@ def adjust_temperature_correction(self): def write_clock_file(self, filename=None, save_nodetrend=False, shift_times=0., highres=False): + """Write a CALDB-compatible clock correction FITS file with the temperature correction. + + + Parameters + ---------- + filename : str, optional + Path to the output FITS file. If None, defaults to 'new_clock_file.fits'. + save_nodetrend : bool, optional + If True, also save the non-detrended temperature correction in the FITS file. + shift_times : float, optional + Time shift (seconds) to apply to the clock offset correction. Default is 0. This can be used to align the correction with the measured clock offsets if needed. + highres : bool, optional + If True, save the clock correction at the full time resolution of the temperature table. If False (default), save a subsampled version with points every 100 seconds, and more points around known clock jump times. + + Returns + ------- + None + Writes the clock correction data to the specified FITS file. + """ from astropy.io.fits import Header, HDUList, BinTableHDU, PrimaryHDU from astropy.time import Time @@ -1237,6 +1688,20 @@ def write_clock_file(self, filename=None, save_nodetrend=False, def interpolate_clock_function(new_clock_table, mets): + """Interpolate the clock correction function at the given METs using the new clock table. + + Parameters + ---------- + new_clock_table : astropy.table.Table + The new clock correction table with columns 'TIME', 'CLOCK_OFF_CORR', and 'CLOCK_FREQ_CORR' + mets : array-like + The METs at which to interpolate the clock correction + + Returns + ------- + tuple + A tuple containing the interpolated clock corrections and a boolean mask indicating which METs were successfully interpolated + """ tab_times = new_clock_table['TIME'] good_mets = (mets > tab_times.min()) & (mets < tab_times.max()) mets = mets[good_mets] @@ -1254,6 +1719,26 @@ def interpolate_clock_function(new_clock_table, mets): def _analyze_residuals(filtered_data, nbins=101, fit_range=[-np.inf, np.inf]): + """Analyze the residuals by fitting a Gaussian to the histogram of the filtered data. + + Parameters + ---------- + filtered_data : array-like + The data to analyze, typically the residuals after detrending. + nbins : int, optional + The number of bins to use for the histogram. Default is 101. + fit_range : list, optional + The range of values to consider for fitting the Gaussian. Default is [-inf, inf]. + + Returns + ------- + edges: array-like + edges of the histogram bins, + frequencies: array-like + the normalized frequencies + info_dict: dict + a dictionary of statistics including the median, MAD, and fitted Gaussian parameters. + """ frequencies, edges = np.histogram( filtered_data, bins=np.linspace(-1000, 1500, nbins) @@ -1477,6 +1962,41 @@ def plot_scatter(new_clock_table, clock_offset_table, shift_times=0, class NuSTARCorr(): + """Apply clock correction to a NuSTAR event file. + + Convenience class that creates a ClockCorrection for a specific + observation and applies it to the event times. + + Parameters + ---------- + events_file : str + Path to input NuSTAR event file (FITS). + outfile : str, optional + Path for output corrected file. Default: input with '_tc' suffix. + adjust : bool, optional + If True, adjust temperature model using clock offsets. Default True. + force_divisor : int, optional + Force a specific clock divisor value. + temperature_file : str, optional + Path to temperature data. If None, uses default. + hdf_dump_file : str, optional + Path for caching intermediate results. + + Attributes + ---------- + clock_correction : ClockCorrection + The underlying correction calculator. + temperature_correction_fun : callable + Interpolation function for the correction. + tstart, tstop : float + Observation time range (MET). + + Examples + -------- + >>> corr = NuSTARCorr('nu12345_01_cl.evt') # doctest: +SKIP + >>> corr.apply_clock_correction() # doctest: +SKIP + 'nu12345_01_cl_tc.evt' + """ def __init__(self, events_file, outfile=None, adjust=True, force_divisor=None, temperature_file=None, hdf_dump_file=None): @@ -1510,11 +2030,19 @@ def __init__(self, events_file, outfile=None, self.clock_correction.temperature_correction_fun(adjust=adjust) def read_observation_info(self): + """Read TSTART and TSTOP from the event file header.""" with fits.open(self.events_file) as hdul: hdr = hdul[1].header self.tstart, self.tstop = hdr['TSTART'], hdr['TSTOP'] def apply_clock_correction(self): + """Apply the clock correction to the event times and write a new FITS file. + + Returns + ------- + str + Path to the output corrected FITS file. + """ events_file = self.events_file outfile = self.outfile @@ -1534,6 +2062,28 @@ def apply_clock_correction(self): def simpcumquad(x, y): + """Cumulative Simpson's rule integration for equally-spaced data. + + Computes the cumulative integral of y(x) using Simpson's rule, + returning the integrated value at each sample point. + + Parameters + ---------- + x : array-like + Independent variable values (must be equally spaced). + y : array-like + Dependent variable values to integrate. + + Returns + ------- + integral : np.ndarray + Cumulative integral values at each x point. + + Raises + ------ + ValueError + If x and y have different lengths or x is empty. + """ n = len(x) if (n != len(y)): @@ -1567,30 +2117,71 @@ def simpcumquad(x, y): def abs_des_fun(x, b0, b1, b2, b3, t0=77509250): - """An absorption-desorption function""" - x = (x - t0) / 86400. / 365.25 + """Absorption-desorption function for long-term clock drift. + + Models the long-term systematic drift of the TCXO due to material + outgassing and other aging effects. The drift follows a sum of + logarithmic terms representing different decay processes. + + Parameters + ---------- + x : array-like + Time values in MET (Mission Elapsed Time in seconds). + b0, b1, b2, b3 : float + Model coefficients for the two logarithmic components. + t0 : float, optional + Reference epoch in MET. Default is ABS_DES_REFERENCE_MET (77509250), + corresponding to early mission. + + Returns + ------- + drift : np.ndarray + Long-term drift contribution in ppm. + + Notes + ----- + The function computes: b0*ln(b1*t + 1) + b2*ln(b3*t + 1) + where t is time in years since t0. + """ + x = (x - t0) / SECONDS_PER_YEAR return b0 * np.log(b1 * x + 1) + b2 * np.log(b3 * x + 1) def clock_ppm_model(nustar_met, temperature, craig_fit=False, version=None, pars=None): - """Improved clock model + """Compute the clock frequency deviation model in parts-per-million. + + The NuSTAR TCXO frequency varies with both temperature and time. + This function combines: + 1. A quadratic temperature dependence around a reference temperature + 2. A long-term drift modeled by absorption-desorption functions Parameters ---------- - time : np.array of floats - Times in MET - temperature : np.array of floats - TCXO Temperatures in C - T0 : float - Reference temperature - ppm_vs_T_pars : list - parameters of the ppm-temperature relation - ppm_vs_T_pars : list - parameters of the ppm-time relation (long-term clock decay) - version : str - Version of the model to use. If None, use the latest - pars : dict - Parameters of the model, forced to be used + nustar_met : array-like + Mission Elapsed Time values in seconds. + temperature : array-like + TCXO temperature values in Celsius. + craig_fit : bool, optional + If True, use the original Craig Markwardt fit parameters. + version : str, optional + Model version identifier. If None, uses the latest version. + pars : dict, optional + Override model parameters. Expected keys: + - 'T0': Reference temperature + - 'offset': Constant ppm offset + - 'ppm_vs_T_pars': [linear_coeff, quadratic_coeff] + - 'ppm_vs_time_pars': [b0, b1, b2, b3] for abs_des_fun + + Returns + ------- + ppm : np.ndarray + Clock frequency deviation in parts-per-million. Positive values + mean the clock runs fast; negative means it runs slow. + + See Also + -------- + abs_des_fun : Long-term drift component + temperature_delay : Integrates ppm to get time delay """ if craig_fit: version = "craig" @@ -1616,6 +2207,37 @@ def temperature_delay(temptable, divisor, met_start=None, met_stop=None, debug=False, craig_fit=False, time_resolution=10, version=None, pars=None): + """Compute cumulative clock delay from temperature-dependent drift. + + Integrates the clock rate correction (from clock_ppm_model) over time + to get the cumulative time delay. Returns an interpolation function + that can evaluate the delay at any MET within the range. + + Parameters + ---------- + temptable : astropy.table.Table + Temperature table with 'met' and 'temperature_smooth' columns. + divisor : int + Clock frequency divisor value (nominally ~24000336). + met_start, met_stop : float, optional + Time range for computation. Defaults to temptable range. + debug : bool, optional + Reserved for debugging output. + craig_fit : bool, optional + Use original Craig Markwardt model parameters. + time_resolution : float, optional + Time step in seconds for integration. Default 10s. + version : str, optional + Model version identifier. + pars : dict, optional + Override model parameters. + + Returns + ------- + delay_function : scipy.interpolate.PchipInterpolator + Function that returns cumulative clock delay (seconds) for any MET. + Can extrapolate beyond the input range. + """ table_times = temptable['met'] if met_start is None: @@ -1660,6 +2282,48 @@ def temperature_correction_table(met_start, met_stop, craig_fit=False, version=None, pars=None): + """Build a table of temperature-based clock corrections. + + This is the main function for computing the temperature-dependent + clock correction. It: + 1. Loads or reads cached correction data if available + 2. Processes each time interval between frequency divisor changes + 3. Integrates the clock drift to get cumulative delay + 4. Handles gaps in temperature data gracefully + + Parameters + ---------- + met_start, met_stop : float + Time range (MET) for the correction table. + temptable : astropy.table.Table or str, optional + Temperature table or path to temperature file. If None, uses default. + + Other parameters + ---------------- + freqchange_file : str, optional + Path to frequency changes file. If None, uses default. + hdf_dump_file : str, optional + Path for caching computed corrections. Set to None to disable caching. + force_divisor : int, optional + If set, use this divisor value for entire range instead of + reading from frequency changes file. + time_resolution : float, optional + Output table time step in seconds. Default 0.5s. + craig_fit : bool, optional + Use original Craig Markwardt model. + version : str, optional + Model version identifier. + pars : dict, optional + Override model parameters. + + Returns + ------- + table : astropy.table.Table + Correction table with columns: + - 'met': Mission Elapsed Time + - 'temp_corr': Clock correction in seconds (add to event times) + - 'divisor': Clock divisor value at each time + """ import six if hdf_dump_file is not None and os.path.exists(hdf_dump_file): log.info(f"Reading cached data from file {hdf_dump_file}") From 350499e24aa0fed691f3343e98f91246e0b8a413 Mon Sep 17 00:00:00 2001 From: Matteo Bachetti Date: Thu, 2 Apr 2026 11:48:06 +0200 Subject: [PATCH 20/34] Add comments about added table columns --- nuclockutils/nustarclock.py | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/nuclockutils/nustarclock.py b/nuclockutils/nustarclock.py index d48b0dc..8b1eec3 100644 --- a/nuclockutils/nustarclock.py +++ b/nuclockutils/nustarclock.py @@ -319,6 +319,8 @@ def load_and_flag_clock_table(clockfile="latest_clock.dat", shift_non_malindi=Fa """ clock_offset_table = load_clock_offset_table(clockfile, shift_non_malindi=shift_non_malindi) + + # Column added/updated: 'flag' (bool) - True for known bad measurements based on db_file clock_offset_table = flag_bad_points( clock_offset_table, db_file=db_file) return clock_offset_table @@ -627,6 +629,7 @@ def eliminate_trends_in_residuals(temptable, clock_offset_table, log.info("Final detrending...") + # Columns added to temptable: 'std', 'temp_corr_trend', 'temp_corr_detrend' table_new = spline_detrending( clock_offset_table, temptable, outlier_cuts=[-0.002, -0.001], @@ -1235,8 +1238,10 @@ def load_temptable(temptable_name): if IS_CSV and os.path.exists(hdf5_name): IS_CSV = False + # Returns table with: 'met', 'temperature', 'mjd', 'temperature_smooth' temptable_raw = read_temptable(hdf5_name, dt=10) else: + # Returns table with: 'met', 'temperature', 'mjd', 'temperature_smooth' temptable_raw = read_temptable(temptable_name, dt=10) if IS_CSV: @@ -1247,12 +1252,14 @@ def load_temptable(temptable_name): @lru_cache(maxsize=64) def load_freq_changes(freq_change_file): + # Returns table with: 'uxt', 'met', 'divisor', 'mjd', 'flag' log.info(f"Reading data from {freq_change_file}") return read_freq_changes_table(freq_change_file) @lru_cache(maxsize=64) def load_clock_offset_table(clock_offset_file, shift_non_malindi=False): + # Returns table with: 'met', 'offset', 'station', 'mjd', 'flag' return read_clock_offset_table(clock_offset_file, shift_non_malindi=shift_non_malindi) @@ -1319,6 +1326,7 @@ def __init__(self, temperature_file, mjdstart=None, mjdstop=None, self.mjdstart = mjdstart self.mjdstop = mjdstop + # Sets self.temptable with: 'met', 'temperature', 'mjd', 'temperature_smooth' self.read_temptable() if mjdstart is None: @@ -1327,6 +1335,7 @@ def __init__(self, temperature_file, mjdstart=None, mjdstop=None, mjdstart = mjdstart - additional_days / 2 self.clock_offset_file = clock_offset_file + # Returns table with: 'met', 'offset', 'station', 'mjd', 'flag' self.clock_offset_table = \ read_clock_offset_table(self.clock_offset_file, shift_non_malindi=True) @@ -1372,6 +1381,7 @@ def __init__(self, temperature_file, mjdstart=None, mjdstop=None, if adjust_absolute_timing: log.info("Adjusting temperature correction") try: + # Adds 'temp_corr_nodetrend'; replaces 'temp_corr' with detrended version self.adjust_temperature_correction() except Exception: logfile = 'adjust_temperature_error.log' @@ -1386,6 +1396,7 @@ def read_temptable(self, cache_temptable_name=None): os.path.exists(cache_temptable_name): self.temptable = Table.read(cache_temptable_name, path='temptable') else: + # Returns table with: 'met', 'temperature', 'mjd', 'temperature_smooth' self.temptable = \ read_temptable(temperature_file=self.temperature_file, mjdstart=self.mjdstart, @@ -1411,6 +1422,7 @@ def temperature_correction_fun(self, adjust=False): def adjust_temperature_correction(self): """Adjust the temperature correction using measured clock offsets.""" + # Adds 'temp_corr_detrend' column to returned table table_new = eliminate_trends_in_residuals( self.temperature_correction_data, load_clock_offset_table( @@ -2347,10 +2359,12 @@ def temperature_correction_table(met_start, met_stop, if temptable is None or isinstance(temptable, six.string_types): mjdstart, mjdstop = sec_to_mjd(met_start), sec_to_mjd(met_stop) + # Returns table with: 'met', 'temperature', 'mjd', 'temperature_smooth' temptable = read_temptable(mjdstart=mjdstart, mjdstop=mjdstop, temperature_file=temptable, dt=10) if force_divisor is None: + # Returns table with: 'uxt', 'met', 'divisor', 'mjd', 'flag' freq_changes_table = \ read_freq_changes_table(freqchange_file=freqchange_file) allfreqtimes = np.array(freq_changes_table['met']) From 1a327d8d7d6162f2f5ce73955e9287f9fa69c5f4 Mon Sep 17 00:00:00 2001 From: Matteo Bachetti Date: Thu, 2 Apr 2026 12:21:08 +0200 Subject: [PATCH 21/34] Use helper function for filtering and logging --- nuclockutils/nustarclock.py | 176 ++++++++++++++++++------------------ 1 file changed, 87 insertions(+), 89 deletions(-) diff --git a/nuclockutils/nustarclock.py b/nuclockutils/nustarclock.py index 8b1eec3..659b7d4 100644 --- a/nuclockutils/nustarclock.py +++ b/nuclockutils/nustarclock.py @@ -144,6 +144,24 @@ datadir = os.path.join(curdir, 'data') +def filter_and_log_table( + table, + mask, + intro_text="Filtering table", + comment_to_point=None, + point_log_func=log.info + ): + n_before = mask.size + n_after = np.count_nonzero(mask) + comment_to_point = f" ({comment_to_point})" if comment_to_point else "" + + if n_after < n_before: + log.info(f"{intro_text} ({n_after} -> {n_before})") + for eliminated in table[~mask]['met']: + point_log_func(f" MET {eliminated}{comment_to_point}") + return table[mask] + return table + def get_bad_points_db(db_file=None): """Load the database of known bad clock offset measurement times. @@ -363,15 +381,12 @@ def spline_detrending(clock_offset_table, temptable, outlier_cuts=None, temperature_is_present = tempcorr_idx < temptable['met'].size tempcorr_idx = tempcorr_idx[temperature_is_present] - n_before = temperature_is_present.size - n_after = np.count_nonzero(temperature_is_present) - - if n_after < n_before: - log.info(f"spline_detrending: Filtering clock_offset_table:") - for eliminated in clock_offset_table[~temperature_is_present]['met']: - log.info(f" MET {eliminated} (beyond temperature table range)") - - clock_offset_table = clock_offset_table[temperature_is_present] + clock_offset_table = filter_and_log_table( + clock_offset_table, + temperature_is_present, + intro_text="spline_detrending: Filtering clock_offset_table to times with temperature data", + comment_to_point="beyond temperature table range" + ) clock_residuals = \ np.array(clock_offset_table['offset'] - @@ -394,16 +409,12 @@ def spline_detrending(clock_offset_table, temptable, outlier_cuts=None, do_not_flag = clock_mets > clock_mets.max() - SECONDS_PER_MONTH better_points[do_not_flag] = True - n_before = better_points.size - n_after = np.count_nonzero(better_points) - - if n_after < n_before: - log.info(f"spline_detrending: Filtering clock_offset_table outliers") - for eliminated in clock_offset_table[~better_points]['met']: - log.info(f" MET {eliminated}") - - clock_offset_table = clock_offset_table[better_points] - clock_residuals = clock_residuals[better_points] + clock_offset_table = filter_and_log_table( + clock_offset_table, + better_points, + "spline_detrending: Filtering clock_offset_table outliers" + ) + clock_residuals = clock_residuals[better_points] detrend_fun = spline_through_data( clock_offset_table['met'], clock_residuals, downsample=10, @@ -475,15 +486,12 @@ def eliminate_trends_in_residuals(temptable, clock_offset_table, temperature_is_present = tempcorr_idx < temptable['met'].size tempcorr_idx = tempcorr_idx[temperature_is_present] - n_before = temperature_is_present.size - n_after = np.count_nonzero(temperature_is_present) - - if n_after < n_before: - log.info(f"eliminate_trends_in_residuals: Filtering clock_offset_table") - for eliminated in clock_offset_table[~temperature_is_present]['met']: - log.info(f" MET {eliminated} (beyond temperature table range)") - - clock_offset_table = clock_offset_table[temperature_is_present] + clock_offset_table = filter_and_log_table( + clock_offset_table, + temperature_is_present, + intro_text="eliminate_trends_in_residuals: Filtering clock_offset_table to times with temperature data", + comment_to_point="beyond temperature table range" + ) clock_residuals = \ clock_offset_table['offset'] - temptable['temp_corr'][tempcorr_idx] @@ -499,17 +507,15 @@ def eliminate_trends_in_residuals(temptable, clock_offset_table, good = (clock_residuals == clock_residuals) & ~clock_offset_table['flag'] & use_for_interpol - n_before = good.size - n_after = np.count_nonzero(good) - - if n_after < n_before: - log.info("eliminate_trends_in_residuals: Filtering clock_offset_table to " - "trustworthy points (Malindi or during Malindi outage)") - for eliminated in clock_offset_table[~good]['met']: - log.debug(f" MET {eliminated} (bad point or not from Malindi)") - - clock_offset_table = clock_offset_table[good] - clock_residuals = clock_residuals[good] + clock_offset_table = filter_and_log_table( + clock_offset_table, + good, + intro_text="eliminate_trends_in_residuals: Filtering clock_offset_table to " + "trustworthy points (Malindi or during Malindi outage)", + comment_to_point="bad point or not from Malindi", + point_log_func=log.debug + ) + clock_residuals = clock_residuals[good] for g in gtis: log.info(f"Treating data from METs {g[0]}--{g[1]}") @@ -960,15 +966,13 @@ def read_freq_changes_table(freqchange_file=None, filter_bad=True): np.abs(freq_changes_table['divisor'] - 2.400034e7) > 20 if filter_bad: - n_before = freq_changes_table['flag'].size - n_after = n_before - np.count_nonzero(freq_changes_table['flag']) - - if n_after < n_before: - log.info(f"read_freq_changes_table: Filtering freq_changes_table") - for eliminated in freq_changes_table['met'][freq_changes_table['flag']]: - log.debug(f" MET {eliminated} (bad frequency change point)") - - freq_changes_table = freq_changes_table[~freq_changes_table['flag']] + freq_changes_table = filter_and_log_table( + freq_changes_table, + ~freq_changes_table['flag'], + intro_text="read_freq_changes_table: Filtering freq_changes_table to remove bad points", + comment_to_point="bad frequency change point", + point_log_func=log.debug + ) return freq_changes_table @@ -1209,13 +1213,13 @@ def read_temptable(temperature_file=None, mjdstart=None, mjdstop=None, good = np.diff(temptable['met']) > 0 good = np.concatenate((good, [True])) - n_before = good.size - n_after = np.count_nonzero(good) - if n_after < n_before: - log.info(f"read_temptable: Filtering temptable for bad time points") - for eliminated in temptable['met'][~good]: - log.debug(f" MET {eliminated} (non increasing time)") - temptable = temptable[good] + temptable = filter_and_log_table( + temptable, + good, + intro_text="read_temptable: Filtering temptable for non-increasing time points", + comment_to_point="non-increasing time point", + point_log_func=log.debug + ) temptable.meta['gti'] = temp_gtis window = np.median(1000 / np.diff(temptable['met'])) @@ -1466,32 +1470,27 @@ def write_clock_file(self, filename=None, save_nodetrend=False, clock_offset_table = self.clock_offset_table good_clock = clock_offset_table['met'] < table_new['met'][-1] - n_before = good_clock.size - n_after = np.count_nonzero(good_clock) - if n_after < n_before: - log.info(f"write_clock_file: Filtering clock_offset_table:") - for eliminated in clock_offset_table[~good_clock]['met']: - log.info(f" MET {eliminated} " - f"(beyond table range)") - - clock_offset_table = clock_offset_table[good_clock] + clock_offset_table = filter_and_log_table( + clock_offset_table, + good_clock, + intro_text="write_clock_file: Filtering clock_offset_table to remove points" + " beyond the temperature table range", + comment_to_point="beyound temperature range", + ) tempcorr_idx = np.searchsorted(table_new['met'], clock_offset_table['met']) temperature_is_present = tempcorr_idx < table_new['met'].size tempcorr_idx = tempcorr_idx[temperature_is_present] - n_before = temperature_is_present.size - n_after = np.count_nonzero(temperature_is_present) - - if n_after < n_before: - log.info(f"write_clock_file: Filtering clock_offset_table:") - for eliminated in clock_offset_table[~temperature_is_present]['met']: - log.info(f" MET {eliminated} " - f"(beyond temperature range)") - - clock_offset_table = clock_offset_table[temperature_is_present] + clock_offset_table = filter_and_log_table( + clock_offset_table, + temperature_is_present, + intro_text="write_clock_file: Filtering clock_offset_table to remove points" + " without temperature information", + comment_to_point="no temperature data present", + ) clock_residuals_detrend = clock_offset_table['offset'] - \ table_new['temp_corr'][tempcorr_idx] @@ -1797,15 +1796,13 @@ def plot_scatter(new_clock_table, clock_offset_table, shift_times=0, plt.show() yint = - yint - n_before = good_mets.size - n_after = np.count_nonzero(good_mets) - if n_after < n_before: - log.info("plot_scatter: Filtering clock_offset_table:") - for eliminated in clock_offset_table[~good_mets]['met']: - log.info(f" MET {eliminated} " - f"(beyond interpolation range)") - clock_offset_table = clock_offset_table[good_mets] + clock_offset_table = filter_and_log_table( + clock_offset_table, + good_mets, + intro_text="plot_scatter: Filtering clock_offset_table to remove points beyond interpolation range", + comment_to_point="beyond interpolation range", + ) clock_offset_table['offset'] -= shift_times clock_mets = clock_offset_table['met'] @@ -2348,13 +2345,14 @@ def temperature_correction_table(met_start, met_stop, "Interval not fully included in cached data. Recalculating.") else: good = (mets >= met_start - 86400) & (mets < met_stop + 86400) - n_before = good.size - n_after = np.count_nonzero(good) - if n_after < n_before: - log.info(f"temperature_correction_table: Filtering result_table:") - for eliminated in result_table[~good]['met']: - log.info(f" MET {eliminated} (beyond requested range)") - filtered_table = result_table[good] + + filtered_table = filter_and_log_table( + result_table, + good, + intro_text="temperature_correction_table: Filtering cached result_table to requested MET range", + comment_to_point="beyond requested range", + point_log_func=log.debug + ) return filtered_table if temptable is None or isinstance(temptable, six.string_types): From 98a49267e995e58c39c681d2e65b13317254c0f1 Mon Sep 17 00:00:00 2001 From: Matteo Bachetti Date: Thu, 2 Apr 2026 14:27:31 +0200 Subject: [PATCH 22/34] Use named variables instead of magic numbers, throughout --- nuclockutils/__init__.py | 6 +++++ nuclockutils/barycorr.py | 8 +++---- nuclockutils/clean_clock/app.py | 4 ++-- nuclockutils/diagnostics/bary_and_fold_all.py | 6 ++--- .../diagnostics/compare_clock_files.py | 6 ++--- nuclockutils/diagnostics/compare_pulses.py | 3 ++- nuclockutils/diagnostics/fold_to_ephemeris.py | 11 +++++---- .../diagnostics/get_crab_ephemeris.py | 3 ++- nuclockutils/nustarclock.py | 23 ++++++++----------- nuclockutils/utils.py | 9 ++++---- 10 files changed, 42 insertions(+), 37 deletions(-) diff --git a/nuclockutils/__init__.py b/nuclockutils/__init__.py index 5e4a97e..1ffb7bb 100644 --- a/nuclockutils/__init__.py +++ b/nuclockutils/__init__.py @@ -22,3 +22,9 @@ version = '0.0.0' __version__ = version + +# Time constants (all in seconds unless noted) +SECONDS_PER_DAY = 86400 +HALF_DAY_SECONDS = 43200 # Minimum GTI duration for valid processing +SECONDS_PER_MONTH = SECONDS_PER_DAY * 30 +SECONDS_PER_YEAR = SECONDS_PER_DAY * 365.25 diff --git a/nuclockutils/barycorr.py b/nuclockutils/barycorr.py index 7c90279..c26aeaa 100644 --- a/nuclockutils/barycorr.py +++ b/nuclockutils/barycorr.py @@ -14,7 +14,7 @@ import nuclockutils from .utils import filter_with_region, high_precision_keyword_read from .nustarclock import interpolate_clock_function - +from . import SECONDS_PER_DAY class OrbitalFunctions(): lat_fun = None @@ -29,7 +29,7 @@ def get_orbital_functions(orbfile): orbtable = Table.read(orbfile) mjdref = high_precision_keyword_read(orbtable.meta, 'MJDREF') - times = Time(np.array(orbtable['TIME'] / 86400 + mjdref), format='mjd') + times = Time(np.array(orbtable['TIME'] / SECONDS_PER_DAY + mjdref), format='mjd') if 'GEODETIC' in orbtable.colnames: geod = np.array(orbtable['GEODETIC']) lat, lon, alt = geod[:, 0] * u.deg, geod[:, 1] * u.deg, geod[:, @@ -82,8 +82,8 @@ def get_barycentric_correction(orbfile, parfile, dt=5, ephem='DE421'): mjdref = high_precision_keyword_read(hdul[1].header, 'MJDREF') knots = no.X.get_knots() - mjds = np.arange(knots[1], knots[-2], dt / 86400) - mets = (mjds - mjdref) * 86400 + mjds = np.arange(knots[1], knots[-2], dt / SECONDS_PER_DAY) + mets = (mjds - mjdref) * SECONDS_PER_DAY obs, scale = 'nustar', "tt" toalist = [None] * len(mjds) diff --git a/nuclockutils/clean_clock/app.py b/nuclockutils/clean_clock/app.py index c84e202..ad1e73e 100644 --- a/nuclockutils/clean_clock/app.py +++ b/nuclockutils/clean_clock/app.py @@ -19,7 +19,7 @@ from nuclockutils.nustarclock import load_temptable, load_freq_changes, \ load_and_flag_clock_table, find_good_time_intervals, calculate_stats, \ - eliminate_trends_in_residuals, _BAD_POINTS_FILE + eliminate_trends_in_residuals, _BAD_POINTS_FILE, FIXED_CONTROL_POINTS import dash from dash import dcc @@ -80,7 +80,7 @@ def recalc(outfile='save_all.pickle'): table_new = eliminate_trends_in_residuals( table_new, clock_offset_table_corr, gtis, - fixed_control_points=np.arange(291e6, 295e6, 86400)) + fixed_control_points=FIXED_CONTROL_POINTS) mets = np.array(table_new['met']) start = mets[0] diff --git a/nuclockutils/diagnostics/bary_and_fold_all.py b/nuclockutils/diagnostics/bary_and_fold_all.py index 25db873..41fe73a 100644 --- a/nuclockutils/diagnostics/bary_and_fold_all.py +++ b/nuclockutils/diagnostics/bary_and_fold_all.py @@ -17,7 +17,7 @@ get_ephemeris_from_parfile from .compare_pulses import main as main_compare_pulses from nuclockutils.barycorr import main_barycorr as nubarycorr -from nuclockutils.utils import high_precision_keyword_read +from nuclockutils.utils import high_precision_keyword_read, SECONDS_PER_DAY def mkdir(folder): @@ -40,7 +40,7 @@ def get_observing_mjd(fname): log.error(f"{fname} might be corrupted") return None - return np.array([tstart, tstop]) / 86400 + mjdref + return np.array([tstart, tstop]) / SECONDS_PER_DAY + mjdref def fold_file_to_ephemeris(fname, parfile, emin=None, emax=None, @@ -61,7 +61,7 @@ def fold_file_to_ephemeris(fname, parfile, emin=None, emax=None, good = good & (events.energy < emax) times = events.time[good] - event_mjds = times / 86400 + events.mjdref + event_mjds = times / SECONDS_PER_DAY + events.mjdref phase = correction_fun(event_mjds) t = calculate_profile_from_phase(phase, nbin=nbin) diff --git a/nuclockutils/diagnostics/compare_clock_files.py b/nuclockutils/diagnostics/compare_clock_files.py index 8125293..3329ce8 100644 --- a/nuclockutils/diagnostics/compare_clock_files.py +++ b/nuclockutils/diagnostics/compare_clock_files.py @@ -3,17 +3,17 @@ import numpy as np from astropy.table import Table import matplotlib.pyplot as plt -from nuclockutils.nustarclock import interpolate_clock_function +from nuclockutils.nustarclock import interpolate_clock_function, SECONDS_PER_DAY def get_todays_met(): t = Time.now() - met = (t.mjd - 55197.00076601852) * 86400 + met = (t.mjd - 55197.00076601852) * SECONDS_PER_DAY return met def get_launch_met(): t = Time('2012-06-12') - met = (t.mjd - 55197.00076601852) * 86400 + met = (t.mjd - 55197.00076601852) * SECONDS_PER_DAY return met diff --git a/nuclockutils/diagnostics/compare_pulses.py b/nuclockutils/diagnostics/compare_pulses.py index 86e914d..44fc9a1 100644 --- a/nuclockutils/diagnostics/compare_pulses.py +++ b/nuclockutils/diagnostics/compare_pulses.py @@ -7,6 +7,7 @@ from uncertainties import ufloat from astropy import log from statsmodels.robust import mad +from nuclockutils.utils import SECONDS_PER_DAY from .fftfit import fftfit @@ -33,7 +34,7 @@ def format_profile_and_get_phase(file, template=None): fddot = 0 if 'F2' in table.meta: fddot = table.meta['F2'] - delta_t = (mjd - pepoch) * 86400 + delta_t = (mjd - pepoch) * SECONDS_PER_DAY f0 = freq + fdot * delta_t + 0.5 * fddot * delta_t ** 2 period_ms = 1 / f0 * 1000 diff --git a/nuclockutils/diagnostics/fold_to_ephemeris.py b/nuclockutils/diagnostics/fold_to_ephemeris.py index 54a938f..0411c85 100644 --- a/nuclockutils/diagnostics/fold_to_ephemeris.py +++ b/nuclockutils/diagnostics/fold_to_ephemeris.py @@ -16,6 +16,7 @@ from pint import models from stingray.pulse.pulsar import _load_and_prepare_TOAs, get_model from scipy.interpolate import PchipInterpolator +from nuclockutils import SECONDS_PER_DAY @@ -147,8 +148,8 @@ def get_exposure_per_bin(event_times, event_priors, parfile, sampling_time = 1 / m.F0.value / nbin / 7.1234514351515132414 chunk_size = np.min((sampling_time * 50000000, 1000)) - mjdstart = (event_times[0] - 10) / 86400 + mjdref - mjdstop = (event_times[-1] + 10) / 86400 + mjdref + mjdstart = (event_times[0] - 10) / SECONDS_PER_DAY + mjdref + mjdstop = (event_times[-1] + 10) / SECONDS_PER_DAY + mjdref phase_correction_fun = get_phase_from_ephemeris_file(mjdstart, mjdstop, parfile) @@ -163,7 +164,7 @@ def get_exposure_per_bin(event_times, event_priors, parfile, dead = idxs_live == idxs_dead - sample_times = sample_times / 86400 + mjdref + sample_times = sample_times / SECONDS_PER_DAY + mjdref # phases_all = _fast_phase_fddot(sample_times, freq, fdot, fddot) # phases_dead = _fast_phase_fddot(sample_times[dead], freq, fdot, fddot) phases_all = phase_correction_fun(sample_times) @@ -272,7 +273,7 @@ def main(args=None): events = get_events_from_fits(evfile) MJDREF = events.mjdref - MJD = events.time.mean() / 86400 + MJDREF + MJD = events.time.mean() / SECONDS_PER_DAY + MJDREF good = np.ones(events.time.size, dtype=bool) if emin is not None: @@ -280,7 +281,7 @@ def main(args=None): if emax is not None: good = good & (events.energy < emax) - events_time = events.time[good] / 86400 + MJDREF + events_time = events.time[good] / SECONDS_PER_DAY + MJDREF mjdstart = events_time[0] - 0.01 mjdstop = events_time[-1] + 0.01 diff --git a/nuclockutils/diagnostics/get_crab_ephemeris.py b/nuclockutils/diagnostics/get_crab_ephemeris.py index 86a9e2d..5bde8ba 100644 --- a/nuclockutils/diagnostics/get_crab_ephemeris.py +++ b/nuclockutils/diagnostics/get_crab_ephemeris.py @@ -7,6 +7,7 @@ from astropy.table import Table from hendrics.io import high_precision_keyword_read +from nuclockutils import SECONDS_PER_DAY from urllib.request import urlopen @@ -85,7 +86,7 @@ def main(args=None): t0 = high_precision_keyword_read(header, 'TSTART') t1 = high_precision_keyword_read(header, 'TSTOP') mjdref = high_precision_keyword_read(header, 'MJDREF') - MJD = (t1 + t0) / 2 / 86400 + mjdref + MJD = (t1 + t0) / 2 / SECONDS_PER_DAY + mjdref get_crab_ephemeris( MJD, fname=os.path.splitext(os.path.basename(evfile))[0] + '.par') diff --git a/nuclockutils/nustarclock.py b/nuclockutils/nustarclock.py index 659b7d4..ec84647 100644 --- a/nuclockutils/nustarclock.py +++ b/nuclockutils/nustarclock.py @@ -99,6 +99,7 @@ from holoviews.operation.datashader import datashade from holoviews import opts +from . import SECONDS_PER_DAY, SECONDS_PER_MONTH, SECONDS_PER_YEAR, HALF_DAY_SECONDS # ============================================================================= # Constants # ============================================================================= @@ -106,12 +107,6 @@ # File paths _BAD_POINTS_FILE = "BAD_POINTS_DB.dat" -# Time constants (all in seconds unless noted) -SECONDS_PER_DAY = 86400 -HALF_DAY_SECONDS = 43200 # Minimum GTI duration for valid processing -SECONDS_PER_MONTH = SECONDS_PER_DAY * 30 -SECONDS_PER_YEAR = SECONDS_PER_DAY * 365.25 - # Station timing offset: non-Malindi stations have a systematic 0.5 ms offset # due to different signal processing pipelines NON_MALINDI_OFFSET_SECONDS = 0.0005 @@ -281,7 +276,7 @@ def find_good_time_intervals(temperature_table, clock_jump_times=None): gtis = cross_two_gtis(temp_gtis, clock_gtis) lengths = gtis[:, 1] - gtis[:, 0] # ensure at least half a day duration for GTIs - good = lengths > 43200 + good = lengths > HALF_DAY_SECONDS if not np.all(good): log.info(f"Some GTIs are too short. cleaning up: {gtis[~good]}") for g, is_good in zip(gtis, good): @@ -1083,7 +1078,7 @@ def read_csv_temptable(mjdstart=None, mjdstop=None, temperature_file=None): in_subfmt="date_hms").mjd log.info("Done.") temptable["mjd"] = np.array(times_mjd) - temptable['met'] = (temptable["mjd"] - NUSTAR_MJDREF) * 86400 + temptable['met'] = (temptable["mjd"] - NUSTAR_MJDREF) * SECONDS_PER_DAY temptable.remove_column('Time') temptable.rename_column('tp_eps_ceu_txco_tmp', 'temperature') temptable["temperature"] = np.array(temptable["temperature"], dtype=float) @@ -1353,8 +1348,8 @@ def __init__(self, temperature_file, mjdstart=None, mjdstop=None, self.mjdstart = mjdstart self.mjdstop = mjdstop - self.met_start = (self.mjdstart - NUSTAR_MJDREF) * 86400 - self.met_stop = (self.mjdstop - NUSTAR_MJDREF) * 86400 + self.met_start = (self.mjdstart - NUSTAR_MJDREF) * SECONDS_PER_DAY + self.met_stop = (self.mjdstop - NUSTAR_MJDREF) * SECONDS_PER_DAY if label is None or label == "": label = f"{self.met_start}-{self.met_stop}" @@ -1532,7 +1527,7 @@ def write_clock_file(self, filename=None, save_nodetrend=False, if not highres: good_for_clockfile = np.zeros(allmets.size, dtype=bool) good_for_clockfile[::100] = True - twodays = 86400 * 2 + twodays = SECONDS_PER_DAY * 2 for jumptime in self.clock_jump_times: idx0, idx1 = np.searchsorted( allmets, [jumptime - twodays, jumptime + twodays]) @@ -2340,11 +2335,11 @@ def temperature_correction_table(met_start, met_stop, result_table = fix_byteorder(Table.read(hdf_dump_file)) mets = np.array(result_table['met']) if (met_start > mets[10] or met_stop < mets[-20]) and ( - met_stop - met_start < 3 * 365 * 86400): + met_stop - met_start < 3 * SECONDS_PER_YEAR): log.warning( "Interval not fully included in cached data. Recalculating.") else: - good = (mets >= met_start - 86400) & (mets < met_stop + 86400) + good = (mets >= met_start - SECONDS_PER_DAY * 2) & (mets < met_stop + SECONDS_PER_DAY * 2) filtered_table = filter_and_log_table( result_table, @@ -2563,7 +2558,7 @@ def main_update_temptable(args=None): if args.outfile is not None and os.path.exists(args.outfile): log.info("Reading existing temperature table") existing_table = Table.read(args.outfile) - last_measurement = existing_table['mjd'][-1] + 0.001 / 86400 + last_measurement = existing_table['mjd'][-1] + 0.001 / SECONDS_PER_DAY log.info("Reading new temperature values") new_table = read_csv_temptable(temperature_file=args.tempfile, diff --git a/nuclockutils/utils.py b/nuclockutils/utils.py index 2fef441..07f5bca 100644 --- a/nuclockutils/utils.py +++ b/nuclockutils/utils.py @@ -19,6 +19,7 @@ except ImportError: pass +from . import SECONDS_PER_DAY NUSTAR_MJDREF = np.longdouble("55197.00076601852") @@ -60,11 +61,11 @@ def fix_byteorder(table): def sec_to_mjd(time, mjdref=NUSTAR_MJDREF, dtype=np.double): - return np.array(np.asarray(time) / 86400 + mjdref, dtype=dtype) + return np.array(np.asarray(time) / SECONDS_PER_DAY + mjdref, dtype=dtype) def mjd_to_sec(mjd, mjdref=NUSTAR_MJDREF, dtype=np.double): - return np.array((np.asarray(mjd) - mjdref) * 86400, dtype=dtype) + return np.array((np.asarray(mjd) - mjdref) * SECONDS_PER_DAY, dtype=dtype) def sec_to_ut(time, mjdref=NUSTAR_MJDREF, dtype=np.double): @@ -73,7 +74,7 @@ def sec_to_ut(time, mjdref=NUSTAR_MJDREF, dtype=np.double): def ut_to_sec(isot, mjdref=NUSTAR_MJDREF, dtype=np.double): time = Time(isot) - return np.array((np.asarray(time.mjd) - mjdref) * 86400, dtype=dtype) + return np.array((np.asarray(time.mjd) - mjdref) * SECONDS_PER_DAY, dtype=dtype) def splitext_improved(path): @@ -215,7 +216,7 @@ def get_obsid_list_from_heasarc(cache_file='heasarc.hdf5'): mjd_ends = Time(np.array(all_nustar_obs['end_time']), format='mjd') # all_nustar_obs = all_nustar_obs[all_nustar_obs["observation_mode"] == 'SCIENCE'] - all_nustar_obs['met'] = np.array(all_nustar_obs['time'] - NUSTAR_MJDREF) * 86400 + all_nustar_obs['met'] = np.array(all_nustar_obs['time'] - NUSTAR_MJDREF) * SECONDS_PER_DAY all_nustar_obs['date'] = mjds.fits all_nustar_obs['date-end'] = mjd_ends.fits From cdbd2dd4e49f7f192b7442ecfb497de3e65df8b2 Mon Sep 17 00:00:00 2001 From: Matteo Bachetti Date: Thu, 2 Apr 2026 15:29:19 +0200 Subject: [PATCH 23/34] Fix __init__ import --- nuclockutils/__init__.py | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/nuclockutils/__init__.py b/nuclockutils/__init__.py index 1ffb7bb..9cb7725 100644 --- a/nuclockutils/__init__.py +++ b/nuclockutils/__init__.py @@ -1,5 +1,14 @@ # Licensed under a 3-clause BSD style license - see LICENSE.rst +# Time constants (all in seconds unless noted) +SECONDS_PER_DAY = 86400 +HALF_DAY_SECONDS = 43200 # Minimum GTI duration for valid processing +SECONDS_PER_MONTH = SECONDS_PER_DAY * 30 +SECONDS_PER_YEAR = SECONDS_PER_DAY * 365.25 + +# Mission timeline constants (MET values) +# MET when science operations began (end of commissioning phase) +SCIENCE_START_MET = 77674700 from .nustarclock import * from .utils import * @@ -22,9 +31,3 @@ version = '0.0.0' __version__ = version - -# Time constants (all in seconds unless noted) -SECONDS_PER_DAY = 86400 -HALF_DAY_SECONDS = 43200 # Minimum GTI duration for valid processing -SECONDS_PER_MONTH = SECONDS_PER_DAY * 30 -SECONDS_PER_YEAR = SECONDS_PER_DAY * 365.25 From 7f2c351ecf1028f010029b5acf1aa338652aab10 Mon Sep 17 00:00:00 2001 From: Hannah Date: Fri, 3 Apr 2026 13:27:44 -0700 Subject: [PATCH 24/34] Catch and warn for invalid MET intervals --- nuclockutils/nustarclock.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/nuclockutils/nustarclock.py b/nuclockutils/nustarclock.py index ec84647..9a0b324 100644 --- a/nuclockutils/nustarclock.py +++ b/nuclockutils/nustarclock.py @@ -2392,6 +2392,9 @@ def temperature_correction_table(met_start, met_stop, continue if met_intv[0] > met_stop: break + if met_intv[0] > met_intv[1]: + log.warning(f"Invalid interval: {met_intv}") + continue start, stop = met_intv From b6e54f8f03f53f0fb516e39aa230fd0ea50ec953 Mon Sep 17 00:00:00 2001 From: Matteo Bachetti Date: Wed, 8 Apr 2026 09:11:00 +0200 Subject: [PATCH 25/34] Re-read and sort temperature table when hdf5 is older --- nuclockutils/clean_clock/app.py | 1 + nuclockutils/nustarclock.py | 8 +++++++- 2 files changed, 8 insertions(+), 1 deletion(-) diff --git a/nuclockutils/clean_clock/app.py b/nuclockutils/clean_clock/app.py index ad1e73e..06a407f 100644 --- a/nuclockutils/clean_clock/app.py +++ b/nuclockutils/clean_clock/app.py @@ -519,6 +519,7 @@ def refresh_output(n_cl_refresh, n_cl_bad, n_cl_good, n_cl_recalc, selectedData) if who_triggered == 'recalculate-button': log.info("Recalculating all") + stored_analysis.cache_clear() # cache.delete_memoized(stored_analysis, 'save_all.pickle') if os.path.exists('save_all.pickle'): diff --git a/nuclockutils/nustarclock.py b/nuclockutils/nustarclock.py index 9a0b324..bee521d 100644 --- a/nuclockutils/nustarclock.py +++ b/nuclockutils/nustarclock.py @@ -1235,7 +1235,12 @@ def load_temptable(temptable_name): IS_CSV = temptable_name.endswith('.csv') hdf5_name = temptable_name.replace('.csv', '.hdf5') - if IS_CSV and os.path.exists(hdf5_name): + h5_exists = os.path.exists(hdf5_name) + h5_newer = h5_exists and os.path.getmtime(hdf5_name) > os.path.getmtime(temptable_name) + if h5_exists and not h5_newer: + log.info(f"HDF5 file {hdf5_name} is older than CSV file {temptable_name}. " + "Re-reading from CSV and overwriting HDF5.") + if IS_CSV and h5_newer: IS_CSV = False # Returns table with: 'met', 'temperature', 'mjd', 'temperature_smooth' temptable_raw = read_temptable(hdf5_name, dt=10) @@ -1243,6 +1248,7 @@ def load_temptable(temptable_name): # Returns table with: 'met', 'temperature', 'mjd', 'temperature_smooth' temptable_raw = read_temptable(temptable_name, dt=10) + temptable_raw.sort("met") if IS_CSV: log.info(f"Saving temperature data to {hdf5_name}") temptable_raw.write(hdf5_name, overwrite=True) From c192f2ca0dc9705a7e6228d6ed101305c18fe53b Mon Sep 17 00:00:00 2001 From: Matteo Bachetti Date: Wed, 8 Apr 2026 11:29:12 +0200 Subject: [PATCH 26/34] Fix condition for re-read of temperature --- nuclockutils/nustarclock.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/nuclockutils/nustarclock.py b/nuclockutils/nustarclock.py index bee521d..e977d97 100644 --- a/nuclockutils/nustarclock.py +++ b/nuclockutils/nustarclock.py @@ -1237,7 +1237,7 @@ def load_temptable(temptable_name): h5_exists = os.path.exists(hdf5_name) h5_newer = h5_exists and os.path.getmtime(hdf5_name) > os.path.getmtime(temptable_name) - if h5_exists and not h5_newer: + if IS_CSV and h5_exists and not h5_newer: log.info(f"HDF5 file {hdf5_name} is older than CSV file {temptable_name}. " "Re-reading from CSV and overwriting HDF5.") if IS_CSV and h5_newer: From f2fe2b359c10060f0c914d55fcf4f89be17841de Mon Sep 17 00:00:00 2001 From: Matteo Bachetti Date: Wed, 8 Apr 2026 11:41:54 +0200 Subject: [PATCH 27/34] Thoroughly eliminate cache files when recalculating --- nuclockutils/clean_clock/app.py | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/nuclockutils/clean_clock/app.py b/nuclockutils/clean_clock/app.py index 06a407f..6c221c5 100644 --- a/nuclockutils/clean_clock/app.py +++ b/nuclockutils/clean_clock/app.py @@ -520,10 +520,14 @@ def refresh_output(n_cl_refresh, n_cl_bad, n_cl_good, n_cl_recalc, selectedData) if who_triggered == 'recalculate-button': log.info("Recalculating all") + for file in ['save_all.pickle', 'dump.hdf5', 'all_data_res.pkl', 'rolling_data.pkl']: + + if os.path.exists(file): + log.info(f"Removing file {file}") + os.unlink(file) + stored_analysis.cache_clear() - # cache.delete_memoized(stored_analysis, 'save_all.pickle') - if os.path.exists('save_all.pickle'): - os.unlink('save_all.pickle') + elif who_triggered == 'actually-good-data-button' and len(NEW_BAD_POINTS) > 0: log.info("Removing point(s) from bad clock offset database") ALL_BAD_POINTS = eliminate_array_from_array( From 36e7d1f06ab2e0dceedff48786999efe9e057925 Mon Sep 17 00:00:00 2001 From: Matteo Bachetti Date: Wed, 8 Apr 2026 11:57:13 +0200 Subject: [PATCH 28/34] Make points more visible --- nuclockutils/clean_clock/app.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/nuclockutils/clean_clock/app.py b/nuclockutils/clean_clock/app.py index 6c221c5..45e7e5a 100644 --- a/nuclockutils/clean_clock/app.py +++ b/nuclockutils/clean_clock/app.py @@ -220,7 +220,7 @@ def plot_dash(all_data, table_new, gti, all_nustar_obs, 'text': text[bad], 'mode': 'markers', 'name': f'Bad clock offset measurements', - 'marker': {'color': 'grey', 'symbol': "x-dot", 'size': 3} + 'marker': {'color': 'grey', 'symbol': "x-dot", 'size': 5} }), 1, 1) all_data_bad = all_data[bad] @@ -233,7 +233,7 @@ def plot_dash(all_data, table_new, gti, all_nustar_obs, 'text': text[bad], 'mode': 'markers', 'showlegend': False, - 'marker': {'color': 'grey', 'symbol': "x-dot", 'size': 3} + 'marker': {'color': 'grey', 'symbol': "x-dot", 'size': 5} }), row, 1) for station, color in zip(['MLD', 'SNG', 'UHI'], ['blue', 'red', 'orange']): @@ -245,7 +245,7 @@ def plot_dash(all_data, table_new, gti, all_nustar_obs, 'hovertemplate': hovertemplate, 'text': text[good], 'mode': 'markers', - 'marker': {'color': color, 'size': 3}, + 'marker': {'color': color, 'size': 5}, 'name': f'Clock offset - {station}' }), 1, 1) for ydata, row in zip(['residual', 'residual_detrend'], [2, 3]): @@ -256,7 +256,7 @@ def plot_dash(all_data, table_new, gti, all_nustar_obs, 'text': text[good], 'showlegend': False, 'mode': 'markers', - 'marker': {'color': color, 'size': 3} + 'marker': {'color': color, 'size': 5} }), row, 1) # bad_intervals = [[0, 77.767e6]] From 5487f0025bdf145318de46b102362a1899bfb2c2 Mon Sep 17 00:00:00 2001 From: Matteo Bachetti Date: Wed, 8 Apr 2026 13:27:37 +0200 Subject: [PATCH 29/34] Improve margins of temperature plots --- nuclockutils/clean_clock/app.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/nuclockutils/clean_clock/app.py b/nuclockutils/clean_clock/app.py index 45e7e5a..ebb956b 100644 --- a/nuclockutils/clean_clock/app.py +++ b/nuclockutils/clean_clock/app.py @@ -563,11 +563,12 @@ def create_temperature_timeseries(x, y, axis_type='linear'): "data": [dict(x=x.astype(float), y=y.astype(float), mode="lines")], "layout": { "height": 300, + "margin": {"t": 20, "b": 40, "l": 50, "r": 20}, "yaxis": { "title": "TCXO Temperature", "type": "linear" if axis_type == "Linear" else "log", }, - "xaxis": {"title": "met", "showgrid": False, "margin": {"t": 20}}, + "xaxis": {"title": "met", "showgrid": False}, }, } @@ -577,8 +578,9 @@ def create_temperature_gradient_timeseries(x, y, axis_type='linear'): "data": [dict(x=x.astype(float), y=y.astype(float), mode="lines")], "layout": { "height": 300, + "margin": {"t": 20, "b": 40, "l": 50, "r": 20}, "yaxis": {"title": "TCXO Temp Gradient", "type": "linear"}, - "xaxis": {"title": "met", "showgrid": False, "margin": {"t": 20}}, + "xaxis": {"title": "met", "showgrid": False}, }, } From d4bb594edf5805e1b8469f738c963cb9085d1923 Mon Sep 17 00:00:00 2001 From: Matteo Bachetti Date: Wed, 8 Apr 2026 13:34:40 +0200 Subject: [PATCH 30/34] Fix axis labels of temperature plots --- nuclockutils/clean_clock/app.py | 44 +++++++++++++++++---------------- 1 file changed, 23 insertions(+), 21 deletions(-) diff --git a/nuclockutils/clean_clock/app.py b/nuclockutils/clean_clock/app.py index ebb956b..c9e2743 100644 --- a/nuclockutils/clean_clock/app.py +++ b/nuclockutils/clean_clock/app.py @@ -559,30 +559,32 @@ def display_selected_data(selectedData): def create_temperature_timeseries(x, y, axis_type='linear'): - return { - "data": [dict(x=x.astype(float), y=y.astype(float), mode="lines")], - "layout": { - "height": 300, - "margin": {"t": 20, "b": 40, "l": 50, "r": 20}, - "yaxis": { - "title": "TCXO Temperature", - "type": "linear" if axis_type == "Linear" else "log", - }, - "xaxis": {"title": "met", "showgrid": False}, - }, - } + import plotly.graph_objects as go + fig = go.Figure() + fig.add_trace(go.Scattergl(x=np.asarray(x, dtype=float), + y=np.asarray(y, dtype=float), mode="lines")) + fig.update_layout( + height=300, + margin=dict(t=20, b=50, l=60, r=20), + ) + fig.update_xaxes(title_text="MET (s)", showgrid=False) + fig.update_yaxes(title_text="TCXO Temperature (°C)", + type="linear" if axis_type == "linear" else "log") + return fig def create_temperature_gradient_timeseries(x, y, axis_type='linear'): - return { - "data": [dict(x=x.astype(float), y=y.astype(float), mode="lines")], - "layout": { - "height": 300, - "margin": {"t": 20, "b": 40, "l": 50, "r": 20}, - "yaxis": {"title": "TCXO Temp Gradient", "type": "linear"}, - "xaxis": {"title": "met", "showgrid": False}, - }, - } + import plotly.graph_objects as go + fig = go.Figure() + fig.add_trace(go.Scattergl(x=np.asarray(x, dtype=float), + y=np.asarray(y, dtype=float), mode="lines")) + fig.update_layout( + height=300, + margin=dict(t=20, b=50, l=60, r=20), + ) + fig.update_xaxes(title_text="MET (s)", showgrid=False) + fig.update_yaxes(title_text="TCXO Temp Gradient (°C/s)", type="linear") + return fig @app.callback( From fdcb06c23cdadddf6edce596b620c6e48111d7c9 Mon Sep 17 00:00:00 2001 From: Matteo Bachetti Date: Wed, 8 Apr 2026 15:29:51 +0200 Subject: [PATCH 31/34] Fix temperature labels --- nuclockutils/clean_clock/app.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/nuclockutils/clean_clock/app.py b/nuclockutils/clean_clock/app.py index c9e2743..7cec961 100644 --- a/nuclockutils/clean_clock/app.py +++ b/nuclockutils/clean_clock/app.py @@ -418,8 +418,8 @@ def stored_analysis(file): html.Div( className="three columns div-for-charts bg-grey", children=[ - dcc.Graph(id='temperature-time-series'), - dcc.Graph(id='temperature-gradient-time-series'), + dcc.Graph(mathjax=True, id='temperature-time-series'), + dcc.Graph(mathjax=True, id='temperature-gradient-time-series'), ] ), html.Div(id='dummy', style={'display': 'none'}), @@ -577,13 +577,13 @@ def create_temperature_gradient_timeseries(x, y, axis_type='linear'): import plotly.graph_objects as go fig = go.Figure() fig.add_trace(go.Scattergl(x=np.asarray(x, dtype=float), - y=np.asarray(y, dtype=float), mode="lines")) + y=np.asarray(y, dtype=float) * 1e3, mode="lines")) fig.update_layout( height=300, margin=dict(t=20, b=50, l=60, r=20), ) fig.update_xaxes(title_text="MET (s)", showgrid=False) - fig.update_yaxes(title_text="TCXO Temp Gradient (°C/s)", type="linear") + fig.update_yaxes(title_text=r"TCXO Temp Gradient (1e-3 °C/s)", type="linear") return fig From 1d2c3526a3942c20b7853c9b6fefbf873484cb5e Mon Sep 17 00:00:00 2001 From: Matteo Bachetti Date: Wed, 8 Apr 2026 23:00:52 +0200 Subject: [PATCH 32/34] Sort MET immediately --- nuclockutils/nustarclock.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/nuclockutils/nustarclock.py b/nuclockutils/nustarclock.py index e977d97..711a83e 100644 --- a/nuclockutils/nustarclock.py +++ b/nuclockutils/nustarclock.py @@ -1080,6 +1080,7 @@ def read_csv_temptable(mjdstart=None, mjdstop=None, temperature_file=None): temptable["mjd"] = np.array(times_mjd) temptable['met'] = (temptable["mjd"] - NUSTAR_MJDREF) * SECONDS_PER_DAY temptable.remove_column('Time') + temptable.sort("met") temptable.rename_column('tp_eps_ceu_txco_tmp', 'temperature') temptable["temperature"] = np.array(temptable["temperature"], dtype=float) if os.path.exists('tmp.csv'): @@ -1248,7 +1249,6 @@ def load_temptable(temptable_name): # Returns table with: 'met', 'temperature', 'mjd', 'temperature_smooth' temptable_raw = read_temptable(temptable_name, dt=10) - temptable_raw.sort("met") if IS_CSV: log.info(f"Saving temperature data to {hdf5_name}") temptable_raw.write(hdf5_name, overwrite=True) From 7e8587607b4c1d081c79f3d24d0af5efa7916f69 Mon Sep 17 00:00:00 2001 From: Matteo Bachetti Date: Mon, 13 Apr 2026 17:12:18 +0200 Subject: [PATCH 33/34] Test up to Python 3.14 --- .github/workflows/ci-test.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/ci-test.yml b/.github/workflows/ci-test.yml index 9f225b7..ff5be5a 100644 --- a/.github/workflows/ci-test.yml +++ b/.github/workflows/ci-test.yml @@ -16,7 +16,7 @@ jobs: strategy: fail-fast: false matrix: - python-version: ["3.9", "3.11", "3.13"] + python-version: ["3.9", "3.12", "3.14"] steps: - uses: actions/checkout@v4 From 09752f5303e40f0fc63788c7fab32f16f514b93d Mon Sep 17 00:00:00 2001 From: Murray Brightman Date: Mon, 20 Apr 2026 15:07:11 -0700 Subject: [PATCH 34/34] Updated bad samples for April 20, 2026 release --- nuclockutils/data/BAD_POINTS_DB.dat | 1 + 1 file changed, 1 insertion(+) diff --git a/nuclockutils/data/BAD_POINTS_DB.dat b/nuclockutils/data/BAD_POINTS_DB.dat index 8dfce4f..64ac16f 100644 --- a/nuclockutils/data/BAD_POINTS_DB.dat +++ b/nuclockutils/data/BAD_POINTS_DB.dat @@ -1694,3 +1694,4 @@ 512698376 512839329 513133512 +513704513