From 3d0919432580510319a16e41341c6b3a3b3c7994 Mon Sep 17 00:00:00 2001 From: Rutger van Haasteren Date: Tue, 16 Jun 2026 06:40:11 +0000 Subject: [PATCH 01/16] Support hierarchical triple systems with two binary components Enable a second (outer) binary delay component in a PINT timing model so a pulsar binary can itself orbit a third body. The outer orbit's delay is computed before, and propagated into, the inner binary calculation. - Make PulsarBinary suffix-aware (param_suffix, binary_param_tag, _bp/_hasbp helpers, _apply_param_suffix) so a single implementation serves as both the inner (canonical params) and outer (_2-suffixed params) component. - Route pb() and change_binary_epoch() through _bp()/_hasbp() so the inherited public methods work on suffixed outer components. - Add BinaryDD2 and BinaryBT2 outer wrappers (category pulsar_system_outer, param_suffix _2, BINARY2 tag). - Register BINARY2, order pulsar_system_outer before pulsar_system, and handle BINARY2 selection / validation in TimingModel and ModelBuilder. - Add triple-system par file and tests. --- CHANGELOG-unreleased.md | 1 + src/pint/models/__init__.py | 4 +- src/pint/models/binary_bt.py | 42 ++++- src/pint/models/binary_dd.py | 58 +++++-- src/pint/models/model_builder.py | 47 +++++- src/pint/models/pulsar_binary.py | 228 +++++++++++++++++++------- src/pint/models/timing_model.py | 50 +++++- tests/datafile/B1855+09_triple_DD.par | 36 ++++ tests/test_binary_generic.py | 8 +- tests/test_triple_binary.py | 179 ++++++++++++++++++++ 10 files changed, 564 insertions(+), 89 deletions(-) create mode 100644 tests/datafile/B1855+09_triple_DD.par create mode 100644 tests/test_triple_binary.py diff --git a/CHANGELOG-unreleased.md b/CHANGELOG-unreleased.md index aeeee5a7dc..ea284c4e7d 100644 --- a/CHANGELOG-unreleased.md +++ b/CHANGELOG-unreleased.md @@ -16,6 +16,7 @@ the released changes. ### 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` and `BinaryBT2` are provided. ### 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 diff --git a/src/pint/models/__init__.py b/src/pint/models/__init__.py index d13ceec8c9..b48670b4da 100644 --- a/src/pint/models/__init__.py +++ b/src/pint/models/__init__.py @@ -21,8 +21,8 @@ # 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.chromatic_model import ChromaticCM, ChromaticCMX diff --git a/src/pint/models/binary_bt.py b/src/pint/models/binary_bt.py index e055bd3993..51a5e81768 100644 --- a/src/pint/models/binary_bt.py +++ b/src/pint/models/binary_bt.py @@ -67,18 +67,46 @@ def validate(self): """Validate BT model parameters""" super().validate() for p in ("T0", "A1"): - if getattr(self, p).value is None: + if self._bp(p).value is None: raise MissingParameter("BT", p, f"{p} is required for BT") # If any *DOT is set, we need T0 for p in ("PBDOT", "OMDOT", "EDOT", "A1DOT"): - if getattr(self, p).value is None: - getattr(self, p).value = "0" - getattr(self, p).frozen = True + 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 39352bbdb8..17420d9682 100644 --- a/src/pint/models/binary_dd.py +++ b/src/pint/models/binary_dd.py @@ -115,21 +115,59 @@ def validate(self): 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 + 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/model_builder.py b/src/pint/models/model_builder.py index 84a8f4c130..9216c415e5 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 79b2800b5a..92ea36c051 100644 --- a/src/pint/models/pulsar_binary.py +++ b/src/pint/models/pulsar_binary.py @@ -88,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 @@ -260,6 +270,68 @@ def __init__(self): # Set up delay function self.delay_funcs_component += [self.binarymodel_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() for bpar in self.params: @@ -379,48 +451,61 @@ 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" @@ -509,21 +594,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) @@ -556,19 +651,25 @@ 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 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) @@ -618,32 +719,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" @@ -651,23 +758,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 @@ -685,22 +792,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 @@ -711,7 +819,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/timing_model.py b/src/pint/models/timing_model.py index fa4484bc2e..2491a3af53 100644 --- a/src/pint/models/timing_model.py +++ b/src/pint/models/timing_model.py @@ -125,6 +125,7 @@ "dispersion_constant", "dispersion_dmx", "dispersion_jump", + "pulsar_system_outer", "pulsar_system", "frequency_dependent", "absolute_phase", @@ -322,6 +323,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 +498,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_binaries_with_tag("BINARY") <= 1 + ), "Model can have at most one inner PulsarBinary component." assert ( - num_components_of_type(PulsarBinary) <= 1 - ), "Model can have at most one PulsarBinary component." + num_binaries_with_tag("BINARY2") <= 1 + ), "Model can have at most one outer PulsarBinary component." from pint.models.solar_wind_dispersion import ( SolarWindDispersion, @@ -3133,7 +3161,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: @@ -4223,13 +4252,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 ------ @@ -4241,7 +4277,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/tests/datafile/B1855+09_triple_DD.par b/tests/datafile/B1855+09_triple_DD.par new file mode 100644 index 0000000000..346917fc19 --- /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 dd7d5e3521..4caac397db 100644 --- a/tests/test_binary_generic.py +++ b/tests/test_binary_generic.py @@ -10,7 +10,13 @@ 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_triple_binary.py b/tests/test_triple_binary.py new file mode 100644 index 0000000000..fb781111a5 --- /dev/null +++ b/tests/test_triple_binary.py @@ -0,0 +1,179 @@ +"""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 + +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.residuals import Residuals +from pinttestdata import datadir + +TRIPLE_PAR = os.path.join(datadir, "B1855+09_triple_DD.par") + + +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) + + +@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): + 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 + + # The inner DD model is unchanged. + inner = BinaryDD() + assert inner.category == "pulsar_system" + assert inner.param_suffix == "" + assert "PB" in inner.params From 88a87a78522d621d746198dec2b96ec9fb7537b8 Mon Sep 17 00:00:00 2001 From: Rutger van Haasteren Date: Tue, 16 Jun 2026 07:48:22 +0000 Subject: [PATCH 02/16] Use the BINARY-tagged orbit in timing-model orbital utilities Select the primary binary component by binary_param_tag so orbital phase, pulsar radial velocity, and conjunction computations follow the inner BINARY orbit when BINARY2 is also present. Add a triple-system regression test that validates these utilities match the inner-only model and differ from the outer-orbit phase. --- src/pint/models/timing_model.py | 36 +++++++++++++++++++---------- tests/test_orbit_phase.py | 40 ++++++++++++++++++++++++++++++++- 2 files changed, 63 insertions(+), 13 deletions(-) diff --git a/src/pint/models/timing_model.py b/src/pint/models/timing_model.py index 2491a3af53..40d9e1c1ac 100644 --- a/src/pint/models/timing_model.py +++ b/src/pint/models/timing_model.py @@ -953,6 +953,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], @@ -991,10 +1011,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 @@ -1062,9 +1079,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 ( @@ -1136,10 +1151,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) @@ -1161,7 +1173,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") diff --git a/tests/test_orbit_phase.py b/tests/test_orbit_phase.py index cc73fc4f32..04307916c4 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) From 7476674b25844c0f6d913f8dc27ce3e30f5bc38f Mon Sep 17 00:00:00 2001 From: Colin Clark Date: Tue, 30 Jun 2026 14:22:19 +0200 Subject: [PATCH 03/16] Account for derivatives of later delay components when computing d_delay_d_param --- src/pint/models/pulsar_binary.py | 6 ++++++ .../stand_alone_psr_binaries/ELL1_model.py | 1 + .../binary_generic.py | 14 +++++++++++++ src/pint/models/timing_model.py | 20 ++++++++++++++++--- 4 files changed, 38 insertions(+), 3 deletions(-) diff --git a/src/pint/models/pulsar_binary.py b/src/pint/models/pulsar_binary.py index 92ea36c051..1a0a25340e 100644 --- a/src/pint/models/pulsar_binary.py +++ b/src/pint/models/pulsar_binary.py @@ -269,6 +269,7 @@ 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_func += [self.d_binary_delay_d_prev_delay] def _apply_param_suffix(self): """Rename this component's PINT-facing parameters with ``param_suffix``. @@ -657,6 +658,11 @@ def d_binary_delay_d_xxxx(self, toas, param, acc_delay): 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 derivative of binary delay w.r.t. previous delays""" + 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 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 f883650190..04fa70090b 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): 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 571a98182b..6659ae02f0 100644 --- a/src/pint/models/stand_alone_psr_binaries/binary_generic.py +++ b/src/pint/models/stand_alone_psr_binaries/binary_generic.py @@ -122,6 +122,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 +263,19 @@ 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 + ------- + float + Derivative of binary delay w.r.t. earlier delays (dimensionless) + """ + + result = self.d_binarydelay_d_par(self.d_binarydelay_d_prev_delay_par).to("") + return result + def prtl_der(self, y, x): """Find the partial derivatives in binary model pdy/pdx diff --git a/src/pint/models/timing_model.py b/src/pint/models/timing_model.py index 40d9e1c1ac..e006a929c2 100644 --- a/src/pint/models/timing_model.py +++ b/src/pint/models/timing_model.py @@ -2262,10 +2262,23 @@ def d_delay_d_param( 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() + + result = np.longdouble(np.zeros(toas.ntoas) << (u.s / par.units)) + for cp in self.DelayComponent_list: + d = self.delay( + toas, cutoff_component=cp.__class__.__name__, include_last=False ) + + d_delay_d_prev_delay = 0 + for f in cp.delay_derivs_wrt_prev_delay: + d_delay_d_prev_delay += f(toas, d)().to("") + result *= 1 + d_delay_d_prev_delay + + if param in cp.deriv_funcs: + for df in cp.deriv_funcs[param]: + result += df(toas, param, d).to( + result.unit, equivalencies=u.dimensionless_angles() + ) return result def d_phase_d_param_num( @@ -4052,6 +4065,7 @@ class DelayComponent(Component): def __init__(self): super().__init__() self.delay_funcs_component = [] + self.delay_derivs_wrt_prev_delay = [] class PhaseComponent(Component): From 240ca35dded0971046a2cc0f306de26d98ea3bcd Mon Sep 17 00:00:00 2001 From: Colin Clark Date: Tue, 30 Jun 2026 14:51:41 +0200 Subject: [PATCH 04/16] Fixing previous commit --- src/pint/models/pulsar_binary.py | 2 +- src/pint/models/timing_model.py | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/pint/models/pulsar_binary.py b/src/pint/models/pulsar_binary.py index 1a0a25340e..14cb217c42 100644 --- a/src/pint/models/pulsar_binary.py +++ b/src/pint/models/pulsar_binary.py @@ -269,7 +269,7 @@ 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_func += [self.d_binary_delay_d_prev_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``. diff --git a/src/pint/models/timing_model.py b/src/pint/models/timing_model.py index e006a929c2..0bdb8e59c0 100644 --- a/src/pint/models/timing_model.py +++ b/src/pint/models/timing_model.py @@ -2270,7 +2270,7 @@ def d_delay_d_param( ) d_delay_d_prev_delay = 0 - for f in cp.delay_derivs_wrt_prev_delay: + for f in cp.delay_deriv_wrt_prev_delay_funcs: d_delay_d_prev_delay += f(toas, d)().to("") result *= 1 + d_delay_d_prev_delay @@ -4065,7 +4065,7 @@ class DelayComponent(Component): def __init__(self): super().__init__() self.delay_funcs_component = [] - self.delay_derivs_wrt_prev_delay = [] + self.delay_deriv_wrt_prev_delay_funcs = [] class PhaseComponent(Component): From d3aeae9233f0bab1ef5f513335d305a249dc9d73 Mon Sep 17 00:00:00 2001 From: Rutger van Haasteren Date: Mon, 13 Jul 2026 17:29:12 +0200 Subject: [PATCH 05/16] Add A1DOT2: second time derivative of the projected semi-major axis A hierarchical triple's outer orbit imposes a slowly varying Doppler factor on the inner orbit; over a short data span this appears as secular changes of the inner A1 beyond the linear A1DOT term. A1DOT2 parameterizes the quadratic term: a1(t) = A1 + A1DOT * (t - T0) + 0.5 * A1DOT2 * (t - T0)**2 Changes: - new A1DOT2 floatParameter (lsec/s^2) on PulsarBinary, frozen at 0 when unset (BinaryDD.validate, like the other *DOT parameters) - a1(), d_a1_d_T0, d_a1_d_A1DOT2 and d_Pobs_d_A1DOT2 in the standalone PSR_BINARY (d_a1_d_T0 now includes the A1DOT2 chain term) - d_beta_d_A1DOT2 in DDmodel for symmetry with d_beta_d_A1DOT (the generic d_beta_d_par path chains through d_a1_d_par either way) - tests: A1DOT2 affects the delay; analytic delay derivatives for A1DOT2/A1DOT/T0 and outer parameters match finite differences --- src/pint/models/binary_dd.py | 2 +- src/pint/models/pulsar_binary.py | 10 +++++ .../stand_alone_psr_binaries/DD_model.py | 12 +++++ .../binary_generic.py | 19 ++++++-- tests/test_triple_binary.py | 45 +++++++++++++++++++ 5 files changed, 83 insertions(+), 5 deletions(-) diff --git a/src/pint/models/binary_dd.py b/src/pint/models/binary_dd.py index 17420d9682..6d39c8d4af 100644 --- a/src/pint/models/binary_dd.py +++ b/src/pint/models/binary_dd.py @@ -114,7 +114,7 @@ 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"): + 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 diff --git a/src/pint/models/pulsar_binary.py b/src/pint/models/pulsar_binary.py index b04cc150e2..09f40861e5 100644 --- a/src/pint/models/pulsar_binary.py +++ b/src/pint/models/pulsar_binary.py @@ -144,6 +144,16 @@ def __init__(self): tcb2tdb_scale_factor=(1 / consts.c), ) ) + self.add_param( + floatParameter( + name="A1DOT2", + aliases=[], + 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", 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 604397dee9..f7ae12cac8 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/binary_generic.py b/src/pint/models/stand_alone_psr_binaries/binary_generic.py index b662768031..0a9e4067d9 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, @@ -408,19 +409,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." @@ -821,6 +829,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/tests/test_triple_binary.py b/tests/test_triple_binary.py index fb781111a5..60e4bf4330 100644 --- a/tests/test_triple_binary.py +++ b/tests/test_triple_binary.py @@ -177,3 +177,48 @@ def test_outer_wrapper_classes(): assert inner.category == "pulsar_system" assert inner.param_suffix == "" assert "PB" in inner.params + + +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("param", ["A1DOT2", "A1DOT", "T0", "A1_2", "T0_2"]) +def test_delay_derivatives_match_numerical(toas, param): + """Analytic delay derivatives (including A1DOT2 and the chain rule through + the outer->inner delay coupling) agree with central finite differences.""" + m = mb.get_model(TRIPLE_PAR) + # Nonzero secular terms so the T0 derivative exercises the A1DOT/A1DOT2 + # contributions to d_a1_d_T0. + m.A1DOT.quantity = 3e-13 * u.lsec / u.s + m.A1DOT2.quantity = 2e-21 * u.lsec / u.s**2 + + steps = { + "A1DOT2": 1e-22, + "A1DOT": 1e-15, + "T0": 1e-6, + "A1_2": 1e-4, + "T0_2": 1e-4, + } + 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 From f8807c608f4fa44daf221edc5dbdf261726a6fd7 Mon Sep 17 00:00:00 2001 From: Colin Clark Date: Mon, 27 Jul 2026 09:43:06 +0200 Subject: [PATCH 06/16] Added A1DOT2 parameter - second derivative of projected semi-major axis --- src/pint/models/binary_dd.py | 2 +- src/pint/models/pulsar_binary.py | 11 +++++++++++ .../stand_alone_psr_binaries/DD_model.py | 12 ++++++++++++ .../binary_generic.py | 19 +++++++++++++++---- 4 files changed, 39 insertions(+), 5 deletions(-) diff --git a/src/pint/models/binary_dd.py b/src/pint/models/binary_dd.py index 17420d9682..6d39c8d4af 100644 --- a/src/pint/models/binary_dd.py +++ b/src/pint/models/binary_dd.py @@ -114,7 +114,7 @@ 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"): + 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 diff --git a/src/pint/models/pulsar_binary.py b/src/pint/models/pulsar_binary.py index b04cc150e2..cf3374a480 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) @@ -144,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", 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 604397dee9..f7ae12cac8 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/binary_generic.py b/src/pint/models/stand_alone_psr_binaries/binary_generic.py index b662768031..0a9e4067d9 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, @@ -408,19 +409,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." @@ -821,6 +829,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) From 0ab9e06908b33583dbb7662b8148a0b4372cc093 Mon Sep 17 00:00:00 2001 From: Rutger van Haasteren Date: Mon, 27 Jul 2026 10:37:49 +0000 Subject: [PATCH 07/16] Add BinaryELL12 and fix BT/ELL1 derivatives for hierarchical triples. MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Complete BINARY2 support for BT and ELL1 so A1DOT2 and the outer→inner prev-delay chain rule match finite-difference tests for all three families. --- CHANGELOG-unreleased.md | 2 +- src/pint/models/__init__.py | 2 +- src/pint/models/binary_bt.py | 2 +- src/pint/models/binary_ell1.py | 38 ++++- .../stand_alone_psr_binaries/BT_model.py | 22 ++- .../stand_alone_psr_binaries/ELL1_model.py | 22 ++- .../binary_generic.py | 5 +- tests/test_triple_binary.py | 147 ++++++++++++++++-- 8 files changed, 209 insertions(+), 31 deletions(-) diff --git a/CHANGELOG-unreleased.md b/CHANGELOG-unreleased.md index 24da258522..dd1cd34b75 100644 --- a/CHANGELOG-unreleased.md +++ b/CHANGELOG-unreleased.md @@ -17,7 +17,7 @@ the released changes. ### 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` and `BinaryBT2` are provided. +- 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 diff --git a/src/pint/models/__init__.py b/src/pint/models/__init__.py index b48670b4da..bb16388baf 100644 --- a/src/pint/models/__init__.py +++ b/src/pint/models/__init__.py @@ -24,7 +24,7 @@ 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 34472603e7..9cc0fea5ac 100644 --- a/src/pint/models/binary_bt.py +++ b/src/pint/models/binary_bt.py @@ -70,7 +70,7 @@ def validate(self): raise MissingParameter("BT", p, f"{p} is required for BT") # If any *DOT is set, we need T0 - for p in ("PBDOT", "OMDOT", "EDOT", "A1DOT"): + 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 diff --git a/src/pint/models/binary_ell1.py b/src/pint/models/binary_ell1.py index d83f46c0b6..d61b97b29d 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/stand_alone_psr_binaries/BT_model.py b/src/pint/models/stand_alone_psr_binaries/BT_model.py index feefcceb16..17ab2a4fc0 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/ELL1_model.py b/src/pint/models/stand_alone_psr_binaries/ELL1_model.py index 04fa70090b..b803180c54 100644 --- a/src/pint/models/stand_alone_psr_binaries/ELL1_model.py +++ b/src/pint/models/stand_alone_psr_binaries/ELL1_model.py @@ -52,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 0a9e4067d9..6d2c8c750e 100644 --- a/src/pint/models/stand_alone_psr_binaries/binary_generic.py +++ b/src/pint/models/stand_alone_psr_binaries/binary_generic.py @@ -274,7 +274,10 @@ def d_binarydelay_d_prevdelay(self): Derivative of binary delay w.r.t. earlier delays (dimensionless) """ - result = self.d_binarydelay_d_par(self.d_binarydelay_d_prev_delay_par).to("") + with u.set_enabled_equivalencies(u.dimensionless_angles()): + result = self.d_binarydelay_d_par(self.d_binarydelay_d_prev_delay_par).to( + "" + ) return result def prtl_der(self, y, x): diff --git a/tests/test_triple_binary.py b/tests/test_triple_binary.py index 60e4bf4330..8fe0de66ef 100644 --- a/tests/test_triple_binary.py +++ b/tests/test_triple_binary.py @@ -14,15 +14,98 @@ 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 -from pinttestdata import datadir 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(): @@ -37,6 +120,24 @@ def _inner_only_par(): 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) @@ -164,7 +265,7 @@ def test_outer_param_derivative(triple_model, toas): def test_outer_wrapper_classes(): """The outer wrappers are configured for the BINARY2 tag and _2 suffix.""" - for cls in (BinaryDD2, BinaryBT2): + for cls in (BinaryDD2, BinaryBT2, BinaryELL12): outer = cls() assert outer.category == "pulsar_system_outer" assert outer.param_suffix == "_2" @@ -172,6 +273,10 @@ def test_outer_wrapper_classes(): 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" @@ -191,26 +296,21 @@ def test_a1dot2_changes_delay(toas): assert np.max(np.abs((d1 - d0).to_value(u.s))) > 1e-9 -@pytest.mark.parametrize("param", ["A1DOT2", "A1DOT", "T0", "A1_2", "T0_2"]) -def test_delay_derivatives_match_numerical(toas, param): +@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 = mb.get_model(TRIPLE_PAR) - # Nonzero secular terms so the T0 derivative exercises the A1DOT/A1DOT2 - # contributions to d_a1_d_T0. + 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 - steps = { - "A1DOT2": 1e-22, - "A1DOT": 1e-15, - "T0": 1e-6, - "A1_2": 1e-4, - "T0_2": 1e-4, - } ana = m.d_delay_d_param(toas, param) q = getattr(m, param) - v0, h = q.value, steps[param] + v0, h = q.value, STEPS[param] q.value = v0 + h dp = m.delay(toas) q.value = v0 - h @@ -222,3 +322,20 @@ def test_delay_derivatives_match_numerical(toas, param): 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 From 748ca967f542f73377e07bf3cdd38919059515d7 Mon Sep 17 00:00:00 2001 From: Rutger van Haasteren Date: Fri, 31 Jul 2026 07:48:28 +0000 Subject: [PATCH 08/16] fix(timing): put solar_windx before binary in DEFAULT_ORDER MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Chain-ruled d_delay_d_param requires SWX ahead of the binary; without solar_windx in DEFAULT_ORDER it sorted last and broke SW↔SWX fits. Also harden A1DOT parfile checks and Fermi observatory overwrite. --- src/pint/models/timing_model.py | 1 + tests/test_fermiphase.py | 2 +- tests/test_parfile_writing_format.py | 19 +++++++++++-------- 3 files changed, 13 insertions(+), 9 deletions(-) diff --git a/src/pint/models/timing_model.py b/src/pint/models/timing_model.py index 9dc5d5abfb..37965f6ce6 100644 --- a/src/pint/models/timing_model.py +++ b/src/pint/models/timing_model.py @@ -122,6 +122,7 @@ "troposphere", "solar_system_shapiro", "solar_wind", + "solar_windx", "dispersion_constant", "dispersion_dmx", "dispersion_jump", diff --git a/tests/test_fermiphase.py b/tests/test_fermiphase.py index 4aa62f3d33..734cc5312c 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_parfile_writing_format.py b/tests/test_parfile_writing_format.py index 222243b3d9..5f3ad589f3 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(): From 3a634bbee8b13b4b49cd3a09928255e20aef5969 Mon Sep 17 00:00:00 2001 From: Rutger van Haasteren Date: Fri, 31 Jul 2026 08:18:48 +0000 Subject: [PATCH 09/16] perf(timing): cache delay-deriv chain once per design matrix Compute parameter-independent acc_delays/chain_factors via delay_deriv_chain and thread them through designmatrix paths, restoring near-master cost while keeping the binary chain rule. Also harden triple-aware convert_binary, derived-param summary, and MissingParameter names. --- CHANGELOG-unreleased.md | 1 + src/pint/binaryconvert.py | 11 +- src/pint/models/binary_bt.py | 3 +- src/pint/models/pulsar_binary.py | 10 +- .../binary_generic.py | 10 +- src/pint/models/timing_model.py | 138 ++++++++++++++---- src/pint/pint_matrix.py | 7 +- tests/test_triple_binary.py | 28 ++++ 8 files changed, 171 insertions(+), 37 deletions(-) diff --git a/CHANGELOG-unreleased.md b/CHANGELOG-unreleased.md index c705a4396a..c2cf7f509a 100644 --- a/CHANGELOG-unreleased.md +++ b/CHANGELOG-unreleased.md @@ -15,6 +15,7 @@ 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 diff --git a/src/pint/binaryconvert.py b/src/pint/binaryconvert.py index f300ab5a32..2eef2820d3 100644 --- a/src/pint/binaryconvert.py +++ b/src/pint/binaryconvert.py @@ -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/binary_bt.py b/src/pint/models/binary_bt.py index 24a302411b..3242254e9c 100644 --- a/src/pint/models/binary_bt.py +++ b/src/pint/models/binary_bt.py @@ -67,7 +67,8 @@ def validate(self): super().validate() for p in ("T0", "A1"): if self._bp(p).value is None: - raise MissingParameter("BT", p, f"{p} is required for BT") + 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"): diff --git a/src/pint/models/pulsar_binary.py b/src/pint/models/pulsar_binary.py index b04cc150e2..52126642a3 100644 --- a/src/pint/models/pulsar_binary.py +++ b/src/pint/models/pulsar_binary.py @@ -512,9 +512,12 @@ def check_required_params(self, required_params): 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() @@ -658,9 +661,10 @@ def d_binary_delay_d_xxxx(self, toas, param, acc_delay): return self.binary_instance.d_binarydelay_d_par(param) def d_binary_delay_d_prev_delay(self, toas, acc_delay): - """Return derivative of binary delay w.r.t. previous delays""" + """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 + return self.binary_instance.d_binarydelay_d_prevdelay() def print_par(self, format="pint"): tag = self.binary_param_tag 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 b662768031..4f22d81a7f 100644 --- a/src/pint/models/stand_alone_psr_binaries/binary_generic.py +++ b/src/pint/models/stand_alone_psr_binaries/binary_generic.py @@ -269,12 +269,14 @@ def d_binarydelay_d_prevdelay(self): Returns ------- - float - Derivative of binary delay w.r.t. earlier delays (dimensionless) + 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. """ - result = self.d_binarydelay_d_par(self.d_binarydelay_d_prev_delay_par).to("") - return result + 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 diff --git a/src/pint/models/timing_model.py b/src/pint/models/timing_model.py index 37965f6ce6..7623f8e0c8 100644 --- a/src/pint/models/timing_model.py +++ b/src/pint/models/timing_model.py @@ -34,7 +34,17 @@ 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 @@ -137,6 +147,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. @@ -2197,7 +2224,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 @@ -2220,6 +2253,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 ------- @@ -2230,7 +2268,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 @@ -2246,40 +2284,89 @@ 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. " ) + if chain is None: + chain = self.delay_deriv_chain(toas) result = np.longdouble(np.zeros(toas.ntoas) << (u.s / par.units)) - for cp in self.DelayComponent_list: - d = self.delay( - toas, cutoff_component=cp.__class__.__name__, include_last=False - ) - - d_delay_d_prev_delay = 0 - for f in cp.delay_deriv_wrt_prev_delay_funcs: - d_delay_d_prev_delay += f(toas, d)().to("") - result *= 1 + d_delay_d_prev_delay - + 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, d).to( + result += df(toas, param, acc_delay).to( result.unit, equivalencies=u.dimensionless_angles() ) return result @@ -2469,7 +2556,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'] @@ -2482,7 +2570,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) @@ -3522,9 +3610,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 diff --git a/src/pint/pint_matrix.py b/src/pint/pint_matrix.py index cca7eace01..5c082189ed 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/test_triple_binary.py b/tests/test_triple_binary.py index fb781111a5..d0df30d9f6 100644 --- a/tests/test_triple_binary.py +++ b/tests/test_triple_binary.py @@ -177,3 +177,31 @@ def test_outer_wrapper_classes(): 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" \ No newline at end of file From fdc7ee5691cb5972d2798df2003c79b1f574816a Mon Sep 17 00:00:00 2001 From: Rutger van Haasteren Date: Fri, 31 Jul 2026 08:33:47 +0000 Subject: [PATCH 10/16] style(docs): black binary_generic and document BINARY2 Satisfy CI black --check, note solar_windx DEFAULT_ORDER in the changelog, and mention BINARY2 in timing-model docs. --- CHANGELOG-unreleased.md | 1 + docs/explanation.rst | 3 +++ docs/timingmodels.rst | 6 ++++-- src/pint/models/stand_alone_psr_binaries/binary_generic.py | 4 +--- 4 files changed, 9 insertions(+), 5 deletions(-) diff --git a/CHANGELOG-unreleased.md b/CHANGELOG-unreleased.md index f0b284cd35..f8724f8b7e 100644 --- a/CHANGELOG-unreleased.md +++ b/CHANGELOG-unreleased.md @@ -38,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 a6bd78c769..3d3d24c8cd 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 252fd4d658..b714cf7f1f 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/models/stand_alone_psr_binaries/binary_generic.py b/src/pint/models/stand_alone_psr_binaries/binary_generic.py index 70589763b6..a7fa071a00 100644 --- a/src/pint/models/stand_alone_psr_binaries/binary_generic.py +++ b/src/pint/models/stand_alone_psr_binaries/binary_generic.py @@ -278,9 +278,7 @@ def d_binarydelay_d_prevdelay(self): """ with u.set_enabled_equivalencies(u.dimensionless_angles()): - return self.d_binarydelay_d_par(self.d_binarydelay_d_prev_delay_par).to( - "" - ) + 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 From 1175030989c936d23185f4c00b790e9c527e57b1 Mon Sep 17 00:00:00 2001 From: Rutger van Haasteren Date: Fri, 31 Jul 2026 08:35:52 +0000 Subject: [PATCH 11/16] chore: drop unused imports in timing_model and binaryconvert Keep MissingParameter/MissingBinaryError/MissingTOAs re-exported from timing_model (listed in __all__) for existing importers. Remove a dead unreachable return in binary_generic. --- src/pint/binaryconvert.py | 2 +- .../models/stand_alone_psr_binaries/binary_generic.py | 1 - src/pint/models/timing_model.py | 8 +++++--- 3 files changed, 6 insertions(+), 5 deletions(-) diff --git a/src/pint/binaryconvert.py b/src/pint/binaryconvert.py index 2eef2820d3..101bd4f4a9 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 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 a7fa071a00..e2adb9bec8 100644 --- a/src/pint/models/stand_alone_psr_binaries/binary_generic.py +++ b/src/pint/models/stand_alone_psr_binaries/binary_generic.py @@ -556,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)""" diff --git a/src/pint/models/timing_model.py b/src/pint/models/timing_model.py index 7623f8e0c8..a377488822 100644 --- a/src/pint/models/timing_model.py +++ b/src/pint/models/timing_model.py @@ -50,8 +50,7 @@ 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 @@ -91,7 +90,6 @@ colorize, open_or_use, split_prefixed_name, - xxxselections, get_unit, ) @@ -100,6 +98,10 @@ "TimingModel", "Component", "AllComponents", + # Re-exported for compatibility with older imports. + "MissingBinaryError", + "MissingParameter", + "MissingTOAs", ] # Parameters or lines in par files we don't understand but shouldn't # complain about. These are still passed to components so that they From 672906443141b9cf578b52b326e7e933f0b8cbc9 Mon Sep 17 00:00:00 2001 From: Rutger van Haasteren Date: Fri, 31 Jul 2026 08:41:43 +0000 Subject: [PATCH 12/16] refactor: import Missing* exceptions from pint.exceptions Stop re-exporting them from timing_model; point remaining call sites at the canonical module. --- src/pint/models/chromatic_model.py | 4 ++-- src/pint/models/solar_wind_dispersion.py | 2 +- src/pint/models/timing_model.py | 6 ------ tests/test_binary_generic.py | 2 +- tests/test_chromatic.py | 2 +- tests/test_ddk.py | 2 +- tests/test_fitter_error_checking.py | 2 +- tests/test_piecewise.py | 2 +- tests/test_wavex.py | 2 +- 9 files changed, 9 insertions(+), 15 deletions(-) diff --git a/src/pint/models/chromatic_model.py b/src/pint/models/chromatic_model.py index e1f3a327c0..e6335e53f8 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/solar_wind_dispersion.py b/src/pint/models/solar_wind_dispersion.py index ff81471172..112b33dc9e 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/timing_model.py b/src/pint/models/timing_model.py index a377488822..c65ebfc407 100644 --- a/src/pint/models/timing_model.py +++ b/src/pint/models/timing_model.py @@ -60,8 +60,6 @@ from pint.derived_quantities import dispersion_slope from pint.exceptions import ( AliasConflict, - MissingBinaryError, - MissingParameter, MissingTOAs, PrefixError, PropertyAttributeError, @@ -98,10 +96,6 @@ "TimingModel", "Component", "AllComponents", - # Re-exported for compatibility with older imports. - "MissingBinaryError", - "MissingParameter", - "MissingTOAs", ] # Parameters or lines in par files we don't understand but shouldn't # complain about. These are still passed to components so that they diff --git a/tests/test_binary_generic.py b/tests/test_binary_generic.py index 4caac397db..c32ba91d7a 100644 --- a/tests/test_binary_generic.py +++ b/tests/test_binary_generic.py @@ -5,7 +5,7 @@ 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 diff --git a/tests/test_chromatic.py b/tests/test_chromatic.py index 227553ab74..b0808484b6 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 e9401cc6d8..d32171ebe8 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_fitter_error_checking.py b/tests/test_fitter_error_checking.py index bcdc7f3acb..b352c0247f 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_piecewise.py b/tests/test_piecewise.py index 0809a63289..955488a412 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_wavex.py b/tests/test_wavex.py index 7c5c7f9f59..2a60decc13 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 From 22636053e8e78610dc8b4df6d50a599b997df487 Mon Sep 17 00:00:00 2001 From: Rutger van Haasteren Date: Fri, 31 Jul 2026 09:14:58 +0000 Subject: [PATCH 13/16] fix(tests): import MissingParameter from pint.exceptions Stops collection failure in test_model_manual after Missing* stopped being re-exported from timing_model. --- tests/test_model_manual.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/tests/test_model_manual.py b/tests/test_model_manual.py index de406b7b4d..8ffc5fb60f 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 From 7aa3f8ff390fc4b30865bdb266c7ed42e5a14305 Mon Sep 17 00:00:00 2001 From: Rutger van Haasteren Date: Fri, 11 Sep 2026 09:52:20 +0200 Subject: [PATCH 14/16] Updated explanation.rst regarding hierarchical triple system --- docs/explanation.rst | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/docs/explanation.rst b/docs/explanation.rst index cdae31d125..d8adcafeda 100644 --- a/docs/explanation.rst +++ b/docs/explanation.rst @@ -597,6 +597,12 @@ components. 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``). + The outer delay is accumulated first, and the inner binary is evaluated at + barycentric time minus that accumulated delay, so the outer Rømer delay + Doppler-shifts the inner orbit. This is the usual hierarchical + light-travel-time coupling of two independent Keplerian orbits; it is not + a dynamical three-body model (there is no companion–companion gravity, + Kozai–Lidov coupling, or mutual Shapiro delay between the companions). - 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 From 4b483dee9ced378410342368b0abe9dd60b03f2d Mon Sep 17 00:00:00 2001 From: Rutger van Haasteren Date: Fri, 11 Sep 2026 10:52:02 +0200 Subject: [PATCH 15/16] test(binary): FB0/FB1 inner orbit in hierarchical triple Check that the inner binary of a hierarchical triple can use the orbital-frequency (FBn) parameterization: FB0/FB1 live on the inner BinaryDD, the outer BinaryDD2 carries no FBn parameters, the delay matches the equivalent PB/PBDOT triple to numerical noise, and the FBn derivatives stay finite and nonzero. --- tests/test_triple_binary.py | 56 +++++++++++++++++++++++++++++++++++++ 1 file changed, 56 insertions(+) diff --git a/tests/test_triple_binary.py b/tests/test_triple_binary.py index 6a31f772bf..8c8c3ee48e 100644 --- a/tests/test_triple_binary.py +++ b/tests/test_triple_binary.py @@ -367,3 +367,59 @@ def test_outer_orbit_affects_delay_family(family): m.A1_2.quantity = 0 * u.lsec d1 = m.delay(toas) assert np.max(np.abs((d1 - d0).to_value(u.s))) > 1e-9 + + +def _fb_triple_par(fb1=None): + """Return the DD triple parfile text with the inner orbit expressed in + orbital frequency (``FB0`` = 1/``PB``, optionally ``FB1``) instead of ``PB``.""" + lines = [] + with open(TRIPLE_PAR) as f: + for line in f: + parts = line.split() + if parts and parts[0] == "PB": + pb_s = np.longdouble(parts[1]) * 86400 + lines.append(f"FB0 {1 / pb_s:.25e} 1\n") + if fb1 is not None: + lines.append(f"FB1 {fb1:.25e} 1\n") + continue + lines.append(line) + return "".join(lines) + + +@pytest.mark.parametrize("fb1", [None, -3.0e-20]) +def test_inner_orbit_fb_parameterization(triple_model, toas, fb1): + """The inner orbit of a hierarchical triple can use the orbital-frequency + (``FBn``) parameterization: ``FB0``/``FB1`` are owned by the inner + ``BinaryDD``, the outer ``BinaryDD2`` has no ``FBn`` parameters, and the + delay matches the equivalent ``PB``/``PBDOT`` triple to numerical noise.""" + m_fb = mb.get_model(io.StringIO(_fb_triple_par(fb1))) + + inner = m_fb.components["BinaryDD"] + outer = m_fb.components["BinaryDD2"] + assert "FB0" in inner.params + assert inner.FB0.value is not None + assert m_fb.PB.value is None + assert not any(p.startswith("FB") for p in outer.params) + assert "PB_2" in outer.params + assert m_fb.PB_2.quantity == triple_model.PB_2.quantity + + d_fb = m_fb.delay(toas) + assert np.all(np.isfinite(d_fb.value)) + + # PB-equivalent reference: PBDOT = -FB1 / FB0**2. + m_pb = mb.get_model(TRIPLE_PAR) + if fb1 is not None: + assert "FB1" in inner.params + m_pb.PBDOT.quantity = (-m_fb.FB1.quantity / m_fb.FB0.quantity**2).to( + u.dimensionless_unscaled + ) + d_pb = m_pb.delay(toas) + + diff = np.max(np.abs((d_fb - d_pb).to_value(u.s))) + assert diff < 1e-11 + + # The FBn parameters of the inner orbit remain fittable inside the triple. + for p in ["FB0"] + (["FB1"] if fb1 is not None else []): + d = m_fb.d_delay_d_param(toas, p) + assert np.all(np.isfinite(d.value)) + assert np.any(d.value != 0) From 7e39259596d8f099555c9d6b6fcccd0bf4cdaa13 Mon Sep 17 00:00:00 2001 From: Rutger van Haasteren Date: Fri, 11 Sep 2026 11:31:11 +0200 Subject: [PATCH 16/16] Added note in explanation.rst that the outer orbit is period-parameterized only. No FBn and ORBWAVE are supported, those are only valid parameters on the inner orbit --- docs/explanation.rst | 2 ++ 1 file changed, 2 insertions(+) diff --git a/docs/explanation.rst b/docs/explanation.rst index d8adcafeda..fef8073681 100644 --- a/docs/explanation.rst +++ b/docs/explanation.rst @@ -603,6 +603,8 @@ components. light-travel-time coupling of two independent Keplerian orbits; it is not a dynamical three-body model (there is no companion–companion gravity, Kozai–Lidov coupling, or mutual Shapiro delay between the companions). + The outer orbit is period-parameterized only: ``FBn`` and ``ORBWAVE`` are + not supported for ``BINARY2`` (the inner orbit may still use them). - 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