diff --git a/CHANGELOG-unreleased.md b/CHANGELOG-unreleased.md index 59a2a78b3..789b2062f 100644 --- a/CHANGELOG-unreleased.md +++ b/CHANGELOG-unreleased.md @@ -16,6 +16,7 @@ the released changes. - Type hinting for most of the `pint.utils` module - `funcParameter`s are no longer listed in the `pintk` interface. - Updated location of CCERA +- New algorithm for TCB <-> TDB conversion ### Added - `bayesian_information_criterion()` function - `dmx_setup` function @@ -24,6 +25,8 @@ the released changes. - Test for `pint.utils.split_swx()` - Custom type definitions for type hints - Added `citation.cff` +- `convert_tcb2tdb`, `tcb2tdb_scale_factor`, and `effective_dimensionality` attributes for `floatParameter`s, `MJDParameter`s, `AngleParameter`s, `maskParameter`s, and `prefixParameter`s. +- Documentation: HOWTO about determining tcb<->tdb scaling factors ### Fixed - `pint.utils.split_swx()` to use updated `SolarWindDispersionX()` parameter naming convention - Fix #1759 by changing order of comparison diff --git a/docs/examples/How_to_build_a_timing_model_component.py b/docs/examples/How_to_build_a_timing_model_component.py index 981c2e48e..45c78d2e5 100644 --- a/docs/examples/How_to_build_a_timing_model_component.py +++ b/docs/examples/How_to_build_a_timing_model_component.py @@ -107,6 +107,7 @@ def __init__(self): units=u.s, description="Spin period", longdouble=True, + tcb2tdb_scale_factor=u.Quantity(1), ) ) # Add spin period derivative P1. Since it is not all rquired, we are setting the @@ -118,6 +119,7 @@ def __init__(self): units=u.s / u.s, description="Spin period derivative", longdouble=True, + tcb2tdb_scale_factor=u.Quantity(1), ) ) # Add reference epoch time. @@ -126,6 +128,7 @@ def __init__(self): name="PEPOCH_P0", description="Reference epoch for spin-down", time_scale="tdb", + tcb2tdb_scale_factor=u.Quantity(1), ) ) # Add spindown phase model function to phase functions @@ -159,6 +162,7 @@ def F0(self): units="Hz", description="Spin-frequency", long_double=True, + tcb2tdb_scale_factor=u.Quantity(1), ) # Defining the derivatives. In the PINT, a common format of derivative naming is @@ -175,6 +179,7 @@ def F1(self): units=u.Hz / u.s, description="Spin down frequency", long_double=True, + tcb2tdb_scale_factor=u.Quantity(1), ) @property diff --git a/docs/examples/build_model_from_scratch.py b/docs/examples/build_model_from_scratch.py index 3c82b36c9..0310fb553 100644 --- a/docs/examples/build_model_from_scratch.py +++ b/docs/examples/build_model_from_scratch.py @@ -280,7 +280,11 @@ # %% # Add prefix parameters dmx_0003 = p.prefixParameter( - parameter_type="float", name="DMX_0003", value=None, units=u.pc / u.cm**3 + parameter_type="float", + name="DMX_0003", + value=None, + units=u.pc / u.cm**3, + tcb2tdb_scale_factor=u.Quantity(1), ) tm.components["DispersionDMX"].add_param(dmx_0003, setup=True) @@ -299,10 +303,18 @@ # %% dmxr1_0003 = p.prefixParameter( - parameter_type="mjd", name="DMXR1_0003", value=None, units=u.day + parameter_type="mjd", + name="DMXR1_0003", + value=None, + units=u.day, + tcb2tdb_scale_factor=u.Quantity(1), ) # DMXR1 is a type of MJD parameter internally. dmxr2_0003 = p.prefixParameter( - parameter_type="mjd", name="DMXR2_0003", value=None, units=u.day + parameter_type="mjd", + name="DMXR2_0003", + value=None, + units=u.day, + tcb2tdb_scale_factor=u.Quantity(1), ) # DMXR1 is a type of MJD parameter internally. tm.components["DispersionDMX"].add_param(dmxr1_0003, setup=True) @@ -346,7 +358,12 @@ # %% f1 = p.prefixParameter( - parameter_type="float", name="F1", value=0.0, units=u.Hz / (u.s), longdouble=True + parameter_type="float", + name="F1", + value=0.0, + units=u.Hz / (u.s), + longdouble=True, + tcb2tdb_scale_factor=u.Quantity(1), ) f2 = p.prefixParameter( @@ -355,6 +372,7 @@ value=0.0, units=u.Hz / (u.s) ** 2, longdouble=True, + tcb2tdb_scale_factor=u.Quantity(1), ) # %% diff --git a/docs/examples/understanding_parameters.py b/docs/examples/understanding_parameters.py index e5d1199d0..365cccfd7 100644 --- a/docs/examples/understanding_parameters.py +++ b/docs/examples/understanding_parameters.py @@ -84,18 +84,34 @@ # %% {"jupyter": {"outputs_hidden": false}} # You can specify the vaue as a number -t = pp.floatParameter(name="TEST", value=100, units="Hz", uncertainty=0.03) +t = pp.floatParameter( + name="TEST", + value=100, + units="Hz", + uncertainty=0.03, + tcb2tdb_scale_factor=u.Quantity(1), +) print(t) # %% {"jupyter": {"outputs_hidden": false}} # Or as a string that will be parsed -t2 = pp.floatParameter(name="TEST", value="200", units="Hz", uncertainty=".04") +t2 = pp.floatParameter( + name="TEST", + value="200", + units="Hz", + uncertainty=".04", + tcb2tdb_scale_factor=u.Quantity(1), +) print(t2) # %% {"jupyter": {"outputs_hidden": false}} # Or as an astropy Quantity with units (this is the preferred method!) t3 = pp.floatParameter( - name="TEST", value=0.3 * u.kHz, units="Hz", uncertainty=4e-5 * u.kHz + name="TEST", + value=0.3 * u.kHz, + units="Hz", + uncertainty=4e-5 * u.kHz, + tcb2tdb_scale_factor=u.Quantity(1), ) print(t3) print(t3.quantity) @@ -103,6 +119,10 @@ print(t3.uncertainty) print(t3.uncertainty_value) +# `tcb2tdb_scale_factor` is a quantity that controls how the parameter transforms +# under the TCB <-> TDB conversion. Please see the Explanation about Timescales and +# the How-To about TCB <-> TDB conversion for more details. + # %% [markdown] # ## Setting Parameters # diff --git a/docs/howto.rst b/docs/howto.rst index 0657a16e0..1a507c9e1 100644 --- a/docs/howto.rst +++ b/docs/howto.rst @@ -19,3 +19,4 @@ Please feel free to change this by writing more; there are also some entries on editing-documentation working-with-notebooks controlling-logging + tcb2tdb-factors diff --git a/docs/tcb2tdb-factors.rst b/docs/tcb2tdb-factors.rst new file mode 100644 index 000000000..65d8dc6b8 --- /dev/null +++ b/docs/tcb2tdb-factors.rst @@ -0,0 +1,71 @@ +.. highlight:: shell + +How to determine the scaling factor for TCB <-> TDB conversion for a parameter +------------------------------------------------------------------------------ + +The TCB and TDB timescales differ by a constant factor in the definition of the second. +This means that the epochs and TOAs must be scaled according to the following expression:: + + t_tdb = (t_tcb - IFTE_MJD0) / IFTE_K + IFTE_MJD0 + +Similarly, a time difference will be scaled as:: + + dt_tdb = dt_tcb / IFTE_K + +Since the definition of the second is changing, all parameters involved in the timing model +must also be transformed. In the simplest case, if a quantity x has dimensions of [T^n], it +will be transformed as:: + + x_tdb = x_tcb / IFTE_K^n + +This rule applies to the majority of parameters. + +However, there are some parameters in pulsar timing which appear in the timing model multiplied +my some constants. Examples include + + 1. DM appears as DMconst * DM + 2. A1 appears as A1 / c + 3. PX appears as PX * c / AU + 4. M2 appears as M2 * G / c^3 + +In these cases, it is customary to keep the multiplication factor numerically constant during +TCB <-> TDB conversion, and to absorb the effect of converting this constant into the parameter +itself. If a parameter has such a factor, it must be specified using the `tcb2tdb_scale_factor` +argument while constructing the `Parameter` object. For example, DM will have +`tcb2tdb_scale_factor=DMconst`. + +Note that the parameter multiplied by the constant in these cases has dimensions of the +form [T^n]. In the above cases, the value of n is as follows. + + 1. DM has n = -1 + 2. A1 has n = 1 + 3. PX has n = -1 + 4. M2 has n = 1 + +In general, if a parameter x appears in the timing model as C*x and if C*x has dimensionality of +the form [T^n], the scaling should be done with the "effective dimensionality" n. + +If a parameter doesn't have a dimensionality of [T^n], a general rule is to reorganize the +factors in the equation such that each group has a dimensionality [T^n]. This is ALWAYS possible +because the timing model components produce either a delay ([T^1]) or a phase ([T^0]). + +A useful trick for doing this type of transformation is to express the parameters in geometrized +units, where everything has dimensions of [T^n]. For example, a mass M will be expressed as M*(G/c^3) +(e.g., M2, MTOT), and a distance L will be expressed as L/c (e.g., A1). Please note that this may not +work in every case, and each parameter should be treated in a case-by-case basis depending on the +delay/phase expression they appear in. + +Exceptions to this are noise parameters. The TOA uncertainties are measured in the observatory +timescale and the are not converted into TCB or TDB before computing the likelihood function/ +chi-squared. Hence, we don't convert the quantities that modify TOA uncertainties, namely EFACs and`` +EQUADs. Since we are not converting TOA variances, it doesn't make sense to convert TOA covariances +either. Hence, ECORRs and red and DM noise parameters are not converted. This means that +the noise parameters must ALWAYS be re-estimated after a TCB <-> TDB conversion. + +Another exception to this are FD parameters and FD jumps, which involve polynomials of logarithms. +Such functions transform in a non-linear fashion during TCB <-> TDB conversion, and we have chosen +not to apply the conversion to such parameters. They too must be re-estimated after a TCB <-> TDB +conversion. + +In such cases, the parameter can be excluded from TCB <-> TDB conversion by passing `convert_tcb2tdb=False` +into the `Parameter` class constructor. \ No newline at end of file diff --git a/src/pint/models/absolute_phase.py b/src/pint/models/absolute_phase.py index 940bea474..2db5ede34 100644 --- a/src/pint/models/absolute_phase.py +++ b/src/pint/models/absolute_phase.py @@ -1,4 +1,5 @@ -"""Timing model absolute phase (TZRMJD, TZRSITE ...)""" +"""Timing model absolute phase""" + import astropy.units as u from loguru import logger as log @@ -33,6 +34,7 @@ def __init__(self): name="TZRMJD", description="Epoch of the zero phase TOA.", time_scale="utc", + convert_tcb2tdb=False, ) ) self.add_param( @@ -45,6 +47,7 @@ def __init__(self): name="TZRFRQ", units=u.MHz, description="The frequency of the zero phase TOA.", + convert_tcb2tdb=False, ) ) self.tz_cache = None diff --git a/src/pint/models/astrometry.py b/src/pint/models/astrometry.py index 005d09349..7c83584ac 100644 --- a/src/pint/models/astrometry.py +++ b/src/pint/models/astrometry.py @@ -8,6 +8,7 @@ import astropy.constants as const import astropy.coordinates as coords import astropy.units as u +import astropy.constants as consts import numpy as np from astropy.time import Time @@ -50,11 +51,18 @@ def __init__(self): name="POSEPOCH", description="Reference epoch for position", time_scale="tdb", + tcb2tdb_scale_factor=u.Quantity(1), ) ) self.add_param( - floatParameter(name="PX", units="mas", value=0.0, description="Parallax") + floatParameter( + name="PX", + units="mas", + value=0.0, + description="Parallax", + tcb2tdb_scale_factor=(consts.c / u.au), + ) ) self.delay_funcs_component += [self.solar_system_geometric_delay] @@ -162,7 +170,9 @@ def solar_system_geometric_delay( re_dot_L = np.sum(tbl["ssb_obs_pos"][c] * L_hat, axis=1) delay[c] = -re_dot_L.to(ls).value if self.PX.value != 0.0: + # This is equivalent to PX * c / AU. L = (1.0 / self.PX.value) * u.kpc + # TODO: np.sum currently loses units in some cases... re_sqr = ( np.sum(tbl["ssb_obs_pos"][c] ** 2, axis=1) @@ -278,6 +288,7 @@ def __init__(self): units="H:M:S", description="Right ascension (J2000)", aliases=["RA"], + tcb2tdb_scale_factor=u.Quantity(1), ) ) @@ -287,6 +298,7 @@ def __init__(self): units="D:M:S", description="Declination (J2000)", aliases=["DEC"], + tcb2tdb_scale_factor=u.Quantity(1), ) ) @@ -296,6 +308,7 @@ def __init__(self): units="mas/year", description="Proper motion in RA", value=0.0, + tcb2tdb_scale_factor=u.Quantity(1), ) ) @@ -305,6 +318,7 @@ def __init__(self): units="mas/year", description="Proper motion in DEC", value=0.0, + tcb2tdb_scale_factor=u.Quantity(1), ) ) self.set_special_params(["RAJ", "DECJ", "PMRA", "PMDEC"]) @@ -755,6 +769,7 @@ def __init__(self): units="deg", description="Ecliptic longitude", aliases=["LAMBDA"], + tcb2tdb_scale_factor=u.Quantity(1), ) ) @@ -764,6 +779,7 @@ def __init__(self): units="deg", description="Ecliptic latitude", aliases=["BETA"], + tcb2tdb_scale_factor=u.Quantity(1), ) ) @@ -774,6 +790,7 @@ def __init__(self): description="Proper motion in ecliptic longitude", aliases=["PMLAMBDA"], value=0.0, + tcb2tdb_scale_factor=u.Quantity(1), ) ) @@ -784,6 +801,7 @@ def __init__(self): description="Proper motion in ecliptic latitude", aliases=["PMBETA"], value=0.0, + tcb2tdb_scale_factor=u.Quantity(1), ) ) diff --git a/src/pint/models/binary_bt.py b/src/pint/models/binary_bt.py index 6a40ae9f4..8c4969c2d 100644 --- a/src/pint/models/binary_bt.py +++ b/src/pint/models/binary_bt.py @@ -1,19 +1,18 @@ """The BT (Blandford & Teukolsky) model.""" + import numpy as np from pint.models.parameter import floatParameter from pint.models.pulsar_binary import PulsarBinary from pint.models.stand_alone_psr_binaries.BT_model import BTmodel from pint.models.stand_alone_psr_binaries.BT_piecewise import BTpiecewise -from pint.models.timing_model import MissingParameter, TimingModel +from pint.models.timing_model import MissingParameter import astropy.units as u -from pint import GMsun, Tsun, ls -from astropy.table import Table +import astropy.constants as consts +from pint import ls from astropy.time import Time from pint.models.parameter import ( - MJDParameter, floatParameter, prefixParameter, - maskParameter, ) from pint.toa_select import TOASelect @@ -58,6 +57,7 @@ def __init__(self): value=0.0, units="second", description="Time dilation & gravitational redshift", + tcb2tdb_scale_factor=u.Quantity(1), ) ) self.remove_param("M2") @@ -81,11 +81,6 @@ def validate(self): self.GAMMA.frozen = True -"""The BT (Blandford & Teukolsky) model with piecewise orbital parameters. -See Blandford & Teukolsky 1976, ApJ, 205, 580. -""" - - class BinaryBTPiecewise(PulsarBinary): """Model implementing the BT model with piecewise orbital parameters A1X and T0X. This model lets the user specify time ranges and fit for a different piecewise orbital parameter in each time range, This is a PINT pulsar binary BTPiecewise model class, a subclass of PulsarBinary. @@ -111,6 +106,7 @@ def __init__(self): value=0.0, units="second", description="Time dilation & gravitational redshift", + tcb2tdb_scale_factor=u.Quantity(1), ) ) self.A1_value_funcs = [] @@ -137,7 +133,16 @@ def add_group_range( piece_index : int Number to label the piece being added. """ - if group_start_mjd is not None and group_end_mjd is not None: + if group_start_mjd is None: + if group_end_mjd is not None: + group_start_mjd = group_end_mjd - 100 + else: + group_start_mjd = 50000 + group_end_mjd = 60000 + + elif group_end_mjd is None: + group_end_mjd = group_start_mjd + 100 + else: if isinstance(group_start_mjd, Time): group_start_mjd = group_start_mjd.mjd elif isinstance(group_start_mjd, u.quantity.Quantity): @@ -147,22 +152,11 @@ def add_group_range( elif isinstance(group_end_mjd, u.quantity.Quantity): group_end_mjd = group_end_mjd.value - elif group_start_mjd is None or group_end_mjd is None: - if group_start_mjd is None and group_end_mjd is not None: - group_start_mjd = group_end_mjd - 100 - elif group_start_mjd is not None and group_end_mjd is None: - group_end_mjd = group_start_mjd + 100 - else: - group_start_mjd = 50000 - group_end_mjd = 60000 - if piece_index is None: dct = self.get_prefix_mapping_component("XR1_") - if len(list(dct.keys())) > 0: - piece_index = np.max(list(dct.keys())) + 1 - else: - piece_index = 0 - + piece_index = ( + np.max(list(dct.keys())) + 1 if len(list(dct.keys())) > 0 else 0 + ) # check the validity of the desired group to add if group_end_mjd is not None and group_start_mjd is not None: @@ -174,7 +168,7 @@ def add_group_range( ) elif piece_index > 9999: raise ValueError( - f"Invalid index for group. Cannot index beyond 9999 (yet?)" + "Invalid index for group. Cannot index beyond 9999 (yet?)" ) i = f"{int(piece_index):04d}" @@ -186,6 +180,7 @@ def add_group_range( parameter_type="MJD", time_scale="utc", value=group_start_mjd, + tcb2tdb_scale_factor=u.Quantity(1), ) ) self.add_param( @@ -196,6 +191,7 @@ def add_group_range( parameter_type="MJD", time_scale="utc", value=group_end_mjd, + tcb2tdb_scale_factor=u.Quantity(1), ) ) self.setup() @@ -208,11 +204,7 @@ def remove_range(self, index): Number or list/array of numbers corresponding to T0X/A1X indices to be removed from model. """ - if ( - isinstance(index, int) - or isinstance(index, float) - or isinstance(index, np.int64) - ): + if isinstance(index, (int, float, np.int64)): indices = [index] elif not isinstance(index, list) or not isinstance(index, np.ndarray): raise TypeError( @@ -223,14 +215,16 @@ def remove_range(self, index): for prefix in ["T0X_", "A1X_", "XR1_", "XR2_"]: if hasattr(self, f"{prefix+index_rf}"): self.remove_param(prefix + index_rf) - if hasattr(self.binary_instance, "param_pieces"): - if len(self.binary_instance.param_pieces) > 0: - temp_piece_information = [] - for item in self.binary_instance.param_pieces: - if item[0] != index_rf: - temp_piece_information.append(item) - self.binary_instance.param_pieces = temp_piece_information - # self.binary_instance.param_pieces = self.binary_instance.param_pieces.remove('index_rf') + if ( + hasattr(self.binary_instance, "param_pieces") + and len(self.binary_instance.param_pieces) > 0 + ): + temp_piece_information = [ + item + for item in self.binary_instance.param_pieces + if item[0] != index_rf + ] + self.binary_instance.param_pieces = temp_piece_information self.validate() self.setup() @@ -251,29 +245,27 @@ def add_piecewise_param(self, piece_index, **kwargs): param = key paramx = kwargs[key] - if key == "T0": - param_unit = u.d + if key == "A1": + param_unit = ls if isinstance(paramx, u.quantity.Quantity): paramx = paramx.value - elif isinstance(paramx, np.float128): + elif isinstance(paramx, np.float64): paramx = paramx - elif isinstance(paramx, Time): - paramx = paramx.mjd else: raise ValueError( - "Unspported data type '%s' for piecewise T0. Ensure the piecewise parameter value is a np.float128, Time or astropy.quantity.Quantity" - % type(paramx) + f"Unspported data type '{type(paramx)}' for piecewise A1. Ensure the piecewise parameter value is a np.float64 or astropy.quantity.Quantity" ) - elif key == "A1": - param_unit = ls + elif key == "T0": + param_unit = u.d if isinstance(paramx, u.quantity.Quantity): paramx = paramx.value - elif isinstance(paramx, np.float64): + elif isinstance(paramx, np.float128): paramx = paramx + elif isinstance(paramx, Time): + paramx = paramx.mjd else: raise ValueError( - "Unspported data type '%s' for piecewise A1. Ensure the piecewise parameter value is a np.float64 or astropy.quantity.Quantity" - % type(paramx) + f"Unspported data type '{type(paramx)}' for piecewise T0. Ensure the piecewise parameter value is a np.float128, Time or astropy.quantity.Quantity" ) key_found = True @@ -283,36 +275,33 @@ def add_piecewise_param(self, piece_index, **kwargs): ) if piece_index is None: - dct = self.get_prefix_mapping_component(param + "X_") - if len(list(dct.keys())) > 0: - piece_index = np.max(list(dct.keys())) + 1 - else: - piece_index = 0 - elif int(piece_index) in self.get_prefix_mapping_component(param + "X_"): + dct = self.get_prefix_mapping_component(f"{param}X_") + piece_index = ( + np.max(list(dct.keys())) + 1 if len(list(dct.keys())) > 0 else 0 + ) + elif int(piece_index) in self.get_prefix_mapping_component(f"{param}X_"): raise ValueError( - "Index '%s' is already in use in this model. Please choose another." - % piece_index + f"Index '{piece_index}' is already in use in this model. Please choose another." ) i = f"{int(piece_index):04d}" # handling if None are passed as arguments - if any(i is None for i in [param, param_unit, paramx]): - if param is not None: - # if parameter value or unit unset, set with default according to param - if param_unit is None: - param_unit = (getattr(self, param)).units - if paramx is None: - paramx = (getattr(self, param)).value + if any(i is None for i in [param, param_unit, paramx]) and param is not None: + if param_unit is None: + param_unit = (getattr(self, param)).units + if paramx is None: + paramx = (getattr(self, param)).value # check if name exists and is currently available self.add_param( prefixParameter( - name=param + f"X_{i}", + name=f"{param}X_{i}", units=param_unit, value=paramx, - description="Parameter" + param + "variation", + description=f"Parameter{param}variation", parameter_type="float", frozen=False, + tcb2tdb_scale_factor=(1 / consts.c if param == "A1" else u.Quantity(1)), ) ) self.setup() @@ -355,36 +344,33 @@ def setup(self): ) if hasattr(self, f"A1X_{piece_index}"): - if hasattr(self, f"A1X_{piece_index}"): - if getattr(self, f"A1X_{piece_index}") is not None: - self.binary_instance.add_binary_params( - f"A1X_{piece_index}", getattr(self, f"A1X_{piece_index}") - ) - else: - self.binary_instance.add_binary_params( - f"A1X_{piece_index}", self.A1.value - ) - - if hasattr(self, f"XR1_{piece_index}"): - if getattr(self, f"XR1_{piece_index}") is not None: + if getattr(self, f"A1X_{piece_index}") is not None: self.binary_instance.add_binary_params( - f"XR1_{piece_index}", getattr(self, f"XR1_{piece_index}") + f"A1X_{piece_index}", getattr(self, f"A1X_{piece_index}") ) else: - raise ValueError(f"No date provided to create a group with") - else: - raise ValueError(f"No name provided to create a group with") - - if hasattr(self, f"XR2_{piece_index}"): - if getattr(self, f"XR2_{piece_index}") is not None: self.binary_instance.add_binary_params( - f"XR2_{piece_index}", getattr(self, f"XR2_{piece_index}") + f"A1X_{piece_index}", self.A1.value ) - else: - raise ValueError(f"No date provided to create a group with") + + if not hasattr(self, f"XR1_{piece_index}"): + raise ValueError("No name provided to create a group with") + + if getattr(self, f"XR1_{piece_index}") is not None: + self.binary_instance.add_binary_params( + f"XR1_{piece_index}", getattr(self, f"XR1_{piece_index}") + ) else: - raise ValueError(f"No name provided to create a group with") + raise ValueError("No date provided to create a group with") + if not hasattr(self, f"XR2_{piece_index}"): + raise ValueError("No name provided to create a group with") + if getattr(self, f"XR2_{piece_index}") is not None: + self.binary_instance.add_binary_params( + f"XR2_{piece_index}", getattr(self, f"XR2_{piece_index}") + ) + else: + raise ValueError("No date provided to create a group with") self.update_binary_object(None) def validate(self): @@ -401,7 +387,7 @@ def validate(self): super().validate() for p in ("T0", "A1"): if getattr(self, p).value is None: - raise MissingParameter("BT", p, "%s is required for BT" % p) + raise MissingParameter("BT", p, f"{p} is required for BT") # If any *DOT is set, we need T0 for p in ("PBDOT", "OMDOT", "EDOT", "A1DOT"): @@ -409,9 +395,8 @@ def validate(self): getattr(self, p).set("0") getattr(self, p).frozen = True - if getattr(self, p).value is not None: - if self.T0.value is None: - raise MissingParameter("BT", "T0", "T0 is required if *DOT is set") + if getattr(self, p).value is not None and self.T0.value is None: + raise MissingParameter("BT", "T0", "T0 is required if *DOT is set") if self.GAMMA.value is None: self.GAMMA.set("0") @@ -438,30 +423,33 @@ def validate(self): ) if len(np.setdiff1d(j_plb, j_pub)) > 0: raise ValueError( - f"Group index mismatch error. Check the indexes of XR1_/XR2_ parameters in the model" + "Group index mismatch error. Check the indexes of XR1_/XR2_ parameters in the model" + ) + if ( + len(ls_A1X) <= 0 + and (len(ls_pub) > 0 and len(ls_T0X) > 0) + and len(np.setdiff1d(j_pub, j_T0X)) > 0 + ): + raise ValueError( + "Group index mismatch error. Check the indexes of T0X groups, make sure they match there are corresponding group ranges (XR1/XR2)" + ) + if ( + len(ls_T0X) <= 0 + and (len(ls_pub) > 0 and len(ls_A1X) > 0) + and len(np.setdiff1d(j_pub, j_A1X)) > 0 + ): + raise ValueError( + "Group index mismatch error. Check the indexes of A1X groups, make sure they match there are corresponding group ranges (/XR2)" ) - if not len(ls_A1X) > 0: - if len(ls_pub) > 0 and len(ls_T0X) > 0: - if len(np.setdiff1d(j_pub, j_T0X)) > 0: - raise ValueError( - f"Group index mismatch error. Check the indexes of T0X groups, make sure they match there are corresponding group ranges (XR1/XR2)" - ) - if not len(ls_T0X) > 0: - if len(ls_pub) > 0 and len(ls_A1X) > 0: - if len(np.setdiff1d(j_pub, j_A1X)) > 0: - raise ValueError( - f"Group index mismatch error. Check the indexes of A1X groups, make sure they match there are corresponding group ranges (/XR2)" - ) lb = [(getattr(self, tup[1])).value for tup in ls_plb] ub = [(getattr(self, tup[1])).value for tup in ls_pub] for i in range(len(lb)): for j in range(len(lb)): - if i != j: - if max(lb[i], lb[j]) < min(ub[i], ub[j]): - raise ValueError( - f"Group boundary overlap detected. Make sure groups are not overlapping" - ) + if i != j and max(lb[i], lb[j]) < min(ub[i], ub[j]): + raise ValueError( + "Group boundary overlap detected. Make sure groups are not overlapping" + ) def paramx_per_toa(self, param_name, toas): """Find the piecewise parameter value each toa will reference during calculations @@ -481,16 +469,15 @@ def paramx_per_toa(self, param_name, toas): XR2_mapping = self.get_prefix_mapping_component("XR2_") if not hasattr(self, "toas_selector"): self.toas_selector = TOASelect(is_range=True) - if param_name[0:2] == "T0": + if param_name[:2] == "T0": paramX_mapping = self.get_prefix_mapping_component("T0X_") param_unit = u.d - elif param_name[0:2] == "A1": + elif param_name[:2] == "A1": paramX_mapping = self.get_prefix_mapping_component("A1X_") param_unit = ls else: raise AttributeError( - "param '%s' not found. Please choose another. Currently implemented: 'T0' or 'A1' " - % param_name + f"param '{param_name}' not found. Please choose another. Currently implemented: 'T0' or 'A1' " ) for piece_index in paramX_mapping.keys(): r1 = getattr(self, XR1_mapping[piece_index]).quantity @@ -502,7 +489,7 @@ def paramx_per_toa(self, param_name, toas): paramx[v] += getattr(self, k).quantity for i in range(len(paramx)): if paramx[i] == 0: - paramx[i] = (getattr(self, param_name[0:2])).value * param_unit + paramx[i] = getattr(self, param_name[:2]).value * param_unit return paramx diff --git a/src/pint/models/binary_dd.py b/src/pint/models/binary_dd.py index b07df462f..cdff58eb7 100644 --- a/src/pint/models/binary_dd.py +++ b/src/pint/models/binary_dd.py @@ -1,9 +1,10 @@ """Damour and Deruelle binary model.""" + import numpy as np from astropy import units as u, constants as c from pint import Tsun -from pint.models.parameter import floatParameter, funcParameter, intParameter +from pint.models.parameter import floatParameter, funcParameter from pint.models.pulsar_binary import PulsarBinary from pint.models.stand_alone_psr_binaries.DD_model import DDmodel from pint.models.stand_alone_psr_binaries.DDS_model import DDSmodel @@ -63,6 +64,7 @@ def __init__( value=0.0, units="s", description="DD model aberration parameter A0", + tcb2tdb_scale_factor=u.Quantity(1), ) ) @@ -72,6 +74,7 @@ def __init__( value=0.0, units="s", description="DD model aberration parameter B0", + tcb2tdb_scale_factor=u.Quantity(1), ) ) @@ -81,6 +84,7 @@ def __init__( value=0.0, units="second", description="Time dilation & gravitational redshift", + tcb2tdb_scale_factor=u.Quantity(1), ) ) @@ -90,6 +94,7 @@ def __init__( value=0.0, units="", description="Relativistic deformation of the orbit", + tcb2tdb_scale_factor=u.Quantity(1), ) ) @@ -100,6 +105,7 @@ def __init__( value=0.0, units="", description="Relativistic deformation of the orbit", + tcb2tdb_scale_factor=u.Quantity(1), ) ) @@ -173,7 +179,11 @@ def __init__( self.add_param( floatParameter( - name="SHAPMAX", value=0.0, description="Function of inclination angle" + name="SHAPMAX", + value=0.0, + units=u.dimensionless_unscaled, + description="Function of inclination angle", + tcb2tdb_scale_factor=u.Quantity(1), ) ) self.remove_param("SINI") @@ -251,6 +261,7 @@ def __init__( name="MTOT", units=u.M_sun, description="Total system mass in units of Solar mass", + tcb2tdb_scale_factor=(c.G / c.c**3), ) ) self.add_param( @@ -259,6 +270,7 @@ def __init__( units="deg/year", description="Excess longitude of periastron advance compared to GR", long_double=True, + tcb2tdb_scale_factor=u.Quantity(1), ) ) self.add_param( @@ -269,6 +281,7 @@ def __init__( unit_scale=True, scale_factor=1e-12, scale_threshold=1e-7, + tcb2tdb_scale_factor=u.Quantity(1), ) ) for p in ["OMDOT", "PBDOT", "GAMMA", "SINI", "DR", "DTH"]: @@ -404,6 +417,7 @@ def __init__(self): units="second", description="Shapiro delay parameter H3 as in Freire and Wex 2010 Eq(20)", long_double=True, + tcb2tdb_scale_factor=u.Quantity(1), ) ) self.add_param( @@ -413,6 +427,7 @@ def __init__(self): description="Shapiro delay parameter STIGMA as in Freire and Wex 2010 Eq(12)", long_double=True, aliases=["VARSIGMA", "STIG"], + tcb2tdb_scale_factor=u.Quantity(1), ) ) self.remove_param("M2") diff --git a/src/pint/models/binary_ddk.py b/src/pint/models/binary_ddk.py index 6b40f5ebc..8a31f6209 100644 --- a/src/pint/models/binary_ddk.py +++ b/src/pint/models/binary_ddk.py @@ -1,4 +1,5 @@ """The DDK model - Damour and Deruelle with kinematics.""" + import warnings import numpy as np from astropy import units as u @@ -113,7 +114,11 @@ def __init__( self.add_param( floatParameter( - name="KIN", value=0.0, units="deg", description="Inclination angle" + name="KIN", + value=0.0, + units="deg", + description="Inclination angle", + tcb2tdb_scale_factor=u.Quantity(1), ) ) self.add_param( @@ -122,6 +127,7 @@ def __init__( value=0.0, units="deg", description="The longitude of the ascending node", + tcb2tdb_scale_factor=u.Quantity(1), ) ) self.add_param( diff --git a/src/pint/models/binary_ell1.py b/src/pint/models/binary_ell1.py index aaa6092e9..a34b3c186 100644 --- a/src/pint/models/binary_ell1.py +++ b/src/pint/models/binary_ell1.py @@ -110,7 +110,10 @@ def __init__(self): self.add_param( MJDParameter( - name="TASC", description="Epoch of ascending node", time_scale="tdb" + name="TASC", + description="Epoch of ascending node", + time_scale="tdb", + tcb2tdb_scale_factor=u.Quantity(1), ) ) @@ -120,6 +123,7 @@ def __init__(self): units="", description="First Laplace-Lagrange parameter, ECC*sin(OM)", long_double=True, + tcb2tdb_scale_factor=u.Quantity(1), ) ) @@ -129,6 +133,7 @@ def __init__(self): units="", description="Second Laplace-Lagrange parameter, ECC*cos(OM)", long_double=True, + tcb2tdb_scale_factor=u.Quantity(1), ) ) @@ -138,6 +143,7 @@ def __init__(self): units="1e-12/s", description="First derivative of first Laplace-Lagrange parameter", long_double=True, + tcb2tdb_scale_factor=u.Quantity(1), ) ) @@ -147,6 +153,7 @@ def __init__(self): units="1e-12/s", description="Second derivative of first Laplace-Lagrange parameter", long_double=True, + tcb2tdb_scale_factor=u.Quantity(1), ) ) self.remove_param("ECC") @@ -342,6 +349,7 @@ def __init__(self): units="second", description="Shapiro delay parameter H3 as in Freire and Wex 2010 Eq(20)", long_double=True, + tcb2tdb_scale_factor=u.Quantity(1), ) ) @@ -351,6 +359,7 @@ def __init__(self): units="second", description="Shapiro delay parameter H4 as in Freire and Wex 2010 Eq(21)", long_double=True, + tcb2tdb_scale_factor=u.Quantity(1), ) ) @@ -361,6 +370,7 @@ def __init__(self): description="Shapiro delay parameter STIGMA as in Freire and Wex 2010 Eq(12)", long_double=True, aliases=["VARSIGMA", "STIG"], + tcb2tdb_scale_factor=u.Quantity(1), ) ) self.add_param( @@ -458,6 +468,7 @@ def __init__(self): units="deg/year", description="Rate of advance of periastron", long_double=True, + tcb2tdb_scale_factor=u.Quantity(1), ) ) @@ -467,6 +478,7 @@ def __init__(self): units="1/year", description="Log-derivative of the eccentricity EDOT/ECC", long_double=True, + tcb2tdb_scale_factor=u.Quantity(1), ) ) diff --git a/src/pint/models/dispersion_model.py b/src/pint/models/dispersion_model.py index 44148ad6c..d626754dc 100644 --- a/src/pint/models/dispersion_model.py +++ b/src/pint/models/dispersion_model.py @@ -20,7 +20,6 @@ split_prefixed_name, taylor_horner, taylor_horner_deriv, - get_prefix_timeranges, ) from pint import DMconst @@ -154,6 +153,7 @@ def __init__(self): value=0.0, description="Dispersion measure", long_double=True, + tcb2tdb_scale_factor=DMconst, ) ) self.add_param( @@ -165,11 +165,15 @@ def __init__(self): description_template=self.DM_dervative_description, type_match="float", long_double=True, + tcb2tdb_scale_factor=DMconst, ) ) self.add_param( MJDParameter( - name="DMEPOCH", description="Epoch of DM measurement", time_scale="tdb" + name="DMEPOCH", + description="Epoch of DM measurement", + time_scale="tdb", + tcb2tdb_scale_factor=u.Quantity(1), ) ) @@ -320,13 +324,16 @@ class DispersionDMX(Dispersion): def __init__(self): super().__init__() + # DMX is for info output right now + # @abhisrkckl: What exactly is the use of this parameter? self.add_param( floatParameter( name="DMX", units="pc cm^-3", value=0.0, description="Dispersion measure", + convert_tcb2tdb=False, ) ) @@ -396,6 +403,7 @@ def add_DMX_range(self, mjd_start, mjd_end, index=None, dmx=0, frozen=True): description="Dispersion measure variation", parameter_type="float", frozen=frozen, + tcb2tdb_scale_factor=DMconst, ) ) self.add_param( @@ -406,6 +414,7 @@ def add_DMX_range(self, mjd_start, mjd_end, index=None, dmx=0, frozen=True): parameter_type="MJD", time_scale="utc", value=mjd_start, + tcb2tdb_scale_factor=u.Quantity(1), ) ) self.add_param( @@ -416,6 +425,7 @@ def add_DMX_range(self, mjd_start, mjd_end, index=None, dmx=0, frozen=True): parameter_type="MJD", time_scale="utc", value=mjd_end, + tcb2tdb_scale_factor=u.Quantity(1), ) ) self.setup() @@ -512,6 +522,7 @@ def add_DMX_ranges(self, mjd_starts, mjd_ends, indices=None, dmxs=0, frozens=Tru description="Dispersion measure variation", parameter_type="float", frozen=frozen, + tcb2tdb_scale_factor=DMconst, ) ) self.add_param( @@ -522,6 +533,7 @@ def add_DMX_ranges(self, mjd_starts, mjd_ends, indices=None, dmxs=0, frozens=Tru parameter_type="MJD", time_scale="utc", value=mjd_start, + tcb2tdb_scale_factor=u.Quantity(1), ) ) self.add_param( @@ -532,6 +544,7 @@ def add_DMX_ranges(self, mjd_starts, mjd_ends, indices=None, dmxs=0, frozens=Tru parameter_type="MJD", time_scale="utc", value=mjd_end, + tcb2tdb_scale_factor=u.Quantity(1), ) ) self.setup() @@ -735,7 +748,8 @@ def __init__(self): name="DMJUMP", units="pc cm^-3", value=None, - description="DM value offset.", + description="Wideband DM value offset.", + convert_tcb2tdb=False, ) ) @@ -825,6 +839,7 @@ def __init__(self): units="pc cm^-3", value=None, description="System-dependent DM offset.", + tcb2tdb_scale_factor=DMconst, ) ) diff --git a/src/pint/models/dmwavex.py b/src/pint/models/dmwavex.py index d61560dd2..258eb3556 100644 --- a/src/pint/models/dmwavex.py +++ b/src/pint/models/dmwavex.py @@ -1,4 +1,5 @@ """DM variations expressed as a sum of sinusoids.""" + import astropy.units as u import numpy as np from loguru import logger as log @@ -7,7 +8,7 @@ from pint.models.parameter import MJDParameter, prefixParameter from pint.models.timing_model import MissingParameter from pint.models.dispersion_model import Dispersion -from pint import dmu +from pint import DMconst, dmu class DMWaveX(Dispersion): @@ -35,6 +36,7 @@ def __init__(self): name="DMWXEPOCH", description="Reference epoch for Fourier representation of DM noise", time_scale="tdb", + tcb2tdb_scale_factor=u.Quantity(1), ) ) self.add_dmwavex_component(0.1, index=1, dmwxsin=0, dmwxcos=0, frozen=False) @@ -93,6 +95,7 @@ def add_dmwavex_component( units="1/d", value=dmwxfreq, parameter_type="float", + tcb2tdb_scale_factor=u.Quantity(1), ) ) self.add_param( @@ -103,6 +106,7 @@ def add_dmwavex_component( value=dmwxsin, frozen=frozen, parameter_type="float", + tcb2tdb_scale_factor=DMconst, ) ) self.add_param( @@ -113,6 +117,7 @@ def add_dmwavex_component( value=dmwxcos, frozen=frozen, parameter_type="float", + tcb2tdb_scale_factor=DMconst, ) ) self.setup() @@ -204,6 +209,7 @@ def add_dmwavex_components( units="1/d", value=dmwxfreq, parameter_type="float", + tcb2tdb_scale_factor=u.Quantity(1), ) ) self.add_param( @@ -214,6 +220,7 @@ def add_dmwavex_components( value=dmwxsin, parameter_type="float", frozen=frozen, + tcb2tdb_scale_factor=DMconst, ) ) self.add_param( @@ -224,6 +231,7 @@ def add_dmwavex_components( value=dmwxcos, parameter_type="float", frozen=frozen, + tcb2tdb_scale_factor=DMconst, ) ) self.setup() diff --git a/src/pint/models/fdjump.py b/src/pint/models/fdjump.py index 9cec3354e..a5258defe 100644 --- a/src/pint/models/fdjump.py +++ b/src/pint/models/fdjump.py @@ -79,6 +79,7 @@ def __init__(self): name=f"FD{j}JUMP", units="second", description=f"System-dependent FD parameter of polynomial index {j}", + convert_tcb2tdb=False, ) ) diff --git a/src/pint/models/frequency_dependent.py b/src/pint/models/frequency_dependent.py index 3094a8a9a..2215b5c08 100644 --- a/src/pint/models/frequency_dependent.py +++ b/src/pint/models/frequency_dependent.py @@ -1,4 +1,5 @@ """Frequency-dependent delays to model profile evolution.""" + from warnings import warn import astropy.units as u @@ -43,6 +44,7 @@ def __init__(self): descriptionTplt=self._description_template, # unitTplt=lambda x: "second", type_match="float", + convert_tcb2tdb=False, ) ) diff --git a/src/pint/models/glitch.py b/src/pint/models/glitch.py index d89eec6b8..f198a27c2 100644 --- a/src/pint/models/glitch.py +++ b/src/pint/models/glitch.py @@ -1,4 +1,5 @@ """Pulsar timing glitches.""" + import astropy.units as u import numpy as np @@ -20,31 +21,31 @@ class Glitch(PhaseComponent): @classmethod def _description_glitch_phase(cls, x): - return "Phase change for glitch %d" % x + return f"Phase change for glitch {x}" @classmethod def _description_glitch_epoch(cls, x): - return "Epoch of glitch %d" % x + return f"Epoch of glitch {x}" @classmethod def _description_glitch_frequencychange(cls, x): - return ("Permanent frequency change for glitch %d" % x,) + return (f"Permanent frequency change for glitch {x}",) @classmethod def _description_glitch_frequencyderivativechange(cls, x): - return ("Permanent frequency-derivative change for glitch %d" % x,) + return (f"Permanent frequency-derivative change for glitch {x}",) @classmethod def _description_glitch_frequencysecondderivativechange(cls, x): - return ("Permanent second frequency-derivative change for glitch %d" % x,) + return (f"Permanent second frequency-derivative change for glitch {x}",) @classmethod def _description_decaying_frequencychange(cls, x): - return "Decaying frequency change for glitch %d " % x + return f"Decaying frequency change for glitch {x}" @classmethod def _description_decaytimeconstant(cls, x): - return "Decay time constant for glitch %d" % x + return f"Decay time constant for glitch {x}" register = True category = "glitch" @@ -59,6 +60,7 @@ def __init__(self): value=0.0, description_template=self._description_glitch_phase, type_match="float", + tcb2tdb_scale_factor=u.Quantity(1), ) ) self.add_param( @@ -68,6 +70,7 @@ def __init__(self): description_template=self._description_glitch_epoch, parameter_type="MJD", time_scale="tdb", + tcb2tdb_scale_factor=u.Quantity(1), ) ) self.add_param( @@ -77,6 +80,7 @@ def __init__(self): value=0.0, description_template=self._description_glitch_frequencychange, type_match="float", + tcb2tdb_scale_factor=u.Quantity(1), ) ) self.add_param( @@ -85,6 +89,7 @@ def __init__(self): units="Hz/s", value=0.0, description_template=self._description_glitch_frequencyderivativechange, + tcb2tdb_scale_factor=u.Quantity(1), ) ) self.add_param( @@ -93,6 +98,7 @@ def __init__(self): units="Hz/s^2", value=0.0, description_template=self._description_glitch_frequencysecondderivativechange, + tcb2tdb_scale_factor=u.Quantity(1), ) ) self.add_param( @@ -102,6 +108,7 @@ def __init__(self): value=0.0, description_template=self._description_decaying_frequencychange, type_match="float", + tcb2tdb_scale_factor=u.Quantity(1), ) ) @@ -112,6 +119,7 @@ def __init__(self): value=0.0, description_template=self._description_decaytimeconstant, type_match="float", + tcb2tdb_scale_factor=u.Quantity(1), ) ) self.phase_funcs_component += [self.glitch_phase] @@ -208,7 +216,7 @@ def glitch_phase(self, toas, delay): tau = getattr(self, "GLTD_%d" % idx).quantity decayterm = dF0D * tau * (1.0 - np.exp(-(dt[affected] / tau))) else: - decayterm = 0.0 * u.Unit("") + decayterm = u.Quantity(0) log.debug(f"Glitch phase for glitch {idx}: {dphs} {dphs.unit}") phs[affected] += ( diff --git a/src/pint/models/ifunc.py b/src/pint/models/ifunc.py index 91a3b79cb..323380ac3 100644 --- a/src/pint/models/ifunc.py +++ b/src/pint/models/ifunc.py @@ -1,4 +1,5 @@ """Tabulated extra delays.""" + import astropy.units as u import numpy as np @@ -56,7 +57,12 @@ def __init__(self): super().__init__() self.add_param( - floatParameter(name="SIFUNC", description="Type of interpolation", units="") + floatParameter( + name="SIFUNC", + description="Type of interpolation", + units="", + convert_tcb2tdb=False, + ) ) self.add_param( prefixParameter( @@ -66,6 +72,7 @@ def __init__(self): type_match="pair", long_double=True, parameter_type="pair", + convert_tcb2tdb=False, ) ) self.phase_funcs_component += [self.ifunc_phase] diff --git a/src/pint/models/jump.py b/src/pint/models/jump.py index b133040d2..75a972acd 100644 --- a/src/pint/models/jump.py +++ b/src/pint/models/jump.py @@ -28,7 +28,11 @@ class DelayJump(DelayComponent): def __init__(self): super().__init__() - self.add_param(maskParameter(name="JUMP", units="second")) + self.add_param( + maskParameter( + name="JUMP", units="second", tcb2tdb_scale_factor=u.Quantity(1) + ) + ) self.delay_funcs_component += [self.jump_delay] def setup(self): @@ -94,6 +98,7 @@ def __init__(self): name="JUMP", units="second", description="Phase jump for selection.", + tcb2tdb_scale_factor=u.Quantity(1), ) ) self.phase_funcs_component += [self.jump_phase] @@ -216,6 +221,7 @@ def add_jump_and_flags(self, toa_table): value=0.0, units="second", frozen=False, + tcb2tdb_scale_factor=u.Quantity(1), ) # otherwise add on jump with next index else: @@ -237,6 +243,7 @@ def add_jump_and_flags(self, toa_table): value=0.0, units="second", frozen=False, + tcb2tdb_scale_factor=u.Quantity(1), ) self.add_param(param) ind = param.index diff --git a/src/pint/models/model_builder.py b/src/pint/models/model_builder.py index 677da6b86..9e3998ce3 100644 --- a/src/pint/models/model_builder.py +++ b/src/pint/models/model_builder.py @@ -991,10 +991,7 @@ def guess_binary_model(parfile_dict): def add_sini(parameters): """If 'KIN' is a model parameter, Tempo2 doesn't really use SINI""" - if "KIN" in parameters: - return list(parameters) + ["SINI"] - else: - return list(parameters) + return list(parameters) + ["SINI"] if "KIN" in parameters else list(parameters) all_components = AllComponents() binary_models = all_components.category_component_map["pulsar_system"] @@ -1006,7 +1003,7 @@ def add_sini(parameters): ) for binary_model in binary_models } - binary_parameters_map.update({"Isolated": []}) + binary_parameters_map["Isolated"] = [] all_binary_parameters = { parname for parnames in binary_parameters_map.values() for parname in parnames } @@ -1056,7 +1053,7 @@ def convert_binary_params_dict( parameters if they exist. """ binary = parfile_dict.get("BINARY", None) - binary = binary if not binary else binary[0] + binary = binary[0] if binary else binary log.debug(f"Requested to convert binary model for BINARY model: {binary}") if binary: @@ -1067,11 +1064,8 @@ def convert_binary_params_dict( ) if not binary_model_guesses: - error_message = f"Unable to determine binary model for this par file" - log_message = ( - f"Unable to determine the binary model based" - f"on the model parameters in the par file." - ) + error_message = "Unable to determine binary model for this par file" + log_message = "Unable to determine the binary model basedon the model parameters in the par file." log.error(log_message) raise UnknownBinaryModel(error_message) @@ -1085,21 +1079,21 @@ def convert_binary_params_dict( # Convert KIN if requested if convert_komkin and "KIN" in parfile_dict: log.info(f"Converting KOM to/from IAU <--> DT96: {parfile_dict['KIN']}") - log.debug(f"Converting KIN to/from IAU <--> DT96") + log.debug("Converting KIN to/from IAU <--> DT96") entries = parfile_dict["KIN"][0].split() new_value = _convert_kin(float(entries[0]) * u.deg).value parfile_dict["KIN"] = [" ".join([repr(new_value)] + entries[1:])] # Convert KOM if requested if convert_komkin and "KOM" in parfile_dict: - log.debug(f"Converting KOM to/from IAU <--> DT96") + log.debug("Converting KOM to/from IAU <--> DT96") entries = parfile_dict["KOM"][0].split() new_value = _convert_kom(float(entries[0]) * u.deg).value parfile_dict["KOM"] = [" ".join([repr(new_value)] + entries[1:])] # Drop SINI if requested if drop_ddk_sini and binary_model_guesses[0] == "DDK": - log.debug(f"Dropping SINI from DDK model") + log.debug("Dropping SINI from DDK model") parfile_dict.pop("SINI", None) return parfile_dict diff --git a/src/pint/models/noise_model.py b/src/pint/models/noise_model.py index debe7efdc..8bc64d379 100644 --- a/src/pint/models/noise_model.py +++ b/src/pint/models/noise_model.py @@ -59,6 +59,7 @@ def __init__( units="", aliases=["T2EFAC", "TNEF"], description="A multiplication factor on the measured TOA uncertainties,", + convert_tcb2tdb=False, ) ) @@ -68,6 +69,7 @@ def __init__( units="us", aliases=["T2EQUAD"], description="An error term added in quadrature to the scaled (by EFAC) TOA uncertainty.", + convert_tcb2tdb=False, ) ) @@ -76,6 +78,7 @@ def __init__( name="TNEQ", units=u.LogUnit(physical_unit=u.second), description="An error term added in quadrature to the scaled (by EFAC) TOA uncertainty in units of log10(second).", + convert_tcb2tdb=False, ) ) self.covariance_matrix_funcs += [self.sigma_scaled_cov_matrix] @@ -238,6 +241,7 @@ def __init__( name="DMEFAC", units="", description="A multiplication factor on the measured DM uncertainties,", + convert_tcb2tdb=False, ) ) @@ -246,6 +250,7 @@ def __init__( name="DMEQUAD", units="pc / cm ^ 3", description="An error term added in quadrature to the scaled (by EFAC) TOA uncertainty.", + convert_tcb2tdb=False, ) ) @@ -347,6 +352,7 @@ def __init__( units="us", aliases=["TNECORR"], description="An error term that is correlated among all TOAs in an observing epoch.", + convert_tcb2tdb=False, ) ) @@ -475,6 +481,7 @@ def __init__( units="", aliases=[], description="Amplitude of powerlaw DM noise in tempo2 format", + convert_tcb2tdb=False, ) ) self.add_param( @@ -483,6 +490,7 @@ def __init__( units="", aliases=[], description="Spectral index of powerlaw " "DM noise in tempo2 format", + convert_tcb2tdb=False, ) ) self.add_param( @@ -491,6 +499,7 @@ def __init__( units="", aliases=[], description="Number of DM noise frequencies.", + convert_tcb2tdb=False, ) ) @@ -585,6 +594,7 @@ def __init__( units="", aliases=[], description="Amplitude of powerlaw red noise.", + convert_tcb2tdb=False, ) ) self.add_param( @@ -593,6 +603,7 @@ def __init__( units="", aliases=[], description="Spectral index of powerlaw red noise.", + convert_tcb2tdb=False, ) ) @@ -602,6 +613,7 @@ def __init__( units="", aliases=[], description="Amplitude of powerlaw red noise in tempo2 format", + convert_tcb2tdb=False, ) ) self.add_param( @@ -610,6 +622,7 @@ def __init__( units="", aliases=[], description="Spectral index of powerlaw red noise in tempo2 format", + convert_tcb2tdb=False, ) ) self.add_param( @@ -618,6 +631,7 @@ def __init__( units="", aliases=[], description="Number of red noise frequencies.", + convert_tcb2tdb=False, ) ) diff --git a/src/pint/models/parameter.py b/src/pint/models/parameter.py index c94e62f6e..a32347347 100644 --- a/src/pint/models/parameter.py +++ b/src/pint/models/parameter.py @@ -21,6 +21,7 @@ parameters PINT understands. """ + import numbers from warnings import warn @@ -577,7 +578,7 @@ def from_parfile_line(self, line): return True def value_as_latex(self): - return f"${self.as_ufloat():.1uSL}$" if not self.frozen else f"{self.value:f}" + return f"{self.value:f}" if self.frozen else f"${self.as_ufloat():.1uSL}$" def as_latex(self): try: @@ -649,6 +650,11 @@ class floatParameter(Parameter): parameter exist. long_double : bool, optional, default False A flag specifying whether value is float or long double. + convert_tcb2tdb: bool + Whether to convert this parameter during TCB <-> TDB conversion. + tcb2tdb_scale_factor: astropy.units.Quantity + The scaling factor to be applied while computing the effective + dimensionality. The default is 1. Example ------- @@ -672,6 +678,8 @@ def __init__( unit_scale=False, scale_factor=None, scale_threshold=None, + convert_tcb2tdb=True, + tcb2tdb_scale_factor=None, **kwargs, ): self.long_double = long_double @@ -699,10 +707,18 @@ def __init__( ] self.unit_scale = unit_scale + assert ( + not convert_tcb2tdb or tcb2tdb_scale_factor is not None + ), "Please specify the tcb2tdb_scale_factor explicitly." + self.convert_tcb2tdb = convert_tcb2tdb + self.tcb2tdb_scale_factor = tcb2tdb_scale_factor + @property def long_double(self): """Whether the parameter has long double precision.""" # FIXME: why not just always keep long double precision? + # AS: Because it is not supported in Mac M* machines, and we + # should avoid long doubles as much as we can. return self._long_double @long_double.setter @@ -754,8 +770,8 @@ def _set_quantity(self, val): 1. Astropy quantity 2. float 3. string - """ + # Check long_double if not self._long_double: setfunc_with_unit = _identity_function @@ -850,6 +866,13 @@ def from_ufloat(self, value, units=None): self.quantity = value.n * units self.uncertainty = value.s * units if value.s > 0 else None + @property + def effective_dimensionality(self) -> int: + """Compute the effective dimensionality for TCB <-> TDB conversion.""" + return compute_effective_dimensionality( + self.quantity, self.tcb2tdb_scale_factor + ) + class strParameter(Parameter): """String-valued parameter. @@ -1088,6 +1111,8 @@ def __init__( continuous=True, aliases=None, time_scale="tdb", + convert_tcb2tdb=True, + tcb2tdb_scale_factor=None, **kwargs, ): self._time_scale = time_scale @@ -1106,6 +1131,12 @@ def __init__( self.paramType = "MJDParameter" self.special_arg += ["time_scale"] + assert ( + not convert_tcb2tdb or tcb2tdb_scale_factor is not None + ), "Please specify the tcb2tdb_scale_factor explicitly." + self.convert_tcb2tdb = convert_tcb2tdb + self.tcb2tdb_scale_factor = tcb2tdb_scale_factor + def str_quantity(self, quan): return time_to_mjd_string(quan) @@ -1212,6 +1243,13 @@ def as_ufloat(self): """ return ufloat(self.value, self.uncertainty_value) + @property + def effective_dimensionality(self) -> int: + """Compute the effective dimensionality for TCB <-> TDB conversion.""" + return compute_effective_dimensionality( + self.value * u.day, self.tcb2tdb_scale_factor + ) + class AngleParameter(Parameter): """Parameter in angle units. @@ -1239,6 +1277,11 @@ class AngleParameter(Parameter): aliases : list, optional An optional list of strings specifying alternate names that can also be accepted for this parameter. + convert_tcb2tdb: bool + Whether to convert this parameter during TCB <-> TDB conversion. + tcb2tdb_scale_factor: astropy.units.Quantity + The scaling factor to be applied while computing the effective + dimensionality. The default is 1. Example ------- @@ -1258,6 +1301,8 @@ def __init__( frozen=True, continuous=True, aliases=None, + convert_tcb2tdb=True, + tcb2tdb_scale_factor=None, **kwargs, ): self._str_unit = units @@ -1287,6 +1332,12 @@ def __init__( aliases=aliases, ) + assert ( + not convert_tcb2tdb or tcb2tdb_scale_factor is not None + ), "Please specify the tcb2tdb_scale_factor explicitly." + self.convert_tcb2tdb = convert_tcb2tdb + self.tcb2tdb_scale_factor = tcb2tdb_scale_factor + def _get_value(self, quan): # return Angle(x * self.unit_identifier[units.lower()][0]) return quan.value @@ -1372,6 +1423,13 @@ def as_ufloat(self, units=None): error = self.uncertainty.to_value(units) if self.uncertainty is not None else 0 return ufloat(value, error) + @property + def effective_dimensionality(self) -> int: + """Compute the effective dimensionality for TCB <-> TDB conversion.""" + return compute_effective_dimensionality( + self.quantity, self.tcb2tdb_scale_factor + ) + class prefixParameter: """Families of parameters identified by a prefix like ``DMX_0123``. @@ -1422,6 +1480,12 @@ class prefixParameter: Set float type quantity and value in numpy long doubles. time_scale : str, optional Time scale for MJDParameter class. + convert_tcb2tdb: bool + Whether to convert this parameter during TCB <-> TDB conversion. + tcb2tdb_scale_factor: astropy.units.Quantity or function + The scaling factor to be applied while computing the effective + dimensionality. If this is a function, it should take the prefix as + argument and return the scaling factor. The default is 1. """ def __init__( @@ -1442,6 +1506,8 @@ def __init__( scale_factor=None, scale_threshold=None, time_scale="utc", + convert_tcb2tdb=True, + tcb2tdb_scale_factor=None, **kwargs, ): # Split prefixed name, if the name is not in the prefixed format, error @@ -1483,6 +1549,18 @@ def __init__( real_description = input_description aliases = [pa + self.idxfmt for pa in self.prefix_aliases] self.long_double = long_double + + # For prefix parameters, the scaling factor can in principle be + # a function of the prefix. + assert ( + not convert_tcb2tdb or tcb2tdb_scale_factor is not None + ), "Please specify the tcb2tdb_scale_factor explicitly." + tcb2tdb_scale_factor_val = ( + tcb2tdb_scale_factor(self.prefix) + if hasattr(tcb2tdb_scale_factor, "__call__") + else tcb2tdb_scale_factor + ) + # initiate parameter class self.param_comp = self.param_class( name=self.name, @@ -1498,6 +1576,8 @@ def __init__( unit_scale=unit_scale, scale_factor=scale_factor, scale_threshold=scale_threshold, + convert_tcb2tdb=convert_tcb2tdb, + tcb2tdb_scale_factor=tcb2tdb_scale_factor_val, ) self.is_prefix = True self.time_scale = time_scale @@ -1598,6 +1678,14 @@ def description(self, val): def special_arg(self): return self.param_comp.special_arg + @property + def convert_tcb2tdb(self): + return self.param_comp.convert_tcb2tdb + + @property + def tcb2tdb_scale_factor(self): + return self.param_comp.tcb2tdb_scale_factor + def __repr__(self): return self.param_comp.__repr__() @@ -1657,6 +1745,8 @@ def new_param(self, index, inheritfrozen=False): "long_double", "time_scale", "parameter_type", + "convert_tcb2tdb", + "tcb2tdb_scale_factor", ] if hasattr(self, key) and (key != "frozen" or inheritfrozen) } @@ -1683,6 +1773,11 @@ def as_ufloat(self, units=None): error = self.uncertainty.to_value(units) if self.uncertainty is not None else 0 return ufloat(value, error) + @property + def effective_dimensionality(self) -> int: + """Compute the effective dimensionality for TCB <-> TDB conversion.""" + return self.param_comp.effective_dimensionality + class maskParameter(floatParameter): """Parameter that applies to a subset of TOAs. @@ -1737,6 +1832,11 @@ class maskParameter(floatParameter): Whether derivatives with respect to this parameter make sense. aliases : list, optional List of aliases for parameter name. + convert_tcb2tdb: bool + Whether to convert this parameter during TCB <-> TDB conversion. + tcb2tdb_scale_factor: astropy.units.Quantity + The scaling factor to be applied while computing the effective + dimensionality. The default is 1. """ # TODO: Is mask parameter provide some other type of parameters other then floatParameter? @@ -1755,6 +1855,8 @@ def __init__( frozen=True, continuous=False, aliases=[], + convert_tcb2tdb=True, + tcb2tdb_scale_factor=None, ): self.is_mask = True # {key_name: (keyvalue parse function, keyvalue length)} @@ -1805,6 +1907,8 @@ def __init__( continuous=continuous, aliases=idx_aliases + aliases, long_double=long_double, + convert_tcb2tdb=convert_tcb2tdb, + tcb2tdb_scale_factor=tcb2tdb_scale_factor, ) # For the first mask parameter, add name to aliases for the reading @@ -2002,6 +2106,8 @@ def new_param(self, index, copy_all=False): frozen=self.frozen, continuous=self.continuous, aliases=self.prefix_aliases, + convert_tcb2tdb=self.convert_tcb2tdb, + tcb2tdb_scale_factor=self.tcb2tdb_scale_factor, ) if copy_all else maskParameter( @@ -2010,6 +2116,8 @@ def new_param(self, index, copy_all=False): long_double=self.long_double, units=self.units, aliases=self.prefix_aliases, + convert_tcb2tdb=self.convert_tcb2tdb, + tcb2tdb_scale_factor=self.tcb2tdb_scale_factor, ) ) @@ -2239,10 +2347,12 @@ def str_quantity(self, quan): try: # Maybe it's a singleton quantity return floatParameter.str_quantity(self, quan) - except AttributeError: + except AttributeError as e: # Not a quantity, let's hope it's a list of length two? if len(quan) != 2: - raise ValueError(f"Don't know how to print this as a pair: {quan}") + raise ValueError( + f"Don't know how to print this as a pair: {quan}" + ) from e v0 = quan[0].to(self.units).value v1 = quan[1].to(self.units).value @@ -2486,3 +2596,21 @@ def as_parfile_line(self, format="pint"): if self.inpar else f"# {super().as_parfile_line(format=format)}" ) + + +def compute_effective_dimensionality( + quantity: u.Quantity, scaling_factor: u.Quantity +) -> int: + """Compute the effective dimensionality for TCB <-> TDB conversion.""" + unit = (quantity * scaling_factor).si.unit + + if len(unit.bases) == 0 or unit.bases == [u.rad]: + return 0 + elif unit.bases == [u.s]: + return unit.powers[0] + elif set(unit.bases) == {u.s, u.rad}: + return unit.powers[unit.bases.index(u.s)] + else: + raise ValueError( + "The scaled quantity has an unsupported unit. Check the scaling_factor.", + ) diff --git a/src/pint/models/phase_offset.py b/src/pint/models/phase_offset.py index 3793b5bca..b044b5867 100644 --- a/src/pint/models/phase_offset.py +++ b/src/pint/models/phase_offset.py @@ -27,6 +27,7 @@ def __init__(self): value=0.0, units="", description="Overall phase offset between physical TOAs and the TZR TOA.", + tcb2tdb_scale_factor=u.Quantity(1), ) ) self.phase_funcs_component += [self.offset_phase] diff --git a/src/pint/models/piecewise.py b/src/pint/models/piecewise.py index 9e78cd2ff..9b9fdd8b3 100644 --- a/src/pint/models/piecewise.py +++ b/src/pint/models/piecewise.py @@ -1,4 +1,5 @@ """Pulsar timing piecewise spin-down solution.""" + import astropy.units as u import numpy as np @@ -57,6 +58,7 @@ def __init__(self): description_template=self._description_solution_epochstart, parameter_type="MJD", time_scale="tdb", + tcb2tdb_scale_factor=u.Quantity(1), ) ) @@ -67,6 +69,7 @@ def __init__(self): description_template=self._description_solution_epochstop, parameter_type="MJD", time_scale="tdb", + tcb2tdb_scale_factor=u.Quantity(1), ) ) self.add_param( @@ -76,6 +79,7 @@ def __init__(self): description_template=self._description_solution_epoch, parameter_type="MJD", time_scale="tdb", + tcb2tdb_scale_factor=u.Quantity(1), ) ) self.add_param( @@ -86,6 +90,7 @@ def __init__(self): description_template=self._description_solution_startphase, type_match="float", uncertainty=1, + tcb2tdb_scale_factor=u.Quantity(1), ) ) self.add_param( @@ -95,6 +100,7 @@ def __init__(self): value=0.0, description_template=self._description_solution_frequency, type_match="float", + tcb2tdb_scale_factor=u.Quantity(1), ) ) self.add_param( @@ -103,6 +109,7 @@ def __init__(self): units="Hz/s", value=0.0, description_template=self._description_solution_frequencyderivative, + tcb2tdb_scale_factor=u.Quantity(1), ) ) self.add_param( @@ -111,6 +118,7 @@ def __init__(self): units="Hz/s^2", value=0.0, description_template=self._description_solution_secondfrequencyderivative, + tcb2tdb_scale_factor=u.Quantity(1), ) ) diff --git a/src/pint/models/pulsar_binary.py b/src/pint/models/pulsar_binary.py index b27dc923c..5f81e29ac 100644 --- a/src/pint/models/pulsar_binary.py +++ b/src/pint/models/pulsar_binary.py @@ -4,11 +4,11 @@ as PINT timing models. """ - import astropy.units as u import contextlib import numpy as np from astropy.coordinates import SkyCoord +import astropy.constants as consts from loguru import logger as log from pint import ls @@ -87,7 +87,11 @@ def __init__(self): self.binary_model_class = None self.add_param( floatParameter( - name="PB", units=u.day, description="Orbital period", long_double=True + name="PB", + units=u.day, + description="Orbital period", + long_double=True, + tcb2tdb_scale_factor=u.Quantity(1), ) ) self.add_param( @@ -98,6 +102,7 @@ def __init__(self): unit_scale=True, scale_factor=1e-12, scale_threshold=1e-7, + tcb2tdb_scale_factor=u.Quantity(1), ) ) self.add_param( @@ -105,6 +110,7 @@ def __init__(self): name="A1", units=ls, description="Projected semi-major axis of pulsar orbit, ap*sin(i)", + tcb2tdb_scale_factor=(1 / consts.c), ) ) # NOTE: the DOT here takes the value and times 1e-12, tempo/tempo2 can @@ -118,11 +124,16 @@ def __init__(self): unit_scale=True, scale_factor=1e-12, scale_threshold=1e-7, + tcb2tdb_scale_factor=(1 / consts.c), ) ) self.add_param( floatParameter( - name="ECC", units="", aliases=["E"], description="Eccentricity" + name="ECC", + units="", + aliases=["E"], + description="Eccentricity", + tcb2tdb_scale_factor=u.Quantity(1), ) ) self.add_param( @@ -133,11 +144,15 @@ def __init__(self): unit_scale=True, scale_factor=1e-12, scale_threshold=1e-7, + tcb2tdb_scale_factor=u.Quantity(1), ) ) self.add_param( MJDParameter( - name="T0", description="Epoch of periastron passage", time_scale="tdb" + name="T0", + description="Epoch of periastron passage", + time_scale="tdb", + tcb2tdb_scale_factor=u.Quantity(1), ) ) self.add_param( @@ -146,6 +161,7 @@ def __init__(self): units=u.deg, description="Longitude of periastron", long_double=True, + tcb2tdb_scale_factor=u.Quantity(1), ) ) self.add_param( @@ -154,6 +170,7 @@ def __init__(self): units="deg/year", description="Rate of advance of periastron", long_double=True, + tcb2tdb_scale_factor=u.Quantity(1), ) ) self.add_param( @@ -161,11 +178,15 @@ def __init__(self): name="M2", units=u.M_sun, description="Companion mass", + tcb2tdb_scale_factor=(consts.G / consts.c**3), ) ) self.add_param( floatParameter( - name="SINI", units="", description="Sine of inclination angle" + name="SINI", + units="", + description="Sine of inclination angle", + tcb2tdb_scale_factor=u.Quantity(1), ) ) self.add_param( @@ -179,6 +200,7 @@ def __init__(self): description_template=self.FBX_description, type_match="float", long_double=True, + tcb2tdb_scale_factor=u.Quantity(1), ) ) self.internal_params = [] @@ -310,11 +332,11 @@ def check_required_params(self, required_params): method_name = f"{p.lower()}_func" try: par_method = getattr(self.binary_instance, method_name) - except AttributeError: + except AttributeError as e: raise MissingParameter( self.binary_model_name, f"{p} is required for '{self.binary_model_name}'.", - ) + ) from e par_method() # With new parameter class set up, do we need this? @@ -402,13 +424,13 @@ def update_binary_object(self, toas, acc_delay=None): if hasattr(self._parent, par) or set(alias).intersection(self.params): try: pint_bin_name = self._parent.match_param_aliases(par) - except UnknownParameter: + except UnknownParameter as e: if par in self.internal_params: pint_bin_name = par else: raise UnknownParameter( f"Unable to find {par} in the parent model" - ) + ) from e binObjpar = getattr(self._parent, pint_bin_name) # make sure we aren't passing along derived parameters to the binary instance diff --git a/src/pint/models/solar_wind_dispersion.py b/src/pint/models/solar_wind_dispersion.py index 0f220e2e6..9ca243d66 100644 --- a/src/pint/models/solar_wind_dispersion.py +++ b/src/pint/models/solar_wind_dispersion.py @@ -1,4 +1,5 @@ """Dispersion due to the solar wind.""" + from warnings import warn import astropy.constants as const @@ -7,8 +8,9 @@ import numpy as np import scipy.special +from pint import DMconst from pint.models.dispersion_model import Dispersion -from pint.models.parameter import floatParameter, prefixParameter +from pint.models.parameter import floatParameter, intParameter, prefixParameter import pint.utils from pint.models.timing_model import MissingTOAs from pint.toa_select import TOASelect @@ -297,6 +299,7 @@ def __init__(self): value=0.0, aliases=["NE1AU", "SOLARN0"], description="Solar Wind density at 1 AU", + tcb2tdb_scale_factor=(const.c * DMconst), ) ) self.add_param( @@ -305,13 +308,13 @@ def __init__(self): value=2.0, units="", description="Solar Wind Model radial power-law index (only for SWM=1)", + tcb2tdb_scale_factor=u.Quantity(1), ) ) self.add_param( - floatParameter( + intParameter( name="SWM", - value=0.0, - units="", + value=0, description="Solar Wind Model (0 is from Edwards+ 2006, 1 is from You+2007,2012/Hazboun+ 2022)", ) ) @@ -740,6 +743,7 @@ def add_swx_range( description="Max Solar Wind DM", parameter_type="float", frozen=frozen, + tcb2tdb_scale_factor=DMconst, ) ) self.add_param( @@ -749,6 +753,7 @@ def add_swx_range( description="Solar wind power-law index", parameter_type="float", frozen=frozen, + tcb2tdb_scale_factor=u.Quantity(1), ) ) self.add_param( @@ -759,6 +764,7 @@ def add_swx_range( parameter_type="MJD", time_scale="utc", value=mjd_start, + tcb2tdb_scale_factor=u.Quantity(1), ) ) self.add_param( @@ -769,6 +775,7 @@ def add_swx_range( parameter_type="MJD", time_scale="utc", value=mjd_end, + tcb2tdb_scale_factor=u.Quantity(1), ) ) self.setup() diff --git a/src/pint/models/spindown.py b/src/pint/models/spindown.py index ff9972968..8a8bce411 100644 --- a/src/pint/models/spindown.py +++ b/src/pint/models/spindown.py @@ -1,4 +1,5 @@ """Polynomial pulsar spindown.""" + # spindown.py # Defines Spindown timing model class import astropy.units as u @@ -56,6 +57,7 @@ def __init__(self): description_template=self.F_description, type_match="float", long_double=True, + tcb2tdb_scale_factor=u.Quantity(1), ) ) # self.add_param( @@ -75,6 +77,7 @@ def __init__(self): name="PEPOCH", description="Reference epoch for spin-down", time_scale="tdb", + tcb2tdb_scale_factor=u.Quantity(1), ) ) diff --git a/src/pint/models/tcb_conversion.py b/src/pint/models/tcb_conversion.py index 13dd65fc6..266e70316 100644 --- a/src/pint/models/tcb_conversion.py +++ b/src/pint/models/tcb_conversion.py @@ -2,7 +2,14 @@ import numpy as np -from pint.models.parameter import MJDParameter +from pint.models.timing_model import TimingModel +from pint.models.parameter import ( + AngleParameter, + MJDParameter, + floatParameter, + maskParameter, + prefixParameter, +) from loguru import logger as log __all__ = [ @@ -19,7 +26,7 @@ IFTE_K = 1 + IFTE_KM1 -def scale_parameter(model, param, n, backwards): +def scale_parameter(model: TimingModel, param: str, n: int, backwards: bool) -> None: """Scale a parameter x by a power of IFTE_K x_tdb = x_tcb * IFTE_K**n @@ -52,14 +59,14 @@ def scale_parameter(model, param, n, backwards): factor = IFTE_K ** (p * n) - if hasattr(model, param) and getattr(model, param).quantity is not None: - par = getattr(model, param) + if (param in model) and model[param].quantity is not None: + par = model[param] par.value *= factor if par.uncertainty_value is not None: par.uncertainty_value *= factor -def transform_mjd_parameter(model, param, backwards): +def transform_mjd_parameter(model: TimingModel, param: str, backwards: bool) -> None: """Convert an MJD from TCB to TDB or vice versa. t_tdb = (t_tcb - IFTE_MJD0) / IFTE_K + IFTE_MJD0 t_tcb = (t_tdb - IFTE_MJD0) * IFTE_K + IFTE_MJD0 @@ -76,42 +83,35 @@ def transform_mjd_parameter(model, param, backwards): factor = IFTE_K if backwards else 1 / IFTE_K tref = IFTE_MJD0 - if hasattr(model, param) and getattr(model, param).quantity is not None: - par = getattr(model, param) - assert isinstance(par, MJDParameter) + if (param in model) and model[param].quantity is not None: + par = model[param] + assert isinstance(par, MJDParameter) or ( + isinstance(par, prefixParameter) + and isinstance(par.param_comp, MJDParameter) + ) par.value = (par.value - tref) * factor + tref if par.uncertainty_value is not None: par.uncertainty_value *= factor -def convert_tcb_tdb(model, backwards=False): +def convert_tcb_tdb(model: TimingModel, backwards: bool = False) -> None: """This function performs a partial conversion of a model specified in TCB to TDB. While this should be sufficient as a starting point, the resulting parameters are only approximate and the model should be re-fit. - This is based on the `transform` plugin of tempo2. - - The following parameters are converted to TDB: - 1. Spin frequency, its derivatives and spin epoch - 2. Sky coordinates, proper motion and the position epoch - 3. DM, DM derivatives and DM epoch - 4. Keplerian binary parameters and FB1 + This is roughly based on the `transform` plugin of tempo2, but uses + a different algorithm and does a more complete conversion. The following parameters are NOT converted although they are in fact affected by the TCB to TDB conversion: - 1. Parallax - 2. TZRMJD and TZRFRQ - 2. DMX parameters - 3. Solar wind parameters - 4. Binary post-Keplerian parameters including Shapiro delay - parameters (except FB1) - 5. Jumps and DM Jumps - 6. FD parameters - 7. EQUADs - 8. Red noise parameters including FITWAVES, powerlaw red noise and - powerlaw DM noise parameters + 1. TZRMJD and TZRFRQ + 2. DMJUMPs (the wideband kind) + 3. FD parameters and FD jumps + 4. EQUADs and ECORRs. + 5. GP Red noise and GP DM noise parameters + 6. Pair parameters such as Wave and IFunc parameters Parameters ---------- @@ -123,8 +123,8 @@ def convert_tcb_tdb(model, backwards=False): target_units = "TCB" if backwards else "TDB" - if model.UNITS.value == target_units or ( - model.UNITS.value is None and not backwards + if model["UNITS"].value == target_units or ( + model["UNITS"].value is None and not backwards ): log.warning("The input par file is already in the target units. Doing nothing.") return @@ -135,38 +135,24 @@ def convert_tcb_tdb(model, backwards=False): "the resulting timing model should be re-fit to get reliable results." ) - if "Spindown" in model.components: - for n, Fn_par in model.get_prefix_mapping("F").items(): - scale_parameter(model, Fn_par, n + 1, backwards) - - transform_mjd_parameter(model, "PEPOCH", backwards) - - if "AstrometryEquatorial" in model.components: - scale_parameter(model, "PMRA", 1, backwards) - scale_parameter(model, "PMDEC", 1, backwards) - elif "AstrometryEcliptic" in model.components: - scale_parameter(model, "PMELAT", 1, backwards) - scale_parameter(model, "PMELONG", 1, backwards) - transform_mjd_parameter(model, "POSEPOCH", backwards) - - # Although DM has the unit pc/cm^3, the quantity that enters - # the timing model is DMconst*DM, which has dimensions - # of frequency. Hence, DM and its derivatives will be - # scaled by IFTE_K**(i+1). - if "DispersionDM" in model.components: - scale_parameter(model, "DM", 1, backwards) - for n, DMn_par in model.get_prefix_mapping("DM").items(): - scale_parameter(model, DMn_par, n + 1, backwards) - transform_mjd_parameter(model, "DMEPOCH", backwards) - - if hasattr(model, "BINARY") and getattr(model, "BINARY").value is not None: - transform_mjd_parameter(model, "T0", backwards) - transform_mjd_parameter(model, "TASC", backwards) - scale_parameter(model, "PB", -1, backwards) - scale_parameter(model, "FB0", 1, backwards) - scale_parameter(model, "FB1", 2, backwards) - scale_parameter(model, "A1", -1, backwards) - - model.UNITS.value = target_units + for par in model.params: + param = model[par] + if ( + param.quantity is not None + and hasattr(param, "convert_tcb2tdb") + and param.convert_tcb2tdb + ): + if isinstance(param, (floatParameter, AngleParameter, maskParameter)) or ( + isinstance(param, prefixParameter) + and isinstance(param.param_comp, (floatParameter, AngleParameter)) + ): + scale_parameter(model, par, -param.effective_dimensionality, backwards) + elif isinstance(param, MJDParameter) or ( + isinstance(param, prefixParameter) + and isinstance(param.param_comp, MJDParameter) + ): + transform_mjd_parameter(model, par, backwards) + + model["UNITS"].value = target_units model.validate(allow_tcb=backwards) diff --git a/src/pint/models/timing_model.py b/src/pint/models/timing_model.py index e2d6cd298..c5216853f 100644 --- a/src/pint/models/timing_model.py +++ b/src/pint/models/timing_model.py @@ -276,14 +276,23 @@ def __init__(self, name="", components=[]): strParameter(name="UNITS", description="Units (TDB assumed)"), "" ) self.add_param_from_top( - MJDParameter(name="START", description="Start MJD for fitting"), "" + MJDParameter( + name="START", description="Start MJD for fitting", convert_tcb2tdb=False + ), + "", ) self.add_param_from_top( - MJDParameter(name="FINISH", description="End MJD for fitting"), "" + MJDParameter( + name="FINISH", description="End MJD for fitting", convert_tcb2tdb=False + ), + "", ) self.add_param_from_top( floatParameter( - name="RM", description="Rotation measure", units=u.radian / u.m**2 + name="RM", + description="Rotation measure", + units=u.radian / u.m**2, + convert_tcb2tdb=False, ), "", ) @@ -343,6 +352,7 @@ def __init__(self, name="", components=[]): name="CHI2", units="", description="Chi-squared value obtained during fitting", + convert_tcb2tdb=False, ), "", ) @@ -351,6 +361,7 @@ def __init__(self, name="", components=[]): name="CHI2R", units="", description="Reduced chi-squared value obtained during fitting", + convert_tcb2tdb=False, ), "", ) @@ -360,6 +371,7 @@ def __init__(self, name="", components=[]): name="TRES", units=u.us, description="TOA residual after fitting", + convert_tcb2tdb=False, ), "", ) @@ -368,6 +380,7 @@ def __init__(self, name="", components=[]): name="DMRES", units=u.pc / u.cm**3, description="DM residual after fitting (wideband only)", + convert_tcb2tdb=False, ), "", ) @@ -1738,6 +1751,7 @@ def jump_flags_to_params(self, toas): value=0.0, units="second", uncertainty=0.0, + tcb2tdb_scale_factor=u.Quantity(1), ) self.add_param_from_top(param, "PhaseJump") getattr(self, param.name).frozen = False @@ -2486,11 +2500,7 @@ def compare( ufloat(otherpar.value, otherpar.uncertainty.value) ) else: - # otherpar must have no uncertainty - if otherpar.value is not None: - value2[pn] = str(otherpar.value) - else: - value2[pn] = "Missing" + value2[pn] = str(otherpar.value) else: value2[pn] = "Missing" if value2[pn] == "Missing": @@ -2641,9 +2651,7 @@ def compare( and not "unc_rat" in m ): sout = colorize(sout, "red") - elif ( - "change" in m or "diff1" in m or "diff2" in m and "unc_rat" in m - ): + elif "diff2" in m: sout = colorize(sout, "red", bg_color="green") elif "unc_rat" in m: sout = colorize(sout, bg_color="green") @@ -2888,10 +2896,10 @@ def setup(self): for cp in self.components.values(): cp.setup() - def __contains__(self, name): + def __contains__(self, name: str) -> bool: return name in self.params - def __getitem__(self, name): + def __getitem__(self, name: str) -> Parameter: if name in self.top_level_params: return getattr(self, name) for cp in self.components.values(): diff --git a/src/pint/models/wave.py b/src/pint/models/wave.py index 194b2b5b8..0a7cb1dfc 100644 --- a/src/pint/models/wave.py +++ b/src/pint/models/wave.py @@ -1,4 +1,5 @@ """Delays expressed as a sum of sinusoids.""" + import astropy.units as u import numpy as np @@ -33,6 +34,7 @@ def __init__(self): name="WAVEEPOCH", description="Reference epoch for wave solution", time_scale="tdb", + convert_tcb2tdb=False, ) ) self.add_param( @@ -40,6 +42,7 @@ def __init__(self): name="WAVE_OM", description="Base frequency of wave solution", units="1/d", + convert_tcb2tdb=False, ) ) self.add_param( @@ -50,6 +53,7 @@ def __init__(self): type_match="pair", long_double=True, parameter_type="pair", + convert_tcb2tdb=False, ) ) self.phase_funcs_component += [self.wave_phase] @@ -134,6 +138,7 @@ def add_wave_component(self, amps, index=None): type_match="pair", long_double=True, parameter_type="pair", + convert_tcb2tdb=False, ) ) self.setup() diff --git a/src/pint/models/wavex.py b/src/pint/models/wavex.py index e6f85cd33..1714f3644 100644 --- a/src/pint/models/wavex.py +++ b/src/pint/models/wavex.py @@ -53,6 +53,7 @@ def __init__(self): name="WXEPOCH", description="Reference epoch for Fourier representation of red noise", time_scale="tdb", + tcb2tdb_scale_factor=u.Quantity(1), ) ) self.add_wavex_component(0.1, index=1, wxsin=0, wxcos=0, frozen=False) @@ -108,6 +109,7 @@ def add_wavex_component(self, wxfreq, index=None, wxsin=0, wxcos=0, frozen=True) units="1/d", value=wxfreq, parameter_type="float", + tcb2tdb_scale_factor=u.Quantity(1), ) ) self.add_param( @@ -118,6 +120,7 @@ def add_wavex_component(self, wxfreq, index=None, wxsin=0, wxcos=0, frozen=True) value=wxsin, frozen=frozen, parameter_type="float", + tcb2tdb_scale_factor=u.Quantity(1), ) ) self.add_param( @@ -128,6 +131,7 @@ def add_wavex_component(self, wxfreq, index=None, wxsin=0, wxcos=0, frozen=True) value=wxcos, frozen=frozen, parameter_type="float", + tcb2tdb_scale_factor=u.Quantity(1), ) ) self.setup() @@ -219,6 +223,7 @@ def add_wavex_components( units="1/d", value=wxfreq, parameter_type="float", + tcb2tdb_scale_factor=u.Quantity(1), ) ) self.add_param( @@ -229,6 +234,7 @@ def add_wavex_components( value=wxsin, parameter_type="float", frozen=frozen, + tcb2tdb_scale_factor=u.Quantity(1), ) ) self.add_param( @@ -239,6 +245,7 @@ def add_wavex_components( value=wxcos, parameter_type="float", frozen=frozen, + tcb2tdb_scale_factor=u.Quantity(1), ) ) self.setup() diff --git a/src/pint/scripts/tcb2tdb.py b/src/pint/scripts/tcb2tdb.py index 505d74d4e..4c427d21a 100644 --- a/src/pint/scripts/tcb2tdb.py +++ b/src/pint/scripts/tcb2tdb.py @@ -17,26 +17,14 @@ def main(argv=None): description="""`tcb2tdb` converts TCB par files to TDB. Please note that this conversion is not exact and the timing model should be re-fit to the TOAs. - - The following parameters are converted to TDB: - 1. Spin frequency, its derivatives and spin epoch - 2. Sky coordinates, proper motion and the position epoch - 3. DM, DM derivatives and DM epoch - 4. Keplerian binary parameters and FB1 - + The following parameters are NOT converted although they are in fact affected by the TCB to TDB conversion: - 1. Parallax - 2. TZRMJD and TZRFRQ - 3. DMX parameters - 4. Solar wind parameters - 5. Binary post-Keplerian parameters including Shapiro delay - parameters (except FB1) - 6. Jumps and DM Jumps - 7. FD parameters - 8. EQUADs - 9. Red noise parameters including FITWAVES, powerlaw red noise and - powerlaw DM noise parameters + 1. TZRMJD and TZRFRQ + 2. DM Jumps (the wideband kind) + 3. FD parameters and FD jumps + 4. EQUADs and ECORRs + 5. GP Red noise parameters and GP DM noise parameters """, formatter_class=argparse.ArgumentDefaultsHelpFormatter, ) diff --git a/src/pint/utils.py b/src/pint/utils.py index c0723cb17..f48b9f05b 100644 --- a/src/pint/utils.py +++ b/src/pint/utils.py @@ -1979,8 +1979,8 @@ def ELL1_check( Checks whether the assumptions that allow ELL1 to be safely used are satisfied. To work properly, we should have: - :math:`asini/c e^4 \ll {\\rm timing precision} / \sqrt N_{\\rm TOA}` - or :math:`A1 E^4 \ll TRES / \sqrt N_{\\rm TOA}` + :math:`asini/c e^4 \\ll {\\rm timing precision} / \\sqrt N_{\\rm TOA}` + or :math:`A1 E^4 \\ll TRES / \\sqrt N_{\\rm TOA}` since the ELL1 model now includes terms up to O(E^3) diff --git a/tests/test_all_component_and_model_builder.py b/tests/test_all_component_and_model_builder.py index 997ca61a3..8dab2f33a 100644 --- a/tests/test_all_component_and_model_builder.py +++ b/tests/test_all_component_and_model_builder.py @@ -17,6 +17,7 @@ from pint.models.parameter import floatParameter from pint.utils import split_prefixed_name, PrefixError from pinttestdata import datadir +from astropy import units as u class SimpleModel(PhaseComponent): @@ -27,7 +28,14 @@ class SimpleModel(PhaseComponent): def __init__(self): super(SimpleModel, self).__init__() - self.add_param(floatParameter(name="TESTPARAM", value=0.0, unit="s")) + self.add_param( + floatParameter( + name="TESTPARAM", + value=0.0, + unit="s", + tcb2tdb_scale_factor=u.Quantity(1), + ) + ) class SubsetModel(PhaseComponent): @@ -38,7 +46,11 @@ class SubsetModel(PhaseComponent): def __init__(self): super(SubsetModel, self).__init__() - self.add_param(floatParameter(name="F0", value=0.0, unit="1/s")) + self.add_param( + floatParameter( + name="F0", value=0.0, unit="1/s", tcb2tdb_scale_factor=u.Quantity(1) + ) + ) # self.add_param(floatParameter(name="F1", value=0.0, unit="1/s^2")) @@ -54,7 +66,15 @@ def simple_model(): @pytest.fixture def simple_model_overlap(simple_model): - simple_model.add_param(floatParameter(name="F0", aliases=[], value=0.0, unit="1/s")) + simple_model.add_param( + floatParameter( + name="F0", + aliases=[], + value=0.0, + unit="1/s", + tcb2tdb_scale_factor=u.Quantity(1), + ) + ) return simple_model @@ -62,10 +82,22 @@ def simple_model_overlap(simple_model): def simple_model_alias_overlap(): simple_model = SimpleModel() simple_model.add_param( - floatParameter(name="TESTPARAM2", aliases=["F0"], value=0.0, unit="1/s") + floatParameter( + name="TESTPARAM2", + aliases=["F0"], + value=0.0, + unit="1/s", + tcb2tdb_scale_factor=u.Quantity(1), + ) ) simple_model.add_param( - floatParameter(name="TESTPARAM3", aliases=["LAMBDA"], value=0.0, unit="deg") + floatParameter( + name="TESTPARAM3", + aliases=["LAMBDA"], + value=0.0, + unit="deg", + tcb2tdb_scale_factor=u.Quantity(1), + ) ) return simple_model @@ -108,7 +140,13 @@ def test_model_builder_class(): assert "SimpleModel" in mb.components simple_comp = mb.components["SimpleModel"] simple_comp.add_param( - floatParameter(name="TESTPARAM2", aliases=["F0"], value=0.0, unit="s") + floatParameter( + name="TESTPARAM2", + aliases=["F0"], + value=0.0, + unit="s", + tcb2tdb_scale_factor=u.Quantity(1), + ) ) @@ -173,7 +211,13 @@ class SimpleModel2(PhaseComponent): def __init__(self): super(SimpleModel2, self).__init__() self.add_param( - floatParameter(name="TESTPARAMF0", aliases=["F0"], value=0.0, unit="s") + floatParameter( + name="TESTPARAMF0", + aliases=["F0"], + value=0.0, + unit="s", + tcb2tdb_scale_factor=u.Quantity(1), + ) ) mb2 = AllComponents() @@ -392,7 +436,11 @@ class SubsetModel2(PhaseComponent): def __init__(self): super(SubsetModel2, self).__init__() - self.add_param(floatParameter(name="F0", value=0.0, unit="1/s")) + self.add_param( + floatParameter( + name="F0", value=0.0, unit="1/s", tcb2tdb_scale_factor=u.Quantity(1) + ) + ) with pytest.raises(ComponentConflict): mb = ModelBuilder() diff --git a/tests/test_fitter.py b/tests/test_fitter.py index dd1409dc2..b22a12ade 100644 --- a/tests/test_fitter.py +++ b/tests/test_fitter.py @@ -103,6 +103,7 @@ def test_ftest_nb(): value=0.0, units=u.Hz / u.s / u.s, frozen=False, + tcb2tdb_scale_factor=u.Quantity(1), ) ft = f.ftest(F2, "Spindown", remove=False) assert isinstance(ft["ft"], (float, bool)) @@ -111,7 +112,12 @@ def test_ftest_nb(): assert isinstance(Ftest_dict["ft"], (float, bool)) # Test removing parameter F1 = param.prefixParameter( - parameter_type="float", name="F1", value=0.0, units=u.Hz / u.s, frozen=False + parameter_type="float", + name="F1", + value=0.0, + units=u.Hz / u.s, + frozen=False, + tcb2tdb_scale_factor=u.Quantity(1), ) ft = f.ftest(F1, "Spindown", remove=True) assert isinstance(ft["ft"], (float, bool)) @@ -129,7 +135,12 @@ def test_ftest_wb(): wb_f.fit_toas() # Parallax PX = param.floatParameter( - parameter_type="float", name="PX", value=0.0, units=u.mas, frozen=False + parameter_type="float", + name="PX", + value=0.0, + units=u.mas, + frozen=False, + tcb2tdb_scale_factor=u.Quantity(1), ) PX_Component = "AstrometryEcliptic" # A1DOT @@ -139,6 +150,7 @@ def test_ftest_wb(): value=0.0, units=ls / u.second, frozen=False, + tcb2tdb_scale_factor=u.Quantity(1), ) A1DOT_Component = "BinaryELL1" # Test adding A1DOT diff --git a/tests/test_flagging_clustering.py b/tests/test_flagging_clustering.py index 21f0e426c..c2584b40b 100644 --- a/tests/test_flagging_clustering.py +++ b/tests/test_flagging_clustering.py @@ -1,4 +1,5 @@ """Tests for clustering and flagging""" + import os import pytest import copy @@ -38,7 +39,12 @@ def test_jump_by_cluster(setup_NGC6440E): setup_NGC6440E.m.add_component(PhaseJump(), validate=False) cp = setup_NGC6440E.m.components["PhaseJump"] par = p.maskParameter( - name="JUMP", key="mjd", value=0.2, key_value=[54099, 54100], units=u.s + name="JUMP", + key="mjd", + value=0.2, + key_value=[54099, 54100], + units=u.s, + tcb2tdb_scale_factor=u.Quantity(1), ) # this should be the last group of any TOAs cp.add_param(par, setup=True) @@ -51,7 +57,12 @@ def test_jump_by_cluster(setup_NGC6440E): m_copy.add_component(PhaseJump(), validate=False) cp_copy = m_copy.components["PhaseJump"] par_copy = p.maskParameter( - name="JUMP", key="-cluster", value=0.2, key_value=41, units=u.s + name="JUMP", + key="-cluster", + value=0.2, + key_value=41, + units=u.s, + tcb2tdb_scale_factor=u.Quantity(1), ) # this should be identical to the cluster above cp_copy.add_param(par_copy, setup=True) diff --git a/tests/test_grid.py b/tests/test_grid.py index a3ffbc652..eb12dc28b 100644 --- a/tests/test_grid.py +++ b/tests/test_grid.py @@ -1,4 +1,5 @@ """Test chi^2 gridding routines""" + import concurrent.futures import pytest @@ -296,6 +297,7 @@ def test_grid_3param_prefix_singleprocessor(): description=modelcomponent.F_description(2), longdouble=True, frozen=False, + tcb2tdb_scale_factor=u.Quantity(1), ) modelcomponent.add_param(p, setup=True) m.validate() @@ -340,6 +342,7 @@ def test_grid_3param_prefix_multiprocessor(): description=modelcomponent.F_description(2), longdouble=True, frozen=False, + tcb2tdb_scale_factor=u.Quantity(1), ) modelcomponent.add_param(p, setup=True) m.validate() diff --git a/tests/test_jump.py b/tests/test_jump.py index be53697fc..e5125e5b4 100644 --- a/tests/test_jump.py +++ b/tests/test_jump.py @@ -1,4 +1,5 @@ """Tests for jump model component """ + import logging import os import pytest @@ -100,7 +101,12 @@ def test_jump_params_to_flags(setup_NGC6440E): cp = setup_NGC6440E.m.components["PhaseJump"] par = p.maskParameter( - name="JUMP", key="freq", value=0.2, key_value=[1440, 1700], units=u.s + name="JUMP", + key="freq", + value=0.2, + key_value=[1440, 1700], + units=u.s, + tcb2tdb_scale_factor=u.Quantity(1), ) # TOAs indexed 48, 49, 54 in NGC6440E are within this frequency range cp.add_param(par, setup=True) @@ -128,7 +134,12 @@ def test_jump_params_to_flags(setup_NGC6440E): # check that adding overlapping jump works par2 = p.maskParameter( - name="JUMP", key="freq", value=0.2, key_value=[1600, 1900], units=u.s + name="JUMP", + key="freq", + value=0.2, + key_value=[1600, 1900], + units=u.s, + tcb2tdb_scale_factor=u.Quantity(1), ) # frequency range overlaps with par, 2nd jump will have common TOAs w/ 1st cp.add_param(par2, setup=True) # add flags based off jumps added to model @@ -146,7 +157,12 @@ def test_multijump_toa(setup_NGC6440E): setup_NGC6440E.m.add_component(PhaseJump(), validate=False) cp = setup_NGC6440E.m.components["PhaseJump"] par = p.maskParameter( - name="JUMP", key="freq", value=0.2, key_value=[1440, 1700], units=u.s + name="JUMP", + key="freq", + value=0.2, + key_value=[1440, 1700], + units=u.s, + tcb2tdb_scale_factor=u.Quantity(1), ) # TOAs indexed 48, 49, 54 in NGC6440E are within this frequency range selected_toa_ind = [48, 49, 54] cp.add_param(par, setup=True) @@ -171,7 +187,12 @@ def test_unfrozen_jump(setup_NGC6440E): setup_NGC6440E.m.add_component(PhaseJump(), validate=False) # this has no TOAs par = p.maskParameter( - name="JUMP", key="freq", value=0.2, key_value=[3000, 3200], units=u.s + name="JUMP", + key="freq", + value=0.2, + key_value=[3000, 3200], + units=u.s, + tcb2tdb_scale_factor=u.Quantity(1), ) setup_NGC6440E.m.components["PhaseJump"].add_param(par, setup=True) setup_NGC6440E.m.JUMP1.frozen = False @@ -183,7 +204,12 @@ def test_find_empty_masks(setup_NGC6440E): setup_NGC6440E.m.add_component(PhaseJump(), validate=False) # this has no TOAs par = p.maskParameter( - name="JUMP", key="freq", value=0.2, key_value=[3000, 3200], units=u.s + name="JUMP", + key="freq", + value=0.2, + key_value=[3000, 3200], + units=u.s, + tcb2tdb_scale_factor=u.Quantity(1), ) setup_NGC6440E.m.components["PhaseJump"].add_param(par, setup=True) setup_NGC6440E.m.JUMP1.frozen = False diff --git a/tests/test_mask_parameter.py b/tests/test_mask_parameter.py index ffb532c7c..245c7264f 100644 --- a/tests/test_mask_parameter.py +++ b/tests/test_mask_parameter.py @@ -20,7 +20,9 @@ def toas(): def test_mjd_mask(toas): - mp = maskParameter("test1", key="mjd", key_value=[54000, 54100]) + mp = maskParameter( + "test1", key="mjd", key_value=[54000, 54100], tcb2tdb_scale_factor=u.Quantity(1) + ) assert mp.key == "mjd" assert mp.key_value == [54000, 54100] assert mp.value is None @@ -32,16 +34,27 @@ def test_mjd_mask(toas): assert np.all(select_toas == raw_selection[0]) assert np.all(toas.table["mjd_float"][select_toas] <= 54100) assert np.all(toas.table["mjd_float"][select_toas] >= 54000) - mp_str_keyval = maskParameter("test2", key="mjd", key_value=["54000", "54100"]) + mp_str_keyval = maskParameter( + "test2", + key="mjd", + key_value=["54000", "54100"], + tcb2tdb_scale_factor=u.Quantity(1), + ) assert mp_str_keyval.key_value == [54000, 54100] - mp_value_switch = maskParameter("test1", key="mjd", key_value=[54100, 54000]) + mp_value_switch = maskParameter( + "test1", key="mjd", key_value=[54100, 54000], tcb2tdb_scale_factor=u.Quantity(1) + ) assert mp_value_switch.key_value == [54000, 54100] with pytest.raises(ValueError): - mp_str_keyval = maskParameter("test2", key="mjd", key_value=["54000"]) + mp_str_keyval = maskParameter( + "test2", key="mjd", key_value=["54000"], tcb2tdb_scale_factor=u.Quantity(1) + ) def test_freq_mask(toas): - mp = maskParameter("test2", key="freq", key_value=[1400, 1430]) + mp = maskParameter( + "test2", key="freq", key_value=[1400, 1430], tcb2tdb_scale_factor=u.Quantity(1) + ) assert mp.key_value == [1400 * u.MHz, 1430 * u.MHz] select_toas = mp.select_toa_mask(toas) assert len(select_toas) > 0 @@ -49,41 +62,69 @@ def test_freq_mask(toas): (toas.table["freq"] <= 1430 * u.MHz) * (toas.table["freq"] >= 1400 * u.MHz) ) assert np.all(select_toas == raw_selection[0]) - mp_str = maskParameter("test2", key="freq", key_value=["1400", "2000"]) + mp_str = maskParameter( + "test2", + key="freq", + key_value=["1400", "2000"], + tcb2tdb_scale_factor=u.Quantity(1), + ) assert mp_str.key_value == [1400 * u.MHz, 2000 * u.MHz] - mp_switch = maskParameter("test2", key="freq", key_value=[2000, 1400]) + mp_switch = maskParameter( + "test2", key="freq", key_value=[2000, 1400], tcb2tdb_scale_factor=u.Quantity(1) + ) assert mp_switch.key_value == [1400 * u.MHz, 2000 * u.MHz] mp_quantity = maskParameter( - "test2", key="freq", key_value=[2000 * u.MHz, 1400 * u.MHz] + "test2", + key="freq", + key_value=[2000 * u.MHz, 1400 * u.MHz], + tcb2tdb_scale_factor=u.Quantity(1), ) assert mp_quantity.key_value == [1400 * u.MHz, 2000 * u.MHz] with pytest.raises(ValueError): - mp_wrong_keyvalue = maskParameter("test2", key="freq", key_value=[1400]) + mp_wrong_keyvalue = maskParameter( + "test2", key="freq", key_value=[1400], tcb2tdb_scale_factor=u.Quantity(1) + ) def test_tel_mask(toas): - mp_ao = maskParameter("test2", key="tel", key_value="ao") + mp_ao = maskParameter( + "test2", key="tel", key_value="ao", tcb2tdb_scale_factor=u.Quantity(1) + ) assert mp_ao.key_value == ["arecibo"] select_toas = mp_ao.select_toa_mask(toas) assert len(select_toas) > 0 raw_selection = np.where(toas.table["obs"] == "arecibo") assert np.all(select_toas == raw_selection[0]) - mp_gbt = maskParameter("test2", key="tel", key_value=["gbt"]) + mp_gbt = maskParameter( + "test2", key="tel", key_value=["gbt"], tcb2tdb_scale_factor=u.Quantity(1) + ) assert mp_gbt.key_value == ["gbt"] select_toas = mp_gbt.select_toa_mask(toas) assert np.all(toas.table["obs"][select_toas] == "gbt") - mp_special_obs1 = maskParameter("test2", key="tel", key_value="@") + mp_special_obs1 = maskParameter( + "test2", key="tel", key_value="@", tcb2tdb_scale_factor=u.Quantity(1) + ) assert mp_special_obs1.key_value == ["barycenter"] - mp_special_obs2 = maskParameter("test2", key="tel", key_value=0) + mp_special_obs2 = maskParameter( + "test2", key="tel", key_value=0, tcb2tdb_scale_factor=u.Quantity(1) + ) assert mp_special_obs2.key_value == ["geocenter"] with pytest.raises(ValueError): - mp_wrong_keyvalue = maskParameter("test2", key="tel", key_value=["gbt", "ao"]) + mp_wrong_keyvalue = maskParameter( + "test2", + key="tel", + key_value=["gbt", "ao"], + tcb2tdb_scale_factor=u.Quantity(1), + ) def test_name_mask(toas): # Not sure about the use case for name mask mp_name = maskParameter( - "test2", key="name", key_value="53393.000009.3.000.000.9y.x.ff" + "test2", + key="name", + key_value="53393.000009.3.000.000.9y.x.ff", + tcb2tdb_scale_factor=u.Quantity(1), ) assert mp_name.key_value == ["53393.000009.3.000.000.9y.x.ff"] select_toas = mp_name.select_toa_mask(toas) @@ -92,18 +133,29 @@ def test_name_mask(toas): assert np.all(select_toas == raw_selection[0]) with pytest.raises(ValueError): mp_wrong_keyvalue = maskParameter( - "test2", key="name", key_value=["name1", "name2"] + "test2", + key="name", + key_value=["name1", "name2"], + tcb2tdb_scale_factor=u.Quantity(1), ) def test_flag_mask(toas): - mp_flag = maskParameter("test2", key="-fe", key_value=430) + mp_flag = maskParameter( + "test2", key="-fe", key_value=430, tcb2tdb_scale_factor=u.Quantity(1) + ) assert mp_flag.key_value == ["430"] - mp_flag2 = maskParameter("test2", key="-fe", key_value="430") + mp_flag2 = maskParameter( + "test2", key="-fe", key_value="430", tcb2tdb_scale_factor=u.Quantity(1) + ) assert mp_flag2.key_value == ["430"] with pytest.raises(ValueError): - mp_wrong_key = maskParameter("test2", key="fe", key_value="430") - mp_flag3 = maskParameter("test2", key="-fe", key_value="L-wide") + mp_wrong_key = maskParameter( + "test2", key="fe", key_value="430", tcb2tdb_scale_factor=u.Quantity(1) + ) + mp_flag3 = maskParameter( + "test2", key="-fe", key_value="L-wide", tcb2tdb_scale_factor=u.Quantity(1) + ) assert mp_flag3.key_value == ["L-wide"] select_toas = mp_flag3.select_toa_mask(toas) assert len(select_toas) > 0 diff --git a/tests/test_parameters.py b/tests/test_parameters.py index 460bd1ca8..a232a7b59 100644 --- a/tests/test_parameters.py +++ b/tests/test_parameters.py @@ -476,7 +476,7 @@ def test_start_finish_notfrozen(self): ], ) def test_parameter_parsing(type_, str_value, value): - p = type_(name="TEST") + p = type_(name="TEST", tcb2tdb_scale_factor=u.Quantity(1)) assert p.from_parfile_line(str_value) assert p.value == value @@ -494,7 +494,7 @@ def test_parameter_parsing(type_, str_value, value): ], ) def test_parameter_parsing_fail(type_, par_line): - p = type_(name="TEST") + p = type_(name="TEST", tcb2tdb_scale_factor=u.Quantity(1)) with pytest.raises(ValueError): p.from_parfile_line(par_line) @@ -514,7 +514,7 @@ def test_parameter_parsing_fail(type_, par_line): ], ) def test_parameter_setting(type_, set_value, value): - p = type_(name="TEST") + p = type_(name="TEST", tcb2tdb_scale_factor=u.Quantity(1)) if isinstance(value, type(Exception)) and issubclass(value, Exception): with pytest.raises(value): p.value = set_value @@ -533,10 +533,18 @@ def test_set_uncertainty_bogus_raises(p): def test_compare_key_value(): - par1 = maskParameter(name="Test1", key="-k1", key_value="kv1") - par2 = maskParameter(name="Test2", key="-k1", key_value="kv1") - par3 = maskParameter(name="Test3", key="-k2", key_value="kv1") - par4 = maskParameter(name="Test4", key="-k1", key_value="kv2") + par1 = maskParameter( + name="Test1", key="-k1", key_value="kv1", tcb2tdb_scale_factor=u.Quantity(1) + ) + par2 = maskParameter( + name="Test2", key="-k1", key_value="kv1", tcb2tdb_scale_factor=u.Quantity(1) + ) + par3 = maskParameter( + name="Test3", key="-k2", key_value="kv1", tcb2tdb_scale_factor=u.Quantity(1) + ) + par4 = maskParameter( + name="Test4", key="-k1", key_value="kv2", tcb2tdb_scale_factor=u.Quantity(1) + ) assert par1.compare_key_value(par2) assert par2.compare_key_value(par1) assert not par1.compare_key_value(par3) @@ -546,18 +554,53 @@ def test_compare_key_value(): def test_compare_key_value_list(): - par1 = maskParameter(name="Test1", key="-k1", key_value=["k", "q"]) - par2 = maskParameter(name="Test2", key="-k1", key_value=["q", "k"]) - par3 = maskParameter(name="Test3", key="-k1", key_value=["k", "q"]) + par1 = maskParameter( + name="Test1", + key="-k1", + key_value=["k", "q"], + tcb2tdb_scale_factor=u.Quantity(1), + ) + par2 = maskParameter( + name="Test2", + key="-k1", + key_value=["q", "k"], + tcb2tdb_scale_factor=u.Quantity(1), + ) + par3 = maskParameter( + name="Test3", + key="-k1", + key_value=["k", "q"], + tcb2tdb_scale_factor=u.Quantity(1), + ) assert par1.compare_key_value(par2) assert par2.compare_key_value(par1) assert par3.compare_key_value(par1) - par4 = maskParameter(name="Test4", key="-k1", key_value=[2000, 3000]) - par5 = maskParameter(name="Test5", key="-k1", key_value=[3000, 2000]) + par4 = maskParameter( + name="Test4", + key="-k1", + key_value=[2000, 3000], + tcb2tdb_scale_factor=u.Quantity(1), + ) + par5 = maskParameter( + name="Test5", + key="-k1", + key_value=[3000, 2000], + tcb2tdb_scale_factor=u.Quantity(1), + ) assert par4.compare_key_value(par5) assert par5.compare_key_value(par4) - par6 = maskParameter(name="Test6", key="-k1", key_value=["430", "guppi", "puppi"]) - par7 = maskParameter(name="Test7", key="-k1", key_value=["puppi", "guppi", "430"]) + par6 = maskParameter( + name="Test6", + key="-k1", + key_value=["430", "guppi", "puppi"], + tcb2tdb_scale_factor=u.Quantity(1), + ) + par7 = maskParameter( + name="Test7", + key="-k1", + key_value=["puppi", "guppi", "430"], + tcb2tdb_scale_factor=u.Quantity(1), + ) assert par6.compare_key_value(par7) assert par7.compare_key_value(par6) @@ -569,13 +612,13 @@ def test_compare_key_value_list(): intParameter(name="FISH"), strParameter(name="FISH"), pytest.param( - maskParameter(name="JUMP"), + maskParameter(name="JUMP", tcb2tdb_scale_factor=u.Quantity(1)), ), pytest.param( - prefixParameter(name="F0"), + prefixParameter(name="F0", tcb2tdb_scale_factor=u.Quantity(1)), ), - pairParameter(name="WEAVE"), - AngleParameter(name="BEND"), + pairParameter(name="WEAVE", convert_tcb2tdb=False), + AngleParameter(name="BEND", tcb2tdb_scale_factor=u.Quantity(1)), ], ) def test_parameter_can_be_pickled(p): @@ -586,7 +629,9 @@ def test_fitter_construction_success_after_remove_param(): """Checks that add_param and remove_param don't require m.setup() to be run prior to constructing a fitter. This addresses issue #1260.""" m = get_model(os.path.join(datadir, "B1855+09_NANOGrav_9yv1.gls.par")) t = get_TOAs(os.path.join(datadir, "B1855+09_NANOGrav_9yv1.tim")) - FD4 = prefixParameter(parameter_type="float", name="FD4", value=0.0, units=u.s) + FD4 = prefixParameter( + parameter_type="float", name="FD4", value=0.0, units=u.s, convert_tcb2tdb=False + ) """Fitter construction used to fail after remove_param without m.setup(). Test this (for both adding and removing, just to be safe):""" m.add_param_from_top(FD4, "FD") f = pint.fitter.GLSFitter(toas=t, model=m) @@ -598,7 +643,9 @@ def test_correct_number_of_params_and_FD_terms_after_add_or_remove_param(): """Checks that the number of parameters in some component after add_param or remove_param is correct. Similarly checks that len(m.get_prefix_mapping('FD')) (i.e. num_FD_terms) is correct after seeing a bug where the length didn't change after FD param removal (see #1260).""" m = get_model(os.path.join(datadir, "B1855+09_NANOGrav_9yv1.gls.par")) t = get_TOAs(os.path.join(datadir, "B1855+09_NANOGrav_9yv1.tim")) - FD4 = prefixParameter(parameter_type="float", name="FD4", value=0.0, units=u.s) + FD4 = prefixParameter( + parameter_type="float", name="FD4", value=0.0, units=u.s, convert_tcb2tdb=False + ) m.add_param_from_top(FD4, "FD") assert len(m.components["FD"].params) == 4 assert len(m.get_prefix_mapping("FD")) == 4 diff --git a/tests/test_pickle.py b/tests/test_pickle.py index 90f4809ff..4e4302663 100644 --- a/tests/test_pickle.py +++ b/tests/test_pickle.py @@ -7,6 +7,7 @@ from pint.models.model_builder import get_model from pinttestdata import datadir import pickle +from astropy import units as u def test_pickle_prefixparameter(): @@ -23,6 +24,7 @@ def test_pickle_prefixparameter(): description=modelcomponent.F_description(2), longdouble=True, frozen=False, + tcb2tdb_scale_factor=u.Quantity(1), ) modelcomponent.add_param(p, setup=True) m.validate() diff --git a/tests/test_solarsystemshapiro.py b/tests/test_solarsystemshapiro.py new file mode 100644 index 000000000..0ef610479 --- /dev/null +++ b/tests/test_solarsystemshapiro.py @@ -0,0 +1,11 @@ +from pint.models import get_model_and_toas +from pinttestdata import datadir +import pytest + + +def test_solarsystemshapiro_exception(): + m, t = get_model_and_toas(datadir / "NGC6440E.par", datadir / "NGC6440E.tim") + m.PLANET_SHAPIRO.value = True + + with pytest.raises(KeyError): + m.components["SolarSystemShapiro"].solar_system_shapiro_delay(t) diff --git a/tests/test_tcb2tdb.py b/tests/test_tcb2tdb.py index 34e7373a1..7c34d9a93 100644 --- a/tests/test_tcb2tdb.py +++ b/tests/test_tcb2tdb.py @@ -7,7 +7,8 @@ import numpy as np import pytest -from pint.models.model_builder import ModelBuilder +from pint import DMconst, dmu +from pint.models.model_builder import ModelBuilder, get_model from pint.models.tcb_conversion import convert_tcb_tdb from pint.scripts import tcb2tdb @@ -27,6 +28,7 @@ ECC 1.0 OM 0.0 PB 10.0 +FD1 1e-3 EPHEM DE436 CLK TT(BIPM2017) UNITS TCB @@ -68,6 +70,60 @@ def test_convert_units_roundtrip(): assert np.isclose(getattr(m, par).value, getattr(m_, par).value) +def test_effective_dimensionality(): + m = ModelBuilder()(StringIO(simplepar), allow_tcb=True) + assert m.PEPOCH.effective_dimensionality == 1 + assert m.F0.effective_dimensionality == -1 + assert m.F1.effective_dimensionality == -2 + + assert m.POSEPOCH.effective_dimensionality == 1 + assert m.RAJ.effective_dimensionality == 0 + assert m.DECJ.effective_dimensionality == 0 + assert m.PMRA.effective_dimensionality == -1 + assert m.PMDEC.effective_dimensionality == -1 + assert m.PX.effective_dimensionality == -1 + + assert m.DM.effective_dimensionality == -1 + + assert m.T0.effective_dimensionality == 1 + assert m.A1.effective_dimensionality == 1 + assert m.ECC.effective_dimensionality == 0 + assert m.OM.effective_dimensionality == 0 + assert m.PB.effective_dimensionality == 1 + + assert m.NE_SW.effective_dimensionality == -2 + + +def test_dm_scaling_factor(): + m = get_model( + StringIO( + """ + PSR TEST + F0 100 + F1 -1e-14 + PEPOCH 55000 + DMEPOCH 55000 + DM 12.5 + DM1 -0.001 + DM2 1e-5 + DMXR1_0001 51000 + DMXR2_0001 51000 + DMX_0001 0.002 + DMWXEPOCH 55000 + DMWXFREQ_0001 0.001 + DMWXSIN_0001 0.0003 + DMWXCOS_0001 0.0002 + """ + ) + ) + + for param in m.params: + par = m[param] + + if hasattr(par, "units") and par.units == dmu: + assert not par.convert_tcb2tdb or par.tcb2tdb_scale_factor == DMconst + + def test_tcb2tdb(tmp_path): tmppar1 = tmp_path / "tmp1.par" tmppar2 = tmp_path / "tmp2.par" diff --git a/tests/test_timing_model.py b/tests/test_timing_model.py index bef73fac1..2cd37645c 100644 --- a/tests/test_timing_model.py +++ b/tests/test_timing_model.py @@ -56,15 +56,9 @@ def test_from_par(self): assert len(tm.DelayComponent_list) == 4 assert len(tm.PhaseComponent_list) == 2 - # Check delay component order - order = [] - for dcp in tm.DelayComponent_list: - order.append(DEFAULT_ORDER.index(dcp.category)) + order = [DEFAULT_ORDER.index(dcp.category) for dcp in tm.DelayComponent_list] assert all(np.diff(np.array(order)) > 0) - # Check phase component order - order = [] - for dcp in tm.PhaseComponent_list: - order.append(DEFAULT_ORDER.index(dcp.category)) + order = [DEFAULT_ORDER.index(dcp.category) for dcp in tm.PhaseComponent_list] assert all(np.diff(np.array(order)) > 0) def test_component_input(self): @@ -77,16 +71,10 @@ def test_component_input(self): # test the link to timing model assert v._parent == tm - # Test Delay order - order = [] - for dcp in tm.DelayComponent_list: - order.append(DEFAULT_ORDER.index(dcp.category)) + order = [DEFAULT_ORDER.index(dcp.category) for dcp in tm.DelayComponent_list] assert all(np.diff(np.array(order)) > 0) - # Test Phase order - order = [] - for dcp in tm.PhaseComponent_list: - order.append(DEFAULT_ORDER.index(dcp.category)) + order = [DEFAULT_ORDER.index(dcp.category) for dcp in tm.PhaseComponent_list] assert all(np.diff(np.array(order)) > 0) def test_add_component(self): @@ -123,6 +111,7 @@ def test_add_component(self): value=p_vals["value"], key_value=p_vals["key_value"], units=u.s, + tcb2tdb_scale_factor=u.Quantity(1), ) print("test", par.name) cp.add_param(par, setup=True) @@ -268,7 +257,7 @@ def test_items(model_0437): def test_iterator(model_0437): - assert [k for k in model_0437] == model_0437.params + assert list(model_0437) == model_0437.params par_base = """