diff --git a/CHANGELOG-unreleased.md b/CHANGELOG-unreleased.md index c199d6083..24da25852 100644 --- a/CHANGELOG-unreleased.md +++ b/CHANGELOG-unreleased.md @@ -17,6 +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. ### 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 d13ceec8c..b48670b4d 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 4eb81fa29..34472603e 100644 --- a/src/pint/models/binary_bt.py +++ b/src/pint/models/binary_bt.py @@ -66,18 +66,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 39352bbdb..17420d968 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 84a8f4c13..9216c415e 100644 --- a/src/pint/models/model_builder.py +++ b/src/pint/models/model_builder.py @@ -271,8 +271,8 @@ def _validate_components(self): f" class, please set register to 'False' in the class" f" of component {k}." ) - if v.category == "pulsar_system": - # The pulsar system will be selected by parameter BINARY + if v.category in ("pulsar_system", "pulsar_system_outer"): + # The pulsar system is selected by parameter BINARY/BINARY2 continue else: raise ComponentConflict(m) @@ -488,6 +488,12 @@ def choose_model(self, param_inpar, force_binary_model=None, allow_T2=False): self.choose_binary_model(param_inpar, force_binary_model, allow_T2) ) + # Outer-orbit binary for hierarchical triple systems. + binary2 = param_inpar.get("BINARY2", None) + if binary2: + binary2 = binary2[0] + selected_components.add(self.choose_outer_binary_model(param_inpar)) + # 2. Get the component list from the parameters in the parfile. # 2.1 Check the aliases of input parameters. # This does not include the repeating parameters, but it should not @@ -531,6 +537,17 @@ def choose_model(self, param_inpar, force_binary_model=None, allow_T2=False): ) else: continue + # Outer-orbit binary component, controlled by the BINARY2 tag. + if self.all_components.components[cps[0]].category == "pulsar_system_outer": + if binary2 is None: + raise MissingBinaryError( + f"The outer-orbit binary model is decided by the" + f" parameter 'BINARY2'. Please indicate the outer" + f" binary model before using parameter {k}, which" + f" is an outer binary model parameter." + ) + else: + continue if len(cps) == 1: # No conflict, parameter only shows in one component. selected_components.add(cps[0]) @@ -645,6 +662,32 @@ def choose_binary_model(self, param_inpar, force_binary_model=None, allow_T2=Fal return binary_cp.__class__.__name__ + def choose_outer_binary_model(self, param_inpar): + """Choose the outer-orbit BINARY2 model for a hierarchical triple. + + Parameters + ---------- + param_inpar: dict + Dictionary of the unique parameters in .par file with the key being + the parfile line. :func:`parse_parfile` returns this dictionary. + + Returns + ------- + str + Name of the outer binary component class. + + Note + ---- + The outer model is taken verbatim from the ``BINARY2`` parameter (no T2 + guessing is performed) and must correspond to a registered component in + the ``pulsar_system_outer`` category (e.g. ``DD`` -> ``BinaryDD2``). + """ + binary2 = param_inpar["BINARY2"][0] + binary_cp = self.all_components.search_binary_components( + binary2, category="pulsar_system_outer" + ) + return binary_cp.__class__.__name__ + def _setup_model( self, timing_model, diff --git a/src/pint/models/pulsar_binary.py b/src/pint/models/pulsar_binary.py index d9b50460b..b04cc150e 100644 --- a/src/pint/models/pulsar_binary.py +++ b/src/pint/models/pulsar_binary.py @@ -87,6 +87,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 @@ -258,6 +268,69 @@ def __init__(self): self.warn_default_params = ["ECC", "OM"] # Set up delay function self.delay_funcs_component += [self.binarymodel_delay] + self.delay_deriv_wrt_prev_delay_funcs += [self.d_binary_delay_d_prev_delay] + + def _apply_param_suffix(self): + """Rename this component's PINT-facing parameters with ``param_suffix``. + + This is used by outer-orbit components (e.g. the outer binary of a + hierarchical triple) so that their parameters appear in the parfile as + ``PB_2``, ``A1_2``, ... while the underlying standalone binary instance + keeps the canonical names (``PB``, ``A1``, ...). + + The canonical attribute (e.g. ``self.PB``) is *removed* from this + component so that it does not leak into the parent + :class:`~pint.models.timing_model.TimingModel` namespace (where it would + collide with the inner binary's identically named parameter). Base-class + methods therefore access parameters through :meth:`_bp` / :meth:`_hasbp`, + which apply the suffix. + + Prefix parameters (``FBn``, ``ORBWAVECn``, ...) are removed entirely: + orbital-frequency and ORBWAVE parameterizations are not supported for a + suffixed (outer) orbit. + """ + suffix = self.param_suffix + if not suffix: + return + for canonical in list(self.params): + par = getattr(self, canonical) + if getattr(par, "is_prefix", False): + # Orbital-frequency / ORBWAVE parameterizations are not + # supported for a suffixed (outer) orbit. + self.remove_param(canonical) + continue + new_name = canonical + suffix + par.aliases = [alias + suffix for alias in par.aliases] + par.name = new_name + # Keep any funcParameter cross-references consistent with the + # renamed parameters (no-op for DD/BT, which have none). + if hasattr(par, "_params"): + par._params = [p + suffix for p in par._params] + # Expose the renamed parameter under its suffixed name only and drop + # the canonical attribute so it cannot leak to the parent model. + setattr(self, new_name, par) + delattr(self, canonical) + self.params[self.params.index(canonical)] = new_name + + def _bp(self, name): + """Return this component's parameter object for canonical ``name``. + + Applies :attr:`param_suffix` so that base-class code written in terms of + canonical names (e.g. ``PB``) resolves to the suffixed parameter (e.g. + ``PB_2``) for an outer-orbit component. For a normal binary + (``param_suffix == ""``) this is just ``getattr(self, name)``. + """ + suffixed = name + self.param_suffix + if self.param_suffix and hasattr(self, suffixed): + return getattr(self, suffixed) + return getattr(self, name) + + def _hasbp(self, name): + """Whether this component has the (possibly suffixed) parameter ``name``.""" + suffixed = name + self.param_suffix + return bool(self.param_suffix and hasattr(self, suffixed)) or hasattr( + self, name + ) def setup(self): super().setup() @@ -378,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" @@ -508,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) @@ -555,19 +651,30 @@ def binarymodel_delay(self, toas, acc_delay=None): def d_binary_delay_d_xxxx(self, toas, param, acc_delay): """Return the binary model delay derivatives.""" self.update_binary_object(toas, acc_delay) + # The standalone binary instance only knows the canonical (unsuffixed) + # parameter names, so strip the suffix before requesting a derivative. + if self.param_suffix and param.endswith(self.param_suffix): + param = param[: -len(self.param_suffix)] return self.binary_instance.d_binarydelay_d_par(param) + def d_binary_delay_d_prev_delay(self, toas, acc_delay): + """Return 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 if self._parent is None: - result = f"BINARY {self.binary_model_name}\n" - elif self._parent.BINARY.value != self.binary_model_name: + result = f"{tag} {self.binary_model_name}\n" + elif tag_par.value != self.binary_model_name: raise TimingModelError( - f"Parameter BINARY {self._parent.BINARY.value}" + f"Parameter {tag} {tag_par.value}" f" does not match the binary" f" component {self.binary_model_name}" ) else: - result = self._parent.BINARY.as_parfile_line(format=format) + result = tag_par.as_parfile_line(format=format) for p in self.params: par = getattr(self, p) @@ -617,32 +724,38 @@ def change_binary_epoch(self, new_epoch): """ new_epoch = parse_time(new_epoch, scale="tdb", precision=9) + # Parameter access goes through _bp() so that this works for both a + # normal (inner) binary and a suffixed outer-orbit component. + PB_par = self._bp("PB") + PBDOT_par = self._bp("PBDOT") + T0_par = self._bp("T0") + # Get PB and PBDOT from model - if self.PB.quantity is not None and not isinstance(self.PB, funcParameter): - PB = self.PB.quantity - if self.PBDOT.quantity is not None: - PBDOT = self.PBDOT.quantity + if PB_par.quantity is not None and not isinstance(PB_par, funcParameter): + PB = PB_par.quantity + if PBDOT_par.quantity is not None: + PBDOT = PBDOT_par.quantity else: PBDOT = 0.0 * u.Unit("") else: - PB = 1.0 / self.FB0.quantity + PB = 1.0 / self._bp("FB0").quantity try: - PBDOT = -self.FB1.quantity / self.FB0.quantity**2 + PBDOT = -self._bp("FB1").quantity / self._bp("FB0").quantity ** 2 except AttributeError: PBDOT = 0.0 * u.Unit("") # Find closest periapsis time and reassign T0 - t0_ld = self.T0.quantity.tdb.mjd_long + t0_ld = T0_par.quantity.tdb.mjd_long dt = (new_epoch.tdb.mjd_long - t0_ld) * u.day d_orbits = dt / PB - PBDOT * dt**2 / (2.0 * PB**2) n_orbits = np.round(d_orbits.to(u.Unit(""))) if n_orbits == 0: return dt_integer_orbits = PB * n_orbits + PB * PBDOT * n_orbits**2 / 2.0 - self.T0.quantity = self.T0.quantity + dt_integer_orbits + T0_par.quantity = T0_par.quantity + dt_integer_orbits with contextlib.suppress(AttributeError): - if self.FB2.quantity is not None: + if self._bp("FB2").quantity is not None: log.warning( "Ignoring orbital frequency derivatives higher than FB1" "in computing new T0; a model fit should resolve this" @@ -650,23 +763,23 @@ def change_binary_epoch(self, new_epoch): # Update PB or FB0, FB1, etc. if isinstance(self.binary_instance.orbits_cls, bo.OrbitPB): dPB = PBDOT * dt_integer_orbits - self.PB.quantity = self.PB.quantity + dPB + PB_par.quantity = PB_par.quantity + dPB else: fbterms = [0.0 * u.Unit("")] + self._parent.get_prefix_list("FB") for n in range(len(fbterms) - 1): - cur_deriv = getattr(self, f"FB{n}") + cur_deriv = self._bp(f"FB{n}") cur_deriv.value = taylor_horner_deriv( dt_integer_orbits.to(u.s), fbterms, deriv_order=n + 1 ) # Update ECC, OM, and A1 - dECC = self.EDOT.quantity * dt_integer_orbits - self.ECC.quantity = self.ECC.quantity + dECC - dOM = self.OMDOT.quantity * dt_integer_orbits - self.OM.quantity = self.OM.quantity + dOM - dA1 = self.A1DOT.quantity * dt_integer_orbits - self.A1.quantity = self.A1.quantity + dA1 + dECC = self._bp("EDOT").quantity * dt_integer_orbits + self._bp("ECC").quantity = self._bp("ECC").quantity + dECC + dOM = self._bp("OMDOT").quantity * dt_integer_orbits + self._bp("OM").quantity = self._bp("OM").quantity + dOM + dA1 = self._bp("A1DOT").quantity * dt_integer_orbits + self._bp("A1").quantity = self._bp("A1").quantity + dA1 def pb(self, t=None): """Return binary period and uncertainty (optionally evaluated at different times) regardless of binary model @@ -684,22 +797,23 @@ def pb(self, t=None): Binary period uncertainty """ + PB_par = self._bp("PB") + PBDOT_par = self._bp("PBDOT") if self.binary_model_name.startswith("ELL1"): - t0 = self.TASC.quantity + t0 = self._bp("TASC").quantity else: - t0 = self.T0.quantity + t0 = self._bp("T0").quantity t = t0 if t is None else parse_time(t) - if self.PB.quantity is not None: - if self.PBDOT.quantity is None and ( - not hasattr(self, "XPBDOT") - or getattr(self, "XPBDOT").quantity is not None + if PB_par.quantity is not None: + if PBDOT_par.quantity is None and ( + not self._hasbp("XPBDOT") or self._bp("XPBDOT").quantity is not None ): - return self.PB.quantity, self.PB.uncertainty - pb = self.PB.as_ufloat(u.d) - if self.PBDOT.quantity is not None: - pbdot = self.PBDOT.as_ufloat(u.s / u.s) - if hasattr(self, "XPBDOT") and self.XPBDOT.quantity is not None: - pbdot += self.XPBDOT.as_ufloat(u.s / u.s) + return PB_par.quantity, PB_par.uncertainty + pb = PB_par.as_ufloat(u.d) + if PBDOT_par.quantity is not None: + pbdot = PBDOT_par.as_ufloat(u.s / u.s) + if self._hasbp("XPBDOT") and self._bp("XPBDOT").quantity is not None: + pbdot += self._bp("XPBDOT").as_ufloat(u.s / u.s) pnew = pb + pbdot * (t - t0).jd if not isinstance(pnew, np.ndarray): return pnew.n * u.d, pnew.s * u.d if pnew.s > 0 else None @@ -710,7 +824,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/stand_alone_psr_binaries/ELL1_model.py b/src/pint/models/stand_alone_psr_binaries/ELL1_model.py index f88365019..04fa70090 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 d6d61b298..b66276803 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 73954157f..c77e94427 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_components_of_type(PulsarBinary) <= 1 - ), "Model can have at most one PulsarBinary component." + num_binaries_with_tag("BINARY") <= 1 + ), "Model can have at most one inner PulsarBinary component." + assert ( + num_binaries_with_tag("BINARY2") <= 1 + ), "Model can have at most one outer PulsarBinary component." from pint.models.solar_wind_dispersion import ( SolarWindDispersion, @@ -925,6 +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], @@ -963,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 @@ -1034,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 ( @@ -1108,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) @@ -1133,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") @@ -2222,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_deriv_wrt_prev_delay_funcs: + 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( @@ -3133,7 +3186,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: @@ -4011,6 +4065,7 @@ class DelayComponent(Component): def __init__(self): super().__init__() self.delay_funcs_component = [] + self.delay_deriv_wrt_prev_delay_funcs = [] class PhaseComponent(Component): @@ -4223,13 +4278,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 +4303,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 000000000..346917fc1 --- /dev/null +++ b/tests/datafile/B1855+09_triple_DD.par @@ -0,0 +1,36 @@ +PSRJ 1855+09 +RAJ 18:57:36.3932884 1 0.00002602730280675029 +DECJ +09:43:17.29196 1 0.00078789485676919773 +F0 186.49408156698235146 1 0.00000000000698911818 +F1 -6.2049547277487420583e-16 1 1.7380934373573401505e-20 +PEPOCH 49453 +POSEPOCH 49453 +DMEPOCH 49453 +DM 13.29709 +PMRA -2.5054345161030380639 1 0.03104958261053317181 +PMDEC -5.4974558631993817232 1 0.06348008663748286318 +PX 1.2288569063263405232 1 0.21243361289239687251 +SINI 0.99741717335200923866 1 0.00182023515130851988 +BINARY DD +PB 12.327171194774200418 1 0.00000000079493185824 +T0 49452.940695077335647 1 0.00169031830532837251 +A1 9.2307804312998001928 1 0.00000036890718667634 +OM 276.55142180589701234 1 0.04936551005019605698 +ECC 2.1745265668236919017e-05 1 0.00000004027191312623 +M2 0.26111312480723428917 1 0.02616161008932908066 +BINARY2 DD +PB_2 1400.0 1 +T0_2 49452.0 1 +A1_2 120.0 1 +OM_2 110.0 1 +ECC_2 0.3 1 +START 53358.726464889485214 +FINISH 55108.922917417192366 +TZRMJD 54177.508359343262555 +TZRFRQ 424 +TZRSITE ao +TRES 0.395 +CLK TT(TAI) +MODE 1 +UNITS TDB +EPHEM DE405 diff --git a/tests/test_binary_generic.py b/tests/test_binary_generic.py index dd7d5e352..4caac397d 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_orbit_phase.py b/tests/test_orbit_phase.py index cc73fc4f3..04307916c 100644 --- a/tests/test_orbit_phase.py +++ b/tests/test_orbit_phase.py @@ -1,5 +1,5 @@ -import pytest import os +import io import pytest import numpy as np @@ -91,3 +91,41 @@ def test_j0737(self): assert len(x) == 2, "conjunction is not returning an array" # make sure true anomaly before T0 is positive assert mJ0737.orbital_phase(52000.0, anom="true").value > 0.0 + + def test_triple_orbital_utilities_use_inner_binary(self): + with open(os.path.join(datadir, "B1855+09_triple_DD.par")) as f: + triple_lines = f.readlines() + triple_model = m.get_model(io.StringIO("".join(triple_lines))) + + inner_only_par = "".join( + line + for line in triple_lines + if line.split() + and line.split()[0] != "BINARY2" + and not line.split()[0].endswith("_2") + ) + inner_model = m.get_model(io.StringIO(inner_only_par)) + + # The timing-model orbital utilities should follow the BINARY-tagged + # (inner) component, not the outer BINARY2 component. + ts = triple_model.T0.value + np.linspace(0, triple_model.PB.value, 32) + triple_phase = triple_model.orbital_phase(ts, anom="mean", radians=False) + inner_phase = inner_model.orbital_phase(ts, anom="mean", radians=False) + assert np.allclose(triple_phase, inner_phase) + assert np.allclose( + triple_model.pulsar_radial_velocity(ts), + inner_model.pulsar_radial_velocity(ts), + ) + assert np.isclose( + triple_model.conjunction(triple_model.T0.value), + inner_model.conjunction(inner_model.T0.value), + ) + + # Sanity check that the outer model's mean anomaly is generally different. + outer = triple_model.components["BinaryDD2"] + outer.update_binary_object(None) + outer.binary_instance.update_input(barycentric_toa=np.asarray(ts)) + outer_phase = np.remainder(outer.binary_instance.M().value, 2 * np.pi) / ( + 2 * np.pi + ) + assert not np.allclose(triple_phase, outer_phase) diff --git a/tests/test_triple_binary.py b/tests/test_triple_binary.py new file mode 100644 index 000000000..fb781111a --- /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