Merge branch 'nanograv:master' into feat/jgb-solar-wind-gp

.github/workflows/ci_test.yml

@@ -24,16 +24,16 @@ jobs:
       matrix:
         include:
           - os: ubuntu-latest
-            python: '3.12'
+            python: '3.13'
             tox_env: 'ephemeris_connection'
           - os: ubuntu-latest
-            python: '3.12'
+            python: '3.13'
             tox_env: 'black'
           - os: ubuntu-latest
-            python: '3.12'
-            tox_env: 'py312-test-cov'
+            python: '3.13'
+            tox_env: 'py313-test-cov'
           - os: ubuntu-latest
-            python: '3.12'
+            python: '3.13'
             tox_env: 'notebooks'
 #          - os: ubuntu-latest
 #            python: '3.8'
@@ -42,13 +42,13 @@ jobs:
 #            python: '3.10'
 #            tox_env: 'py310-test-alldeps-cov'
           - os: macos-12
-            python: '3.12'
-            tox_env: 'py312-test'
+            python: '3.13'
+            tox_env: 'py313-test'
 #          - os: windows-latest
 #            python: '3.8'
 #            tox_env: 'py38-test'
           - os: ubuntu-latest
-            python: '3.8'
+            python: '3.9'
             tox_env: 'oldestdeps'
 #          - os: ubuntu-latest
 #            python: '3.8'

.pre-commit-config.yaml

@@ -6,6 +6,6 @@ repos:
     -   id: check-merge-conflict
     -   id: check-symlinks
 -   repo: https://github.com/psf/black
-    rev: 24.2.0
+    rev: 24.10.0
     hooks:
     -   id: black

CHANGELOG-unreleased.md

@@ -9,6 +9,16 @@ the released changes.
 
 ## Unreleased
 ### Changed
+- Command line scripts now automatically do `allow_tcb` and `allow_T2` while reading par files.
+- Updated the `plot_chains` function in `event_optimize` so that the subplots are a fixed size to prevent the subplots from being condensed in the case of many fit parameters.
 ### Added
+- Added an option `nbin` to `photonphase` to decide how many phase bins to use for the phaseogram
+- Added an option `linearize_model` to speed up the photon phases calculation within `event_optimize` through the designmatrix.
+- Added AIC and BIC calculation to be written in the post fit parfile from `event_optimize`
+- When TCB->TDB conversion info is missing, will print parameter name
+- `add_param` returns the name of the parameter (useful for numbered parameters)
 ### Fixed
+- Changed WAVE_OM units from 1/d to rad/d.
+- When EQUAD is created from TNEQ, has proper TCB->TDB conversion info
+- TOA selection masks will work when only TOA is the first one
 ### Removed

CHANGELOG.md

