diff --git a/CHANGELOG-unreleased.md b/CHANGELOG-unreleased.md index 25f3d6283..286321f2c 100644 --- a/CHANGELOG-unreleased.md +++ b/CHANGELOG-unreleased.md @@ -19,6 +19,10 @@ the released changes. - 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 @@ -36,5 +40,5 @@ the released changes. - 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/src/pint/models/noise_model.py b/src/pint/models/noise_model.py index b53ebdb2d..7fa364bb9 100644 --- a/src/pint/models/noise_model.py +++ b/src/pint/models/noise_model.py @@ -8,7 +8,7 @@ from loguru import logger as log -from pint.models.parameter import floatParameter, maskParameter +from pint.models.parameter import floatParameter, intParameter, maskParameter from pint.models.timing_model import Component @@ -480,7 +480,7 @@ def __init__( name="TNDMAMP", units="", aliases=[], - description="Amplitude of powerlaw DM noise in tempo2 format", + description="Log-amplitude of powerlaw DM noise in tempo2 format", convert_tcb2tdb=False, ) ) @@ -494,20 +494,40 @@ def __init__( ) ) self.add_param( - floatParameter( + intParameter( name="TNDMC", - units="", + value=30, aliases=[], description="Number of DM noise frequencies.", convert_tcb2tdb=False, ) ) + self.add_param( + floatParameter( + name="TNDMFLOW", + value=0, + description="Fundamental log-frequency of the DM noise Fourier basis as a multiple of the inverse data span (tempo2 format).", + aliases=[], + frozen=True, + convert_tcb2tdb=False, + ) + ) + self.add_param( + floatParameter( + name="TNDMCORNER", + value=0, + description="Corner frequency of the DM noise as a multiple of the inverse data span (tempo2 format).", + aliases=[], + frozen=True, + convert_tcb2tdb=False, + ) + ) self.covariance_matrix_funcs += [self.pl_dm_cov_matrix] self.basis_funcs += [self.pl_dm_basis_weight_pair] def get_pl_vals(self): - nf = int(self.TNDMC.value) if self.TNDMC.value is not None else 30 + nf = self.TNDMC.value amp, gam = 10**self.TNDMAMP.value, self.TNDMGAM.value return (amp, gam, nf) @@ -522,9 +542,20 @@ def get_noise_basis(self, toas): fref = 1400 * u.MHz D = (fref.value / freqs.value) ** 2 nf = self.get_pl_vals()[2] - Fmat = create_fourier_design_matrix(t, nf) + + f_1 = self.get_fundamental_frequency(t) + + Fmat = create_fourier_design_matrix(t, nf, f_1) return Fmat * D[:, None] + def get_fundamental_frequency(self, t): + f_low_factor = 10**self.TNDMFLOW.value + return f_low_factor / (t.max() - t.min()) + + def get_corner_frequency(self, t): + f_c_factor = self.TNDMCORNER.value + return f_c_factor / (t.max() - t.min()) + def get_noise_weights(self, toas): """Return power law DM noise weights. @@ -533,8 +564,14 @@ def get_noise_weights(self, toas): 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] + f_1 = self.get_fundamental_frequency(t) + f_c = self.get_corner_frequency(t) + Ffreqs = get_rednoise_freqs(nf, f_1) + return ( + powerlaw(Ffreqs, amp, gam) + if f_c == 0 + else powerlaw_corner(Ffreqs, amp, gam, f_c) + ) * Ffreqs[0] def pl_dm_basis_weight_pair(self, toas): """Return a Fourier design matrix and power law DM noise weights. @@ -732,7 +769,7 @@ def __init__( name="TNREDAMP", units="", aliases=[], - description="Amplitude of powerlaw red noise in tempo2 format", + description="Log-amplitude of powerlaw red noise in tempo2 format", convert_tcb2tdb=False, ) ) @@ -746,11 +783,31 @@ def __init__( ) ) self.add_param( - floatParameter( + intParameter( name="TNREDC", - units="", - aliases=[], + value=30, description="Number of red noise frequencies.", + aliases=[], + convert_tcb2tdb=False, + ) + ) + self.add_param( + floatParameter( + name="TNREDFLOW", + value=0, + description="Fundamental log-frequency of the red noise Fourier basis as a multiple of the inverse data span (tempo2 format).", + aliases=[], + frozen=True, + convert_tcb2tdb=False, + ) + ) + self.add_param( + floatParameter( + name="TNREDCORNER", + value=0, + description="Corner frequency of the red noise as a multiple of the inverse data span (tempo2 format).", + aliases=[], + frozen=True, convert_tcb2tdb=False, ) ) @@ -759,7 +816,7 @@ def __init__( self.basis_funcs += [self.pl_rn_basis_weight_pair] def get_pl_vals(self): - nf = int(self.TNREDC.value) if self.TNREDC.value is not None else 30 + nf = self.TNREDC.value if self.TNREDAMP.value is not None and self.TNREDGAM.value is not None: amp, gam = 10**self.TNREDAMP.value, self.TNREDGAM.value elif self.RNAMP.value is not None and self.RNIDX is not None: @@ -767,6 +824,14 @@ def get_pl_vals(self): amp, gam = self.RNAMP.value / fac, -1 * self.RNIDX.value return (amp, gam, nf) + def get_fundamental_frequency(self, t): + f_low_factor = 10**self.TNREDFLOW.value + return f_low_factor / (t.max() - t.min()) + + def get_corner_frequency(self, t): + f_c_factor = self.TNREDCORNER.value + return f_c_factor / (t.max() - t.min()) + def get_noise_basis(self, toas): """Return a Fourier design matrix for red noise. @@ -774,8 +839,12 @@ def get_noise_basis(self, toas): tbl = toas.table t = (tbl["tdbld"].quantity * u.day).to(u.s).value + + f_1 = self.get_fundamental_frequency(t) + nf = self.get_pl_vals()[2] - return create_fourier_design_matrix(t, nf) + + return create_fourier_design_matrix(t, nf, f_1) def get_noise_weights(self, toas): """Return power law red noise weights. @@ -785,8 +854,16 @@ def get_noise_weights(self, toas): 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] + + f_1 = self.get_fundamental_frequency(t) + f_c = self.get_corner_frequency(t) + + Ffreqs = get_rednoise_freqs(nf, f_1) + return ( + powerlaw(Ffreqs, amp, gam) + if f_c == 0 + else powerlaw_corner(Ffreqs, amp, gam, f_c) + ) * Ffreqs[0] def pl_rn_basis_weight_pair(self, toas): """Return a Fourier design matrix and power law red noise weights. @@ -844,12 +921,10 @@ def create_ecorr_quantization_matrix(toas_table, dt=1, nmin=2): return U -def get_rednoise_freqs(t, nmodes, Tspan=None): - """Frequency components for creating the red noise basis matrix.""" - - T = Tspan if Tspan is not None else t.max() - t.min() +def get_rednoise_freqs(nmodes, f_1): + """Frequency components for creating the red/DM noise basis matrix.""" - f = np.linspace(1 / T, nmodes / T, nmodes) + f = np.linspace(f_1, nmodes * f_1, nmodes) Ffreqs = np.zeros(2 * nmodes) Ffreqs[::2] = f @@ -858,13 +933,13 @@ def get_rednoise_freqs(t, nmodes, Tspan=None): return Ffreqs -def create_fourier_design_matrix(t, nmodes, Tspan=None): +def create_fourier_design_matrix(t, nmodes, f_1): """ Construct fourier design matrix from eq 11 of Lentati et al, 2013 :param t: vector of time series in seconds :param nmodes: number of fourier coefficients to use - :param Tspan: option to some other Tspan + :param f_1: fundamental frequency of the Fourier basis :return: F: fourier design matrix :return: f: Sampling frequencies """ @@ -872,7 +947,7 @@ def create_fourier_design_matrix(t, nmodes, Tspan=None): N = len(t) F = np.zeros((N, 2 * nmodes)) - Ffreqs = get_rednoise_freqs(t, nmodes, Tspan=Tspan) + Ffreqs = get_rednoise_freqs(nmodes, f_1) F[:, ::2] = np.sin(2 * np.pi * t[:, None] * Ffreqs[::2]) F[:, 1::2] = np.cos(2 * np.pi * t[:, None] * Ffreqs[1::2]) @@ -880,7 +955,7 @@ def create_fourier_design_matrix(t, nmodes, Tspan=None): return F -def powerlaw(f, A=1e-16, gamma=5): +def powerlaw(f, A, gamma): """Power-law PSD. :param f: Sampling frequencies @@ -890,3 +965,23 @@ def powerlaw(f, A=1e-16, gamma=5): fyr = 1 / 3.16e7 return A**2 / 12.0 / np.pi**2 * fyr ** (gamma - 3) * f ** (-gamma) + + +def powerlaw_corner(f, A, gamma, fc): + """Power-law PSD that saturates at low frequencies. + + :param f: Sampling frequencies + :param A: Amplitude of red noise [GW units] + :param gamma: Spectral index of red noise process + :param fc: Corner frequency + """ + + fyr = 1 / 3.16e7 + + return ( + A**2 + / 12.0 + / np.pi**2 + / fyr**3 + * ((1 + (fyr / fc) ** (gamma / 2)) / (1 + (f / fc) ** (gamma / 2))) ** 2 + ) diff --git a/src/pint/utils.py b/src/pint/utils.py index 758307558..00e9c4540 100644 --- a/src/pint/utils.py +++ b/src/pint/utils.py @@ -3253,12 +3253,21 @@ def plrednoise_from_wavex( tnredc = len(model.components["WaveX"].get_indices()) + tnredflow = ( + 0 + if (model.START.quantity is None or model.FINISH.quantity is None) + else np.log10( + (model.FINISH.quantity - model.START.quantity) * model.WXFREQ_0001.quantity + ).si.value + ) + model1 = deepcopy(model) model1.remove_component("WaveX") model1.add_component(PLRedNoise()) model1.TNREDAMP.value = log10_A_val model1.TNREDGAM.value = gamma_val model1.TNREDC.value = tnredc + model1.TNREDFLOW.value = tnredflow model1.TNREDAMP.uncertainty_value = log10_A_err model1.TNREDGAM.uncertainty_value = gamma_err @@ -3307,12 +3316,22 @@ def pldmnoise_from_dmwavex( tndmc = len(model.components["DMWaveX"].get_indices()) + tndmflow = ( + 0 + if (model.START.quantity is None or model.FINISH.quantity is None) + else np.log10( + (model.FINISH.quantity - model.START.quantity) + * model.DMWXFREQ_0001.quantity + ).si.value + ) + model1 = deepcopy(model) model1.remove_component("DMWaveX") model1.add_component(PLDMNoise()) model1.TNDMAMP.value = log10_A_val model1.TNDMGAM.value = gamma_val model1.TNDMC.value = tndmc + model1.TNDMFLOW.value = tndmflow model1.TNDMAMP.uncertainty_value = log10_A_err model1.TNDMGAM.uncertainty_value = gamma_err diff --git a/tests/test_gls_fitter.py b/tests/test_gls_fitter.py index 9ad629d36..76ee7158e 100644 --- a/tests/test_gls_fitter.py +++ b/tests/test_gls_fitter.py @@ -1,6 +1,5 @@ import json import os -import pytest import astropy.units as u import numpy as np