Brute force method to find the optimal number of harmonics for WaveX, DMWaveX, and CMWaveX

CHANGELOG-unreleased.md

@@ -14,13 +14,19 @@ the released changes.
 - `maskParameter.__repr__()` output now includes the frozen attribute.
 - Changed default value of `FDJUMPLOG` to `Y`
 - Bumped `black` version to 24.x
+- Replaced `pint.utils.find_optimal_nharms` by a more general function `pint.noise_analysis.find_optimal_nharms` which optimizes the `Nharms` for multiple noise components simultaneously. 
+- Updated the example `rednoise-fit-example.py`
 ### 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()`
 ### 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

docs/examples/rednoise-fit-example.py

@@ -18,12 +18,12 @@
 # %%
 from pint import DMconst
 from pint.models import get_model
+from pint.noise_analysis import find_optimal_nharms
 from pint.simulation import make_fake_toas_uniform
 from pint.logging import setup as setup_log
 from pint.fitter import WLSFitter
 from pint.utils import (
     dmwavex_setup,
-    find_optimal_nharms,
     wavex_setup,
     plrednoise_from_wavex,
     pldmnoise_from_dmwavex,
@@ -100,7 +100,11 @@
 m1 = deepcopy(m)
 m1.remove_component("PLRedNoise")
 
-nharm_opt, d_aics = find_optimal_nharms(m1, t, "WaveX", 30)
+aics, nharm_opt = find_optimal_nharms(
+    m1, t, include_components=["WaveX"], nharms_max=30
+)
+d_aics = aics - np.min(aics)
+nharm_opt = nharm_opt[0]
 
 print("Optimum no of harmonics = ", nharm_opt)
 
@@ -118,7 +122,7 @@
 plt.ylabel("AIC - AIC$_\\min{} + 1$")
 plt.legend()
 plt.yscale("log")
-# plt.savefig("sim3-aic.pdf")
+plt.show()
 
 # %%
 # Now create a new model with the optimum number of harmonics
@@ -181,7 +185,7 @@
 )
 plt.axvline(fyr, color="black", ls="dotted")
 plt.axhline(0, color="grey", ls="--")
-plt.ylabel("Fourier coeffs ($\mu$s)")
+plt.ylabel("Fourier coeffs ($\\mu$s)")
 plt.xscale("log")
 plt.legend(fontsize=8)
 
@@ -200,6 +204,7 @@
 plt.xlabel("Frequency (Hz)")
 plt.axvline(fyr, color="black", ls="dotted", label="1 yr$^{-1}$")
 plt.legend()
+plt.show()
 
 # %% [markdown]
 # Note the outlier in the 1 year^-1 bin. This is caused by the covariance with RA and DEC, which introduce a delay with the same frequency.
@@ -259,7 +264,12 @@
 
 m2 = deepcopy(m1)
 
-nharm_opt, d_aics = find_optimal_nharms(m2, t, "DMWaveX", 30)
+aics, nharm_opt = find_optimal_nharms(
+    m2, t, include_components=["DMWaveX"], nharms_max=30
+)
+d_aics = aics - np.min(aics)
+nharm_opt = nharm_opt[0]
+
 print("Optimum no of harmonics = ", nharm_opt)
 
 # %%
@@ -273,7 +283,7 @@
 plt.ylabel("AIC - AIC$_\\min{} + 1$")
 plt.legend()
 plt.yscale("log")
-# plt.savefig("sim3-aic.pdf")
+plt.show()
 
 # %%
 # Now create a new model with the optimum number of
@@ -351,7 +361,7 @@
 )
 plt.axvline(fyr, color="black", ls="dotted")
 plt.axhline(0, color="grey", ls="--")
-plt.ylabel("Fourier coeffs ($\mu$s)")
+plt.ylabel("Fourier coeffs ($\\mu$s)")
 plt.xscale("log")
 plt.legend(fontsize=8)
 
@@ -370,3 +380,4 @@
 plt.xlabel("Frequency (Hz)")
 plt.axvline(fyr, color="black", ls="dotted", label="1 yr$^{-1}$")
 plt.legend()
+plt.show()

requirements.txt

@@ -10,3 +10,4 @@ uncertainties
 loguru
 nestle>=0.2.0
 numdifftools
+joblib

setup.cfg

@@ -44,6 +44,7 @@ install_requires =
     loguru
     nestle>=0.2.0
     numdifftools
+    joblib
 
 [options.packages.find]
 where = src

