From 8a9eef34561a24c4722d37c35f9cd9f8da2aae08 Mon Sep 17 00:00:00 2001 From: Hazboun6 Date: Thu, 10 Oct 2019 14:47:00 -0700 Subject: [PATCH 01/13] Reorganizing gp bases and priors. --- enterprise/signals/gp_bases.py | 230 ++++++++++++++++++++++++++++++++ enterprise/signals/gp_priors.py | 145 ++++++++++++++++++++ enterprise/signals/utils.py | 218 ++---------------------------- 3 files changed, 385 insertions(+), 208 deletions(-) create mode 100644 enterprise/signals/gp_bases.py create mode 100644 enterprise/signals/gp_priors.py diff --git a/enterprise/signals/gp_bases.py b/enterprise/signals/gp_bases.py new file mode 100644 index 00000000..84e7446a --- /dev/null +++ b/enterprise/signals/gp_bases.py @@ -0,0 +1,230 @@ +# gp_bases.py +"""Utilities module containing various useful +functions for use in other modules. +""" + +from __future__ import (absolute_import, division, + print_function, unicode_literals) + +import numpy as np +from enterprise.signals.parameter import function +###################################### +# Fourier-basis signal functions ##### +###################################### + + +@function +def createfourierdesignmatrix_red(toas, nmodes=30, Tspan=None, + logf=False, fmin=None, fmax=None, + pshift=False, modes=None): + """ + Construct fourier design matrix from eq 11 of Lentati et al, 2013 + :param toas: vector of time series in seconds + :param nmodes: number of fourier coefficients to use + :param freq: option to output frequencies + :param Tspan: option to some other Tspan + :param logf: use log frequency spacing + :param fmin: lower sampling frequency + :param fmax: upper sampling frequency + :param pshift: option to add random phase shift + :param modes: option to provide explicit list or array of + sampling frequencies + + :return: F: fourier design matrix + :return: f: Sampling frequencies + """ + + T = Tspan if Tspan is not None else toas.max() - toas.min() + + # define sampling frequencies + if modes is not None: + nmodes = len(modes) + f = modes + elif fmin is None and fmax is None and not logf: + # make sure partially overlapping sets of modes + # have identical frequencies + f = 1.0 * np.arange(1, nmodes + 1) / T + else: + # more general case + + if fmin is None: + fmin = 1 / T + + if fmax is None: + fmax = nmodes / T + + if logf: + f = np.logspace(np.log10(fmin), np.log10(fmax), nmodes) + else: + f = np.linspace(fmin, fmax, nmodes) + + # add random phase shift to basis functions + ranphase = (np.random.uniform(0.0, 2 * np.pi, nmodes) + if pshift else np.zeros(nmodes)) + + Ffreqs = np.repeat(f, 2) + + N = len(toas) + F = np.zeros((N, 2 * nmodes)) + + # The sine/cosine modes + F[:, ::2] = np.sin(2*np.pi*toas[:, None]*f[None, :] + + ranphase[None, :]) + F[:, 1::2] = np.cos(2*np.pi*toas[:, None]*f[None, :] + + ranphase[None, :]) + + return F, Ffreqs + + +@function +def createfourierdesignmatrix_dm(toas, freqs, nmodes=30, Tspan=None, + pshift=False, fref=1400, logf=False, + fmin=None, fmax=None, modes=None): + """ + Construct DM-variation fourier design matrix. Current + normalization expresses DM signal as a deviation [seconds] + at fref [MHz] + + :param toas: vector of time series in seconds + :param freqs: radio frequencies of observations [MHz] + :param nmodes: number of fourier coefficients to use + :param Tspan: option to some other Tspan + :param pshift: option to add random phase shift + :param fref: reference frequency [MHz] + :param logf: use log frequency spacing + :param fmin: lower sampling frequency + :param fmax: upper sampling frequency + :param modes: option to provide explicit list or array of + sampling frequencies + + :return: F: DM-variation fourier design matrix + :return: f: Sampling frequencies + """ + + # get base fourier design matrix and frequencies + F, Ffreqs = createfourierdesignmatrix_red( + toas, nmodes=nmodes, Tspan=Tspan, logf=logf, + fmin=fmin, fmax=fmax, pshift=pshift, modes=modes) + + # compute the DM-variation vectors + Dm = (fref/freqs)**2 + + return F * Dm[:, None], Ffreqs + + +@function +def createfourierdesignmatrix_env(toas, log10_Amp=-7, log10_Q=np.log10(300), + t0=53000*86400, nmodes=30, Tspan=None, + logf=False, fmin=None, fmax=None, + modes=None): + """ + Construct fourier design matrix with gaussian envelope. + + :param toas: vector of time series in seconds + :param nmodes: number of fourier coefficients to use + :param freqs: radio frequencies of observations [MHz] + :param freq: option to output frequencies + :param Tspan: option to some other Tspan + :param logf: use log frequency spacing + :param fmin: lower sampling frequency + :param fmax: upper sampling frequency + :param log10_Amp: log10 of the Amplitude [s] + :param t0: mean of gaussian envelope [s] + :param log10_Q: log10 of standard deviation of gaussian envelope [days] + :param modes: option to provide explicit list or array of + sampling frequencies + + :return: F: fourier design matrix with gaussian envelope + :return: f: Sampling frequencies + """ + + # get base fourier design matrix and frequencies + F, Ffreqs = createfourierdesignmatrix_red( + toas, nmodes=nmodes, Tspan=Tspan, logf=logf, + fmin=fmin, fmax=fmax, modes=modes) + + # compute gaussian envelope + A = 10**log10_Amp + Q = 10**log10_Q * 86400 + env = A * np.exp(-(toas-t0)**2/2/Q**2) + return F * env[:, None], Ffreqs + + +@function +def createfourierdesignmatrix_ephem(toas, pos, nmodes=30, Tspan=None): + """ + Construct ephemeris perturbation Fourier design matrix and frequencies. + The matrix contains nmodes*6 columns, ordered as by frequency first, + Cartesian coordinate second: + + sin(f0) [x], sin(f0) [y], sin(f0) [z], + cos(f0) [x], cos(f0) [y], cos(f0) [z], + sin(f1) [x], sin(f1) [y], sin(f1) [z], ... + + The corresponding frequency vector repeats every entry six times. + This design matrix should be used with monopole_orf and with + a powerlaw that specifies components=6. + + :param toas: vector of time series in seconds + :param pos: pulsar position as Cartesian vector + :param nmodes: number of Fourier coefficients + :param Tspan: Tspan used to define Fourier bins + + :return: F: Fourier design matrix of shape (len(toas),6*nmodes) + :return: f: Sampling frequencies (6*nmodes) + """ + + F0, F0f = createfourierdesignmatrix_red( + toas, nmodes=nmodes, Tspan=Tspan) + + F1 = np.zeros((len(toas), nmodes, 2, 3), 'd') + F1[:, :, 0, :] = F0[:, 0::2, np.newaxis] + F1[:, :, 1, :] = F0[:, 1::2, np.newaxis] + + # verify this is the scalar product we want + F1 *= pos + + F1f = np.zeros((nmodes, 2, 3), 'd') + F1f[:, :, :] = F0f[::2, np.newaxis, np.newaxis] + + return F1.reshape((len(toas), nmodes*6)), F1f.reshape((nmodes*6, )) + + +def createfourierdesignmatrix_eph(t, nmodes, phi, theta, freq=False, + Tspan=None, logf=False, fmin=None, + fmax=None, modes=None): + raise NotImplementedError( + "createfourierdesignmatrix_eph was removed, " + + "and replaced with createfourierdesignmatrix_ephem") + + +@function +def createfourierdesignmatrix_chromatic(toas, freqs, nmodes=30, Tspan=None, + logf=False, fmin=None, fmax=None, + idx=4): + + """ + Construct Scattering-variation fourier design matrix. + + :param toas: vector of time series in seconds + :param freqs: radio frequencies of observations [MHz] + :param nmodes: number of fourier coefficients to use + :param freq: option to output frequencies + :param Tspan: option to some other Tspan + :param logf: use log frequency spacing + :param fmin: lower sampling frequency + :param fmax: upper sampling frequency + :param idx: Index of chromatic effects + + :return: F: Chromatic-variation fourier design matrix + :return: f: Sampling frequencies + """ + + # get base fourier design matrix and frequencies + F, Ffreqs = createfourierdesignmatrix_red(toas, nmodes=nmodes, Tspan=Tspan, + logf=logf, fmin=fmin, fmax=fmax) + + # compute the DM-variation vectors + Dm = (1400/freqs) ** idx + + return F * Dm[:, None], Ffreqs diff --git a/enterprise/signals/gp_priors.py b/enterprise/signals/gp_priors.py new file mode 100644 index 00000000..76d421d9 --- /dev/null +++ b/enterprise/signals/gp_priors.py @@ -0,0 +1,145 @@ +# gp_priors.py +"""Utilities module containing various useful +functions for use in other modules. +""" + +from __future__ import (absolute_import, division, + print_function, unicode_literals) + +import numpy as np +import scipy.stats + +from enterprise.signals import parameter +from enterprise.signals.parameter import function +import enterprise.constants as const + + +@function +def powerlaw(f, log10_A=-16, gamma=5, components=2): + df = np.diff(np.concatenate((np.array([0]), f[::components]))) + return ((10**log10_A)**2 / 12.0 / np.pi**2 * + const.fyr**(gamma-3) * f**(-gamma) * np.repeat(df, components)) + + +@function +def turnover(f, log10_A=-15, gamma=4.33, lf0=-8.5, kappa=10/3, beta=0.5): + df = np.diff(np.concatenate((np.array([0]), f[::2]))) + hcf = (10**log10_A * (f / const.fyr) ** ((3-gamma) / 2) / + (1 + (10**lf0 / f) ** kappa) ** beta) + return hcf**2/12/np.pi**2/f**3*np.repeat(df, 2) + + +@function +def free_spectrum(f, log10_rho=None): + """ + Free spectral model. PSD amplitude at each frequency + is a free parameter. Model is parameterized by + S(f_i) = \rho_i^2 * T, + where \rho_i is the free parameter and T is the observation + length. + """ + return np.repeat(10**(2*np.array(log10_rho)), 2) + + +@function +def t_process(f, log10_A=-15, gamma=4.33, alphas=None): + """ + t-process model. PSD amplitude at each frequency + is a fuzzy power-law. + """ + alphas = np.ones_like(f) if alphas is None else np.repeat(alphas, 2) + return powerlaw(f, log10_A=log10_A, gamma=gamma) * alphas + + +@function +def t_process_adapt(f, log10_A=-15, gamma=4.33, alphas_adapt=None, nfreq=None): + """ + t-process model. PSD amplitude at each frequency + is a fuzzy power-law. + """ + if alphas_adapt is None: + alpha_model = np.ones_like(f) + else: + if nfreq is None: + alpha_model = np.repeat(alphas_adapt, 2) + else: + alpha_model = np.ones_like(f) + alpha_model[2*int(np.rint(nfreq))] = alphas_adapt + alpha_model[2*int(np.rint(nfreq))+1] = alphas_adapt + + return powerlaw(f, log10_A=log10_A, gamma=gamma) * alpha_model + + +def InvGammaPrior(value, alpha=1, gamma=1): + """Prior function for InvGamma parameters.""" + return scipy.stats.invgamma.pdf(value, alpha, scale=gamma) + + +def InvGammaSampler(alpha=1, gamma=1, size=None): + """Sampling function for Uniform parameters.""" + return scipy.stats.invgamma.rvs(alpha, scale=gamma, size=size) + + +def InvGamma(alpha=1, gamma=1, size=None): + """Class factory for Inverse Gamma parameters.""" + class InvGamma(parameter.Parameter): + _size = size + _prior = parameter.Function(InvGammaPrior, alpha=alpha, gamma=gamma) + _sampler = staticmethod(InvGammaSampler) + _alpha = alpha + _gamma = gamma + + def __repr__(self): + return '"{}": InvGamma({},{})'.format(self.name, alpha, gamma) \ + + ('' if self._size is None else '[{}]'.format(self._size)) + + return InvGamma + + +@function +def turnover_knee(f, log10_A, gamma, lfb, lfk, kappa, delta): + """ + Generic turnover spectrum with a high-frequency knee. + :param f: sampling frequencies of GWB + :param A: characteristic strain amplitude at f=1/yr + :param gamma: negative slope of PSD around f=1/yr (usually 13/3) + :param lfb: log10 transition frequency at which environment dominates GWs + :param lfk: log10 knee frequency due to population finiteness + :param kappa: smoothness of turnover (10/3 for 3-body stellar scattering) + :param delta: slope at higher frequencies + """ + df = np.diff(np.concatenate((np.array([0]), f[::2]))) + hcf = (10**log10_A * (f / const.fyr) ** ((3-gamma) / 2) * + (1.0 + (f / 10**lfk)) ** delta / + np.sqrt(1 + (10**lfb / f) ** kappa)) + return hcf**2 / 12 / np.pi**2 / f**3 * np.repeat(df, 2) + + +@function +def broken_powerlaw(f, log10_A, gamma, delta, log10_fb, kappa=0.1): + """ + Generic broken powerlaw spectrum. + :param f: sampling frequencies + :param A: characteristic strain amplitude [set for gamma at f=1/yr] + :param gamma: negative slope of PSD for f > f_break [set for comparison + at f=1/yr (default 13/3)] + :param delta: slope for frequencies < f_break + :param log10_fb: log10 transition frequency at which slope switches from + gamma to delta + :param kappa: smoothness of transition (Default = 0.1) + """ + df = np.diff(np.concatenate((np.array([0]), f[::2]))) + hcf = (10**log10_A * (f / const.fyr) ** ((3-gamma) / 2) * + (1 + (f / 10**log10_fb) ** (1/kappa)) ** + (kappa * (gamma - delta) / 2)) + return hcf**2 / 12 / np.pi**2 / f**3 * np.repeat(df, 2) + + +@function +def powerlaw_genmodes(f, log10_A=-16, gamma=5, components=2, wgts=None): + if wgts is not None: + df = wgts**2 + else: + df = np.diff(np.concatenate((np.array([0]), f[::components]))) + return ((10**log10_A)**2 / 12.0 / np.pi**2 * + const.fyr**(gamma-3) * f**(-gamma) * np.repeat(df, components)) diff --git a/enterprise/signals/utils.py b/enterprise/signals/utils.py index ec4ffc7b..910019f6 100644 --- a/enterprise/signals/utils.py +++ b/enterprise/signals/utils.py @@ -1,5 +1,6 @@ -#utils.py -"""Utilities module containing various useful +# utils.py +""" +Utilities module containing various useful functions for use in other modules. """ @@ -19,6 +20,12 @@ import enterprise import enterprise.constants as const from enterprise.signals.parameter import function +from enterprise.signals.gp_priors import powerlaw, turnover +from enterprise.signals.gp_bases import (createfourierdesignmatrix_red, + createfourierdesignmatrix_dm, + createfourierdesignmatrix_env, + createfourierdesignmatrix_ephem, + createfourierdesignmatrix_eph) try: from sksparse.cholmod import cholesky @@ -226,196 +233,6 @@ def create_stabletimingdesignmatrix(designmat, fastDesign=True): return Mm -###################################### -# Fourier-basis signal functions ##### -###################################### - - -@function -def createfourierdesignmatrix_red(toas, nmodes=30, Tspan=None, - logf=False, fmin=None, fmax=None, - pshift=False, modes=None): - """ - Construct fourier design matrix from eq 11 of Lentati et al, 2013 - :param toas: vector of time series in seconds - :param nmodes: number of fourier coefficients to use - :param freq: option to output frequencies - :param Tspan: option to some other Tspan - :param logf: use log frequency spacing - :param fmin: lower sampling frequency - :param fmax: upper sampling frequency - :param pshift: option to add random phase shift - :param modes: option to provide explicit list or array of - sampling frequencies - - :return: F: fourier design matrix - :return: f: Sampling frequencies - """ - - T = Tspan if Tspan is not None else toas.max() - toas.min() - - # define sampling frequencies - if modes is not None: - nmodes = len(modes) - f = modes - elif fmin is None and fmax is None and not logf: - # make sure partially overlapping sets of modes - # have identical frequencies - f = 1.0 * np.arange(1, nmodes + 1) / T - else: - # more general case - - if fmin is None: - fmin = 1 / T - - if fmax is None: - fmax = nmodes / T - - if logf: - f = np.logspace(np.log10(fmin), np.log10(fmax), nmodes) - else: - f = np.linspace(fmin, fmax, nmodes) - - # add random phase shift to basis functions - ranphase = (np.random.uniform(0.0, 2 * np.pi, nmodes) - if pshift else np.zeros(nmodes)) - - Ffreqs = np.repeat(f, 2) - - N = len(toas) - F = np.zeros((N, 2 * nmodes)) - - # The sine/cosine modes - F[:,::2] = np.sin(2*np.pi*toas[:,None]*f[None,:] + - ranphase[None,:]) - F[:,1::2] = np.cos(2*np.pi*toas[:,None]*f[None,:] + - ranphase[None,:]) - - return F, Ffreqs - - -@function -def createfourierdesignmatrix_dm(toas, freqs, nmodes=30, Tspan=None, - pshift=False, fref=1400, logf=False, - fmin=None, fmax=None, modes=None): - """ - Construct DM-variation fourier design matrix. Current - normalization expresses DM signal as a deviation [seconds] - at fref [MHz] - - :param toas: vector of time series in seconds - :param freqs: radio frequencies of observations [MHz] - :param nmodes: number of fourier coefficients to use - :param Tspan: option to some other Tspan - :param pshift: option to add random phase shift - :param fref: reference frequency [MHz] - :param logf: use log frequency spacing - :param fmin: lower sampling frequency - :param fmax: upper sampling frequency - :param modes: option to provide explicit list or array of - sampling frequencies - - :return: F: DM-variation fourier design matrix - :return: f: Sampling frequencies - """ - - # get base fourier design matrix and frequencies - F, Ffreqs = createfourierdesignmatrix_red( - toas, nmodes=nmodes, Tspan=Tspan, logf=logf, - fmin=fmin, fmax=fmax, pshift=pshift, modes=modes) - - # compute the DM-variation vectors - Dm = (fref/freqs)**2 - - return F * Dm[:, None], Ffreqs - - -@function -def createfourierdesignmatrix_env(toas, log10_Amp=-7, log10_Q=np.log10(300), - t0=53000*86400, nmodes=30, Tspan=None, - logf=False, fmin=None, fmax=None, - modes=None): - """ - Construct fourier design matrix with gaussian envelope. - - :param toas: vector of time series in seconds - :param nmodes: number of fourier coefficients to use - :param freqs: radio frequencies of observations [MHz] - :param freq: option to output frequencies - :param Tspan: option to some other Tspan - :param logf: use log frequency spacing - :param fmin: lower sampling frequency - :param fmax: upper sampling frequency - :param log10_Amp: log10 of the Amplitude [s] - :param t0: mean of gaussian envelope [s] - :param log10_Q: log10 of standard deviation of gaussian envelope [days] - :param modes: option to provide explicit list or array of - sampling frequencies - - :return: F: fourier design matrix with gaussian envelope - :return: f: Sampling frequencies - """ - - # get base fourier design matrix and frequencies - F, Ffreqs = createfourierdesignmatrix_red( - toas, nmodes=nmodes, Tspan=Tspan, logf=logf, - fmin=fmin, fmax=fmax, modes=modes) - - # compute gaussian envelope - A = 10**log10_Amp - Q = 10**log10_Q * 86400 - env = A * np.exp(-(toas-t0)**2/2/Q**2) - return F * env[:, None], Ffreqs - - -@function -def createfourierdesignmatrix_ephem(toas, pos, nmodes=30, Tspan=None): - """ - Construct ephemeris perturbation Fourier design matrix and frequencies. - The matrix contains nmodes*6 columns, ordered as by frequency first, - Cartesian coordinate second: - - sin(f0) [x], sin(f0) [y], sin(f0) [z], - cos(f0) [x], cos(f0) [y], cos(f0) [z], - sin(f1) [x], sin(f1) [y], sin(f1) [z], ... - - The corresponding frequency vector repeats every entry six times. - This design matrix should be used with monopole_orf and with - a powerlaw that specifies components=6. - - :param toas: vector of time series in seconds - :param pos: pulsar position as Cartesian vector - :param nmodes: number of Fourier coefficients - :param Tspan: Tspan used to define Fourier bins - - :return: F: Fourier design matrix of shape (len(toas),6*nmodes) - :return: f: Sampling frequencies (6*nmodes) - """ - - F0, F0f = createfourierdesignmatrix_red( - toas, nmodes=nmodes, Tspan=Tspan) - - F1 = np.zeros((len(toas),nmodes,2,3), 'd') - F1[:,:,0,:] = F0[:,0::2,np.newaxis] - F1[:,:,1,:] = F0[:,1::2,np.newaxis] - - # verify this is the scalar product we want - F1 *= pos - - F1f = np.zeros((nmodes,2,3), 'd') - F1f[:,:,:] = F0f[::2,np.newaxis,np.newaxis] - - return F1.reshape((len(toas),nmodes*6)), F1f.reshape((nmodes*6,)) - - -def createfourierdesignmatrix_eph(t, nmodes, phi, theta, freq=False, - Tspan=None, logf=False, fmin=None, - fmax=None, modes=None): - raise NotImplementedError( - "createfourierdesignmatrix_eph was removed, " + - "and replaced with createfourierdesignmatrix_ephem") - - ################################### # Deterministic GW signal functions ################################### @@ -921,21 +738,6 @@ def linear_interp_basis(toas, dt=30*86400): return M[:, idx], x[idx] -@function -def powerlaw(f, log10_A=-16, gamma=5, components=2): - df = np.diff(np.concatenate((np.array([0]), f[::components]))) - return ((10**log10_A)**2 / 12.0 / np.pi**2 * - const.fyr**(gamma-3) * f**(-gamma) * np.repeat(df, components)) - - -@function -def turnover(f, log10_A=-15, gamma=4.33, lf0=-8.5, kappa=10/3, beta=0.5): - df = np.diff(np.concatenate((np.array([0]), f[::2]))) - hcf = (10**log10_A * (f / const.fyr) ** ((3-gamma) / 2) / - (1 + (10**lf0 / f) ** kappa) ** beta) - return hcf**2/12/np.pi**2/f**3*np.repeat(df, 2) - - # overlap reduction functions @function @@ -1170,7 +972,7 @@ def createfourierdesignmatrix_physicalephem(toas, planetssb, pos_t, c = np.zeros(6) c[i] = dpar - #Fl.append(physical_ephem_delay(toas, planetssb, pos_t, + # Fl.append(physical_ephem_delay(toas, planetssb, pos_t, # **{parname: c}, **oa)/dpar) kwarg_dict = {parname: c} kwarg_dict.update(oa) From be495c659a7039b90b810c97689b3ebebff78c91 Mon Sep 17 00:00:00 2001 From: Hazboun6 Date: Thu, 10 Oct 2019 16:08:29 -0700 Subject: [PATCH 02/13] Add lint exceptions. --- enterprise/signals/gp_bases.py | 5 +++++ enterprise/signals/utils.py | 9 +++------ 2 files changed, 8 insertions(+), 6 deletions(-) diff --git a/enterprise/signals/gp_bases.py b/enterprise/signals/gp_bases.py index 84e7446a..d7b2598d 100644 --- a/enterprise/signals/gp_bases.py +++ b/enterprise/signals/gp_bases.py @@ -12,6 +12,11 @@ # Fourier-basis signal functions ##### ###################################### +__all__ = ['createfourierdesignmatrix_red', + 'createfourierdesignmatrix_dm', + 'createfourierdesignmatrix_env', + 'createfourierdesignmatrix_ephem', + 'createfourierdesignmatrix_eph'] @function def createfourierdesignmatrix_red(toas, nmodes=30, Tspan=None, diff --git a/enterprise/signals/utils.py b/enterprise/signals/utils.py index 910019f6..fb6ce299 100644 --- a/enterprise/signals/utils.py +++ b/enterprise/signals/utils.py @@ -20,12 +20,9 @@ import enterprise import enterprise.constants as const from enterprise.signals.parameter import function -from enterprise.signals.gp_priors import powerlaw, turnover -from enterprise.signals.gp_bases import (createfourierdesignmatrix_red, - createfourierdesignmatrix_dm, - createfourierdesignmatrix_env, - createfourierdesignmatrix_ephem, - createfourierdesignmatrix_eph) +from enterprise.signals.gp_priors import powerlaw, turnover # noqa: F401 +from enterprise.signals.gp_bases import * # noqa: F401, F403 + try: from sksparse.cholmod import cholesky From 020dcf884c425d6b138ce297817e51e0d2ebef07 Mon Sep 17 00:00:00 2001 From: Hazboun6 Date: Thu, 10 Oct 2019 20:06:27 -0700 Subject: [PATCH 03/13] Fix lint. --- enterprise/signals/gp_bases.py | 1 + 1 file changed, 1 insertion(+) diff --git a/enterprise/signals/gp_bases.py b/enterprise/signals/gp_bases.py index d7b2598d..63559123 100644 --- a/enterprise/signals/gp_bases.py +++ b/enterprise/signals/gp_bases.py @@ -18,6 +18,7 @@ 'createfourierdesignmatrix_ephem', 'createfourierdesignmatrix_eph'] + @function def createfourierdesignmatrix_red(toas, nmodes=30, Tspan=None, logf=False, fmin=None, fmax=None, From 838483652df82ef46de96c88b2cd91e256adf45d Mon Sep 17 00:00:00 2001 From: Hazboun6 Date: Tue, 22 Oct 2019 20:41:11 -0700 Subject: [PATCH 04/13] Made some import and lint changes. --- enterprise/signals/utils.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/enterprise/signals/utils.py b/enterprise/signals/utils.py index fb6ce299..cce28777 100644 --- a/enterprise/signals/utils.py +++ b/enterprise/signals/utils.py @@ -21,7 +21,13 @@ import enterprise.constants as const from enterprise.signals.parameter import function from enterprise.signals.gp_priors import powerlaw, turnover # noqa: F401 -from enterprise.signals.gp_bases import * # noqa: F401, F403 +import enterprise.signals as sigs # noqa: F401 +from enterprise.signals.gp_bases import ( # noqa: F401 + createfourierdesignmatrix_red, + createfourierdesignmatrix_dm, + createfourierdesignmatrix_env, + createfourierdesignmatrix_ephem, + createfourierdesignmatrix_eph) try: From dee22b749964a7afcdfa4097d875e28f00f600fd Mon Sep 17 00:00:00 2001 From: Hazboun6 Date: Tue, 22 Oct 2019 20:41:48 -0700 Subject: [PATCH 05/13] Adding tests for new priors. --- tests/test_gp_priors.py | 317 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 317 insertions(+) create mode 100644 tests/test_gp_priors.py diff --git a/tests/test_gp_priors.py b/tests/test_gp_priors.py new file mode 100644 index 00000000..666745df --- /dev/null +++ b/tests/test_gp_priors.py @@ -0,0 +1,317 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- + +""" +test_gp_priors +---------------------------------- + +Tests for GP priors and bases. +""" + + +import unittest +import numpy as np + +from tests.enterprise_test_data import datadir +from enterprise.pulsar import Pulsar +from enterprise.signals import parameter +from enterprise.signals import gp_signals +from enterprise.signals import gp_priors +from enterprise.signals import gp_bases +import scipy.stats + + +class TestGPSignals(unittest.TestCase): + + @classmethod + def setUpClass(cls): + """Setup the Pulsar object.""" + + # initialize Pulsar class + cls.psr = Pulsar(datadir + '/B1855+09_NANOGrav_9yv1.gls.par', + datadir + '/B1855+09_NANOGrav_9yv1.tim') + + def test_turnover_prior(self): + """Test that red noise signal returns correct values.""" + # set up signal parameter + pr = gp_priors.turnover(log10_A=parameter.Uniform(-18,-12), + gamma=parameter.Uniform(1,7), + lf0=parameter.Uniform(-9,-7.5), + kappa=parameter.Uniform(2.5,5), + beta=parameter.Uniform(0.01,1),) + basis = gp_bases.createfourierdesignmatrix_red(nmodes=30) + rn = gp_signals.BasisGP(priorFunction=pr, basisFunction=basis, + name='red_noise') + rnm = rn(self.psr) + + # parameters + log10_A, gamma, lf0, kappa, beta = -14.5, 4.33, -8.5, 3, 0.5 + params = {'B1855+09_red_noise_log10_A': log10_A, + 'B1855+09_red_noise_gamma': gamma, + 'B1855+09_red_noise_lf0': lf0, + 'B1855+09_red_noise_kappa': kappa, + 'B1855+09_red_noise_beta': beta} + + # basis matrix test + F, f2 = gp_bases.createfourierdesignmatrix_red(self.psr.toas, + nmodes=30) + msg = 'F matrix incorrect for turnover.' + assert np.allclose(F, rnm.get_basis(params)), msg + + # spectrum test + phi = gp_priors.turnover(f2, + log10_A=log10_A, + gamma=gamma, + lf0=lf0, + kappa=kappa, + beta=beta) + msg = 'Spectrum incorrect for turnover.' + assert np.all(rnm.get_phi(params) == phi), msg + + # inverse spectrum test + msg = 'Spectrum inverse incorrect for turnover.' + assert np.all(rnm.get_phiinv(params) == 1/phi), msg + + # test shape + msg = 'F matrix shape incorrect' + assert rnm.get_basis(params).shape == F.shape, msg + + def test_free_spec_prior(self): + """Test that red noise signal returns correct values.""" + # set up signal parameter + pr = gp_priors.free_spectrum(log10_rho=parameter.Uniform(-10, -4, + size=30)) + basis = gp_bases.createfourierdesignmatrix_red(nmodes=30) + rn = gp_signals.BasisGP(priorFunction=pr, basisFunction=basis, + name='red_noise') + rnm = rn(self.psr) + # parameters + rhos = np.random.uniform(-10, -4, size=30) + params = {'B1855+09_red_noise_log10_rho': rhos} + # basis matrix test + F, f2 = gp_bases.createfourierdesignmatrix_red(self.psr.toas, + nmodes=30) + msg = 'F matrix incorrect for free spectrum.' + assert np.allclose(F, rnm.get_basis(params)), msg + + # spectrum test + phi = gp_priors.free_spectrum(f2, log10_rho=rhos) + + msg = 'Spectrum incorrect for free spectrum.' + assert np.all(rnm.get_phi(params) == phi), msg + + # inverse spectrum test + msg = 'Spectrum inverse incorrect for free spectrum.' + assert np.all(rnm.get_phiinv(params) == 1/phi), msg + + # test shape + msg = 'F matrix shape incorrect' + assert rnm.get_basis(params).shape == F.shape, msg + + def test_t_process_prior(self): + """Test that red noise signal returns correct values.""" + # set up signal parameter + pr = gp_priors.t_process(log10_A=parameter.Uniform(-18,-12), + gamma=parameter.Uniform(1,7), + alphas=gp_priors.InvGamma(alpha=1, gamma=1, + size=30)) + basis = gp_bases.createfourierdesignmatrix_red(nmodes=30) + rn = gp_signals.BasisGP(priorFunction=pr, basisFunction=basis, + name='red_noise') + rnm = rn(self.psr) + # parameters + alphas = scipy.stats.invgamma.rvs(1, scale=1, size=30) + log10_A, gamma = -15, 4.33 + params = {'B1855+09_red_noise_log10_A': log10_A, + 'B1855+09_red_noise_gamma': gamma, + 'B1855+09_red_noise_alphas': alphas} + # basis matrix test + F, f2 = gp_bases.createfourierdesignmatrix_red(self.psr.toas, + nmodes=30) + msg = 'F matrix incorrect for free spectrum.' + assert np.allclose(F, rnm.get_basis(params)), msg + + # spectrum test + phi = gp_priors.t_process(f2, log10_A=log10_A, gamma=gamma, + alphas=alphas) + + msg = 'Spectrum incorrect for free spectrum.' + assert np.all(rnm.get_phi(params) == phi), msg + + # inverse spectrum test + msg = 'Spectrum inverse incorrect for free spectrum.' + assert np.all(rnm.get_phiinv(params) == 1/phi), msg + + # test shape + msg = 'F matrix shape incorrect' + assert rnm.get_basis(params).shape == F.shape, msg + + def test_adapt_t_process_prior(self): + """Test that red noise signal returns correct values.""" + # set up signal parameter + pr = gp_priors.t_process_adapt(log10_A=parameter.Uniform(-18,-12), + gamma=parameter.Uniform(1,7), + alphas_adapt=gp_priors.InvGamma(), + nfreq=parameter.Uniform(5,25)) + basis = gp_bases.createfourierdesignmatrix_red(nmodes=30) + rn = gp_signals.BasisGP(priorFunction=pr, basisFunction=basis, + name='red_noise') + rnm = rn(self.psr) + # parameters + alphas = scipy.stats.invgamma.rvs(1, scale=1, size=1) + log10_A, gamma, nfreq = -15, 4.33, 12 + params = {'B1855+09_red_noise_log10_A': log10_A, + 'B1855+09_red_noise_gamma': gamma, + 'B1855+09_red_noise_alphas_adapt': alphas, + 'B1855+09_red_noise_nfreq': nfreq} + # basis matrix test + F, f2 = gp_bases.createfourierdesignmatrix_red(self.psr.toas, + nmodes=30) + msg = 'F matrix incorrect for free spectrum.' + assert np.allclose(F, rnm.get_basis(params)), msg + + # spectrum test + phi = gp_priors.t_process_adapt(f2, log10_A=log10_A, gamma=gamma, + alphas_adapt=alphas, nfreq=nfreq) + + msg = 'Spectrum incorrect for free spectrum.' + assert np.all(rnm.get_phi(params) == phi), msg + + # inverse spectrum test + msg = 'Spectrum inverse incorrect for free spectrum.' + assert np.all(rnm.get_phiinv(params) == 1/phi), msg + + # test shape + msg = 'F matrix shape incorrect' + assert rnm.get_basis(params).shape == F.shape, msg + + def test_turnover_knee_prior(self): + """Test that red noise signal returns correct values.""" + # set up signal parameter + pr = gp_priors.turnover_knee(log10_A=parameter.Uniform(-18,-12), + gamma=parameter.Uniform(1,7), + lfb=parameter.Uniform(-9,-7.5), + lfk=parameter.Uniform(-9,-7.5), + kappa=parameter.Uniform(2.5,5), + delta=parameter.Uniform(0.01,1),) + basis = gp_bases.createfourierdesignmatrix_red(nmodes=30) + rn = gp_signals.BasisGP(priorFunction=pr, basisFunction=basis, + name='red_noise') + rnm = rn(self.psr) + + # parameters + log10_A, gamma, lfb = -14.5, 4.33, -8.5 + lfk, kappa, delta = -8.5, 3, 0.5 + params = {'B1855+09_red_noise_log10_A': log10_A, + 'B1855+09_red_noise_gamma': gamma, + 'B1855+09_red_noise_lfb': lfb, + 'B1855+09_red_noise_lfk': lfk, + 'B1855+09_red_noise_kappa': kappa, + 'B1855+09_red_noise_delta': delta} + + # basis matrix test + F, f2 = gp_bases.createfourierdesignmatrix_red(self.psr.toas, + nmodes=30) + msg = 'F matrix incorrect for turnover.' + assert np.allclose(F, rnm.get_basis(params)), msg + + # spectrum test + phi = gp_priors.turnover_knee(f2, + log10_A=log10_A, + gamma=gamma, + lfb=lfb, + lfk=lfk, + kappa=kappa, + delta=delta) + msg = 'Spectrum incorrect for turnover.' + assert np.all(rnm.get_phi(params) == phi), msg + + # inverse spectrum test + msg = 'Spectrum inverse incorrect for turnover.' + assert np.all(rnm.get_phiinv(params) == 1/phi), msg + + # test shape + msg = 'F matrix shape incorrect' + assert rnm.get_basis(params).shape == F.shape, msg + + def test_broken_powerlaw_prior(self): + """Test that red noise signal returns correct values.""" + # set up signal parameter + pr = gp_priors.broken_powerlaw(log10_A=parameter.Uniform(-18,-12), + gamma=parameter.Uniform(1,7), + log10_fb=parameter.Uniform(-9,-7.5), + kappa=parameter.Uniform(0.1,1.0), + delta=parameter.Uniform(0.01,1),) + basis = gp_bases.createfourierdesignmatrix_red(nmodes=30) + rn = gp_signals.BasisGP(priorFunction=pr, basisFunction=basis, + name='red_noise') + rnm = rn(self.psr) + + # parameters + log10_A, gamma, log10_fb, kappa, delta = -14.5, 4.33, -8.5, 1, 0.5 + params = {'B1855+09_red_noise_log10_A': log10_A, + 'B1855+09_red_noise_gamma': gamma, + 'B1855+09_red_noise_log10_fb': log10_fb, + 'B1855+09_red_noise_kappa': kappa, + 'B1855+09_red_noise_delta': delta} + + # basis matrix test + F, f2 = gp_bases.createfourierdesignmatrix_red(self.psr.toas, + nmodes=30) + msg = 'F matrix incorrect for turnover.' + assert np.allclose(F, rnm.get_basis(params)), msg + + # spectrum test + phi = gp_priors.broken_powerlaw(f2, + log10_A=log10_A, + gamma=gamma, + log10_fb=log10_fb, + kappa=kappa, + delta=delta) + msg = 'Spectrum incorrect for turnover.' + assert np.all(rnm.get_phi(params) == phi), msg + + # inverse spectrum test + msg = 'Spectrum inverse incorrect for turnover.' + assert np.all(rnm.get_phiinv(params) == 1/phi), msg + + # test shape + msg = 'F matrix shape incorrect' + assert rnm.get_basis(params).shape == F.shape, msg + + def test_powerlaw_genmodes_prior(self): + """Test that red noise signal returns correct values.""" + # set up signal parameter + pr = gp_priors.powerlaw_genmodes(log10_A=parameter.Uniform(-18,-12), + gamma=parameter.Uniform(1,7)) + basis = gp_bases.createfourierdesignmatrix_chromatic(nmodes=30) + rn = gp_signals.BasisGP(priorFunction=pr, basisFunction=basis, + name='red_noise') + rnm = rn(self.psr) + + # parameters + log10_A, gamma = -14.5, 4.33 + params = {'B1855+09_red_noise_log10_A': log10_A, + 'B1855+09_red_noise_gamma': gamma} + + # basis matrix test + F, f2 = gp_bases.createfourierdesignmatrix_chromatic(self.psr.toas, + self.psr.freqs, + nmodes=30) + msg = 'F matrix incorrect for turnover.' + assert np.allclose(F, rnm.get_basis(params)), msg + + # spectrum test + phi = gp_priors.powerlaw_genmodes(f2, + log10_A=log10_A, + gamma=gamma) + msg = 'Spectrum incorrect for turnover.' + assert np.all(rnm.get_phi(params) == phi), msg + + # inverse spectrum test + msg = 'Spectrum inverse incorrect for turnover.' + assert np.all(rnm.get_phiinv(params) == 1/phi), msg + + # test shape + msg = 'F matrix shape incorrect' + assert rnm.get_basis(params).shape == F.shape, msg From 9170a7ddb2080bb16df947c391b3ec3389825250 Mon Sep 17 00:00:00 2001 From: Hazboun6 Date: Thu, 27 Feb 2020 15:36:57 -0800 Subject: [PATCH 06/13] Added sunssb to TempoPulsar and tests --- enterprise/pulsar.py | 24 ++++++++++++++++++++++++ tests/test_pulsar.py | 8 ++++++++ 2 files changed, 32 insertions(+) diff --git a/enterprise/pulsar.py b/enterprise/pulsar.py index 8aa29a08..6fa1d534 100644 --- a/enterprise/pulsar.py +++ b/enterprise/pulsar.py @@ -280,6 +280,11 @@ def planetssb(self): """Return planetary position vectors at all timestamps""" return self._planetssb[self._isort, :, :] + @property + def sunssb(self): + """Return sun position vector at all timestamps""" + return self._sunssb[self._isort, :] + class PintPulsar(BasePulsar): @@ -378,6 +383,7 @@ def __init__(self, t2pulsar, sort=True, self._raj, self._decj = self._get_radec(t2pulsar) self._pos = self._get_pos() self._planetssb = self._get_planetssb(t2pulsar) + self._sunssb = self._get_sunssb(t2pulsar) # gather DM/DMX information if available self._set_dm(t2pulsar) @@ -447,6 +453,24 @@ def _get_planetssb(self, t2pulsar): planetssb[:,ii,3:] = utils.ecl2eq_vec(planetssb[:,ii,3:]) return planetssb + def _get_sunssb(self, t2pulsar): + sunssb = None + if self.planets: + # for ii in range(1, 10): + # tag = 'DMASSPLANET' + str(ii) + # self.t2pulsar[tag].val = 0.0 + self.t2pulsar.formbats() + sunssb = np.zeros((len(self._toas), 6)) + sunssb[:, :] = self.t2pulsar.sun_ssb + + + if 'ELONG' and 'ELAT' in np.concatenate((t2pulsar.pars(), + t2pulsar.pars( + which='set'))): + sunssb[:,:3] = utils.ecl2eq_vec(sunssb[:,:3]) + sunssb[:,3:] = utils.ecl2eq_vec(sunssb[:,3:]) + return sunssb + def Pulsar(*args, **kwargs): diff --git a/tests/test_pulsar.py b/tests/test_pulsar.py index 7f510da3..0eeb4475 100644 --- a/tests/test_pulsar.py +++ b/tests/test_pulsar.py @@ -104,6 +104,14 @@ def test_filter_data(self): """Place holder for filter_data tests.""" assert self.psr.filter_data() is None + def test_planetssb(self): + """Place holder for filter_data tests.""" + assert hasattr(self.psr, 'planetssb') + + def test_sunssb(self): + """Place holder for filter_data tests.""" + assert hasattr(self.psr, 'sunssb') + def test_to_pickle(self): """Place holder for to_pickle tests.""" self.psr.to_pickle() From ceb8962b6e0b7adbf2635cc422dc05a237559928 Mon Sep 17 00:00:00 2001 From: Hazboun6 Date: Thu, 27 Feb 2020 15:45:54 -0800 Subject: [PATCH 07/13] Merged sunssb changes 2 --- enterprise/signals/gp_bases.py | 104 ------------- enterprise/signals/gp_priors.py | 68 --------- tests/test_gp_priors.py | 250 -------------------------------- 3 files changed, 422 deletions(-) diff --git a/enterprise/signals/gp_bases.py b/enterprise/signals/gp_bases.py index cce9f5e1..63559123 100644 --- a/enterprise/signals/gp_bases.py +++ b/enterprise/signals/gp_bases.py @@ -3,24 +3,15 @@ functions for use in other modules. """ -<<<<<<< HEAD from __future__ import (absolute_import, division, print_function, unicode_literals) import numpy as np from enterprise.signals.parameter import function -======= -from __future__ import absolute_import, division, print_function, unicode_literals - -import numpy as np -from enterprise.signals.parameter import function - ->>>>>>> 7e6931cbc5a819ee0bf223984eef09d963644927 ###################################### # Fourier-basis signal functions ##### ###################################### -<<<<<<< HEAD __all__ = ['createfourierdesignmatrix_red', 'createfourierdesignmatrix_dm', 'createfourierdesignmatrix_env', @@ -32,21 +23,6 @@ def createfourierdesignmatrix_red(toas, nmodes=30, Tspan=None, logf=False, fmin=None, fmax=None, pshift=False, modes=None): -======= -__all__ = [ - "createfourierdesignmatrix_red", - "createfourierdesignmatrix_dm", - "createfourierdesignmatrix_env", - "createfourierdesignmatrix_ephem", - "createfourierdesignmatrix_eph", -] - - -@function -def createfourierdesignmatrix_red( - toas, nmodes=30, Tspan=None, logf=False, fmin=None, fmax=None, pshift=False, modes=None -): ->>>>>>> 7e6931cbc5a819ee0bf223984eef09d963644927 """ Construct fourier design matrix from eq 11 of Lentati et al, 2013 :param toas: vector of time series in seconds @@ -89,12 +65,8 @@ def createfourierdesignmatrix_red( f = np.linspace(fmin, fmax, nmodes) # add random phase shift to basis functions -<<<<<<< HEAD ranphase = (np.random.uniform(0.0, 2 * np.pi, nmodes) if pshift else np.zeros(nmodes)) -======= - ranphase = np.random.uniform(0.0, 2 * np.pi, nmodes) if pshift else np.zeros(nmodes) ->>>>>>> 7e6931cbc5a819ee0bf223984eef09d963644927 Ffreqs = np.repeat(f, 2) @@ -102,29 +74,18 @@ def createfourierdesignmatrix_red( F = np.zeros((N, 2 * nmodes)) # The sine/cosine modes -<<<<<<< HEAD F[:, ::2] = np.sin(2*np.pi*toas[:, None]*f[None, :] + ranphase[None, :]) F[:, 1::2] = np.cos(2*np.pi*toas[:, None]*f[None, :] + ranphase[None, :]) -======= - F[:, ::2] = np.sin(2 * np.pi * toas[:, None] * f[None, :] + ranphase[None, :]) - F[:, 1::2] = np.cos(2 * np.pi * toas[:, None] * f[None, :] + ranphase[None, :]) ->>>>>>> 7e6931cbc5a819ee0bf223984eef09d963644927 return F, Ffreqs @function -<<<<<<< HEAD def createfourierdesignmatrix_dm(toas, freqs, nmodes=30, Tspan=None, pshift=False, fref=1400, logf=False, fmin=None, fmax=None, modes=None): -======= -def createfourierdesignmatrix_dm( - toas, freqs, nmodes=30, Tspan=None, pshift=False, fref=1400, logf=False, fmin=None, fmax=None, modes=None -): ->>>>>>> 7e6931cbc5a819ee0bf223984eef09d963644927 """ Construct DM-variation fourier design matrix. Current normalization expresses DM signal as a deviation [seconds] @@ -148,43 +109,20 @@ def createfourierdesignmatrix_dm( # get base fourier design matrix and frequencies F, Ffreqs = createfourierdesignmatrix_red( -<<<<<<< HEAD toas, nmodes=nmodes, Tspan=Tspan, logf=logf, fmin=fmin, fmax=fmax, pshift=pshift, modes=modes) # compute the DM-variation vectors Dm = (fref/freqs)**2 -======= - toas, nmodes=nmodes, Tspan=Tspan, logf=logf, fmin=fmin, fmax=fmax, pshift=pshift, modes=modes - ) - - # compute the DM-variation vectors - Dm = (fref / freqs) ** 2 ->>>>>>> 7e6931cbc5a819ee0bf223984eef09d963644927 return F * Dm[:, None], Ffreqs @function -<<<<<<< HEAD def createfourierdesignmatrix_env(toas, log10_Amp=-7, log10_Q=np.log10(300), t0=53000*86400, nmodes=30, Tspan=None, logf=False, fmin=None, fmax=None, modes=None): -======= -def createfourierdesignmatrix_env( - toas, - log10_Amp=-7, - log10_Q=np.log10(300), - t0=53000 * 86400, - nmodes=30, - Tspan=None, - logf=False, - fmin=None, - fmax=None, - modes=None, -): ->>>>>>> 7e6931cbc5a819ee0bf223984eef09d963644927 """ Construct fourier design matrix with gaussian envelope. @@ -208,7 +146,6 @@ def createfourierdesignmatrix_env( # get base fourier design matrix and frequencies F, Ffreqs = createfourierdesignmatrix_red( -<<<<<<< HEAD toas, nmodes=nmodes, Tspan=Tspan, logf=logf, fmin=fmin, fmax=fmax, modes=modes) @@ -216,15 +153,6 @@ def createfourierdesignmatrix_env( A = 10**log10_Amp Q = 10**log10_Q * 86400 env = A * np.exp(-(toas-t0)**2/2/Q**2) -======= - toas, nmodes=nmodes, Tspan=Tspan, logf=logf, fmin=fmin, fmax=fmax, modes=modes - ) - - # compute gaussian envelope - A = 10 ** log10_Amp - Q = 10 ** log10_Q * 86400 - env = A * np.exp(-((toas - t0) ** 2) / 2 / Q ** 2) ->>>>>>> 7e6931cbc5a819ee0bf223984eef09d963644927 return F * env[:, None], Ffreqs @@ -252,23 +180,16 @@ def createfourierdesignmatrix_ephem(toas, pos, nmodes=30, Tspan=None): :return: f: Sampling frequencies (6*nmodes) """ -<<<<<<< HEAD F0, F0f = createfourierdesignmatrix_red( toas, nmodes=nmodes, Tspan=Tspan) F1 = np.zeros((len(toas), nmodes, 2, 3), 'd') -======= - F0, F0f = createfourierdesignmatrix_red(toas, nmodes=nmodes, Tspan=Tspan) - - F1 = np.zeros((len(toas), nmodes, 2, 3), "d") ->>>>>>> 7e6931cbc5a819ee0bf223984eef09d963644927 F1[:, :, 0, :] = F0[:, 0::2, np.newaxis] F1[:, :, 1, :] = F0[:, 1::2, np.newaxis] # verify this is the scalar product we want F1 *= pos -<<<<<<< HEAD F1f = np.zeros((nmodes, 2, 3), 'd') F1f[:, :, :] = F0f[::2, np.newaxis, np.newaxis] @@ -287,24 +208,6 @@ def createfourierdesignmatrix_eph(t, nmodes, phi, theta, freq=False, def createfourierdesignmatrix_chromatic(toas, freqs, nmodes=30, Tspan=None, logf=False, fmin=None, fmax=None, idx=4): -======= - F1f = np.zeros((nmodes, 2, 3), "d") - F1f[:, :, :] = F0f[::2, np.newaxis, np.newaxis] - - return F1.reshape((len(toas), nmodes * 6)), F1f.reshape((nmodes * 6,)) - - -def createfourierdesignmatrix_eph( - t, nmodes, phi, theta, freq=False, Tspan=None, logf=False, fmin=None, fmax=None, modes=None -): - raise NotImplementedError( - "createfourierdesignmatrix_eph was removed, " + "and replaced with createfourierdesignmatrix_ephem" - ) - - -@function -def createfourierdesignmatrix_chromatic(toas, freqs, nmodes=30, Tspan=None, logf=False, fmin=None, fmax=None, idx=4): ->>>>>>> 7e6931cbc5a819ee0bf223984eef09d963644927 """ Construct Scattering-variation fourier design matrix. @@ -324,17 +227,10 @@ def createfourierdesignmatrix_chromatic(toas, freqs, nmodes=30, Tspan=None, logf """ # get base fourier design matrix and frequencies -<<<<<<< HEAD F, Ffreqs = createfourierdesignmatrix_red(toas, nmodes=nmodes, Tspan=Tspan, logf=logf, fmin=fmin, fmax=fmax) # compute the DM-variation vectors Dm = (1400/freqs) ** idx -======= - F, Ffreqs = createfourierdesignmatrix_red(toas, nmodes=nmodes, Tspan=Tspan, logf=logf, fmin=fmin, fmax=fmax) - - # compute the DM-variation vectors - Dm = (1400 / freqs) ** idx ->>>>>>> 7e6931cbc5a819ee0bf223984eef09d963644927 return F * Dm[:, None], Ffreqs diff --git a/enterprise/signals/gp_priors.py b/enterprise/signals/gp_priors.py index ca3407a2..76d421d9 100644 --- a/enterprise/signals/gp_priors.py +++ b/enterprise/signals/gp_priors.py @@ -3,29 +3,20 @@ functions for use in other modules. """ -<<<<<<< HEAD from __future__ import (absolute_import, division, print_function, unicode_literals) -======= -from __future__ import absolute_import, division, print_function, unicode_literals ->>>>>>> 7e6931cbc5a819ee0bf223984eef09d963644927 import numpy as np import scipy.stats from enterprise.signals import parameter from enterprise.signals.parameter import function -<<<<<<< HEAD import enterprise.constants as const -======= -from enterprise import constants as const ->>>>>>> 7e6931cbc5a819ee0bf223984eef09d963644927 @function def powerlaw(f, log10_A=-16, gamma=5, components=2): df = np.diff(np.concatenate((np.array([0]), f[::components]))) -<<<<<<< HEAD return ((10**log10_A)**2 / 12.0 / np.pi**2 * const.fyr**(gamma-3) * f**(-gamma) * np.repeat(df, components)) @@ -36,18 +27,6 @@ def turnover(f, log10_A=-15, gamma=4.33, lf0=-8.5, kappa=10/3, beta=0.5): hcf = (10**log10_A * (f / const.fyr) ** ((3-gamma) / 2) / (1 + (10**lf0 / f) ** kappa) ** beta) return hcf**2/12/np.pi**2/f**3*np.repeat(df, 2) -======= - return ( - (10 ** log10_A) ** 2 / 12.0 / np.pi ** 2 * const.fyr ** (gamma - 3) * f ** (-gamma) * np.repeat(df, components) - ) - - -@function -def turnover(f, log10_A=-15, gamma=4.33, lf0=-8.5, kappa=10 / 3, beta=0.5): - df = np.diff(np.concatenate((np.array([0]), f[::2]))) - hcf = 10 ** log10_A * (f / const.fyr) ** ((3 - gamma) / 2) / (1 + (10 ** lf0 / f) ** kappa) ** beta - return hcf ** 2 / 12 / np.pi ** 2 / f ** 3 * np.repeat(df, 2) ->>>>>>> 7e6931cbc5a819ee0bf223984eef09d963644927 @function @@ -59,11 +38,7 @@ def free_spectrum(f, log10_rho=None): where \rho_i is the free parameter and T is the observation length. """ -<<<<<<< HEAD return np.repeat(10**(2*np.array(log10_rho)), 2) -======= - return np.repeat(10 ** (2 * np.array(log10_rho)), 2) ->>>>>>> 7e6931cbc5a819ee0bf223984eef09d963644927 @function @@ -89,13 +64,8 @@ def t_process_adapt(f, log10_A=-15, gamma=4.33, alphas_adapt=None, nfreq=None): alpha_model = np.repeat(alphas_adapt, 2) else: alpha_model = np.ones_like(f) -<<<<<<< HEAD alpha_model[2*int(np.rint(nfreq))] = alphas_adapt alpha_model[2*int(np.rint(nfreq))+1] = alphas_adapt -======= - alpha_model[2 * int(np.rint(nfreq))] = alphas_adapt - alpha_model[2 * int(np.rint(nfreq)) + 1] = alphas_adapt ->>>>>>> 7e6931cbc5a819ee0bf223984eef09d963644927 return powerlaw(f, log10_A=log10_A, gamma=gamma) * alpha_model @@ -112,10 +82,6 @@ def InvGammaSampler(alpha=1, gamma=1, size=None): def InvGamma(alpha=1, gamma=1, size=None): """Class factory for Inverse Gamma parameters.""" -<<<<<<< HEAD -======= - ->>>>>>> 7e6931cbc5a819ee0bf223984eef09d963644927 class InvGamma(parameter.Parameter): _size = size _prior = parameter.Function(InvGammaPrior, alpha=alpha, gamma=gamma) @@ -124,14 +90,8 @@ class InvGamma(parameter.Parameter): _gamma = gamma def __repr__(self): -<<<<<<< HEAD return '"{}": InvGamma({},{})'.format(self.name, alpha, gamma) \ + ('' if self._size is None else '[{}]'.format(self._size)) -======= - return '"{}": InvGamma({},{})'.format(self.name, alpha, gamma) + ( - "" if self._size is None else "[{}]".format(self._size) - ) ->>>>>>> 7e6931cbc5a819ee0bf223984eef09d963644927 return InvGamma @@ -149,20 +109,10 @@ def turnover_knee(f, log10_A, gamma, lfb, lfk, kappa, delta): :param delta: slope at higher frequencies """ df = np.diff(np.concatenate((np.array([0]), f[::2]))) -<<<<<<< HEAD hcf = (10**log10_A * (f / const.fyr) ** ((3-gamma) / 2) * (1.0 + (f / 10**lfk)) ** delta / np.sqrt(1 + (10**lfb / f) ** kappa)) return hcf**2 / 12 / np.pi**2 / f**3 * np.repeat(df, 2) -======= - hcf = ( - 10 ** log10_A - * (f / const.fyr) ** ((3 - gamma) / 2) - * (1.0 + (f / 10 ** lfk)) ** delta - / np.sqrt(1 + (10 ** lfb / f) ** kappa) - ) - return hcf ** 2 / 12 / np.pi ** 2 / f ** 3 * np.repeat(df, 2) ->>>>>>> 7e6931cbc5a819ee0bf223984eef09d963644927 @function @@ -179,35 +129,17 @@ def broken_powerlaw(f, log10_A, gamma, delta, log10_fb, kappa=0.1): :param kappa: smoothness of transition (Default = 0.1) """ df = np.diff(np.concatenate((np.array([0]), f[::2]))) -<<<<<<< HEAD hcf = (10**log10_A * (f / const.fyr) ** ((3-gamma) / 2) * (1 + (f / 10**log10_fb) ** (1/kappa)) ** (kappa * (gamma - delta) / 2)) return hcf**2 / 12 / np.pi**2 / f**3 * np.repeat(df, 2) -======= - hcf = ( - 10 ** log10_A - * (f / const.fyr) ** ((3 - gamma) / 2) - * (1 + (f / 10 ** log10_fb) ** (1 / kappa)) ** (kappa * (gamma - delta) / 2) - ) - return hcf ** 2 / 12 / np.pi ** 2 / f ** 3 * np.repeat(df, 2) ->>>>>>> 7e6931cbc5a819ee0bf223984eef09d963644927 @function def powerlaw_genmodes(f, log10_A=-16, gamma=5, components=2, wgts=None): if wgts is not None: -<<<<<<< HEAD df = wgts**2 else: df = np.diff(np.concatenate((np.array([0]), f[::components]))) return ((10**log10_A)**2 / 12.0 / np.pi**2 * const.fyr**(gamma-3) * f**(-gamma) * np.repeat(df, components)) -======= - df = wgts ** 2 - else: - df = np.diff(np.concatenate((np.array([0]), f[::components]))) - return ( - (10 ** log10_A) ** 2 / 12.0 / np.pi ** 2 * const.fyr ** (gamma - 3) * f ** (-gamma) * np.repeat(df, components) - ) ->>>>>>> 7e6931cbc5a819ee0bf223984eef09d963644927 diff --git a/tests/test_gp_priors.py b/tests/test_gp_priors.py index 81768d48..666745df 100644 --- a/tests/test_gp_priors.py +++ b/tests/test_gp_priors.py @@ -22,26 +22,18 @@ class TestGPSignals(unittest.TestCase): -<<<<<<< HEAD -======= ->>>>>>> 7e6931cbc5a819ee0bf223984eef09d963644927 @classmethod def setUpClass(cls): """Setup the Pulsar object.""" # initialize Pulsar class -<<<<<<< HEAD cls.psr = Pulsar(datadir + '/B1855+09_NANOGrav_9yv1.gls.par', datadir + '/B1855+09_NANOGrav_9yv1.tim') -======= - cls.psr = Pulsar(datadir + "/B1855+09_NANOGrav_9yv1.gls.par", datadir + "/B1855+09_NANOGrav_9yv1.tim") ->>>>>>> 7e6931cbc5a819ee0bf223984eef09d963644927 def test_turnover_prior(self): """Test that red noise signal returns correct values.""" # set up signal parameter -<<<<<<< HEAD pr = gp_priors.turnover(log10_A=parameter.Uniform(-18,-12), gamma=parameter.Uniform(1,7), lf0=parameter.Uniform(-9,-7.5), @@ -50,22 +42,10 @@ def test_turnover_prior(self): basis = gp_bases.createfourierdesignmatrix_red(nmodes=30) rn = gp_signals.BasisGP(priorFunction=pr, basisFunction=basis, name='red_noise') -======= - pr = gp_priors.turnover( - log10_A=parameter.Uniform(-18, -12), - gamma=parameter.Uniform(1, 7), - lf0=parameter.Uniform(-9, -7.5), - kappa=parameter.Uniform(2.5, 5), - beta=parameter.Uniform(0.01, 1), - ) - basis = gp_bases.createfourierdesignmatrix_red(nmodes=30) - rn = gp_signals.BasisGP(priorFunction=pr, basisFunction=basis, name="red_noise") ->>>>>>> 7e6931cbc5a819ee0bf223984eef09d963644927 rnm = rn(self.psr) # parameters log10_A, gamma, lf0, kappa, beta = -14.5, 4.33, -8.5, 3, 0.5 -<<<<<<< HEAD params = {'B1855+09_red_noise_log10_A': log10_A, 'B1855+09_red_noise_gamma': gamma, 'B1855+09_red_noise_lf0': lf0, @@ -94,38 +74,11 @@ def test_turnover_prior(self): # test shape msg = 'F matrix shape incorrect' -======= - params = { - "B1855+09_red_noise_log10_A": log10_A, - "B1855+09_red_noise_gamma": gamma, - "B1855+09_red_noise_lf0": lf0, - "B1855+09_red_noise_kappa": kappa, - "B1855+09_red_noise_beta": beta, - } - - # basis matrix test - F, f2 = gp_bases.createfourierdesignmatrix_red(self.psr.toas, nmodes=30) - msg = "F matrix incorrect for turnover." - assert np.allclose(F, rnm.get_basis(params)), msg - - # spectrum test - phi = gp_priors.turnover(f2, log10_A=log10_A, gamma=gamma, lf0=lf0, kappa=kappa, beta=beta) - msg = "Spectrum incorrect for turnover." - assert np.all(rnm.get_phi(params) == phi), msg - - # inverse spectrum test - msg = "Spectrum inverse incorrect for turnover." - assert np.all(rnm.get_phiinv(params) == 1 / phi), msg - - # test shape - msg = "F matrix shape incorrect" ->>>>>>> 7e6931cbc5a819ee0bf223984eef09d963644927 assert rnm.get_basis(params).shape == F.shape, msg def test_free_spec_prior(self): """Test that red noise signal returns correct values.""" # set up signal parameter -<<<<<<< HEAD pr = gp_priors.free_spectrum(log10_rho=parameter.Uniform(-10, -4, size=30)) basis = gp_bases.createfourierdesignmatrix_red(nmodes=30) @@ -139,24 +92,11 @@ def test_free_spec_prior(self): F, f2 = gp_bases.createfourierdesignmatrix_red(self.psr.toas, nmodes=30) msg = 'F matrix incorrect for free spectrum.' -======= - pr = gp_priors.free_spectrum(log10_rho=parameter.Uniform(-10, -4, size=30)) - basis = gp_bases.createfourierdesignmatrix_red(nmodes=30) - rn = gp_signals.BasisGP(priorFunction=pr, basisFunction=basis, name="red_noise") - rnm = rn(self.psr) - # parameters - rhos = np.random.uniform(-10, -4, size=30) - params = {"B1855+09_red_noise_log10_rho": rhos} - # basis matrix test - F, f2 = gp_bases.createfourierdesignmatrix_red(self.psr.toas, nmodes=30) - msg = "F matrix incorrect for free spectrum." ->>>>>>> 7e6931cbc5a819ee0bf223984eef09d963644927 assert np.allclose(F, rnm.get_basis(params)), msg # spectrum test phi = gp_priors.free_spectrum(f2, log10_rho=rhos) -<<<<<<< HEAD msg = 'Spectrum incorrect for free spectrum.' assert np.all(rnm.get_phi(params) == phi), msg @@ -166,23 +106,11 @@ def test_free_spec_prior(self): # test shape msg = 'F matrix shape incorrect' -======= - msg = "Spectrum incorrect for free spectrum." - assert np.all(rnm.get_phi(params) == phi), msg - - # inverse spectrum test - msg = "Spectrum inverse incorrect for free spectrum." - assert np.all(rnm.get_phiinv(params) == 1 / phi), msg - - # test shape - msg = "F matrix shape incorrect" ->>>>>>> 7e6931cbc5a819ee0bf223984eef09d963644927 assert rnm.get_basis(params).shape == F.shape, msg def test_t_process_prior(self): """Test that red noise signal returns correct values.""" # set up signal parameter -<<<<<<< HEAD pr = gp_priors.t_process(log10_A=parameter.Uniform(-18,-12), gamma=parameter.Uniform(1,7), alphas=gp_priors.InvGamma(alpha=1, gamma=1, @@ -190,20 +118,10 @@ def test_t_process_prior(self): basis = gp_bases.createfourierdesignmatrix_red(nmodes=30) rn = gp_signals.BasisGP(priorFunction=pr, basisFunction=basis, name='red_noise') -======= - pr = gp_priors.t_process( - log10_A=parameter.Uniform(-18, -12), - gamma=parameter.Uniform(1, 7), - alphas=gp_priors.InvGamma(alpha=1, gamma=1, size=30), - ) - basis = gp_bases.createfourierdesignmatrix_red(nmodes=30) - rn = gp_signals.BasisGP(priorFunction=pr, basisFunction=basis, name="red_noise") ->>>>>>> 7e6931cbc5a819ee0bf223984eef09d963644927 rnm = rn(self.psr) # parameters alphas = scipy.stats.invgamma.rvs(1, scale=1, size=30) log10_A, gamma = -15, 4.33 -<<<<<<< HEAD params = {'B1855+09_red_noise_log10_A': log10_A, 'B1855+09_red_noise_gamma': gamma, 'B1855+09_red_noise_alphas': alphas} @@ -226,36 +144,11 @@ def test_t_process_prior(self): # test shape msg = 'F matrix shape incorrect' -======= - params = { - "B1855+09_red_noise_log10_A": log10_A, - "B1855+09_red_noise_gamma": gamma, - "B1855+09_red_noise_alphas": alphas, - } - # basis matrix test - F, f2 = gp_bases.createfourierdesignmatrix_red(self.psr.toas, nmodes=30) - msg = "F matrix incorrect for free spectrum." - assert np.allclose(F, rnm.get_basis(params)), msg - - # spectrum test - phi = gp_priors.t_process(f2, log10_A=log10_A, gamma=gamma, alphas=alphas) - - msg = "Spectrum incorrect for free spectrum." - assert np.all(rnm.get_phi(params) == phi), msg - - # inverse spectrum test - msg = "Spectrum inverse incorrect for free spectrum." - assert np.all(rnm.get_phiinv(params) == 1 / phi), msg - - # test shape - msg = "F matrix shape incorrect" ->>>>>>> 7e6931cbc5a819ee0bf223984eef09d963644927 assert rnm.get_basis(params).shape == F.shape, msg def test_adapt_t_process_prior(self): """Test that red noise signal returns correct values.""" # set up signal parameter -<<<<<<< HEAD pr = gp_priors.t_process_adapt(log10_A=parameter.Uniform(-18,-12), gamma=parameter.Uniform(1,7), alphas_adapt=gp_priors.InvGamma(), @@ -263,21 +156,10 @@ def test_adapt_t_process_prior(self): basis = gp_bases.createfourierdesignmatrix_red(nmodes=30) rn = gp_signals.BasisGP(priorFunction=pr, basisFunction=basis, name='red_noise') -======= - pr = gp_priors.t_process_adapt( - log10_A=parameter.Uniform(-18, -12), - gamma=parameter.Uniform(1, 7), - alphas_adapt=gp_priors.InvGamma(), - nfreq=parameter.Uniform(5, 25), - ) - basis = gp_bases.createfourierdesignmatrix_red(nmodes=30) - rn = gp_signals.BasisGP(priorFunction=pr, basisFunction=basis, name="red_noise") ->>>>>>> 7e6931cbc5a819ee0bf223984eef09d963644927 rnm = rn(self.psr) # parameters alphas = scipy.stats.invgamma.rvs(1, scale=1, size=1) log10_A, gamma, nfreq = -15, 4.33, 12 -<<<<<<< HEAD params = {'B1855+09_red_noise_log10_A': log10_A, 'B1855+09_red_noise_gamma': gamma, 'B1855+09_red_noise_alphas_adapt': alphas, @@ -301,37 +183,11 @@ def test_adapt_t_process_prior(self): # test shape msg = 'F matrix shape incorrect' -======= - params = { - "B1855+09_red_noise_log10_A": log10_A, - "B1855+09_red_noise_gamma": gamma, - "B1855+09_red_noise_alphas_adapt": alphas, - "B1855+09_red_noise_nfreq": nfreq, - } - # basis matrix test - F, f2 = gp_bases.createfourierdesignmatrix_red(self.psr.toas, nmodes=30) - msg = "F matrix incorrect for free spectrum." - assert np.allclose(F, rnm.get_basis(params)), msg - - # spectrum test - phi = gp_priors.t_process_adapt(f2, log10_A=log10_A, gamma=gamma, alphas_adapt=alphas, nfreq=nfreq) - - msg = "Spectrum incorrect for free spectrum." - assert np.all(rnm.get_phi(params) == phi), msg - - # inverse spectrum test - msg = "Spectrum inverse incorrect for free spectrum." - assert np.all(rnm.get_phiinv(params) == 1 / phi), msg - - # test shape - msg = "F matrix shape incorrect" ->>>>>>> 7e6931cbc5a819ee0bf223984eef09d963644927 assert rnm.get_basis(params).shape == F.shape, msg def test_turnover_knee_prior(self): """Test that red noise signal returns correct values.""" # set up signal parameter -<<<<<<< HEAD pr = gp_priors.turnover_knee(log10_A=parameter.Uniform(-18,-12), gamma=parameter.Uniform(1,7), lfb=parameter.Uniform(-9,-7.5), @@ -341,24 +197,11 @@ def test_turnover_knee_prior(self): basis = gp_bases.createfourierdesignmatrix_red(nmodes=30) rn = gp_signals.BasisGP(priorFunction=pr, basisFunction=basis, name='red_noise') -======= - pr = gp_priors.turnover_knee( - log10_A=parameter.Uniform(-18, -12), - gamma=parameter.Uniform(1, 7), - lfb=parameter.Uniform(-9, -7.5), - lfk=parameter.Uniform(-9, -7.5), - kappa=parameter.Uniform(2.5, 5), - delta=parameter.Uniform(0.01, 1), - ) - basis = gp_bases.createfourierdesignmatrix_red(nmodes=30) - rn = gp_signals.BasisGP(priorFunction=pr, basisFunction=basis, name="red_noise") ->>>>>>> 7e6931cbc5a819ee0bf223984eef09d963644927 rnm = rn(self.psr) # parameters log10_A, gamma, lfb = -14.5, 4.33, -8.5 lfk, kappa, delta = -8.5, 3, 0.5 -<<<<<<< HEAD params = {'B1855+09_red_noise_log10_A': log10_A, 'B1855+09_red_noise_gamma': gamma, 'B1855+09_red_noise_lfb': lfb, @@ -389,39 +232,11 @@ def test_turnover_knee_prior(self): # test shape msg = 'F matrix shape incorrect' -======= - params = { - "B1855+09_red_noise_log10_A": log10_A, - "B1855+09_red_noise_gamma": gamma, - "B1855+09_red_noise_lfb": lfb, - "B1855+09_red_noise_lfk": lfk, - "B1855+09_red_noise_kappa": kappa, - "B1855+09_red_noise_delta": delta, - } - - # basis matrix test - F, f2 = gp_bases.createfourierdesignmatrix_red(self.psr.toas, nmodes=30) - msg = "F matrix incorrect for turnover." - assert np.allclose(F, rnm.get_basis(params)), msg - - # spectrum test - phi = gp_priors.turnover_knee(f2, log10_A=log10_A, gamma=gamma, lfb=lfb, lfk=lfk, kappa=kappa, delta=delta) - msg = "Spectrum incorrect for turnover." - assert np.all(rnm.get_phi(params) == phi), msg - - # inverse spectrum test - msg = "Spectrum inverse incorrect for turnover." - assert np.all(rnm.get_phiinv(params) == 1 / phi), msg - - # test shape - msg = "F matrix shape incorrect" ->>>>>>> 7e6931cbc5a819ee0bf223984eef09d963644927 assert rnm.get_basis(params).shape == F.shape, msg def test_broken_powerlaw_prior(self): """Test that red noise signal returns correct values.""" # set up signal parameter -<<<<<<< HEAD pr = gp_priors.broken_powerlaw(log10_A=parameter.Uniform(-18,-12), gamma=parameter.Uniform(1,7), log10_fb=parameter.Uniform(-9,-7.5), @@ -430,22 +245,10 @@ def test_broken_powerlaw_prior(self): basis = gp_bases.createfourierdesignmatrix_red(nmodes=30) rn = gp_signals.BasisGP(priorFunction=pr, basisFunction=basis, name='red_noise') -======= - pr = gp_priors.broken_powerlaw( - log10_A=parameter.Uniform(-18, -12), - gamma=parameter.Uniform(1, 7), - log10_fb=parameter.Uniform(-9, -7.5), - kappa=parameter.Uniform(0.1, 1.0), - delta=parameter.Uniform(0.01, 1), - ) - basis = gp_bases.createfourierdesignmatrix_red(nmodes=30) - rn = gp_signals.BasisGP(priorFunction=pr, basisFunction=basis, name="red_noise") ->>>>>>> 7e6931cbc5a819ee0bf223984eef09d963644927 rnm = rn(self.psr) # parameters log10_A, gamma, log10_fb, kappa, delta = -14.5, 4.33, -8.5, 1, 0.5 -<<<<<<< HEAD params = {'B1855+09_red_noise_log10_A': log10_A, 'B1855+09_red_noise_gamma': gamma, 'B1855+09_red_noise_log10_fb': log10_fb, @@ -474,53 +277,20 @@ def test_broken_powerlaw_prior(self): # test shape msg = 'F matrix shape incorrect' -======= - params = { - "B1855+09_red_noise_log10_A": log10_A, - "B1855+09_red_noise_gamma": gamma, - "B1855+09_red_noise_log10_fb": log10_fb, - "B1855+09_red_noise_kappa": kappa, - "B1855+09_red_noise_delta": delta, - } - - # basis matrix test - F, f2 = gp_bases.createfourierdesignmatrix_red(self.psr.toas, nmodes=30) - msg = "F matrix incorrect for turnover." - assert np.allclose(F, rnm.get_basis(params)), msg - - # spectrum test - phi = gp_priors.broken_powerlaw(f2, log10_A=log10_A, gamma=gamma, log10_fb=log10_fb, kappa=kappa, delta=delta) - msg = "Spectrum incorrect for turnover." - assert np.all(rnm.get_phi(params) == phi), msg - - # inverse spectrum test - msg = "Spectrum inverse incorrect for turnover." - assert np.all(rnm.get_phiinv(params) == 1 / phi), msg - - # test shape - msg = "F matrix shape incorrect" ->>>>>>> 7e6931cbc5a819ee0bf223984eef09d963644927 assert rnm.get_basis(params).shape == F.shape, msg def test_powerlaw_genmodes_prior(self): """Test that red noise signal returns correct values.""" # set up signal parameter -<<<<<<< HEAD pr = gp_priors.powerlaw_genmodes(log10_A=parameter.Uniform(-18,-12), gamma=parameter.Uniform(1,7)) basis = gp_bases.createfourierdesignmatrix_chromatic(nmodes=30) rn = gp_signals.BasisGP(priorFunction=pr, basisFunction=basis, name='red_noise') -======= - pr = gp_priors.powerlaw_genmodes(log10_A=parameter.Uniform(-18, -12), gamma=parameter.Uniform(1, 7)) - basis = gp_bases.createfourierdesignmatrix_chromatic(nmodes=30) - rn = gp_signals.BasisGP(priorFunction=pr, basisFunction=basis, name="red_noise") ->>>>>>> 7e6931cbc5a819ee0bf223984eef09d963644927 rnm = rn(self.psr) # parameters log10_A, gamma = -14.5, 4.33 -<<<<<<< HEAD params = {'B1855+09_red_noise_log10_A': log10_A, 'B1855+09_red_noise_gamma': gamma} @@ -544,24 +314,4 @@ def test_powerlaw_genmodes_prior(self): # test shape msg = 'F matrix shape incorrect' -======= - params = {"B1855+09_red_noise_log10_A": log10_A, "B1855+09_red_noise_gamma": gamma} - - # basis matrix test - F, f2 = gp_bases.createfourierdesignmatrix_chromatic(self.psr.toas, self.psr.freqs, nmodes=30) - msg = "F matrix incorrect for turnover." - assert np.allclose(F, rnm.get_basis(params)), msg - - # spectrum test - phi = gp_priors.powerlaw_genmodes(f2, log10_A=log10_A, gamma=gamma) - msg = "Spectrum incorrect for turnover." - assert np.all(rnm.get_phi(params) == phi), msg - - # inverse spectrum test - msg = "Spectrum inverse incorrect for turnover." - assert np.all(rnm.get_phiinv(params) == 1 / phi), msg - - # test shape - msg = "F matrix shape incorrect" ->>>>>>> 7e6931cbc5a819ee0bf223984eef09d963644927 assert rnm.get_basis(params).shape == F.shape, msg From c6c5dc42542708900cc673a0bf7988633c6b251a Mon Sep 17 00:00:00 2001 From: Hazboun6 Date: Thu, 27 Feb 2020 16:20:09 -0800 Subject: [PATCH 08/13] Removed test of planetssb and sunssb from PintPulsar --- tests/test_pulsar.py | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/tests/test_pulsar.py b/tests/test_pulsar.py index 4871f4ce..40e59bec 100644 --- a/tests/test_pulsar.py +++ b/tests/test_pulsar.py @@ -155,3 +155,9 @@ def test_stoas(self): def test_dm(self): pass + + def test_planetssb(self): + pass + + def test_sunssb(self): + pass From c788f1010a6bc5d3e7d13e73e8c43071e584f25c Mon Sep 17 00:00:00 2001 From: Hazboun6 Date: Thu, 27 Feb 2020 16:22:49 -0800 Subject: [PATCH 09/13] PintPulsar sets _get_sunssb to zero --- enterprise/pulsar.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/enterprise/pulsar.py b/enterprise/pulsar.py index de1165c5..fb00e352 100644 --- a/enterprise/pulsar.py +++ b/enterprise/pulsar.py @@ -344,6 +344,9 @@ def _get_radec(self, model): def _get_planetssb(self): return np.zeros((len(self._toas), 9, 6)) + def _get_sunssb(self): + return np.zeros((len(self._toas), 6)) + class Tempo2Pulsar(BasePulsar): def __init__(self, t2pulsar, sort=True, drop_t2pulsar=True, planets=True): From 76c3800e620879af23626eef6cacab76dfe6f16c Mon Sep 17 00:00:00 2001 From: Hazboun6 Date: Thu, 27 Feb 2020 16:57:21 -0800 Subject: [PATCH 10/13] Ran black to reformat --- enterprise/pulsar.py | 9 +- enterprise/signals/gp_bases.py | 97 ++++++------ enterprise/signals/gp_priors.py | 55 ++++--- tests/test_gp_priors.py | 273 +++++++++++++++----------------- tests/test_pulsar.py | 4 +- 5 files changed, 212 insertions(+), 226 deletions(-) diff --git a/enterprise/pulsar.py b/enterprise/pulsar.py index fb00e352..72e9b4cb 100644 --- a/enterprise/pulsar.py +++ b/enterprise/pulsar.py @@ -459,12 +459,9 @@ def _get_sunssb(self, t2pulsar): sunssb = np.zeros((len(self._toas), 6)) sunssb[:, :] = self.t2pulsar.sun_ssb - - if 'ELONG' and 'ELAT' in np.concatenate((t2pulsar.pars(), - t2pulsar.pars( - which='set'))): - sunssb[:,:3] = utils.ecl2eq_vec(sunssb[:,:3]) - sunssb[:,3:] = utils.ecl2eq_vec(sunssb[:,3:]) + if "ELONG" and "ELAT" in np.concatenate((t2pulsar.pars(), t2pulsar.pars(which="set"))): + sunssb[:, :3] = utils.ecl2eq_vec(sunssb[:, :3]) + sunssb[:, 3:] = utils.ecl2eq_vec(sunssb[:, 3:]) return sunssb diff --git a/enterprise/signals/gp_bases.py b/enterprise/signals/gp_bases.py index 63559123..031f5d80 100644 --- a/enterprise/signals/gp_bases.py +++ b/enterprise/signals/gp_bases.py @@ -3,26 +3,28 @@ functions for use in other modules. """ -from __future__ import (absolute_import, division, - print_function, unicode_literals) +from __future__ import absolute_import, division, print_function, unicode_literals import numpy as np from enterprise.signals.parameter import function + ###################################### # Fourier-basis signal functions ##### ###################################### -__all__ = ['createfourierdesignmatrix_red', - 'createfourierdesignmatrix_dm', - 'createfourierdesignmatrix_env', - 'createfourierdesignmatrix_ephem', - 'createfourierdesignmatrix_eph'] +__all__ = [ + "createfourierdesignmatrix_red", + "createfourierdesignmatrix_dm", + "createfourierdesignmatrix_env", + "createfourierdesignmatrix_ephem", + "createfourierdesignmatrix_eph", +] @function -def createfourierdesignmatrix_red(toas, nmodes=30, Tspan=None, - logf=False, fmin=None, fmax=None, - pshift=False, modes=None): +def createfourierdesignmatrix_red( + toas, nmodes=30, Tspan=None, logf=False, fmin=None, fmax=None, pshift=False, modes=None +): """ Construct fourier design matrix from eq 11 of Lentati et al, 2013 :param toas: vector of time series in seconds @@ -65,8 +67,7 @@ def createfourierdesignmatrix_red(toas, nmodes=30, Tspan=None, f = np.linspace(fmin, fmax, nmodes) # add random phase shift to basis functions - ranphase = (np.random.uniform(0.0, 2 * np.pi, nmodes) - if pshift else np.zeros(nmodes)) + ranphase = np.random.uniform(0.0, 2 * np.pi, nmodes) if pshift else np.zeros(nmodes) Ffreqs = np.repeat(f, 2) @@ -74,18 +75,16 @@ def createfourierdesignmatrix_red(toas, nmodes=30, Tspan=None, F = np.zeros((N, 2 * nmodes)) # The sine/cosine modes - F[:, ::2] = np.sin(2*np.pi*toas[:, None]*f[None, :] + - ranphase[None, :]) - F[:, 1::2] = np.cos(2*np.pi*toas[:, None]*f[None, :] + - ranphase[None, :]) + F[:, ::2] = np.sin(2 * np.pi * toas[:, None] * f[None, :] + ranphase[None, :]) + F[:, 1::2] = np.cos(2 * np.pi * toas[:, None] * f[None, :] + ranphase[None, :]) return F, Ffreqs @function -def createfourierdesignmatrix_dm(toas, freqs, nmodes=30, Tspan=None, - pshift=False, fref=1400, logf=False, - fmin=None, fmax=None, modes=None): +def createfourierdesignmatrix_dm( + toas, freqs, nmodes=30, Tspan=None, pshift=False, fref=1400, logf=False, fmin=None, fmax=None, modes=None +): """ Construct DM-variation fourier design matrix. Current normalization expresses DM signal as a deviation [seconds] @@ -109,20 +108,28 @@ def createfourierdesignmatrix_dm(toas, freqs, nmodes=30, Tspan=None, # get base fourier design matrix and frequencies F, Ffreqs = createfourierdesignmatrix_red( - toas, nmodes=nmodes, Tspan=Tspan, logf=logf, - fmin=fmin, fmax=fmax, pshift=pshift, modes=modes) + toas, nmodes=nmodes, Tspan=Tspan, logf=logf, fmin=fmin, fmax=fmax, pshift=pshift, modes=modes + ) # compute the DM-variation vectors - Dm = (fref/freqs)**2 + Dm = (fref / freqs) ** 2 return F * Dm[:, None], Ffreqs @function -def createfourierdesignmatrix_env(toas, log10_Amp=-7, log10_Q=np.log10(300), - t0=53000*86400, nmodes=30, Tspan=None, - logf=False, fmin=None, fmax=None, - modes=None): +def createfourierdesignmatrix_env( + toas, + log10_Amp=-7, + log10_Q=np.log10(300), + t0=53000 * 86400, + nmodes=30, + Tspan=None, + logf=False, + fmin=None, + fmax=None, + modes=None, +): """ Construct fourier design matrix with gaussian envelope. @@ -146,13 +153,13 @@ def createfourierdesignmatrix_env(toas, log10_Amp=-7, log10_Q=np.log10(300), # get base fourier design matrix and frequencies F, Ffreqs = createfourierdesignmatrix_red( - toas, nmodes=nmodes, Tspan=Tspan, logf=logf, - fmin=fmin, fmax=fmax, modes=modes) + toas, nmodes=nmodes, Tspan=Tspan, logf=logf, fmin=fmin, fmax=fmax, modes=modes + ) # compute gaussian envelope - A = 10**log10_Amp - Q = 10**log10_Q * 86400 - env = A * np.exp(-(toas-t0)**2/2/Q**2) + A = 10 ** log10_Amp + Q = 10 ** log10_Q * 86400 + env = A * np.exp(-((toas - t0) ** 2) / 2 / Q ** 2) return F * env[:, None], Ffreqs @@ -180,34 +187,31 @@ def createfourierdesignmatrix_ephem(toas, pos, nmodes=30, Tspan=None): :return: f: Sampling frequencies (6*nmodes) """ - F0, F0f = createfourierdesignmatrix_red( - toas, nmodes=nmodes, Tspan=Tspan) + F0, F0f = createfourierdesignmatrix_red(toas, nmodes=nmodes, Tspan=Tspan) - F1 = np.zeros((len(toas), nmodes, 2, 3), 'd') + F1 = np.zeros((len(toas), nmodes, 2, 3), "d") F1[:, :, 0, :] = F0[:, 0::2, np.newaxis] F1[:, :, 1, :] = F0[:, 1::2, np.newaxis] # verify this is the scalar product we want F1 *= pos - F1f = np.zeros((nmodes, 2, 3), 'd') + F1f = np.zeros((nmodes, 2, 3), "d") F1f[:, :, :] = F0f[::2, np.newaxis, np.newaxis] - return F1.reshape((len(toas), nmodes*6)), F1f.reshape((nmodes*6, )) + return F1.reshape((len(toas), nmodes * 6)), F1f.reshape((nmodes * 6,)) -def createfourierdesignmatrix_eph(t, nmodes, phi, theta, freq=False, - Tspan=None, logf=False, fmin=None, - fmax=None, modes=None): +def createfourierdesignmatrix_eph( + t, nmodes, phi, theta, freq=False, Tspan=None, logf=False, fmin=None, fmax=None, modes=None +): raise NotImplementedError( - "createfourierdesignmatrix_eph was removed, " + - "and replaced with createfourierdesignmatrix_ephem") + "createfourierdesignmatrix_eph was removed, " + "and replaced with createfourierdesignmatrix_ephem" + ) @function -def createfourierdesignmatrix_chromatic(toas, freqs, nmodes=30, Tspan=None, - logf=False, fmin=None, fmax=None, - idx=4): +def createfourierdesignmatrix_chromatic(toas, freqs, nmodes=30, Tspan=None, logf=False, fmin=None, fmax=None, idx=4): """ Construct Scattering-variation fourier design matrix. @@ -227,10 +231,9 @@ def createfourierdesignmatrix_chromatic(toas, freqs, nmodes=30, Tspan=None, """ # get base fourier design matrix and frequencies - F, Ffreqs = createfourierdesignmatrix_red(toas, nmodes=nmodes, Tspan=Tspan, - logf=logf, fmin=fmin, fmax=fmax) + F, Ffreqs = createfourierdesignmatrix_red(toas, nmodes=nmodes, Tspan=Tspan, logf=logf, fmin=fmin, fmax=fmax) # compute the DM-variation vectors - Dm = (1400/freqs) ** idx + Dm = (1400 / freqs) ** idx return F * Dm[:, None], Ffreqs diff --git a/enterprise/signals/gp_priors.py b/enterprise/signals/gp_priors.py index 76d421d9..89417f45 100644 --- a/enterprise/signals/gp_priors.py +++ b/enterprise/signals/gp_priors.py @@ -3,8 +3,7 @@ functions for use in other modules. """ -from __future__ import (absolute_import, division, - print_function, unicode_literals) +from __future__ import absolute_import, division, print_function, unicode_literals import numpy as np import scipy.stats @@ -17,16 +16,16 @@ @function def powerlaw(f, log10_A=-16, gamma=5, components=2): df = np.diff(np.concatenate((np.array([0]), f[::components]))) - return ((10**log10_A)**2 / 12.0 / np.pi**2 * - const.fyr**(gamma-3) * f**(-gamma) * np.repeat(df, components)) + return ( + (10 ** log10_A) ** 2 / 12.0 / np.pi ** 2 * const.fyr ** (gamma - 3) * f ** (-gamma) * np.repeat(df, components) + ) @function -def turnover(f, log10_A=-15, gamma=4.33, lf0=-8.5, kappa=10/3, beta=0.5): +def turnover(f, log10_A=-15, gamma=4.33, lf0=-8.5, kappa=10 / 3, beta=0.5): df = np.diff(np.concatenate((np.array([0]), f[::2]))) - hcf = (10**log10_A * (f / const.fyr) ** ((3-gamma) / 2) / - (1 + (10**lf0 / f) ** kappa) ** beta) - return hcf**2/12/np.pi**2/f**3*np.repeat(df, 2) + hcf = 10 ** log10_A * (f / const.fyr) ** ((3 - gamma) / 2) / (1 + (10 ** lf0 / f) ** kappa) ** beta + return hcf ** 2 / 12 / np.pi ** 2 / f ** 3 * np.repeat(df, 2) @function @@ -38,7 +37,7 @@ def free_spectrum(f, log10_rho=None): where \rho_i is the free parameter and T is the observation length. """ - return np.repeat(10**(2*np.array(log10_rho)), 2) + return np.repeat(10 ** (2 * np.array(log10_rho)), 2) @function @@ -64,8 +63,8 @@ def t_process_adapt(f, log10_A=-15, gamma=4.33, alphas_adapt=None, nfreq=None): alpha_model = np.repeat(alphas_adapt, 2) else: alpha_model = np.ones_like(f) - alpha_model[2*int(np.rint(nfreq))] = alphas_adapt - alpha_model[2*int(np.rint(nfreq))+1] = alphas_adapt + alpha_model[2 * int(np.rint(nfreq))] = alphas_adapt + alpha_model[2 * int(np.rint(nfreq)) + 1] = alphas_adapt return powerlaw(f, log10_A=log10_A, gamma=gamma) * alpha_model @@ -82,6 +81,7 @@ def InvGammaSampler(alpha=1, gamma=1, size=None): def InvGamma(alpha=1, gamma=1, size=None): """Class factory for Inverse Gamma parameters.""" + class InvGamma(parameter.Parameter): _size = size _prior = parameter.Function(InvGammaPrior, alpha=alpha, gamma=gamma) @@ -90,8 +90,9 @@ class InvGamma(parameter.Parameter): _gamma = gamma def __repr__(self): - return '"{}": InvGamma({},{})'.format(self.name, alpha, gamma) \ - + ('' if self._size is None else '[{}]'.format(self._size)) + return '"{}": InvGamma({},{})'.format(self.name, alpha, gamma) + ( + "" if self._size is None else "[{}]".format(self._size) + ) return InvGamma @@ -109,10 +110,13 @@ def turnover_knee(f, log10_A, gamma, lfb, lfk, kappa, delta): :param delta: slope at higher frequencies """ df = np.diff(np.concatenate((np.array([0]), f[::2]))) - hcf = (10**log10_A * (f / const.fyr) ** ((3-gamma) / 2) * - (1.0 + (f / 10**lfk)) ** delta / - np.sqrt(1 + (10**lfb / f) ** kappa)) - return hcf**2 / 12 / np.pi**2 / f**3 * np.repeat(df, 2) + hcf = ( + 10 ** log10_A + * (f / const.fyr) ** ((3 - gamma) / 2) + * (1.0 + (f / 10 ** lfk)) ** delta + / np.sqrt(1 + (10 ** lfb / f) ** kappa) + ) + return hcf ** 2 / 12 / np.pi ** 2 / f ** 3 * np.repeat(df, 2) @function @@ -129,17 +133,20 @@ def broken_powerlaw(f, log10_A, gamma, delta, log10_fb, kappa=0.1): :param kappa: smoothness of transition (Default = 0.1) """ df = np.diff(np.concatenate((np.array([0]), f[::2]))) - hcf = (10**log10_A * (f / const.fyr) ** ((3-gamma) / 2) * - (1 + (f / 10**log10_fb) ** (1/kappa)) ** - (kappa * (gamma - delta) / 2)) - return hcf**2 / 12 / np.pi**2 / f**3 * np.repeat(df, 2) + hcf = ( + 10 ** log10_A + * (f / const.fyr) ** ((3 - gamma) / 2) + * (1 + (f / 10 ** log10_fb) ** (1 / kappa)) ** (kappa * (gamma - delta) / 2) + ) + return hcf ** 2 / 12 / np.pi ** 2 / f ** 3 * np.repeat(df, 2) @function def powerlaw_genmodes(f, log10_A=-16, gamma=5, components=2, wgts=None): if wgts is not None: - df = wgts**2 + df = wgts ** 2 else: df = np.diff(np.concatenate((np.array([0]), f[::components]))) - return ((10**log10_A)**2 / 12.0 / np.pi**2 * - const.fyr**(gamma-3) * f**(-gamma) * np.repeat(df, components)) + return ( + (10 ** log10_A) ** 2 / 12.0 / np.pi ** 2 * const.fyr ** (gamma - 3) * f ** (-gamma) * np.repeat(df, components) + ) diff --git a/tests/test_gp_priors.py b/tests/test_gp_priors.py index 666745df..b60b8c9b 100644 --- a/tests/test_gp_priors.py +++ b/tests/test_gp_priors.py @@ -22,296 +22,275 @@ class TestGPSignals(unittest.TestCase): - @classmethod def setUpClass(cls): """Setup the Pulsar object.""" # initialize Pulsar class - cls.psr = Pulsar(datadir + '/B1855+09_NANOGrav_9yv1.gls.par', - datadir + '/B1855+09_NANOGrav_9yv1.tim') + cls.psr = Pulsar(datadir + "/B1855+09_NANOGrav_9yv1.gls.par", datadir + "/B1855+09_NANOGrav_9yv1.tim") def test_turnover_prior(self): """Test that red noise signal returns correct values.""" # set up signal parameter - pr = gp_priors.turnover(log10_A=parameter.Uniform(-18,-12), - gamma=parameter.Uniform(1,7), - lf0=parameter.Uniform(-9,-7.5), - kappa=parameter.Uniform(2.5,5), - beta=parameter.Uniform(0.01,1),) + pr = gp_priors.turnover( + log10_A=parameter.Uniform(-18, -12), + gamma=parameter.Uniform(1, 7), + lf0=parameter.Uniform(-9, -7.5), + kappa=parameter.Uniform(2.5, 5), + beta=parameter.Uniform(0.01, 1), + ) basis = gp_bases.createfourierdesignmatrix_red(nmodes=30) - rn = gp_signals.BasisGP(priorFunction=pr, basisFunction=basis, - name='red_noise') + rn = gp_signals.BasisGP(priorFunction=pr, basisFunction=basis, name="red_noise") rnm = rn(self.psr) # parameters log10_A, gamma, lf0, kappa, beta = -14.5, 4.33, -8.5, 3, 0.5 - params = {'B1855+09_red_noise_log10_A': log10_A, - 'B1855+09_red_noise_gamma': gamma, - 'B1855+09_red_noise_lf0': lf0, - 'B1855+09_red_noise_kappa': kappa, - 'B1855+09_red_noise_beta': beta} + params = { + "B1855+09_red_noise_log10_A": log10_A, + "B1855+09_red_noise_gamma": gamma, + "B1855+09_red_noise_lf0": lf0, + "B1855+09_red_noise_kappa": kappa, + "B1855+09_red_noise_beta": beta, + } # basis matrix test - F, f2 = gp_bases.createfourierdesignmatrix_red(self.psr.toas, - nmodes=30) - msg = 'F matrix incorrect for turnover.' + F, f2 = gp_bases.createfourierdesignmatrix_red(self.psr.toas, nmodes=30) + msg = "F matrix incorrect for turnover." assert np.allclose(F, rnm.get_basis(params)), msg # spectrum test - phi = gp_priors.turnover(f2, - log10_A=log10_A, - gamma=gamma, - lf0=lf0, - kappa=kappa, - beta=beta) - msg = 'Spectrum incorrect for turnover.' + phi = gp_priors.turnover(f2, log10_A=log10_A, gamma=gamma, lf0=lf0, kappa=kappa, beta=beta) + msg = "Spectrum incorrect for turnover." assert np.all(rnm.get_phi(params) == phi), msg # inverse spectrum test - msg = 'Spectrum inverse incorrect for turnover.' - assert np.all(rnm.get_phiinv(params) == 1/phi), msg + msg = "Spectrum inverse incorrect for turnover." + assert np.all(rnm.get_phiinv(params) == 1 / phi), msg # test shape - msg = 'F matrix shape incorrect' + msg = "F matrix shape incorrect" assert rnm.get_basis(params).shape == F.shape, msg def test_free_spec_prior(self): """Test that red noise signal returns correct values.""" # set up signal parameter - pr = gp_priors.free_spectrum(log10_rho=parameter.Uniform(-10, -4, - size=30)) + pr = gp_priors.free_spectrum(log10_rho=parameter.Uniform(-10, -4, size=30)) basis = gp_bases.createfourierdesignmatrix_red(nmodes=30) - rn = gp_signals.BasisGP(priorFunction=pr, basisFunction=basis, - name='red_noise') + rn = gp_signals.BasisGP(priorFunction=pr, basisFunction=basis, name="red_noise") rnm = rn(self.psr) # parameters rhos = np.random.uniform(-10, -4, size=30) - params = {'B1855+09_red_noise_log10_rho': rhos} + params = {"B1855+09_red_noise_log10_rho": rhos} # basis matrix test - F, f2 = gp_bases.createfourierdesignmatrix_red(self.psr.toas, - nmodes=30) - msg = 'F matrix incorrect for free spectrum.' + F, f2 = gp_bases.createfourierdesignmatrix_red(self.psr.toas, nmodes=30) + msg = "F matrix incorrect for free spectrum." assert np.allclose(F, rnm.get_basis(params)), msg # spectrum test phi = gp_priors.free_spectrum(f2, log10_rho=rhos) - msg = 'Spectrum incorrect for free spectrum.' + msg = "Spectrum incorrect for free spectrum." assert np.all(rnm.get_phi(params) == phi), msg # inverse spectrum test - msg = 'Spectrum inverse incorrect for free spectrum.' - assert np.all(rnm.get_phiinv(params) == 1/phi), msg + msg = "Spectrum inverse incorrect for free spectrum." + assert np.all(rnm.get_phiinv(params) == 1 / phi), msg # test shape - msg = 'F matrix shape incorrect' + msg = "F matrix shape incorrect" assert rnm.get_basis(params).shape == F.shape, msg def test_t_process_prior(self): """Test that red noise signal returns correct values.""" # set up signal parameter - pr = gp_priors.t_process(log10_A=parameter.Uniform(-18,-12), - gamma=parameter.Uniform(1,7), - alphas=gp_priors.InvGamma(alpha=1, gamma=1, - size=30)) + pr = gp_priors.t_process( + log10_A=parameter.Uniform(-18, -12), + gamma=parameter.Uniform(1, 7), + alphas=gp_priors.InvGamma(alpha=1, gamma=1, size=30), + ) basis = gp_bases.createfourierdesignmatrix_red(nmodes=30) - rn = gp_signals.BasisGP(priorFunction=pr, basisFunction=basis, - name='red_noise') + rn = gp_signals.BasisGP(priorFunction=pr, basisFunction=basis, name="red_noise") rnm = rn(self.psr) # parameters alphas = scipy.stats.invgamma.rvs(1, scale=1, size=30) log10_A, gamma = -15, 4.33 - params = {'B1855+09_red_noise_log10_A': log10_A, - 'B1855+09_red_noise_gamma': gamma, - 'B1855+09_red_noise_alphas': alphas} + params = { + "B1855+09_red_noise_log10_A": log10_A, + "B1855+09_red_noise_gamma": gamma, + "B1855+09_red_noise_alphas": alphas, + } # basis matrix test - F, f2 = gp_bases.createfourierdesignmatrix_red(self.psr.toas, - nmodes=30) - msg = 'F matrix incorrect for free spectrum.' + F, f2 = gp_bases.createfourierdesignmatrix_red(self.psr.toas, nmodes=30) + msg = "F matrix incorrect for free spectrum." assert np.allclose(F, rnm.get_basis(params)), msg # spectrum test - phi = gp_priors.t_process(f2, log10_A=log10_A, gamma=gamma, - alphas=alphas) + phi = gp_priors.t_process(f2, log10_A=log10_A, gamma=gamma, alphas=alphas) - msg = 'Spectrum incorrect for free spectrum.' + msg = "Spectrum incorrect for free spectrum." assert np.all(rnm.get_phi(params) == phi), msg # inverse spectrum test - msg = 'Spectrum inverse incorrect for free spectrum.' - assert np.all(rnm.get_phiinv(params) == 1/phi), msg + msg = "Spectrum inverse incorrect for free spectrum." + assert np.all(rnm.get_phiinv(params) == 1 / phi), msg # test shape - msg = 'F matrix shape incorrect' + msg = "F matrix shape incorrect" assert rnm.get_basis(params).shape == F.shape, msg def test_adapt_t_process_prior(self): """Test that red noise signal returns correct values.""" # set up signal parameter - pr = gp_priors.t_process_adapt(log10_A=parameter.Uniform(-18,-12), - gamma=parameter.Uniform(1,7), - alphas_adapt=gp_priors.InvGamma(), - nfreq=parameter.Uniform(5,25)) + pr = gp_priors.t_process_adapt( + log10_A=parameter.Uniform(-18, -12), + gamma=parameter.Uniform(1, 7), + alphas_adapt=gp_priors.InvGamma(), + nfreq=parameter.Uniform(5, 25), + ) basis = gp_bases.createfourierdesignmatrix_red(nmodes=30) - rn = gp_signals.BasisGP(priorFunction=pr, basisFunction=basis, - name='red_noise') + rn = gp_signals.BasisGP(priorFunction=pr, basisFunction=basis, name="red_noise") rnm = rn(self.psr) # parameters alphas = scipy.stats.invgamma.rvs(1, scale=1, size=1) log10_A, gamma, nfreq = -15, 4.33, 12 - params = {'B1855+09_red_noise_log10_A': log10_A, - 'B1855+09_red_noise_gamma': gamma, - 'B1855+09_red_noise_alphas_adapt': alphas, - 'B1855+09_red_noise_nfreq': nfreq} + params = { + "B1855+09_red_noise_log10_A": log10_A, + "B1855+09_red_noise_gamma": gamma, + "B1855+09_red_noise_alphas_adapt": alphas, + "B1855+09_red_noise_nfreq": nfreq, + } # basis matrix test - F, f2 = gp_bases.createfourierdesignmatrix_red(self.psr.toas, - nmodes=30) - msg = 'F matrix incorrect for free spectrum.' + F, f2 = gp_bases.createfourierdesignmatrix_red(self.psr.toas, nmodes=30) + msg = "F matrix incorrect for free spectrum." assert np.allclose(F, rnm.get_basis(params)), msg # spectrum test - phi = gp_priors.t_process_adapt(f2, log10_A=log10_A, gamma=gamma, - alphas_adapt=alphas, nfreq=nfreq) + phi = gp_priors.t_process_adapt(f2, log10_A=log10_A, gamma=gamma, alphas_adapt=alphas, nfreq=nfreq) - msg = 'Spectrum incorrect for free spectrum.' + msg = "Spectrum incorrect for free spectrum." assert np.all(rnm.get_phi(params) == phi), msg # inverse spectrum test - msg = 'Spectrum inverse incorrect for free spectrum.' - assert np.all(rnm.get_phiinv(params) == 1/phi), msg + msg = "Spectrum inverse incorrect for free spectrum." + assert np.all(rnm.get_phiinv(params) == 1 / phi), msg # test shape - msg = 'F matrix shape incorrect' + msg = "F matrix shape incorrect" assert rnm.get_basis(params).shape == F.shape, msg def test_turnover_knee_prior(self): """Test that red noise signal returns correct values.""" # set up signal parameter - pr = gp_priors.turnover_knee(log10_A=parameter.Uniform(-18,-12), - gamma=parameter.Uniform(1,7), - lfb=parameter.Uniform(-9,-7.5), - lfk=parameter.Uniform(-9,-7.5), - kappa=parameter.Uniform(2.5,5), - delta=parameter.Uniform(0.01,1),) + pr = gp_priors.turnover_knee( + log10_A=parameter.Uniform(-18, -12), + gamma=parameter.Uniform(1, 7), + lfb=parameter.Uniform(-9, -7.5), + lfk=parameter.Uniform(-9, -7.5), + kappa=parameter.Uniform(2.5, 5), + delta=parameter.Uniform(0.01, 1), + ) basis = gp_bases.createfourierdesignmatrix_red(nmodes=30) - rn = gp_signals.BasisGP(priorFunction=pr, basisFunction=basis, - name='red_noise') + rn = gp_signals.BasisGP(priorFunction=pr, basisFunction=basis, name="red_noise") rnm = rn(self.psr) # parameters log10_A, gamma, lfb = -14.5, 4.33, -8.5 lfk, kappa, delta = -8.5, 3, 0.5 - params = {'B1855+09_red_noise_log10_A': log10_A, - 'B1855+09_red_noise_gamma': gamma, - 'B1855+09_red_noise_lfb': lfb, - 'B1855+09_red_noise_lfk': lfk, - 'B1855+09_red_noise_kappa': kappa, - 'B1855+09_red_noise_delta': delta} + params = { + "B1855+09_red_noise_log10_A": log10_A, + "B1855+09_red_noise_gamma": gamma, + "B1855+09_red_noise_lfb": lfb, + "B1855+09_red_noise_lfk": lfk, + "B1855+09_red_noise_kappa": kappa, + "B1855+09_red_noise_delta": delta, + } # basis matrix test - F, f2 = gp_bases.createfourierdesignmatrix_red(self.psr.toas, - nmodes=30) - msg = 'F matrix incorrect for turnover.' + F, f2 = gp_bases.createfourierdesignmatrix_red(self.psr.toas, nmodes=30) + msg = "F matrix incorrect for turnover." assert np.allclose(F, rnm.get_basis(params)), msg # spectrum test - phi = gp_priors.turnover_knee(f2, - log10_A=log10_A, - gamma=gamma, - lfb=lfb, - lfk=lfk, - kappa=kappa, - delta=delta) - msg = 'Spectrum incorrect for turnover.' + phi = gp_priors.turnover_knee(f2, log10_A=log10_A, gamma=gamma, lfb=lfb, lfk=lfk, kappa=kappa, delta=delta) + msg = "Spectrum incorrect for turnover." assert np.all(rnm.get_phi(params) == phi), msg # inverse spectrum test - msg = 'Spectrum inverse incorrect for turnover.' - assert np.all(rnm.get_phiinv(params) == 1/phi), msg + msg = "Spectrum inverse incorrect for turnover." + assert np.all(rnm.get_phiinv(params) == 1 / phi), msg # test shape - msg = 'F matrix shape incorrect' + msg = "F matrix shape incorrect" assert rnm.get_basis(params).shape == F.shape, msg def test_broken_powerlaw_prior(self): """Test that red noise signal returns correct values.""" # set up signal parameter - pr = gp_priors.broken_powerlaw(log10_A=parameter.Uniform(-18,-12), - gamma=parameter.Uniform(1,7), - log10_fb=parameter.Uniform(-9,-7.5), - kappa=parameter.Uniform(0.1,1.0), - delta=parameter.Uniform(0.01,1),) + pr = gp_priors.broken_powerlaw( + log10_A=parameter.Uniform(-18, -12), + gamma=parameter.Uniform(1, 7), + log10_fb=parameter.Uniform(-9, -7.5), + kappa=parameter.Uniform(0.1, 1.0), + delta=parameter.Uniform(0.01, 1), + ) basis = gp_bases.createfourierdesignmatrix_red(nmodes=30) - rn = gp_signals.BasisGP(priorFunction=pr, basisFunction=basis, - name='red_noise') + rn = gp_signals.BasisGP(priorFunction=pr, basisFunction=basis, name="red_noise") rnm = rn(self.psr) # parameters log10_A, gamma, log10_fb, kappa, delta = -14.5, 4.33, -8.5, 1, 0.5 - params = {'B1855+09_red_noise_log10_A': log10_A, - 'B1855+09_red_noise_gamma': gamma, - 'B1855+09_red_noise_log10_fb': log10_fb, - 'B1855+09_red_noise_kappa': kappa, - 'B1855+09_red_noise_delta': delta} + params = { + "B1855+09_red_noise_log10_A": log10_A, + "B1855+09_red_noise_gamma": gamma, + "B1855+09_red_noise_log10_fb": log10_fb, + "B1855+09_red_noise_kappa": kappa, + "B1855+09_red_noise_delta": delta, + } # basis matrix test - F, f2 = gp_bases.createfourierdesignmatrix_red(self.psr.toas, - nmodes=30) - msg = 'F matrix incorrect for turnover.' + F, f2 = gp_bases.createfourierdesignmatrix_red(self.psr.toas, nmodes=30) + msg = "F matrix incorrect for turnover." assert np.allclose(F, rnm.get_basis(params)), msg # spectrum test - phi = gp_priors.broken_powerlaw(f2, - log10_A=log10_A, - gamma=gamma, - log10_fb=log10_fb, - kappa=kappa, - delta=delta) - msg = 'Spectrum incorrect for turnover.' + phi = gp_priors.broken_powerlaw(f2, log10_A=log10_A, gamma=gamma, log10_fb=log10_fb, kappa=kappa, delta=delta) + msg = "Spectrum incorrect for turnover." assert np.all(rnm.get_phi(params) == phi), msg # inverse spectrum test - msg = 'Spectrum inverse incorrect for turnover.' - assert np.all(rnm.get_phiinv(params) == 1/phi), msg + msg = "Spectrum inverse incorrect for turnover." + assert np.all(rnm.get_phiinv(params) == 1 / phi), msg # test shape - msg = 'F matrix shape incorrect' + msg = "F matrix shape incorrect" assert rnm.get_basis(params).shape == F.shape, msg def test_powerlaw_genmodes_prior(self): """Test that red noise signal returns correct values.""" # set up signal parameter - pr = gp_priors.powerlaw_genmodes(log10_A=parameter.Uniform(-18,-12), - gamma=parameter.Uniform(1,7)) + pr = gp_priors.powerlaw_genmodes(log10_A=parameter.Uniform(-18, -12), gamma=parameter.Uniform(1, 7)) basis = gp_bases.createfourierdesignmatrix_chromatic(nmodes=30) - rn = gp_signals.BasisGP(priorFunction=pr, basisFunction=basis, - name='red_noise') + rn = gp_signals.BasisGP(priorFunction=pr, basisFunction=basis, name="red_noise") rnm = rn(self.psr) # parameters log10_A, gamma = -14.5, 4.33 - params = {'B1855+09_red_noise_log10_A': log10_A, - 'B1855+09_red_noise_gamma': gamma} + params = {"B1855+09_red_noise_log10_A": log10_A, "B1855+09_red_noise_gamma": gamma} # basis matrix test - F, f2 = gp_bases.createfourierdesignmatrix_chromatic(self.psr.toas, - self.psr.freqs, - nmodes=30) - msg = 'F matrix incorrect for turnover.' + F, f2 = gp_bases.createfourierdesignmatrix_chromatic(self.psr.toas, self.psr.freqs, nmodes=30) + msg = "F matrix incorrect for turnover." assert np.allclose(F, rnm.get_basis(params)), msg # spectrum test - phi = gp_priors.powerlaw_genmodes(f2, - log10_A=log10_A, - gamma=gamma) - msg = 'Spectrum incorrect for turnover.' + phi = gp_priors.powerlaw_genmodes(f2, log10_A=log10_A, gamma=gamma) + msg = "Spectrum incorrect for turnover." assert np.all(rnm.get_phi(params) == phi), msg # inverse spectrum test - msg = 'Spectrum inverse incorrect for turnover.' - assert np.all(rnm.get_phiinv(params) == 1/phi), msg + msg = "Spectrum inverse incorrect for turnover." + assert np.all(rnm.get_phiinv(params) == 1 / phi), msg # test shape - msg = 'F matrix shape incorrect' + msg = "F matrix shape incorrect" assert rnm.get_basis(params).shape == F.shape, msg diff --git a/tests/test_pulsar.py b/tests/test_pulsar.py index 40e59bec..d4128d90 100644 --- a/tests/test_pulsar.py +++ b/tests/test_pulsar.py @@ -105,11 +105,11 @@ def test_filter_data(self): def test_planetssb(self): """Place holder for filter_data tests.""" - assert hasattr(self.psr, 'planetssb') + assert hasattr(self.psr, "planetssb") def test_sunssb(self): """Place holder for filter_data tests.""" - assert hasattr(self.psr, 'sunssb') + assert hasattr(self.psr, "sunssb") def test_to_pickle(self): """Place holder for to_pickle tests.""" From f59783ac6c99bf65de6ba61b599a84eace3ddf77 Mon Sep 17 00:00:00 2001 From: Hazboun6 Date: Mon, 18 Jan 2021 11:21:48 -0800 Subject: [PATCH 11/13] start adding SSB vec to PintPulsar --- enterprise/pulsar.py | 41 ++++++++++++++++++++++++++++++++++++++++- 1 file changed, 40 insertions(+), 1 deletion(-) diff --git a/enterprise/pulsar.py b/enterprise/pulsar.py index 1921c70e..6d219a91 100644 --- a/enterprise/pulsar.py +++ b/enterprise/pulsar.py @@ -7,6 +7,7 @@ import os import astropy.units as u +import astropy.constants as const import numpy as np from ephem import Ecliptic, Equatorial @@ -322,7 +323,8 @@ def __init__(self, toas, model, sort=True, planets=True): self._pdist = self._get_pdist() self._raj, self._decj = self._get_radec(model) self._pos = self._get_pos() - self._planetssb = self._get_planetssb() + self._planetssb = self._get_planetssb(toas) + self._sunssb = self._get_sunssb(toas) # TODO: pos_t not currently implemented self._pos_t = np.zeros((len(self._toas), 3)) @@ -338,12 +340,49 @@ def _get_radec(self, model): elong, elat = model.ELONG.value, model.ELAT.value return self._get_radec_from_ecliptic(elong * d2r, elat * d2r) + def _get_ssb_lsec(toas,obs_pl,ssb_obs): + """Get the planet to SSB vector in lightseconds from Pint table""" + vec = toas.table[obs_pl]-toas.table[ssb_obs] + return (vec/const.c).to('s').value + def _get_planetssb(self): return np.zeros((len(self._toas), 9, 6)) def _get_sunssb(self): return np.zeros((len(self._toas), 6)) + def _get_planetssb(self, toas): + planetssb = None + if self.planets: + + planetssb = np.zeros((len(self._toas), 9, 6)) + # planetssb[:, 0, :] = self.t2pulsar.mercury_ssb + # planetssb[:, 1, :] = self.t2pulsar.venus_ssb + planetssb[:, 2, :3] = _get_ssb_lsec(toas,'obs_earth_pos','ssb_obs_pos') + # planetssb[:, 3, :] = self.t2pulsar.mars_ssb + planetssb[:, 4, :3] = _get_ssb_lsec(toas,'obs_jupiter_pos','ssb_obs_pos') + planetssb[:, 5, :3] = _get_ssb_lsec(toas,'obs_saturn_pos','ssb_obs_pos') + planetssb[:, 6, :3] = _get_ssb_lsec(toas,'obs_uranus_pos','ssb_obs_pos') + planetssb[:, 7, :3] = _get_ssb_lsec(toas,'obs_neptune_pos','ssb_obs_pos') + # planetssb[:, 8, :] = self.t2pulsar.pluto_ssb + + # if "ELONG" and "ELAT" in np.concatenate((t2pulsar.pars(), t2pulsar.pars(which="set"))): + # for ii in range(9): + # planetssb[:, ii, :3] = utils.ecl2eq_vec(planetssb[:, ii, :3]) + # planetssb[:, ii, 3:] = utils.ecl2eq_vec(planetssb[:, ii, 3:]) + return planetssb + + def _get_sunssb(self, toas): + sunssb = None + if self.planets: + sunssb = np.zeros((len(self._toas), 6)) + sunssb[:, :3] = _get_ssb_lsec(toas,'obs_sun_pos','ssb_obs_pos') + + # if "ELONG" and "ELAT" in np.concatenate((t2pulsar.pars(), t2pulsar.pars(which="set"))): + # sunssb[:, :3] = utils.ecl2eq_vec(sunssb[:, :3]) + # sunssb[:, 3:] = utils.ecl2eq_vec(sunssb[:, 3:]) + return sunssb + class Tempo2Pulsar(BasePulsar): def __init__(self, t2pulsar, sort=True, drop_t2pulsar=True, planets=True): From 399fcd0cdc1691a71830920967c780a7f37e18ac Mon Sep 17 00:00:00 2001 From: Hazboun6 Date: Mon, 18 Jan 2021 16:52:15 -0800 Subject: [PATCH 12/13] black wanted to change this --- travis_pypi_setup.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/travis_pypi_setup.py b/travis_pypi_setup.py index e787045d..ca8e0fc6 100644 --- a/travis_pypi_setup.py +++ b/travis_pypi_setup.py @@ -66,8 +66,7 @@ def fetch_public_key(repo): def prepend_line(filepath, line): - """Rewrite a file adding a line to its beginning. - """ + """Rewrite a file adding a line to its beginning.""" with open(filepath) as f: lines = f.readlines() From fdd7b59752407c5a01f7d73cac2c7074bd820556 Mon Sep 17 00:00:00 2001 From: Hazboun6 Date: Mon, 18 Jan 2021 16:52:40 -0800 Subject: [PATCH 13/13] first changes of PintPulsar for ssb vectors --- enterprise/pulsar.py | 45 ++++++++++++++++++++------------------------ tests/test_pulsar.py | 10 ++++++---- 2 files changed, 26 insertions(+), 29 deletions(-) diff --git a/enterprise/pulsar.py b/enterprise/pulsar.py index 4dae3153..20a8933e 100644 --- a/enterprise/pulsar.py +++ b/enterprise/pulsar.py @@ -323,8 +323,8 @@ def __init__(self, toas, model, sort=True, planets=True): self._pdist = self._get_pdist() self._raj, self._decj = self._get_radec(model) self._pos = self._get_pos() - self._planetssb = self._get_planetssb(toas) - self._sunssb = self._get_sunssb(toas) + self._planetssb = self._get_planetssb(toas, model) + self._sunssb = self._get_sunssb(toas, model) # TODO: pos_t not currently implemented self._pos_t = np.zeros((len(self._toas), 3)) @@ -361,47 +361,40 @@ def _get_radec(self, model): elong, elat = model.ELONG.value, model.ELAT.value return self._get_radec_from_ecliptic(elong * d2r, elat * d2r) - def _get_ssb_lsec(toas,obs_pl,ssb_obs): + def _get_ssb_lsec(self, toas, obs_planet): """Get the planet to SSB vector in lightseconds from Pint table""" - vec = toas.table[obs_pl]-toas.table[ssb_obs] - return (vec/const.c).to('s').value + vec = toas.table[obs_planet] + toas.table["ssb_obs_pos"] + return (vec / const.c).to("s").value - def _get_planetssb(self): - return np.zeros((len(self._toas), 9, 6)) - - def _get_sunssb(self): - return np.zeros((len(self._toas), 6)) - - def _get_planetssb(self, toas): + def _get_planetssb(self, toas, model): planetssb = None if self.planets: - planetssb = np.zeros((len(self._toas), 9, 6)) # planetssb[:, 0, :] = self.t2pulsar.mercury_ssb # planetssb[:, 1, :] = self.t2pulsar.venus_ssb - planetssb[:, 2, :3] = _get_ssb_lsec(toas,'obs_earth_pos','ssb_obs_pos') + planetssb[:, 2, :3] = self._get_ssb_lsec(toas, "obs_earth_pos") # planetssb[:, 3, :] = self.t2pulsar.mars_ssb - planetssb[:, 4, :3] = _get_ssb_lsec(toas,'obs_jupiter_pos','ssb_obs_pos') - planetssb[:, 5, :3] = _get_ssb_lsec(toas,'obs_saturn_pos','ssb_obs_pos') - planetssb[:, 6, :3] = _get_ssb_lsec(toas,'obs_uranus_pos','ssb_obs_pos') - planetssb[:, 7, :3] = _get_ssb_lsec(toas,'obs_neptune_pos','ssb_obs_pos') + planetssb[:, 4, :3] = self._get_ssb_lsec(toas, "obs_jupiter_pos") + planetssb[:, 5, :3] = self._get_ssb_lsec(toas, "obs_saturn_pos") + planetssb[:, 6, :3] = self._get_ssb_lsec(toas, "obs_uranus_pos") + planetssb[:, 7, :3] = self._get_ssb_lsec(toas, "obs_neptune_pos") # planetssb[:, 8, :] = self.t2pulsar.pluto_ssb - # if "ELONG" and "ELAT" in np.concatenate((t2pulsar.pars(), t2pulsar.pars(which="set"))): + # if hasattr(model, "ELAT") and hasattr(model, "ELONG"): # for ii in range(9): # planetssb[:, ii, :3] = utils.ecl2eq_vec(planetssb[:, ii, :3]) - # planetssb[:, ii, 3:] = utils.ecl2eq_vec(planetssb[:, ii, 3:]) + # # planetssb[:, ii, 3:] = utils.ecl2eq_vec(planetssb[:, ii, 3:]) return planetssb - def _get_sunssb(self, toas): + def _get_sunssb(self, toas, model): sunssb = None if self.planets: sunssb = np.zeros((len(self._toas), 6)) - sunssb[:, :3] = _get_ssb_lsec(toas,'obs_sun_pos','ssb_obs_pos') + sunssb[:, :3] = self._get_ssb_lsec(toas, "obs_sun_pos") - # if "ELONG" and "ELAT" in np.concatenate((t2pulsar.pars(), t2pulsar.pars(which="set"))): + # if hasattr(model, "ELAT") and hasattr(model, "ELONG"): # sunssb[:, :3] = utils.ecl2eq_vec(sunssb[:, :3]) - # sunssb[:, 3:] = utils.ecl2eq_vec(sunssb[:, 3:]) + # # sunssb[:, 3:] = utils.ecl2eq_vec(sunssb[:, 3:]) return sunssb @@ -568,7 +561,9 @@ def Pulsar(*args, **kwargs): if timing_package.lower() == "pint": if (clk is not None) and (bipm_version is None): bipm_version = clk.split("(")[1][:-1] - model, toas = get_model_and_toas(relparfile, reltimfile, ephem=ephem, bipm_version=bipm_version) + model, toas = get_model_and_toas( + relparfile, reltimfile, ephem=ephem, bipm_version=bipm_version, planets=planets + ) os.chdir(cwd) return PintPulsar(toas, model, sort=sort, planets=planets) diff --git a/tests/test_pulsar.py b/tests/test_pulsar.py index d4128d90..95cad51e 100644 --- a/tests/test_pulsar.py +++ b/tests/test_pulsar.py @@ -151,13 +151,15 @@ def setUpClass(cls): # exclude tests pending implementation of .stoas, .dm, .dmx in PintPulsar def test_stoas(self): - pass + assert hasattr(self.psr, "stoas") def test_dm(self): - pass + assert hasattr(self.psr, "dm") def test_planetssb(self): - pass + """Place holder for filter_data tests.""" + assert hasattr(self.psr, "planetssb") def test_sunssb(self): - pass + """Place holder for filter_data tests.""" + assert hasattr(self.psr, "sunssb")