diff --git a/CHANGELOG-unreleased.md b/CHANGELOG-unreleased.md index 6c43bdf65..f8724f8b7 100644 --- a/CHANGELOG-unreleased.md +++ b/CHANGELOG-unreleased.md @@ -15,9 +15,11 @@ the released changes. - Updated GMRT coordinates. - Replaced custom ``pint.ls`` with astropy ``u.lsec`` - Updated code to remove deprecation warnings during CI +- `d_delay_d_param` now applies the chain rule for delay components that respond to delays accumulated from earlier components (e.g. a binary delay's dependence on its evaluation epoch). Design-matrix entries for binary pulsars change at the ~1e-5 relative level. The parameter-independent ingredients are computed once per design matrix via the new `TimingModel.delay_deriv_chain`; the unused `acc_delay` argument of `d_delay_d_param` was replaced by the optional `chain` argument. ### Added - Plot whitened DM residuals in pintk. - `ssb_to_psb_xyz_ECL` and `ssb_to_psb_xyz_ICRS` are now cached +- Support for hierarchical triple systems: a second (outer) binary component can be added via a `BINARY2` line with `_2`-suffixed orbital parameters (e.g. `PB_2`, `A1_2`). Outer orbit delay is computed before, and propagated into, the inner binary. Outer wrappers `BinaryDD2`, `BinaryBT2`, and `BinaryELL12` are provided. Delay derivatives account for the outer→inner coupling (chain rule through the previous delay). The projected semi-major axis includes a second time derivative `A1DOT2` (alias `X2DOT`). ### Fixed - `WidebandTOAFitter` raises a warning if the model has correlated errors (It used to give wrong results before). - Fixed bug where "include_bipm" flag was being ignored when loading Fermi TOAs with weights, now defaults to using EPHEM, CLOCK and PLANET_SHAPIRO from the timing model @@ -36,4 +38,5 @@ the released changes. - Make VLBI frame rotation work correctly when proper motion is present. - Changed some API to pass Mac CI - Log-separated frequency computation for red noise components. +- Place ``solar_windx`` before the binary in ``DEFAULT_ORDER`` so SolarWindDispersionX delays and derivatives chain-rule through the binary the same way as ``solar_wind``. ### Removed diff --git a/docs/explanation.rst b/docs/explanation.rst index a6bd78c76..3d3d24c8c 100644 --- a/docs/explanation.rst +++ b/docs/explanation.rst @@ -588,6 +588,9 @@ components. - If the ``BINARY`` line is present in the parameter file, its value determines which binary model to use; if not, no binary model is used. + An optional ``BINARY2`` line selects an outer-orbit component for a + hierarchical triple (parameters are suffixed with ``_2``, e.g. ``PB_2``; + outer wrappers include ``BinaryDD2``, ``BinaryBT2``, and ``BinaryELL12``). - Each model component has one or more "special parameters" or families of parameters identified by a common prefix. If a par file contains a special parameter, or a known alias of one, then the timing model uses the diff --git a/docs/timingmodels.rst b/docs/timingmodels.rst index 252fd4d65..b714cf7f1 100644 --- a/docs/timingmodels.rst +++ b/docs/timingmodels.rst @@ -13,8 +13,10 @@ file, these are selected based on the parameters present. Binary models are selected explicitly using the ``BINARY`` parameter, while each non-binary component is selected if some parameter unique to it is included (for example if ``ELAT`` is present, :class:`~pint.models.astrometry.AstrometryEcliptic` is -selected). Ambiguous or contradictory parameter files are possible, and for -these PINT raises an exception. +selected). Hierarchical triples may also include an outer orbit selected with +``BINARY2`` (parameters carry a ``_2`` suffix, e.g. ``PB_2``). Ambiguous or +contradictory parameter files are possible, and for these PINT raises an +exception. .. componentlist:: diff --git a/src/pint/binaryconvert.py b/src/pint/binaryconvert.py index f300ab5a3..101bd4f4a 100644 --- a/src/pint/binaryconvert.py +++ b/src/pint/binaryconvert.py @@ -7,7 +7,7 @@ """ import copy -from typing import List, Optional, Tuple, Union +from typing import Tuple import numpy as np from astropy import units as u @@ -598,9 +598,16 @@ def convert_binary( if not model.is_binary: raise AttributeError("Input model is not a binary") - binary_component_name = [ + binary_component_names = [ x for x in model.components.keys() if x.startswith("Binary") - ][0] + ] + if len(binary_component_names) > 1: + raise ValueError( + "convert_binary does not support hierarchical triple systems " + f"with multiple binary components ({binary_component_names}); " + "convert each orbit separately or remove BINARY2 first." + ) + binary_component_name = binary_component_names[0] binary_component = model.components[binary_component_name] if binary_component.binary_model_name == output: log.debug( diff --git a/src/pint/models/__init__.py b/src/pint/models/__init__.py index d13ceec8c..bb16388ba 100644 --- a/src/pint/models/__init__.py +++ b/src/pint/models/__init__.py @@ -21,10 +21,10 @@ # Import all standard model components here from pint.models.astrometry import AstrometryEcliptic, AstrometryEquatorial -from pint.models.binary_bt import BinaryBT, BinaryBTPiecewise -from pint.models.binary_dd import BinaryDD, BinaryDDGR, BinaryDDH, BinaryDDS +from pint.models.binary_bt import BinaryBT, BinaryBT2, BinaryBTPiecewise +from pint.models.binary_dd import BinaryDD, BinaryDD2, BinaryDDGR, BinaryDDH, BinaryDDS from pint.models.binary_ddk import BinaryDDK -from pint.models.binary_ell1 import BinaryELL1, BinaryELL1H, BinaryELL1k +from pint.models.binary_ell1 import BinaryELL1, BinaryELL1H, BinaryELL1k, BinaryELL12 from pint.models.chromatic_model import ChromaticCM, ChromaticCMX from pint.models.cmwavex import CMWaveX from pint.models.dispersion_model import ( diff --git a/src/pint/models/binary_bt.py b/src/pint/models/binary_bt.py index 83651fe5d..68fa609f7 100644 --- a/src/pint/models/binary_bt.py +++ b/src/pint/models/binary_bt.py @@ -66,18 +66,47 @@ def validate(self): """Validate BT model parameters""" super().validate() for p in ("T0", "A1"): - if getattr(self, p).value is None: - raise MissingParameter("BT", p, f"{p} is required for BT") + if self._bp(p).value is None: + pname = self._bp(p).name + raise MissingParameter("BT", pname, f"{pname} is required for BT") # If any *DOT is set, we need T0 - for p in ("PBDOT", "OMDOT", "EDOT", "A1DOT"): - if getattr(self, p).value is None: - getattr(self, p).value = "0" - getattr(self, p).frozen = True + for p in ("PBDOT", "OMDOT", "EDOT", "A1DOT", "A1DOT2"): + if self._bp(p).value is None: + self._bp(p).value = "0" + self._bp(p).frozen = True - if self.GAMMA.value is None: - self.GAMMA.value = "0" - self.GAMMA.frozen = True + if self._bp("GAMMA").value is None: + self._bp("GAMMA").value = "0" + self._bp("GAMMA").frozen = True + + +class BinaryBT2(BinaryBT): + """Outer-orbit Blandford and Teukolsky model for a hierarchical triple. + + This is identical to :class:`pint.models.binary_bt.BinaryBT` except that all + of its parameters carry a ``_2`` suffix (``PB_2``, ``A1_2``, ``T0_2``, ...) + and it is selected with the ``BINARY2`` parfile parameter instead of + ``BINARY``. See :class:`pint.models.binary_dd.BinaryDD2` for a description of + how the outer orbit couples into the inner binary. + + Orbital-frequency (``FBn``) and ``ORBWAVE`` parameterizations are not + supported for the outer orbit. + + Parameters supported: + + .. paramtable:: + :class: pint.models.binary_bt.BinaryBT2 + """ + + register = True + category = "pulsar_system_outer" + param_suffix = "_2" + binary_param_tag = "BINARY2" + + def __init__(self): + super().__init__() + self._apply_param_suffix() class BinaryBTPiecewise(PulsarBinary): diff --git a/src/pint/models/binary_dd.py b/src/pint/models/binary_dd.py index 39352bbdb..6d39c8d4a 100644 --- a/src/pint/models/binary_dd.py +++ b/src/pint/models/binary_dd.py @@ -114,22 +114,60 @@ def validate(self): super().validate() self.check_required_params(["T0", "A1"]) # If any *DOT is set, we need T0 - for p in ("PBDOT", "OMDOT", "EDOT", "A1DOT"): - if hasattr(self, p) and getattr(self, p).value is None: - getattr(self, p).value = 0.0 - getattr(self, p).frozen = True + for p in ("PBDOT", "OMDOT", "EDOT", "A1DOT", "A1DOT2"): + if self._hasbp(p) and self._bp(p).value is None: + self._bp(p).value = 0.0 + self._bp(p).frozen = True - if hasattr(self, "GAMMA") and self.GAMMA.value is None: - self.GAMMA.value = 0.0 - self.GAMMA.frozen = True + if self._hasbp("GAMMA") and self._bp("GAMMA").value is None: + self._bp("GAMMA").value = 0.0 + self._bp("GAMMA").frozen = True # If eccentricity is zero, freeze some parameters to 0 # OM = 0 -> T0 = TASC - if self.ECC.value == 0 or self.ECC.value is None: + if self._bp("ECC").value == 0 or self._bp("ECC").value is None: for p in ("ECC", "OM", "OMDOT", "EDOT"): - if hasattr(self, p): - getattr(self, p).value = 0.0 - getattr(self, p).frozen = True + if self._hasbp(p): + self._bp(p).value = 0.0 + self._bp(p).frozen = True + + +class BinaryDD2(BinaryDD): + """Outer-orbit Damour and Deruelle model for a hierarchical triple system. + + This is identical to :class:`pint.models.binary_dd.BinaryDD` except that all + of its parameters carry a ``_2`` suffix (``PB_2``, ``A1_2``, ``T0_2``, ...) + and it is selected with the ``BINARY2`` parfile parameter instead of + ``BINARY``. It is intended to model the *outer* orbit of a hierarchical + triple, alongside a normal inner binary component. + + Because this component belongs to the ``pulsar_system_outer`` category, + which is ordered before ``pulsar_system`` in + :data:`pint.models.timing_model.DEFAULT_ORDER`, its delay is accumulated + before the inner binary's delay. PINT evaluates each binary at + ``barycentric time - accumulated delay``, so the outer orbit's light-travel + delay automatically shifts the epoch at which the inner orbit is evaluated. + This reproduces the physical coupling of a hierarchical triple (the wide + outer orbit Doppler-shifting the inner orbit), rather than naively adding two + independent binary delays. + + Orbital-frequency (``FBn``) and ``ORBWAVE`` parameterizations are not + supported for the outer orbit. + + Parameters supported: + + .. paramtable:: + :class: pint.models.binary_dd.BinaryDD2 + """ + + register = True + category = "pulsar_system_outer" + param_suffix = "_2" + binary_param_tag = "BINARY2" + + def __init__(self): + super().__init__() + self._apply_param_suffix() class BinaryDDS(BinaryDD): diff --git a/src/pint/models/binary_ell1.py b/src/pint/models/binary_ell1.py index d83f46c0b..d61b97b29 100644 --- a/src/pint/models/binary_ell1.py +++ b/src/pint/models/binary_ell1.py @@ -218,12 +218,16 @@ def validate(self): """Validate parameters.""" super().validate() - if self.TASC.value is None: + if self._bp("TASC").value is None: raise MissingParameter("ELL1", "TASC", "TASC is required for ELL1 model.") - for p in ["EPS1", "EPS2"]: - pm = getattr(self, p) + for p in ("EPS1", "EPS2"): + pm = self._bp(p) if pm.value is None: pm.value = 0 + for p in ("A1DOT", "A1DOT2"): + if self._hasbp(p) and self._bp(p).value is None: + self._bp(p).value = 0.0 + self._bp(p).frozen = True def change_binary_epoch(self, new_epoch): """Change the epoch for this binary model. @@ -307,6 +311,34 @@ def change_binary_epoch(self, new_epoch): return dt_integer_orbits +class BinaryELL12(BinaryELL1): + """Outer-orbit ELL1 model for a hierarchical triple system. + + This is identical to :class:`pint.models.binary_ell1.BinaryELL1` except that + all of its parameters carry a ``_2`` suffix (``PB_2``, ``A1_2``, + ``TASC_2``, ...) and it is selected with the ``BINARY2`` parfile parameter + instead of ``BINARY``. See :class:`pint.models.binary_dd.BinaryDD2` for a + description of how the outer orbit couples into the inner binary. + + Orbital-frequency (``FBn``) and ``ORBWAVE`` parameterizations are not + supported for the outer orbit. + + Parameters supported: + + .. paramtable:: + :class: pint.models.binary_ell1.BinaryELL12 + """ + + register = True + category = "pulsar_system_outer" + param_suffix = "_2" + binary_param_tag = "BINARY2" + + def __init__(self): + super().__init__() + self._apply_param_suffix() + + class BinaryELL1H(BinaryELL1): """ELL1 modified to use H3 parameter for Shapiro delay. diff --git a/src/pint/models/chromatic_model.py b/src/pint/models/chromatic_model.py index e1f3a327c..e6335e53f 100644 --- a/src/pint/models/chromatic_model.py +++ b/src/pint/models/chromatic_model.py @@ -7,9 +7,9 @@ from loguru import logger as log from pint import DMconst -from pint.exceptions import MissingParameter +from pint.exceptions import MissingParameter, MissingTOAs from pint.models.parameter import MJDParameter, floatParameter, prefixParameter -from pint.models.timing_model import DelayComponent, MissingParameter, MissingTOAs +from pint.models.timing_model import DelayComponent from pint.toa_select import TOASelect from pint.utils import split_prefixed_name, taylor_horner, taylor_horner_deriv diff --git a/src/pint/models/model_builder.py b/src/pint/models/model_builder.py index 84a8f4c13..9216c415e 100644 --- a/src/pint/models/model_builder.py +++ b/src/pint/models/model_builder.py @@ -271,8 +271,8 @@ def _validate_components(self): f" class, please set register to 'False' in the class" f" of component {k}." ) - if v.category == "pulsar_system": - # The pulsar system will be selected by parameter BINARY + if v.category in ("pulsar_system", "pulsar_system_outer"): + # The pulsar system is selected by parameter BINARY/BINARY2 continue else: raise ComponentConflict(m) @@ -488,6 +488,12 @@ def choose_model(self, param_inpar, force_binary_model=None, allow_T2=False): self.choose_binary_model(param_inpar, force_binary_model, allow_T2) ) + # Outer-orbit binary for hierarchical triple systems. + binary2 = param_inpar.get("BINARY2", None) + if binary2: + binary2 = binary2[0] + selected_components.add(self.choose_outer_binary_model(param_inpar)) + # 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 @@ -531,6 +537,17 @@ def choose_model(self, param_inpar, force_binary_model=None, allow_T2=False): ) else: continue + # Outer-orbit binary component, controlled by the BINARY2 tag. + if self.all_components.components[cps[0]].category == "pulsar_system_outer": + if binary2 is None: + raise MissingBinaryError( + f"The outer-orbit binary model is decided by the" + f" parameter 'BINARY2'. Please indicate the outer" + f" binary model before using parameter {k}, which" + f" is an outer binary model parameter." + ) + else: + continue if len(cps) == 1: # No conflict, parameter only shows in one component. selected_components.add(cps[0]) @@ -645,6 +662,32 @@ def choose_binary_model(self, param_inpar, force_binary_model=None, allow_T2=Fal return binary_cp.__class__.__name__ + def choose_outer_binary_model(self, param_inpar): + """Choose the outer-orbit BINARY2 model for a hierarchical triple. + + Parameters + ---------- + param_inpar: dict + Dictionary of the unique parameters in .par file with the key being + the parfile line. :func:`parse_parfile` returns this dictionary. + + Returns + ------- + str + Name of the outer binary component class. + + Note + ---- + The outer model is taken verbatim from the ``BINARY2`` parameter (no T2 + guessing is performed) and must correspond to a registered component in + the ``pulsar_system_outer`` category (e.g. ``DD`` -> ``BinaryDD2``). + """ + binary2 = param_inpar["BINARY2"][0] + binary_cp = self.all_components.search_binary_components( + binary2, category="pulsar_system_outer" + ) + return binary_cp.__class__.__name__ + def _setup_model( self, timing_model, diff --git a/src/pint/models/pulsar_binary.py b/src/pint/models/pulsar_binary.py index d9b50460b..01472e84f 100644 --- a/src/pint/models/pulsar_binary.py +++ b/src/pint/models/pulsar_binary.py @@ -49,6 +49,7 @@ class PulsarBinary(DelayComponent): - PBDOT - time derivative of binary period (s/s) - A1 - projected orbital amplitude, $a \sin i$ (ls, non-negative) - A1DOT - time derivative of projected orbital amplitude (ls/s) + - A1DOT2 - Second time derivative of projected orbital amplitude (ls/s^2) - ECC (or E) - eccentricity (no units, 0<=ECC<1) - EDOT - time derivative of eccentricity (1/s) - OM - longitude of periastron (deg) @@ -87,6 +88,16 @@ class PulsarBinary(DelayComponent): category = "pulsar_system" + # Suffix appended to the PINT-facing parameter names of this component + # (e.g. ``"_2"`` for the outer orbit of a hierarchical triple). The + # standalone binary instance always uses the canonical (unsuffixed) names. + # An empty string means the normal single-binary behaviour. + param_suffix = "" + + # The top-level parameter that selects this binary model in the parfile. + # Outer-orbit components override this with ``"BINARY2"``. + binary_param_tag = "BINARY" + def __init__(self): super().__init__() self.binary_model_name = None @@ -134,6 +145,16 @@ def __init__(self): tcb2tdb_scale_factor=(1 / consts.c), ) ) + self.add_param( + floatParameter( + name="A1DOT2", + aliases=["X2DOT"], + units=u.lsec / u.s**2, + description="Second Derivative of projected semi-major axis, d2[ap*sin(i)]/dt2", + unit_scale=False, + tcb2tdb_scale_factor=(1 / consts.c), + ) + ) self.add_param( floatParameter( name="ECC", @@ -258,6 +279,69 @@ def __init__(self): self.warn_default_params = ["ECC", "OM"] # Set up delay function self.delay_funcs_component += [self.binarymodel_delay] + self.delay_deriv_wrt_prev_delay_funcs += [self.d_binary_delay_d_prev_delay] + + def _apply_param_suffix(self): + """Rename this component's PINT-facing parameters with ``param_suffix``. + + This is used by outer-orbit components (e.g. the outer binary of a + hierarchical triple) so that their parameters appear in the parfile as + ``PB_2``, ``A1_2``, ... while the underlying standalone binary instance + keeps the canonical names (``PB``, ``A1``, ...). + + The canonical attribute (e.g. ``self.PB``) is *removed* from this + component so that it does not leak into the parent + :class:`~pint.models.timing_model.TimingModel` namespace (where it would + collide with the inner binary's identically named parameter). Base-class + methods therefore access parameters through :meth:`_bp` / :meth:`_hasbp`, + which apply the suffix. + + Prefix parameters (``FBn``, ``ORBWAVECn``, ...) are removed entirely: + orbital-frequency and ORBWAVE parameterizations are not supported for a + suffixed (outer) orbit. + """ + suffix = self.param_suffix + if not suffix: + return + for canonical in list(self.params): + par = getattr(self, canonical) + if getattr(par, "is_prefix", False): + # Orbital-frequency / ORBWAVE parameterizations are not + # supported for a suffixed (outer) orbit. + self.remove_param(canonical) + continue + new_name = canonical + suffix + par.aliases = [alias + suffix for alias in par.aliases] + par.name = new_name + # Keep any funcParameter cross-references consistent with the + # renamed parameters (no-op for DD/BT, which have none). + if hasattr(par, "_params"): + par._params = [p + suffix for p in par._params] + # Expose the renamed parameter under its suffixed name only and drop + # the canonical attribute so it cannot leak to the parent model. + setattr(self, new_name, par) + delattr(self, canonical) + self.params[self.params.index(canonical)] = new_name + + def _bp(self, name): + """Return this component's parameter object for canonical ``name``. + + Applies :attr:`param_suffix` so that base-class code written in terms of + canonical names (e.g. ``PB``) resolves to the suffixed parameter (e.g. + ``PB_2``) for an outer-orbit component. For a normal binary + (``param_suffix == ""``) this is just ``getattr(self, name)``. + """ + suffixed = name + self.param_suffix + if self.param_suffix and hasattr(self, suffixed): + return getattr(self, suffixed) + return getattr(self, name) + + def _hasbp(self, name): + """Whether this component has the (possibly suffixed) parameter ``name``.""" + suffixed = name + self.param_suffix + return bool(self.param_suffix and hasattr(self, suffixed)) or hasattr( + self, name + ) def setup(self): super().setup() @@ -378,57 +462,73 @@ def setup(self): def validate(self): super().validate() if ( - hasattr(self, "SINI") - and self.SINI.value is not None - and not 0 <= self.SINI.value <= 1 + self._hasbp("SINI") + and self._bp("SINI").value is not None + and not 0 <= self._bp("SINI").value <= 1 ): raise ValueError( - f"Sine of inclination angle must be between zero and one ({self.SINI.quantity})" + f"Sine of inclination angle must be between zero and one ({self._bp('SINI').quantity})" ) - if hasattr(self, "M2") and self.M2.value is not None and self.M2.value < 0: + if ( + self._hasbp("M2") + and self._bp("M2").value is not None + and self._bp("M2").value < 0 + ): raise ValueError( - f"Companion mass M2 cannot be negative ({self.M2.quantity})" + f"Companion mass M2 cannot be negative ({self._bp('M2').quantity})" ) if ( - hasattr(self, "ECC") - and self.ECC.value is not None - and not 0 <= self.ECC.value <= 1 + self._hasbp("ECC") + and self._bp("ECC").value is not None + and not 0 <= self._bp("ECC").value <= 1 ): raise ValueError( - f"Eccentricity ECC must be between zero and one ({self.ECC.quantity})" + f"Eccentricity ECC must be between zero and one ({self._bp('ECC').quantity})" ) - if self.A1.value is not None and self.A1.value < 0: + if ( + self._hasbp("A1") + and self._bp("A1").value is not None + and self._bp("A1").value < 0 + ): raise ValueError( - f"Projected semi-major axis A1 cannot be negative ({self.A1.quantity})" + f"Projected semi-major axis A1 cannot be negative ({self._bp('A1').quantity})" ) - if self.PB.value is not None: - if self.PB.value <= 0: + has_fb0 = self._hasbp("FB0") + if self._hasbp("PB") and self._bp("PB").value is not None: + if self._bp("PB").value <= 0: raise ValueError( - f"Binary period PB must be non-negative ({self.PB.quantity})" + f"Binary period PB must be non-negative ({self._bp('PB').quantity})" + ) + if ( + has_fb0 + and self._bp("FB0").value is not None + and not ( + isinstance(self._bp("FB0"), funcParameter) + or isinstance(self._bp("PB"), funcParameter) ) - if self.FB0.value is not None and not ( - isinstance(self.FB0, funcParameter) - or isinstance(self.PB, funcParameter) ): raise ValueError("Model cannot have values for both FB0 and PB") - if self.FB0.value is not None and self.FB0.value <= 0: + if has_fb0 and self._bp("FB0").value is not None and self._bp("FB0").value <= 0: raise ValueError( - f"Binary frequency FB0 must be non-negative ({self.FB0.quantity})" + f"Binary frequency FB0 must be non-negative ({self._bp('FB0').quantity})" ) def check_required_params(self, required_params): # search for all the possible to get the parameters. for p in required_params: - par = getattr(self, p) + par = self._bp(p) if par.value is None: # try to search if there is any class method that computes it method_name = f"{p.lower()}_func" try: par_method = getattr(self.binary_instance, method_name) except AttributeError as e: + # Use the (possibly suffixed) parameter name so outer-orbit + # components report e.g. T0_2 rather than the canonical T0. raise MissingParameter( self.binary_model_name, - f"{p} is required for '{self.binary_model_name}'.", + par.name, + f"{par.name} is required for '{self.binary_model_name}'.", ) from e par_method() @@ -508,21 +608,31 @@ def update_binary_object(self, toas, acc_delay=None): epoch=tbl["tdbld"].astype(np.float64), ecl=self._parent.ECL.value ) for par in self.binary_instance.binary_params: + # The standalone binary instance uses canonical names (``par``), + # while this component's PINT-facing parameter may carry a suffix + # (e.g. ``PB_2`` for an outer orbit). ``param_suffix`` is empty for + # a normal single binary, in which case this is a no-op. + pint_par_name = par + self.param_suffix if par in self.binary_instance.param_aliases.keys(): - alias = self.binary_instance.param_aliases[par] + alias = [ + a + self.param_suffix + for a in self.binary_instance.param_aliases[par] + ] else: alias = [] # the _parent attribute should give access to all the components - if hasattr(self._parent, par) or set(alias).intersection(self.params): + if hasattr(self._parent, pint_par_name) or set(alias).intersection( + self.params + ): try: - pint_bin_name = self._parent.match_param_aliases(par) + pint_bin_name = self._parent.match_param_aliases(pint_par_name) except UnknownParameter as e: if par in self.internal_params: - pint_bin_name = par + pint_bin_name = pint_par_name else: raise UnknownParameter( - f"Unable to find {par} in the parent model" + f"Unable to find {pint_par_name} in the parent model" ) from e binObjpar = getattr(self._parent, pint_bin_name) @@ -555,19 +665,31 @@ def binarymodel_delay(self, toas, acc_delay=None): def d_binary_delay_d_xxxx(self, toas, param, acc_delay): """Return the binary model delay derivatives.""" self.update_binary_object(toas, acc_delay) + # The standalone binary instance only knows the canonical (unsuffixed) + # parameter names, so strip the suffix before requesting a derivative. + if self.param_suffix and param.endswith(self.param_suffix): + param = param[: -len(self.param_suffix)] return self.binary_instance.d_binarydelay_d_par(param) + def d_binary_delay_d_prev_delay(self, toas, acc_delay): + """Return the derivative of the binary delay w.r.t. the accumulated + delay from preceding components (dimensionless).""" + self.update_binary_object(toas, acc_delay) + return self.binary_instance.d_binarydelay_d_prevdelay() + def print_par(self, format="pint"): + tag = self.binary_param_tag + tag_par = getattr(self._parent, tag) if self._parent is not None else None if self._parent is None: - result = f"BINARY {self.binary_model_name}\n" - elif self._parent.BINARY.value != self.binary_model_name: + result = f"{tag} {self.binary_model_name}\n" + elif tag_par.value != self.binary_model_name: raise TimingModelError( - f"Parameter BINARY {self._parent.BINARY.value}" + f"Parameter {tag} {tag_par.value}" f" does not match the binary" f" component {self.binary_model_name}" ) else: - result = self._parent.BINARY.as_parfile_line(format=format) + result = tag_par.as_parfile_line(format=format) for p in self.params: par = getattr(self, p) @@ -617,32 +739,38 @@ def change_binary_epoch(self, new_epoch): """ new_epoch = parse_time(new_epoch, scale="tdb", precision=9) + # Parameter access goes through _bp() so that this works for both a + # normal (inner) binary and a suffixed outer-orbit component. + PB_par = self._bp("PB") + PBDOT_par = self._bp("PBDOT") + T0_par = self._bp("T0") + # Get PB and PBDOT from model - if self.PB.quantity is not None and not isinstance(self.PB, funcParameter): - PB = self.PB.quantity - if self.PBDOT.quantity is not None: - PBDOT = self.PBDOT.quantity + if PB_par.quantity is not None and not isinstance(PB_par, funcParameter): + PB = PB_par.quantity + if PBDOT_par.quantity is not None: + PBDOT = PBDOT_par.quantity else: PBDOT = 0.0 * u.Unit("") else: - PB = 1.0 / self.FB0.quantity + PB = 1.0 / self._bp("FB0").quantity try: - PBDOT = -self.FB1.quantity / self.FB0.quantity**2 + PBDOT = -self._bp("FB1").quantity / self._bp("FB0").quantity ** 2 except AttributeError: PBDOT = 0.0 * u.Unit("") # Find closest periapsis time and reassign T0 - t0_ld = self.T0.quantity.tdb.mjd_long + t0_ld = T0_par.quantity.tdb.mjd_long dt = (new_epoch.tdb.mjd_long - t0_ld) * u.day d_orbits = dt / PB - PBDOT * dt**2 / (2.0 * PB**2) n_orbits = np.round(d_orbits.to(u.Unit(""))) if n_orbits == 0: return dt_integer_orbits = PB * n_orbits + PB * PBDOT * n_orbits**2 / 2.0 - self.T0.quantity = self.T0.quantity + dt_integer_orbits + T0_par.quantity = T0_par.quantity + dt_integer_orbits with contextlib.suppress(AttributeError): - if self.FB2.quantity is not None: + if self._bp("FB2").quantity is not None: log.warning( "Ignoring orbital frequency derivatives higher than FB1" "in computing new T0; a model fit should resolve this" @@ -650,23 +778,23 @@ def change_binary_epoch(self, new_epoch): # Update PB or FB0, FB1, etc. if isinstance(self.binary_instance.orbits_cls, bo.OrbitPB): dPB = PBDOT * dt_integer_orbits - self.PB.quantity = self.PB.quantity + dPB + PB_par.quantity = PB_par.quantity + dPB else: fbterms = [0.0 * u.Unit("")] + self._parent.get_prefix_list("FB") for n in range(len(fbterms) - 1): - cur_deriv = getattr(self, f"FB{n}") + cur_deriv = self._bp(f"FB{n}") cur_deriv.value = taylor_horner_deriv( dt_integer_orbits.to(u.s), fbterms, deriv_order=n + 1 ) # Update ECC, OM, and A1 - dECC = self.EDOT.quantity * dt_integer_orbits - self.ECC.quantity = self.ECC.quantity + dECC - dOM = self.OMDOT.quantity * dt_integer_orbits - self.OM.quantity = self.OM.quantity + dOM - dA1 = self.A1DOT.quantity * dt_integer_orbits - self.A1.quantity = self.A1.quantity + dA1 + dECC = self._bp("EDOT").quantity * dt_integer_orbits + self._bp("ECC").quantity = self._bp("ECC").quantity + dECC + dOM = self._bp("OMDOT").quantity * dt_integer_orbits + self._bp("OM").quantity = self._bp("OM").quantity + dOM + dA1 = self._bp("A1DOT").quantity * dt_integer_orbits + self._bp("A1").quantity = self._bp("A1").quantity + dA1 def pb(self, t=None): """Return binary period and uncertainty (optionally evaluated at different times) regardless of binary model @@ -684,22 +812,23 @@ def pb(self, t=None): Binary period uncertainty """ + PB_par = self._bp("PB") + PBDOT_par = self._bp("PBDOT") if self.binary_model_name.startswith("ELL1"): - t0 = self.TASC.quantity + t0 = self._bp("TASC").quantity else: - t0 = self.T0.quantity + t0 = self._bp("T0").quantity t = t0 if t is None else parse_time(t) - if self.PB.quantity is not None: - if self.PBDOT.quantity is None and ( - not hasattr(self, "XPBDOT") - or getattr(self, "XPBDOT").quantity is not None + if PB_par.quantity is not None: + if PBDOT_par.quantity is None and ( + not self._hasbp("XPBDOT") or self._bp("XPBDOT").quantity is not None ): - return self.PB.quantity, self.PB.uncertainty - pb = self.PB.as_ufloat(u.d) - if self.PBDOT.quantity is not None: - pbdot = self.PBDOT.as_ufloat(u.s / u.s) - if hasattr(self, "XPBDOT") and self.XPBDOT.quantity is not None: - pbdot += self.XPBDOT.as_ufloat(u.s / u.s) + return PB_par.quantity, PB_par.uncertainty + pb = PB_par.as_ufloat(u.d) + if PBDOT_par.quantity is not None: + pbdot = PBDOT_par.as_ufloat(u.s / u.s) + if self._hasbp("XPBDOT") and self._bp("XPBDOT").quantity is not None: + pbdot += self._bp("XPBDOT").as_ufloat(u.s / u.s) pnew = pb + pbdot * (t - t0).jd if not isinstance(pnew, np.ndarray): return pnew.n * u.d, pnew.s * u.d if pnew.s > 0 else None @@ -710,7 +839,7 @@ def pb(self, t=None): uncertainties.unumpy.std_devs(pnew) * u.d, ) - elif self.FB0.quantity is not None: + elif self._hasbp("FB0") and self._bp("FB0").quantity is not None: # assume FB terms dt = (t - t0).sec coeffs = [] diff --git a/src/pint/models/solar_wind_dispersion.py b/src/pint/models/solar_wind_dispersion.py index ff8147117..112b33dc9 100644 --- a/src/pint/models/solar_wind_dispersion.py +++ b/src/pint/models/solar_wind_dispersion.py @@ -10,6 +10,7 @@ import pint.utils from pint import DMconst +from pint.exceptions import MissingTOAs from pint.models.dispersion_model import Dispersion from pint.models.parameter import ( MJDParameter, @@ -17,7 +18,6 @@ intParameter, prefixParameter, ) -from pint.models.timing_model import MissingTOAs from pint.toa_select import TOASelect diff --git a/src/pint/models/stand_alone_psr_binaries/BT_model.py b/src/pint/models/stand_alone_psr_binaries/BT_model.py index feefcceb1..17ab2a4fc 100644 --- a/src/pint/models/stand_alone_psr_binaries/BT_model.py +++ b/src/pint/models/stand_alone_psr_binaries/BT_model.py @@ -168,6 +168,12 @@ def d_delayL2_d_A1(self): def d_delayL2_d_A1DOT(self): return self.tt0 * self.d_delayL2_d_A1() + def d_delayL1_d_A1DOT2(self): + return 0.5 * self.tt0**2 * self.d_delayL1_d_A1() + + def d_delayL2_d_A1DOT2(self): + return 0.5 * self.tt0**2 * self.d_delayL2_d_A1() + def d_delayL1_d_OM(self): a1 = self.a1() / c.c return a1 * np.cos(self.omega()) * (np.cos(self.E()) - self.ecc()) @@ -207,10 +213,22 @@ def d_delayL2_d_GAMMA(self): return np.sin(self.E()) def d_delayL1_d_T0(self): - return self.d_delayL1_d_E() * self.d_E_d_T0() + # Include a1(t) secular terms (A1DOT, A1DOT2): required for the + # outer→inner prev-delay chain rule in hierarchical triples. + with u.set_enabled_equivalencies(u.dimensionless_angles()): + d_a1_d_T0 = self.prtl_der("a1", "T0") + return ( + self.d_delayL1_d_E() * self.d_E_d_T0() + + self.d_delayL1_d_A1() * d_a1_d_T0 + ) def d_delayL2_d_T0(self): - return self.d_delayL2_d_E() * self.d_E_d_T0() + with u.set_enabled_equivalencies(u.dimensionless_angles()): + d_a1_d_T0 = self.prtl_der("a1", "T0") + return ( + self.d_delayL2_d_E() * self.d_E_d_T0() + + self.d_delayL2_d_A1() * d_a1_d_T0 + ) def d_delayL1_d_par(self, par): if par not in self.binary_params: diff --git a/src/pint/models/stand_alone_psr_binaries/DD_model.py b/src/pint/models/stand_alone_psr_binaries/DD_model.py index 604397dee..f7ae12cac 100644 --- a/src/pint/models/stand_alone_psr_binaries/DD_model.py +++ b/src/pint/models/stand_alone_psr_binaries/DD_model.py @@ -349,6 +349,18 @@ def d_beta_d_A1DOT(self): d_a1_d_A1DOT = self.d_a1_d_A1DOT() return d_a1_d_A1DOT / c.c * (1 - eTheta**2) ** 0.5 * cosOmg + def d_beta_d_A1DOT2(self): + """Derivative. + + Computes:: + + dBeta/dA1DOT2 = * d_a1_d_A1DOT2/c*(1-eTheta**2)**0.5*cos(omega) + """ + eTheta = self.eTheta() + cosOmg = np.cos(self.omega()) + d_a1_d_A1DOT2 = self.d_a1_d_A1DOT2() + return d_a1_d_A1DOT2 / c.c * (1 - eTheta**2) ** 0.5 * cosOmg + def d_beta_d_T0(self): """Derivative. diff --git a/src/pint/models/stand_alone_psr_binaries/ELL1_model.py b/src/pint/models/stand_alone_psr_binaries/ELL1_model.py index f88365019..b803180c5 100644 --- a/src/pint/models/stand_alone_psr_binaries/ELL1_model.py +++ b/src/pint/models/stand_alone_psr_binaries/ELL1_model.py @@ -34,6 +34,7 @@ def __init__(self): self.ELL1_interVars = ["eps1", "eps2", "Phi", "Dre", "Drep", "Drepp", "nhat"] self.add_inter_vars(self.ELL1_interVars) self.orbits_func = self.orbits_ELL1 + self.d_binarydelay_d_prev_delay_par = "TASC" @property def tt0(self): @@ -51,24 +52,32 @@ def ttasc(self): return (t - self.TASC.value * u.day).to("second") def a1(self): - """ELL1 model a1 calculation. + """ELL1 projected semi-major axis a1(t). - This method overrides the a1() method in pulsar_binary.py. Instead of tt0, - it uses ttasc. + Overrides the generic ``a1()`` to use ``ttasc`` instead of ``tt0``. + Includes A1DOT2 so hierarchical-triple prev-delay derivatives are correct. """ - return self.A1 + self.ttasc() * self.A1DOT + return ( + self.A1 + self.ttasc() * self.A1DOT + 0.5 * self.ttasc() ** 2 * self.A1DOT2 + ) def d_a1_d_A1(self): return np.longdouble(np.ones(len(self.ttasc()))) * u.Unit("") + def d_a1_d_TASC(self): + # Critical for d_binarydelay_d_prevdelay (par == "TASC"). + return (-self.A1DOT - self.A1DOT2 * self.ttasc()).to(self.A1DOT.unit) + def d_a1_d_T0(self): - result = np.empty(len(self.ttasc())) - result.fill(-self.A1DOT.value) - return result * u.Unit(self.A1DOT.unit) + # Kept for compatibility; ELL1's epoch parameter is TASC. + return self.d_a1_d_TASC() def d_a1_d_A1DOT(self): return self.ttasc() + def d_a1_d_A1DOT2(self): + return 0.5 * self.ttasc() ** 2 + def eps1(self): return self.EPS1 + self.ttasc() * self.EPS1DOT diff --git a/src/pint/models/stand_alone_psr_binaries/binary_generic.py b/src/pint/models/stand_alone_psr_binaries/binary_generic.py index d6d61b298..e2adb9bec 100644 --- a/src/pint/models/stand_alone_psr_binaries/binary_generic.py +++ b/src/pint/models/stand_alone_psr_binaries/binary_generic.py @@ -98,6 +98,7 @@ def __init__( "EDOT": 0.0 / u.second, "A1": 10.0 * u.lsec, "A1DOT": 0.0 * u.lsec / u.second, + "A1DOT2": 0.0 * u.lsec / u.second**2, "T0": np.longdouble(54000.0) * u.day, "OM": 0.0 * u.deg, "OMDOT": 0.0 * u.deg / u.year, @@ -122,6 +123,7 @@ def __init__( self.cache_vars = ["E", "nu"] self.binary_delay_funcs = [] self.d_binarydelay_d_par_funcs = [] + self.d_binarydelay_d_prev_delay_par = "T0" self.orbits_cls = OrbitPB(self, ["PB", "PBDOT", "XPBDOT", "T0"]) @property @@ -262,6 +264,22 @@ def d_binarydelay_d_par(self, par): return result + def d_binarydelay_d_prevdelay(self): + """Get the derivative of the binary delay w.r.t. delays accumulated + from previous components + + Returns + ------- + astropy.units.Quantity + Dimensionless derivative of the binary delay w.r.t. earlier delays. + The binary delay depends on time only through ``t - T0`` (or + ``t - TASC``), so the epoch-parameter derivative equals + ``-d(delay)/dt``, which is the response to an earlier delay. + """ + + with u.set_enabled_equivalencies(u.dimensionless_angles()): + return self.d_binarydelay_d_par(self.d_binarydelay_d_prev_delay_par).to("") + def prtl_der(self, y, x): """Find the partial derivatives in binary model pdy/pdx @@ -394,19 +412,26 @@ def d_ecc_d_EDOT(self): return self.tt0 def a1(self): - return self.A1 + self.tt0 * self.A1DOT if hasattr(self, "_tt0") else self.A1 + + if hasattr(self, "_tt0"): + return self.A1 + self.tt0 * self.A1DOT + 0.5 * self.tt0**2 * self.A1DOT2 + else: + return self.A1 def d_a1_d_A1(self): return np.longdouble(np.ones(len(self.tt0))) * u.Unit("") def d_a1_d_T0(self): - result = np.empty(len(self.tt0)) - result.fill(-self.A1DOT.value) - return result * u.Unit(self.A1DOT.unit) + + result = -self.A1DOT - self.A1DOT2 * self.tt0 + return result.to(self.A1DOT.unit) def d_a1_d_A1DOT(self): return self.tt0 + def d_a1_d_A1DOT2(self): + return 0.5 * self.tt0**2 + def d_a1_d_par(self, par): if par not in self.binary_params: errorMesg = f"{par} is not in binary parameter list." @@ -531,7 +556,6 @@ def d_E_d_par(self, par): else: E = self.E() return np.zeros(len(self.tt0)) * E.unit / par_obj.unit - return func() def nu(self): """True anomaly (Ae)""" @@ -807,6 +831,9 @@ def d_Pobs_d_OMDOT(self): def d_Pobs_d_A1DOT(self): return self.tt0 * self.d_Pobs_d_A1() + def d_Pobs_d_A1DOT2(self): + return 0.5 * self.tt0**2 * self.d_Pobs_d_A1() + ############## Calculation for design matrix ################ def Pobs_designmatrix(self, params): npars = len(params) diff --git a/src/pint/models/timing_model.py b/src/pint/models/timing_model.py index e8930d039..c65ebfc40 100644 --- a/src/pint/models/timing_model.py +++ b/src/pint/models/timing_model.py @@ -34,14 +34,23 @@ import inspect from collections import OrderedDict, defaultdict from functools import wraps -from typing import Callable, Dict, List, Literal, Optional, Set, Tuple, Union +from typing import ( + Callable, + Dict, + List, + Literal, + NamedTuple, + Optional, + Set, + Tuple, + Union, +) from warnings import warn import astropy.coordinates as coords import astropy.time as time import numpy as np -from astropy import constants as c, units as u -from astropy.table import Table +from astropy import units as u from astropy.utils.decorators import lazyproperty from loguru import logger as log from scipy.optimize import brentq @@ -51,8 +60,6 @@ from pint.derived_quantities import dispersion_slope from pint.exceptions import ( AliasConflict, - MissingBinaryError, - MissingParameter, MissingTOAs, PrefixError, PropertyAttributeError, @@ -81,7 +88,6 @@ colorize, open_or_use, split_prefixed_name, - xxxselections, get_unit, ) @@ -122,9 +128,11 @@ "troposphere", "solar_system_shapiro", "solar_wind", + "solar_windx", "dispersion_constant", "dispersion_dmx", "dispersion_jump", + "pulsar_system_outer", "pulsar_system", "frequency_dependent", "absolute_phase", @@ -135,6 +143,23 @@ ] +class DelayDerivChain(NamedTuple): + """Parameter-independent ingredients of delay-derivative calculations. + + Produced by :meth:`TimingModel.delay_deriv_chain`; all lists are ordered + like ``TimingModel.DelayComponent_list``. + """ + + acc_delays: List[u.Quantity] + """Delay accumulated from all components preceding each component.""" + + chain_factors: List[u.Quantity] + """Per component, ``1 + d(component delay)/d(accumulated delay)``.""" + + total_delay: u.Quantity + """The total delay, equal to ``TimingModel.delay(toas)``.""" + + def property_exists(f): """Mark a function as a property but handle AttributeErrors. @@ -322,6 +347,14 @@ def __init__(self, name: str = "", components: List["Component"] = []): ), "", ) + self.add_param_from_top( + strParameter( + name="BINARY2", + description="Outer-orbit binary model for a hierarchical triple system", + value=None, + ), + "", + ) self.add_param_from_top( boolParameter( name="DILATEFREQ", @@ -489,14 +522,33 @@ def num_components_of_type(type): from pint.models.pulsar_binary import PulsarBinary + def num_binaries_with_tag(tag): + return len( + list( + filter( + lambda c: isinstance(c, PulsarBinary) + and getattr(c, "binary_param_tag", "BINARY") == tag, + self.components.values(), + ) + ) + ) + has_binary_attr = hasattr(self, "BINARY") and self.BINARY.value if has_binary_attr: assert ( - num_components_of_type(PulsarBinary) == 1 - ), "BINARY attribute is set but no PulsarBinary component found." + num_binaries_with_tag("BINARY") == 1 + ), "BINARY attribute is set but no (inner) PulsarBinary component found." + has_binary2_attr = hasattr(self, "BINARY2") and self.BINARY2.value + if has_binary2_attr: + assert ( + num_binaries_with_tag("BINARY2") == 1 + ), "BINARY2 attribute is set but no outer PulsarBinary component found." assert ( - num_components_of_type(PulsarBinary) <= 1 - ), "Model can have at most one PulsarBinary component." + num_binaries_with_tag("BINARY") <= 1 + ), "Model can have at most one inner PulsarBinary component." + assert ( + num_binaries_with_tag("BINARY2") <= 1 + ), "Model can have at most one outer PulsarBinary component." from pint.models.solar_wind_dispersion import ( SolarWindDispersion, @@ -925,6 +977,26 @@ def is_binary(self) -> bool: return any(isinstance(x, PulsarBinary) for x in self.components.values()) + def _get_primary_binary_component(self): + """Return the binary component tagged by ``BINARY`` when present. + + For hierarchical triples we want orbital utility methods to operate on + the inner orbit (the component selected by the ``BINARY`` line), not the + outer ``BINARY2`` component. + """ + from pint.models.pulsar_binary import PulsarBinary + + binaries = [c for c in self.components.values() if isinstance(c, PulsarBinary)] + if not binaries: + return None + + for b in binaries: + if getattr(b, "binary_param_tag", "BINARY") == "BINARY": + return b + + # Fallback for legacy/single-binary behavior. + return binaries[0] + def orbital_phase( self, barytimes: Union[time.Time, TOAs, np.ndarray, float, MJDParameter], @@ -963,10 +1035,7 @@ def orbital_phase( """ if not self.is_binary: # punt if not a binary return None - # Find the binary model - b = self.components[ - [x for x in self.components.keys() if x.startswith("Binary")][0] - ] + b = self._get_primary_binary_component() # Make sure that the binary instance has the binary params b.update_binary_object(None) # Handle input times and update them in stand-alone binary models @@ -1034,9 +1103,7 @@ def pulsar_radial_velocity( """ # this should also update the binary instance nu = self.orbital_phase(barytimes, anom="true") - b = self.components[ - [x for x in self.components.keys() if x.startswith("Binary")][0] - ] + b = self._get_primary_binary_component() bbi = b.binary_instance # shorthand psi = nu + bbi.omega() return ( @@ -1108,10 +1175,7 @@ def conjunction(self, baryMJD: Union[float, time.Time]) -> Union[float, np.ndarr """ if not self.is_binary: # punt if not a binary return None - # Find the binary model - b = self.components[ - [x for x in self.components.keys() if x.startswith("Binary")][0] - ] + b = self._get_primary_binary_component() bbi = b.binary_instance # shorthand # Superior conjunction occurs when true anomaly + omega == 90 deg # We will need to solve for this using a root finder (brentq) @@ -1133,7 +1197,7 @@ def funct(t): scs = [] for bt in bts: # Make 11 times over one orbit after bt - pb = self.pb()[0].to_value("day") + pb = b.pb()[0].to_value("day") ts = np.linspace(bt, bt + pb, 11) # Compute the true anomalies and omegas for those times nus = self.orbital_phase(ts, anom="true") @@ -2156,7 +2220,13 @@ def d_phase_d_tpulsar(self, toas: TOAs): """ raise NotImplementedError - def d_phase_d_param(self, toas: TOAs, delay: u.Quantity, param: str) -> u.Quantity: + def d_phase_d_param( + self, + toas: TOAs, + delay: u.Quantity, + param: str, + chain: Optional[DelayDerivChain] = None, + ) -> u.Quantity: """Return the derivative of phase with respect to the parameter. This is the derivative of the phase observed at each TOA with @@ -2179,6 +2249,11 @@ def d_phase_d_param(self, toas: TOAs, delay: u.Quantity, param: str) -> u.Quanti the value should be ``self.delay(toas)``. param : str The name of the parameter to differentiate with respect to. + chain : DelayDerivChain, optional + Precomputed output of :meth:`delay_deriv_chain` for these TOAs, + forwarded to :meth:`d_delay_d_param`. It is parameter-independent, + so callers evaluating derivatives for many parameters should + compute it once and pass it in. Returns ------- @@ -2189,7 +2264,7 @@ def d_phase_d_param(self, toas: TOAs, delay: u.Quantity, param: str) -> u.Quanti # Is it safe to assume that any param affecting delay only affects # phase indirectly (and vice-versa)?? if delay is None: - delay = self.delay(toas) + delay = chain.total_delay if chain is not None else self.delay(toas) par = getattr(self, param) result = np.longdouble(np.zeros(toas.ntoas)) / par.units phase_derivs = self.phase_deriv_funcs @@ -2205,29 +2280,91 @@ def d_phase_d_param(self, toas: TOAs, delay: u.Quantity, param: str) -> u.Quanti # d_Phase2/d_delay*d_delay/d_param # = (d_Phase1/d_delay + d_Phase2/d_delay) * # d_delay_d_param - d_delay_d_p = self.d_delay_d_param(toas, param) + d_delay_d_p = self.d_delay_d_param(toas, param, chain=chain) dpdd_result = np.longdouble(np.zeros(toas.ntoas)) / u.second for dpddf in self.d_phase_d_delay_funcs: dpdd_result += dpddf(toas, delay) result = dpdd_result * d_delay_d_p return result.to(result.unit, equivalencies=u.dimensionless_angles()) + def delay_deriv_chain(self, toas: TOAs) -> DelayDerivChain: + """Compute the parameter-independent parts of delay derivatives. + + In a single pass over ``self.DelayComponent_list`` this computes, for + each delay component, the delay accumulated from all preceding + components (which sets the epoch at which the component is evaluated) + and the chain-rule factor ``1 + d(component delay)/d(accumulated + delay)`` describing how the component's delay responds to a change in + the preceding delays (e.g. a binary delay responds to a shift of its + evaluation epoch; most components do not depend on it at all and get a + factor of one). + + These quantities do not depend on the parameter being differentiated, + so callers evaluating derivatives for many parameters (notably + :meth:`designmatrix`) should compute this once and pass it to + :meth:`d_delay_d_param` / :meth:`d_phase_d_param`. + + Parameters + ---------- + toas : pint.toa.TOAs + The TOAs at which delays and derivatives are evaluated. + + Returns + ------- + DelayDerivChain + """ + acc_delays = [] + chain_factors = [] + delay = np.zeros(toas.ntoas) * u.second + for cp in self.DelayComponent_list: + acc_delays.append(delay) + factor = 1 + for f in cp.delay_deriv_wrt_prev_delay_funcs: + factor = factor + f(toas, delay).to(u.dimensionless_unscaled) + chain_factors.append(factor) + # Accumulate out of place so the stored acc_delays stay intact. + for df in cp.delay_funcs_component: + delay = delay + df(toas, delay) + return DelayDerivChain(acc_delays, chain_factors, delay) + def d_delay_d_param( - self, toas: TOAs, param: str, acc_delay: Optional[u.Quantity] = None + self, toas: TOAs, param: str, chain: Optional[DelayDerivChain] = None ) -> u.Quantity: - """Return the derivative of delay with respect to the parameter.""" + """Return the derivative of delay with respect to the parameter. + + Parameters + ---------- + toas : pint.toa.TOAs + The TOAs at which the derivative should be evaluated. + param : str + The name of the parameter to differentiate with respect to. + chain : DelayDerivChain, optional + Precomputed output of :meth:`delay_deriv_chain` for these TOAs. + It is parameter-independent, so callers evaluating derivatives + for many parameters should compute it once and pass it in; + it is computed internally when not provided. + """ par = getattr(self, param) - result = np.longdouble(np.zeros(toas.ntoas) << (u.s / par.units)) - delay_derivs = self.delay_deriv_funcs - if param not in list(delay_derivs.keys()): + if param not in self.delay_deriv_funcs: raise AttributeError( f"Derivative function for '{param}' is not provided" f" or not registered; parameter '{param}' may not be fittable. " ) - for df in delay_derivs[param]: - result += df(toas, param, acc_delay).to( - result.unit, equivalencies=u.dimensionless_angles() - ) + if chain is None: + chain = self.delay_deriv_chain(toas) + + result = np.longdouble(np.zeros(toas.ntoas) << (u.s / par.units)) + for cp, acc_delay, factor in zip( + self.DelayComponent_list, chain.acc_delays, chain.chain_factors + ): + # Derivatives of the delays accumulated so far propagate through + # this component's dependence on its evaluation epoch. + result *= factor + if param in cp.deriv_funcs: + for df in cp.deriv_funcs[param]: + result += df(toas, param, acc_delay).to( + result.unit, equivalencies=u.dimensionless_angles() + ) return result def d_phase_d_param_num( @@ -2415,7 +2552,8 @@ def designmatrix( F0 = self.F0.quantity # 1/sec ntoas = len(toas) nparams = len(params) - delay = self.delay(toas) + chain = self.delay_deriv_chain(toas) + delay = chain.total_delay units = [] # Apply all delays ? # tt = toas['tdbld'] @@ -2428,7 +2566,7 @@ def designmatrix( M[:, ii] = 1.0 / F0.value units.append(u.s / u.s) else: - q = -self.d_phase_d_param(toas, delay, param) + q = -self.d_phase_d_param(toas, delay, param, chain=chain) the_unit = u.Unit("") / getattr(self, param).units M[:, ii] = q.to_value(the_unit) / F0.value units.append(the_unit / F0.unit) @@ -3135,7 +3273,8 @@ def as_parfile( if format.lower() == "tempo2": result_begin += "MODE 1\n" for p in self.top_level_params: - if p == "BINARY": # Will print the Binary model name in the binary section + # Will print the binary model name in the (outer) binary section + if p in ("BINARY", "BINARY2"): continue result_begin += getattr(self, p).as_parfile_line(format=format) for cat in start_order: @@ -3467,9 +3606,9 @@ def get_derived_params( outdict["Dist (pc)"] = 1.0 / px # Now binary system derived parameters if self.is_binary: - for x in self.components: - if x.startswith("Binary"): - binary = x + # Prefer the BINARY-tagged (inner) component so hierarchical triples + # do not report BinaryDD2 / BinaryBT2 here. + binary = self._get_primary_binary_component().__class__.__name__ s += f"\nBinary model {binary}\n" outdict["Binary"] = binary @@ -4013,6 +4152,7 @@ class DelayComponent(Component): def __init__(self): super().__init__() self.delay_funcs_component = [] + self.delay_deriv_wrt_prev_delay_funcs = [] class PhaseComponent(Component): @@ -4225,13 +4365,20 @@ def component_unique_params(self) -> Dict[str, List[str]]: component_special_params[cps[0]].append(param) return component_special_params - def search_binary_components(self, system_name: str) -> "Component": + def search_binary_components( + self, system_name: str, category: str = "pulsar_system" + ) -> "Component": """Search the pulsar binary component based on given name. Parameters ---------- system_name : str Searching name for the pulsar binary/system + category : str, optional + The component category to search within. Defaults to + ``"pulsar_system"`` (the inner binary). Use + ``"pulsar_system_outer"`` to find the outer-orbit component of a + hierarchical triple. Return ------ @@ -4243,7 +4390,7 @@ def search_binary_components(self, system_name: str) -> "Component": If the input binary model name does not match any PINT defined binary model. """ - all_systems = self.category_component_map["pulsar_system"] + all_systems = self.category_component_map[category] if system_name in all_systems: return self.components[system_name] for cp_name in all_systems: diff --git a/src/pint/pint_matrix.py b/src/pint/pint_matrix.py index cca7eace0..5c082189e 100644 --- a/src/pint/pint_matrix.py +++ b/src/pint/pint_matrix.py @@ -449,7 +449,8 @@ def __call__(self, data, model, derivative_params, offset=True, offset_padding=1 M = np.zeros((data.ntoas, len(params))) labels = [{self.derivative_quantity: (0, M.shape[0], self.quantity_unit)}] labels_dim2 = {} - delay = model.delay(data) + chain = model.delay_deriv_chain(data) + delay = chain.total_delay for ii, param in enumerate(params): if param == "Offset": M[:, ii] = offset_padding @@ -457,7 +458,9 @@ def __call__(self, data, model, derivative_params, offset=True, offset_padding=1 else: param_unit = getattr(model, param).units # Since this is the phase derivative, we know the quantity unit. - q = deriv_func(data, delay, param).to(u.Unit("") / param_unit) + q = deriv_func(data, delay, param, chain=chain).to( + u.Unit("") / param_unit + ) # NOTE Here we have negative sign here. Since in pulsar timing # the residuals are calculated as (Phase - int(Phase)), which is different diff --git a/tests/datafile/B1855+09_triple_DD.par b/tests/datafile/B1855+09_triple_DD.par new file mode 100644 index 000000000..346917fc1 --- /dev/null +++ b/tests/datafile/B1855+09_triple_DD.par @@ -0,0 +1,36 @@ +PSRJ 1855+09 +RAJ 18:57:36.3932884 1 0.00002602730280675029 +DECJ +09:43:17.29196 1 0.00078789485676919773 +F0 186.49408156698235146 1 0.00000000000698911818 +F1 -6.2049547277487420583e-16 1 1.7380934373573401505e-20 +PEPOCH 49453 +POSEPOCH 49453 +DMEPOCH 49453 +DM 13.29709 +PMRA -2.5054345161030380639 1 0.03104958261053317181 +PMDEC -5.4974558631993817232 1 0.06348008663748286318 +PX 1.2288569063263405232 1 0.21243361289239687251 +SINI 0.99741717335200923866 1 0.00182023515130851988 +BINARY DD +PB 12.327171194774200418 1 0.00000000079493185824 +T0 49452.940695077335647 1 0.00169031830532837251 +A1 9.2307804312998001928 1 0.00000036890718667634 +OM 276.55142180589701234 1 0.04936551005019605698 +ECC 2.1745265668236919017e-05 1 0.00000004027191312623 +M2 0.26111312480723428917 1 0.02616161008932908066 +BINARY2 DD +PB_2 1400.0 1 +T0_2 49452.0 1 +A1_2 120.0 1 +OM_2 110.0 1 +ECC_2 0.3 1 +START 53358.726464889485214 +FINISH 55108.922917417192366 +TZRMJD 54177.508359343262555 +TZRFRQ 424 +TZRSITE ao +TRES 0.395 +CLK TT(TAI) +MODE 1 +UNITS TDB +EPHEM DE405 diff --git a/tests/test_binary_generic.py b/tests/test_binary_generic.py index dd7d5e352..c32ba91d7 100644 --- a/tests/test_binary_generic.py +++ b/tests/test_binary_generic.py @@ -5,12 +5,18 @@ import pytest from pint.models.model_builder import get_model -from pint.models.timing_model import MissingParameter +from pint.exceptions import MissingParameter from utils import verify_stand_alone_binary_parameter_updates from pinttestdata import datadir -bad_trouble = ["J1923+2515_NANOGrav_9yv1.gls.par", "J1744-1134.basic.ecliptic.par"] +bad_trouble = [ + "J1923+2515_NANOGrav_9yv1.gls.par", + "J1744-1134.basic.ecliptic.par", + # Hierarchical triple (two binary components); this generic single-binary + # check cannot disambiguate the inner vs outer parameter sets. + "B1855+09_triple_DD.par", +] @pytest.mark.parametrize("parfile", glob(join(datadir, "*.par"))) diff --git a/tests/test_chromatic.py b/tests/test_chromatic.py index 227553ab7..b0808484b 100644 --- a/tests/test_chromatic.py +++ b/tests/test_chromatic.py @@ -3,7 +3,7 @@ import pytest from pint.models import get_model, get_model_and_toas from pint.models.chromatic_model import ChromaticCM -from pint.models.timing_model import MissingParameter +from pint.exceptions import MissingParameter from pint.simulation import make_fake_toas_uniform from pint.fitter import WLSFitter import astropy.units as u diff --git a/tests/test_ddk.py b/tests/test_ddk.py index e9401cc6d..d32171ebe 100644 --- a/tests/test_ddk.py +++ b/tests/test_ddk.py @@ -19,7 +19,7 @@ import pint.simulation import pint.toa as toa from pint.models.parameter import boolParameter -from pint.models.timing_model import MissingParameter +from pint.exceptions import MissingParameter from pint.residuals import Residuals import pint.fitter diff --git a/tests/test_fermiphase.py b/tests/test_fermiphase.py index 4aa62f3d3..734cc5312 100644 --- a/tests/test_fermiphase.py +++ b/tests/test_fermiphase.py @@ -55,7 +55,7 @@ def test_process_and_accuracy(): # level by comparison with stored Tempo2 "Fermi plugin" results. modelin = pint.models.get_model(parfile) - get_satellite_observatory("Fermi", ft2file) + get_satellite_observatory("Fermi", ft2file, overwrite=True) ts = get_Fermi_TOAs( eventfileraw, weightcolumn="PSRJ0030+0451", diff --git a/tests/test_fitter_error_checking.py b/tests/test_fitter_error_checking.py index bcdc7f3ac..b352c0247 100644 --- a/tests/test_fitter_error_checking.py +++ b/tests/test_fitter_error_checking.py @@ -7,7 +7,7 @@ import pint.fitter from pint.models import get_model -from pint.models.timing_model import MissingTOAs +from pint.exceptions import MissingTOAs from pint.simulation import make_fake_toas_uniform par_base = """ diff --git a/tests/test_model_manual.py b/tests/test_model_manual.py index de406b7b4..8ffc5fb60 100644 --- a/tests/test_model_manual.py +++ b/tests/test_model_manual.py @@ -5,15 +5,14 @@ import pytest import astropy.units as u +from pint.exceptions import MissingParameter, UnknownBinaryModel from pint.models.astrometry import AstrometryEquatorial from pint.models.dispersion_model import DispersionDM, DispersionDMX from pint.models.spindown import Spindown from pint.models.model_builder import get_model from pint.models.timing_model import ( - MissingParameter, TimingModel, Component, - UnknownBinaryModel, ) from pinttestdata import datadir diff --git a/tests/test_orbit_phase.py b/tests/test_orbit_phase.py index cc73fc4f3..04307916c 100644 --- a/tests/test_orbit_phase.py +++ b/tests/test_orbit_phase.py @@ -1,5 +1,5 @@ -import pytest import os +import io import pytest import numpy as np @@ -91,3 +91,41 @@ def test_j0737(self): assert len(x) == 2, "conjunction is not returning an array" # make sure true anomaly before T0 is positive assert mJ0737.orbital_phase(52000.0, anom="true").value > 0.0 + + def test_triple_orbital_utilities_use_inner_binary(self): + with open(os.path.join(datadir, "B1855+09_triple_DD.par")) as f: + triple_lines = f.readlines() + triple_model = m.get_model(io.StringIO("".join(triple_lines))) + + inner_only_par = "".join( + line + for line in triple_lines + if line.split() + and line.split()[0] != "BINARY2" + and not line.split()[0].endswith("_2") + ) + inner_model = m.get_model(io.StringIO(inner_only_par)) + + # The timing-model orbital utilities should follow the BINARY-tagged + # (inner) component, not the outer BINARY2 component. + ts = triple_model.T0.value + np.linspace(0, triple_model.PB.value, 32) + triple_phase = triple_model.orbital_phase(ts, anom="mean", radians=False) + inner_phase = inner_model.orbital_phase(ts, anom="mean", radians=False) + assert np.allclose(triple_phase, inner_phase) + assert np.allclose( + triple_model.pulsar_radial_velocity(ts), + inner_model.pulsar_radial_velocity(ts), + ) + assert np.isclose( + triple_model.conjunction(triple_model.T0.value), + inner_model.conjunction(inner_model.T0.value), + ) + + # Sanity check that the outer model's mean anomaly is generally different. + outer = triple_model.components["BinaryDD2"] + outer.update_binary_object(None) + outer.binary_instance.update_input(barycentric_toa=np.asarray(ts)) + outer_phase = np.remainder(outer.binary_instance.M().value, 2 * np.pi) / ( + 2 * np.pi + ) + assert not np.allclose(triple_phase, outer_phase) diff --git a/tests/test_parfile_writing_format.py b/tests/test_parfile_writing_format.py index 222243b3d..5f3ad589f 100644 --- a/tests/test_parfile_writing_format.py +++ b/tests/test_parfile_writing_format.py @@ -75,14 +75,17 @@ def test_STIGMA(): def test_A1DOT(): """Should get changed to XDOT for TEMPO/TEMPO2""" m = get_model(os.path.join(datadir, "J1600-3053_test.par")) - assert ( - "A1DOT" in m.as_parfile() - and "XDOT" not in m.as_parfile() - and "A1DOT" not in m.as_parfile(format="tempo") - and "XDOT" in m.as_parfile(format="tempo") - and "A1DOT" not in m.as_parfile(format="tempo2") - and "XDOT" in m.as_parfile(format="tempo2") - ) + + def names(fmt="pint"): + return { + line.split()[0] + for line in m.as_parfile(format=fmt).splitlines() + if line.strip() and not line.startswith("#") + } + + assert "A1DOT" in names() and "XDOT" not in names() + assert "A1DOT" not in names("tempo") and "XDOT" in names("tempo") + assert "A1DOT" not in names("tempo2") and "XDOT" in names("tempo2") def test_ECL(): diff --git a/tests/test_piecewise.py b/tests/test_piecewise.py index 0809a6328..955488a41 100644 --- a/tests/test_piecewise.py +++ b/tests/test_piecewise.py @@ -9,7 +9,7 @@ import pint.residuals import pint.toa from pinttestdata import datadir -from pint.models.timing_model import MissingParameter +from pint.exceptions import MissingParameter from pint import fitter parfile = os.path.join(datadir, "piecewise.par") diff --git a/tests/test_triple_binary.py b/tests/test_triple_binary.py new file mode 100644 index 000000000..6a31f772b --- /dev/null +++ b/tests/test_triple_binary.py @@ -0,0 +1,369 @@ +"""Tests for hierarchical triple systems (two binary components). + +A hierarchical triple is modelled with a normal inner binary (``BINARY``, +parameters ``PB``, ``A1``, ...) plus an outer-orbit binary (``BINARY2``, +parameters ``PB_2``, ``A1_2``, ...). The outer-orbit component belongs to the +``pulsar_system_outer`` category, which is ordered before ``pulsar_system`` in +:data:`pint.models.timing_model.DEFAULT_ORDER`, so its delay is accumulated +first and propagated into the inner binary's evaluation epoch. +""" + +import io +import os + +import astropy.units as u +import numpy as np +import pytest +from pinttestdata import datadir + +import pint.models.model_builder as mb +import pint.simulation as sim +from pint.models.binary_bt import BinaryBT2 +from pint.models.binary_dd import BinaryDD, BinaryDD2 +from pint.models.binary_ell1 import BinaryELL12 +from pint.residuals import Residuals + +TRIPLE_PAR = os.path.join(datadir, "B1855+09_triple_DD.par") +TRIPLE_PAR_DD = TRIPLE_PAR + +TRIPLE_PAR_BT = """\ +PSRJ J1737_triple_BT +RAJ 17:37:47.11235 +DECJ -08:11:08.887 +F0 239.51996484444 +F1 -4.55E-16 +PEPOCH 54987 +DM 55.311 +BINARY BT +PB 79.517379 +ECC 5.38E-5 +A1 9.332791 +T0 54696.879781933 +OM 49.8 +BINARY2 BT +PB_2 1400.0 +T0_2 54696.0 +A1_2 120.0 +OM_2 110.0 +ECC_2 0.3 +TZRMJD 54987 +TZRFRQ 1400 +TZRSITE @ +CLK TT(TAI) +UNITS TDB +EPHEM DE405 +""" + +TRIPLE_PAR_ELL1 = """\ +PSRJ J0023_triple_ELL1 +ELONG 9.07039380 +ELAT 6.30910853 +F0 327.8470205906107 +F1 -1.22783E-15 +PEPOCH 56567 +DM 14.32810 +BINARY ELL1 +PB 0.138799 +A1 0.03484142 +TASC 56567.02609362 +EPS1 7.2E-6 +EPS2 -4.0E-6 +BINARY2 ELL1 +PB_2 1400.0 +A1_2 120.0 +TASC_2 56567.0 +EPS1_2 0.01 +EPS2_2 0.02 +TZRMJD 56567 +TZRFRQ 1400 +TZRSITE @ +CLK TT(TAI) +UNITS TDB +EPHEM DE436 +""" + +FAMILY = { + "DD": { + "par": TRIPLE_PAR_DD, + "params": ["A1DOT2", "A1DOT", "T0", "A1_2", "T0_2"], + }, + "BT": { + "par": TRIPLE_PAR_BT, + "params": ["A1DOT2", "A1DOT", "T0", "A1_2", "T0_2"], + }, + "ELL1": { + "par": TRIPLE_PAR_ELL1, + "params": ["A1DOT2", "A1DOT", "TASC", "A1_2", "TASC_2"], + }, +} + +STEPS = { + "A1DOT2": 1e-22, + "A1DOT": 1e-15, + "T0": 1e-6, + "TASC": 1e-6, + "A1_2": 1e-4, + "T0_2": 1e-4, + "TASC_2": 1e-4, +} + + +def _inner_only_par(): + """Return the triple parfile text with the BINARY2/outer lines removed.""" + lines = [] + with open(TRIPLE_PAR) as f: + for line in f: + key = line.split()[0] if line.split() else "" + if key == "BINARY2" or key.endswith("_2"): + continue + lines.append(line) + return "".join(lines) + + +def _load_family_model(family): + par = FAMILY[family]["par"] + if family == "DD": + return mb.get_model(par) + return mb.get_model(io.StringIO(par)) + + +def _family_toas(model): + return sim.make_fake_toas_uniform( + model.PEPOCH.value - 200, + model.PEPOCH.value + 800, + 50, + model, + freq=1400 * u.MHz, + add_noise=False, + ) + + +@pytest.fixture(scope="module") +def triple_model(): + return mb.get_model(TRIPLE_PAR) + + +@pytest.fixture(scope="module") +def toas(triple_model): + return sim.make_fake_toas_uniform( + 53400, 55000, 50, triple_model, freq=1400 * u.MHz, add_noise=False + ) + + +def test_two_binary_components_built(triple_model): + """Both an inner and an outer binary component are present.""" + assert "BinaryDD" in triple_model.components + assert "BinaryDD2" in triple_model.components + assert triple_model.components["BinaryDD"].category == "pulsar_system" + assert triple_model.components["BinaryDD2"].category == "pulsar_system_outer" + assert triple_model.BINARY.value == "DD" + assert triple_model.BINARY2.value == "DD" + + +def test_outer_ordered_before_inner(triple_model): + """The outer binary must be evaluated before the inner one so that its + delay propagates into the inner orbit.""" + order = [c.__class__.__name__ for c in triple_model.DelayComponent_list] + assert order.index("BinaryDD2") < order.index("BinaryDD") + + +def test_parameters_resolve_to_correct_component(triple_model): + """Canonical names resolve to the inner binary and ``_2`` names to the + outer binary, with no cross-contamination.""" + assert np.isclose(triple_model.PB.quantity.to_value(u.day), 12.327171194774200418) + assert triple_model.PB_2.quantity == 1400.0 * u.day + assert np.isclose(triple_model.A1_2.value, 120.0) + assert triple_model.OM_2.quantity == 110.0 * u.deg + # The outer component exposes suffixed names only (canonical names removed). + outer = triple_model.components["BinaryDD2"] + assert "PB_2" in outer.params + assert "PB" not in outer.params + assert not hasattr(outer, "PB") + + +def test_parfile_roundtrip(triple_model): + """``BINARY2`` and the ``_2`` parameters survive a parfile round-trip.""" + s = triple_model.as_parfile() + par_lines = [line.split() for line in s.splitlines() if line.split()] + assert ["BINARY", "DD"] in par_lines + assert ["BINARY2", "DD"] in par_lines + assert any(parts[0] == "PB_2" for parts in par_lines) + assert any(parts[0] == "A1_2" for parts in par_lines) + # BINARY/BINARY2 should each appear exactly once. + assert sum(parts[0] == "BINARY" for parts in par_lines) == 1 + assert sum(parts[0] == "BINARY2" for parts in par_lines) == 1 + + m2 = mb.get_model(io.StringIO(s)) + assert m2.PB_2.quantity == triple_model.PB_2.quantity + assert m2.A1_2.quantity == triple_model.A1_2.quantity + assert m2.BINARY2.value == "DD" + + +def test_residuals_finite(triple_model, toas): + res = Residuals(toas, triple_model).time_resids + assert np.all(np.isfinite(res.value)) + + +def test_outer_orbit_affects_delay(triple_model, toas): + """Switching the outer orbit on/off changes the total delay substantially.""" + d_on = triple_model.delay(toas) + m_off = mb.get_model(TRIPLE_PAR) + m_off.A1_2.value = 0.0 + d_off = m_off.delay(toas) + # A1_2 = 120 ls means the outer Roemer delay reaches ~100 s. + assert np.max(np.abs((d_on - d_off).to_value(u.s))) > 1.0 + + +def test_outer_delay_propagates_into_inner(triple_model, toas): + """The defining feature of a hierarchical triple: the inner binary is + evaluated at an epoch shifted by the outer orbit's light-travel delay, + rather than the two binaries being treated independently.""" + m_inner = mb.get_model(io.StringIO(_inner_only_par())) + + # Trigger a delay computation so each inner binary instance caches the + # barycentric time it was evaluated at. + triple_model.delay(toas) + m_inner.delay(toas) + + t_triple = triple_model.components["BinaryDD"].binary_instance.t + t_alone = m_inner.components["BinaryDD"].binary_instance.t + + shift = (t_triple - t_alone).to_value(u.s) + # The inner orbit's evaluation epoch is shifted by the outer delay (seconds + # scale), which is exactly the coupling that cures the apparent PBDOT etc. + assert np.max(np.abs(shift)) > 1.0 + + +def test_naive_sum_differs_from_coupled(triple_model, toas): + """The coupled triple delay differs from naively adding an independent + inner-binary delay and an independent outer-binary delay.""" + inner_comp = triple_model.components["BinaryDD"] + outer_comp = triple_model.components["BinaryDD2"] + + # Coupled: inner sees the accumulated outer delay. + coupled_total = triple_model.delay(toas) + + # Naive sum: evaluate each binary at the same (outer-free) accumulated delay. + acc_before = triple_model.delay( + toas, cutoff_component="BinaryDD2", include_last=False + ) + outer_only = outer_comp.binarymodel_delay(toas, acc_before) + inner_only = inner_comp.binarymodel_delay(toas, acc_before) + naive_total = acc_before + outer_only + inner_only + + diff = np.max(np.abs((coupled_total - naive_total).to_value(u.s))) + assert diff > 0.0 + + +def test_outer_param_derivative(triple_model, toas): + """Derivatives with respect to an outer (suffixed) parameter are available + and non-trivial, so the outer orbit is fittable.""" + d = triple_model.d_delay_d_param(toas, "A1_2") + assert np.all(np.isfinite(d.value)) + assert np.any(d.value != 0) + + +def test_outer_wrapper_classes(): + """The outer wrappers are configured for the BINARY2 tag and _2 suffix.""" + for cls in (BinaryDD2, BinaryBT2, BinaryELL12): + outer = cls() + assert outer.category == "pulsar_system_outer" + assert outer.param_suffix == "_2" + assert outer.binary_param_tag == "BINARY2" + assert "PB_2" in outer.params + assert "PB" not in outer.params + + ell = BinaryELL12() + assert "TASC_2" in ell.params + assert "TASC" not in ell.params + + # The inner DD model is unchanged. + inner = BinaryDD() + assert inner.category == "pulsar_system" + assert inner.param_suffix == "" + assert "PB" in inner.params + + +def test_derived_params_report_inner_binary(triple_model): + """Summary / derived-parameter text should name the inner BINARY component.""" + text, info = triple_model.get_derived_params(returndict=True) + assert info["Binary"] == "BinaryDD" + assert "Binary model BinaryDD" in text + assert "BinaryDD2" not in text.split("Binary model")[1].splitlines()[0] + + +def test_convert_binary_rejects_triple(triple_model): + """convert_binary must not silently operate on a hierarchical triple.""" + import pint.binaryconvert + + with pytest.raises(ValueError, match="multiple binary components"): + pint.binaryconvert.convert_binary(triple_model, "BT") + + +def test_outer_missing_parameter_uses_suffixed_name(): + """MissingParameter messages for outer orbits should cite T0_2, not T0.""" + from pint.exceptions import MissingParameter + + outer = BinaryDD2() + # Leave T0_2 unset; give A1_2 a value so only T0_2 is reported missing. + outer.A1_2.value = 1.0 + with pytest.raises(MissingParameter, match="T0_2") as exc_info: + outer.validate() + assert exc_info.value.param == "T0_2" + + +def test_a1dot2_changes_delay(toas): + """A nonzero A1DOT2 (second derivative of the projected semi-major axis) + changes the inner-binary delay.""" + m = mb.get_model(TRIPLE_PAR) + d0 = m.delay(toas) + m.A1DOT2.quantity = 1e-18 * u.lsec / u.s**2 + d1 = m.delay(toas) + # 0.5 * A1DOT2 * tt0**2 with tt0 up to ~1e8 s gives a delay change of + # order milliseconds; just require a clearly nonzero effect. + assert np.max(np.abs((d1 - d0).to_value(u.s))) > 1e-9 + + +@pytest.mark.parametrize( + "family,param", + [(fam, p) for fam, cfg in FAMILY.items() for p in cfg["params"]], +) +def test_delay_derivatives_match_numerical(family, param): + """Analytic delay derivatives (including A1DOT2 and the chain rule through + the outer->inner delay coupling) agree with central finite differences.""" + m = _load_family_model(family) + toas = _family_toas(m) + m.A1DOT.quantity = 3e-13 * u.lsec / u.s + m.A1DOT2.quantity = 2e-21 * u.lsec / u.s**2 + + ana = m.d_delay_d_param(toas, param) + q = getattr(m, param) + v0, h = q.value, STEPS[param] + q.value = v0 + h + dp = m.delay(toas) + q.value = v0 - h + dm = m.delay(toas) + q.value = v0 + num = (dp - dm) / (2 * h * q.units) + a = ana.to_value(num.unit) + n = num.value + scale = np.max(np.abs(n)) + assert scale > 0 + assert np.max(np.abs(a - n)) / scale < 1e-4 + + +@pytest.mark.parametrize("family", ["BT", "ELL1"]) +def test_outer_orbit_affects_delay_family(family): + par = FAMILY[family]["par"] + m = mb.get_model(io.StringIO(par)) + toas = sim.make_fake_toas_uniform( + m.PEPOCH.value - 200, + m.PEPOCH.value + 800, + 40, + m, + freq=1400 * u.MHz, + ) + d0 = m.delay(toas) + m.A1_2.quantity = 0 * u.lsec + d1 = m.delay(toas) + assert np.max(np.abs((d1 - d0).to_value(u.s))) > 1e-9 diff --git a/tests/test_wavex.py b/tests/test_wavex.py index 7c5c7f9f5..2a60decc1 100644 --- a/tests/test_wavex.py +++ b/tests/test_wavex.py @@ -6,7 +6,7 @@ from astropy import units as u from pint.models import get_model from pint.models import model_builder as mb -from pint.models.timing_model import MissingParameter +from pint.exceptions import MissingParameter from pint.fitter import Fitter from pint.residuals import Residuals from pint.simulation import make_fake_toas_uniform