src/pint/logging.py

@@ -141,8 +141,7 @@ def __init__(
         never: Optional[List[str]] = None,
         onlyonce_level: str = "INFO",
     ) -> None:
-        r"""
-        Define regexs for messages that will only be seen once.  Use ``\S+`` for a variable that might change.
+        r"""Define regexes for messages that will only be seen once.  Use ``\S+`` for a variable that might change.
         If a message comes through with a new value for that variable, it will be seen.
 
         Make sure to escape other regex commands like ``()``.

src/pint/models/__init__.py

@@ -26,6 +26,7 @@
 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.cmwavex import CMWaveX
 from pint.models.dispersion_model import (
     DispersionDM,
     DispersionDMX,

src/pint/models/binary_ddk.py

@@ -42,21 +42,26 @@ def _convert_kom(kom):
 
 
 class BinaryDDK(BinaryDD):
-    """Damour and Deruelle model with kinematics.
+    r"""Damour and Deruelle model with kinematics.
 
     This extends the :class:`pint.models.binary_dd.BinaryDD` model with
     "Shklovskii" and "Kopeikin" terms that account for the finite distance
     of the system from Earth, the finite size of the system, and the
     interaction of these with the proper motion.
 
-    From Kopeikin (1995) this includes :math:`\Delta_{\pi M}` (Equation 17), the mixed annual-orbital parallax term, which changes :math:`a_1` and :math:`\omega`
-    (:meth:`~pint.models.stand_alone_psr_binaries.DDK_model.DDKmodel.delta_a1_parallax` and :meth:`~pint.models.stand_alone_psr_binaries.DDK_model.DDKmodel.delta_omega_parallax`).
+    From Kopeikin (1995) this includes :math:`\Delta_{\pi M}` (Equation 17),
+    the mixed annual-orbital parallax term, which changes :math:`a_1` and :math:`\omega`
+    (:meth:`~pint.models.stand_alone_psr_binaries.DDK_model.DDKmodel.delta_a1_parallax`
+    and :meth:`~pint.models.stand_alone_psr_binaries.DDK_model.DDKmodel.delta_omega_parallax`).
 
-    It does not include :math:`\Delta_{\pi P}`, the pure pulsar orbital parallax term (Equation 14).
+    It does not include :math:`\Delta_{\pi P}`, the pure pulsar orbital parallax term
+    (Equation 14).
 
-    From Kopeikin (1996) this includes apparent changes in :math:`\omega`, :math:`a_1`, and :math:`i` due to the proper motion
-    (:meth:`~pint.models.stand_alone_psr_binaries.DDK_model.DDKmodel.delta_omega_proper_motion`, :meth:`~pint.models.stand_alone_psr_binaries.DDK_model.DDKmodel.delta_a1_proper_motion`,
-    :meth:`~pint.models.stand_alone_psr_binaries.DDK_model.DDKmodel.delta_kin_proper_motion`) (Equations 8, 9, 10).
+    From Kopeikin (1996) this includes apparent changes in :math:`\omega`, :math:`a_1`, and
+    :math:`i` due to the proper motion (:meth:`~pint.models.stand_alone_psr_binaries.DDK_model.DDKmodel.delta_omega_proper_motion`,
+    :meth:`~pint.models.stand_alone_psr_binaries.DDK_model.DDKmodel.delta_a1_proper_motion`,
+    :meth:`~pint.models.stand_alone_psr_binaries.DDK_model.DDKmodel.delta_kin_proper_motion`)
+    (Equations 8, 9, 10).
 
     The actual calculations for this are done in
     :class:`pint.models.stand_alone_psr_binaries.DDK_model.DDKmodel`.
@@ -226,14 +231,14 @@ def validate(self):
             warnings.warn("Using A1DOT with a DDK model is not advised.")
 
     def alternative_solutions(self):
-        """Alternative Kopeikin solutions (potential local minima)
+        r"""Alternative Kopeikin solutions (potential local minima)
 
         There are 4 potential local minima for a DDK model where a1dot is the same
         These are given by where Eqn. 8 in Kopeikin (1996) is equal to the best-fit value.
 
         We first define the symmetry point where a1dot is zero (in equatorial coordinates):
 
-        :math:`KOM_0 = \\tan^{-1} (\mu_{\delta} / \mu_{\\alpha})`
+        :math:`KOM_0 = \tan^{-1} (\mu_{\delta} / \mu_{\alpha})`
 
         The solutions are then: