Merge branch 'nanograv:master' into chromatic

.gitignore

@@ -119,4 +119,10 @@ docs/examples/*.ipynb
 docs/examples-rendered/*.py
 
 # VSCode wants to put virtualenvs here
-.env
+.env
+
+# pintpublish output
+*.tex
+*.aux
+*.log
+*.pdf

.readthedocs.yml

@@ -17,7 +17,8 @@ build:
 
 # Build documentation in the docs/ directory with Sphinx
 sphinx:
-   configuration: docs/conf.py
+  builder: html
+  configuration: docs/conf.py
 
 # If using Sphinx, optionally build your docs in additional formats such as PDF
 formats:
@@ -29,4 +30,4 @@ python:
     - requirements: requirements_dev.txt
     - requirements: requirements.txt
     - method: pip
-      path: .
+      path: .

AUTHORS.rst

@@ -16,6 +16,7 @@ Contributors
 
 Active developers are indicated by (*). Authors of the PINT paper are indicated by (#).
 
+* Gabriella Agazie (*)
 * Akash Anumarlapudi (*)
 * Anne Archibald (#*)
 * Matteo Bachetti (#)

CHANGELOG-unreleased.md

@@ -9,7 +9,21 @@ the released changes.
 
 ## Unreleased
 ### Changed
+- `WAVE` parameters can be added to a `Wave` model with `add_wave_component()` in `wave.py` 
+- Moved design matrix normalization code from `pint.fitter` to the new `pint.utils.normalize_designmatrix()` function.
+- Made `Residuals` independent of `GLSFitter` (GLS chi2 is now computed using the new function `Residuals._calc_gls_chi2()`).
 ### Added
+- Added `WaveX` model as a `DelayComponent` with Fourier amplitudes as fitted parameters
+- `Parameter.as_latex` method for latex representation of a parameter.
+- `pint.output.publish` module and `pintpublish` script for generating publication (LaTeX) output.
+- Added radial velocity methods for binary models
+- Added `DMWaveX` model (Fourier representation of DM noise)
 ### Fixed
+- Wave model `validate()` can correctly use PEPOCH to assign WAVEEPOCH parameter
 - Fixed RTD by specifying theme explicitly.
+- `.value()` now works for pairParameters
+- Setting `model.PARAM1 = model.PARAM2` no longer overrides the name of `PARAM1`
+- Fixed an incorrect docstring in `pbprime()` functions. 
+- Fix ICRS -> ECL conversion when parameter uncertainties are not set.
+- `get_TOAs` raises an exception upon finding mixed narrowband and wideband TOAs in a tim file. `TOAs.is_wideband` returns True only if *ALL* TOAs have the -pp_dm flag.
 ### Removed

docs/examples/understanding_fitters.py

@@ -144,7 +144,7 @@
 # straightforward but the full covariance matrix may be enormous.
 # If False, an algorithm is used that takes advantage of the structure
 # of the covariance matrix, based on information provided by the noise
-# model. The two algorithms should give the same result to numerical
+# model. The two algorithms should give the same result up to numerical
 # accuracy where they both can be applied.
 
 # %% [markdown]
@@ -178,11 +178,6 @@
 # %%
 glsfit.fit_toas(maxiter=1)
 
-# %%
-# Not sure how to do this properly yet.
-# glsfit2 = pint.fitter.GLSFitter(toas=t, model=glsfit.model, residuals=glsfit.resids)
-# glsfit2.fit_toas(maxiter=0)
-
 # %%
 glsfit.print_summary()
 
@@ -206,7 +201,8 @@
 # %% [markdown]
 # ## Choosing fitters
 #
-# You can use the automatic fitter selection to help you choose between `WLSFitter`, `GLSFitter`, and their wideband variants.  The default `Downhill` fitters generally have better performance than the plain variants.
+# You can use the automatic fitter selection to help you choose between `WLSFitter`, `GLSFitter`, and their wideband variants.
+# The default `Downhill` fitters generally have better performance than the plain variants.
 
 # %%
 autofit = pint.fitter.Fitter.auto(toas=ts1855, model=m1855)
@@ -221,4 +217,5 @@
 # The results are (thankfully) identical.
 
 # %% [markdown]
-# The MCMC fitter is considerably more complicated, so it has its own dedicated walkthroughs in `MCMC_walkthrough.ipynb` (for photon data) and `examples/fit_NGC6440E_MCMC.py` (for fitting TOAs).
+# The MCMC fitter is considerably more complicated, so it has its own dedicated walkthroughs in `MCMC_walkthrough.ipynb`
+# (for photon data) and `examples/fit_NGC6440E_MCMC.py` (for fitting TOAs).

docs/explanation.rst

@@ -716,13 +716,16 @@ model and data require different calculations - narrowband (TOA-only) versus
 wideband (TOA and DM measurements) and uncorrelated errors versus correlated
 errors.
 
-The TEMPO/TEMPO2 and default PINT fitting algorithms (:class:`pint.fitter.WidebandTOAFitter` for example), leaving aside the rank-reduced case, proceed like:
+The TEMPO/TEMPO2 and default PINT fitting algorithms (:class:`pint.fitter.WidebandTOAFitter`, for example), 
+leaving aside the rank-reduced case, proceed like:
 
 1. Evaluate the model and its derivatives at the starting point :math:`x`, producing a set of residuals :math:`\delta y` and a Jacobian `M`.
 2. Compute :math:`\delta x` to minimize :math:`\left| M\delta x - \delta y \right|_C`, where :math:`\left| \cdot \right|_C` is the squared amplitude of a vector with respect to the data uncertainties/covariance :math:`C`.
 3. Update the starting point by :math:`\delta x`.
 
-TEMPO and TEMPO2 can check whether the predicted improvement of chi-squared, assuming the linear model is correct, is enough to warrant continuing; if so, they jump back to step 1 unless the maximum number of iterations is reached. PINT does not contain this check.
+TEMPO and TEMPO2 can check whether the predicted improvement of chi-squared, assuming 
+the linear model is correct, is enough to warrant continuing; if so, they jump back to 
+step 1 unless the maximum number of iterations is reached. PINT does not contain this check.
 
 This algorithm is the Gauss-Newton_algorithm_ for solving nonlinear
 least-squares problems, and even in one-complex-dimensional cases can exhibit
@@ -744,9 +747,16 @@ PINT contains a slightly more sophisticated algorithm, implemented in
 4. Evaluate the model at the starting point plus :math:`\lambda \delta x`. If this is invalid or worse than the starting point, divide :math:`\lambda` by two and repeat this step. If :math:`\lambda` is too small, accept the best point seen to date and exit without convergence.
 5. If the model improved but only slightly with :math:`\lambda=1`, exit with convergence. If the maximum number of iterations was reached, exit without convergence. Otherwise update the starting point and return to step 1.
 
-This ensures that PINT tries taking smaller steps if problems arise, and claims convergence only if a normal step worked. It does not solve the problems that arise if some parameters are nearly degenerate, enough to cause problems with the numerical linear algebra.
+This ensures that PINT tries taking smaller steps if problems arise, and claims convergence 
+only if a normal step worked. It does not solve the problems that arise if some parameters are 
+nearly degenerate, enough to cause problems with the numerical linear algebra.
 
-As a rule, this kind of problem is addressed with the Levenberg-Marquardt algorithm, which operates on the same principle of taking reduced steps when the derivative appears not to match the function, but does so in a way that also reduces issues with degenerate parameters; unfortunately it is not clear how to adapt this problem to the rank-reduced case. Nevertheless PINT contains an implementation, in :class:`pint.fitter.WidebandLMFitter`, but it does not perform as well as one might hope in practice and must be considered experimental.
+As a rule, this kind of problem is addressed with the Levenberg-Marquardt algorithm, which 
+operates on the same principle of taking reduced steps when the derivative appears not to 
+match the function, but does so in a way that also reduces issues with degenerate parameters; 
+unfortunately it is not clear how to adapt this problem to the rank-reduced case. Nevertheless, 
+PINT contains an implementation in :class:`pint.fitter.WidebandLMFitter`, but it does not perform as 
+well as one might hope in practice and must be considered experimental.
 
 .. _Gauss-Newton_algorithm: https://en.wikipedia.org/wiki/Gauss%E2%80%93Newton_algorithm
 .. _chaotic: https://en.wikipedia.org/wiki/Newton_fractal

setup.cfg

@@ -61,6 +61,7 @@ console_scripts =
     convert_parfile = pint.scripts.convert_parfile:main
     compare_parfiles = pint.scripts.compare_parfiles:main
     tcb2tdb = pint.scripts.tcb2tdb:main
+    pintpublish = pint.scripts.pintpublish:main
 
 
 # See the docstring in versioneer.py for instructions. Note that you must

src/pint/fitter.py

@@ -83,12 +83,11 @@
 )
 from pint.residuals import Residuals, WidebandTOAResiduals
 from pint.toa import TOAs
-from pint.utils import FTest
+from pint.utils import FTest, normalize_designmatrix
 
 
 __all__ = [
     "Fitter",
-    "auto",
     "WLSFitter",
     "GLSFitter",
     "WidebandTOAFitter",
@@ -1303,9 +1302,7 @@ def step(self):
         # NOTE, We remove subtract mean value here, since it did not give us a
         # fast converge fitting.
         # M[:,1:] -= M[:,1:].mean(axis=0)
-        fac = np.sqrt((M**2).mean(axis=0))
-        fac[fac == 0] = 1.0
-        M /= fac
+        M, fac = normalize_designmatrix(M, params)
         # Singular value decomp of design matrix:
         #   M = U s V^T
         # Dimensions:
@@ -1457,15 +1454,7 @@ def step(self):
                 M = np.hstack((M, Mn))
 
         # normalize the design matrix
-        norm = np.sqrt(np.sum(M**2, axis=0))
-        for c in np.where(norm == 0)[0]:
-            warn(
-                f"Parameter degeneracy; the following parameter yields "
-                f"almost no change: {params[c]}",
-                DegeneracyWarning,
-            )
-        norm[norm == 0] = 1
-        M /= norm
+        M, norm = normalize_designmatrix(M, params)
         self.M = M
         self.fac = norm
 
@@ -2014,10 +2003,7 @@ def fit_toas(self, maxiter=1, threshold=None, debug=False):
             # NOTE, We remove subtract mean value here, since it did not give us a
             # fast converge fitting.
             # M[:,1:] -= M[:,1:].mean(axis=0)
-            fac = np.sqrt((M**2).mean(axis=0))
-            # fac[0] = 1.0
-            fac[fac == 0] = 1.0
-            M /= fac
+            M, fac = normalize_designmatrix(M, params)
             # Singular value decomp of design matrix:
             #   M = U s V^T
             # Dimensions:
@@ -2165,7 +2151,11 @@ def fit_toas(self, maxiter=1, threshold=0, full_cov=False, debug=False):
             fitperrs = self.model.get_params_dict("free", "uncertainty")
 
             # Define the linear system
+            # normalize the design matrix
             M, params, units = self.get_designmatrix()
+            # M /= norm
+
+            ntmpar = len(fitp)
 
             # Get residuals and TOA uncertainties in seconds
             if i == 0:
@@ -2182,18 +2172,11 @@ def fit_toas(self, maxiter=1, threshold=0, full_cov=False, debug=False):
                     phiinv = np.concatenate((phiinv, 1 / phi))
                     M = np.hstack((M, Mn))
 
-            # normalize the design matrix
-            norm = np.sqrt(np.sum(M**2, axis=0))
             ntmpar = len(fitp)
-            for c in np.where(norm == 0)[0]:
-                warn(
-                    f"Parameter degeneracy; the following parameter yields "
-                    f"almost no change: {params[c]}",
-                    DegeneracyWarning,
-                )
-            norm[norm == 0] = 1
+
+            # normalize the design matrix
+            M, norm = normalize_designmatrix(M, params)
             self.fac = norm
-            M /= norm
 
             # compute covariance matrices
             if full_cov:
@@ -2537,17 +2520,10 @@ def fit_toas(self, maxiter=1, threshold=0, full_cov=False, debug=False):
                         new_d_matrix.param_units,
                     )
 
-            # normalize the design matrix
-            norm = np.sqrt(np.sum(M**2, axis=0))
             ntmpar = len(fitp)
-            for c in np.where(norm == 0)[0]:
-                warn(
-                    f"Parameter degeneracy; the following parameter yields "
-                    f"almost no change: {params[c]}",
-                    DegeneracyWarning,
-                )
-            norm[norm == 0] = 1
-            M /= norm
+
+            # normalize the design matrix
+            M, norm = normalize_designmatrix(M, params)
             self.fac = norm
 
             # compute covariance matrices

src/pint/models/__init__.py

@@ -26,6 +26,7 @@
 from pint.models.binary_ell1 import BinaryELL1, BinaryELL1H, BinaryELL1k
 from pint.models.chromatic_model import ChromaticCM
 from pint.models.dispersion_model import DispersionDM, DispersionDMX
+from pint.models.dmwavex import DMWaveX
 from pint.models.frequency_dependent import FD
 from pint.models.glitch import Glitch
 from pint.models.phase_offset import PhaseOffset
@@ -43,6 +44,7 @@
 from pint.models.timing_model import DEFAULT_ORDER, TimingModel
 from pint.models.troposphere_delay import TroposphereDelay
 from pint.models.wave import Wave
+from pint.models.wavex import WaveX
 
 # Define a standard basic model
 StandardTimingModel = TimingModel(

src/pint/models/absolute_phase.py

@@ -30,19 +30,21 @@ def __init__(self):
         super().__init__()
         self.add_param(
             MJDParameter(
-                name="TZRMJD", description="Epoch of the zero phase.", time_scale="utc"
+                name="TZRMJD",
+                description="Epoch of the zero phase TOA.",
+                time_scale="utc",
             )
         )
         self.add_param(
             strParameter(
-                name="TZRSITE", description="Observatory of the zero phase measured."
+                name="TZRSITE", description="Observatory of the zero phase TOA."
             )
         )
         self.add_param(
             floatParameter(
                 name="TZRFRQ",
                 units=u.MHz,
-                description="The frequency of the zero phase measured.",
+                description="The frequency of the zero phase TOA.",
             )
         )
         self.tz_cache = None

src/pint/models/astrometry.py

@@ -550,8 +550,16 @@ def as_ECL(self, epoch=None, ecl="IERS2010"):
             ra=self.RAJ.quantity,
             dec=self.DECJ.quantity,
             obstime=self.POSEPOCH.quantity,
-            pm_ra_cosdec=self.RAJ.uncertainty * np.cos(self.DECJ.quantity) / dt,
-            pm_dec=self.DECJ.uncertainty / dt,
+            pm_ra_cosdec=(
+                self.RAJ.uncertainty * np.cos(self.DECJ.quantity) / dt
+                if self.RAJ.uncertainty is not None
+                else 0 * self.RAJ.units / dt
+            ),
+            pm_dec=(
+                self.DECJ.uncertainty / dt
+                if self.DECJ.uncertainty is not None
+                else 0 * self.DECJ.units / dt
+            ),
             frame=coords.ICRS,
         )
         c_ECL = c.transform_to(PulsarEcliptic(ecl=ecl))
@@ -563,8 +571,16 @@ def as_ECL(self, epoch=None, ecl="IERS2010"):
             ra=self.RAJ.quantity,
             dec=self.DECJ.quantity,
             obstime=self.POSEPOCH.quantity,
-            pm_ra_cosdec=self.PMRA.uncertainty,
-            pm_dec=self.PMDEC.uncertainty,
+            pm_ra_cosdec=(
+                self.PMRA.uncertainty
+                if self.PMRA.uncertainty is not None
+                else 0 * self.PMRA.units
+            ),
+            pm_dec=(
+                self.PMDEC.uncertainty
+                if self.PMDEC.uncertainty is not None
+                else 0 * self.PMDEC.units
+            ),
             frame=coords.ICRS,
         )
         c_ECL = c.transform_to(PulsarEcliptic(ecl=ecl))
@@ -943,8 +959,16 @@ def as_ECL(self, epoch=None, ecl="IERS2010"):
             lat=self.ELAT.quantity,
             obliquity=OBL[self.ECL.value],
             obstime=self.POSEPOCH.quantity,
-            pm_lon_coslat=self.ELONG.uncertainty * np.cos(self.ELAT.quantity) / dt,
-            pm_lat=self.ELAT.uncertainty / dt,
+            pm_lon_coslat=(
+                self.ELONG.uncertainty * np.cos(self.ELAT.quantity) / dt
+                if self.ELONG.uncertainty is not None
+                else 0 * self.ELONG.units / dt
+            ),
+            pm_lat=(
+                self.ELAT.uncertainty / dt
+                if self.ELAT.uncertainty is not None
+                else 0 * self.ELAT.units / dt
+            ),
             frame=PulsarEcliptic,
         )
         c_ECL = c.transform_to(PulsarEcliptic(ecl=ecl))
@@ -957,12 +981,16 @@ def as_ECL(self, epoch=None, ecl="IERS2010"):
             lat=self.ELAT.quantity,
             obliquity=OBL[self.ECL.value],
             obstime=self.POSEPOCH.quantity,
-            pm_lon_coslat=self.PMELONG.uncertainty
-            if self.PMELONG.uncertainty is not None
-            else 0 * self.PMELONG.units,
-            pm_lat=self.PMELAT.uncertainty
-            if self.PMELAT.uncertainty is not None
-            else 0 * self.PMELAT.units,
+            pm_lon_coslat=(
+                self.PMELONG.uncertainty
+                if self.PMELONG.uncertainty is not None
+                else 0 * self.PMELONG.units
+            ),
+            pm_lat=(
+                self.PMELAT.uncertainty
+                if self.PMELAT.uncertainty is not None
+                else 0 * self.PMELAT.units
+            ),
             frame=PulsarEcliptic,
         )
         c_ECL = c.transform_to(PulsarEcliptic(ecl=ecl))

src/pint/models/binary_ell1.py

@@ -117,7 +117,7 @@ def __init__(self):
             floatParameter(
                 name="EPS1",
                 units="",
-                description="First Laplace-Lagrange parameter, ECC x sin(OM) for ELL1 model",
+                description="First Laplace-Lagrange parameter, ECC*sin(OM)",
                 long_double=True,
             )
         )
@@ -126,7 +126,7 @@ def __init__(self):
             floatParameter(
                 name="EPS2",
                 units="",
-                description="Second Laplace-Lagrange parameter, ECC x cos(OM) for ELL1 model",
+                description="Second Laplace-Lagrange parameter, ECC*cos(OM)",
                 long_double=True,
             )
         )