From 9f822fc559c9ce0b037d35720cc12cf3cebcfdab Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Fri, 10 May 2024 11:24:37 +0200 Subject: [PATCH 01/19] TNRedC as an intParameter --- src/pint/models/noise_model.py | 20 ++++++++++++++++---- 1 file changed, 16 insertions(+), 4 deletions(-) diff --git a/src/pint/models/noise_model.py b/src/pint/models/noise_model.py index debe7efdc..1434a68bc 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 @@ -613,11 +613,10 @@ def __init__( ) ) self.add_param( - floatParameter( + intParameter( name="TNREDC", - units="", - aliases=[], description="Number of red noise frequencies.", + aliases=[], ) ) @@ -756,3 +755,16 @@ 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=1e-16, gamma=5, fc=1e9): + """Power-law PSD. + + :param f: Sampling frequencies + :param A: Amplitude of red noise [GW units] + :param gamma: Spectral index of red noise process + :param fc: The corner frequency for spectral turnover + """ + + fyr = 1 / 3.16e7 + return A**2 / 12.0 / np.pi**2 * fyr ** (gamma - 3) * f ** (-gamma) From c95499f00b714e4b921bce00b70d2a00773d8dc0 Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Fri, 10 May 2024 11:27:48 +0200 Subject: [PATCH 02/19] changelog --- CHANGELOG-unreleased.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGELOG-unreleased.md b/CHANGELOG-unreleased.md index 73b881a60..293c47fce 100644 --- a/CHANGELOG-unreleased.md +++ b/CHANGELOG-unreleased.md @@ -15,6 +15,7 @@ the released changes. - Better documentation for `akaike_information_criterion()` - `funcParameter`s are no longer listed in the `pintk` interface. - Updated location of CCERA +- `TNRedC` to an `intParameter` ### Added - `bayesian_information_criterion()` function - `pintk` now reads and automatically converts TCB par files and par files with `BINARY T2`. From a0b02605f4cdad93dd3b3e733773f7bf0baa1e3d Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Fri, 10 May 2024 11:29:46 +0200 Subject: [PATCH 03/19] default tnredc and tndmc --- src/pint/models/noise_model.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/src/pint/models/noise_model.py b/src/pint/models/noise_model.py index 1434a68bc..b0dc16e55 100644 --- a/src/pint/models/noise_model.py +++ b/src/pint/models/noise_model.py @@ -486,9 +486,9 @@ def __init__( ) ) self.add_param( - floatParameter( + intParameter( name="TNDMC", - units="", + value=30, aliases=[], description="Number of DM noise frequencies.", ) @@ -498,7 +498,7 @@ def __init__( 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) @@ -615,6 +615,7 @@ def __init__( self.add_param( intParameter( name="TNREDC", + value=30, description="Number of red noise frequencies.", aliases=[], ) @@ -624,7 +625,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: From a1f8cb44265949246947d00f7b16e4ba9c7af9e8 Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Fri, 10 May 2024 11:30:16 +0200 Subject: [PATCH 04/19] CHANGELOG --- CHANGELOG-unreleased.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CHANGELOG-unreleased.md b/CHANGELOG-unreleased.md index 293c47fce..c4078b8c7 100644 --- a/CHANGELOG-unreleased.md +++ b/CHANGELOG-unreleased.md @@ -15,7 +15,7 @@ the released changes. - Better documentation for `akaike_information_criterion()` - `funcParameter`s are no longer listed in the `pintk` interface. - Updated location of CCERA -- `TNRedC` to an `intParameter` +- `TNRedC` and `TNDMC` into `intParameter`s ### Added - `bayesian_information_criterion()` function - `pintk` now reads and automatically converts TCB par files and par files with `BINARY T2`. From 9fee879d330df047263b9438761be4de1b37d17d Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Fri, 10 May 2024 11:34:17 +0200 Subject: [PATCH 05/19] TNREDFLOW & TNDMFLOW --- src/pint/models/noise_model.py | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/src/pint/models/noise_model.py b/src/pint/models/noise_model.py index b0dc16e55..8ef811a34 100644 --- a/src/pint/models/noise_model.py +++ b/src/pint/models/noise_model.py @@ -493,6 +493,13 @@ def __init__( description="Number of DM noise frequencies.", ) ) + self.add_param( + floatParameter( + name="TNDMFLOW", + description="Fundamental log-frequency of the DM noise Fourier basis.", + aliases=[], + ) + ) self.covariance_matrix_funcs += [self.pl_dm_cov_matrix] self.basis_funcs += [self.pl_dm_basis_weight_pair] @@ -620,6 +627,13 @@ def __init__( aliases=[], ) ) + self.add_param( + floatParameter( + name="TNREDFLOW", + description="Fundamental log-frequency of the red noise Fourier basis.", + aliases=[], + ) + ) self.covariance_matrix_funcs += [self.pl_rn_cov_matrix] self.basis_funcs += [self.pl_rn_basis_weight_pair] From d58cffc85cc09a78a11b527812d0e5c4ce75c3bd Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Fri, 10 May 2024 12:28:00 +0200 Subject: [PATCH 06/19] get_Tspan --- CHANGELOG-unreleased.md | 2 +- src/pint/models/noise_model.py | 63 +++++++++++++++++++--------------- 2 files changed, 36 insertions(+), 29 deletions(-) diff --git a/CHANGELOG-unreleased.md b/CHANGELOG-unreleased.md index c4078b8c7..9134d7130 100644 --- a/CHANGELOG-unreleased.md +++ b/CHANGELOG-unreleased.md @@ -15,7 +15,7 @@ the released changes. - Better documentation for `akaike_information_criterion()` - `funcParameter`s are no longer listed in the `pintk` interface. - Updated location of CCERA -- `TNRedC` and `TNDMC` into `intParameter`s +- Changed `TNRedC` and `TNDMC` into `intParameter`s; made their default values explicit. ### Added - `bayesian_information_criterion()` function - `pintk` now reads and automatically converts TCB par files and par files with `BINARY T2`. diff --git a/src/pint/models/noise_model.py b/src/pint/models/noise_model.py index 8ef811a34..48d552172 100644 --- a/src/pint/models/noise_model.py +++ b/src/pint/models/noise_model.py @@ -496,7 +496,8 @@ def __init__( self.add_param( floatParameter( name="TNDMFLOW", - description="Fundamental log-frequency of the DM noise Fourier basis.", + value=0, + description="Fundamental log-frequency of the DM noise Fourier basis as a multiple of the inverse data span (tempo2 format).", aliases=[], ) ) @@ -520,9 +521,17 @@ 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_low_factor = 10**self.TNDMFLOW.value + Tspan = (t.max() - t.min()) / f_low_factor + + Fmat = create_fourier_design_matrix(t, nf, Tspan) return Fmat * D[:, None] + def get_Tspan(self, t): + f_low_factor = 10**self.TNDMFLOW.value + return (t.max() - t.min()) / f_low_factor + def get_noise_weights(self, toas): """Return power law DM noise weights. @@ -531,7 +540,8 @@ 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) + Tspan = self.get_Tspan(t) + Ffreqs = get_rednoise_freqs(nf, Tspan) return powerlaw(Ffreqs, amp, gam) * Ffreqs[0] def pl_dm_basis_weight_pair(self, toas): @@ -608,7 +618,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", ) ) self.add_param( @@ -630,7 +640,8 @@ def __init__( self.add_param( floatParameter( name="TNREDFLOW", - description="Fundamental log-frequency of the red noise Fourier basis.", + value=0, + description="Fundamental log-frequency of the red noise Fourier basis as a multiple of the inverse data span (tempo2 format).", aliases=[], ) ) @@ -647,6 +658,10 @@ def get_pl_vals(self): amp, gam = self.RNAMP.value / fac, -1 * self.RNIDX.value return (amp, gam, nf) + def get_Tspan(self, t): + f_low_factor = 10**self.TNREDFLOW.value + return (t.max() - t.min()) / f_low_factor + def get_noise_basis(self, toas): """Return a Fourier design matrix for red noise. @@ -654,8 +669,12 @@ def get_noise_basis(self, toas): tbl = toas.table t = (tbl["tdbld"].quantity * u.day).to(u.s).value + + Tspan = self.get_Tspan(t) + nf = self.get_pl_vals()[2] - return create_fourier_design_matrix(t, nf) + + return create_fourier_design_matrix(t, nf, Tspan) def get_noise_weights(self, toas): """Return power law red noise weights. @@ -665,7 +684,10 @@ 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) + + Tspan = self.get_Tspan(t) + + Ffreqs = get_rednoise_freqs(nf, Tspan) return powerlaw(Ffreqs, amp, gam) * Ffreqs[0] def pl_rn_basis_weight_pair(self, toas): @@ -724,12 +746,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.""" +def get_rednoise_freqs(nmodes, Tspan): + """Frequency components for creating the red/DM noise basis matrix.""" - T = Tspan if Tspan is not None else t.max() - t.min() - - f = np.linspace(1 / T, nmodes / T, nmodes) + f = np.linspace(1 / Tspan, nmodes / Tspan, nmodes) Ffreqs = np.zeros(2 * nmodes) Ffreqs[::2] = f @@ -738,13 +758,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, Tspan): """ 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 Tspan: fundamental period of the Fourier basis :return: F: fourier design matrix :return: f: Sampling frequencies """ @@ -752,7 +772,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, Tspan) F[:, ::2] = np.sin(2 * np.pi * t[:, None] * Ffreqs[::2]) F[:, 1::2] = np.cos(2 * np.pi * t[:, None] * Ffreqs[1::2]) @@ -770,16 +790,3 @@ 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=1e-16, gamma=5, fc=1e9): - """Power-law PSD. - - :param f: Sampling frequencies - :param A: Amplitude of red noise [GW units] - :param gamma: Spectral index of red noise process - :param fc: The corner frequency for spectral turnover - """ - - fyr = 1 / 3.16e7 - return A**2 / 12.0 / np.pi**2 * fyr ** (gamma - 3) * f ** (-gamma) From c1e653ffbb14e12920a72f3c4acb1e56cec84f1d Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Fri, 10 May 2024 12:41:54 +0200 Subject: [PATCH 07/19] plrednoise_from_wavex and pldmnoise_from_dmwavex --- src/pint/models/noise_model.py | 2 ++ src/pint/utils.py | 19 +++++++++++++++++++ 2 files changed, 21 insertions(+) diff --git a/src/pint/models/noise_model.py b/src/pint/models/noise_model.py index 48d552172..74e92a081 100644 --- a/src/pint/models/noise_model.py +++ b/src/pint/models/noise_model.py @@ -499,6 +499,7 @@ def __init__( 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, ) ) @@ -643,6 +644,7 @@ def __init__( 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, ) ) diff --git a/src/pint/utils.py b/src/pint/utils.py index 629a33a70..ff41773a3 100644 --- a/src/pint/utils.py +++ b/src/pint/utils.py @@ -2927,12 +2927,21 @@ def plrednoise_from_wavex(model, ignore_fyr=True): 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 @@ -2979,12 +2988,22 @@ def pldmnoise_from_dmwavex(model, ignore_fyr=False): 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 From 6ecc3db6b215971c4e4a2e72d138994bfc64761e Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Fri, 10 May 2024 12:44:13 +0200 Subject: [PATCH 08/19] CHANGELOG --- CHANGELOG-unreleased.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGELOG-unreleased.md b/CHANGELOG-unreleased.md index 9134d7130..92217def0 100644 --- a/CHANGELOG-unreleased.md +++ b/CHANGELOG-unreleased.md @@ -21,6 +21,7 @@ the released changes. - `pintk` now reads and automatically converts TCB par files and par files with `BINARY T2`. - Test for `pint.utils.split_swx()` - Custom type definitions for type hints +- `plrednoise_from_wavex()` and `pldmnoise_from_dmwavex()` functions now compute `TNRedFLow` and `TNDMFLow` ### Fixed - `pint.utils.split_swx()` to use updated `SolarWindDispersionX()` parameter naming convention - Fix #1759 by changing order of comparison From bdb12a8b722822fe7252359ecdb175416a22d7bf Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Fri, 10 May 2024 12:57:42 +0200 Subject: [PATCH 09/19] -- --- src/pint/models/noise_model.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/pint/models/noise_model.py b/src/pint/models/noise_model.py index 74e92a081..b4cf957bf 100644 --- a/src/pint/models/noise_model.py +++ b/src/pint/models/noise_model.py @@ -523,8 +523,7 @@ def get_noise_basis(self, toas): D = (fref.value / freqs.value) ** 2 nf = self.get_pl_vals()[2] - f_low_factor = 10**self.TNDMFLOW.value - Tspan = (t.max() - t.min()) / f_low_factor + Tspan = self.get_Tspan(t) Fmat = create_fourier_design_matrix(t, nf, Tspan) return Fmat * D[:, None] From 210a9c1c336c0da927ba80f711ccb58231f31b88 Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Fri, 10 May 2024 14:03:43 +0200 Subject: [PATCH 10/19] get_F1 --- src/pint/models/noise_model.py | 34 +++++++++++++++++----------------- 1 file changed, 17 insertions(+), 17 deletions(-) diff --git a/src/pint/models/noise_model.py b/src/pint/models/noise_model.py index b4cf957bf..e43944bf9 100644 --- a/src/pint/models/noise_model.py +++ b/src/pint/models/noise_model.py @@ -523,14 +523,14 @@ def get_noise_basis(self, toas): D = (fref.value / freqs.value) ** 2 nf = self.get_pl_vals()[2] - Tspan = self.get_Tspan(t) + F1 = self.get_F1(t) - Fmat = create_fourier_design_matrix(t, nf, Tspan) + Fmat = create_fourier_design_matrix(t, nf, F1) return Fmat * D[:, None] - def get_Tspan(self, t): + def get_F1(self, t): f_low_factor = 10**self.TNDMFLOW.value - return (t.max() - t.min()) / f_low_factor + return f_low_factor / (t.max() - t.min()) def get_noise_weights(self, toas): """Return power law DM noise weights. @@ -540,8 +540,8 @@ 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() - Tspan = self.get_Tspan(t) - Ffreqs = get_rednoise_freqs(nf, Tspan) + F1 = self.get_F1(t) + Ffreqs = get_rednoise_freqs(nf, F1) return powerlaw(Ffreqs, amp, gam) * Ffreqs[0] def pl_dm_basis_weight_pair(self, toas): @@ -659,9 +659,9 @@ def get_pl_vals(self): amp, gam = self.RNAMP.value / fac, -1 * self.RNIDX.value return (amp, gam, nf) - def get_Tspan(self, t): + def get_F1(self, t): f_low_factor = 10**self.TNREDFLOW.value - return (t.max() - t.min()) / f_low_factor + return f_low_factor / (t.max() - t.min()) def get_noise_basis(self, toas): """Return a Fourier design matrix for red noise. @@ -671,11 +671,11 @@ def get_noise_basis(self, toas): tbl = toas.table t = (tbl["tdbld"].quantity * u.day).to(u.s).value - Tspan = self.get_Tspan(t) + F1 = self.get_F1(t) nf = self.get_pl_vals()[2] - return create_fourier_design_matrix(t, nf, Tspan) + return create_fourier_design_matrix(t, nf, F1) def get_noise_weights(self, toas): """Return power law red noise weights. @@ -686,9 +686,9 @@ def get_noise_weights(self, toas): t = (tbl["tdbld"].quantity * u.day).to(u.s).value amp, gam, nf = self.get_pl_vals() - Tspan = self.get_Tspan(t) + F1 = self.get_F1(t) - Ffreqs = get_rednoise_freqs(nf, Tspan) + Ffreqs = get_rednoise_freqs(nf, F1) return powerlaw(Ffreqs, amp, gam) * Ffreqs[0] def pl_rn_basis_weight_pair(self, toas): @@ -747,10 +747,10 @@ def create_ecorr_quantization_matrix(toas_table, dt=1, nmin=2): return U -def get_rednoise_freqs(nmodes, Tspan): +def get_rednoise_freqs(nmodes, F1): """Frequency components for creating the red/DM noise basis matrix.""" - f = np.linspace(1 / Tspan, nmodes / Tspan, nmodes) + f = np.linspace(F1, nmodes * F1, nmodes) Ffreqs = np.zeros(2 * nmodes) Ffreqs[::2] = f @@ -759,13 +759,13 @@ def get_rednoise_freqs(nmodes, Tspan): return Ffreqs -def create_fourier_design_matrix(t, nmodes, Tspan): +def create_fourier_design_matrix(t, nmodes, F1): """ 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: fundamental period of the Fourier basis + :param F1: fundamental frequency of the Fourier basis :return: F: fourier design matrix :return: f: Sampling frequencies """ @@ -773,7 +773,7 @@ def create_fourier_design_matrix(t, nmodes, Tspan): N = len(t) F = np.zeros((N, 2 * nmodes)) - Ffreqs = get_rednoise_freqs(nmodes, Tspan) + Ffreqs = get_rednoise_freqs(nmodes, F1) F[:, ::2] = np.sin(2 * np.pi * t[:, None] * Ffreqs[::2]) F[:, 1::2] = np.cos(2 * np.pi * t[:, None] * Ffreqs[1::2]) From 89722482f6a0cb29b675a4c7ee553ff681d6242d Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Fri, 10 May 2024 14:06:11 +0200 Subject: [PATCH 11/19] f_1 --- CHANGELOG-unreleased.md | 1 + src/pint/models/noise_model.py | 10 +++++----- 2 files changed, 6 insertions(+), 5 deletions(-) diff --git a/CHANGELOG-unreleased.md b/CHANGELOG-unreleased.md index 92217def0..2d2000532 100644 --- a/CHANGELOG-unreleased.md +++ b/CHANGELOG-unreleased.md @@ -16,6 +16,7 @@ the released changes. - `funcParameter`s are no longer listed in the `pintk` interface. - Updated location of CCERA - Changed `TNRedC` and `TNDMC` into `intParameter`s; made their default values explicit. +- `get_rednoise_freqs()` and `get_rednoise_freqs()` now take `f_1` as input instead of `Tspan`. ### Added - `bayesian_information_criterion()` function - `pintk` now reads and automatically converts TCB par files and par files with `BINARY T2`. diff --git a/src/pint/models/noise_model.py b/src/pint/models/noise_model.py index e43944bf9..a226ef5d4 100644 --- a/src/pint/models/noise_model.py +++ b/src/pint/models/noise_model.py @@ -747,10 +747,10 @@ def create_ecorr_quantization_matrix(toas_table, dt=1, nmin=2): return U -def get_rednoise_freqs(nmodes, F1): +def get_rednoise_freqs(nmodes, f_1): """Frequency components for creating the red/DM noise basis matrix.""" - f = np.linspace(F1, nmodes * F1, nmodes) + f = np.linspace(f_1, nmodes * f_1, nmodes) Ffreqs = np.zeros(2 * nmodes) Ffreqs[::2] = f @@ -759,13 +759,13 @@ def get_rednoise_freqs(nmodes, F1): return Ffreqs -def create_fourier_design_matrix(t, nmodes, F1): +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 F1: fundamental frequency of the Fourier basis + :param f_1: fundamental frequency of the Fourier basis :return: F: fourier design matrix :return: f: Sampling frequencies """ @@ -773,7 +773,7 @@ def create_fourier_design_matrix(t, nmodes, F1): N = len(t) F = np.zeros((N, 2 * nmodes)) - Ffreqs = get_rednoise_freqs(nmodes, F1) + 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]) From c470d1662ea2f444fc8ea798d255af1b14476063 Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Fri, 10 May 2024 14:07:51 +0200 Subject: [PATCH 12/19] get_fundamental_frequency --- src/pint/models/noise_model.py | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/src/pint/models/noise_model.py b/src/pint/models/noise_model.py index a226ef5d4..393be12a0 100644 --- a/src/pint/models/noise_model.py +++ b/src/pint/models/noise_model.py @@ -523,12 +523,12 @@ def get_noise_basis(self, toas): D = (fref.value / freqs.value) ** 2 nf = self.get_pl_vals()[2] - F1 = self.get_F1(t) + f_1 = self.get_fundamental_frequency(t) - Fmat = create_fourier_design_matrix(t, nf, F1) + Fmat = create_fourier_design_matrix(t, nf, f_1) return Fmat * D[:, None] - def get_F1(self, t): + def get_fundamental_frequency(self, t): f_low_factor = 10**self.TNDMFLOW.value return f_low_factor / (t.max() - t.min()) @@ -540,7 +540,7 @@ 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() - F1 = self.get_F1(t) + F1 = self.get_fundamental_frequency(t) Ffreqs = get_rednoise_freqs(nf, F1) return powerlaw(Ffreqs, amp, gam) * Ffreqs[0] @@ -659,7 +659,7 @@ def get_pl_vals(self): amp, gam = self.RNAMP.value / fac, -1 * self.RNIDX.value return (amp, gam, nf) - def get_F1(self, t): + def get_fundamental_frequency(self, t): f_low_factor = 10**self.TNREDFLOW.value return f_low_factor / (t.max() - t.min()) @@ -671,11 +671,11 @@ def get_noise_basis(self, toas): tbl = toas.table t = (tbl["tdbld"].quantity * u.day).to(u.s).value - F1 = self.get_F1(t) + f_1 = self.get_fundamental_frequency(t) nf = self.get_pl_vals()[2] - return create_fourier_design_matrix(t, nf, F1) + return create_fourier_design_matrix(t, nf, f_1) def get_noise_weights(self, toas): """Return power law red noise weights. @@ -686,9 +686,9 @@ def get_noise_weights(self, toas): t = (tbl["tdbld"].quantity * u.day).to(u.s).value amp, gam, nf = self.get_pl_vals() - F1 = self.get_F1(t) + f_1 = self.get_fundamental_frequency(t) - Ffreqs = get_rednoise_freqs(nf, F1) + Ffreqs = get_rednoise_freqs(nf, f_1) return powerlaw(Ffreqs, amp, gam) * Ffreqs[0] def pl_rn_basis_weight_pair(self, toas): From f66b3c6f8205d8947f0ac3190b3cfa719bbd6ac1 Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Fri, 10 May 2024 14:08:24 +0200 Subject: [PATCH 13/19] -- --- src/pint/models/noise_model.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/pint/models/noise_model.py b/src/pint/models/noise_model.py index 393be12a0..ad5f209ad 100644 --- a/src/pint/models/noise_model.py +++ b/src/pint/models/noise_model.py @@ -540,8 +540,8 @@ 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() - F1 = self.get_fundamental_frequency(t) - Ffreqs = get_rednoise_freqs(nf, F1) + f_1 = self.get_fundamental_frequency(t) + Ffreqs = get_rednoise_freqs(nf, f_1) return powerlaw(Ffreqs, amp, gam) * Ffreqs[0] def pl_dm_basis_weight_pair(self, toas): From 5021e5c0f559d137a2d6d16bad57eb2dfafad091 Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Fri, 10 May 2024 14:09:19 +0200 Subject: [PATCH 14/19] amp --- src/pint/models/noise_model.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/pint/models/noise_model.py b/src/pint/models/noise_model.py index ad5f209ad..eb05317fb 100644 --- a/src/pint/models/noise_model.py +++ b/src/pint/models/noise_model.py @@ -474,7 +474,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", ) ) self.add_param( From 0b7a39f81457d070d7124ba6710e3468cc431cea Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Fri, 10 May 2024 14:34:12 +0200 Subject: [PATCH 15/19] powerlaw_corner --- CHANGELOG-unreleased.md | 2 ++ src/pint/models/noise_model.py | 22 +++++++++++++++++++++- 2 files changed, 23 insertions(+), 1 deletion(-) diff --git a/CHANGELOG-unreleased.md b/CHANGELOG-unreleased.md index 2d2000532..341115c37 100644 --- a/CHANGELOG-unreleased.md +++ b/CHANGELOG-unreleased.md @@ -23,7 +23,9 @@ the released changes. - Test for `pint.utils.split_swx()` - Custom type definitions for type hints - `plrednoise_from_wavex()` and `pldmnoise_from_dmwavex()` functions now compute `TNRedFLow` and `TNDMFLow` +- `powerlaw_corner` function ### Fixed - `pint.utils.split_swx()` to use updated `SolarWindDispersionX()` parameter naming convention - Fix #1759 by changing order of comparison ### Removed +- Unnecessary default arguments from the `powerlaw()` function. \ No newline at end of file diff --git a/src/pint/models/noise_model.py b/src/pint/models/noise_model.py index eb05317fb..12784cb1e 100644 --- a/src/pint/models/noise_model.py +++ b/src/pint/models/noise_model.py @@ -781,7 +781,7 @@ def create_fourier_design_matrix(t, nmodes, f_1): return F -def powerlaw(f, A=1e-16, gamma=5): +def powerlaw(f, A, gamma): """Power-law PSD. :param f: Sampling frequencies @@ -791,3 +791,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 + ) From fb1ed740205cc5903c6f7fa5857450879c2206fa Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Fri, 10 May 2024 14:36:08 +0200 Subject: [PATCH 16/19] corner params --- src/pint/models/noise_model.py | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/src/pint/models/noise_model.py b/src/pint/models/noise_model.py index 12784cb1e..d31668d53 100644 --- a/src/pint/models/noise_model.py +++ b/src/pint/models/noise_model.py @@ -502,6 +502,15 @@ def __init__( frozen=True, ) ) + 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, + ) + ) self.covariance_matrix_funcs += [self.pl_dm_cov_matrix] self.basis_funcs += [self.pl_dm_basis_weight_pair] @@ -646,6 +655,15 @@ def __init__( frozen=True, ) ) + 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, + ) + ) self.covariance_matrix_funcs += [self.pl_rn_cov_matrix] self.basis_funcs += [self.pl_rn_basis_weight_pair] From e9a208d4a0e499364b6cfe663699ed9ff0616ce2 Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Fri, 10 May 2024 14:40:46 +0200 Subject: [PATCH 17/19] apply powerlaw_corner --- src/pint/models/noise_model.py | 22 ++++++++++++++++++++-- 1 file changed, 20 insertions(+), 2 deletions(-) diff --git a/src/pint/models/noise_model.py b/src/pint/models/noise_model.py index d31668d53..89f017342 100644 --- a/src/pint/models/noise_model.py +++ b/src/pint/models/noise_model.py @@ -541,6 +541,10 @@ 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. @@ -550,8 +554,13 @@ def get_noise_weights(self, toas): t = (tbl["tdbld"].quantity * u.day).to(u.s).value amp, gam, nf = self.get_pl_vals() 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) * Ffreqs[0] + 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. @@ -681,6 +690,10 @@ 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. @@ -705,9 +718,14 @@ def get_noise_weights(self, toas): amp, gam, nf = self.get_pl_vals() f_1 = self.get_fundamental_frequency(t) + f_c = self.get_fundamental_frequency(t) Ffreqs = get_rednoise_freqs(nf, f_1) - return powerlaw(Ffreqs, amp, gam) * Ffreqs[0] + 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. From 0f9b7154001986f79044ec5d9ac9e4633608e99a Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Mon, 13 May 2024 10:54:24 +0200 Subject: [PATCH 18/19] fix test --- src/pint/models/noise_model.py | 2 +- tests/test_gls_fitter.py | 1 - 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/src/pint/models/noise_model.py b/src/pint/models/noise_model.py index 89f017342..c9ef9b9bd 100644 --- a/src/pint/models/noise_model.py +++ b/src/pint/models/noise_model.py @@ -718,7 +718,7 @@ def get_noise_weights(self, toas): amp, gam, nf = self.get_pl_vals() f_1 = self.get_fundamental_frequency(t) - f_c = self.get_fundamental_frequency(t) + f_c = self.get_corner_frequency(t) Ffreqs = get_rednoise_freqs(nf, f_1) return ( 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 From ed31b69bd83cb1456268c4bf27033b27c376996a Mon Sep 17 00:00:00 2001 From: Abhimanyu Susobhanan Date: Mon, 13 May 2024 11:05:08 +0200 Subject: [PATCH 19/19] CHANGELOG --- CHANGELOG-unreleased.md | 2 ++ 1 file changed, 2 insertions(+) diff --git a/CHANGELOG-unreleased.md b/CHANGELOG-unreleased.md index 341115c37..15296c797 100644 --- a/CHANGELOG-unreleased.md +++ b/CHANGELOG-unreleased.md @@ -24,6 +24,8 @@ the released changes. - Custom type definitions for type hints - `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` ### Fixed - `pint.utils.split_swx()` to use updated `SolarWindDispersionX()` parameter naming convention - Fix #1759 by changing order of comparison