diff --git a/CHANGELOG-unreleased.md b/CHANGELOG-unreleased.md index 8817082c1..2663d3d68 100644 --- a/CHANGELOG-unreleased.md +++ b/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. diff --git a/docs/_ext/componentlist.py b/docs/_ext/componentlist.py index 63bffd6c3..95e1115c4 100644 --- a/docs/_ext/componentlist.py +++ b/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}" diff --git a/src/pint/models/__init__.py b/src/pint/models/__init__.py index ab7820027..4567ba882 100644 --- a/src/pint/models/__init__.py +++ b/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", diff --git a/src/pint/models/noise_model.py b/src/pint/models/noise_model.py index c2f0cfc61..953450074 100644 --- a/src/pint/models/noise_model.py +++ b/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 """ diff --git a/src/pint/utils.py b/src/pint/utils.py index f2bccac4e..c6e67f5b3 100644 --- a/src/pint/utils.py +++ b/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 diff --git a/tests/test_noise_models.py b/tests/test_noise_models.py index 7165c772f..fea45e9c7 100644 --- a/tests/test_noise_models.py +++ b/tests/test_noise_models.py @@ -25,14 +25,24 @@ 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() @@ -40,7 +50,8 @@ 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 diff --git a/tests/test_plchromnoise.py b/tests/test_plchromnoise.py new file mode 100644 index 000000000..46c8db376 --- /dev/null +++ b/tests/test_plchromnoise.py @@ -0,0 +1,166 @@ +"""Tests to ensure that ModelBuilder is able to read PLChromNoise +component properly. + +Tests that PLChromNoise reproduces same results as PLRedNoise using +monochromatic data +""" + +import astropy.units as u +import contextlib +import io +import numpy as np +import pint.fitter as fitters +import pint.models.model_builder as mb +from pint.models.model_builder import get_model +from pint.models.timing_model import Component +from pint.residuals import Residuals +from pint.simulation import make_fake_toas_fromMJDs +import pytest + +parfile_contents = """ + PSRJ J0023+0923 + RAJ 00:23:16.8790858 1 0.00002408141295805134 + DECJ +09:23:23.86936 1 0.00082010713730773120 + F0 327.84702062954047136 1 0.00000000000295205483 + F1 -1.2278326306812866375e-15 1 3.8219244605614075223e-19 + PEPOCH 56199.999797564144902 + POSEPOCH 56199.999797564144902 + DMEPOCH 56200 + DM 14.327978186774068347 1 0.00006751663559857748 + BINARY ELL1 + PB 0.13879914244858396754 1 0.00000000003514075083 + A1 0.034841158415224894973 1 0.00000012173038389247 + TASC 56178.804891768506529 1 0.00000007765191894742 + EPS1 1.6508830631753595232e-05 1 0.00000477568412215803 + EPS2 3.9656838708709247373e-06 1 0.00000458951091435993 + CLK TT(BIPM2015) + MODE 1 + UNITS TDB + TIMEEPH FB90 + DILATEFREQ N + PLANET_SHAPIRO N + CORRECT_TROPOSPHERE N + EPHEM DE436 + JUMP -fe 430 0 0 + TNEF -group 430_ASP 1.0389 + TNEQ -group 430_ASP -8.77109 + TNECORR -group 430_PUPPI 0.00490558 + TNRedAmp -13.3087 + TNRedGam 0.125393 + TNRedC 14 + TNChromAMP -14.2 + TNChromGAM 0.624 + TNChromC 70 + TNChromIdx 4 + CM0 0 1 + CM1 0 1 +""" + + +@pytest.fixture() +def modelJ0023p0923(): + return mb.get_model(io.StringIO(parfile_contents)) + + +def test_read_PLChromNoise_component(modelJ0023p0923): + assert "PLChromNoise" in modelJ0023p0923.components + + +def test_read_PLChromNoise_component_type(modelJ0023p0923): + assert ( + modelJ0023p0923.components["PLChromNoise"] + in modelJ0023p0923.NoiseComponent_list + ) + + +def test_read_PLChromNoise_params(modelJ0023p0923): + params = ["TNCHROMAMP", "TNCHROMGAM", "TNCHROMC", "TNCHROMIDX"] + for param in params: + assert ( + hasattr(modelJ0023p0923, param) + and getattr(modelJ0023p0923, param).quantity is not None + ) + + +def test_PLRedNoise_recovery(): + # basic model, no EFAC or EQUAD + model = get_model( + io.StringIO( + """ + PSRJ J1234+5678 + ELAT 0 1 + ELONG 0 1 + DM 10 0 + F0 1 1 + PEPOCH 58000 + POSEPOCH 58000 + PMELONG 0 1 + PMELAT 0 1 + PX 0 1 + TNEF mjd 57000 58000 1.0389 + TNEQ mjd 57000 58000 -8.77109 + TNRedAmp -11 + TNRedGam 3 + TNRedC 30 + """ + ) + ) + + # simulate toas + MJDs = np.linspace(57001, 58000, 150, dtype=np.longdouble) * u.d + toas = make_fake_toas_fromMJDs(MJDs, model=model, error=1 * u.us, add_noise=True) + + # get red noise basis and weights + rn_basis = model.components["PLRedNoise"].get_noise_basis(toas) + rn_weights = model.components["PLRedNoise"].get_noise_weights(toas) + + # fit model with red noise + f1 = fitters.DownhillGLSFitter(toas, model) + f1.fit_toas() + f1.model.validate() + f1.model.validate_toas(toas) + r1 = Residuals(toas, f1.model) + + # remove red noise + A = model["TNREDAMP"].value + gam = model["TNREDGAM"].value + c = model["TNREDC"].value + model.remove_component("PLRedNoise") + + # create and add PLChromNoise component + all_components = Component.component_types + cm_class = all_components["ChromaticCM"] + model.add_component(cm_class(), validate=False) + model["TNCHROMIDX"].value = 4 + noise_class = all_components["PLChromNoise"] + model.add_component(noise_class(), validate=False) + model["TNCHROMAMP"].quantity = A + model["TNCHROMGAM"].quantity = gam + model["TNCHROMC"].value = c + model.validate() + + # get chromatic noise basis and weights + chrom_basis = model.components["PLChromNoise"].get_noise_basis(toas) + # chromatic basis will be off by a constant factor, depending on frequency + # of the test data + D = (1400 / toas.get_freqs().value) ** 2 + chrom_basis_scaled = chrom_basis / D[:, None] + chrom_weights = model.components["PLChromNoise"].get_noise_weights(toas) + + # refit model + f2 = fitters.DownhillGLSFitter(toas, model) + with contextlib.suppress(fitters.InvalidModelParameters, fitters.StepProblem): + f2.fit_toas() + f2.model.validate() + f2.model.validate_toas(toas) + r2 = Residuals(toas, f2.model) + + # check weights and basis are equivalent within error + basis_diff = rn_basis - chrom_basis_scaled + weights_diff = rn_weights - chrom_weights + assert np.all(np.isclose(basis_diff, 0, atol=1e-3)) + assert np.all(np.isclose(weights_diff, 0)) + + # check residuals are equivalent within error + rs_diff = r2.time_resids.value - r1.time_resids.value + assert np.all(np.isclose(rs_diff, 0, atol=1e-6))