diff --git a/CHANGELOG-unreleased.md b/CHANGELOG-unreleased.md index d68f68b91..cbb483399 100644 --- a/CHANGELOG-unreleased.md +++ b/CHANGELOG-unreleased.md @@ -37,6 +37,7 @@ the released changes. - `freeze_params` option in `wavex_setup` and `dmwavex_setup` - `plrednoise_from_wavex`, `pldmnoise_from_dmwavex`, and `find_optimal_nharms` functions - fake TOAs can be created with `subtract_mean=False`, to maintain phase coherence between different data sets +- Binary models can be guessed by the `ModelBuilder`. Options and script are added to allow reading/conversion of the T2 binary model - Better explanation of ELL1H behavior when H3/H4/STIGMA supplied and when NHARMS is used - FDJumpDM component for System-dependent DM offsets - Documentation: Explanation for FDJUMP and FDJUMPDM diff --git a/setup.cfg b/setup.cfg index ad2570a5f..28dfa2974 100644 --- a/setup.cfg +++ b/setup.cfg @@ -64,6 +64,7 @@ console_scripts = convert_parfile = pint.scripts.convert_parfile:main compare_parfiles = pint.scripts.compare_parfiles:main tcb2tdb = pint.scripts.tcb2tdb:main + t2binary2pint = pint.scripts.t2binary2pint:main pintpublish = pint.scripts.pintpublish:main diff --git a/src/pint/models/model_builder.py b/src/pint/models/model_builder.py index f81d78fd8..677da6b86 100644 --- a/src/pint/models/model_builder.py +++ b/src/pint/models/model_builder.py @@ -34,10 +34,24 @@ get_unit, ) from pint.models.tcb_conversion import convert_tcb_tdb +from pint.models.binary_ddk import _convert_kin, _convert_kom __all__ = ["ModelBuilder", "get_model", "get_model_and_toas"] default_models = ["StandardTimingModel"] +_binary_model_priority = [ + "Isolated", + "BT", + "BT_piecewise", + "ELL1", + "ELL1H", + "ELL1k", + "DD", + "DDK", + "DDGR", + "DDS", + "DDH", +] class ComponentConflict(ValueError): @@ -112,6 +126,8 @@ def __call__( parfile, allow_name_mixing=False, allow_tcb=False, + allow_T2=False, + force_binary_model=None, toas_for_tzr=None, **kwargs, ): @@ -134,6 +150,16 @@ def __call__( converted to TDB upon read. If "raw", an unconverted malformed TCB TimingModel object will be returned. + allow_T2 : bool, optional + Whether to convert a T2 binary model to an appropriate underlying + binary model. Default is False, and will throw an error upon + encountering the T2 binary model. If True, the binary model will be + converted to the most appropriate PINT-compatible binary model. + + force_binary_model : str, optional + When set to some binary model, like force_binary_model="DD", this + will override the binary model set in the parfile. Defaults to None + toas_for_tzr : TOAs or None, optional If this is not None, a TZR TOA (AbsPhase) will be created using the given TOAs object. @@ -151,6 +177,8 @@ def __call__( convert_tcb = allow_tcb == True allow_tcb_ = allow_tcb in [True, "raw"] + assert isinstance(allow_T2, bool) + pint_param_dict, original_name, unknown_param = self._pintify_parfile( parfile, allow_name_mixing ) @@ -168,7 +196,9 @@ def __call__( original_name[k] = k else: remaining_args[k] = v - selected, conflict, param_not_in_pint = self.choose_model(pint_param_dict) + selected, conflict, param_not_in_pint = self.choose_model( + pint_param_dict, force_binary_model=force_binary_model, allow_T2=allow_T2 + ) selected.update(set(self.default_components)) # Add SolarSystemShapiro only if an Astrometry component is present. @@ -404,7 +434,7 @@ def _pintify_parfile(self, parfile, allow_name_mixing=False): return pint_param_dict, original_name_map, unknown_param - def choose_model(self, param_inpar): + def choose_model(self, param_inpar, force_binary_model=None, allow_T2=False): """Choose the model components based on the parfile. Parameters @@ -413,6 +443,16 @@ def choose_model(self, param_inpar): Dictionary of the unique parameters in .par file with the key is the parfile line. :func:`parse_parfile` returns this dictionary. + allow_T2 : bool, optional + Whether to convert a T2 binary model to an appropriate underlying + binary model. Default is False, and will throw an error upon + encountering the T2 binary model. If True, the binary model will be + converted to the most appropriate PINT-compatible binary model. + + force_binary_model : str, optional + When set to some binary model, like force_binary_model="DD", this + will override the binary model set in the parfile. Defaults to None + Returns ------- list @@ -445,10 +485,13 @@ def choose_model(self, param_inpar): # build the base fo the timing model # pint_param_dict, unknown_param = self._pintify_parfile(param_inpar) binary = param_inpar.get("BINARY", None) - if binary is not None: + + if binary: binary = binary[0] - binary_cp = self.all_components.search_binary_components(binary) - selected_components.add(binary_cp.__class__.__name__) + selected_components.add( + self.choose_binary_model(param_inpar, force_binary_model, allow_T2) + ) + # 2. Get the component list from the parameters in the parfile. # 2.1 Check the aliases of input parameters. # This does not include the repeating parameters, but it should not @@ -532,6 +575,80 @@ def choose_model(self, param_inpar): selected_cates[cate] = cp return selected_components, conflict_components, param_not_in_pint + def choose_binary_model(self, param_inpar, force_binary_model=None, allow_T2=False): + """Choose the BINARY model based on the parfile. + + Parameters + ---------- + param_inpar: dict + Dictionary of the unique parameters in .par file with the key is the + parfile line. :func:`parse_parfile` returns this dictionary. + + force_binary_model : str, optional + When set to some binary model, like force_binary_model="DD", this + will override the binary model set in the parfile. Defaults to None + + allow_T2 : bool, optional + Whether to convert a T2 binary model to an appropriate underlying + binary model. Default is False, and will throw an error upon + encountering the T2 binary model. If True, the binary model will be + converted to the most appropriate PINT-compatible binary model. + + Returns + ------- + str + Name of the binary component + + Note + ---- + If the binary model does not have a PINT model (e.g. the T2 model), an + error is thrown with the suggested model that could replace it. If + allow_T2 is set to True, the most appropriate binary model is guessed + and used. If an appropriate model cannot be found, no suggestion is + given and an error is thrown. + """ + binary = param_inpar["BINARY"][0] + + # Guess what the binary model should be, regardless of BINARY parameter + try: + binary_model_guesses = guess_binary_model(param_inpar) + except UnknownBinaryModel as e: + log.error( + "Unable to find suitable binary model that has all the" + "parameters in the parfile. Please fix the par file." + ) + + # Allow for T2 model, gracefully + if force_binary_model is not None and binary != "T2": + binary = force_binary_model + elif binary == "T2" and allow_T2: + binary = binary_model_guesses[0] + log.warning( + f"Found T2 binary model. Gracefully converting T2 to: {binary}." + ) + + # Make sure that DDK parameters are properly converted + convert_binary_params_dict(param_inpar, force_binary_model=binary) + + try: + binary_cp = self.all_components.search_binary_components(binary) + + except UnknownBinaryModel as e: + log.error(f"Could not find binary model {binary}") + + log.info( + f"Compatible models with these parameters: {', '.join(binary_model_guesses)}." + ) + + # Re-raise the error, with an added guess for the binary model if we have one + if binary_model_guesses: + raise UnknownBinaryModel( + str(e), suggestion=binary_model_guesses[0] + ) from None + raise + + return binary_cp.__class__.__name__ + def _setup_model( self, timing_model, @@ -660,7 +777,13 @@ def _report_conflict(self, conflict_graph): def get_model( - parfile, allow_name_mixing=False, allow_tcb=False, toas_for_tzr=None, **kwargs + parfile, + allow_name_mixing=False, + allow_tcb=False, + allow_T2=False, + force_binary_model=None, + toas_for_tzr=None, + **kwargs, ): """A one step function to build model from a parfile. @@ -681,6 +804,16 @@ def get_model( converted to TDB upon read. If "raw", an unconverted malformed TCB TimingModel object will be returned. + allow_T2 : bool, optional + Whether to convert a T2 binary model to an appropriate underlying + binary model. Default is False, and will throw an error upon + encountering the T2 binary model. If True, the binary model will be + converted to the most appropriate PINT-compatible binary model. + + force_binary_model : str, optional + When set to some binary model, like force_binary_model="DD", this will + override the binary model set in the parfile. Defaults to None + toas_for_tzr : TOAs or None, optional If this is not None, a TZR TOA (AbsPhase) will be created using the given TOAs object. @@ -702,6 +835,8 @@ def get_model( StringIO(contents), allow_name_mixing, allow_tcb=allow_tcb, + allow_T2=allow_T2, + force_binary_model=force_binary_model, toas_for_tzr=toas_for_tzr, **kwargs, ) @@ -713,6 +848,8 @@ def get_model( parfile, allow_name_mixing, allow_tcb=allow_tcb, + allow_T2=allow_T2, + force_binary_model=force_binary_model, toas_for_tzr=toas_for_tzr, **kwargs, ) @@ -736,6 +873,8 @@ def get_model_and_toas( allow_name_mixing=False, limits="warn", allow_tcb=False, + allow_T2=False, + force_binary_model=None, add_tzr_to_model=True, **kwargs, ): @@ -784,6 +923,14 @@ def get_model_and_toas( 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. + allow_T2 : bool, optional + Whether to convert a T2 binary model to an appropriate underlying + binary model. Default is False, and will throw an error upon + encountering the T2 binary model. If True, the binary model will be + converted to the most appropriate PINT-compatible binary model. + force_binary_model : str, optional + When set to some binary model, like force_binary_model="DD", this + will override the binary model set in the parfile. Defaults to None add_tzr_to_model : bool, optional Create a TZR TOA in the timing model using the created TOAs object. Default is True. @@ -795,7 +942,14 @@ def get_model_and_toas( A tuple with (model instance, TOAs instance) """ - mm = get_model(parfile, allow_name_mixing, allow_tcb=allow_tcb, **kwargs) + mm = get_model( + parfile, + allow_name_mixing, + allow_tcb=allow_tcb, + allow_T2=allow_T2, + force_binary_model=force_binary_model, + **kwargs, + ) tt = get_TOAs( timfile, @@ -817,3 +971,135 @@ def get_model_and_toas( mm.add_tzr_toa(tt) return mm, tt + + +def guess_binary_model(parfile_dict): + """Based on the PINT parameter dictionary, guess the binary model + + Parameters + ---------- + parfile_dict + The parameter dictionary as read-in by parse_parfile + + Returns + ------- + list: + A priority-ordered list of possible binary models. The first one is the + best-guess + + """ + + 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) + + all_components = AllComponents() + binary_models = all_components.category_component_map["pulsar_system"] + + # Find all binary parameters + binary_parameters_map = { + all_components.components[binary_model].binary_model_name: add_sini( + all_components.search_binary_components(binary_model).aliases_map.keys() + ) + for binary_model in binary_models + } + binary_parameters_map.update({"Isolated": []}) + all_binary_parameters = { + parname for parnames in binary_parameters_map.values() for parname in parnames + } + + # Find all parfile parameters + all_parfile_parameters = set(parfile_dict.keys()) + + # All binary parameters in the parfile + parfile_binary_parameters = all_parfile_parameters & all_binary_parameters + + # Find which binary models include those + allowed_binary_models = { + binary_model + for (binary_model, bmc) in binary_parameters_map.items() + if len(parfile_binary_parameters - set(bmc)) == 0 + } + + # Now select the best-guess binary model + priority = [bm for bm in _binary_model_priority if bm in allowed_binary_models] + omitted = allowed_binary_models - set(priority) + + return priority + list(omitted) + + +def convert_binary_params_dict( + parfile_dict, convert_komkin=True, drop_ddk_sini=True, force_binary_model=None +): + """Convert the PINT parameter dictionary to include the best-guess binary + + Parameters + ---------- + parfile_dict + The parameter dictionary as read-in by parse_parfile + convert_komkin + Whether or not to convert the KOM and KIN parameters + drop_ddk_sini + Whether to drop SINI when converting to the DDK model + + force_binary_model : str, optional + When set to some binary model, like force_binary_model="DD", this will + override the binary model set in the parfile. Defaults to None + + Returns + ------- + A new parfile dictionary with the binary model replaced with the best-guess + model. For a conversion to DDK, this function also converts the KOM/KIN + parameters if they exist. + """ + binary = parfile_dict.get("BINARY", None) + binary = binary if not binary else binary[0] + log.debug(f"Requested to convert binary model for BINARY model: {binary}") + + if binary: + if not force_binary_model: + binary_model_guesses = guess_binary_model(parfile_dict) + log.info( + f"Compatible models with these parameters: {', '.join(binary_model_guesses)}. Using {binary_model_guesses[0]}" + ) + + 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." + ) + + log.error(log_message) + raise UnknownBinaryModel(error_message) + + else: + binary_model_guesses = [force_binary_model] + + # Select the best-guess binary model + parfile_dict["BINARY"] = [binary_model_guesses[0]] + + # 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") + 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") + 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") + parfile_dict.pop("SINI", None) + + return parfile_dict diff --git a/src/pint/models/timing_model.py b/src/pint/models/timing_model.py index bcd1c55c1..e2d6cd298 100644 --- a/src/pint/models/timing_model.py +++ b/src/pint/models/timing_model.py @@ -4070,9 +4070,17 @@ class UnknownParameter(TimingModelError): class UnknownBinaryModel(TimingModelError): - """Signal that the par file requested a binary model no in PINT.""" + """Signal that the par file requested a binary model not in PINT.""" - pass + def __init__(self, message, suggestion=None): + super().__init__(message) + self.suggestion = suggestion + + def __str__(self): + base_message = super().__str__() + if self.suggestion: + return f"{base_message} Perhaps use {self.suggestion}?" + return base_message class MissingBinaryError(TimingModelError): diff --git a/src/pint/scripts/t2binary2pint.py b/src/pint/scripts/t2binary2pint.py new file mode 100644 index 000000000..42d70f89a --- /dev/null +++ b/src/pint/scripts/t2binary2pint.py @@ -0,0 +1,57 @@ +"""PINT-based tool for converting T2 par files to PINT.""" + +import argparse + +from loguru import logger as log + +import pint.logging +from pint.models.model_builder import ModelBuilder + +pint.logging.setup(level="INFO") + +__all__ = ["main"] + + +def main(argv=None): + parser = argparse.ArgumentParser( + description="""`t2binary2pint` converts par with binary models for + Tempo2 to par files with binary models compatible with PINT. + + To guess the binary model, the model parameters in the par file are + compared to all the binary models that PINT knows about. Then, the + simplest binary model that contains all these parameters is chosen. + + The following parameters are optionally converted: + 1. KOM (AIU <--> DT96 convention) + 2. KIN (AIU <--> DT96 convention) + + The following parameters are optionally ignored: + 1. SINI (if the binary model is DDK) + + """, + 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)") + parser.add_argument( + "--convert_komkin", + type=bool, + default=True, + help="Whether to convert KOM/KIN parameters (True)", + ) + parser.add_argument( + "--drop_ddk_sini", + type=bool, + default=True, + help="Whether to drop SINI if the model is DDK (True)", + ) + + args = parser.parse_args(argv) + + mb = ModelBuilder() + + model = mb(args.input_par, allow_T2=True, allow_tcb=True) + model.write_parfile(args.output_par) + print(f"Output written to {args.output_par}") + + log.info(f"Output written to {args.output_par}.") diff --git a/tests/test_process_parfile.py b/tests/test_process_parfile.py index 8aa570250..daaa23d20 100644 --- a/tests/test_process_parfile.py +++ b/tests/test_process_parfile.py @@ -5,6 +5,7 @@ from io import StringIO from pint.models.model_builder import ModelBuilder, parse_parfile, TimingModelError +from pint.models.model_builder import guess_binary_model from pint.models import get_model @@ -157,3 +158,83 @@ def test_pintify_parfile_duplicate_binary(): m = ModelBuilder() with pytest.raises(TimingModelError): m._pintify_parfile(StringIO(dup_par)) + + +ddk_par = """ +PSR J1234+5678 +ELAT 0 1 2 +ELONG 0 1 3 +F0 1 1 5 +DM 10 0 3 +PEPOCH 57000 +BINARY T2 +PB 60 +T0 54300 +A1 30 +OM 150 +ECC 8.0e-05 +PBDOT 2.0e-14 +OMDOT 0.0001 +M2 0.3 +KOM 90 +KIN 70 +SINI KIN +#SINI 0.95 +""" + +ell1h_par = """ +PSR J1234+5678 +ELAT 0 1 2 +ELONG 0 1 3 +F0 1 1 5 +DM 10 0 3 +PEPOCH 57000 +BINARY T2 +PB 60 +A1 1 +PBDOT 2.0e-14 +TASC 54300 +EPS1 1.0e-06 +EPS2 1.0e-6 +H3 1.0e-07 +H4 1.0e-07 +""" + +nobinarymodel_par = """ +PSR J1234+5678 +ELAT 0 1 2 +ELONG 0 1 3 +F0 1 1 5 +DM 10 0 3 +PEPOCH 57000 +BINARY T2 +PB 60 +A1 1 +PBDOT 2.0e-14 +TASC 54300 +T0 54300 +EPS1 1.0e-06 +EPS2 1.0e-6 +H3 1.0e-07 +H4 1.0e-07 +""" + + +def test_guess_binary_model(): + parfiles = [base_par, ddk_par, ell1h_par, nobinarymodel_par] + binary_models = ["Isolated", "DDK", "ELL1H", None] + trip_raise = [False, True, True, True] + + for parfile, binary_model, trip in zip(parfiles, binary_models, trip_raise): + par_dict = parse_parfile(StringIO(parfile)) + + binary_model_guesses = guess_binary_model(par_dict) + + if binary_model: + assert binary_model_guesses[0] == binary_model + else: + assert binary_model_guesses == [] + + if trip: + with pytest.raises(TimingModelError): + m = get_model(StringIO(parfile)) diff --git a/tests/test_t2binary2pint.py b/tests/test_t2binary2pint.py new file mode 100644 index 000000000..36c81690c --- /dev/null +++ b/tests/test_t2binary2pint.py @@ -0,0 +1,117 @@ +"""Tests for `pint.models.tcb_conversion` and the `tcb2tdb` script.""" + +import os +import numpy as np + +import pytest + +from pint.models.model_builder import ModelBuilder +from pint.models.model_builder import _binary_model_priority +from pint.models.timing_model import UnknownBinaryModel +from pint.models.timing_model import AllComponents +from pint.scripts import t2binary2pint + + +ddkpar = """ +PSR PSRTEST +ELONG 256.66869638 1 +ELAT 30.700359811 1 +F0 218.81184039 1 +F1 -4.08278e-16 1 +PEPOCH 55636 +POSEPOCH 55636 +DMEPOCH 58983 +DM 15.988527064 1 +PMELONG 5.2676082172 1 +PMELAT -3.441494631 1 +PX 0.8981169694 1 +SINI KIN +#SINI 0.9509308076719 +BINARY T2 +PB 67.825136364 1 +T0 54303.634515 1 +A1 32.342422610 1 +OM 176.19952646 1 +ECC 7.494087e-05 1 +PBDOT 1.958980e-14 1 +OMDOT 0.0001539624 1 +M2 0.3016563544 1 +KOM 92.824444242 1 +KIN 70.946081736 1 +EPHEM DE436 +EPHVER 5 +CLK TT(BIPM2020) +T2CMETHOD IAU2000B +NE_SW 0.000 +TIMEEPH FB90 +CORRECT_TROPOSPHERE Y +""" + +incompatiblepar = """ +PSR PSRTEST +ELONG 256.66869638 1 +ELAT 30.700359811 1 +F0 218.81184039 1 +F1 -4.08278e-16 1 +PEPOCH 55636 +POSEPOCH 55636 +SINI KIN +#SINI 0.9509308076719 +BINARY T2 +PB 67.825136364 1 +T0 54303.634515 1 +A1 32.342422610 1 +OM 176.19952646 1 +EPS1 0.1237542123 1 +M2 0.3016563544 1 +KOM 92.824444242 1 +KIN 70.946081736 1 +EPHEM DE436 +EPHVER 5 +CLK TT(BIPM2020) +T2CMETHOD IAU2000B +NE_SW 0.000 +TIMEEPH FB90 +CORRECT_TROPOSPHERE Y +""" + + +def test_t2binary2pint(tmp_path): + tmppar1 = tmp_path / "tmp1.par" + tmppar2 = tmp_path / "tmp2.par" + with open(tmppar1, "w") as f: + f.write(ddkpar) + + cmd = f"{tmppar1} {tmppar2}" + t2binary2pint.main(cmd.split()) + + assert os.path.isfile(tmppar2) + + m2 = ModelBuilder()(tmppar2) + assert m2.BINARY.value == "DDK" + assert np.isclose(m2.KIN.value, 109.053918264) + assert np.isclose(m2.KOM.value, -2.8244442419) + + +def test_binary_model_priority(): + all_components = AllComponents() + binary_models = all_components.category_component_map["pulsar_system"] + binary_models = [ + all_components.components[m].binary_model_name for m in binary_models + ] + + msg = ( + "Please update _binary_model_priority in model_builder.py " + "to include all binary models" + ) + + assert set(binary_models) - set(_binary_model_priority) == set(), msg + + +def test_binary_exception(tmp_path): + tmppar = tmp_path / "tmp1.par" + with open(tmppar, "w") as f: + f.write(incompatiblepar) + + with pytest.raises(UnknownBinaryModel): + m = ModelBuilder()(tmppar)