@@ -7,6 +7,12 @@ and this project, at least loosely, adheres to [Semantic Versioning](https://sem
 This file contains the released changes to the codebase. See CHANGELOG-unreleased.md for
 the unreleased changes. This file should only be changed while tagging a new version.
 
+## [1.1] 2024-11-05
+### Changed
+* Bump oldest python to 3.9
+* Change oldest dependencies: `numpy` 1.18.5 to 1.23.0; `astropy` 4.0 to 5.0.5; `scipy` 1.4.1 to 1.9.0; `matplotlib` 3.2.0 to 3.4.3
+* Update CI testing to use python 3.13
+
 ## [1.0.2] 2024-10-18
 ### Changed
 - Moved the events -> TOAs and photon weights code into the function `load_events_weights` within `event_optimize`.

docs/installation.rst

@@ -14,7 +14,7 @@ is more complicated (but not too much).
 Prerequisites
 -------------
 
-PINT requires Python 3.8+ [1]_
+PINT requires Python 3.9+ [1]_
 
 Your Python must have the package installation tool pip_ installed.  Also make sure your ``setuptools`` are up to date (e.g. ``pip install -U setuptools``).
 

pyproject.toml

@@ -19,22 +19,22 @@ classifiers = [
     "Operating System :: OS Independent",
     "Programming Language :: Python",
     "Programming Language :: Python :: 3",
-    "Programming Language :: Python :: 3.8",
     "Programming Language :: Python :: 3.9",
     "Programming Language :: Python :: 3.10",
     "Programming Language :: Python :: 3.11",
     "Programming Language :: Python :: 3.12",
+    "Programming Language :: Python :: 3.13",
     "Topic :: Scientific/Engineering :: Astronomy",
     "Topic :: Software Development :: Libraries :: Python Modules",
 ]
-requires-python = ">=3.8"
+requires-python = ">=3.9"
 dependencies = [
-    "numpy>=1.18.5",
-    "astropy>=4.0,!=4.0.1,!=4.0.1.post1",
+    "numpy>=1.23.0",
+    "astropy>=5.0.5",
     "pyerfa",
-    "scipy>=1.4.1",
+    "scipy>=1.9.0",
     "jplephem>=2.6",
-    "matplotlib>=3.2.0",
+    "matplotlib>=3.4.3",
     "emcee>=3.0.1",
     "corner>=2.0.1",
     "uncertainties",

requirements.txt

@@ -1,9 +1,9 @@
-numpy>=1.18.5
-astropy>=4.0,!=4.0.1,!=4.0.1.post1
+numpy>=1.23.0
+astropy>=5.0.5
 pyerfa
-scipy>=1.4.1
+scipy>=1.9.0
 jplephem>=2.6
-matplotlib>=3.2.0
+matplotlib>=3.4.3
 emcee>=3.0.1
 corner>=2.0.1
 uncertainties

requirements_dev.txt

@@ -1,10 +1,9 @@
 attrs>=19.2
 pip>=9.0.1
-setuptools>=41.0
+setuptools>=61.2
 coverage>=4.3.4
 traitlets
 Sphinx==6.2.1
-astropy>=4.0,!=4.0.1,!=4.0.1.post1
 #astropy-helpers>=1.3
 sphinx-rtd-theme==1.2.2
 coveralls>=1.1
@@ -29,8 +28,7 @@ pdbpp
 tox
 pre-commit
 typed-ast>=1.5.0
-#black>=19.0a0,<20.0a0
-black~=23.0
+black~=24.0
 pygments
 ipython
 pathlib2

src/pint/__init__.py

@@ -16,7 +16,6 @@
 import astropy.time as time
 import astropy.units as u
 import numpy as np
-import pkg_resources
 from astropy.units import si
 from pathlib import Path
 
@@ -106,6 +105,9 @@
     "hourangle_second": hourangle_second,
 }
 
+# define a units equivalency for gauss in cgs
+gauss_equiv = [u.Gauss, u.Hz * (u.g / u.cm) ** (1 / 2), lambda x: x, lambda x: x]
+
 import astropy.version
 
 if astropy.version.major < 4:

src/pint/config.py

@@ -1,7 +1,7 @@
 """Functions related to PINT configuration."""
 
 import os
-import pkg_resources
+import importlib.resources
 
 __all__ = ["datadir", "examplefile", "runtimefile"]
 
@@ -15,7 +15,7 @@ def datadir() -> str:
     str
         Directory of PINT data files
     """
-    return pkg_resources.resource_filename(__name__, "data/")
+    return os.path.join(importlib.resources.files("pint"), "data/")
 
 
 def examplefile(filename: str) -> str:
@@ -35,9 +35,7 @@ def examplefile(filename: str) -> str:
     This is **not** for files needed at runtime. Those are located by :func:`pint.config.runtimefile`.
     This is for files needed for the example notebooks.
     """
-    return pkg_resources.resource_filename(
-        __name__, os.path.join("data/examples/", filename)
-    )
+    return os.path.join(importlib.resources.files("pint"), f"data/examples/{filename}")
 
 
 def runtimefile(filename: str) -> str:
@@ -57,6 +55,4 @@ def runtimefile(filename: str) -> str:
     This **is**  for files needed at runtime. Files needed for the example notebooks
     are found via :func:`pint.config.examplefile`.
     """
-    return pkg_resources.resource_filename(
-        __name__, os.path.join("data/runtime/", filename)
-    )
+    return os.path.join(importlib.resources.files("pint"), f"data/runtime/{filename}")

src/pint/derived_quantities.py

@@ -136,6 +136,14 @@ def pferrs(
     return (forp, forperr, fdorpd, fdorpderr)
 
 
+def _to_gauss(B: u.Quantity) -> u.G:
+    """Convert quantity with mass, length, and time units to Gauss.
+
+    In cgs units, magnetic field is has units (mass/length)^(1/2) / time.
+    """
+    return B.to(u.Gauss, equivalencies=[pint.gauss_equiv])
+
+
 @u.quantity_input(f=u.Hz, fdot=u.Hz / u.s, fo=u.Hz)
 def pulsar_age(
     f: u.Quantity, fdot: u.Quantity, n: int = 3, fo: u.Quantity = 1e99 * u.Hz
@@ -219,19 +227,27 @@ def pulsar_edot(
     return (-4.0 * np.pi**2 * I * f * fdot).to(u.erg / u.s)
 
 
-@u.quantity_input(f=u.Hz, fdot=u.Hz / u.s)
-def pulsar_B(f: u.Quantity, fdot: u.Quantity) -> u.G:
+@u.quantity_input(f=u.Hz, fdot=u.Hz / u.s, I=u.g * u.cm**2, R=u.km)
+def pulsar_B(
+    f: u.Quantity,
+    fdot: u.Quantity,
+    I: u.Quantity = 1.0e45 * u.g * u.cm**2,
+    R: u.Quantity = 10 * u.km,
+) -> u.G:
     r"""Compute pulsar surface magnetic field
 
-    Return the estimated pulsar surface magnetic field strength
-    given the spin frequency and frequency derivative.
+    Return the pulsar surface magnetic field strength given the spin frequency `f` and frequency derivative `fdot`.
 
     Parameters
     ----------
     f : astropy.units.Quantity
         pulsar frequency
     fdot : astropy.units.Quantity
         frequency derivative :math:`\dot f`
+    I : astropy.units.Quantity, optional
+        pulsar moment of inertia, default of 1e45 g*cm**2
+    R : astropy.units.Quantity, optional
+        pulsar radius, default of 10 km
 
     Returns
     -------
@@ -247,15 +263,19 @@ def pulsar_B(f: u.Quantity, fdot: u.Quantity) -> u.G:
 
     Notes
     -----
-    Calculates :math:`B=3.2\times 10^{19}\,{\rm  G}\sqrt{ f \dot f^{-3}}`
+    Calculates :math:`B=\sqrt{\frac{3\,I\,c^3}{8\pi^2\,R^6}\times\frac{-\dot{f}}{f^3}}`
     """
-    # This is a hack to use the traditional formula by stripping the units.
-    # It would be nice to improve this to a  proper formula with units
-    return 3.2e19 * u.G * np.sqrt(-fdot.to_value(u.Hz / u.s) / f.to_value(u.Hz) ** 3.0)
+    factor = (3.0 * I * const.c**3) / (8.0 * np.pi**2 * R**6)
+    return _to_gauss((factor * (-fdot) / f**3) ** 0.5)
 
 
-@u.quantity_input(f=u.Hz, fdot=u.Hz / u.s)
-def pulsar_B_lightcyl(f: u.Quantity, fdot: u.Quantity) -> u.G:
+@u.quantity_input(f=u.Hz, fdot=u.Hz / u.s, I=u.g * u.cm**2, R=u.km)
+def pulsar_B_lightcyl(
+    f: u.Quantity,
+    fdot: u.Quantity,
+    I: u.Quantity = 1.0e45 * u.g * u.cm**2,
+    R=10 * u.km,
+) -> u.G:
     r"""Compute pulsar magnetic field at the light cylinder
 
     Return the estimated pulsar magnetic field strength at the
@@ -268,6 +288,10 @@ def pulsar_B_lightcyl(f: u.Quantity, fdot: u.Quantity) -> u.G:
         pulsar frequency
     fdot : astropy.units.Quantity
         frequency derivative :math:`\dot f`
+    I : astropy.units.Quantity, optional
+        pulsar moment of inertia, default of 1e45 g*cm**2
+    R : astropy.units.Quantity, optional
+        pulsar radius, default of 10 km
 
     Returns
     -------
@@ -283,17 +307,10 @@ def pulsar_B_lightcyl(f: u.Quantity, fdot: u.Quantity) -> u.G:
 
     Notes
     -----
-    Calculates :math:`B_{LC} = 2.9\times 10^8\,{\rm G} P^{-5/2} \dot P^{1/2}`
+    Calculates :math:`B_{LC} = \sqrt{\frac{-24\pi^4\,I}{c^3}\dot{f}f^3}`
     """
-    p, pd = p_to_f(f, fdot)
-    # This is a hack to use the traditional formula by stripping the units.
-    # It would be nice to improve this to a  proper formula with units
-    return (
-        2.9e8
-        * u.G
-        * p.to_value(u.s) ** (-5.0 / 2.0)
-        * np.sqrt(pd.to(u.dimensionless_unscaled).value)
-    )
+    factor = 24.0 * np.pi**4.0 * I / const.c**3.0
+    return _to_gauss((factor * (-fdot) * f**3.0) ** 0.5)
 
 
 @u.quantity_input(pb=u.d, x=u.cm)

src/pint/models/noise_model.py

@@ -126,6 +126,7 @@ def setup(self):
                             index=tneq_par.index,
                             aliases=["T2EQUAD"],
                             description="An error term added in quadrature to the scaled (by EFAC) TOA uncertainty.",
+                            convert_tcb2tdb=False,
                         )
                     )
                 EQUAD_par = getattr(self, EQUAD_name)
@@ -167,14 +168,14 @@ def scale_toa_sigma(self, toas, warn=True):
             if equad.quantity is None:
                 continue
             mask = equad.select_toa_mask(toas)
-            if np.any(mask):
+            if len(mask) > 0:
                 sigma_scaled[mask] = np.hypot(sigma_scaled[mask], equad.quantity)
             elif warn:
                 warnings.warn(f"EQUAD {equad} has no TOAs")
         for efac_name in self.EFACs:
             efac = getattr(self, efac_name)
             mask = efac.select_toa_mask(toas)
-            if np.any(mask):
+            if len(mask) > 0:
                 sigma_scaled[mask] *= efac.quantity
             elif warn:
                 warnings.warn(f"EFAC {efac} has no TOAs")

src/pint/models/parameter.py

@@ -709,7 +709,7 @@ def __init__(
 
         assert (
             not convert_tcb2tdb or tcb2tdb_scale_factor is not None
-        ), "Please specify the tcb2tdb_scale_factor explicitly."
+        ), f"Please specify the tcb2tdb_scale_factor explicitly for {name}."
         self.convert_tcb2tdb = convert_tcb2tdb
         self.tcb2tdb_scale_factor = tcb2tdb_scale_factor
 
@@ -1133,7 +1133,7 @@ def __init__(
 
         assert (
             not convert_tcb2tdb or tcb2tdb_scale_factor is not None
-        ), "Please specify the tcb2tdb_scale_factor explicitly."
+        ), f"Please specify the tcb2tdb_scale_factor explicitly for {name}."
         self.convert_tcb2tdb = convert_tcb2tdb
         self.tcb2tdb_scale_factor = tcb2tdb_scale_factor
 
@@ -1334,7 +1334,7 @@ def __init__(
 
         assert (
             not convert_tcb2tdb or tcb2tdb_scale_factor is not None
-        ), "Please specify the tcb2tdb_scale_factor explicitly."
+        ), f"Please specify the tcb2tdb_scale_factor explicitly for {name}."
         self.convert_tcb2tdb = convert_tcb2tdb
         self.tcb2tdb_scale_factor = tcb2tdb_scale_factor
 
@@ -1554,7 +1554,7 @@ def __init__(
         # a function of the prefix.
         assert (
             not convert_tcb2tdb or tcb2tdb_scale_factor is not None
-        ), "Please specify the tcb2tdb_scale_factor explicitly."
+        ), f"Please specify the tcb2tdb_scale_factor explicitly for {name}."
         tcb2tdb_scale_factor_val = (
             tcb2tdb_scale_factor(self.prefix)
             if hasattr(tcb2tdb_scale_factor, "__call__")

src/pint/models/timing_model.py

@@ -3425,6 +3425,7 @@ def add_param(self, param, deriv_func=None, setup=False):
         if deriv_func is not None:
             self.register_deriv_funcs(deriv_func, param.name)
         param._parent = self
+        return param.name
 
     def remove_param(self, param):
         """Remove a parameter from the Component.