diff --git a/.github/workflows/ci_test.yml b/.github/workflows/ci_test.yml index 65f603dba..dc096a29c 100644 --- a/.github/workflows/ci_test.yml +++ b/.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' @@ -41,18 +41,23 @@ jobs: # - os: ubuntu-latest # python: '3.10' # tox_env: 'py310-test-alldeps-cov' - - os: macos-12 - python: '3.12' - tox_env: 'py312-test' + # - os: macos-12 + # 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' # tox_env: 'codestyle' +# uncomment this if needed to run a single test quickly +# can specify OS and python versions + # - os: macos-12 + # python: '3.12' + # tox_env: 'singletest' steps: - name: Check out repository @@ -89,3 +94,47 @@ jobs: # with: # name: documentation # path: .tox/docs_out/ + macos-latest-py313: + runs-on: macos-latest + steps: + - uses: actions/checkout@v4 + - name: Set-Up Python Env + uses: mamba-org/setup-micromamba@v1 + with: + init-shell: bash + environment-name: pint + cache-environment: true + cache-downloads: true + create-args: >- + --platform osx-64 + -c conda-forge + python=3.13 + astropy + git + - name: Install base dependencies + shell: bash -el {0} + run: | + python -m pip install --upgrade pip + python -m pip install tox pytest hypothesis numdifftools pathos setuptools + - name: Print OS, machine info + shell: bash -el {0} + run: | + python -c "import os; print(f'os {os.uname()}')" + python -c "import platform; print(f'processor {platform.processor()}')" + python -c "import numpy as np; print(f'eps: {np.finfo(np.longdouble).eps}')" + - name: Print Python, pip, and tox versions + shell: bash -el {0} + run: | + python -c "import sys; print(f'Python {sys.version}')" + python -c "import pip; print(f'pip {pip.__version__}')" + python -c "import tox; print(f'tox {tox.__version__}')" + - name: Install PINT and requirements + shell: bash -el {0} + run: | + python -m pip install -r requirements.txt + python -m pip install -r requirements_dev.txt + python -m pip install --force-reinstall --no-deps . + - name: Run tests + shell: bash -el {0} + run: pytest -v --pyargs tests --reruns 5 + \ No newline at end of file diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 6c2c4b8cc..82fcce346 100644 --- a/.pre-commit-config.yaml +++ b/.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 diff --git a/AUTHORS.rst b/AUTHORS.rst index 1e2ed1489..00342b0b2 100644 --- a/AUTHORS.rst +++ b/AUTHORS.rst @@ -23,6 +23,7 @@ Active developers are indicated by (*). Authors of the PINT paper are indicated * Bastian Beranek * Deven Bhakta (*) * Chloe Champagne (#) +* Colin Clark * Jonathan Colen (#) * H Thankful Cromartie * Christoph Deil diff --git a/CHANGELOG-unreleased.md b/CHANGELOG-unreleased.md index 286321f2c..9169dc1f2 100644 --- a/CHANGELOG-unreleased.md +++ b/CHANGELOG-unreleased.md @@ -9,36 +9,10 @@ the released changes. ## Unreleased ### Changed -- Moved the events -> TOAs and photon weights code into the function `load_events_weights` within `event_optimize`. -- Updated the `maxMJD` argument in `event_optimize` to default to the current mjd -- `maskParameter.__repr__()` output now includes the frozen attribute. -- Changed default value of `FDJUMPLOG` to `Y` -- Bumped `black` version to 24.x ### Added -- arXiv link of PINT noise paper in README -- Type hints in `pint.derived_quantities`, `pint.modelutils`, `pint.binaryconvert`, `pint.config`, -`pint.erfautils`, `pint.fits_utils`, `pint.logging` and `pint.residuals` -- Doing `model.par = something` will try to assign to `par.quantity` or `par.value` but will give warning - `plrednoise_from_wavex()` and `pldmnoise_from_dmwavex()` functions now compute `TNRedFLow` and `TNDMFLow` - `powerlaw_corner` function - `TNREDFLOW` and `TNREDCORNER` parameters in `PLRedNoise` - `TNDMFLOW` and `TNDMCORNER` parameters in `PLDMNoise` -- `PLChromNoise` component to model chromatic red noise with a power law spectrum -- Fourier series representation of chromatic noise (`CMWaveX`) -- `pint.utils.cmwavex_setup` and `pint.utils.plchromnoise_from_cmwavex` functions -- More validation for correlated noise components in `TimingModel.validate_component_types()` ### Fixed -- Bug in `DMWaveX.get_indices()` function -- Explicit type conversion in `woodbury_dot()` function -- Documentation: Fixed empty descriptions in the timing model components table -- BIC implementation -- `event_optimize`: Fixed a bug that was causing the results.txt file to be written without the median values. -- SWX model now has SWXP_0001 frozen by default, and new segments should also have SWXP frozen -- Can now properly use local files for ephemeris -- Typos in `explanation.rst` regarding local ephemeris. ### Removed -- Removed the argument `--usepickle` in `event_optimize` as the `load_events_weights` function checks the events file type to see if the -file is a pickle file. -- Removed obsolete code, such as manually tracking the progress of the MCMC run within `event_optimize` -- Unnecessary default arguments from the `powerlaw()` function. -- `download_data.sh` script and `de432s.bsp` ephemeris file diff --git a/CHANGELOG.md b/CHANGELOG.md index 80348c541..8131c21b4 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,6 +7,73 @@ 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.1] 2024-012-18 +### 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 +- Time derivatives of NE_SW in `SolarWindDispersion` +- New prefix pattern for `split_prefixed_name` to handle derivatives of NE_SW +- 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 +- Piecewise-constant model for chromatic variations (CMX) +- `add_param` returns the name of the parameter (useful for numbered parameters) +- Rerun intermittent failures in CI +- micromamba CI environment for testing macOS-latest, without tox +- models now have metadata dictionary +### 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 +- Condense code in Glitch model and add test coverage. +- `find_empty_masks` will now search through `CMX` parameters +- Fixed some docstrings for binary models. +### Removed +- macOS 12 CI + +## [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`. +- Updated the `maxMJD` argument in `event_optimize` to default to the current mjd +- `maskParameter.__repr__()` output now includes the frozen attribute. +- Changed default value of `FDJUMPLOG` to `Y` +- Bumped `black` version to 24.x +- Moved all custom exceptions and warnings to a single module `pint.exceptions` +- Changed from `setup.cfg` to `pyproject.toml` +### Added +- arXiv link of PINT noise paper in README +- Type hints in `pint.derived_quantities`, `pint.modelutils`, `pint.binaryconvert`, `pint.config`, +`pint.erfautils`, `pint.fits_utils`, `pint.logging` and `pint.residuals` +- Doing `model.par = something` will try to assign to `par.quantity` or `par.value` but will give warning +- `PLChromNoise` component to model chromatic red noise with a power law spectrum +- Fourier series representation of chromatic noise (`CMWaveX`) +- `pint.utils.cmwavex_setup` and `pint.utils.plchromnoise_from_cmwavex` functions +- More validation for correlated noise components in `TimingModel.validate_component_types()` +- ORBWAVEs model for modelling binary orbital period variations in the fourier domain +### Fixed +- Bug in `DMWaveX.get_indices()` function +- Explicit type conversion in `woodbury_dot()` function +- Documentation: Fixed empty descriptions in the timing model components table +- BIC implementation +- `event_optimize`: Fixed a bug that was causing the results.txt file to be written without the median values. +- SWX model now has SWXP_0001 frozen by default, and new segments should also have SWXP frozen +- Can now properly use local files for ephemeris +- Typos in `explanation.rst` regarding local ephemeris. +- DD/ELL1 models will check for valid SINI and raise exception if it strays, which will tell the fitter to try elsewhere +- Don't try to print AIC for wideband data in `pintk` (it's not yet implemented) +### Removed +- Removed the argument `--usepickle` in `event_optimize` as the `load_events_weights` function checks the events file type to see if the file is a pickle file. +- Removed obsolete code, such as manually tracking the progress of the MCMC run within `event_optimize` +- `download_data.sh` script and `de432s.bsp` ephemeris file + ## [1.0.1] 2024-07-01 ### Changed - Avoided unnecessary creation of `SkyCoord` objects in `AstrometryEquatorial` and `AstrometryEcliptic`. diff --git a/README.rst b/README.rst index 76085fe66..75b7862a1 100644 --- a/README.rst +++ b/README.rst @@ -127,6 +127,7 @@ email pint@nanograv.org or one of the people below: * Scott Ransom (sransom@nrao.edu) * Paul Ray (paul.s.ray3.civ@us.navy.mil) * David Kaplan (kaplan@uwm.edu) +* Abhimanyu Susobhanan (abhimanyu.susobhanan@nanograv.org) Want to do something new? Submit a github `issue `_. diff --git a/docs/installation.rst b/docs/installation.rst index eb780c9ea..c6e23b549 100644 --- a/docs/installation.rst +++ b/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``). diff --git a/pyproject.toml b/pyproject.toml new file mode 100644 index 000000000..fa5ebdbe7 --- /dev/null +++ b/pyproject.toml @@ -0,0 +1,144 @@ +[build-system] +requires = ["setuptools>=61.2", "versioneer"] +build-backend = "setuptools.build_meta" + +[project] +name = "pint-pulsar" +description = "A Pulsar Timing Package, written in Python from scratch" +authors = [ + {name = "Luo Jing", email = "sransom@nrao.edu"}, + {name = "Scott Ransom"}, + {name = "Paul Demorest"}, + {name = "Paul Ray"}, + {name = "et al."}, +] +license = {text = "License :: OSI Approved :: BSD License"} +classifiers = [ + "Intended Audience :: Science/Research", + "License :: OSI Approved :: BSD License", + "Operating System :: OS Independent", + "Programming Language :: Python", + "Programming Language :: Python :: 3", + "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.9" +dependencies = [ + "numpy>=1.23.0", + "astropy>=5.0.5", + "pyerfa", + "scipy>=1.9.0", + "jplephem>=2.6", + "matplotlib>=3.4.3", + "emcee>=3.0.1", + "corner>=2.0.1", + "uncertainties", + "loguru", + "nestle>=0.2.0", + "numdifftools", +] +dynamic = ["version"] + +[project.readme] +file = "README.rst" +content-type = "text/x-rst" + +[project.urls] +Homepage = "https://github.com/nanograv/PINT" +Documentation = "https://nanograv-pint.readthedocs.io/" + +[project.entry-points] +# See the docstring in versioneer.py for instructions. Note that you must +# re-run 'versioneer.py setup' after changing this section, and commit the +# resulting files. + +[project.scripts] +photonphase = "pint.scripts.photonphase:main" +event_optimize = "pint.scripts.event_optimize:main" +event_optimize_multiple = "pint.scripts.event_optimize_multiple:main" +pintempo = "pint.scripts.pintempo:main" +zima = "pint.scripts.zima:main" +pintbary = "pint.scripts.pintbary:main" +fermiphase = "pint.scripts.fermiphase:main" +pintk = "pint.scripts.pintk:main" +convert_parfile = "pint.scripts.convert_parfile:main" +compare_parfiles = "pint.scripts.compare_parfiles:main" +tcb2tdb = "pint.scripts.tcb2tdb:main" +t2binary2pint = "pint.scripts.t2binary2pint:main" +pintpublish = "pint.scripts.pintpublish:main" + +[tool.setuptools] +zip-safe = false +package-dir = {"" = "src"} +include-package-data = true +# These should match requirements.txt + +[tool.setuptools.packages.find] +where = ["src"] +namespaces = false + +[tool.setuptools.package-data] +"*" = ["*.*"] + +[tool.versioneer] +VCS = "git" +style = "pep440" +versionfile_source = "src/pint/extern/_version.py" +versionfile_build = "pint/extern/_version.py" +tag_prefix = "''" +parentdir_prefix = "'pint-'" + +[tool.distutils.bdist_wheel] +universal = 0 + +[tool.aliases] +test = "pytest" + +[tool.flake8] +max-line-length = "100" +# This is an inappropriate non-error for slicing +extend-ignore = """ +E203, +E265 +# __init__ doesn't need a docstring, should be in the class +D107 +# Other magic methods don't necessarily need docstings, what they do is well-defined +D105 +# Style issues, suppress these for a full flake8 run +# E111,E114,E115,E116,E122,E123,E124,E125,E126,E127,E128,E129,E131 +# E201,E202,E203,E221,E222,E225,E226,E227,E228,E231,E241,E251,E261,E262,E265,E266,E271,E272 +# E301,E302,E303,E305,E306 +# E401,E501,E502 +# E701,E702,E703,E704,E741 +# W291,W293,W391,W503,W504 +# D100,D101,D102,D103,D104,D105 +# D200,D202,D204,D205,D207,D208,D209,D210 +# D300 +# D400,D401,D402,D403,D412,D413 +# RST201,RST202,RST203,RST210,RST212,RST299 +# RST301,RST304,RST306 +# Ugh people want to break these rules +N802 # Function names should be lowercase +N803 # argument name should be lowercase +N806 # variable should be lowercase""" +statistics = "True" +exclude = """ +docs/conf.py +versioneer.py +pint/mcmc_fitter.py""" +rst-roles = """ +class, +module, +func,""" + +[tool.isort] +multi_line_output = 3 +line_length = 88 +skip_glob = ["src/pint/extern/*"] +include_trailing_comma = true +combine_as_imports = true diff --git a/requirements.txt b/requirements.txt index 16f76401e..293e3d07a 100644 --- a/requirements.txt +++ b/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 diff --git a/requirements_dev.txt b/requirements_dev.txt index a00511840..be6046e27 100644 --- a/requirements_dev.txt +++ b/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 @@ -12,6 +11,8 @@ wheel>=0.29.0 pytest>=4.3 pytest-cov>=2.7.1 pytest-runner>=5.1 +pytest-xdist +pytest-rerunfailures flake8>=3.7 pep8-naming>=0.8.2 flake8-docstrings>=1.4 @@ -28,9 +29,7 @@ jupytext pdbpp tox pre-commit -typed-ast>=1.5.0 -#black>=19.0a0,<20.0a0 -black~=23.0 +black~=24.0 pygments ipython pathlib2 @@ -40,4 +39,3 @@ loguru # click<=8.0.4 gprof2dot py-cpuinfo -pytest-xdist diff --git a/setup.cfg b/setup.cfg deleted file mode 100644 index 28dfa2974..000000000 --- a/setup.cfg +++ /dev/null @@ -1,132 +0,0 @@ -[metadata] -name = pint-pulsar -description = A Pulsar Timing Package, written in Python from scratch -long_description = file: README.rst -long_description_content_type = text/x-rst -author = Luo Jing, Scott Ransom, Paul Demorest, Paul Ray, et al. -author_email = sransom@nrao.edu -url = https://github.com/nanograv/PINT -project_urls = - Documentation = https://nanograv-pint.readthedocs.io/ -license = License :: OSI Approved :: BSD License -classifier = - Intended Audience :: Science/Research - License :: OSI Approved :: BSD License - 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 - Topic :: Scientific/Engineering :: Astronomy - Topic :: Software Development :: Libraries :: Python Modules - -[options] -zip_safe = False -packages = find: -package_dir = - = src -include_package_data = True -python_requires = >=3.8 -# These should match requirements.txt -install_requires = - numpy>=1.18.5 - astropy>=4.0,!=4.0.1,!=4.0.1.post1 - pyerfa - scipy>=1.4.1 - jplephem>=2.6 - matplotlib>=3.2.0 - emcee>=3.0.1 - corner>=2.0.1 - uncertainties - loguru - nestle>=0.2.0 - numdifftools - -[options.packages.find] -where = src - -[options.package_data] -* = *.* - -[options.entry_points] -console_scripts = - photonphase = pint.scripts.photonphase:main - event_optimize = pint.scripts.event_optimize:main - event_optimize_multiple = pint.scripts.event_optimize_multiple:main - pintempo = pint.scripts.pintempo:main - zima = pint.scripts.zima:main - pintbary = pint.scripts.pintbary:main - fermiphase = pint.scripts.fermiphase:main - pintk = pint.scripts.pintk:main - convert_parfile = pint.scripts.convert_parfile:main - compare_parfiles = pint.scripts.compare_parfiles:main - tcb2tdb = pint.scripts.tcb2tdb:main - t2binary2pint = pint.scripts.t2binary2pint:main - pintpublish = pint.scripts.pintpublish:main - - -# See the docstring in versioneer.py for instructions. Note that you must -# re-run 'versioneer.py setup' after changing this section, and commit the -# resulting files. - -[versioneer] -VCS = git -style = pep440 -versionfile_source = src/pint/extern/_version.py -versionfile_build = pint/extern/_version.py -tag_prefix = '' -parentdir_prefix = 'pint-' - -[bdist_wheel] -universal = 0 - -[aliases] -test=pytest - -[flake8] -max-line-length = 100 -# This is an inappropriate non-error for slicing -extend-ignore = E203, - E265 -# __init__ doesn't need a docstring, should be in the class - D107 -# Other magic methods don't necessarily need docstings, what they do is well-defined - D105 -# Style issues, suppress these for a full flake8 run -# E111,E114,E115,E116,E122,E123,E124,E125,E126,E127,E128,E129,E131 -# E201,E202,E203,E221,E222,E225,E226,E227,E228,E231,E241,E251,E261,E262,E265,E266,E271,E272 -# E301,E302,E303,E305,E306 -# E401,E501,E502 -# E701,E702,E703,E704,E741 -# W291,W293,W391,W503,W504 -# D100,D101,D102,D103,D104,D105 -# D200,D202,D204,D205,D207,D208,D209,D210 -# D300 -# D400,D401,D402,D403,D412,D413 -# RST201,RST202,RST203,RST210,RST212,RST299 -# RST301,RST304,RST306 -# Ugh people want to break these rules - N802 # Function names should be lowercase - N803 # argument name should be lowercase - N806 # variable should be lowercase - -statistics = True -exclude = - docs/conf.py - versioneer.py - pint/mcmc_fitter.py -rst-roles = - class, - module, - func, - -[isort] -multi_line_output = 3 -line_length = 88 -skip_glob = - src/pint/extern/* -include_trailing_comma = True -combine_as_imports = True diff --git a/src/pint/__init__.py b/src/pint/__init__.py index c0d511e7c..7906be01f 100644 --- a/src/pint/__init__.py +++ b/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: diff --git a/src/pint/config.py b/src/pint/config.py index a24bec249..344bc71b1 100644 --- a/src/pint/config.py +++ b/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}") diff --git a/src/pint/derived_quantities.py b/src/pint/derived_quantities.py index d5852ae4b..24dcb1934 100644 --- a/src/pint/derived_quantities.py +++ b/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,12 +227,16 @@ 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 ---------- @@ -232,6 +244,10 @@ def pulsar_B(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 ------- @@ -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) diff --git a/src/pint/exceptions.py b/src/pint/exceptions.py new file mode 100644 index 000000000..81b11a606 --- /dev/null +++ b/src/pint/exceptions.py @@ -0,0 +1,177 @@ +__all__ = [ + "DegeneracyWarning", + "ConvergenceFailure", + "MaxiterReached", + "StepProblem", + "CorrelatedErrors", + "MissingTOAs", + "PropertyAttributeError", + "TimingModelError", + "MissingParameter", + "AliasConflict", + "UnknownParameter", + "UnknownBinaryModel", + "MissingBinaryError", + "PINTPrecisionError", + "PrefixError", + "InvalidModelParameters", + "ComponentConflict", + "ClockCorrectionError", + "NoClockCorrections", + "ClockCorrectionOutOfRange", +] + + +# originally from fitter.py +class DegeneracyWarning(UserWarning): + pass + + +class ConvergenceFailure(ValueError): + pass + + +class MaxiterReached(ConvergenceFailure): + pass + + +class StepProblem(ConvergenceFailure): + pass + + +class CorrelatedErrors(ValueError): + def __init__(self, model): + trouble_components = [ + c.__class__.__name__ + for c in model.NoiseComponent_list + if c.introduces_correlated_errors + ] + super().__init__( + f"Model has correlated errors and requires a GLS-based fitter; " + f"remove {trouble_components} if you want to use WLS" + ) + self.trouble_components = trouble_components + + +# from timing_model.py +class MissingTOAs(ValueError): + """Some parameter does not describe any TOAs.""" + + def __init__(self, parameter_names): + if isinstance(parameter_names, str): + parameter_names = [parameter_names] + if len(parameter_names) == 1: + msg = f"Parameter {parameter_names[0]} does not correspond to any TOAs: you might need to run `model.find_empty_masks(toas, freeze=True)`" + elif len(parameter_names) > 1: + msg = f"Parameters {' '.join(parameter_names)} do not correspond to any TOAs: you might need to run `model.find_empty_masks(toas, freeze=True)`" + else: + raise ValueError("Incorrect attempt to construct MissingTOAs") + super().__init__(msg) + self.parameter_names = parameter_names + + +class PropertyAttributeError(ValueError): + pass + + +class TimingModelError(ValueError): + """Generic base class for timing model errors.""" + + pass + + +class MissingParameter(TimingModelError): + """A required model parameter was not included. + + Parameters + ---------- + module + name of the model class that raised the error + param + name of the missing parameter + msg + additional message + + """ + + def __init__(self, module, param, msg=None): + super().__init__(msg) + self.module = module + self.param = param + self.msg = msg + + def __str__(self): + result = f"{self.module}.{self.param}" + if self.msg is not None: + result += "\n " + self.msg + return result + + +class AliasConflict(TimingModelError): + """If the same alias is used for different parameters.""" + + pass + + +class UnknownParameter(TimingModelError): + """Signal that a parameter name does not match any PINT parameters and their aliases.""" + + pass + + +class UnknownBinaryModel(TimingModelError): + """Signal that the par file requested a binary model not in PINT.""" + + def __init__(self, message, suggestion=None): + super().__init__(message) + self.suggestion = suggestion + + def __str__(self): + base_message = super().__str__() + if self.suggestion: + return f"{base_message} Perhaps use {self.suggestion}?" + return base_message + + +class MissingBinaryError(TimingModelError): + """Error for missing BINARY parameter.""" + + pass + + +# from utils.py +class PINTPrecisionError(RuntimeError): + pass + + +class PrefixError(ValueError): + pass + + +# from models.parameter.py +class InvalidModelParameters(ValueError): + pass + + +# models.model_builder.py +class ComponentConflict(ValueError): + """Error for multiple components can be select but no other indications.""" + + +# observatories.__init__.py +class ClockCorrectionError(RuntimeError): + """Unspecified error doing clock correction.""" + + pass + + +class NoClockCorrections(ClockCorrectionError): + """Clock corrections are expected but none are available.""" + + pass + + +class ClockCorrectionOutOfRange(ClockCorrectionError): + """Clock corrections are available but the requested time is not covered.""" + + pass diff --git a/src/pint/fitter.py b/src/pint/fitter.py index 1474b6b44..b187fad65 100644 --- a/src/pint/fitter.py +++ b/src/pint/fitter.py @@ -76,6 +76,7 @@ AngleParameter, boolParameter, strParameter, + InvalidModelParameters, ) from pint.pint_matrix import ( CorrelationMatrix, @@ -89,6 +90,13 @@ from pint.residuals import Residuals, WidebandTOAResiduals from pint.toa import TOAs from pint.utils import FTest, normalize_designmatrix +from pint.exceptions import ( + DegeneracyWarning, + ConvergenceFailure, + MaxiterReached, + StepProblem, + CorrelatedErrors, +) __all__ = [ @@ -166,22 +174,6 @@ def __get__(self, instance, owner=None): return val -class DegeneracyWarning(UserWarning): - pass - - -class ConvergenceFailure(ValueError): - pass - - -class MaxiterReached(ConvergenceFailure): - pass - - -class StepProblem(ConvergenceFailure): - pass - - class Fitter: """Base class for objects encapsulating fitting problems. @@ -887,24 +879,6 @@ def covariance_matrix(self): return self.parameter_covariance_matrix -class InvalidModelParameters(ValueError): - pass - - -class CorrelatedErrors(ValueError): - def __init__(self, model): - trouble_components = [ - c.__class__.__name__ - for c in model.NoiseComponent_list - if c.introduces_correlated_errors - ] - super().__init__( - f"Model has correlated errors and requires a GLS-based fitter; " - f"remove {trouble_components} if you want to use WLS" - ) - self.trouble_components = trouble_components - - class ModelState: """Record a model state and cache calculations diff --git a/src/pint/models/__init__.py b/src/pint/models/__init__.py index 67201a259..27c78a4f4 100644 --- a/src/pint/models/__init__.py +++ b/src/pint/models/__init__.py @@ -25,7 +25,7 @@ from pint.models.binary_dd import BinaryDD, BinaryDDS, BinaryDDGR, BinaryDDH from pint.models.binary_ddk import BinaryDDK from pint.models.binary_ell1 import BinaryELL1, BinaryELL1H, BinaryELL1k -from pint.models.chromatic_model import ChromaticCM +from pint.models.chromatic_model import ChromaticCM, ChromaticCMX from pint.models.cmwavex import CMWaveX from pint.models.dispersion_model import ( DispersionDM, diff --git a/src/pint/models/absolute_phase.py b/src/pint/models/absolute_phase.py index 7228c255c..fa1acb741 100644 --- a/src/pint/models/absolute_phase.py +++ b/src/pint/models/absolute_phase.py @@ -5,7 +5,8 @@ import pint.toa as toa from pint.models.parameter import MJDParameter, floatParameter, strParameter -from pint.models.timing_model import MissingParameter, PhaseComponent +from pint.models.timing_model import PhaseComponent +from pint.exceptions import MissingParameter class AbsPhase(PhaseComponent): diff --git a/src/pint/models/astrometry.py b/src/pint/models/astrometry.py index 1ac7e4a90..2f1d4feee 100644 --- a/src/pint/models/astrometry.py +++ b/src/pint/models/astrometry.py @@ -24,9 +24,10 @@ strParameter, ) import pint.toa -from pint.models.timing_model import DelayComponent, MissingParameter +from pint.models.timing_model import DelayComponent from pint.pulsar_ecliptic import OBL, PulsarEcliptic from pint.utils import add_dummy_distance, remove_dummy_distance +from pint.exceptions import MissingParameter astropy_version = sys.modules["astropy"].__version__ mas_yr = u.mas / u.yr diff --git a/src/pint/models/binary_bt.py b/src/pint/models/binary_bt.py index 8c4969c2d..03ad52583 100644 --- a/src/pint/models/binary_bt.py +++ b/src/pint/models/binary_bt.py @@ -5,7 +5,7 @@ from pint.models.pulsar_binary import PulsarBinary from pint.models.stand_alone_psr_binaries.BT_model import BTmodel from pint.models.stand_alone_psr_binaries.BT_piecewise import BTpiecewise -from pint.models.timing_model import MissingParameter +from pint.exceptions import MissingParameter import astropy.units as u import astropy.constants as consts from pint import ls @@ -21,7 +21,7 @@ class BinaryBT(PulsarBinary): """Blandford and Teukolsky binary model. - This binary model is described in Blandford and Teukolshy 1976. It is + This binary model is described in Blandford and Teukolsky (1976). It is a relatively simple parametrized post-Keplerian model that does not support Shapiro delay calculations. @@ -36,12 +36,16 @@ class BinaryBT(PulsarBinary): Notes ----- Because PINT's binary models all support specification of multiple orbital - frequency derivatives FBn, this is capable of behaving like the model called - BTX in tempo2. The model called BTX in tempo instead supports multiple - (non-interacting) companions, and that is not supported here. Neither can - PINT accept "BTX" as an alias for this model. + frequency derivatives ``FBn``, this is capable of behaving like the model called + ``BTX`` in ``tempo2``. The model called ``BTX`` in ``tempo`` instead supports multiple + (non-interacting) companions, and that is not supported here. + + References + ---------- + - Blandford & Teukolsky 1976, ApJ, 205, 580 [1]_ + + .. [1] https://ui.adsabs.harvard.edu/abs/1976ApJ...205..580B/abstract - See Blandford & Teukolsky 1976, ApJ, 205, 580. """ register = True @@ -82,16 +86,23 @@ def validate(self): class BinaryBTPiecewise(PulsarBinary): - """Model implementing the BT model with piecewise orbital parameters A1X and T0X. This model lets the user specify time ranges and fit for a different piecewise orbital parameter in each time range, - This is a PINT pulsar binary BTPiecewise model class, a subclass of PulsarBinary. - It is a wrapper for stand alone BTPiecewise class defined in - ./stand_alone_psr_binary/BT_piecewise.py - The aim for this class is to connect the stand alone binary model with the PINT platform. - BTpiecewise special parameters, where xxxx denotes the 4-digit index of the piece: - T0X_xxxx Piecewise T0 values for piece - A1X_xxxx Piecewise A1 values for piece - XR1_xxxx Lower time boundary of piece - XR2_xxxx Upper time boundary of piece + """BT model with piecewise orbital parameters ``A1X`` and ``T0X``. This model lets the user specify time ranges and fit for a different piecewise orbital parameter in each time range. + + ``BTpiecewise`` special parameters, where xxxx denotes the 4-digit index of the piece: + + - ``T0X_xxxx``: Piecewise ``T0`` values for piece + - ``A1X_xxxx``: Piecewise ``A1`` values for piece + - ``XR1_xxxx``: Lower time boundary of piece + - ``XR2_xxxx``: Upper time boundary of piece + + The actual calculations for this are done in + :class:`pint.models.stand_alone_psr_binaries.BT_piecewise.BTpiecewise` + + Parameters supported: + + .. paramtable:: + :class: pint.models.binary_bt.BinaryBTPiecewise + """ register = True diff --git a/src/pint/models/binary_ddk.py b/src/pint/models/binary_ddk.py index 37e31856c..9d83403a1 100644 --- a/src/pint/models/binary_ddk.py +++ b/src/pint/models/binary_ddk.py @@ -8,7 +8,7 @@ from pint.models.binary_dd import BinaryDD from pint.models.parameter import boolParameter, floatParameter, funcParameter from pint.models.stand_alone_psr_binaries.DDK_model import DDKmodel -from pint.models.timing_model import MissingParameter, TimingModelError +from pint.exceptions import MissingParameter, TimingModelError def _convert_kin(kin): diff --git a/src/pint/models/binary_ell1.py b/src/pint/models/binary_ell1.py index a34b3c186..bec23f0b3 100644 --- a/src/pint/models/binary_ell1.py +++ b/src/pint/models/binary_ell1.py @@ -17,7 +17,7 @@ from pint.models.stand_alone_psr_binaries.ELL1_model import ELL1model from pint.models.stand_alone_psr_binaries.ELL1H_model import ELL1Hmodel from pint.models.stand_alone_psr_binaries.ELL1k_model import ELL1kmodel -from pint.models.timing_model import MissingParameter +from pint.exceptions import MissingParameter from pint.utils import taylor_horner_deriv from pint import Tsun diff --git a/src/pint/models/chromatic_model.py b/src/pint/models/chromatic_model.py index 49a5c152a..8f428fea2 100644 --- a/src/pint/models/chromatic_model.py +++ b/src/pint/models/chromatic_model.py @@ -2,17 +2,20 @@ from warnings import warn import numpy as np import astropy.units as u -from pint.models.timing_model import DelayComponent, MissingParameter +from pint.models.timing_model import DelayComponent, MissingParameter, MissingTOAs from pint.models.parameter import floatParameter, prefixParameter, MJDParameter +from pint.toa_select import TOASelect from pint.utils import split_prefixed_name, taylor_horner, taylor_horner_deriv from pint import DMconst +from pint.exceptions import MissingParameter from astropy.time import Time +from loguru import logger as log cmu = u.pc / u.cm**3 / u.MHz**2 class Chromatic(DelayComponent): - """A base chromatic timing model.""" + """A base chromatic timing model with a constant chromatic index.""" def __init__(self): super().__init__() @@ -111,11 +114,14 @@ def register_cm_deriv_funcs(self, func, param): class ChromaticCM(Chromatic): - """Simple chromatic delay model. + """Simple chromatic delay model with a constant chromatic index. This model uses Taylor expansion to model CM variation over time. It can also be used for a constant CM. + Fitting for the chromatic index is not supported because the fit is too + unstable when fit simultaneously with the DM. + Parameters supported: .. paramtable:: @@ -300,3 +306,401 @@ def change_cmepoch(self, new_epoch): dt.to(u.yr), cmterms, deriv_order=n + 1 ) self.CMEPOCH.value = new_epoch + + +class ChromaticCMX(Chromatic): + """This class provides a CMX model - piecewise-constant chromatic variations with constant + chromatic index. + + This model lets the user specify time ranges and fit for a different CMX value in each time range. + + It should be used in combination with the `ChromaticCM` model. Specifically, TNCHROMIDX must be + set. + + Parameters supported: + + .. paramtable:: + :class: pint.models.chromatic_model.ChromaticCMX + """ + + register = True + category = "chromatic_cmx" + + def __init__(self): + super().__init__() + + self.add_CMX_range(None, None, cmx=0, frozen=False, index=1) + + self.cm_value_funcs += [self.cmx_cm] + self.set_special_params(["CMX_0001", "CMXR1_0001", "CMXR2_0001"]) + self.delay_funcs_component += [self.CMX_chromatic_delay] + + def add_CMX_range(self, mjd_start, mjd_end, index=None, cmx=0, frozen=True): + """Add CMX range to a chromatic model with specified start/end MJDs and CMX value. + + Parameters + ---------- + + mjd_start : float or astropy.quantity.Quantity or astropy.time.Time + MJD for beginning of CMX event. + mjd_end : float or astropy.quantity.Quantity or astropy.time.Time + MJD for end of CMX event. + index : int, None + Integer label for CMX event. If None, will increment largest used index by 1. + cmx : float or astropy.quantity.Quantity + Change in CM during CMX event. + frozen : bool + Indicates whether CMX will be fit. + + Returns + ------- + + index : int + Index that has been assigned to new CMX event. + """ + + #### Setting up the CMX title convention. If index is None, want to increment the current max CMX index by 1. + if index is None: + dct = self.get_prefix_mapping_component("CMX_") + index = np.max(list(dct.keys())) + 1 + i = f"{int(index):04d}" + + if mjd_end is not None and mjd_start is not None: + if mjd_end < mjd_start: + raise ValueError("Starting MJD is greater than ending MJD.") + elif mjd_start != mjd_end: + raise ValueError("Only one MJD bound is set.") + + if int(index) in self.get_prefix_mapping_component("CMX_"): + raise ValueError( + f"Index '{index}' is already in use in this model. Please choose another." + ) + + if isinstance(cmx, u.quantity.Quantity): + cmx = cmx.to_value(cmu) + + if isinstance(mjd_start, Time): + mjd_start = mjd_start.mjd + elif isinstance(mjd_start, u.quantity.Quantity): + mjd_start = mjd_start.value + if isinstance(mjd_end, Time): + mjd_end = mjd_end.mjd + elif isinstance(mjd_end, u.quantity.Quantity): + mjd_end = mjd_end.value + self.add_param( + prefixParameter( + name=f"CMX_{i}", + units=cmu, + value=cmx, + description="Dispersion measure variation", + parameter_type="float", + frozen=frozen, + convert_tcb2tdb=False, + ) + ) + self.add_param( + prefixParameter( + name=f"CMXR1_{i}", + units="MJD", + description="Beginning of CMX interval", + parameter_type="MJD", + time_scale="utc", + value=mjd_start, + convert_tcb2tdb=False, + ) + ) + self.add_param( + prefixParameter( + name=f"CMXR2_{i}", + units="MJD", + description="End of CMX interval", + parameter_type="MJD", + time_scale="utc", + value=mjd_end, + convert_tcb2tdb=False, + ) + ) + self.setup() + self.validate() + return index + + def add_CMX_ranges(self, mjd_starts, mjd_ends, indices=None, cmxs=0, frozens=True): + """Add CMX ranges to a dispersion model with specified start/end MJDs and CMXs. + + Parameters + ---------- + + mjd_starts : iterable of float or astropy.quantity.Quantity or astropy.time.Time + MJD for beginning of CMX event. + mjd_end : iterable of float or astropy.quantity.Quantity or astropy.time.Time + MJD for end of CMX event. + indices : iterable of int, None + Integer label for CMX event. If None, will increment largest used index by 1. + cmxs : iterable of float or astropy.quantity.Quantity, or float or astropy.quantity.Quantity + Change in CM during CMX event. + frozens : iterable of bool or bool + Indicates whether CMX will be fit. + + Returns + ------- + + indices : list + Indices that has been assigned to new CMX events + """ + if len(mjd_starts) != len(mjd_ends): + raise ValueError( + f"Number of mjd_start values {len(mjd_starts)} must match number of mjd_end values {len(mjd_ends)}" + ) + if indices is None: + indices = [None] * len(mjd_starts) + cmxs = np.atleast_1d(cmxs) + if len(cmxs) == 1: + cmxs = np.repeat(cmxs, len(mjd_starts)) + if len(cmxs) != len(mjd_starts): + raise ValueError( + f"Number of mjd_start values {len(mjd_starts)} must match number of cmx values {len(cmxs)}" + ) + frozens = np.atleast_1d(frozens) + if len(frozens) == 1: + frozens = np.repeat(frozens, len(mjd_starts)) + if len(frozens) != len(mjd_starts): + raise ValueError( + f"Number of mjd_start values {len(mjd_starts)} must match number of frozen values {len(frozens)}" + ) + + #### Setting up the CMX title convention. If index is None, want to increment the current max CMX index by 1. + dct = self.get_prefix_mapping_component("CMX_") + last_index = np.max(list(dct.keys())) + added_indices = [] + for mjd_start, mjd_end, index, cmx, frozen in zip( + mjd_starts, mjd_ends, indices, cmxs, frozens + ): + if index is None: + index = last_index + 1 + last_index += 1 + elif index in list(dct.keys()): + raise ValueError( + f"Attempting to insert CMX_{index:04d} but it already exists" + ) + added_indices.append(index) + i = f"{int(index):04d}" + + if mjd_end is not None and mjd_start is not None: + if mjd_end < mjd_start: + raise ValueError("Starting MJD is greater than ending MJD.") + elif mjd_start != mjd_end: + raise ValueError("Only one MJD bound is set.") + if int(index) in dct: + raise ValueError( + f"Index '{index}' is already in use in this model. Please choose another." + ) + if isinstance(cmx, u.quantity.Quantity): + cmx = cmx.to_value(u.pc / u.cm**3) + if isinstance(mjd_start, Time): + mjd_start = mjd_start.mjd + elif isinstance(mjd_start, u.quantity.Quantity): + mjd_start = mjd_start.value + if isinstance(mjd_end, Time): + mjd_end = mjd_end.mjd + elif isinstance(mjd_end, u.quantity.Quantity): + mjd_end = mjd_end.value + log.trace(f"Adding CMX_{i} from MJD {mjd_start} to MJD {mjd_end}") + self.add_param( + prefixParameter( + name=f"CMX_{i}", + units=cmu, + value=cmx, + description="Dispersion measure variation", + parameter_type="float", + frozen=frozen, + convert_tcb2tdb=False, + ) + ) + self.add_param( + prefixParameter( + name=f"CMXR1_{i}", + units="MJD", + description="Beginning of CMX interval", + parameter_type="MJD", + time_scale="utc", + value=mjd_start, + convert_tcb2tdb=False, + ) + ) + self.add_param( + prefixParameter( + name=f"CMXR2_{i}", + units="MJD", + description="End of CMX interval", + parameter_type="MJD", + time_scale="utc", + value=mjd_end, + convert_tcb2tdb=False, + ) + ) + self.setup() + self.validate() + return added_indices + + def remove_CMX_range(self, index): + """Removes all CMX parameters associated with a given index/list of indices. + + Parameters + ---------- + + index : float, int, list, np.ndarray + Number or list/array of numbers corresponding to CMX indices to be removed from model. + """ + + if isinstance(index, (int, float, np.int64)): + indices = [index] + elif isinstance(index, (list, set, np.ndarray)): + indices = index + else: + raise TypeError( + f"index must be a float, int, set, list, or array - not {type(index)}" + ) + for index in indices: + index_rf = f"{int(index):04d}" + for prefix in ["CMX_", "CMXR1_", "CMXR2_"]: + self.remove_param(prefix + index_rf) + self.validate() + + def get_indices(self): + """Returns an array of integers corresponding to CMX parameters. + + Returns + ------- + inds : np.ndarray + Array of CMX indices in model. + """ + inds = [int(p.split("_")[-1]) for p in self.params if "CMX_" in p] + return np.array(inds) + + def setup(self): + super().setup() + # Get CMX mapping. + # Register the CMX derivatives + for prefix_par in self.get_params_of_type("prefixParameter"): + if prefix_par.startswith("CMX_"): + self.register_deriv_funcs(self.d_delay_d_cmparam, prefix_par) + self.register_cm_deriv_funcs(self.d_cm_d_CMX, prefix_par) + + def validate(self): + """Validate the CMX parameters.""" + super().validate() + CMX_mapping = self.get_prefix_mapping_component("CMX_") + CMXR1_mapping = self.get_prefix_mapping_component("CMXR1_") + CMXR2_mapping = self.get_prefix_mapping_component("CMXR2_") + if CMX_mapping.keys() != CMXR1_mapping.keys(): + # FIXME: report mismatch + raise ValueError( + "CMX_ parameters do not " + "match CMXR1_ parameters. " + "Please check your prefixed parameters." + ) + if CMX_mapping.keys() != CMXR2_mapping.keys(): + raise ValueError( + "CMX_ parameters do not " + "match CMXR2_ parameters. " + "Please check your prefixed parameters." + ) + r1 = np.zeros(len(CMX_mapping)) + r2 = np.zeros(len(CMX_mapping)) + indices = np.zeros(len(CMX_mapping), dtype=np.int32) + for j, index in enumerate(CMX_mapping): + if ( + getattr(self, f"CMXR1_{index:04d}").quantity is not None + and getattr(self, f"CMXR2_{index:04d}").quantity is not None + ): + r1[j] = getattr(self, f"CMXR1_{index:04d}").quantity.mjd + r2[j] = getattr(self, f"CMXR2_{index:04d}").quantity.mjd + indices[j] = index + for j, index in enumerate(CMXR1_mapping): + if np.any((r1[j] > r1) & (r1[j] < r2)): + k = np.where((r1[j] > r1) & (r1[j] < r2))[0] + for kk in k.flatten(): + log.warning( + f"Start of CMX_{index:04d} ({r1[j]}-{r2[j]}) overlaps with CMX_{indices[kk]:04d} ({r1[kk]}-{r2[kk]})" + ) + if np.any((r2[j] > r1) & (r2[j] < r2)): + k = np.where((r2[j] > r1) & (r2[j] < r2))[0] + for kk in k.flatten(): + log.warning( + f"End of CMX_{index:04d} ({r1[j]}-{r2[j]}) overlaps with CMX_{indices[kk]:04d} ({r1[kk]}-{r2[kk]})" + ) + + def validate_toas(self, toas): + CMX_mapping = self.get_prefix_mapping_component("CMX_") + CMXR1_mapping = self.get_prefix_mapping_component("CMXR1_") + CMXR2_mapping = self.get_prefix_mapping_component("CMXR2_") + bad_parameters = [] + for k in CMXR1_mapping.keys(): + if self._parent[CMX_mapping[k]].frozen: + continue + b = self._parent[CMXR1_mapping[k]].quantity.mjd * u.d + e = self._parent[CMXR2_mapping[k]].quantity.mjd * u.d + mjds = toas.get_mjds() + n = np.sum((b <= mjds) & (mjds < e)) + if n == 0: + bad_parameters.append(CMX_mapping[k]) + if bad_parameters: + raise MissingTOAs(bad_parameters) + + def cmx_cm(self, toas): + condition = {} + tbl = toas.table + if not hasattr(self, "cmx_toas_selector"): + self.cmx_toas_selector = TOASelect(is_range=True) + CMX_mapping = self.get_prefix_mapping_component("CMX_") + CMXR1_mapping = self.get_prefix_mapping_component("CMXR1_") + CMXR2_mapping = self.get_prefix_mapping_component("CMXR2_") + for epoch_ind in CMX_mapping.keys(): + r1 = getattr(self, CMXR1_mapping[epoch_ind]).quantity + r2 = getattr(self, CMXR2_mapping[epoch_ind]).quantity + condition[CMX_mapping[epoch_ind]] = (r1.mjd, r2.mjd) + select_idx = self.cmx_toas_selector.get_select_index( + condition, tbl["mjd_float"] + ) + # Get CMX delays + cm = np.zeros(len(tbl)) * self._parent.CM.units + for k, v in select_idx.items(): + cm[v] += getattr(self, k).quantity + return cm + + def CMX_chromatic_delay(self, toas, acc_delay=None): + """This is a wrapper function for interacting with the TimingModel class""" + return self.chromatic_type_delay(toas) + + def d_cm_d_CMX(self, toas, param_name, acc_delay=None): + condition = {} + tbl = toas.table + if not hasattr(self, "cmx_toas_selector"): + self.cmx_toas_selector = TOASelect(is_range=True) + param = getattr(self, param_name) + cmx_index = param.index + CMXR1_mapping = self.get_prefix_mapping_component("CMXR1_") + CMXR2_mapping = self.get_prefix_mapping_component("CMXR2_") + r1 = getattr(self, CMXR1_mapping[cmx_index]).quantity + r2 = getattr(self, CMXR2_mapping[cmx_index]).quantity + condition = {param_name: (r1.mjd, r2.mjd)} + select_idx = self.cmx_toas_selector.get_select_index( + condition, tbl["mjd_float"] + ) + + cmx = np.zeros(len(tbl)) + for k, v in select_idx.items(): + cmx[v] = 1.0 + return cmx * (u.pc / u.cm**3) / (u.pc / u.cm**3) + + def print_par(self, format="pint"): + result = "" + CMX_mapping = self.get_prefix_mapping_component("CMX_") + CMXR1_mapping = self.get_prefix_mapping_component("CMXR1_") + CMXR2_mapping = self.get_prefix_mapping_component("CMXR2_") + sorted_list = sorted(CMX_mapping.keys()) + for ii in sorted_list: + result += getattr(self, CMX_mapping[ii]).as_parfile_line(format=format) + result += getattr(self, CMXR1_mapping[ii]).as_parfile_line(format=format) + result += getattr(self, CMXR2_mapping[ii]).as_parfile_line(format=format) + return result diff --git a/src/pint/models/cmwavex.py b/src/pint/models/cmwavex.py index 465700565..bd1ebca9a 100644 --- a/src/pint/models/cmwavex.py +++ b/src/pint/models/cmwavex.py @@ -6,7 +6,7 @@ from warnings import warn from pint.models.parameter import MJDParameter, prefixParameter -from pint.models.timing_model import MissingParameter +from pint.exceptions import MissingParameter from pint.models.chromatic_model import Chromatic, cmu from pint import DMconst diff --git a/src/pint/models/dispersion_model.py b/src/pint/models/dispersion_model.py index e5317b06c..595379717 100644 --- a/src/pint/models/dispersion_model.py +++ b/src/pint/models/dispersion_model.py @@ -14,14 +14,15 @@ prefixParameter, maskParameter, ) -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, ) -from pint import DMconst +from pint import DMconst, dmu +from pint.exceptions import MissingParameter, MissingTOAs # This value is cited from Duncan Lorimer, Michael Kramer, Handbook of Pulsar # Astronomy, Second edition, Page 86, Note 1 @@ -75,7 +76,7 @@ def dm_value(self, toas): DM values at given TOAs in the unit of DM. """ toas_table = toas if isinstance(toas, Table) else toas.table - dm = np.zeros(len(toas_table)) * self._parent.DM.units + dm = np.zeros(len(toas_table)) * dmu for dm_f in self.dm_value_funcs: dm += dm_f(toas) @@ -182,8 +183,8 @@ def __init__(self): def setup(self): super().setup() - base_dms = list(self.get_prefix_mapping_component("DM").values()) - base_dms += ["DM"] + + base_dms = ["DM"] + list(self.get_prefix_mapping_component("DM").values()) for dm_name in base_dms: self.register_deriv_funcs(self.d_delay_d_dmparam, dm_name) @@ -205,10 +206,10 @@ def validate(self): ) def DM_derivative_unit(self, n): - return "pc cm^-3/yr^%d" % n if n else "pc cm^-3" + return f"pc cm^-3/yr^{n}" if n != 0 else "pc cm^-3" def DM_derivative_description(self, n): - return "%d'th time derivative of the dispersion measure" % n + return f"{n}'th time derivative of the dispersion measure" def get_DM_terms(self): """Return a list of the DM term values in the model: [DM, DM1, ..., DMn]""" @@ -222,7 +223,7 @@ def base_dm(self, toas): if DMEPOCH is None: # Should be ruled out by validate() raise ValueError( - f"DMEPOCH not set but some derivatives are not zero: {dm_terms}" + f"DMEPOCH not set but some DM derivatives are not zero: {dm_terms}" ) else: dt = (toas["tdbld"] - DMEPOCH) * u.day @@ -308,7 +309,7 @@ def change_dmepoch(self, new_epoch): class DispersionDMX(Dispersion): - """This class provides a DMX model - multiple DM values. + """This class provides a DMX model - piecewise-constant DM variations. This model lets the user specify time ranges and fit for a different DM value in each time range. @@ -344,7 +345,7 @@ def __init__(self): self.delay_funcs_component += [self.DMX_dispersion_delay] def add_DMX_range(self, mjd_start, mjd_end, index=None, dmx=0, frozen=True): - """Add DMX range to a dispersion model with specified start/end MJDs and DMX. + """Add DMX range to a dispersion model with specified start/end MJDs and DMX value. Parameters ---------- @@ -626,6 +627,10 @@ def validate(self): r2[j] = getattr(self, f"DMXR2_{index:04d}").quantity.mjd indices[j] = index for j, index in enumerate(DMXR1_mapping): + if (r1[j] == r2[j]) and (r1[j] > 0): + log.warning( + f"Start of DMX_{index:04d} ({r1[j]}) equal to end of DMX_{index:04d} ({r2[j]})" + ) if np.any((r1[j] > r1) & (r1[j] < r2)): k = np.where((r1[j] > r1) & (r1[j] < r2))[0] for kk in k.flatten(): @@ -638,6 +643,8 @@ def validate(self): log.warning( f"End of DMX_{index:04d} ({r1[j]}-{r2[j]}) overlaps with DMX_{indices[kk]:04d} ({r1[kk]}-{r2[kk]})" ) + if not hasattr(self, "dmx_toas_selector"): + self.dmx_toas_selector = TOASelect(is_range=True) def validate_toas(self, toas): DMX_mapping = self.get_prefix_mapping_component("DMX_") @@ -650,7 +657,7 @@ def validate_toas(self, toas): b = self._parent[DMXR1_mapping[k]].quantity.mjd * u.d e = self._parent[DMXR2_mapping[k]].quantity.mjd * u.d mjds = toas.get_mjds() - n = np.sum((b <= mjds) & (mjds < e)) + n = np.sum((b <= mjds) & (mjds <= e)) if n == 0: bad_parameters.append(DMX_mapping[k]) if bad_parameters: @@ -672,7 +679,7 @@ def dmx_dm(self, toas): condition, tbl["mjd_float"] ) # Get DMX delays - dm = np.zeros(len(tbl)) * self._parent.DM.units + dm = np.zeros(len(tbl)) * dmu for k, v in select_idx.items(): dm[v] += getattr(self, k).quantity return dm diff --git a/src/pint/models/dmwavex.py b/src/pint/models/dmwavex.py index 368862391..9d3fa614f 100644 --- a/src/pint/models/dmwavex.py +++ b/src/pint/models/dmwavex.py @@ -6,7 +6,7 @@ from warnings import warn from pint.models.parameter import MJDParameter, prefixParameter -from pint.models.timing_model import MissingParameter +from pint.exceptions import MissingParameter from pint.models.dispersion_model import Dispersion from pint import DMconst, dmu diff --git a/src/pint/models/frequency_dependent.py b/src/pint/models/frequency_dependent.py index 2215b5c08..a5d37b7a8 100644 --- a/src/pint/models/frequency_dependent.py +++ b/src/pint/models/frequency_dependent.py @@ -6,7 +6,8 @@ import numpy as np from pint.models.parameter import prefixParameter -from pint.models.timing_model import DelayComponent, MissingParameter +from pint.models.timing_model import DelayComponent +from pint.exceptions import MissingParameter class FD(DelayComponent): diff --git a/src/pint/models/glitch.py b/src/pint/models/glitch.py index f198a27c2..33d90e9f1 100644 --- a/src/pint/models/glitch.py +++ b/src/pint/models/glitch.py @@ -6,8 +6,9 @@ from loguru import logger as log from pint.models.parameter import prefixParameter -from pint.models.timing_model import MissingParameter, PhaseComponent +from pint.models.timing_model import PhaseComponent from pint.utils import split_prefixed_name +from pint.exceptions import MissingParameter class Glitch(PhaseComponent): @@ -144,49 +145,48 @@ def setup(self): ] for idx in set(self.glitch_indices): for param in self.glitch_prop: - if not hasattr(self, param + "%d" % idx): + check = f"{param}{idx}" + if not hasattr(self, check): param0 = getattr(self, f"{param}1") self.add_param(param0.new_param(idx)) - getattr(self, param + "%d" % idx).value = 0.0 + getattr(self, check).value = 0.0 self.register_deriv_funcs( - getattr(self, f"d_phase_d_{param[:-1]}"), param + "%d" % idx + getattr(self, f"d_phase_d_{param[:-1]}"), check ) def validate(self): """Validate parameters input.""" super().validate() for idx in set(self.glitch_indices): - if not hasattr(self, "GLEP_%d" % idx): - msg = "Glitch Epoch is needed for Glitch %d." % idx - raise MissingParameter("Glitch", "GLEP_%d" % idx, msg) - else: # Check to see if both the epoch and phase are to be fit - if ( - hasattr(self, "GLPH_%d" % idx) - and (not getattr(self, "GLEP_%d" % idx).frozen) - and (not getattr(self, "GLPH_%d" % idx).frozen) - ): - raise ValueError( - "Both the glitch epoch and phase cannot be fit for Glitch %d." - % idx - ) + glep = f"GLEP_{idx}" + glph = f"GLPH_{idx}" + if (not hasattr(self, glep)) or (getattr(self, glep).quantity is None): + msg = f"Glitch Epoch is needed for Glitch {idx}" + raise MissingParameter("Glitch", glep, msg) + # Check to see if both the epoch and phase are to be fit + if ( + hasattr(self, glph) + and (not getattr(self, glep).frozen) + and (not getattr(self, glph).frozen) + ): + raise ValueError( + f"Both the glitch epoch and phase cannot be fit for Glitch {idx}." + ) # Check the Decay Term. glf0dparams = [x for x in self.params if x.startswith("GLF0D_")] for glf0dnm in glf0dparams: glf0d = getattr(self, glf0dnm) idx = glf0d.index - if glf0d.value != 0.0 and getattr(self, "GLTD_%d" % idx).value == 0.0: - msg = ( - "None zero GLF0D_%d parameter needs a none" - " zero GLTD_%d parameter" % (idx, idx) - ) - raise MissingParameter("Glitch", "GLTD_%d" % idx, msg) + if glf0d.value != 0.0 and getattr(self, f"GLTD_{idx}").value == 0.0: + msg = f"Non-zero GLF0D_{idx} parameter needs a non-zero GLTD_{idx} parameter" + raise MissingParameter("Glitch", f"GLTD_{idx}", msg) def print_par(self, format="pint"): result = "" for idx in set(self.glitch_indices): for param in self.glitch_prop: - par = getattr(self, param + "%d" % idx) + par = getattr(self, f"{param}{idx}") result += par.as_parfile_line(format=format) return result @@ -203,22 +203,21 @@ def glitch_phase(self, toas, delay): glep = getattr(self, glepnm) idx = glep.index eph = glep.value - dphs = getattr(self, "GLPH_%d" % idx).quantity - dF0 = getattr(self, "GLF0_%d" % idx).quantity - dF1 = getattr(self, "GLF1_%d" % idx).quantity - dF2 = getattr(self, "GLF2_%d" % idx).quantity + dphs = getattr(self, f"GLPH_{idx}").quantity + dF0 = getattr(self, f"GLF0_{idx}").quantity + dF1 = getattr(self, f"GLF1_{idx}").quantity + dF2 = getattr(self, f"GLF2_{idx}").quantity dt = (tbl["tdbld"] - eph) * u.day - delay dt = dt.to(u.second) affected = dt > 0.0 # TOAs affected by glitch # decay term - dF0D = getattr(self, "GLF0D_%d" % idx).quantity + dF0D = getattr(self, f"GLF0D_{idx}").quantity if dF0D != 0.0: - tau = getattr(self, "GLTD_%d" % idx).quantity + tau = getattr(self, f"GLTD_{idx}").quantity decayterm = dF0D * tau * (1.0 - np.exp(-(dt[affected] / tau))) else: decayterm = u.Quantity(0) - log.debug(f"Glitch phase for glitch {idx}: {dphs} {dphs.unit}") phs[affected] += ( dphs + dt[affected] @@ -231,114 +230,91 @@ def glitch_phase(self, toas, delay): ) return phs - def deriv_prep(self, toas, param, delay): + def deriv_prep(self, toas, param, delay, check_param): """Get the things we need for any of the derivative calcs""" - tbl = toas.table p, ids, idv = split_prefixed_name(param) + if p != f"{check_param}_": + raise ValueError( + f"Can not calculate d_phase_d_{check_param} with respect to {param}." + ) + tbl = toas.table eph = getattr(self, f"GLEP_{ids}").value dt = (tbl["tdbld"] - eph) * u.day - delay dt = dt.to(u.second) affected = np.where(dt > 0.0)[0] - return tbl, p, ids, idv, dt, affected + par = getattr(self, param) + zeros = np.zeros(len(tbl), dtype=np.longdouble) << 1 / par.units + return tbl, p, ids, idv, dt, affected, par, zeros def d_phase_d_GLPH(self, toas, param, delay): """Calculate the derivative wrt GLPH""" - tbl, p, ids, idv, dt, affected = self.deriv_prep(toas, param, delay) - if p != "GLPH_": - raise ValueError( - f"Can not calculate d_phase_d_GLPH with respect to {param}." - ) - par_GLPH = getattr(self, param) - dpdGLPH = np.zeros(len(tbl), dtype=np.longdouble) / par_GLPH.units - dpdGLPH[affected] += 1.0 / par_GLPH.units + tbl, p, ids, idv, dt, affected, par_GLPH, dpdGLPH = self.deriv_prep( + toas, param, delay, "GLPH" + ) + dpdGLPH[affected] = 1.0 / par_GLPH.units return dpdGLPH def d_phase_d_GLF0(self, toas, param, delay): """Calculate the derivative wrt GLF0""" - tbl, p, ids, idv, dt, affected = self.deriv_prep(toas, param, delay) - if p != "GLF0_": - raise ValueError( - f"Can not calculate d_phase_d_GLF0 with respect to {param}." - ) - par_GLF0 = getattr(self, param) - dpdGLF0 = np.zeros(len(tbl), dtype=np.longdouble) / par_GLF0.units + tbl, p, ids, idv, dt, affected, par_GLF0, dpdGLF0 = self.deriv_prep( + toas, param, delay, "GLF0" + ) dpdGLF0[affected] = dt[affected] return dpdGLF0 def d_phase_d_GLF1(self, toas, param, delay): """Calculate the derivative wrt GLF1""" - tbl, p, ids, idv, dt, affected = self.deriv_prep(toas, param, delay) - if p != "GLF1_": - raise ValueError( - f"Can not calculate d_phase_d_GLF1 with respect to {param}." - ) - par_GLF1 = getattr(self, param) - dpdGLF1 = np.zeros(len(tbl), dtype=np.longdouble) / par_GLF1.units - dpdGLF1[affected] += np.longdouble(0.5) * dt[affected] * dt[affected] + tbl, p, ids, idv, dt, affected, par_GLF1, dpdGLF1 = self.deriv_prep( + toas, param, delay, "GLF1" + ) + dpdGLF1[affected] = 0.5 * dt[affected] ** 2 return dpdGLF1 def d_phase_d_GLF2(self, toas, param, delay): """Calculate the derivative wrt GLF1""" - tbl, p, ids, idv, dt, affected = self.deriv_prep(toas, param, delay) - if p != "GLF2_": - raise ValueError( - f"Can not calculate d_phase_d_GLF2 with respect to {param}." - ) - par_GLF2 = getattr(self, param) - dpdGLF2 = np.zeros(len(tbl), dtype=np.longdouble) / par_GLF2.units - dpdGLF2[affected] += ( - np.longdouble(1.0) / 6.0 * dt[affected] * dt[affected] * dt[affected] + tbl, p, ids, idv, dt, affected, par_GLF2, dpdGLF2 = self.deriv_prep( + toas, param, delay, "GLF2" ) + dpdGLF2[affected] = (1.0 / 6.0) * dt[affected] ** 3 return dpdGLF2 def d_phase_d_GLF0D(self, toas, param, delay): """Calculate the derivative wrt GLF0D""" - tbl, p, ids, idv, dt, affected = self.deriv_prep(toas, param, delay) - if p != "GLF0D_": - raise ValueError( - f"Can not calculate d_phase_d_GLF0D with respect to {param}." - ) - par_GLF0D = getattr(self, param) - tau = getattr(self, "GLTD_%d" % idv).quantity - dpdGLF0D = np.zeros(len(tbl), dtype=np.longdouble) / par_GLF0D.units - dpdGLF0D[affected] += tau * (np.longdouble(1.0) - np.exp(-dt[affected] / tau)) + tbl, p, ids, idv, dt, affected, par_GLF0D, dpdGLF0D = self.deriv_prep( + toas, param, delay, "GLF0D" + ) + tau = getattr(self, f"GLTD_{ids}").quantity + dpdGLF0D[affected] = tau * (1.0 - np.exp(-dt[affected] / tau)) return dpdGLF0D def d_phase_d_GLTD(self, toas, param, delay): """Calculate the derivative wrt GLTD""" - tbl, p, ids, idv, dt, affected = self.deriv_prep(toas, param, delay) - if p != "GLTD_": - raise ValueError( - f"Can not calculate d_phase_d_GLTD with respect to {param}." - ) - par_GLTD = getattr(self, param) + tbl, p, ids, idv, dt, affected, par_GLTD, dpdGLTD = self.deriv_prep( + toas, param, delay, "GLTD" + ) if par_GLTD.value == 0.0: - return np.zeros(len(tbl), dtype=np.longdouble) / par_GLTD.units + return dpdGLTD glf0d = getattr(self, f"GLF0D_{ids}").quantity tau = par_GLTD.quantity - dpdGLTD = np.zeros(len(tbl), dtype=np.longdouble) / par_GLTD.units - dpdGLTD[affected] += glf0d * ( - np.longdouble(1.0) - np.exp(-dt[affected] / tau) - ) + glf0d * tau * (-np.exp(-dt[affected] / tau)) * dt[affected] / (tau * tau) + et = np.exp(-dt[affected] / tau) + dpdGLTD[affected] = glf0d * ( + 1.0 - np.exp(-dt[affected] / tau) * (1.0 + dt[affected] / tau) + ) return dpdGLTD def d_phase_d_GLEP(self, toas, param, delay): """Calculate the derivative wrt GLEP""" - tbl, p, ids, idv, dt, affected = self.deriv_prep(toas, param, delay) - if p != "GLEP_": - raise ValueError( - f"Can not calculate d_phase_d_GLEP with respect to {param}." - ) - par_GLEP = getattr(self, param) + tbl, p, ids, idv, dt, affected, par_GLEP, dpdGLEP = self.deriv_prep( + toas, param, delay, "GLEP" + ) glf0 = getattr(self, f"GLF0_{ids}").quantity glf1 = getattr(self, f"GLF1_{ids}").quantity glf2 = getattr(self, f"GLF2_{ids}").quantity glf0d = getattr(self, f"GLF0D_{ids}").quantity tau = getattr(self, f"GLTD_{ids}").quantity - dpdGLEP = np.zeros(len(tbl), dtype=np.longdouble) / par_GLEP.units dpdGLEP[affected] += ( -glf0 + -glf1 * dt[affected] + -0.5 * glf2 * dt[affected] ** 2 ) if tau.value != 0.0: - dpdGLEP[affected] += -glf0d / np.exp(dt[affected] / tau) + dpdGLEP[affected] -= glf0d * np.exp(-dt[affected] / tau) return dpdGLEP diff --git a/src/pint/models/ifunc.py b/src/pint/models/ifunc.py index 323380ac3..465dd2c87 100644 --- a/src/pint/models/ifunc.py +++ b/src/pint/models/ifunc.py @@ -4,7 +4,8 @@ import numpy as np from pint.models.parameter import floatParameter, prefixParameter -from pint.models.timing_model import PhaseComponent, MissingParameter +from pint.models.timing_model import PhaseComponent +from pint.exceptions import MissingParameter class IFunc(PhaseComponent): diff --git a/src/pint/models/model_builder.py b/src/pint/models/model_builder.py index 088989967..36a7ac5dc 100644 --- a/src/pint/models/model_builder.py +++ b/src/pint/models/model_builder.py @@ -15,16 +15,11 @@ AllComponents, TimingModel, ignore_prefix, - UnknownBinaryModel, - UnknownParameter, - TimingModelError, - MissingBinaryError, ignore_params, ignore_prefix, ) from pint.toa import TOAs, get_TOAs from pint.utils import ( - PrefixError, interesting_lines, lines_of, split_prefixed_name, @@ -33,6 +28,14 @@ from pint.models.tcb_conversion import convert_tcb_tdb from pint.models.binary_ddk import _convert_kin, _convert_kom from pint.types import file_like, quantity_like +from pint.exceptions import ( + PrefixError, + ComponentConflict, + UnknownParameter, + UnknownBinaryModel, + TimingModelError, + MissingBinaryError, +) __all__ = ["ModelBuilder", "get_model", "get_model_and_toas"] @@ -52,10 +55,6 @@ ] -class ComponentConflict(ValueError): - """Error for multiple components can be select but no other indications.""" - - def parse_parfile(parfile): """Function for parsing .par file or .par style StringIO. @@ -252,6 +251,10 @@ def __call__( if not hasattr(tm, "NoiseComponent_list"): setattr(tm, "NoiseComponent_list", []) + tm.meta["allow_tcb"] = allow_tcb_ + tm.meta["convert_tcb"] = convert_tcb + tm.meta["allow_T2"] = allow_T2 + return tm def _validate_components(self): @@ -852,6 +855,7 @@ def get_model( **kwargs, ) model.name = parfile + model.meta["original_name"] = parfile return model diff --git a/src/pint/models/noise_model.py b/src/pint/models/noise_model.py index 7fa364bb9..f5491167a 100644 --- a/src/pint/models/noise_model.py +++ b/src/pint/models/noise_model.py @@ -122,6 +122,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) @@ -163,14 +164,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") @@ -441,8 +442,7 @@ def ecorr_cov_matrix(self, toas): class PLDMNoise(NoiseComponent): - """Model of DM variations as radio frequency-dependent noise with a - power-law spectrum. + """Model of DM variations as radio frequency-dependent noise with a power-law spectrum. Variations in DM over time result from both the proper motion of the pulsar and the changing electron number density along the line of sight @@ -458,9 +458,11 @@ class PLDMNoise(NoiseComponent): .. paramtable:: :class: pint.models.noise_model.PLDMNoise - Note - ---- - Ref: Lentati et al. 2014, MNRAS 437(3), 3004-3023 + References + ---------- + - Lentati et al. 2014, MNRAS 437(3), 3004-3023 [1]_ + + .. [1] https://ui.adsabs.harvard.edu/abs/2014MNRAS.437.3004L/abstract """ @@ -594,8 +596,7 @@ def pl_dm_cov_matrix(self, toas): class PLChromNoise(NoiseComponent): - """Model of a radio frequency-dependent noise with a power-law spectrum and - arbitrary chromatic index. + """Model of a radio frequency-dependent noise with a power-law spectrum and arbitrary chromatic index. Such variations are usually attributed to time-variable scattering in the ISM. Scattering smears/broadens the shape of the pulse profile by convolving it with @@ -615,9 +616,11 @@ class PLChromNoise(NoiseComponent): .. paramtable:: :class: pint.models.noise_model.PLChromNoise - Note - ---- - Ref: Lentati et al. 2014, MNRAS 437(3), 3004-3023 + References + ---------- + - Lentati et al. 2014, MNRAS 437(3), 3004-3023 [1]_ + + .. [1] https://ui.adsabs.harvard.edu/abs/2014MNRAS.437.3004L/abstract """ register = True @@ -728,10 +731,11 @@ class PLRedNoise(NoiseComponent): .. paramtable:: :class: pint.models.noise_model.PLRedNoise - Note - ---- - Ref: Lentati et al. 2014, MNRAS 437(3), 3004-3023 + References + ---------- + - Lentati et al. 2014, MNRAS 437(3), 3004-3023 [1]_ + .. [1] https://ui.adsabs.harvard.edu/abs/2014MNRAS.437.3004L/abstract """ register = True diff --git a/src/pint/models/parameter.py b/src/pint/models/parameter.py index bb7ad588c..5e3d9e607 100644 --- a/src/pint/models/parameter.py +++ b/src/pint/models/parameter.py @@ -48,7 +48,7 @@ ) from pint.toa_select import TOASelect from pint.utils import split_prefixed_name - +from pint.exceptions import InvalidModelParameters # potential parfile formats # in one place for consistency @@ -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__") diff --git a/src/pint/models/piecewise.py b/src/pint/models/piecewise.py index 9b9fdd8b3..93c05487a 100644 --- a/src/pint/models/piecewise.py +++ b/src/pint/models/piecewise.py @@ -4,8 +4,9 @@ import numpy as np from pint.models.parameter import prefixParameter -from pint.models.timing_model import MissingParameter, PhaseComponent +from pint.models.timing_model import PhaseComponent from pint.utils import split_prefixed_name, taylor_horner +from pint.exceptions import MissingParameter class PiecewiseSpindown(PhaseComponent): diff --git a/src/pint/models/pulsar_binary.py b/src/pint/models/pulsar_binary.py index 5f81e29ac..26d77bda2 100644 --- a/src/pint/models/pulsar_binary.py +++ b/src/pint/models/pulsar_binary.py @@ -21,12 +21,10 @@ from pint.models.stand_alone_psr_binaries import binary_orbits as bo from pint.models.timing_model import ( DelayComponent, - MissingParameter, - TimingModelError, - UnknownParameter, ) from pint.utils import taylor_horner_deriv, parse_time from pint.pulsar_ecliptic import PulsarEcliptic +from pint.exceptions import MissingParameter, TimingModelError, UnknownParameter # def _p_to_f(p): @@ -63,6 +61,18 @@ class PulsarBinary(DelayComponent): - FB0 - orbital frequency (1/s, alternative to PB, non-negative) - FBn - time derivatives of orbital frequency (1/s**(n+1)) + The following ORBWAVEs parameters define a Fourier series model for orbital phase + variations, as an alternative to the FBn Taylor series expansion: + + - ORBWAVE_OM - base angular frequency for ORBWAVEs expansion (rad / s) + - ORBWAVE_EPOCH - reference epoch for ORBWAVEs model (MJD) + - ORBWAVECn/ORBWAVESn - coefficients for cosine/sine components (dimensionless) + + The orbital phase is then given by: + orbits(t) = (t - T0) / PB + + \sum_{n=0} (ORBWAVECn cos(ORBWAVE_OM * (n + 1) * (t - ORBWAVE_EPOCH)) + + ORBWAVESn sin(ORBWAVE_OM * (n + 1) * (t - ORBWAVE_EPOCH)) + The internal calculation code uses different names for some parameters: - Eccentric Anomaly: E (not parameter ECC) @@ -203,6 +213,50 @@ def __init__(self): tcb2tdb_scale_factor=u.Quantity(1), ) ) + + self.add_param( + prefixParameter( + name="ORBWAVEC0", + value=None, + units="", + description="Amplitude of cosine wave in Tasc-shift function", + aliases=["ORBWAVEC"], + description_template=self.ORBWAVEC_description, + type_match="float", + tcb2tdb_scale_factor=u.Quantity(1), + ) + ) + + self.add_param( + prefixParameter( + name="ORBWAVES0", + value=None, + units="", + description="Amplitude of sine wave in Tasc-shift function", + aliases=["ORBWAVES"], + description_template=self.ORBWAVES_description, + type_match="float", + tcb2tdb_scale_factor=u.Quantity(1), + ) + ) + + self.add_param( + floatParameter( + name="ORBWAVE_OM", + units="rad/s", + description="Base frequency for ORBWAVEs model", + tcb2tdb_scale_factor=u.Quantity(1), + ) + ) + self.add_param( + MJDParameter( + name="ORBWAVE_EPOCH", + description="Reference epoch for ORBWAVEs model", + time_scale="tdb", + tcb2tdb_scale_factor=u.Quantity(1), + ) + ) + self.internal_params = [] self.warn_default_params = ["ECC", "OM"] # Set up delay function @@ -253,6 +307,48 @@ def setup(self): # func=_pdot_to_fdot, # ) # ) + + ORBWAVES_mapping = self.get_prefix_mapping_component("ORBWAVES") + ORBWAVES = { + ows: getattr(self, ows).quantity for ows in ORBWAVES_mapping.values() + } + ORBWAVEC_mapping = self.get_prefix_mapping_component("ORBWAVEC") + ORBWAVEC = { + owc: getattr(self, owc).quantity for owc in ORBWAVEC_mapping.values() + } + + if any(v is not None for v in ORBWAVES.values()): + for k in ["ORBWAVE_OM", "ORBWAVE_EPOCH"]: + self.binary_instance.add_binary_params(k, getattr(self, k).value) + + for k in ORBWAVES.keys(): + self.binary_instance.add_binary_params(k, ORBWAVES[k]) + + for k in ORBWAVEC.keys(): + self.binary_instance.add_binary_params(k, ORBWAVEC[k]) + + using_FBX = any(v is not None for v in FBXs.values()) + if using_FBX: + fbx = sorted(list(FBXs.keys())) + if len(fbx) > 2: + raise ValueError("Only FB0/FB1 are supported.") + if (len(fbx) == 2) and (fbx[1] != "FB1"): + raise ValueError("Only FB0/FB1 are supported.") + self.binary_instance.orbits_cls = bo.OrbitWavesFBX( + self.binary_instance, + fbx + + ["TASC", "ORBWAVE_OM", "ORBWAVE_EPOCH"] + + list(ORBWAVES.keys()) + + list(ORBWAVEC.keys()), + ) + else: + self.binary_instance.orbits_cls = bo.OrbitWaves( + self.binary_instance, + ["PB", "TASC", "ORBWAVE_OM", "ORBWAVE_EPOCH"] + + list(ORBWAVES.keys()) + + list(ORBWAVEC.keys()), + ) + # Note: if we are happy to use these to show alternate parameterizations then this can be uncommented # else: # # remove the FB parameterization, replace with functions @@ -489,6 +585,18 @@ def FBX_unit(self, n): def FBX_description(self, n): return "%dth time derivative of frequency of orbit" % n + def ORBWAVES_description(self, n): + return ( + "Coefficient of the %dth sine wave in Fourier series model of Tasc variations" + % n + ) + + def ORBWAVEC_description(self, n): + return ( + "Coefficient of the %dth cosine wave in Fourier series model of Tasc variations" + % n + ) + def change_binary_epoch(self, new_epoch): """Change the epoch for this binary model. diff --git a/src/pint/models/solar_wind_dispersion.py b/src/pint/models/solar_wind_dispersion.py index dddaeb0f2..8a282fee3 100644 --- a/src/pint/models/solar_wind_dispersion.py +++ b/src/pint/models/solar_wind_dispersion.py @@ -10,7 +10,12 @@ from pint import DMconst from pint.models.dispersion_model import Dispersion -from pint.models.parameter import floatParameter, intParameter, prefixParameter +from pint.models.parameter import ( + MJDParameter, + floatParameter, + intParameter, + prefixParameter, +) import pint.utils from pint.models.timing_model import MissingTOAs from pint.toa_select import TOASelect @@ -302,6 +307,26 @@ def __init__(self): tcb2tdb_scale_factor=(const.c * DMconst), ) ) + self.add_param( + prefixParameter( + name="NE_SW1", + units="cm^-3 / yr", + value=0.0, + description="Derivative of Solar Wind density at 1 AU", + unit_template=self.NE_SW_derivative_unit, + description_template=self.NE_SW_derivative_description, + type_match="float", + tcb2tdb_scale_factor=(const.c * DMconst), + ) + ) + self.add_param( + MJDParameter( + name="SWEPOCH", + description="Epoch of NE_SW measurement", + time_scale="tdb", + tcb2tdb_scale_factor=u.Quantity(1), + ) + ) self.add_param( floatParameter( name="SWP", @@ -324,11 +349,24 @@ def __init__(self): def setup(self): super().setup() - self.register_dm_deriv_funcs(self.d_dm_d_ne_sw, "NE_SW") - self.register_deriv_funcs(self.d_delay_d_ne_sw, "NE_SW") + self.register_dm_deriv_funcs(self.d_dm_d_swp, "SWP") self.register_deriv_funcs(self.d_delay_d_swp, "SWP") + base_nesws = ["NE_SW"] + list( + self.get_prefix_mapping_component("NE_SW").values() + ) + + for nesw_name in base_nesws: + self.register_deriv_funcs(self.d_delay_d_ne_sw, nesw_name) + self.register_dm_deriv_funcs(self.d_dm_d_ne_sw, nesw_name) + + def NE_SW_derivative_unit(self, n): + return f"cm^-3/yr^{n}" if n != 0 else "cm^-3" + + def NE_SW_derivative_description(self, n): + return f"{n}'th time derivative of the Solar Wind density at 1 AU" + def solar_wind_geometry(self, toas): """Return the geometry of solar wind dispersion. @@ -375,22 +413,51 @@ def solar_wind_dm(self, toas): SWM==1: Hazboun et al. 2022 """ - if self.NE_SW.value == 0: + ne_sw_terms = self.get_NE_SW_terms() + + if len(ne_sw_terms) == 1: + ne_sw = self.NE_SW.quantity * np.ones(len(toas)) + else: + if any(t.value != 0 for t in ne_sw_terms[1:]): + SWEPOCH = self.SWEPOCH.value + if SWEPOCH is None: + # Should be ruled out by validate() + raise ValueError( + f"SWEPOCH not set but some NE_SW derivatives are not zero: {ne_sw_terms}" + ) + else: + dt = (toas["tdbld"] - SWEPOCH) * u.day + dt_value = dt.to_value(u.yr) + else: + dt_value = np.zeros(len(toas), dtype=np.longdouble) + + ne_sw_terms_value = [d.value for d in ne_sw_terms] + ne_sw = ( + pint.utils.taylor_horner(dt_value, ne_sw_terms_value) * self.NE_SW.units + ) + + if np.all(ne_sw.value == 0): return np.zeros(len(toas)) * u.pc / u.cm**3 + if self.SWM.value not in [0, 1]: raise NotImplementedError( f"Solar Dispersion Delay not implemented for SWM {self.SWM.value}" ) + solar_wind_geometry = self.solar_wind_geometry(toas) - solar_wind_dm = self.NE_SW.quantity * solar_wind_geometry + solar_wind_dm = ne_sw * solar_wind_geometry return solar_wind_dm.to(u.pc / u.cm**3) def solar_wind_delay(self, toas, acc_delay=None): """This is a wrapper function to compute solar wind dispersion delay.""" - if self.NE_SW.value == 0: - return np.zeros(len(toas)) * u.s return self.dispersion_type_delay(toas) + def get_NE_SW_terms(self): + """Return a list of the NE_SW term values in the model: [NE_SW, NE_SW1, ..., NE_SWn]""" + return [self.NE_SW.quantity] + self._parent.get_prefix_list( + "NE_SW", start_index=1 + ) + def d_dm_d_ne_sw(self, toas, param_name, acc_delay=None): """Derivative of of DM wrt the solar wind ne amplitude.""" if self.SWM.value in [0, 1]: @@ -399,17 +466,31 @@ def d_dm_d_ne_sw(self, toas, param_name, acc_delay=None): raise NotImplementedError( "Solar Dispersion Delay not implemented for SWM %d" % self.SWM.value ) - return solar_wind_geometry + + par = getattr(self, param_name) + order = 0 if param_name == "NE_SW" else par.index + + ne_sws = self.get_NE_SW_terms() + + ne_sw_terms = np.longdouble(np.zeros(len(ne_sws))) + ne_sw_terms[order] = np.longdouble(1.0) + + if self.SWEPOCH.value is None: + if any(t.value != 0 for t in ne_sws[1:]): + # Should be ruled out by validate() + raise ValueError(f"SWEPOCH is not set but {param_name} is not zero") + SWEPOCH = 0 + else: + SWEPOCH = self.SWEPOCH.value + + dt = (toas["tdbld"] - SWEPOCH) * u.day + + return (solar_wind_geometry * dt**order / scipy.special.factorial(order)).to( + u.pc / u.cm**3 / par.units + ) def d_delay_d_ne_sw(self, toas, param_name, acc_delay=None): - try: - bfreq = self._parent.barycentric_radio_freq(toas) - except AttributeError: - warn("Using topocentric frequency for dedispersion!") - bfreq = toas.table["freq"] - deriv = self.d_delay_d_dmparam(toas, "NE_SW") - deriv[bfreq < 1.0 * u.MHz] = 0.0 - return deriv + return self.d_delay_d_dmparam(toas, param_name) def d_solar_wind_geometry_d_swp(self, toas, param_name, acc_delay=None): """Derivative of solar_wind_geometry (path length) wrt power-law index p @@ -452,10 +533,12 @@ def d_delay_d_swp(self, toas, param_name, acc_delay=None): def print_par(self, format="pint"): result = "" - result += getattr(self, "NE_SW").as_parfile_line(format=format) - result += getattr(self, "SWM").as_parfile_line(format=format) - if self.SWM.value == 1: - result += getattr(self, "SWP").as_parfile_line(format=format) + for par in self.params: + if par == "SWP" and self.SWM.value != 1: + continue + else: + result += getattr(self, par).as_parfile_line(format=format) + return result def get_max_dm(self): diff --git a/src/pint/models/spindown.py b/src/pint/models/spindown.py index 8a8bce411..0bc42f996 100644 --- a/src/pint/models/spindown.py +++ b/src/pint/models/spindown.py @@ -6,9 +6,10 @@ import numpy from pint.models.parameter import MJDParameter, prefixParameter -from pint.models.timing_model import MissingParameter, PhaseComponent +from pint.models.timing_model import PhaseComponent from pint.pulsar_mjd import Time from pint.utils import split_prefixed_name, taylor_horner, taylor_horner_deriv +from pint.exceptions import MissingParameter class SpindownBase(PhaseComponent): 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 6e925fd58..828e49610 100644 --- a/src/pint/models/stand_alone_psr_binaries/DD_model.py +++ b/src/pint/models/stand_alone_psr_binaries/DD_model.py @@ -6,6 +6,7 @@ from loguru import logger as log from pint import Tsun +from pint.models.parameter import InvalidModelParameters from .binary_generic import PSR_BINARY @@ -709,6 +710,9 @@ def delayS(self): cOmega = np.cos(self.omega()) TM2 = self.M2.value * Tsun + if np.any(self.SINI < 0) or np.any(self.SINI > 1): + raise InvalidModelParameters("SINI parameter must be between 0 and 1") + return ( -2 * TM2 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 d4b9cd8a0..85ec55b47 100644 --- a/src/pint/models/stand_alone_psr_binaries/ELL1_model.py +++ b/src/pint/models/stand_alone_psr_binaries/ELL1_model.py @@ -4,6 +4,7 @@ import astropy.units as u import numpy as np +from pint.models.parameter import InvalidModelParameters from .binary_generic import PSR_BINARY @@ -600,6 +601,8 @@ def delayS(self): """ELL1 Shapiro delay. Lange et al 2001 eq. A16""" TM2 = self.TM2() Phi = self.Phi() + if np.any(self.SINI < 0) or np.any(self.SINI > 1): + raise InvalidModelParameters("SINI must be between 0 and 1") return -2 * TM2 * np.log(1 - self.SINI * np.sin(Phi)) def d_delayS_d_par(self, par): diff --git a/src/pint/models/stand_alone_psr_binaries/binary_orbits.py b/src/pint/models/stand_alone_psr_binaries/binary_orbits.py index 079607c05..2b3987d6a 100644 --- a/src/pint/models/stand_alone_psr_binaries/binary_orbits.py +++ b/src/pint/models/stand_alone_psr_binaries/binary_orbits.py @@ -237,3 +237,280 @@ def d_pbprime_d_par(self, par): if re.match(r"FB\d+", par) is not None else np.zeros(len(self.tt0)) * u.second / par_obj.unit ) + + +class OrbitWaves(Orbit): + """Orbit with orbital phase variations described by a Fourier series""" + + def __init__( + self, + parent, + orbit_params=[ + "PB", + "TASC", + "ORBWAVE_OM", + "ORBWAVE_EPOCH", + "ORBWAVEC0", + "ORBWAVES0", + ], + ): + # Generalize this to allow for instantiation with OrbitWavesFBX + label = self.__class__.__name__ + label = label[0].lower() + label[1:] + super().__init__(label, parent, orbit_params) + + Cindices = set() + Sindices = set() + + nc = 0 + ns = 0 + + for k in self.binary_params: + if re.match(r"ORBWAVEC\d+", k) is not None: + if k not in self.orbit_params: + self.orbit_params += [k] + Cindices.add(int(k[8:])) + nc += 1 + if re.match(r"ORBWAVES\d+", k) is not None: + if k not in self.orbit_params: + self.orbit_params += [k] + Sindices.add(int(k[8:])) + ns += 1 + + if Cindices != set(range(len(Cindices))) or Sindices != set( + range(len(Sindices)) + ): + raise ValueError( + f"Orbwave indices must be 0 up to some number k without gaps " + f"but are {Cindices} and {Sindices}." + ) + + if nc != ns: + raise ValueError( + f"Must have equal number of sine/cosine ORBWAVE pairs " + f"but have {ns} and {nc}." + ) + + if self.ORBWAVEC0 is None or self.ORBWAVES0 is None: + raise ValueError("The ORBWAVEC0 or ORBWAVES0 parameter is unset") + + self.nwaves = ns + + def _ORBWAVEs(self): + ORBWAVEs = np.zeros((self.nwaves, 2)) * u.Unit("") + + for ii in range(self.nwaves): + ORBWAVEs[ii, 0] = getattr(self, f"ORBWAVEC{ii}") + ORBWAVEs[ii, 1] = getattr(self, f"ORBWAVES{ii}") + + return ORBWAVEs + + def _tw(self): + tw = self.t - self.ORBWAVE_EPOCH.value * u.d + return tw + + def _deltaPhi(self): + tw = self._tw() + waveamps = self._ORBWAVEs() + OM = self.ORBWAVE_OM.to("radian/second") + + nh = np.arange(self.nwaves) + 1 + Cwaves = waveamps[:, 0, None] * np.cos(OM * nh[:, None] * tw[None, :]) + Swaves = waveamps[:, 1, None] * np.sin(OM * nh[:, None] * tw[None, :]) + + delta_Phi = np.sum(Cwaves + Swaves, axis=0) + + return delta_Phi + + def _d_deltaPhi_dt(self): + tw = self._tw() + waveamps = self._ORBWAVEs() + OM = self.ORBWAVE_OM.to("radian/second") + + nh = np.arange(self.nwaves) + 1 + Cwaves = ( + -OM + * nh[:, None] + * waveamps[:, 0, None] + * np.sin(OM * nh[:, None] * tw[None, :]) + ) + Swaves = ( + OM + * nh[:, None] + * waveamps[:, 1, None] + * np.cos(OM * nh[:, None] * tw[None, :]) + ) + + d_deltaPhi_dt = np.sum(Cwaves + Swaves, axis=0) + + return d_deltaPhi_dt.to(u.Unit("1/s"), equivalencies=u.dimensionless_angles()) + + def _d2_deltaPhi_dt2(self): + tw = self._tw() + waveamps = self._ORBWAVEs() + OM = self.ORBWAVE_OM.to("radian/second") + + nh = np.arange(self.nwaves) + 1 + Cwaves = ( + -((OM * nh[:, None]) ** 2) + * waveamps[:, 0, None] + * np.cos(OM * nh[:, None] * tw[None, :]) + ) + Swaves = ( + -((OM * nh[:, None]) ** 2) + * waveamps[:, 1, None] + * np.sin(OM * nh[:, None] * tw[None, :]) + ) + + d2_deltaPhi_dt2 = np.sum(Cwaves + Swaves, axis=0) + + return d2_deltaPhi_dt2.to( + u.Unit("1/s^2"), equivalencies=u.dimensionless_angles() + ) + + def orbits(self): + """Orbital phase (number of orbits since TASC).""" + PB = self.PB.to("second") + orbits = self.tt0 / PB + + dphi = self._deltaPhi() + + orbits += dphi + return orbits.decompose() + + def pbprime(self): + """Instantaneous binary period as a function of time.""" + PB = self.PB.to("second") + FB0 = 1.0 / PB + + FB0_shift = self._d_deltaPhi_dt() + + return 1.0 / (FB0 + FB0_shift).decompose() + + def pbdot_orbit(self): + """Reported value of PBDOT.""" + FB1 = self._d2_deltaPhi_dt2() + + return -(self.pbprime() ** 2) * FB1 + + def d_orbits_d_TASC(self): + PB = self.PB.to("second") + return -(1 / PB) * 2 * np.pi * u.rad + + def d_orbits_d_PB(self): + PB = self.PB.to("second") + return -(self.tt0 / PB**2).decompose() * 2 * np.pi * u.rad + + def d_orbits_d_orbwave(self, par): + tw = self._tw() + WOM = self.ORBWAVE_OM.to("radian/second") + nh = int(par[8:]) + 1 + + return ( + (np.cos(WOM * nh * tw) if par[7] == "C" else np.sin(WOM * nh * tw)) + * 2 + * np.pi + * u.rad + ) + + def d_pbprime_d_PB(self): + PB = self.PB.to("second") + FB0 = 1.0 / PB + + FB0_shift = self._d_deltaPhi_dt() + + FB0_prime = FB0 + FB0_shift + + return (1.0 / ((FB0_prime * PB) ** 2)).decompose() + + def d_pbprime_d_orbwave(self, par): + tw = self._tw() + WOM = self.ORBWAVE_OM.to("radian/second") + + nh = int(par[8:]) + 1 + if par[7] == "C": + d_deltaFB0_d_orbwave = -WOM * nh * np.sin(WOM * nh * tw) + else: + d_deltaFB0_d_orbwave = WOM * nh * np.cos(WOM * nh * tw) + + # use of pbprime should account for both Taylor and Fourier terms + return (-self.pbprime() ** 2 * d_deltaFB0_d_orbwave).to( + "d", equivalencies=u.dimensionless_angles() + ) + + def d_orbits_d_par(self, par): + return ( + self.d_orbits_d_orbwave(par) + if re.match(r"ORBWAVE[CS]\d+", par) is not None + else super().d_orbits_d_par(par) + ) + + def d_pbprime_d_par(self, par): + par_obj = getattr(self, par) + + if re.match(r"ORBWAVE[CS]\d+", par) is not None: + return self.d_pbprime_d_orbwave(par) + + return super().d_pbprime_d_par(par) + + +class OrbitWavesFBX(OrbitWaves): + """Orbit with orbital phase variations described both by FBX terms and + by a Fourier series""" + + def __init__( + self, + parent, + orbit_params=[ + "FB0", + "TASC", + "ORBWAVE_OM", + "ORBWAVE_EPOCH", + "ORBWAVEC0", + "ORBWAVES0", + ], + ): + if not hasattr(parent, "FB1"): + parent.add_param( + floatParameter(name="FB1", units=u.Hz**2, long_double=True, value=0) + ) + parent.binary_instance.add_binary_params("FB1", parent.FB1) + orbit_params += ["FB1"] + super().__init__(parent, orbit_params) + + def orbits(self): + """Orbital phase (number of orbits since TASC).""" + tt = self.tt0 + orbits = self.FB0 * tt * (1 + (0.5 * self.FB1 / self.FB0) * tt) + orbits += self._deltaPhi() + return orbits.decompose() + + def pbprime(self): + """Instantaneous binary period as a function of time.""" + + orbit_freq = self.FB0 + self.FB1 * self.tt0 + orbit_freq += self._d_deltaPhi_dt() + return 1.0 / orbit_freq.decompose() + + def pbdot_orbit(self): + """Reported value of PBDOT.""" + orbit_freq_dot = self._d2_deltaPhi_dt2() + self.FB1 + return -(self.pbprime() ** 2) * orbit_freq_dot + + def d_orbits_d_TASC(self): + return -self.FB0.to("Hz") * 2 * np.pi * u.rad + + def d_orbits_d_FB0(self): + return self.tt0.decompose() * (2 * np.pi * u.rad) + + def d_orbits_d_FB1(self): + return (0.5 * self.tt0**2).decompose() * (2 * np.pi * u.rad) + + def d_pbprime_d_FB0(self): + return -1 * self.pbprime() ** 2 + + def d_pbprime_d_FB1(self): + return -self.tt0 * self.pbprime() ** 2 + + def d_pbprime_d_TASC(self): + return self.FB1 * self.pbprime() ** 2 diff --git a/src/pint/models/timing_model.py b/src/pint/models/timing_model.py index 3c8fd217e..d6aeb0525 100644 --- a/src/pint/models/timing_model.py +++ b/src/pint/models/timing_model.py @@ -29,6 +29,7 @@ import abc import copy +import datetime import inspect import contextlib from collections import OrderedDict, defaultdict @@ -62,13 +63,23 @@ from pint.phase import Phase from pint.toa import TOAs from pint.utils import ( - PrefixError, split_prefixed_name, open_or_use, colorize, xxxselections, ) from pint.derived_quantities import dispersion_slope +from pint.exceptions import ( + PrefixError, + MissingTOAs, + PropertyAttributeError, + TimingModelError, + MissingParameter, + AliasConflict, + UnknownParameter, + UnknownBinaryModel, + MissingBinaryError, +) __all__ = [ @@ -76,11 +87,6 @@ "TimingModel", "Component", "AllComponents", - "TimingModelError", - "MissingParameter", - "MissingTOAs", - "MissingBinaryError", - "UnknownBinaryModel", ] # Parameters or lines in par files we don't understand but shouldn't # complain about. These are still passed to components so that they @@ -104,6 +110,9 @@ ignore_prefix = {"DMXF1_", "DMXF2_", "DMXEP_"} +# prefixes of parameters that may need to be checked for empty ranges +prefixes = ["DM", "SW", "CM"] + DEFAULT_ORDER = [ "astrometry", "jump_delay", @@ -123,26 +132,6 @@ ] -class MissingTOAs(ValueError): - """Some parameter does not describe any TOAs.""" - - def __init__(self, parameter_names): - if isinstance(parameter_names, str): - parameter_names = [parameter_names] - if len(parameter_names) == 1: - msg = f"Parameter {parameter_names[0]} does not correspond to any TOAs: you might need to run `model.find_empty_masks(toas, freeze=True)`" - elif len(parameter_names) > 1: - msg = f"Parameters {' '.join(parameter_names)} do not correspond to any TOAs: you might need to run `model.find_empty_masks(toas, freeze=True)`" - else: - raise ValueError("Incorrect attempt to construct MissingTOAs") - super().__init__(msg) - self.parameter_names = parameter_names - - -class PropertyAttributeError(ValueError): - pass - - def property_exists(f): """Mark a function as a property but handle AttributeErrors. @@ -240,6 +229,8 @@ class TimingModel: ---------- name : str The name of the timing model + meta : dict + A dictionary of metadata component_types : list A list of the distinct categories of component. For example, delay components will be register as 'DelayComponent'. @@ -254,6 +245,9 @@ def __init__(self, name="", components=[]): "First parameter should be the model name, was {!r}".format(name) ) self.name = name + self.meta = { + "read_time": f"{datetime.datetime.now().isoformat()}", + } self.component_types = [] self.top_level_params = [] self.add_param_from_top( @@ -2778,6 +2772,7 @@ def as_parfile( if include_info: info_string = pint.utils.info_string(prefix_string="# ", comment=comment) info_string += f"\n# Format: {format.lower()}" + info_string += "".join([f"\n# {x}: {self.meta[x]}" for x in self.meta]) result_begin = info_string + "\n" else: result_begin = "" @@ -2916,7 +2911,7 @@ def find_empty_masks(self, toas, freeze=False): if freeze: log.info(f"'{maskpar}' has no TOAs so freezing") getattr(self, maskpar).frozen = True - for prefix in ["DM", "SW"]: + for prefix in prefixes: mapping = pint.utils.xxxselections(self, toas, prefix=prefix) for k in mapping: if len(mapping[k]) == 0: @@ -3440,6 +3435,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. @@ -4023,68 +4019,3 @@ def param_to_unit(self, name): if getattr(self.components[component], firstname).unit_template is None: return self._param_unit_map[firstname] return u.Unit(getattr(self.components[component], firstname).unit_template(idx)) - - -class TimingModelError(ValueError): - """Generic base class for timing model errors.""" - - pass - - -class MissingParameter(TimingModelError): - """A required model parameter was not included. - - Parameters - ---------- - module - name of the model class that raised the error - param - name of the missing parameter - msg - additional message - - """ - - def __init__(self, module, param, msg=None): - super().__init__(msg) - self.module = module - self.param = param - self.msg = msg - - def __str__(self): - result = f"{self.module}.{self.param}" - if self.msg is not None: - result += "\n " + self.msg - return result - - -class AliasConflict(TimingModelError): - """If the same alias is used for different parameters.""" - - pass - - -class UnknownParameter(TimingModelError): - """Signal that a parameter name does not match any PINT parameters and their aliases.""" - - pass - - -class UnknownBinaryModel(TimingModelError): - """Signal that the par file requested a binary model not in PINT.""" - - def __init__(self, message, suggestion=None): - super().__init__(message) - self.suggestion = suggestion - - def __str__(self): - base_message = super().__str__() - if self.suggestion: - return f"{base_message} Perhaps use {self.suggestion}?" - return base_message - - -class MissingBinaryError(TimingModelError): - """Error for missing BINARY parameter.""" - - pass diff --git a/src/pint/models/wave.py b/src/pint/models/wave.py index 0a7cb1dfc..64e6da264 100644 --- a/src/pint/models/wave.py +++ b/src/pint/models/wave.py @@ -4,7 +4,8 @@ import numpy as np from pint.models.parameter import MJDParameter, floatParameter, prefixParameter -from pint.models.timing_model import PhaseComponent, MissingParameter +from pint.models.timing_model import PhaseComponent +from pint.exceptions import MissingParameter class Wave(PhaseComponent): @@ -14,9 +15,8 @@ class Wave(PhaseComponent): sine/cosine components. For consistency with the implementation in tempo2, this signal is treated - as a time series, but trivially converted into phase by multiplication by - F0, which could makes changes to PEPOCH fragile if there is strong spin - frequency evolution. + as a time series and trivially converted into phase by multiplication by + F0. Parameters supported: @@ -41,7 +41,7 @@ def __init__(self): floatParameter( name="WAVE_OM", description="Base frequency of wave solution", - units="1/d", + units="rad/d", convert_tcb2tdb=False, ) ) diff --git a/src/pint/models/wavex.py b/src/pint/models/wavex.py index 1714f3644..ecaf34f5a 100644 --- a/src/pint/models/wavex.py +++ b/src/pint/models/wavex.py @@ -6,14 +6,13 @@ from warnings import warn from pint.models.parameter import MJDParameter, prefixParameter -from pint.models.timing_model import DelayComponent, MissingParameter +from pint.models.timing_model import DelayComponent +from pint.exceptions import MissingParameter class WaveX(DelayComponent): """ - Implementation of the wave model as a delay correction - - Delays are expressed as a sum of sinusoids. + Implementation of the wave model as a delay correction, with delays are expressed as a sum of sinusoids. Used for decomposition of timing noise into a series of sine/cosine components with the amplitudes as fitted parameters. @@ -22,25 +21,34 @@ class WaveX(DelayComponent): .. paramtable:: :class: pint.models.wavex.WaveX - This is an extension of the L13 method described in Lentati et al., 2013 doi: 10.1103/PhysRevD.87.104021 - This model is similar to the TEMPO2 WAVE model parameters and users can convert a `TimingModel` with a Wave model - to a WaveX model and produce the same results. The main differences are that the WaveX frequencies are explicitly stated, + This is an extension of the method described in Lentati et al. (2013). + This model is similar to the TEMPO2 WAVE model parameters and users can convert a :class:`~pint.models/timing_model.TimingModel` + with a :class:`~pint.models.wave.Wave` model to a ``WaveX`` model and produce the same results. + The main differences are that the ``WaveX`` frequencies are explicitly stated, they do not necessarily need to be harmonics of some base frequency, the wave amplitudes are fittable parameters, and the - sine and cosine amplutides are reported as separate `prefixParameter`s rather than as a single `pairParameter`. + sine and cosine amplutides are reported as separate :class:`~pint.models.parameter.prefixParameter` rather than as a + single :class:`pint.models.parameter.pairParameter`. Analogous parameters in both models have the same units: - WAVEEPOCH is the same as WXEPOCH - WAVEOM and WXFREQ_000N have units of 1/d - WAVEN and WXSIN_000N/WXCOS_000N have units of seconds - The `pint.utils` functions `translate_wave_to_wavex()` and `translate_wavex_to_wave()` can be used to go back and forth between - two model. + - ``WAVEEPOCH`` is the same as ``WXEPOCH`` + - ``WAVEOM`` and ``WXFREQ_000N`` have units of 1/d + - ``WAVEN`` and ``WXSIN_000N/WXCOS_000N`` have units of seconds + + The :mod:`pint.utils` functions :func:`~pint.utils.translate_wave_to_wavex` and :func:`~pint.utils.translate_wavex_to_wave` + can be used to go back and forth between two model. + + WARNING: If the choice of ``WaveX`` frequencies in a :class:`~pint.models/timing_model.TimingModel` doesn't correspond to harmonics of some base + freqeuncy, it will not be possible to convert it to a :class:`~pint.models.wave.Wave` model. + + To set up a ``WaveX`` model, users can use the :mod:`pint.utils` function :func:`~pint.utils.wavex_setup` with either a list of frequencies or a choice + of harmonics of a base frequency determined by ``2 * pi /Timespan`` - WARNING: If the choice of WaveX frequencies in a `TimingModel` doesn't correspond to harmonics of some base - freqeuncy, it will not be possible to convert it to a Wave model. + References + ---------- + - Lentati et al. (2013), PRD, 87, 104021 [1]_ - To set up a WaveX model, users can use the `pint.utils` function `wavex_setup()` with either a list of frequencies or a choice - of harmonics of a base frequency determined by 2 * pi /Timespan + .. [1] https://ui.adsabs.harvard.edu/abs/2013PhRvD..87j4021L/abstract """ register = True diff --git a/src/pint/observatory/__init__.py b/src/pint/observatory/__init__.py index 8d98e28d7..6f9c3fbcb 100644 --- a/src/pint/observatory/__init__.py +++ b/src/pint/observatory/__init__.py @@ -37,6 +37,11 @@ from pint.config import runtimefile from pint.pulsar_mjd import Time from pint.utils import interesting_lines +from pint.exceptions import ( + ClockCorrectionError, + NoClockCorrections, + ClockCorrectionOutOfRange, +) # Include any files that define observatories here. This will start # with the standard distribution files, then will read any system- or @@ -93,24 +98,6 @@ def find_latest_bipm(): pint_clock_env_var = "PINT_CLOCK_OVERRIDE" -class ClockCorrectionError(RuntimeError): - """Unspecified error doing clock correction.""" - - pass - - -class NoClockCorrections(ClockCorrectionError): - """Clock corrections are expected but none are available.""" - - pass - - -class ClockCorrectionOutOfRange(ClockCorrectionError): - """Clock corrections are available but the requested time is not covered.""" - - pass - - # Global clock files shared by all observatories _gps_clock = None _bipm_clock_versions = {} diff --git a/src/pint/pintk/pulsar.py b/src/pint/pintk/pulsar.py index 7e8f3a80e..ce6aa8560 100644 --- a/src/pint/pintk/pulsar.py +++ b/src/pint/pintk/pulsar.py @@ -151,7 +151,9 @@ def __contains__(self, key): return key in self.prefit_model.params def reset_model(self): - self.prefit_model = pint.models.get_model(self.parfile) + self.prefit_model = pint.models.get_model( + self.parfile, allow_T2=True, allow_tcb=True + ) self.add_model_params() self.postfit_model = None self.postfit_resids = None @@ -172,7 +174,9 @@ def reset_TOAs(self): self.update_resids() def resetAll(self): - self.prefit_model = pint.models.get_model(self.parfile) + self.prefit_model = pint.models.get_model( + self.parfile, allow_T2=True, allow_tcb=True + ) self.postfit_model = None self.postfit_resids = None self.fitted = False @@ -603,9 +607,10 @@ def fit(self, selected, iters=4, compute_random=False): # adds extra prefix params for fitting self.add_model_params() - print( - f"Akaike information criterion = {akaike_information_criterion(self.fitter.model, self.fitter.toas)}" - ) + if not self.all_toas.is_wideband(): + print( + f"Akaike information criterion = {akaike_information_criterion(self.fitter.model, self.fitter.toas)}" + ) def random_models(self, selected): """Compute and plot random models""" diff --git a/src/pint/plot_utils.py b/src/pint/plot_utils.py index 45eba3977..289655f6d 100644 --- a/src/pint/plot_utils.py +++ b/src/pint/plot_utils.py @@ -283,7 +283,9 @@ def plot_priors( a, x = np.histogram(values[i], bins=bins, density=True) counts.append(a) - fig, axs = plt.subplots(len(keys), figsize=(8, 11), constrained_layout=True) + fig, axs = plt.subplots( + len(keys), figsize=(8, len(keys) * 1.5), constrained_layout=True + ) for i, p in enumerate(keys): if i != len(keys[:-1]): diff --git a/src/pint/scripts/compare_parfiles.py b/src/pint/scripts/compare_parfiles.py index 0c4ca5da2..d65785041 100644 --- a/src/pint/scripts/compare_parfiles.py +++ b/src/pint/scripts/compare_parfiles.py @@ -83,14 +83,25 @@ def main(argv=None): parser.add_argument( "-q", "--quiet", default=0, action="count", help="Decrease output verbosity" ) + parser.add_argument( + "--allow_tcb", + action="store_true", + help="Convert TCB par files to TDB automatically", + ) + parser.add_argument( + "--allow_T2", + action="store_true", + help="Guess the underlying binary model when T2 is given", + ) args = parser.parse_args(argv) pint.logging.setup( level=pint.logging.get_level(args.loglevel, args.verbosity, args.quiet) ) - m1 = get_model(args.input1) - m2 = get_model(args.input2) + m1 = get_model(args.input1, allow_T2=args.allow_T2, allow_tcb=args.allow_tcb) + m2 = get_model(args.input2, allow_T2=args.allow_T2, allow_tcb=args.allow_tcb) + print( m1.compare( m2, diff --git a/src/pint/scripts/convert_parfile.py b/src/pint/scripts/convert_parfile.py index 55b45415b..2975d5acc 100644 --- a/src/pint/scripts/convert_parfile.py +++ b/src/pint/scripts/convert_parfile.py @@ -73,6 +73,16 @@ def main(argv=None): parser.add_argument( "-q", "--quiet", default=0, action="count", help="Decrease output verbosity" ) + parser.add_argument( + "--allow_tcb", + action="store_true", + help="Convert TCB par files to TDB automatically", + ) + parser.add_argument( + "--allow_T2", + action="store_true", + help="Guess the underlying binary model when T2 is given", + ) args = parser.parse_args(argv) pint.logging.setup( @@ -83,7 +93,9 @@ def main(argv=None): return log.info(f"Reading '{args.input}'") - model = get_model(args.input) + + model = get_model(args.input, allow_T2=args.allow_T2, allow_tcb=args.allow_tcb) + if hasattr(model, "BINARY") and args.binary is not None: log.info(f"Converting from {model.BINARY.value} to {args.binary}") if args.binary == "ELL1H": diff --git a/src/pint/scripts/event_optimize.py b/src/pint/scripts/event_optimize.py index e4c9dd4ff..6e26ba784 100755 --- a/src/pint/scripts/event_optimize.py +++ b/src/pint/scripts/event_optimize.py @@ -414,6 +414,16 @@ def __init__( self.model, phs, phserr ) self.n_fit_params = len(self.fitvals) + self.M, _, _ = self.model.designmatrix(self.toas) + self.M = self.M.transpose() * -self.model.F0.value + self.phases = self.get_event_phases() + self.linearize_model = False + + def calc_phase_matrix(self, theta): + d_phs = np.zeros(len(self.toas)) + for i in range(len(theta) - 1): + d_phs += self.M[i + 1] * (self.fitvals[i] - theta[i]) + return (self.phases - d_phs) % 1 def get_event_phases(self): """ @@ -446,7 +456,10 @@ def lnposterior(self, theta): return -np.inf, -np.inf, -np.inf # Call PINT to compute the phases - phases = self.get_event_phases() + if self.linearize_model: + phases = self.calc_phase_matrix(theta) + else: + phases = self.get_event_phases() lnlikelihood = profile_likelihood( theta[-1], self.xtemp, phases, self.template, self.weights ) @@ -686,6 +699,23 @@ def main(argv=None): action="store_true", dest="noautocorr", ) + parser.add_argument( + "--linearize_model", + help="Calculates the phase at each MCMC step using the designmatrix", + default=False, + action="store_true", + dest="linearize_model", + ) + parser.add_argument( + "--allow_tcb", + action="store_true", + help="Convert TCB par files to TDB automatically", + ) + parser.add_argument( + "--allow_T2", + action="store_true", + help="Guess the underlying binary model when T2 is given", + ) args = parser.parse_args(argv) pint.logging.setup( @@ -719,7 +749,9 @@ def main(argv=None): ncores = args.ncores # Read in initial model - modelin = pint.models.get_model(parfile) + modelin = pint.models.get_model( + parfile, allow_T2=args.allow_T2, allow_tcb=args.allow_tcb + ) # File name setup and clobber file check filepath = args.filepath or os.getcwd() @@ -862,6 +894,9 @@ def main(argv=None): # This way, one walker should always be in a good position pos[0] = ftr.fitvals + # How phase will be calculated at each step (either with the designmatrix orthe exact phase calculation) + ftr.linearize_model = args.linearize_model + import emcee # Setting up a backend to save the chains into an h5 file @@ -925,7 +960,7 @@ def chains_to_dict(names, sampler): def plot_chains(chain_dict, file=False): npts = len(chain_dict) - fig, axes = plt.subplots(npts, 1, sharex=True, figsize=(8, 9)) + fig, axes = plt.subplots(npts, 1, sharex=True, figsize=(8, npts * 1.5)) for ii, name in enumerate(chain_dict.keys()): axes[ii].plot(chain_dict[name], color="k", alpha=0.3) axes[ii].set_ylabel(name) @@ -950,6 +985,7 @@ def plot_chains(chain_dict, file=False): lnprior_samps = blobs["lnprior"] lnlikelihood_samps = blobs["lnlikelihood"] lnpost_samps = lnprior_samps + lnlikelihood_samps + maxpost = lnpost_samps[:][burnin:].max() ind = np.unravel_index( np.argmax(lnpost_samps[:][burnin:]), lnpost_samps[:][burnin:].shape ) @@ -1000,8 +1036,15 @@ def plot_chains(chain_dict, file=False): ] ftr.set_param_uncertainties(dict(zip(ftr.fitkeys[:-1], errors[:-1]))) + # Calculating the AIC and BIC + n_params = len(ftr.model.free_params) + AIC = 2 * (n_params - maxpost) + BIC = n_params * np.log(len(ts)) - 2 * maxpost + ftr.model.NTOA.value = ts.ntoas f = open(filename + "_post.par", "w") f.write(ftr.model.as_parfile()) + f.write(f"\n#The AIC is {AIC}") + f.write(f"\n#The BIC is {BIC}") f.close() # Print the best MCMC values and ranges diff --git a/src/pint/scripts/event_optimize_MCMCFitter.py b/src/pint/scripts/event_optimize_MCMCFitter.py index bbab8ccee..7357b1028 100755 --- a/src/pint/scripts/event_optimize_MCMCFitter.py +++ b/src/pint/scripts/event_optimize_MCMCFitter.py @@ -128,6 +128,17 @@ def main(argv=None): help="Logging level", dest="loglevel", ) + parser.add_argument( + "--allow_tcb", + action="store_true", + help="Convert TCB par files to TDB automatically", + ) + parser.add_argument( + "--allow_T2", + action="store_true", + help="Guess the underlying binary model when T2 is given", + ) + global nwalkers, nsteps, ftr args = parser.parse_args(argv) @@ -164,7 +175,9 @@ def main(argv=None): wgtexp = args.wgtexp # Read in initial model - modelin = pint.models.get_model(parfile) + modelin = pint.models.get_model( + parfile, allow_T2=args.allow_T2, allow_tcb=args.allow_tcb + ) # The custom_timing version below is to manually construct the TimingModel # class, which allows it to be pickled. This is needed for parallelizing diff --git a/src/pint/scripts/event_optimize_multiple.py b/src/pint/scripts/event_optimize_multiple.py index 41316a000..0a3dc7336 100755 --- a/src/pint/scripts/event_optimize_multiple.py +++ b/src/pint/scripts/event_optimize_multiple.py @@ -228,6 +228,16 @@ def main(argv=None): help="Logging level", dest="loglevel", ) + parser.add_argument( + "--allow_tcb", + action="store_true", + help="Convert TCB par files to TDB automatically", + ) + parser.add_argument( + "--allow_T2", + action="store_true", + help="Guess the underlying binary model when T2 is given", + ) global nwalkers, nsteps, ftr @@ -261,7 +271,9 @@ def main(argv=None): wgtexp = args.wgtexp # Read in initial model - modelin = pint.models.get_model(parfile) + modelin = pint.models.get_model( + parfile, allow_T2=args.allow_T2, allow_tcb=args.allow_tcb + ) # Set the target coords for automatic weighting if necessary if "ELONG" in modelin.params: diff --git a/src/pint/scripts/fermiphase.py b/src/pint/scripts/fermiphase.py index b1427e61f..c5e43bb93 100755 --- a/src/pint/scripts/fermiphase.py +++ b/src/pint/scripts/fermiphase.py @@ -77,6 +77,16 @@ def main(argv=None): parser.add_argument( "-q", "--quiet", default=0, action="count", help="Decrease output verbosity" ) + parser.add_argument( + "--allow_tcb", + action="store_true", + help="Convert TCB par files to TDB automatically", + ) + parser.add_argument( + "--allow_T2", + action="store_true", + help="Guess the underlying binary model when T2 is given", + ) args = parser.parse_args(argv) pint.logging.setup( @@ -88,7 +98,10 @@ def main(argv=None): args.addphase = True # Read in model - modelin = pint.models.get_model(args.parfile) + modelin = pint.models.get_model( + args.parfile, allow_T2=args.allow_T2, allow_tcb=args.allow_tcb + ) + if "ELONG" in modelin.params: tc = SkyCoord( modelin.ELONG.quantity, diff --git a/src/pint/scripts/photonphase.py b/src/pint/scripts/photonphase.py index 7d72eec18..2798e0a7f 100755 --- a/src/pint/scripts/photonphase.py +++ b/src/pint/scripts/photonphase.py @@ -106,12 +106,25 @@ def main(argv=None): help="Logging level", dest="loglevel", ) + parser.add_argument( + "--nbin", help="Number of phase bins in the phaseogram", default=100, type=int + ) parser.add_argument( "-v", "--verbosity", default=0, action="count", help="Increase output verbosity" ) parser.add_argument( "-q", "--quiet", default=0, action="count", help="Decrease output verbosity" ) + parser.add_argument( + "--allow_tcb", + action="store_true", + help="Convert TCB par files to TDB automatically", + ) + parser.add_argument( + "--allow_T2", + action="store_true", + help="Guess the underlying binary model when T2 is given", + ) args = parser.parse_args(argv) pint.logging.setup( @@ -153,7 +166,10 @@ def main(argv=None): "Please barycenter the event file using the official mission tools before processing with PINT" ) # Read in model - modelin = pint.models.get_model(args.parfile) + modelin = pint.models.get_model( + args.parfile, allow_T2=args.allow_T2, allow_tcb=args.allow_tcb + ) + use_planets = False if "PLANET_SHAPIRO" in modelin.params: if modelin.PLANET_SHAPIRO.value: @@ -254,7 +270,7 @@ def main(argv=None): print("Htest : {0:.2f} ({1:.2f} sigma)".format(h, h2sig(h))) if args.plot: - phaseogram_binned(mjds, phases, bins=100, plotfile=args.plotfile) + phaseogram_binned(mjds, phases, bins=args.nbin, plotfile=args.plotfile) # Compute orbital phases for each photon TOA if args.addorbphase: diff --git a/src/pint/scripts/pintbary.py b/src/pint/scripts/pintbary.py index 4876474d8..61cdcd0a0 100755 --- a/src/pint/scripts/pintbary.py +++ b/src/pint/scripts/pintbary.py @@ -74,6 +74,16 @@ def main(argv=None): parser.add_argument( "-q", "--quiet", default=0, action="count", help="Decrease output verbosity" ) + parser.add_argument( + "--allow_tcb", + action="store_true", + help="Convert TCB par files to TDB automatically", + ) + parser.add_argument( + "--allow_T2", + action="store_true", + help="Guess the underlying binary model when T2 is given", + ) args = parser.parse_args(argv) pint.logging.setup( @@ -105,7 +115,9 @@ def main(argv=None): ) if args.parfile is not None: - m = pint.models.get_model(args.parfile) + m = pint.models.get_model( + args.parfile, allow_T2=args.allow_T2, allow_tcb=args.allow_tcb + ) else: # Construct model by hand m = pint.models.StandardTimingModel diff --git a/src/pint/scripts/pintempo.py b/src/pint/scripts/pintempo.py index d99c772ba..ac9099008 100755 --- a/src/pint/scripts/pintempo.py +++ b/src/pint/scripts/pintempo.py @@ -62,6 +62,16 @@ def main(argv=None): parser.add_argument( "-q", "--quiet", default=0, action="count", help="Decrease output verbosity" ) + parser.add_argument( + "--allow_tcb", + action="store_true", + help="Convert TCB par files to TDB automatically", + ) + parser.add_argument( + "--allow_T2", + action="store_true", + help="Guess the underlying binary model when T2 is given", + ) args = parser.parse_args(argv) pint.logging.setup( @@ -69,7 +79,9 @@ def main(argv=None): ) log.info("Reading model from {0}".format(args.parfile)) - m = pint.models.get_model(args.parfile) + m = pint.models.get_model( + args.parfile, allow_T2=args.allow_T2, allow_tcb=args.allow_tcb + ) log.warning(m.params) diff --git a/src/pint/scripts/pintpublish.py b/src/pint/scripts/pintpublish.py index 6d5b0142f..1828a05c6 100644 --- a/src/pint/scripts/pintpublish.py +++ b/src/pint/scripts/pintpublish.py @@ -59,10 +59,22 @@ def main(argv=None): action="store_true", default=False, ) + parser.add_argument( + "--allow_tcb", + action="store_true", + help="Convert TCB par files to TDB automatically", + ) + parser.add_argument( + "--allow_T2", + action="store_true", + help="Guess the underlying binary model when T2 is given", + ) args = parser.parse_args(argv) - model, toas = get_model_and_toas(args.parfile, args.timfile) + model, toas = get_model_and_toas( + args.parfile, args.timfile, allow_T2=args.allow_T2, allow_tcb=args.allow_tcb + ) output = publish( model, diff --git a/src/pint/scripts/t2binary2pint.py b/src/pint/scripts/t2binary2pint.py index 42d70f89a..bab43eef4 100644 --- a/src/pint/scripts/t2binary2pint.py +++ b/src/pint/scripts/t2binary2pint.py @@ -45,12 +45,17 @@ def main(argv=None): default=True, help="Whether to drop SINI if the model is DDK (True)", ) + parser.add_argument( + "--allow_tcb", + action="store_true", + help="Convert TCB par files to TDB automatically", + ) args = parser.parse_args(argv) mb = ModelBuilder() - model = mb(args.input_par, allow_T2=True, allow_tcb=True) + model = mb(args.input_par, allow_T2=True, allow_tcb=args.allow_tcb) model.write_parfile(args.output_par) print(f"Output written to {args.output_par}") diff --git a/src/pint/scripts/tcb2tdb.py b/src/pint/scripts/tcb2tdb.py index 4c427d21a..02a485d3c 100644 --- a/src/pint/scripts/tcb2tdb.py +++ b/src/pint/scripts/tcb2tdb.py @@ -30,11 +30,16 @@ def main(argv=None): ) parser.add_argument("input_par", help="Input par file name (TCB)") parser.add_argument("output_par", help="Output par file name (TDB)") + parser.add_argument( + "--allow_T2", + action="store_true", + help="Guess the underlying binary model when T2 is given", + ) args = parser.parse_args(argv) mb = ModelBuilder() - model = mb(args.input_par, allow_tcb=True) + model = mb(args.input_par, allow_tcb=True, allow_T2=args.allow_T2) model.write_parfile(args.output_par) log.info(f"Output written to {args.output_par}.") diff --git a/src/pint/scripts/zima.py b/src/pint/scripts/zima.py index 3305129f7..29b3e354c 100755 --- a/src/pint/scripts/zima.py +++ b/src/pint/scripts/zima.py @@ -115,6 +115,16 @@ def main(argv=None): parser.add_argument( "-q", "--quiet", default=0, action="count", help="Decrease output verbosity" ) + parser.add_argument( + "--allow_tcb", + action="store_true", + help="Convert TCB par files to TDB automatically", + ) + parser.add_argument( + "--allow_T2", + action="store_true", + help="Guess the underlying binary model when T2 is given", + ) args = parser.parse_args(argv) pint.logging.setup( @@ -122,7 +132,9 @@ def main(argv=None): ) log.info("Reading model from {0}".format(args.parfile)) - m = pint.models.get_model(args.parfile) + m = pint.models.get_model( + args.parfile, allow_T2=args.allow_T2, allow_tcb=args.allow_tcb + ) out_format = args.format error = args.error * u.microsecond diff --git a/src/pint/utils.py b/src/pint/utils.py index 00e9c4540..2b78bce6d 100644 --- a/src/pint/utils.py +++ b/src/pint/utils.py @@ -75,9 +75,9 @@ import pint.pulsar_ecliptic from pint.toa_select import TOASelect from pint.types import file_like, quantity_like +from pint.exceptions import PINTPrecisionError, PrefixError __all__ = [ - "PINTPrecisionError", "check_longdouble_precision", "require_longdouble_precision", "PosVel", @@ -154,10 +154,6 @@ # Actual exported tools -class PINTPrecisionError(RuntimeError): - pass - - # A warning is emitted in pint.pulsar_mjd if sufficient precision is not available @@ -361,14 +357,10 @@ def has_astropy_unit(x: Any) -> bool: re.compile(r"^([a-zA-Z]*\d+[a-zA-Z]+)(\d+)$"), # For the prefix like T2EFAC2 re.compile(r"^([a-zA-Z]+)0*(\d+)$"), # For the prefix like F12 re.compile(r"^([a-zA-Z0-9]+_)(\d+)$"), # For the prefix like DMXR1_3 - # re.compile(r'([a-zA-Z]\d[a-zA-Z]+)(\d+)'), # for prefixes like PLANET_SHAPIRO2? + re.compile(r"([a-zA-Z]+_[a-zA-Z]+)(\d+)$"), # for prefixes like NE_SW2? ] -class PrefixError(ValueError): - pass - - def split_prefixed_name(name: str) -> Tuple[str, str, int]: """Split a prefixed name. @@ -1753,11 +1745,8 @@ def _translate_wave_freqs(om: Union[float, u.Quantity], k: int) -> u.Quantity: astropy.units.Quantity WXFREQ_ quantity in units 1/d that can be used in WaveX model """ - if isinstance(om, u.quantity.Quantity): - om.to(u.d**-1) - else: - om *= u.d**-1 - return (om * (k + 1)) / (2.0 * np.pi) + om <<= u.rad / u.d + return (om * (k + 1)) / (2.0 * np.pi * u.rad) def _translate_wavex_freqs(wxfreq: Union[float, u.Quantity], k: int) -> u.Quantity: @@ -1777,13 +1766,12 @@ def _translate_wavex_freqs(wxfreq: Union[float, u.Quantity], k: int) -> u.Quanti astropy.units.Quantity WAVEOM quantity in units 1/d that can be used in Wave model """ - if isinstance(wxfreq, u.quantity.Quantity): - wxfreq.to(u.d**-1) - else: - wxfreq *= u.d**-1 + wxfreq <<= u.d**-1 if len(wxfreq) == 1: - return (2.0 * np.pi * wxfreq) / (k + 1.0) - wave_om = [((2.0 * np.pi * wxfreq[i]) / (k[i] + 1.0)) for i in range(len(wxfreq))] + return (2.0 * np.pi * u.rad * wxfreq) / (k + 1.0) + wave_om = [ + ((2.0 * np.pi * u.rad * wxfreq[i]) / (k[i] + 1.0)) for i in range(len(wxfreq)) + ] return ( sum(wave_om) / len(wave_om) if np.allclose(wave_om, wave_om[0], atol=1e-3) diff --git a/tests/datafile/J1048+2339_3PC_fake.tim b/tests/datafile/J1048+2339_3PC_fake.tim new file mode 100644 index 000000000..33d3f61b0 --- /dev/null +++ b/tests/datafile/J1048+2339_3PC_fake.tim @@ -0,0 +1,107 @@ +FORMAT 1 +C Created: 2024-01-17T17:35:49.743600 +C PINT_version: 0.9.8+44.g1d25d7ec.dirty +C User: cclark +C Host: HNB775 +C OS: macOS-11.7.10-x86_64-i386-64bit +C Python: 3.11.0 | packaged by conda-forge | (main, Jan 15 2023, 05:44:48) [Clang 14.0.6 ] +fake 1400.000000 56499.6542924265793634 20.000 arecibo -pn -14556992621.0 +fake 1400.000000 56515.4714585239121180 20.000 arecibo -pn -14264065295.0 +fake 1400.000000 56530.1658492425925694 20.000 arecibo -pn -13991925153.0 +fake 1400.000000 56552.8508442386096066 20.000 arecibo -pn -13571785893.0 +fake 1400.000000 56573.1875466751208218 20.000 arecibo -pn -13195126377.0 +fake 1400.000000 56566.0911291819698843 20.000 arecibo -pn -13326561640.0 +fake 1400.000000 56570.6109755583843172 20.000 arecibo -pn -13242848100.0 +fake 1400.000000 56593.0993371540435069 20.000 arecibo -pn -12826327673.0 +fake 1400.000000 56597.6169268134667246 20.000 arecibo -pn -12742653507.0 +fake 1400.000000 56584.9166539062826851 20.000 arecibo -pn -12977884987.0 +fake 1400.000000 56628.4981058130401389 20.000 arecibo -pn -12170671808.0 +fake 1400.000000 56629.6491682273526967 20.000 arecibo -pn -12149352121.0 +fake 1400.000000 56631.3948759159077430 20.000 arecibo -pn -12117018055.0 +fake 1400.000000 56647.7043342506005207 20.000 arecibo -pn -11814934720.0 +fake 1400.000000 56664.0323435136898380 20.000 arecibo -pn -11512511204.0 +fake 1400.000000 56662.6202116121144444 20.000 arecibo -pn -11538666096.0 +fake 1400.000000 56672.5319303522428934 20.000 arecibo -pn -11355086157.0 +fake 1400.000000 56684.7216506798818518 20.000 arecibo -pn -11129317426.0 +fake 1400.000000 56707.5113283280297106 20.000 arecibo -pn -10707235585.0 +fake 1400.000000 56718.1584645803764815 20.000 arecibo -pn -10510047952.0 +fake 1400.000000 56724.6996789532842129 20.000 arecibo -pn -10388904838.0 +fake 1400.000000 56736.8561543310538079 20.000 arecibo -pn -10163771090.0 +fake 1400.000000 56742.9772745043780441 20.000 arecibo -pn -10050411113.0 +fake 1400.000000 56744.8663333106221991 20.000 arecibo -pn -10015427356.0 +fake 1400.000000 56768.2681409742178935 20.000 arecibo -pn -9582048339.0 +fake 1400.000000 56772.0303959887401968 20.000 arecibo -pn -9512376009.0 +fake 1400.000000 56791.3104227062363658 20.000 arecibo -pn -9155336606.0 +fake 1400.000000 56797.7828604962445255 20.000 arecibo -pn -9035476833.0 +fake 1400.000000 56828.4317284397404398 20.000 arecibo -pn -8467904208.0 +fake 1400.000000 56812.0591817380739699 20.000 arecibo -pn -8771100758.0 +fake 1400.000000 56815.4905032654365741 20.000 arecibo -pn -8707557734.0 +fake 1400.000000 56841.7902342958811921 20.000 arecibo -pn -8220521613.0 +fake 1400.000000 56837.2661549845360417 20.000 arecibo -pn -8304302163.0 +fake 1400.000000 56848.5214910022805440 20.000 arecibo -pn -8095865857.0 +fake 1400.000000 56883.4384623373841320 20.000 arecibo -pn -7449224561.0 +fake 1400.000000 56880.8078014433801619 20.000 arecibo -pn -7497943625.0 +fake 1400.000000 56863.8792153396387384 20.000 arecibo -pn -7811454003.0 +fake 1400.000000 56888.1516839050684954 20.000 arecibo -pn -7361935924.0 +fake 1400.000000 56900.3285856223282407 20.000 arecibo -pn -7136417053.0 +fake 1400.000000 56926.5954650101420255 20.000 arecibo -pn -6649934542.0 +fake 1400.000000 56940.6035905419990163 20.000 arecibo -pn -6390485783.0 +fake 1400.000000 56936.8112658763676505 20.000 arecibo -pn -6460725217.0 +fake 1400.000000 56952.1893387214877662 20.000 arecibo -pn -6175898993.0 +fake 1400.000000 56955.1061180325763079 20.000 arecibo -pn -6121875452.0 +fake 1400.000000 56955.3718469374112500 20.000 arecibo -pn -6116953617.0 +fake 1400.000000 56980.1841145056263194 20.000 arecibo -pn -5657382463.0 +fake 1400.000000 56999.6866740820557408 20.000 arecibo -pn -5296155314.0 +fake 1400.000000 57003.1191636566178009 20.000 arecibo -pn -5232578657.0 +fake 1400.000000 57008.9625675219576159 20.000 arecibo -pn -5124347386.0 +fake 1400.000000 57039.6211866433689237 20.000 arecibo -pn -4556497914.0 +fake 1400.000000 57035.2024855226981944 20.000 arecibo -pn -4638338817.0 +fake 1400.000000 57030.6871973368739236 20.000 arecibo -pn -4721968796.0 +fake 1400.000000 57061.5047784739957523 20.000 arecibo -pn -4151190878.0 +fake 1400.000000 57057.6151107016963889 20.000 arecibo -pn -4223230403.0 +fake 1400.000000 57077.0368741394684028 20.000 arecibo -pn -3863528613.0 +fake 1400.000000 57070.3456692022659491 20.000 arecibo -pn -3987452131.0 +fake 1400.000000 57110.1177331390579977 20.000 arecibo -pn -3250878928.0 +fake 1400.000000 57112.4934863738438657 20.000 arecibo -pn -3206881410.0 +fake 1400.000000 57113.2766723567323727 20.000 arecibo -pn -3192377388.0 +fake 1400.000000 57134.0700569265716551 20.000 arecibo -pn -2807304956.0 +fake 1400.000000 57146.6515469834735648 20.000 arecibo -pn -2574312382.0 +fake 1400.000000 57149.7566600342336226 20.000 arecibo -pn -2516810137.0 +fake 1400.000000 57165.2526190755023611 20.000 arecibo -pn -2229847804.0 +fake 1400.000000 57187.2688830774209259 20.000 arecibo -pn -1822139031.0 +fake 1400.000000 57190.6849732852207292 20.000 arecibo -pn -1758877503.0 +fake 1400.000000 57187.9589681610874190 20.000 arecibo -pn -1809359372.0 +fake 1400.000000 57200.3717603740499769 20.000 arecibo -pn -1579491245.0 +fake 1400.000000 57210.6059847836613542 20.000 arecibo -pn -1389965223.0 +fake 1400.000000 57227.0308377842580208 20.000 arecibo -pn -1085791426.0 +fake 1400.000000 57236.8596851869279052 20.000 arecibo -pn -903767168.0 +fake 1400.000000 57261.2487812880531366 20.000 arecibo -pn -452083881.0 +fake 1400.000000 57260.4442328007430556 20.000 arecibo -pn -466984383.0 +fake 1400.000000 57268.5769349500900232 20.000 arecibo -pn -316363616.0 +fake 1400.000000 57279.9979853937153241 20.000 arecibo -pn -104838855.0 +fake 1400.000000 57295.6360907784355671 20.000 arecibo -pn 184795178.0 +fake 1400.000000 57301.3161497624851158 20.000 arecibo -pn 289996962.0 +fake 1400.000000 57317.2132137646926157 20.000 arecibo -pn 584435547.0 +fake 1400.000000 57305.1445115133814120 20.000 arecibo -pn 360903928.0 +fake 1400.000000 57331.7799537955628588 20.000 arecibo -pn 854238088.0 +fake 1400.000000 57329.0529708326988542 20.000 arecibo -pn 803729018.0 +fake 1400.000000 57337.5205917926200810 20.000 arecibo -pn 960566046.0 +fake 1400.000000 57365.2644749378174190 20.000 arecibo -pn 1474439290.0 +fake 1400.000000 57357.4298283792541783 20.000 arecibo -pn 1329325411.0 +fake 1400.000000 57373.5424751155939351 20.000 arecibo -pn 1627764289.0 +fake 1400.000000 57388.6568674133486228 20.000 arecibo -pn 1907710417.0 +fake 1400.000000 57393.9312912733891782 20.000 arecibo -pn 2005401471.0 +fake 1400.000000 57410.9067652986654399 20.000 arecibo -pn 2319811514.0 +fake 1400.000000 57421.2851809244780440 20.000 arecibo -pn 2512030434.0 +fake 1400.000000 57432.9383376233199190 20.000 arecibo -pn 2727855790.0 +fake 1400.000000 57439.2307955621883102 20.000 arecibo -pn 2844394914.0 +fake 1400.000000 57465.9557426709758101 20.000 arecibo -pn 3339339546.0 +fake 1400.000000 57473.2128096176505093 20.000 arecibo -pn 3473736409.0 +fake 1400.000000 57469.1446187069599884 20.000 arecibo -pn 3398395700.0 +fake 1400.000000 57473.7402821754510069 20.000 arecibo -pn 3483504941.0 +fake 1400.000000 57494.4992201782944907 20.000 arecibo -pn 3867941489.0 +fake 1400.000000 57516.7439985751001505 20.000 arecibo -pn 4279885696.0 +fake 1400.000000 57498.9767069652037038 20.000 arecibo -pn 3950859181.0 +fake 1400.000000 57540.7899028004051968 20.000 arecibo -pn 4725180159.0 +fake 1400.000000 57518.1608066433601852 20.000 arecibo -pn 4306123159.0 +fake 1400.000000 57539.6812402803568286 20.000 arecibo -pn 4704649696.0 diff --git a/tests/datafile/J1048+2339_orbwaves.par b/tests/datafile/J1048+2339_orbwaves.par new file mode 100644 index 000000000..a1f230cd8 --- /dev/null +++ b/tests/datafile/J1048+2339_orbwaves.par @@ -0,0 +1,48 @@ +PSR J1048+2339 +RAJ 10:48:43.41835421 1 0.00007631 +DECJ 23:39:53.4043452 1 0.0020013 +PMRA -18.7000 +PMDEC -9.4000 +POSEPOCH 56700.0000 +F0 214.3547853411294852 1 0.0000000000170721 +F1 -1.381753144359e-15 1 1.531201039458e-18 +PEPOCH 57285.000000 +START 56508.798 +FINISH 57548.936 +#DM 16.654332 +DMEPOCH 56700.0000 +SOLARN0 10.00 +EPHEM DE421 +CLK TT(BIPM2021) +UNITS TDB +TIMEEPH FB90 +T2CMETHOD IAU2000B +CORRECT_TROPOSPHERE N +PLANET_SHAPIRO N +DILATEFREQ N +NTOA 104 +TRES 8.84 +TZRMJD 57285.65854990785506 +TZRFRQ 326.988 +TZRSITE 3 +NITS 1 +BINARY ELL1 +A1 0.836120368 1 0.000003165 +EPS1 0 +EPS2 0 +# ATNF DM info (from drc+16): +DM 16.6544 0.0001 +ORBWAVE_OM 3.495788643739122e-08 +ORBWAVE_EPOCH 57028.867 +TASC 56637.61958641038 1 +PB 0.250506696660119 1 +ORBWAVEC0 0.01359122306797746 1 +ORBWAVES0 -0.1138721359160273 1 +ORBWAVEC1 -0.009093465738757574 1 +ORBWAVES1 0.03712548529066605 1 +ORBWAVEC2 0.00414066930830874 1 +ORBWAVES2 -0.011358879259751396 1 +ORBWAVEC3 -0.001400204000025873 1 +ORBWAVES3 0.002749601047126059 1 +ORBWAVEC4 0.00019385501703957642 1 +ORBWAVES4 -0.00028698933945954427 1 \ No newline at end of file diff --git a/tests/datafile/J1048+2339_orbwaves_DD.par b/tests/datafile/J1048+2339_orbwaves_DD.par new file mode 100644 index 000000000..acbc3dd11 --- /dev/null +++ b/tests/datafile/J1048+2339_orbwaves_DD.par @@ -0,0 +1,48 @@ +PSR J1048+2339 +RAJ 10:48:43.41835421 1 0.00007631 +DECJ 23:39:53.4043452 1 0.0020013 +PMRA -18.7000 +PMDEC -9.4000 +POSEPOCH 56700.0000 +F0 214.3547853411294852 1 0.0000000000170721 +F1 -1.381753144359e-15 1 1.531201039458e-18 +PEPOCH 57285.000000 +START 56508.798 +FINISH 57548.936 +#DM 16.654332 +DMEPOCH 56700.0000 +SOLARN0 10.00 +EPHEM DE421 +CLK TT(BIPM2021) +UNITS TDB +TIMEEPH FB90 +T2CMETHOD IAU2000B +CORRECT_TROPOSPHERE N +PLANET_SHAPIRO N +DILATEFREQ N +NTOA 104 +TRES 8.84 +TZRMJD 57285.65854990785506 +TZRFRQ 326.988 +TZRSITE 3 +NITS 1 +BINARY DD +A1 0.836120368 1 0.000003165 +E 0 +OM 0 +# ATNF DM info (from drc+16): +DM 16.6544 0.0001 +ORBWAVE_OM 3.495788643739122e-08 +ORBWAVE_EPOCH 57028.867 +T0 56637.61958641038 1 +PB 0.250506696660119 1 +ORBWAVEC0 0.01359122306797746 1 +ORBWAVES0 -0.1138721359160273 1 +ORBWAVEC1 -0.009093465738757574 1 +ORBWAVES1 0.03712548529066605 1 +ORBWAVEC2 0.00414066930830874 1 +ORBWAVES2 -0.011358879259751396 1 +ORBWAVEC3 -0.001400204000025873 1 +ORBWAVES3 0.002749601047126059 1 +ORBWAVEC4 0.00019385501703957642 1 +ORBWAVES4 -0.00028698933945954427 1 \ No newline at end of file diff --git a/tests/datafile/vela_wave.par b/tests/datafile/vela_wave.par index e8aba5a4c..889361b01 100644 --- a/tests/datafile/vela_wave.par +++ b/tests/datafile/vela_wave.par @@ -28,45 +28,18 @@ TZRMJD 55896.55312516020475 TZRFRQ 1370.2919919999999365 TZRSITE pks TRES 310.453 -EPHVER 2 CLK TT(TAI) MODE 1 UNITS TDB TIMEEPH FB90 DILATEFREQ N PLANET_SHAPIRO N -T2CMETHOD TEMPO NE_SW 9.961 CORRECT_TROPOSPHERE N EPHEM DE405 NITS 1 NTOA 339 CHI2R 328088.1294 298 -JUMP -B 20CM 0 0 -JUMP -B 40CM 0 0 -JUMP -B 50CM 0 0 -JUMP -dfb3_J0437_55319_56160 1 2.2e-07 0 -JUMP -dfb3_J0437_56160_60000 1 4.5e-07 0 -JUMP -pdfb1_128_ch 1 -3.5547e-06 0 -JUMP -pdfb1_2048_ch 1 -4.68829e-05 0 -JUMP -pdfb1_512_ch 1 -1.22813e-05 0 -JUMP -pdfb1_post_2006 1 -1.3e-07 0 -JUMP -pdfb1_pre_2006 1 -1.13e-06 0 -JUMP -pdfb2_1024_MHz 1 -5.435e-06 0 -JUMP -pdfb2_256MHz_1024_ch 1 -1.1395e-05 0 -JUMP -pdfb2_256MHz_512ch 1 -4.75e-06 0 -JUMP -pdfb3_1024_256_512 1 2.45e-06 0 -JUMP -pdfb3_1024_MHz 1 1.03e-06 0 -JUMP -pdfb3_256MHz_1024ch 1 4.295e-06 0 -JUMP -pdfb3_64MHz_1024ch 1 1.494e-05 0 -JUMP -pdfb3_64MHz_512ch 1 8.9e-06 0 -JUMP -pdfb4.*_1024_[1,2]... 1 2.23e-06 0 -JUMP -pdfb4_256MHz_1024ch 1 5.05e-06 0 -JUMP -pdfb4_55319_56055_cals 1 9.27e-07 0 -JUMP -pdfb4_56110_56160_cals 1 5.41e-07 0 -JUMP -pdfb4_56160_60000_cals 1 4.25e-07 0 -JUMP -wbb256_512_128_3p_b 1 -6.2e-07 0 -JUMP -wbb_c_config 1 3.8e-07 0 WAVEEPOCH 55305 WAVE_OM 0.0015182579855022 0 WAVE1 -0.21573979911255 -0.049063841960712 diff --git a/tests/test_all_component_and_model_builder.py b/tests/test_all_component_and_model_builder.py index 8dab2f33a..893b7bfd9 100644 --- a/tests/test_all_component_and_model_builder.py +++ b/tests/test_all_component_and_model_builder.py @@ -448,7 +448,10 @@ def __init__(self): del Component.component_types["SubsetModel2"] -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", +] # Test all the parameters. diff --git a/tests/test_clockcorr.py b/tests/test_clockcorr.py index a3871cd0d..b232676c4 100644 --- a/tests/test_clockcorr.py +++ b/tests/test_clockcorr.py @@ -13,8 +13,8 @@ from pint.observatory import ( Observatory, get_observatory, - ClockCorrectionOutOfRange, ) +from pint.exceptions import ClockCorrectionOutOfRange from pint.observatory.clock_file import ClockFile from pint.toa import get_TOAs diff --git a/tests/test_cmx.py b/tests/test_cmx.py new file mode 100644 index 000000000..976581880 --- /dev/null +++ b/tests/test_cmx.py @@ -0,0 +1,132 @@ +from io import StringIO +from pint.models import get_model +from pint.simulation import make_fake_toas_uniform +from pint.fitter import Fitter + +import pytest +import numpy as np +import astropy.units as u + + +@pytest.fixture +def model_and_toas(): + par = """ + RAJ 05:00:00 + DECJ 12:00:00 + F0 101.0 + F1 -1.1e-14 + PEPOCH 55000 + DM 12 + DMEPOCH 55000 + CMEPOCH 55000 + CM 3.5 + TNCHROMIDX 4.0 + EPHEM DE440 + CLOCK TT(BIPM2021) + UNITS TDB + TZRMJD 55000 + TZRFRQ 1400 + TZRSITE gmrt + CMXR1_0001 53999.9 + CMXR2_0001 54500 + CMX_0001 0.5 1 + CMXR1_0002 54500.1 + CMXR2_0002 55000 + CMX_0002 -0.5 1 + CMXR1_0003 55000.1 + CMXR2_0003 55500 + CMX_0003 -0.1 1 + """ + + model = get_model(StringIO(par)) + + freqs = np.linspace(300, 1600, 16) * u.MHz + toas = make_fake_toas_uniform( + startMJD=54000, + endMJD=56000, + ntoas=2000, + model=model, + freq=freqs, + add_noise=True, + obs="gmrt", + include_bipm=True, + multi_freqs_in_epoch=True, + ) + + return model, toas + + +def test_cmx(model_and_toas): + model, toas = model_and_toas + + assert "ChromaticCMX" in model.components + + ftr = Fitter.auto(toas, model) + ftr.fit_toas() + + assert np.abs(model.CMX_0001.value - ftr.model.CMX_0001.value) / ( + 3 * ftr.model.CMX_0001.uncertainty_value + ) + assert np.abs(model.CMX_0002.value - ftr.model.CMX_0002.value) / ( + 3 * ftr.model.CMX_0002.uncertainty_value + ) + assert np.abs(model.CMX_0003.value - ftr.model.CMX_0003.value) / ( + 3 * ftr.model.CMX_0003.uncertainty_value + ) + + assert ftr.resids.chi2_reduced < 1.6 + + assert "CMX_0001" in str(ftr.model) + + +def test_cmx_delay(model_and_toas): + model, toas = model_and_toas + + # Zero delay outside CMX ranges + nocmx_mask = toas.get_mjds().value > 55500 + assert all( + model.components["ChromaticCMX"].CMX_chromatic_delay(toas)[nocmx_mask] == 0 + ) + + # The delay is consistent + cmx1_mask = np.logical_and( + toas.get_mjds().value >= model.CMXR1_0001.value, + toas.get_mjds().value <= model.CMXR2_0001.value, + ) + cmx1_freqs = toas.get_freqs()[cmx1_mask] + assert all( + np.isclose( + model.components["ChromaticCMX"].chromatic_time_delay( + model.CMX_0001.quantity, model.TNCHROMIDX.quantity, cmx1_freqs + ), + model.components["ChromaticCMX"].CMX_chromatic_delay(toas)[cmx1_mask], + ) + ) + + cmx2_mask = np.logical_and( + toas.get_mjds().value >= model.CMXR1_0002.value, + toas.get_mjds().value <= model.CMXR2_0002.value, + ) + cmx2_freqs = toas.get_freqs()[cmx2_mask] + assert all( + np.isclose( + model.components["ChromaticCMX"].chromatic_time_delay( + model.CMX_0002.quantity, model.TNCHROMIDX.quantity, cmx2_freqs + ), + model.components["ChromaticCMX"].CMX_chromatic_delay(toas)[cmx2_mask], + ) + ) + + cmx3_mask = np.logical_and( + toas.get_mjds().value >= model.CMXR1_0003.value, + toas.get_mjds().value <= model.CMXR2_0003.value, + ) + cmx3_freqs = toas.get_freqs()[cmx3_mask] + assert all( + np.isclose( + model.components["ChromaticCMX"].chromatic_time_delay( + model.CMX_0003.quantity, model.TNCHROMIDX.quantity, cmx3_freqs + ), + model.components["ChromaticCMX"].CMX_chromatic_delay(toas)[cmx3_mask], + ) + ) diff --git a/tests/test_derived_quantities.py b/tests/test_derived_quantities.py index cf5e845b0..833e98a86 100644 --- a/tests/test_derived_quantities.py +++ b/tests/test_derived_quantities.py @@ -67,7 +67,7 @@ def test_Edot(): def test_Bfield(): # B assert np.isclose( - pulsar_B(0.033 * u.Hz, -2.0e-15 * u.Hz / u.s), 238722891596281.66 * u.G + pulsar_B(0.033 * u.Hz, -2.0e-15 * u.Hz / u.s), 238693670891966.22 * u.G ) @@ -75,7 +75,7 @@ def test_Blc(): # B_lc assert np.isclose( pulsar_B_lightcyl(0.033 * u.Hz, -2.0e-15 * u.Hz / u.s), - 0.07774704753236616 * u.G, + 0.07896965114785195 * u.G, ) diff --git a/tests/test_glitch.py b/tests/test_glitch.py index 4d7eac8ad..150879120 100644 --- a/tests/test_glitch.py +++ b/tests/test_glitch.py @@ -1,3 +1,4 @@ +from io import StringIO import os import pytest @@ -5,8 +6,9 @@ import numpy as np import pytest +from pint.exceptions import MissingParameter import pint.fitter -import pint.models +from pint.models import get_model import pint.residuals import pint.toa from pinttestdata import datadir @@ -14,11 +16,54 @@ parfile = os.path.join(datadir, "prefixtest.par") timfile = os.path.join(datadir, "prefixtest.tim") +basepar = """ + PSRJ J0835-4510 + RAJ 08:35:20.61149 + DECJ -45:10:34.8751 + F0 11.18965156782 + PEPOCH 55305 + DM 67.99 + UNITS TDB +""" + +good = """ + GLEP_1 55555 + GLPH_1 0 + GLF0_1 1.0e-6 + GLF1_1 -1.0e-12 + GLF0D_1 1.0e-6 + GLTD_1 10.0 +""" + +# no exponential decay set +bad1 = """ + GLEP_1 55555 + GLF0_1 1.0e-6 + GLF1_1 -1.0e-12 + GLF0D_1 1.0e-6 + GLTD_1 0.0 +""" + +# no epoch set +bad2 = """ + GLPH_1 0 + GLF0_1 1.0e-6 + GLF1_1 -1.0e-12 +""" + +# fitting both epoch and glitch phase +bad3 = """ + GLEP_1 55555 1 0 + GLPH_1 0 1 0 + GLF0_1 1.0e-6 + GLF1_1 -1.0e-12 +""" + class TestGlitch: @classmethod def setup_class(cls): - cls.m = pint.models.get_model(parfile) + cls.m = get_model(parfile) cls.t = pint.toa.get_TOAs(timfile, ephem="DE405", include_bipm=False) cls.f = pint.fitter.WLSFitter(cls.t, cls.m) @@ -37,11 +82,11 @@ def test_glitch_der(self): for pf in self.m.glitch_prop: for idx in set(self.m.glitch_indices): if pf in ["GLF0D_", "GLTD_"]: - getattr(self.m, "GLF0D_%d" % idx).value = 1.0 - getattr(self.m, "GLTD_%d" % idx).value = 100 + getattr(self.m, f"GLF0D_{idx}").value = 1.0 + getattr(self.m, f"GLTD_{idx}").value = 100 else: - getattr(self.m, "GLF0D_%d" % idx).value = 0.0 - getattr(self.m, "GLTD_%d" % idx).value = 0.0 + getattr(self.m, f"GLF0D_{idx}").value = 0.0 + getattr(self.m, f"GLTD_{idx}").value = 0.0 param = pf + str(idx) adf = self.m.d_phase_d_param(self.t, delay, param) param_obj = getattr(self.m, param) @@ -54,3 +99,13 @@ def test_glitch_der(self): errormsg = f"Derivatives for {param} is not accurate, max relative difference is" errormsg += " %lf" % np.nanmax(np.abs(r_diff.value)) assert np.nanmax(np.abs(r_diff.value)) < 1e-3, errormsg + + def test_bad_input(self): + get_model(StringIO(basepar + good)) + with pytest.raises(MissingParameter): + get_model(StringIO(basepar + bad1)) + with pytest.raises(MissingParameter): + m = get_model(StringIO(basepar + bad2)) + print(m) + with pytest.raises(ValueError): + get_model(StringIO(basepar + bad3)) diff --git a/tests/test_model_derivatives.py b/tests/test_model_derivatives.py index a3bcb28a4..74737eb78 100644 --- a/tests/test_model_derivatives.py +++ b/tests/test_model_derivatives.py @@ -153,7 +153,7 @@ def f(value): try: dphase = m.phase(toas, abs_phase=False) - phase except ValueError: - return np.nan * np.zeros_like(phase.frac) + return np.nan * np.zeros_like(phase.frac).astype(np.float64) return np.float64(dphase.int + dphase.frac) if param == "ECC": @@ -178,6 +178,9 @@ def f(value): stepgen = numdifftools.MaxStepGenerator( np.abs(model.FB3.value.astype(np.float64)) * 1e7 ) + elif param == "SINI": + sini = model.SINI.value.astype(np.float64) + stepgen = numdifftools.MaxStepGenerator(abs(sini) / 10) else: stepgen = None df = numdifftools.Derivative(f, step=stepgen) diff --git a/tests/test_model_wave.py b/tests/test_model_wave.py index 1e57d640f..ab1707088 100644 --- a/tests/test_model_wave.py +++ b/tests/test_model_wave.py @@ -1,28 +1,81 @@ +from io import StringIO import os import pytest import astropy.units as u +import numpy as np -import pint.fitter -import pint.models +from pint.models import get_model +from pint import toa import pint.residuals -import pint.toa from pinttestdata import datadir -# Not included in the test here, but as a sanity check I used this same -# ephemeris to phase up Fermi data, and it looks good. +par_nowave = """ + PSRJ J0835-4510 + RAJ 08:35:20.61149 + DECJ -45:10:34.8751 + F0 11.18965156782 + PEPOCH 55305 + DM 67.99 + UNITS TDB +""" -parfile = os.path.join(datadir, "vela_wave.par") -timfile = os.path.join(datadir, "vela_wave.tim") +wave_terms = """ + WAVEEPOCH 55305 + WAVE_OM 0.0015182579855022 + WAVE1 -0.21573979911255 -0.049063841960712 + WAVE2 0.62795320246729 -0.11984954655989 + WAVE3 0.099618608456845 0.28530756908788 + WAVE4 -0.21537340649058 0.18849486610628 + WAVE5 0.021980474493165 -0.23566696662127 +""" -class TestWave: - @classmethod - def setup_class(cls): - cls.m = pint.models.get_model(parfile) - cls.t = pint.toa.get_TOAs(timfile, ephem="DE405", include_bipm=False) +def test_wave_ephem(): + parfile = os.path.join(datadir, "vela_wave.par") + timfile = os.path.join(datadir, "vela_wave.tim") + m = get_model(parfile) + t = toa.get_TOAs(timfile, ephem="DE405", include_bipm=False) + rs = pint.residuals.Residuals(t, m).time_resids + assert rs.std() < 350.0 * u.us - def test_vela_rms_is_small_enough(self): - rs = pint.residuals.Residuals(self.t, self.m).time_resids - rms = rs.to(u.us).std() - assert rms < 350.0 * u.us + +def test_wave_construction(): + m = get_model(StringIO(par_nowave + wave_terms)) + assert np.allclose(m.WAVE_OM.quantity, 0.0015182579855022 * u.rad / u.day) + assert np.allclose(m.WAVE1.quantity[0], -0.21573979911255 * u.s) + assert np.allclose(m.WAVE2.quantity[1], -0.1198495465598 * u.s) + + +def test_wave_computation(): + m0 = get_model(StringIO(par_nowave)) + m1 = get_model(StringIO(par_nowave + wave_terms)) + # make some barycentric TOAs + tdbs = np.linspace(54500, 60000, 10) + ts = toa.TOAs( + toalist=[ + toa.TOA(t, obs="@", freq=np.inf, error=1 * u.ms, scale="tdb") for t in tdbs + ] + ) + ts.compute_TDBs(ephem="DE421") + ts.compute_posvels(ephem="DE421") + ph0 = m0.phase(ts) + ph1 = m1.phase(ts) + dphi = (ph1.int - ph0.int) + (ph1.frac - ph0.frac) + test_phase = np.zeros(len(tdbs)) + WAVEEPOCH = 55305 + WAVE_OM = 0.0015182579855022 + WAVE = [ + [-0.21573979911255, -0.049063841960712], + [0.62795320246729, -0.11984954655989], + [0.099618608456845, 0.28530756908788], + [-0.21537340649058, 0.18849486610628], + [0.021980474493165, -0.23566696662127], + ] + ph = (tdbs - WAVEEPOCH) * WAVE_OM + for i in range(5): + test_phase += WAVE[i][0] * np.sin((i + 1) * ph) + WAVE[i][1] * np.cos( + (i + 1) * ph + ) + test_phase *= m0.F0.quantity.to(u.Hz).value + assert np.allclose(test_phase, dphi) diff --git a/tests/test_orbwaves.py b/tests/test_orbwaves.py new file mode 100644 index 000000000..06ce0141d --- /dev/null +++ b/tests/test_orbwaves.py @@ -0,0 +1,70 @@ +from collections import deque +from io import StringIO +import os +import pytest + +import astropy.units as u +import pint.models as tm +from pint import fitter, toa, simulation +from pinttestdata import datadir +import pint.models.parameter as param +from pint import ls +from pint.models import get_model, get_model_and_toas + + +def test_orbwaves_fit(): + m, t = get_model_and_toas( + datadir / "J1048+2339_orbwaves.par", datadir / "J1048+2339_3PC_fake.tim" + ) + + f = fitter.WLSFitter(toas=t, model=m) + + rms = f.resids.rms_weighted().to_value("us") + + assert rms > 29.2 and rms < 29.3 + + f.fit_toas() + rms = f.resids.rms_weighted().to_value("us") + + assert rms > 21.1 and rms < 21.2 + + +def test_orbwaves_DD_fit(): + m, t = get_model_and_toas( + datadir / "J1048+2339_orbwaves_DD.par", datadir / "J1048+2339_3PC_fake.tim" + ) + + f = fitter.WLSFitter(toas=t, model=m) + + rms = f.resids.rms_weighted().to_value("us") + + assert rms > 29.2 and rms < 29.3 + + f.fit_toas() + rms = f.resids.rms_weighted().to_value("us") + + assert rms > 21.1 and rms < 21.2 + + +def test_invalid_parfiles(): + parlines = open(datadir / "J1048+2339_orbwaves.par").readlines() + + def _delete_line(labels): + lines = deque() + for line in parlines: + label = line.split()[0] + if label in labels: + continue + lines.append(line) + return StringIO(lines) + + deleted_lines = [["ORBWAVEC0", "ORBWAVES0"], ["ORBWAVEC3"], ["ORBWAVES4"]] + + for lines in deleted_lines: + with pytest.raises(Exception): + m = get_model(_delete_line(lines)) + + +if __name__ == "__main__": + test_orbwaves_fit() + test_invalid_parfiles() diff --git a/tests/test_phase_offset.py b/tests/test_phase_offset.py index d28fc26c0..61d46bd5f 100644 --- a/tests/test_phase_offset.py +++ b/tests/test_phase_offset.py @@ -26,7 +26,7 @@ def test_phase_offset(): m = get_model(io.StringIO(simplepar)) assert hasattr(m, "PHOFF") and m.PHOFF.value == 0.2 - + np.random.seed(1929) t = make_fake_toas_uniform( startMJD=50000, endMJD=50500, diff --git a/tests/test_process_parfile.py b/tests/test_process_parfile.py index daaa23d20..da3dbf9e7 100644 --- a/tests/test_process_parfile.py +++ b/tests/test_process_parfile.py @@ -4,9 +4,10 @@ import pytest from io import StringIO -from pint.models.model_builder import ModelBuilder, parse_parfile, TimingModelError +from pint.models.model_builder import ModelBuilder, parse_parfile from pint.models.model_builder import guess_binary_model +from pint.exceptions import TimingModelError from pint.models import get_model base_par = """ diff --git a/tests/test_sini.py b/tests/test_sini.py new file mode 100644 index 000000000..7f6f2b9a6 --- /dev/null +++ b/tests/test_sini.py @@ -0,0 +1,69 @@ +import pytest +import io +from astropy import units as u +import pint.logging +from pint.models import get_model +import pint.fitter +import pint.simulation + +pint.logging.setup("ERROR") + + +def test_sini_fit(): + s = """PSRJ J1853+1303 + EPHEM DE440 + CLK TT(TAI) + UNITS TDB + START 55731.1937186050359860 + FINISH 59051.1334073910116910 + TIMEEPH FB90 + T2CMETHOD IAU2000B + DILATEFREQ N + DMDATA N + NTOA 4570 + CHI2R 0.9873 0 4455.0 + TRES 1.297 + ELONG 286.257303932061347 1 0.00000000931838521883 + ELAT 35.743349195626379 1 0.00000001342269042838 + PMELONG -1.9106838589844408 1 0.011774688289719207 + PMELAT -2.7055634767966588 1 0.02499026250312041 + PX 0.6166118941596491 1 0.13495439191146907 + ECL IERS2010 + POSEPOCH 57391.0000000000000000 + F0 244.39137771284777388 1 4.741681125242900587e-13 + F1 -5.2068158193285555695e-16 1 1.657252200776189767e-20 + PEPOCH 57391.0000000000000000 + CORRECT_TROPOSPHERE Y + PLANET_SHAPIRO Y + NE_SW 0.0 + SWM 0.0 + DM 30.570200097819536894 + DM1 0.0 + DMEPOCH 57391.0000000000000000 + BINARY DD + PB 115.65378643873621578 1 3.2368969907e-09 + PBDOT 0.0 + A1 40.769522249658145 1 1.53854961349644e-06 + XDOT 1.564794972862985e-14 1 8.617009015429654e-16 + ECC 2.3701427108198803e-05 1 3.467883031e-09 + EDOT 0.0 + T0 57400.7367807004642730 1 0.00999271094974152688 + OM 346.60040313716600618 1 0.03110468777178034688 + OMDOT 0.0 + M2 0.18373369620702712 1 0.19097347492091507 + SINI 0.8630519453652643 1 0.14457246899831822 + A0 0.0 + B0 0.0 + GAMMA 0.0 + DR 0.0 + DTH 0.0 + TZRMJD 57390.6587742196275694 + TZRSITE arecibo + TZRFRQ 426.875 + """ + m = get_model(io.StringIO(s)) + t = pint.simulation.make_fake_toas_uniform( + 55731, 59051, 100, m, obs="axis", error=1 * u.us, add_noise=True + ) + f = pint.fitter.Fitter.auto(t, m) + f.fit_toas() diff --git a/tests/test_solar_wind.py b/tests/test_solar_wind.py index bfa8e0eea..69a79fd52 100644 --- a/tests/test_solar_wind.py +++ b/tests/test_solar_wind.py @@ -359,3 +359,83 @@ def test_overlapping_swx(): dm3 = model3.components["SolarWindDispersionX"].swx_dm(toas) dm3[toas.get_mjds() >= 54180 * u.d] = 0 assert np.allclose(dm2 - dm, dm3) + + +def test_nesw_derivatives(): + par = """ + PSRJ J1744-1134 + ELONG 266.119458498 6.000e-09 + ELAT 11.80517508 3.000e-08 + DM 3.1379 4.000e-04 + PEPOCH 55000 + F0 245.4261196898081 5.000e-13 + F1 -5.38156E-16 3.000e-21 + POSEPOCH 59150 + DMEPOCH 55000 + CLK TT(BIPM2018) + EPHEM DE436 + RM 1.62 9.000e-02 + PX 2.61 9.000e-02 + DM1 -7.0E-5 3.000e-05 + PMELONG 19.11 2.000e-02 + DM2 1.7E-5 4.000e-06 + PMELAT -9.21 1.300e-01 + UNITS TDB + NE_SW 8 1 + NE_SW1 0 1 + NE_SW2 0 1 + SWEPOCH 55000 + """ + m = get_model(StringIO(par)) + t = make_fake_toas_uniform(54500, 55500, 1000, m, add_noise=True) + ftr = Fitter.auto(t, m) + ftr.fit_toas() + + assert ( + m.NE_SW.value - ftr.model.NE_SW.value + ) / ftr.model.NE_SW.uncertainty_value < 3 + assert ( + m.NE_SW1.value - ftr.model.NE_SW1.value + ) / ftr.model.NE_SW1.uncertainty_value < 3 + assert ( + m.NE_SW2.value - ftr.model.NE_SW2.value + ) / ftr.model.NE_SW2.uncertainty_value < 3 + + m.components["SolarWindDispersion"] + + +def test_expression(): + par = """ + PSRJ J1744-1134 + ELONG 266.119458498 6.000e-09 + ELAT 11.80517508 3.000e-08 + DM 3.1379 4.000e-04 + PEPOCH 55000 + F0 245.4261196898081 5.000e-13 + F1 -5.38156E-16 3.000e-21 + POSEPOCH 59150 + DMEPOCH 55000 + CLK TT(BIPM2018) + EPHEM DE436 + PX 2.61 9.000e-02 + PMELONG 19.11 2.000e-02 + PMELAT -9.21 1.300e-01 + UNITS TDB + NE_SW 8 1 + NE_SW1 1 1 + SWEPOCH 55000 + """ + m = get_model(StringIO(par)) + t = make_fake_toas_uniform(54500, 55500, 10, m, add_noise=True) + + t0 = m.SWEPOCH.value * u.day + t1 = t["tdbld"][-1] * u.day + dt = t1 - t0 + + assert np.isclose( + ( + m.components["SolarWindDispersion"].solar_wind_dm(t) + / m.components["SolarWindDispersion"].solar_wind_geometry(t) + )[-1], + (m.NE_SW.quantity + dt * m.NE_SW1.quantity), + ) diff --git a/tests/test_variety_parfiles.py b/tests/test_variety_parfiles.py index c137d59a1..bac88e9ac 100644 --- a/tests/test_variety_parfiles.py +++ b/tests/test_variety_parfiles.py @@ -3,7 +3,7 @@ import pytest from io import StringIO -from pint.models.timing_model import ( +from pint.exceptions import ( TimingModelError, UnknownBinaryModel, MissingBinaryError, diff --git a/tox.ini b/tox.ini index 87a498e01..4e19b0fda 100644 --- a/tox.ini +++ b/tox.ini @@ -13,7 +13,8 @@ envlist = report codestyle black - py{38,39,310,311,312}-test{,-alldeps,-devdeps}{,-cov} + singletest + py{38,39,310,311,312,313}-test{,-alldeps,-devdeps}{,-cov} skip_missing_interpreters = True @@ -36,21 +37,39 @@ deps = cov: coverage cov: pytest-cov cov: pytest-remotedata + pytest-rerunfailures hypothesis numdifftools pathos setuptools commands = pip freeze - !cov: pytest - cov: pytest -v --pyargs tests --cov=pint --cov-config={toxinidir}/.coveragerc {posargs} + !cov: pytest --reruns 5 + cov: pytest -v --pyargs tests --cov=pint --cov-config={toxinidir}/.coveragerc {posargs} --reruns 5 cov: coverage xml -o {toxinidir}/coverage.xml depends = - {py38,py39,py310,py311,py312}: clean - report: py38,py39,py310,py312 + {py39,py310,py311,py312,313}: clean + report: py39,py310,py312,py313 docs: notebooks +[testenv:singletest] +description = + Try a simple run with a single test +deps = + numpy + numdifftools + astropy + matplotlib + scipy + pytest + coverage + hypothesis<=6.72.0 + setuptools +# can change this as needed for a single test run +commands = pytest tests/test_precision.py + + [testenv:ephemeris_connection] description = Check whether PINT can obtain the DE440 ephemeris (usually from a server) @@ -62,17 +81,19 @@ deps = [testenv:oldestdeps] description = Run tests on Python 3 with minimum supported versions of astropy, numpy -basepython = python3.8 +basepython = python3.9 deps = - numpy==1.18.5 + numpy==1.23.0 numdifftools==0.9.39 - astropy==4.0 - matplotlib==3.2.0 - scipy==1.4.1 + astropy==5.0.5 + matplotlib==3.4.3 + scipy==1.9.0 pytest + pytest-rerunfailures coverage hypothesis<=6.72.0 -commands = {posargs:pytest} +commands = + pytest --reruns 5 [testenv:report] skip_install = true @@ -84,7 +105,7 @@ commands = [testenv:notebooks] description = update the notebooks -basepython = python3.12 +basepython = python3.13 deps = traitlets sphinx >= 2.2, != 5.1.0 @@ -112,7 +133,7 @@ commands = coverage erase [testenv:docs] changedir = {toxinidir}/docs description = invoke sphinx-build to build the HTML docs -basepython = python3.12 +basepython = python3.13 deps = traitlets sphinx >= 2.2, != 5.1.0