From dbdef6d25a54e932281f0edd7862120d1e19b9ed Mon Sep 17 00:00:00 2001 From: Matthew Kerr Date: Thu, 14 Nov 2024 10:45:49 -0500 Subject: [PATCH 01/10] Fix units on WAVE_OM. --- src/pint/models/wave.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/src/pint/models/wave.py b/src/pint/models/wave.py index 8a902a5ca..64e6da264 100644 --- a/src/pint/models/wave.py +++ b/src/pint/models/wave.py @@ -15,9 +15,8 @@ class Wave(PhaseComponent): sine/cosine components. For consistency with the implementation in tempo2, this signal is treated - as a time series, but trivially converted into phase by multiplication by - F0, which could makes changes to PEPOCH fragile if there is strong spin - frequency evolution. + as a time series and trivially converted into phase by multiplication by + F0. Parameters supported: @@ -42,7 +41,7 @@ def __init__(self): floatParameter( name="WAVE_OM", description="Base frequency of wave solution", - units="1/d", + units="rad/d", convert_tcb2tdb=False, ) ) From e5ce925376eeef08e3054047724bedf18321ae70 Mon Sep 17 00:00:00 2001 From: Matthew Kerr Date: Thu, 14 Nov 2024 10:46:08 -0500 Subject: [PATCH 02/10] Add explicit test for WAVE calculation. --- tests/test_model_wave.py | 85 ++++++++++++++++++++++++++++++++-------- 1 file changed, 69 insertions(+), 16 deletions(-) diff --git a/tests/test_model_wave.py b/tests/test_model_wave.py index 1e57d640f..ab1707088 100644 --- a/tests/test_model_wave.py +++ b/tests/test_model_wave.py @@ -1,28 +1,81 @@ +from io import StringIO import os import pytest import astropy.units as u +import numpy as np -import pint.fitter -import pint.models +from pint.models import get_model +from pint import toa import pint.residuals -import pint.toa from pinttestdata import datadir -# Not included in the test here, but as a sanity check I used this same -# ephemeris to phase up Fermi data, and it looks good. +par_nowave = """ + PSRJ J0835-4510 + RAJ 08:35:20.61149 + DECJ -45:10:34.8751 + F0 11.18965156782 + PEPOCH 55305 + DM 67.99 + UNITS TDB +""" -parfile = os.path.join(datadir, "vela_wave.par") -timfile = os.path.join(datadir, "vela_wave.tim") +wave_terms = """ + WAVEEPOCH 55305 + WAVE_OM 0.0015182579855022 + WAVE1 -0.21573979911255 -0.049063841960712 + WAVE2 0.62795320246729 -0.11984954655989 + WAVE3 0.099618608456845 0.28530756908788 + WAVE4 -0.21537340649058 0.18849486610628 + WAVE5 0.021980474493165 -0.23566696662127 +""" -class TestWave: - @classmethod - def setup_class(cls): - cls.m = pint.models.get_model(parfile) - cls.t = pint.toa.get_TOAs(timfile, ephem="DE405", include_bipm=False) +def test_wave_ephem(): + parfile = os.path.join(datadir, "vela_wave.par") + timfile = os.path.join(datadir, "vela_wave.tim") + m = get_model(parfile) + t = toa.get_TOAs(timfile, ephem="DE405", include_bipm=False) + rs = pint.residuals.Residuals(t, m).time_resids + assert rs.std() < 350.0 * u.us - def test_vela_rms_is_small_enough(self): - rs = pint.residuals.Residuals(self.t, self.m).time_resids - rms = rs.to(u.us).std() - assert rms < 350.0 * u.us + +def test_wave_construction(): + m = get_model(StringIO(par_nowave + wave_terms)) + assert np.allclose(m.WAVE_OM.quantity, 0.0015182579855022 * u.rad / u.day) + assert np.allclose(m.WAVE1.quantity[0], -0.21573979911255 * u.s) + assert np.allclose(m.WAVE2.quantity[1], -0.1198495465598 * u.s) + + +def test_wave_computation(): + m0 = get_model(StringIO(par_nowave)) + m1 = get_model(StringIO(par_nowave + wave_terms)) + # make some barycentric TOAs + tdbs = np.linspace(54500, 60000, 10) + ts = toa.TOAs( + toalist=[ + toa.TOA(t, obs="@", freq=np.inf, error=1 * u.ms, scale="tdb") for t in tdbs + ] + ) + ts.compute_TDBs(ephem="DE421") + ts.compute_posvels(ephem="DE421") + ph0 = m0.phase(ts) + ph1 = m1.phase(ts) + dphi = (ph1.int - ph0.int) + (ph1.frac - ph0.frac) + test_phase = np.zeros(len(tdbs)) + WAVEEPOCH = 55305 + WAVE_OM = 0.0015182579855022 + WAVE = [ + [-0.21573979911255, -0.049063841960712], + [0.62795320246729, -0.11984954655989], + [0.099618608456845, 0.28530756908788], + [-0.21537340649058, 0.18849486610628], + [0.021980474493165, -0.23566696662127], + ] + ph = (tdbs - WAVEEPOCH) * WAVE_OM + for i in range(5): + test_phase += WAVE[i][0] * np.sin((i + 1) * ph) + WAVE[i][1] * np.cos( + (i + 1) * ph + ) + test_phase *= m0.F0.quantity.to(u.Hz).value + assert np.allclose(test_phase, dphi) From 4b2c184a3c0e44eaf05dc7003f3262ae3e27e2b3 Mon Sep 17 00:00:00 2001 From: Matthew Kerr Date: Thu, 14 Nov 2024 10:46:27 -0500 Subject: [PATCH 03/10] Remove extraneous JUMPs from test .par file. --- tests/datafile/vela_wave.par | 27 --------------------------- 1 file changed, 27 deletions(-) diff --git a/tests/datafile/vela_wave.par b/tests/datafile/vela_wave.par index e8aba5a4c..889361b01 100644 --- a/tests/datafile/vela_wave.par +++ b/tests/datafile/vela_wave.par @@ -28,45 +28,18 @@ TZRMJD 55896.55312516020475 TZRFRQ 1370.2919919999999365 TZRSITE pks TRES 310.453 -EPHVER 2 CLK TT(TAI) MODE 1 UNITS TDB TIMEEPH FB90 DILATEFREQ N PLANET_SHAPIRO N -T2CMETHOD TEMPO NE_SW 9.961 CORRECT_TROPOSPHERE N EPHEM DE405 NITS 1 NTOA 339 CHI2R 328088.1294 298 -JUMP -B 20CM 0 0 -JUMP -B 40CM 0 0 -JUMP -B 50CM 0 0 -JUMP -dfb3_J0437_55319_56160 1 2.2e-07 0 -JUMP -dfb3_J0437_56160_60000 1 4.5e-07 0 -JUMP -pdfb1_128_ch 1 -3.5547e-06 0 -JUMP -pdfb1_2048_ch 1 -4.68829e-05 0 -JUMP -pdfb1_512_ch 1 -1.22813e-05 0 -JUMP -pdfb1_post_2006 1 -1.3e-07 0 -JUMP -pdfb1_pre_2006 1 -1.13e-06 0 -JUMP -pdfb2_1024_MHz 1 -5.435e-06 0 -JUMP -pdfb2_256MHz_1024_ch 1 -1.1395e-05 0 -JUMP -pdfb2_256MHz_512ch 1 -4.75e-06 0 -JUMP -pdfb3_1024_256_512 1 2.45e-06 0 -JUMP -pdfb3_1024_MHz 1 1.03e-06 0 -JUMP -pdfb3_256MHz_1024ch 1 4.295e-06 0 -JUMP -pdfb3_64MHz_1024ch 1 1.494e-05 0 -JUMP -pdfb3_64MHz_512ch 1 8.9e-06 0 -JUMP -pdfb4.*_1024_[1,2]... 1 2.23e-06 0 -JUMP -pdfb4_256MHz_1024ch 1 5.05e-06 0 -JUMP -pdfb4_55319_56055_cals 1 9.27e-07 0 -JUMP -pdfb4_56110_56160_cals 1 5.41e-07 0 -JUMP -pdfb4_56160_60000_cals 1 4.25e-07 0 -JUMP -wbb256_512_128_3p_b 1 -6.2e-07 0 -JUMP -wbb_c_config 1 3.8e-07 0 WAVEEPOCH 55305 WAVE_OM 0.0015182579855022 0 WAVE1 -0.21573979911255 -0.049063841960712 From 87939cd468d00cc8760080d6742e19d07f8270ce Mon Sep 17 00:00:00 2001 From: Matthew Kerr Date: Thu, 14 Nov 2024 12:44:48 -0500 Subject: [PATCH 04/10] Fix WAVE<-->WAVEX conversion. --- src/pint/utils.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/pint/utils.py b/src/pint/utils.py index da23092fc..4f6523455 100644 --- a/src/pint/utils.py +++ b/src/pint/utils.py @@ -1746,10 +1746,10 @@ def _translate_wave_freqs(om: Union[float, u.Quantity], k: int) -> u.Quantity: WXFREQ_ quantity in units 1/d that can be used in WaveX model """ if isinstance(om, u.quantity.Quantity): - om.to(u.d**-1) + om.to(u.rad*u.d**-1) else: - om *= u.d**-1 - return (om * (k + 1)) / (2.0 * np.pi) + om *= u.rad*u.d**-1 + return (om * (k + 1)) / (2.0 * np.pi * u.rad) def _translate_wavex_freqs(wxfreq: Union[float, u.Quantity], k: int) -> u.Quantity: @@ -1774,8 +1774,8 @@ def _translate_wavex_freqs(wxfreq: Union[float, u.Quantity], k: int) -> u.Quanti else: wxfreq *= u.d**-1 if len(wxfreq) == 1: - return (2.0 * np.pi * wxfreq) / (k + 1.0) - wave_om = [((2.0 * np.pi * wxfreq[i]) / (k[i] + 1.0)) for i in range(len(wxfreq))] + return (2.0 * np.pi * u.rad * wxfreq) / (k + 1.0) + wave_om = [((2.0 * np.pi * u.rad * wxfreq[i]) / (k[i] + 1.0)) for i in range(len(wxfreq))] return ( sum(wave_om) / len(wave_om) if np.allclose(wave_om, wave_om[0], atol=1e-3) From 4219a9b6294b8af91b4c4e838bc4440df287e598 Mon Sep 17 00:00:00 2001 From: Matthew Kerr Date: Thu, 14 Nov 2024 12:54:09 -0500 Subject: [PATCH 05/10] Simply astropy units syntax. --- src/pint/utils.py | 10 ++-------- 1 file changed, 2 insertions(+), 8 deletions(-) diff --git a/src/pint/utils.py b/src/pint/utils.py index 4f6523455..27af72c85 100644 --- a/src/pint/utils.py +++ b/src/pint/utils.py @@ -1745,10 +1745,7 @@ def _translate_wave_freqs(om: Union[float, u.Quantity], k: int) -> u.Quantity: astropy.units.Quantity WXFREQ_ quantity in units 1/d that can be used in WaveX model """ - if isinstance(om, u.quantity.Quantity): - om.to(u.rad*u.d**-1) - else: - om *= u.rad*u.d**-1 + om <<= u.rad/u.d return (om * (k + 1)) / (2.0 * np.pi * u.rad) @@ -1769,10 +1766,7 @@ def _translate_wavex_freqs(wxfreq: Union[float, u.Quantity], k: int) -> u.Quanti astropy.units.Quantity WAVEOM quantity in units 1/d that can be used in Wave model """ - if isinstance(wxfreq, u.quantity.Quantity): - wxfreq.to(u.d**-1) - else: - wxfreq *= u.d**-1 + wxfreq <<= u.d**-1 if len(wxfreq) == 1: return (2.0 * np.pi * u.rad * wxfreq) / (k + 1.0) wave_om = [((2.0 * np.pi * u.rad * wxfreq[i]) / (k[i] + 1.0)) for i in range(len(wxfreq))] From ba48280d2a0f20e10cd7e6bb13f8f6c35f49401a Mon Sep 17 00:00:00 2001 From: Matthew Kerr Date: Thu, 14 Nov 2024 15:10:11 -0500 Subject: [PATCH 06/10] black is going to make me run screaming into the night --- src/pint/utils.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/pint/utils.py b/src/pint/utils.py index 27af72c85..59f0d3076 100644 --- a/src/pint/utils.py +++ b/src/pint/utils.py @@ -1745,7 +1745,7 @@ def _translate_wave_freqs(om: Union[float, u.Quantity], k: int) -> u.Quantity: astropy.units.Quantity WXFREQ_ quantity in units 1/d that can be used in WaveX model """ - om <<= u.rad/u.d + om <<= u.rad / u.d return (om * (k + 1)) / (2.0 * np.pi * u.rad) @@ -1769,7 +1769,9 @@ def _translate_wavex_freqs(wxfreq: Union[float, u.Quantity], k: int) -> u.Quanti wxfreq <<= u.d**-1 if len(wxfreq) == 1: return (2.0 * np.pi * u.rad * wxfreq) / (k + 1.0) - wave_om = [((2.0 * np.pi * u.rad * wxfreq[i]) / (k[i] + 1.0)) for i in range(len(wxfreq))] + wave_om = [ + ((2.0 * np.pi * u.rad * wxfreq[i]) / (k[i] + 1.0)) for i in range(len(wxfreq)) + ] return ( sum(wave_om) / len(wave_om) if np.allclose(wave_om, wave_om[0], atol=1e-3) From ffc5a2d1cf2a5c1650fa4f0775508d307ca7685c Mon Sep 17 00:00:00 2001 From: Matthew Kerr Date: Fri, 15 Nov 2024 07:57:34 -0500 Subject: [PATCH 07/10] Add changelog entry. --- CHANGELOG-unreleased.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGELOG-unreleased.md b/CHANGELOG-unreleased.md index 2899b99f8..8f6c7f411 100644 --- a/CHANGELOG-unreleased.md +++ b/CHANGELOG-unreleased.md @@ -11,4 +11,5 @@ the released changes. ### Changed ### Added ### Fixed +- Changed WAVE\_OM units from 1/d to rad/d. ### Removed From ab4a22e923bdf2aadcdc1600a7d93bdeefdae7d5 Mon Sep 17 00:00:00 2001 From: David Kaplan Date: Fri, 15 Nov 2024 16:39:55 -0600 Subject: [PATCH 08/10] return name of added param, helpful when adding maskParameter --- src/pint/models/timing_model.py | 1 + 1 file changed, 1 insertion(+) diff --git a/src/pint/models/timing_model.py b/src/pint/models/timing_model.py index a0b27b090..68fe886bf 100644 --- a/src/pint/models/timing_model.py +++ b/src/pint/models/timing_model.py @@ -3425,6 +3425,7 @@ def add_param(self, param, deriv_func=None, setup=False): if deriv_func is not None: self.register_deriv_funcs(deriv_func, param.name) param._parent = self + return param.name def remove_param(self, param): """Remove a parameter from the Component. From 7c2bf806bff01136fb7e7df298539932b1ad0325 Mon Sep 17 00:00:00 2001 From: David Kaplan Date: Fri, 15 Nov 2024 16:40:16 -0600 Subject: [PATCH 09/10] fix logic when selection mask is array([0]) --- src/pint/models/noise_model.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/pint/models/noise_model.py b/src/pint/models/noise_model.py index b1af8e10f..299787fe5 100644 --- a/src/pint/models/noise_model.py +++ b/src/pint/models/noise_model.py @@ -164,14 +164,14 @@ def scale_toa_sigma(self, toas, warn=True): if equad.quantity is None: continue mask = equad.select_toa_mask(toas) - if np.any(mask): + if len(mask) > 0: sigma_scaled[mask] = np.hypot(sigma_scaled[mask], equad.quantity) elif warn: warnings.warn(f"EQUAD {equad} has no TOAs") for efac_name in self.EFACs: efac = getattr(self, efac_name) mask = efac.select_toa_mask(toas) - if np.any(mask): + if len(mask) > 0: sigma_scaled[mask] *= efac.quantity elif warn: warnings.warn(f"EFAC {efac} has no TOAs") From a17f4aabce38b60502d1c0f350980d884cfa9d9c Mon Sep 17 00:00:00 2001 From: David Kaplan Date: Mon, 18 Nov 2024 12:36:53 -0600 Subject: [PATCH 10/10] changelog --- CHANGELOG-unreleased.md | 2 ++ 1 file changed, 2 insertions(+) diff --git a/CHANGELOG-unreleased.md b/CHANGELOG-unreleased.md index d15e5e781..8e7b884e3 100644 --- a/CHANGELOG-unreleased.md +++ b/CHANGELOG-unreleased.md @@ -11,6 +11,8 @@ the released changes. ### Changed ### Added - When TCB->TDB conversion info is missing, will print parameter name +- `add_param` returns the name of the parameter (useful for numbered parameters) ### Fixed - When EQUAD is created from TNEQ, has proper TCB->TDB conversion info +- TOA selection masks will work when only TOA is the first one ### Removed