diff --git a/CHANGELOG.md b/CHANGELOG.md index 755b0d86b..429267750 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -36,6 +36,7 @@ and this project, at least loosely, adheres to [Semantic Versioning](https://sem - Module for converting between binary models also included in `convert_parfile` - Method to get a parameter as a `uncertainties.ufloat` for doing math - Method to get current binary period and uncertainty at a given time regardless of binary model +- TCB to TDB conversion on read, and conversion script (`tcb2tdb`) ### Fixed - Broken notebooks CI test - BIPM correction for simulated TOAs diff --git a/docs/command-line.rst b/docs/command-line.rst index 04ba8ef76..fa04b56c8 100644 --- a/docs/command-line.rst +++ b/docs/command-line.rst @@ -111,3 +111,11 @@ the examples subdirectory of the PINT distro. event_optimize J0030+0451_P8_15.0deg_239557517_458611204_ft1weights_GEO_wt.gt.0.4.fits PSRJ0030+0451_psrcat.par templateJ0030.3gauss --weightcol=PSRJ0030+0451 --minWeight=0.9 --nwalkers=100 --nsteps=500 +tcb2tdb +------- + +A command line tool that converts par files from TCB timescale to TDB timescale. + +:: + + tcb2tdb J0030+0451_tcb.par J0030+0451_tdb.par \ No newline at end of file diff --git a/docs/explanation.rst b/docs/explanation.rst index 35e37ac54..a005c9bc5 100644 --- a/docs/explanation.rst +++ b/docs/explanation.rst @@ -147,12 +147,14 @@ in, and what kind of time you're asking for:: The conventional time scale for working with pulsars, and the one PINT uses, is Barycentric Dynamical Time (TDB). You should be aware that there -is another time scale, not yet supported in PINT, called Barycentric -Coordinate Time (TCB), and that because of different handling of -relativistic corrections, it does not advance at the same rate as TDB +is another time scale, not yet fully supported in PINT, called Barycentric +Coordinate Time (TCB). Because of different handling of relativistic +corrections, the TCB timescale does not advance at the same rate as TDB (there is also a many-second offset). TEMPO2 uses TCB by default, so you may encounter pulsar timing models or even measurements that use -TCB. PINT will attempt to detect this and let you know. +TCB. PINT provides a command line tool `tcb2tdb` to approximately convert +TCB timing models to TDB. PINT can also optionally convert TCB timing models +to TDB (approximately) upon read. Note that the need for leap seconds is because the Earth's rotation is somewhat erratic - no, we're not about to be thrown off, but its diff --git a/setup.cfg b/setup.cfg index c62af7a06..80a0c6408 100644 --- a/setup.cfg +++ b/setup.cfg @@ -60,6 +60,7 @@ console_scripts = pintk = pint.scripts.pintk:main convert_parfile = pint.scripts.convert_parfile:main compare_parfiles = pint.scripts.compare_parfiles:main + tcb2tdb = pint.scripts.tcb2tdb:main # See the docstring in versioneer.py for instructions. Note that you must diff --git a/src/pint/models/model_builder.py b/src/pint/models/model_builder.py index a90383f85..2a14bd73c 100644 --- a/src/pint/models/model_builder.py +++ b/src/pint/models/model_builder.py @@ -1,3 +1,5 @@ +"""Building a timing model from a par file.""" + import copy import warnings from io import StringIO @@ -21,7 +23,7 @@ ) from pint.toa import get_TOAs from pint.utils import PrefixError, interesting_lines, lines_of, split_prefixed_name - +from pint.models.tcb_conversion import convert_tcb_tdb __all__ = ["ModelBuilder", "get_model", "get_model_and_toas"] @@ -73,7 +75,7 @@ def __init__(self): self._validate_components() self.default_components = ["SolarSystemShapiro"] - def __call__(self, parfile, allow_name_mixing=False): + def __call__(self, parfile, allow_name_mixing=False, allow_tcb=False): """Callable object for making a timing model from .par file. Parameters @@ -87,11 +89,22 @@ def __call__(self, parfile, allow_name_mixing=False): T2EFAC and EFAC, both of them maps to PINT parameter EFAC, present in the parfile at the same time. + allow_tcb : True, False, or "raw", optional + Whether to read TCB par files. Default is False, and will throw an + error upon encountering TCB par files. If True, the par file will be + converted to TDB upon read. If "raw", an unconverted malformed TCB + TimingModel object will be returned. + Returns ------- pint.models.timing_model.TimingModel The result timing model based on the input .parfile or file object. """ + + assert allow_tcb in [True, False, "raw"] + convert_tcb = allow_tcb == True + allow_tcb_ = allow_tcb in [True, "raw"] + pint_param_dict, original_name, unknown_param = self._pintify_parfile( parfile, allow_name_mixing ) @@ -103,12 +116,23 @@ def __call__(self, parfile, allow_name_mixing=False): # Make timing model cps = [self.all_components.components[c] for c in selected] tm = TimingModel(components=cps) - self._setup_model(tm, pint_param_dict, original_name, setup=True, validate=True) + self._setup_model( + tm, + pint_param_dict, + original_name, + setup=True, + validate=True, + allow_tcb=allow_tcb_, + ) # Report unknown line for k, v in unknown_param.items(): p_line = " ".join([k] + v) warnings.warn(f"Unrecognized parfile line '{p_line}'", UserWarning) # log.warning(f"Unrecognized parfile line '{p_line}'") + + if tm.UNITS.value == "TCB" and convert_tcb: + convert_tcb_tdb(tm) + return tm def _validate_components(self): @@ -166,9 +190,9 @@ def _get_component_param_overlap(self, component): # Add aliases compare overlap = in_param & cpm_param # translate to PINT parameter - overlap_pint_par = set( - [self.all_components.alias_to_pint_param(ovlp)[0] for ovlp in overlap] - ) + overlap_pint_par = { + self.all_components.alias_to_pint_param(ovlp)[0] for ovlp in overlap + } # The degree of overlapping for input component and compared component overlap_deg_in = len(component.params) - len(overlap_pint_par) overlap_deg_cpm = len(cp.params) - len(overlap_pint_par) @@ -246,38 +270,40 @@ def _pintify_parfile(self, parfile, allow_name_mixing=False): try: pint_name, init0 = self.all_components.alias_to_pint_param(k) except UnknownParameter: - if k in ignore_params: # Parameter is known but in the ingore list + if k in ignore_params: + # Parameter is known but in the ignore list continue - else: # Check ignored prefix - try: - pfx, idxs, idx = split_prefixed_name(k) - if pfx in ignore_prefix: # It is an ignored prefix. - continue - else: - unknown_param[k] += v - except PrefixError: + # Check ignored prefix + try: + pfx, idxs, idx = split_prefixed_name(k) + if pfx in ignore_prefix: # It is an ignored prefix. + continue + else: unknown_param[k] += v + except PrefixError: + unknown_param[k] += v continue pint_param_dict[pint_name] += v original_name_map[pint_name].append(k) repeating[pint_name] += len(v) # Check if this parameter is allowed to be repeated by PINT - if len(pint_param_dict[pint_name]) > 1: - if pint_name not in self.all_components.repeatable_param: - raise TimingModelError( - f"Parameter {pint_name} is not a repeatable parameter. " - f"However, multiple line use it." - ) + if ( + len(pint_param_dict[pint_name]) > 1 + and pint_name not in self.all_components.repeatable_param + ): + raise TimingModelError( + f"Parameter {pint_name} is not a repeatable parameter. " + f"However, multiple line use it." + ) # Check if the name is mixed for p_n, o_n in original_name_map.items(): - if len(o_n) > 1: - if not allow_name_mixing: - raise TimingModelError( - f"Parameter {p_n} have mixed input names/alias " - f"{o_n}. If you want to have mixing names, please use" - f" 'allow_name_mixing=True', and the output .par file " - f"will use '{original_name_map[pint_name][0]}'." - ) + if len(o_n) > 1 and not allow_name_mixing: + raise TimingModelError( + f"Parameter {p_n} have mixed input names/alias " + f"{o_n}. If you want to have mixing names, please use" + f" 'allow_name_mixing=True', and the output .par file " + f"will use '{original_name_map[pint_name][0]}'." + ) original_name_map[p_n] = o_n[0] return pint_param_dict, original_name_map, unknown_param @@ -347,12 +373,11 @@ def choose_model(self, param_inpar): if p_name != first_init: param_not_in_pint.append(pp) - p_cp = self.all_components.param_component_map.get(first_init, None) - if p_cp: + if p_cp := self.all_components.param_component_map.get(first_init, None): param_components_inpar[p_name] = p_cp # Back map the possible_components and the parameters in the parfile # This will remove the duplicate components. - conflict_components = defaultdict(set) # graph for confilict + conflict_components = defaultdict(set) # graph for conflict for k, cps in param_components_inpar.items(): # If `timing_model` in param --> component mapping skip # Timing model is the base. @@ -386,10 +411,10 @@ def choose_model(self, param_inpar): temp_cf_cp.remove(cp) conflict_components[cp].update(set(temp_cf_cp)) continue - # Check if the selected component in the confilict graph. If it is - # remove the selected componens with its conflict components. + # Check if the selected component in the conflict graph. If it is + # remove the selected components with its conflict components. for ps_cp in selected_components: - cf_cps = conflict_components.get(ps_cp, None) + cf_cps = conflict_components.get(ps_cp) if cf_cps is not None: # Had conflict, but resolved. for cf_cp in cf_cps: del conflict_components[cf_cp] @@ -398,17 +423,17 @@ def choose_model(self, param_inpar): selected_cates = {} for cp in selected_components: cate = self.all_components.component_category_map[cp] - if cate not in selected_cates.keys(): - selected_cates[cate] = cp - else: - exisit_cp = selected_cates[cate] + if cate in selected_cates: + exist_cp = selected_cates[cate] raise TimingModelError( - f"Component '{cp}' and '{exisit_cp}' belong to the" + f"Component '{cp}' and '{exist_cp}' belong to the" f" same category '{cate}'. Only one component from" f" the same category can be used for a timing model." f" Please check your input (e.g., .par file)." ) + else: + selected_cates[cate] = cp return selected_components, conflict_components, param_not_in_pint def _setup_model( @@ -418,6 +443,7 @@ def _setup_model( original_name=None, setup=True, validate=True, + allow_tcb=False, ): """Fill up a timing model with parameter values and then setup the model. @@ -443,29 +469,28 @@ def _setup_model( Whether to run the setup function in the timing model. validate : bool, optional Whether to run the validate function in the timing model. + allow_tcb : bool, optional + Whether to allow reading TCB par files """ - if original_name is not None: - use_alias = True - else: - use_alias = False + use_alias = original_name is not None for pp, v in pint_param_dict.items(): try: par = getattr(timing_model, pp) except AttributeError: - # since the input is pintfied, it should be an uninitized indexed parameter + # since the input is pintfied, it should be an uninitialized indexed parameter # double check if the missing parameter an indexed parameter. pint_par, first_init = self.all_components.alias_to_pint_param(pp) try: prefix, _, index = split_prefixed_name(pint_par) - except PrefixError: + except PrefixError as e: par_hosts = self.all_components.param_component_map[pint_par] - currnt_cp = timing_model.components.keys() + current_cp = timing_model.components.keys() raise TimingModelError( f"Parameter {pint_par} is recognized" f" by PINT, but not used in the current" f" timing model. It is used in {par_hosts}," - f" but the current timing model uses {currnt_cp}." - ) + f" but the current timing model uses {current_cp}." + ) from e # TODO need to create a better API for _locate_param_host host_component = timing_model._locate_param_host(first_init) timing_model.add_param_from_top( @@ -477,10 +502,7 @@ def _setup_model( # Fill up the values param_line = len(v) if param_line < 2: - if use_alias: - name = original_name[pp] - else: - name = pp + name = original_name[pp] if use_alias else pp par.from_parfile_line(" ".join([name] + v)) else: # For the repeatable parameters lines = copy.deepcopy(v) # Line queue. @@ -516,20 +538,20 @@ def _setup_model( # There is no current repeatable parameter matching the new line # First try to fill up an empty space. - if empty_repeat_param != []: - emt_par = empty_repeat_param.pop(0) - emt_par.from_parfile_line(" ".join([emt_par.name, li])) - if use_alias: # Use the input alias as input - emt_par.use_alias = original_name[pp] - else: + if not empty_repeat_param: # No empty space, add a new parameter to the timing model. host_component = timing_model._locate_param_host(pp) timing_model.add_param_from_top(temp_par, host_component[0][0]) + else: + emt_par = empty_repeat_param.pop(0) + emt_par.from_parfile_line(" ".join([emt_par.name, li])) + if use_alias: # Use the input alias as input + emt_par.use_alias = original_name[pp] if setup: timing_model.setup() if validate: - timing_model.validate() + timing_model.validate(allow_tcb=allow_tcb) return timing_model def _report_conflict(self, conflict_graph): @@ -538,12 +560,10 @@ def _report_conflict(self, conflict_graph): # Put all the conflict components together from the graph cf_cps = list(v) cf_cps.append(k) - raise ComponentConflict( - "Can not decide the one component from:" " {}".format(cf_cps) - ) + raise ComponentConflict(f"Can not decide the one component from: {cf_cps}") -def get_model(parfile, allow_name_mixing=False): +def get_model(parfile, allow_name_mixing=False, allow_tcb=False): """A one step function to build model from a parfile. Parameters @@ -557,25 +577,31 @@ def get_model(parfile, allow_name_mixing=False): T2EFAC and EFAC, both of them maps to PINT parameter EFAC, present in the parfile at the same time. + allow_tcb : True, False, or "raw", optional + Whether to read TCB par files. Default is False, and will throw an + error upon encountering TCB par files. If True, the par file will be + converted to TDB upon read. If "raw", an unconverted malformed TCB + TimingModel object will be returned. + Returns ------- Model instance get from parfile. """ + model_builder = ModelBuilder() try: contents = parfile.read() except AttributeError: contents = None - if contents is None: - # # parfile is a filename and can be handled by ModelBuilder - # if _model_builder is None: - # _model_builder = ModelBuilder() - model = model_builder(parfile, allow_name_mixing) - model.name = parfile - return model - else: - tm = model_builder(StringIO(contents), allow_name_mixing) - return tm + if contents is not None: + return model_builder(StringIO(contents), allow_name_mixing, allow_tcb=allow_tcb) + + # # parfile is a filename and can be handled by ModelBuilder + # if _model_builder is None: + # _model_builder = ModelBuilder() + model = model_builder(parfile, allow_name_mixing, allow_tcb=allow_tcb) + model.name = parfile + return model def get_model_and_toas( @@ -592,6 +618,7 @@ def get_model_and_toas( picklefilename=None, allow_name_mixing=False, limits="warn", + allow_tcb=False, ): """Load a timing model and a related TOAs, using model commands as needed @@ -630,12 +657,17 @@ def get_model_and_toas( in the parfile at the same time. limits : "warn" or "error" What to do when encountering TOAs for which clock corrections are not available. + allow_tcb : True, False, or "raw", optional + Whether to read TCB par files. Default is False, and will throw an + error upon encountering TCB par files. If True, the par file will be + converted to TDB upon read. If "raw", an unconverted malformed TCB + TimingModel object will be returned. Returns ------- A tuple with (model instance, TOAs instance) """ - mm = get_model(parfile, allow_name_mixing) + mm = get_model(parfile, allow_name_mixing, allow_tcb=allow_tcb) tt = get_TOAs( timfile, include_pn=include_pn, diff --git a/src/pint/models/solar_system_shapiro.py b/src/pint/models/solar_system_shapiro.py index c06d5a366..1cd850065 100644 --- a/src/pint/models/solar_system_shapiro.py +++ b/src/pint/models/solar_system_shapiro.py @@ -111,13 +111,13 @@ def solar_system_shapiro_delay(self, toas, acc_delay=None): if self.PLANET_SHAPIRO.value: for pl in ("jupiter", "saturn", "venus", "uranus", "neptune"): delay[grp] += self.ss_obj_shapiro_delay( - tbl[grp]["obs_" + pl + "_pos"], + tbl[grp][f"obs_{pl}_pos"], psr_dir, self._ss_mass_sec[pl], ) - except KeyError: + except KeyError as e: raise KeyError( "Planet positions not found when trying to compute Solar System Shapiro delay. " "Make sure that you include `planets=True` in your `get_TOAs()` call, or use `get_model_and_toas()`." - ) + ) from e return delay * u.second diff --git a/src/pint/models/tcb_conversion.py b/src/pint/models/tcb_conversion.py new file mode 100644 index 000000000..c4adccfed --- /dev/null +++ b/src/pint/models/tcb_conversion.py @@ -0,0 +1,172 @@ +"""TCB to TDB conversion of a timing model.""" + +import numpy as np + +from pint.models.parameter import MJDParameter +from loguru import logger as log + +__all__ = [ + "IFTE_K", + "scale_parameter", + "transform_mjd_parameter", + "convert_tcb_to_tdb", +] + +# These constants are taken from Irwin & Fukushima 1999. +# These are the same as the constants used in tempo2 as of 10 Feb 2023. +IFTE_MJD0 = np.longdouble("43144.0003725") +IFTE_KM1 = np.longdouble("1.55051979176e-8") +IFTE_K = 1 + IFTE_KM1 + + +def scale_parameter(model, param, n, backwards): + """Scale a parameter x by a power of IFTE_K + x_tdb = x_tcb * IFTE_K**n + + The power n depends on the "effective dimensionality" of + the parameter as it appears in the timing model. Some examples + are given bellow: + + 1. F0 has effective dimensionality of frequency and n = 1 + 2. F1 has effective dimensionality of frequency^2 and n = 2 + 3. A1 has effective dimensionality of time because it appears as + A1/c in the timing model. Therefore, its n = -1 + 4. DM has effective dimensionality of frequency because it appears + as DM*DMconst in the timing model. Therefore, its n = 1 + 5. PBDOT is dimensionless and has n = 0. i.e., it is not scaled. + + Parameter + --------- + model : pint.models.timing_model.TimingModel + The timing model + param : str + The parameter name to be converted + n : int + The power of IFTE_K in the scaling factor + backwards : bool + Whether to do TDB to TCB conversion. + """ + assert isinstance(n, int), "The power must be an integer." + + p = -1 if backwards else 1 + + factor = IFTE_K ** (p * n) + + if hasattr(model, param) and getattr(model, param).quantity is not None: + par = getattr(model, param) + par.value *= factor + if par.uncertainty_value is not None: + par.uncertainty_value *= factor + + +def transform_mjd_parameter(model, param, backwards): + """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 + + Parameters + ---------- + model : pint.models.timing_model.TimingModel + The timing model + param : str + The parameter name to be converted + backwards : bool + Whether to do TDB to TCB conversion. + """ + 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) + + 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): + """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 + + 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 + + Parameters + ---------- + model : pint.models.timing_model.TimingModel + Timing model to be converted. + backwards : bool + Whether to do TDB to TCB conversion. The default is TCB to TDB. + """ + + target_units = "TCB" if backwards else "TDB" + + 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 + + log.warning( + "Converting this timing model from TCB to TDB. " + "Please note that the TCB to TDB conversion is only approximate and " + "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 + + model.validate(allow_tcb=backwards) diff --git a/src/pint/models/timing_model.py b/src/pint/models/timing_model.py index 16d119df9..d119e4119 100644 --- a/src/pint/models/timing_model.py +++ b/src/pint/models/timing_model.py @@ -26,6 +26,8 @@ See :ref:`Timing Models` for more details on how PINT's timing models work. """ + + import abc import copy import inspect @@ -78,7 +80,7 @@ "MissingBinaryError", "UnknownBinaryModel", ] -# Parameters or lines in parfiles we don't understand but shouldn't +# Parameters or lines in par files we don't understand but shouldn't # complain about. These are still passed to components so that they # can use them if they want to. # @@ -86,22 +88,19 @@ # errors in the par file. # # Comparisons with keywords in par file lines is done in a case insensitive way. -ignore_params = set( - [ - "TRES", - "TZRMJD", - "TZRFRQ", - "TZRSITE", - "NITS", - "IBOOT", - "CHI2R", - "MODE", - "PLANET_SHAPIRO2", - # 'NE_SW', 'NE_SW2', - ] -) - -ignore_prefix = set(["DMXF1_", "DMXF2_", "DMXEP_"]) # DMXEP_ for now. +ignore_params = { + "TRES", + "TZRMJD", + "TZRFRQ", + "TZRSITE", + "NITS", + "IBOOT", + "CHI2R", + "MODE", + "PLANET_SHAPIRO2", +} + +ignore_prefix = {"DMXF1_", "DMXF2_", "DMXEP_"} DEFAULT_ORDER = [ "astrometry", @@ -359,7 +358,7 @@ def __repr__(self): def __str__(self): return self.as_parfile() - def validate(self): + def validate(self, allow_tcb=False): """Validate component setup. The checks include required parameters and parameter values. @@ -373,22 +372,34 @@ def validate(self): if self.T2CMETHOD.value not in [None, "IAU2000B"]: # FIXME: really? warn("PINT only supports 'T2CMETHOD IAU2000B'") self.T2CMETHOD.value = "IAU2000B" - if self.UNITS.value not in [None, "TDB"]: - if self.UNITS.value == "TCB": - error_message = """The TCB timescale is not supported by PINT. (PINT only supports 'UNITS TDB'.) - See https://nanograv-pint.readthedocs.io/en/latest/explanation.html#time-scales for an explanation - on different timescales. The par file can be converted from TCB to TDB using the `transform` - plugin of TEMPO2 like so: - $ tempo2 -gr transform J1234+6789_tcb.par J1234+6789_tdb.par tdb + + if self.UNITS.value not in [None, "TDB", "TCB"]: + error_message = f"PINT only supports 'UNITS TDB'. The given timescale '{self.UNITS.value}' is invalid." + raise ValueError(error_message) + elif self.UNITS.value == "TCB": + if not allow_tcb: + error_message = """The TCB timescale is not fully supported by PINT. + PINT only supports 'UNITS TDB' internally. See https://nanograv-pint.readthedocs.io/en/latest/explanation.html#time-scales + for an explanation on different timescales. A TCB par file can be + converted to TDB using the `tcb2tdb` command like so: + + $ tcb2tdb J1234+6789_tcb.par J1234+6789_tdb.par + + However, this conversion is not exact and a fit must be performed to obtain + reliable results. Note that PINT only supports writing TDB par files. """ + raise ValueError(error_message) else: - error_message = f"PINT only supports 'UNITS TDB'. The given timescale '{self.UNITS.value}' is invalid." - raise ValueError(error_message) + log.warning( + "PINT does not support 'UNITS TCB' internally. Reading this par file nevertheless " + "because the `allow_tcb` option was given. This `TimingModel` object should not be " + "used for anything except converting to TDB." + ) if not self.START.frozen: - warn("START cannot be unfrozen...setting START.frozen to True") + warn("START cannot be unfrozen... Setting START.frozen to True") self.START.frozen = True if not self.FINISH.frozen: - warn("FINISH cannot be unfrozen...setting FINISH.frozen to True") + warn("FINISH cannot be unfrozen... Setting FINISH.frozen to True") self.FINISH.frozen = True for cp in self.components.values(): @@ -438,34 +449,30 @@ def params_ordered(self): used_cats = set() pstart = copy.copy(self.top_level_params) for cat in start_order: - if cat in compdict: - cp = compdict[cat] - for cpp in cp: - pstart += cpp.params - used_cats.add(cat) - else: + if cat not in compdict: continue - + cp = compdict[cat] + for cpp in cp: + pstart += cpp.params + used_cats.add(cat) pend = [] for cat in last_order: - if cat in compdict: - cp = compdict[cat] - for cpp in cp: - pend += cpp.parms - used_cats.add(cat) - else: + if cat not in compdict: continue + cp = compdict[cat] + for cpp in cp: + pend += cpp.parms + used_cats.add(cat) # Now collect any components that haven't already been included in the list pmid = [] for cat in compdict: if cat in used_cats: continue - else: - cp = compdict[cat] - for cpp in cp: - pmid += cpp.params - used_cats.add(cat) + cp = compdict[cat] + for cpp in cp: + pmid += cpp.params + used_cats.add(cat) return pstart + pmid + pend @@ -569,7 +576,7 @@ def get_params_of_component_type(self, component_type): ------- list """ - component_type_list_str = "{}_list".format(component_type) + component_type_list_str = f"{component_type}_list" if hasattr(self, component_type_list_str): component_type_list = getattr(self, component_type_list_str) return [ @@ -600,17 +607,14 @@ def set_param_uncertainties(self, fitp): """Set the model parameters to the value contained in the input dict.""" for k, v in fitp.items(): p = getattr(self, k) - if isinstance(v, u.Quantity): - p.uncertainty = v - else: - p.uncertainty = v * p.units + p.uncertainty = v if isinstance(v, u.Quantity) else v * p.units @property_exists def components(self): """All the components in a dictionary indexed by name.""" comps = {} for ct in self.component_types: - for cp in getattr(self, ct + "_list"): + for cp in getattr(self, f"{ct}_list"): comps[cp.__class__.__name__] = cp return comps @@ -698,13 +702,11 @@ def orbital_phase(self, barytimes, anom="mean", radians=True): elif anom.lower() == "true": anoms = bbi.nu() # can be negative else: - raise ValueError("anom='%s' is not a recognized type of anomaly" % anom) + raise ValueError(f"anom='{anom}' is not a recognized type of anomaly") # Make sure all angles are between 0-2*pi anoms = np.remainder(anoms.value, 2 * np.pi) - if radians: # return with radian units - return anoms * u.rad - else: # return as unitless cycles from 0-1 - return anoms / (2 * np.pi) + # return with radian units or return as unitless cycles from 0-1 + return anoms * u.rad if radians else anoms / (2 * np.pi) def conjunction(self, baryMJD): """Return the time(s) of the first superior conjunction(s) after baryMJD. @@ -860,13 +862,13 @@ def d_phase_d_delay_funcs(self): def get_deriv_funcs(self, component_type, derivative_type=""): """Return dictionary of derivative functions.""" - # TODO, this function can be a more generical function collector. + # TODO, this function can be a more generic function collector. deriv_funcs = defaultdict(list) - if not derivative_type == "": + if derivative_type != "": derivative_type += "_" - for cp in getattr(self, component_type + "_list"): + for cp in getattr(self, f"{component_type}_list"): try: - df = getattr(cp, derivative_type + "deriv_funcs") + df = getattr(cp, f"{derivative_type}deriv_funcs") except AttributeError: continue for k, v in df.items(): @@ -909,13 +911,11 @@ def get_component_type(self, component): comp_base = inspect.getmro(component.__class__) if comp_base[-2].__name__ != "Component": raise TypeError( - "Class '%s' is not a Component type class." - % component.__class__.__name__ + f"Class '{component.__class__.__name__}' is not a Component type class." ) elif len(comp_base) < 3: raise TypeError( - "'%s' class is not a subclass of 'Component' class." - % component.__class__.__name__ + f"'{component.__class__.__name__}' class is not a subclass of 'Component' class." ) else: comp_type = comp_base[-3].__name__ @@ -943,17 +943,16 @@ def map_component(self, component): comps = self.components if isinstance(component, str): if component not in list(comps.keys()): - raise AttributeError("No '%s' in the timing model." % component) + raise AttributeError(f"No '{component}' in the timing model.") comp = comps[component] - else: # When component is an component instance. - if component not in list(comps.values()): - raise AttributeError( - "No '%s' in the timing model." % component.__class__.__name__ - ) - else: - comp = component + elif component in list(comps.values()): + comp = component + else: + raise AttributeError( + f"No '{component.__class__.__name__}' in the timing model." + ) comp_type = self.get_component_type(comp) - host_list = getattr(self, comp_type + "_list") + host_list = getattr(self, f"{comp_type}_list") order = host_list.index(comp) return comp, order, host_list, comp_type @@ -972,9 +971,9 @@ def add_component( If true, add a duplicate component. Default is False. """ comp_type = self.get_component_type(component) + cur_cps = [] if comp_type in self.component_types: - comp_list = getattr(self, comp_type + "_list") - cur_cps = [] + comp_list = getattr(self, f"{comp_type}_list") for cp in comp_list: # If component order is not defined. cp_order = ( @@ -984,32 +983,28 @@ def add_component( # Check if the component has been added already. if component.__class__ in (x.__class__ for x in comp_list): log.warning( - "Component '%s' is already present but was added again." - % component.__class__.__name__ + f"Component '{component.__class__.__name__}' is already present but was added again." ) if not force: raise ValueError( - "Component '%s' is already present and will not be " - "added again. To force add it, use force=True option." - % component.__class__.__name__ + f"Component '{component.__class__.__name__}' is already present and will not be " + f"added again. To force add it, use force=True option." ) else: self.component_types.append(comp_type) - cur_cps = [] - # link new component to TimingModel component._parent = self - # If the categore is not in the order list, it will be added to the end. + # If the category is not in the order list, it will be added to the end. if component.category not in order: - new_cp = tuple((len(order) + 1, component)) + new_cp = len(order) + 1, component else: - new_cp = tuple((order.index(component.category), component)) + new_cp = order.index(component.category), component # add new component cur_cps.append(new_cp) cur_cps.sort(key=lambda x: x[0]) new_comp_list = [c[1] for c in cur_cps] - setattr(self, comp_type + "_list", new_comp_list) + setattr(self, f"{comp_type}_list", new_comp_list) # Set up components if setup: self.setup() @@ -1041,8 +1036,8 @@ def _locate_param_host(self, param): list of tuples All possible components that host the target parameter. The first element is the component object that have the target parameter, the - second one is the parameter object. If it is a prefix-style parameter - , it will return one example of such parameter. + second one is the parameter object. If it is a prefix-style parameter, + it will return one example of such parameter. """ result_comp = [] for cp_name, cp in self.components.items(): @@ -1085,12 +1080,12 @@ def add_param_from_top(self, param, target_component, setup=False): if target_component == "": setattr(self, param.name, param) self.top_level_params += [param.name] - else: - if target_component not in list(self.components.keys()): - raise AttributeError( - "Can not find component '%s' in " "timing model." % target_component - ) + elif target_component in list(self.components.keys()): self.components[target_component].add_param(param, setup=setup) + else: + raise AttributeError( + f"Can not find component '{target_component}' in " "timing model." + ) def remove_param(self, param): """Remove a parameter from timing model. @@ -1102,7 +1097,7 @@ def remove_param(self, param): """ param_map = self.get_params_mapping() if param not in param_map: - raise AttributeError("Can not find '%s' in timing model." % param) + raise AttributeError(f"Can not find '{param}' in timing model.") if param_map[param] == "timing_model": delattr(self, param) self.top_level_params.remove(param) @@ -1112,10 +1107,8 @@ def remove_param(self, param): self.setup() def get_params_mapping(self): - """Report whick component each parameter name comes from.""" - param_mapping = {} - for p in self.top_level_params: - param_mapping[p] = "timing_model" + """Report which component each parameter name comes from.""" + param_mapping = {p: "timing_model" for p in self.top_level_params} for cp in list(self.components.values()): for pp in cp.params: param_mapping[pp] = cp.__class__.__name__ @@ -1138,7 +1131,7 @@ def get_prefix_mapping(self, prefix): Returns ------- dict - A dictionary with prefix pararameter real index as key and parameter + A dictionary with prefix parameter real index as key and parameter name as value. """ for cp in self.components.values(): @@ -1224,13 +1217,12 @@ def delay(self, toas, cutoff_component="", include_last=True): idx = len(self.DelayComponent_list) else: delay_names = [x.__class__.__name__ for x in self.DelayComponent_list] - if cutoff_component in delay_names: - idx = delay_names.index(cutoff_component) - if include_last: - idx += 1 - else: - raise KeyError("No delay component named '%s'." % cutoff_component) + if cutoff_component not in delay_names: + raise KeyError(f"No delay component named '{cutoff_component}'.") + idx = delay_names.index(cutoff_component) + if include_last: + idx += 1 # Do NOT cycle through delay_funcs - cycle through components until cutoff for dc in self.DelayComponent_list[:idx]: for df in dc.delay_funcs_component: @@ -1253,30 +1245,24 @@ def phase(self, toas, abs_phase=None): # abs_phase defaults to True if AbsPhase is in the model, otherwise to # False. Of course, if you manually set it, it will use that setting. if abs_phase is None: - if "AbsPhase" in list(self.components.keys()): - abs_phase = True - else: - abs_phase = False - # If the absolute phase flag is on, use the TZR parameters to compute - # the absolute phase. - if abs_phase: - if "AbsPhase" not in list(self.components.keys()): - # if no absolute phase (TZRMJD), add the component to the model and calculate it - from pint.models import absolute_phase - - self.add_component(absolute_phase.AbsPhase(), validate=False) - self.make_TZR_toa( - toas - ) # TODO:needs timfile to get all toas, but model doesn't have access to timfile. different place for this? - self.validate() - tz_toa = self.get_TZR_toa(toas) - tz_delay = self.delay(tz_toa) - tz_phase = Phase(np.zeros(len(toas.table)), np.zeros(len(toas.table))) - for pf in self.phase_funcs: - tz_phase += Phase(pf(tz_toa, tz_delay)) - return phase - tz_phase - else: + abs_phase = "AbsPhase" in list(self.components.keys()) + if not abs_phase: return phase + if "AbsPhase" not in list(self.components.keys()): + # if no absolute phase (TZRMJD), add the component to the model and calculate it + from pint.models import absolute_phase + + self.add_component(absolute_phase.AbsPhase(), validate=False) + self.make_TZR_toa( + toas + ) # TODO:needs timfile to get all toas, but model doesn't have access to timfile. different place for this? + self.validate() + tz_toa = self.get_TZR_toa(toas) + tz_delay = self.delay(tz_toa) + tz_phase = Phase(np.zeros(len(toas.table)), np.zeros(len(toas.table))) + for pf in self.phase_funcs: + tz_phase += Phase(pf(tz_toa, tz_delay)) + return phase - tz_phase def total_dm(self, toas): """Calculate dispersion measure from all the dispersion type of components.""" @@ -1391,8 +1377,8 @@ def noise_model_dimensions(self, toas): """Number of basis functions for each noise model component. Returns a dictionary of correlated-noise components in the noise - model. Each entry contains a tuple (offset, size) where size is the - number of basis funtions for the component, and offset is their + model. Each entry contains a tuple (offset, size) where size is the + number of basis functions for the component, and offset is their starting location in the design matrix and weights vector. """ result = {} @@ -1456,13 +1442,9 @@ def jump_flags_to_params(self, toas): if tjv in tim_jump_values: log.info(f"JUMP -tim_jump {tjv} already exists") tim_jump_values.remove(tjv) - if used_indices: - num = max(used_indices) + 1 - else: - num = 1 - + num = max(used_indices) + 1 if used_indices else 1 if not tim_jump_values: - log.info(f"All tim_jump values have corresponding JUMPs") + log.info("All tim_jump values have corresponding JUMPs") return # FIXME: arrange for these to be in a sensible order (might not be integers @@ -1508,7 +1490,7 @@ def delete_jump_and_flags(self, toa_table, jump_num): Specifies the index of the jump to be deleted. """ # remove jump of specified index - self.remove_param("JUMP" + str(jump_num)) + self.remove_param(f"JUMP{jump_num}") # remove jump flags from selected TOA tables if toa_table is not None: @@ -1573,7 +1555,7 @@ def d_phase_d_toa(self, toas, sample_step=None): if sample_step is None: pulse_period = 1.0 / (self.F0.quantity) sample_step = pulse_period * 2 - # Note that sample_dt is applied cumulatively, so this evaulates phase at TOA-dt and TOA+dt + # Note that sample_dt is applied cumulatively, so this evaluates phase at TOA-dt and TOA+dt sample_dt = [-sample_step, 2 * sample_step] sample_phase = [] @@ -1661,8 +1643,8 @@ def d_delay_d_param(self, toas, param, acc_delay=None): delay_derivs = self.delay_deriv_funcs if param not in list(delay_derivs.keys()): raise AttributeError( - "Derivative function for '%s' is not provided" - " or not registered. " % param + "Derivative function for '{param}' is not provided" + " or not registered. " ) for df in delay_derivs[param]: result += df(toas, param, acc_delay).to( @@ -1679,10 +1661,7 @@ def d_phase_d_param_num(self, toas, param, step=1e-2): par = getattr(self, param) ori_value = par.value unit = par.units - if ori_value == 0: - h = 1.0 * step - else: - h = ori_value * step + h = 1.0 * step if ori_value == 0 else ori_value * step parv = [par.value - h, par.value + h] phase_i = ( @@ -1713,13 +1692,10 @@ def d_delay_d_param_num(self, toas, param, step=1e-2): ori_value = par.value if ori_value is None: # A parameter did not get to use in the model - log.warning("Parameter '%s' is not used by timing model." % param) + log.warning(f"Parameter '{param}' is not used by timing model.") return np.zeros(toas.ntoas) * (u.second / par.units) unit = par.units - if ori_value == 0: - h = 1.0 * step - else: - h = ori_value * step + h = 1.0 * step if ori_value == 0 else ori_value * step parv = [par.value - h, par.value + h] delay = np.zeros((toas.ntoas, 2)) for ii, val in enumerate(parv): @@ -1739,7 +1715,7 @@ def d_dm_d_param(self, data, param): result = np.zeros(len(data)) << (u.pc / u.cm**3 / par.units) dm_df = self.dm_derivs.get(param, None) if dm_df is None: - if param not in self.params: # Maybe add differentitable params + if param not in self.params: # Maybe add differentiable params raise AttributeError(f"Parameter {param} does not exist") else: return result diff --git a/src/pint/scripts/tcb2tdb.py b/src/pint/scripts/tcb2tdb.py new file mode 100644 index 000000000..a920c28bc --- /dev/null +++ b/src/pint/scripts/tcb2tdb.py @@ -0,0 +1,53 @@ +"""PINT-based tool for converting TCB par files to TDB.""" + +import argparse + +from loguru import logger as log + +import pint.logging +from pint.models.model_builder import ModelBuilder +from pint.models.tcb_conversion import convert_tcb_tdb + +pint.logging.setup(level="INFO") + +__all__ = ["main"] + + +def main(argv=None): + parser = argparse.ArgumentParser( + 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 + """, + formatter_class=argparse.ArgumentDefaultsHelpFormatter, + ) + parser.add_argument("input_par", help="Input par file name (TCB)") + parser.add_argument("output_par", help="Output par file name (TDB)") + + args = parser.parse_args(argv) + + mb = ModelBuilder() + model = mb(args.input_par, allow_tcb=True) + model.write_parfile(args.output_par) + + log.info(f"Output written to {args.output_par}.") diff --git a/tests/test_tcb2tdb.py b/tests/test_tcb2tdb.py new file mode 100644 index 000000000..34e7373a1 --- /dev/null +++ b/tests/test_tcb2tdb.py @@ -0,0 +1,83 @@ +"""Tests for `pint.models.tcb_conversion` and the `tcb2tdb` script.""" + +import os +from copy import deepcopy +from io import StringIO + +import numpy as np +import pytest + +from pint.models.model_builder import ModelBuilder +from pint.models.tcb_conversion import convert_tcb_tdb +from pint.scripts import tcb2tdb + +simplepar = """ +PSR PSRTEST +RAJ 17:48:52.75 1 +DECJ -20:21:29.0 1 +F0 61.485476554 1 +F1 -1.181D-15 1 +PEPOCH 53750.000000 +POSEPOCH 53750.000000 +DM 223.9 1 +SOLARN0 0.00 +BINARY BT +T0 53750 +A1 100.0 1 0.1 +ECC 1.0 +OM 0.0 +PB 10.0 +EPHEM DE436 +CLK TT(BIPM2017) +UNITS TCB +TIMEEPH FB90 +T2CMETHOD TEMPO +CORRECT_TROPOSPHERE N +PLANET_SHAPIRO Y +DILATEFREQ N +""" + + +@pytest.mark.parametrize("backwards", [True, False]) +def test_convert_units(backwards): + with pytest.raises(ValueError): + m = ModelBuilder()(StringIO(simplepar)) + + m = ModelBuilder()(StringIO(simplepar), allow_tcb="raw") + f0_tcb = m.F0.value + pb_tcb = m.PB.value + convert_tcb_tdb(m, backwards=backwards) + assert m.UNITS.value == ("TCB" if backwards else "TDB") + assert np.isclose(m.F0.value / f0_tcb, pb_tcb / m.PB.value) + + +def test_convert_units_roundtrip(): + m = ModelBuilder()(StringIO(simplepar), allow_tcb="raw") + m_ = deepcopy(m) + convert_tcb_tdb(m, backwards=False) + convert_tcb_tdb(m, backwards=True) + + for par in m.params: + p = getattr(m, par) + p_ = getattr(m_, par) + if p.value is None: + assert p_.value is None + elif isinstance(p.value, str): + assert getattr(m, par).value == getattr(m_, par).value + else: + assert np.isclose(getattr(m, par).value, getattr(m_, par).value) + + +def test_tcb2tdb(tmp_path): + tmppar1 = tmp_path / "tmp1.par" + tmppar2 = tmp_path / "tmp2.par" + with open(tmppar1, "w") as f: + f.write(simplepar) + + cmd = f"{tmppar1} {tmppar2}" + tcb2tdb.main(cmd.split()) + + assert os.path.isfile(tmppar2) + + m2 = ModelBuilder()(tmppar2) + assert m2.UNITS.value == "TDB"