Merge branch 'master' into redcorner

CHANGELOG-unreleased.md

@@ -18,7 +18,10 @@ the released changes.
 - `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
 ### Fixed
+- Explicit type conversion in `woodbury_dot()` function
+- Documentation: Fixed empty descriptions in the timing model components table
 ### 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.

docs/_ext/componentlist.py

@@ -19,7 +19,7 @@ def run(self):
             class_ = d[k]
             full_name = f"{class_.__module__}.{class_.__name__}"
             if hasattr(class_, "__doc__") and class_.__doc__ is not None:
-                doc = class_.__doc__.split("\n")[0].strip()
+                doc = class_.__doc__.strip().split("\n")[0].strip()
             else:
                 doc = ""
             msg = f"* :class:`~{full_name}` - {doc}"

src/pint/models/__init__.py

@@ -40,7 +40,13 @@
 from pint.models.ifunc import IFunc
 from pint.models.jump import DelayJump, PhaseJump
 from pint.models.model_builder import get_model, get_model_and_toas
-from pint.models.noise_model import EcorrNoise, PLRedNoise, ScaleToaError
+from pint.models.noise_model import (
+    EcorrNoise,
+    PLRedNoise,
+    PLDMNoise,
+    PLChromNoise,
+    ScaleToaError,
+)
 from pint.models.solar_system_shapiro import SolarSystemShapiro
 from pint.models.solar_wind_dispersion import SolarWindDispersion, SolarWindDispersionX
 from pint.models.spindown import Spindown
@@ -57,6 +63,7 @@
     "StandardTimingModel",
     [AstrometryEquatorial(), Spindown(), DispersionDM(), SolarSystemShapiro()],
 )
+
 # BTTimingModel = generate_timing_model("BTTimingModel",
 #         (Astrometry, Spindown, Dispersion, SolarSystemShapiro, BT))
 # DDTimingModel = generate_timing_model("DDTimingModel",

src/pint/models/noise_model.py

@@ -460,7 +460,7 @@ class PLDMNoise(NoiseComponent):
 
     Note
     ----
-    Ref: Lentati et al. 2014
+    Ref: Lentati et al. 2014, MNRAS 437(3), 3004-3023
 
     """
 
@@ -593,6 +593,127 @@ def pl_dm_cov_matrix(self, toas):
         return np.dot(Fmat * phi[None, :], Fmat.T)
 
 
+class PLChromNoise(NoiseComponent):
+    """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
+    a transfer function that is determined by the geometry and electron distribution
+    in the scattering screen(s). The scattering timescale is typically a decreasing
+    function of the observing frequency.
+
+    Scatter broadening causes systematic offsets in the TOA measurements due to the
+    pulse shape mismatch. While this offset need not be a simple function of frequency,
+    it has been often modeled using a delay that is proportional to f^-alpha where alpha
+    is known as the chromatic index.
+
+    This model should be used in combination with the ChromaticCM model.
+
+    Parameters supported:
+
+    .. paramtable::
+        :class: pint.models.noise_model.PLChromNoise
+
+    Note
+    ----
+    Ref: Lentati et al. 2014, MNRAS 437(3), 3004-3023
+
+    """
+
+    register = True
+    category = "pl_chrom_noise"
+
+    introduces_correlated_errors = True
+    is_time_correlated = True
+
+    def __init__(
+        self,
+    ):
+        super().__init__()
+
+        self.add_param(
+            floatParameter(
+                name="TNCHROMAMP",
+                units="",
+                aliases=[],
+                description="Amplitude of powerlaw chromatic noise in tempo2 format",
+                convert_tcb2tdb=False,
+            )
+        )
+        self.add_param(
+            floatParameter(
+                name="TNCHROMGAM",
+                units="",
+                aliases=[],
+                description="Spectral index of powerlaw chromatic noise in tempo2 format",
+                convert_tcb2tdb=False,
+            )
+        )
+        self.add_param(
+            floatParameter(
+                name="TNCHROMC",
+                units="",
+                aliases=[],
+                description="Number of chromatic noise frequencies.",
+                convert_tcb2tdb=False,
+            )
+        )
+
+        self.covariance_matrix_funcs += [self.pl_chrom_cov_matrix]
+        self.basis_funcs += [self.pl_chrom_basis_weight_pair]
+
+    def get_pl_vals(self):
+        nf = int(self.TNCHROMC.value) if self.TNCHROMC.value is not None else 30
+        amp, gam = 10**self.TNCHROMAMP.value, self.TNCHROMGAM.value
+        return (amp, gam, nf)
+
+    def get_noise_basis(self, toas):
+        """Return a Fourier design matrix for chromatic noise.
+
+        See the documentation for pl_chrom_basis_weight_pair function for details."""
+
+        tbl = toas.table
+        t = (tbl["tdbld"].quantity * u.day).to(u.s).value
+        freqs = self._parent.barycentric_radio_freq(toas).to(u.MHz)
+        fref = 1400 * u.MHz
+        alpha = self._parent.TNCHROMIDX.value
+        D = (fref.value / freqs.value) ** alpha
+        nf = self.get_pl_vals()[2]
+        Fmat = create_fourier_design_matrix(t, nf)
+        return Fmat * D[:, None]
+
+    def get_noise_weights(self, toas):
+        """Return power law chromatic noise weights.
+
+        See the documentation for pl_chrom_basis_weight_pair for details."""
+
+        tbl = toas.table
+        t = (tbl["tdbld"].quantity * u.day).to(u.s).value
+        amp, gam, nf = self.get_pl_vals()
+        Ffreqs = get_rednoise_freqs(t, nf)
+        return powerlaw(Ffreqs, amp, gam) * Ffreqs[0]
+
+    def pl_chrom_basis_weight_pair(self, toas):
+        """Return a Fourier design matrix and power law chromatic noise weights.
+
+        A Fourier design matrix contains the sine and cosine basis_functions
+        in a Fourier series expansion. Here we scale the design matrix by
+        (fref/f)**2, where fref = 1400 MHz to match the convention used in
+        enterprise.
+
+        The weights used are the power-law PSD values at frequencies n/T,
+        where n is in [1, TNCHROMC] and T is the total observing duration of
+        the dataset.
+
+        """
+        return (self.get_noise_basis(toas), self.get_noise_weights(toas))
+
+    def pl_chrom_cov_matrix(self, toas):
+        Fmat, phi = self.pl_chrom_basis_weight_pair(toas)
+        return np.dot(Fmat * phi[None, :], Fmat.T)
+
+
 class PLRedNoise(NoiseComponent):
     """Timing noise with a power-law spectrum.
 
@@ -610,7 +731,7 @@ class PLRedNoise(NoiseComponent):
 
     Note
     ----
-    Ref: NANOGrav 11 yrs data
+    Ref: Lentati et al. 2014, MNRAS 437(3), 3004-3023
 
     """
 

src/pint/utils.py

@@ -3026,7 +3026,7 @@ def woodbury_dot(
 
     logdet_N = np.sum(np.log(Ndiag))
     logdet_Phi = np.sum(np.log(Phidiag))
-    _, logdet_Sigma = np.linalg.slogdet(Sigma)
+    _, logdet_Sigma = np.linalg.slogdet(Sigma.astype(float))
 
     logdet_C = logdet_N + logdet_Phi + logdet_Sigma
 

tests/test_noise_models.py

@@ -25,22 +25,33 @@
 
 def add_DM_noise_to_model(model):
     all_components = Component.component_types
-    noise_class = all_components["PLDMNoise"]
-    noise = noise_class()  # Make the DM noise instance.
-    model.add_component(noise, validate=False)
+    model.add_component(all_components["PLDMNoise"](), validate=False)
     model["TNDMAMP"].quantity = 1e-13
     model["TNDMGAM"].quantity = 1.2
     model["TNDMC"].value = 30
     model.validate()
-    return model
+
+
+def add_chrom_noise_to_model(model):
+    all_components = Component.component_types
+    model.add_component(all_components["PLChromNoise"](), validate=False)
+    model["TNCHROMAMP"].quantity = 1e-14
+    model["TNCHROMGAM"].quantity = 1.2
+    model["TNCHROMC"].value = 30
+
+    model.add_component(all_components["ChromaticCM"](), validate=False)
+    model["TNCHROMIDX"].value = 4
+
+    model.validate()
 
 
 @pytest.fixture()
 def model_and_toas():
     parfile = examplefile("B1855+09_NANOGrav_9yv1.gls.par")
     timfile = examplefile("B1855+09_NANOGrav_9yv1.tim")
     model, toas = get_model_and_toas(parfile, timfile)
-    model = add_DM_noise_to_model(model)
+    add_DM_noise_to_model(model)
+    add_chrom_noise_to_model(model)
     return model, toas