From 186b718d1754280c268706bc3bf671483d39fd9c Mon Sep 17 00:00:00 2001 From: Jeremy Baier Date: Wed, 5 Apr 2023 08:51:03 -0700 Subject: [PATCH 01/32] initial commit --- hasasia/sensitivity.py | 25 +++++++++++++++++++++++++ 1 file changed, 25 insertions(+) diff --git a/hasasia/sensitivity.py b/hasasia/sensitivity.py index b60378f..ef28ddf 100644 --- a/hasasia/sensitivity.py +++ b/hasasia/sensitivity.py @@ -1189,6 +1189,31 @@ def red_noise_powerlaw(A, freqs, gamma=None, alpha=None): return A**2*(freqs/fyr)**(-gamma)/(12*np.pi**2) * yr_sec**3 +def red_noise_realistic_curve(A, freqs, gamma=None, alpha=None): + r""" + jeremy + Add power law red noise to the prefit residual power spectral density. + As :math:`P=A^2(f/fyr)^{-\gamma}` + + Parameters + ---------- + A : float + Amplitude of red noise. + + gamma : float + Spectral index of red noise powerlaw. + + freqs : array + Frequencies at which to calculate the red noise power law. + """ + if gamma is None and alpha is not None: + gamma = 3-2*alpha + elif ((gamma is None and alpha is None) + or (gamma is not None and alpha is not None)): + ValueError('Must specify one version of spectral index.') + + return A**2*(freqs/fyr)**(-gamma)/(12*np.pi**2) * yr_sec**3 + def S_h(A, alpha, freqs): r""" Add power law red noise to the prefit residual power spectral density. From 3c79ceb77a0940b32e2c3b22aa9e12e18c16d1f1 Mon Sep 17 00:00:00 2001 From: Jeremy Baier Date: Wed, 31 May 2023 01:59:41 -0700 Subject: [PATCH 02/32] adding gitignore and new function --- .gitignore | 2 +- hasasia/sensitivity.py | 27 +++++++++++---------------- 2 files changed, 12 insertions(+), 17 deletions(-) diff --git a/.gitignore b/.gitignore index 84229f4..a06d623 100644 --- a/.gitignore +++ b/.gitignore @@ -2,7 +2,7 @@ __pycache__/ *.py[cod] *$py.class - +**/.DS_Store # C extensions *.so diff --git a/hasasia/sensitivity.py b/hasasia/sensitivity.py index ef28ddf..e7a0c85 100644 --- a/hasasia/sensitivity.py +++ b/hasasia/sensitivity.py @@ -1189,30 +1189,25 @@ def red_noise_powerlaw(A, freqs, gamma=None, alpha=None): return A**2*(freqs/fyr)**(-gamma)/(12*np.pi**2) * yr_sec**3 -def red_noise_realistic_curve(A, freqs, gamma=None, alpha=None): +def psd_from_holodeck_realization(hc_bg, freqs): r""" - jeremy - Add power law red noise to the prefit residual power spectral density. - As :math:`P=A^2(f/fyr)^{-\gamma}` + Calculate the power spectral density with given strain and frequency. Parameters ---------- - A : float - Amplitude of red noise. - - gamma : float - Spectral index of red noise powerlaw. + hc_bg : array + characteristic strain of background at each frequency freqs : array - Frequencies at which to calculate the red noise power law. + Frequency + + Returns + ------- + S_h : array + the power spectral density from the background """ - if gamma is None and alpha is not None: - gamma = 3-2*alpha - elif ((gamma is None and alpha is None) - or (gamma is not None and alpha is not None)): - ValueError('Must specify one version of spectral index.') - return A**2*(freqs/fyr)**(-gamma)/(12*np.pi**2) * yr_sec**3 + return hc_bg**2 / (12 * np.pi**2 * freqs[:,np.newaxis]**3) def S_h(A, alpha, freqs): r""" From ad3c628a81407f4408951a7c0a5392ee436dbe61 Mon Sep 17 00:00:00 2001 From: Jeremy Baier Date: Sat, 3 Jun 2023 13:38:39 -0700 Subject: [PATCH 03/32] updating psd function; adding detection probab ;) --- hasasia/sensitivity.py | 14 +++++++------- hasasia/utils.py | 18 ++++++++++++------ 2 files changed, 19 insertions(+), 13 deletions(-) diff --git a/hasasia/sensitivity.py b/hasasia/sensitivity.py index e7a0c85..76b37d3 100644 --- a/hasasia/sensitivity.py +++ b/hasasia/sensitivity.py @@ -1189,25 +1189,25 @@ def red_noise_powerlaw(A, freqs, gamma=None, alpha=None): return A**2*(freqs/fyr)**(-gamma)/(12*np.pi**2) * yr_sec**3 -def psd_from_holodeck_realization(hc_bg, freqs): +def psd_from_background_realization(background_hc, freqs): r""" - Calculate the power spectral density with given strain and frequency. + Calculate the power spectral density with given background strain and frequency binning. Parameters ---------- - hc_bg : array - characteristic strain of background at each frequency + background_hc : array + Characteristic strain (hc) of background at each frequency. freqs : array - Frequency + Frequency bins over which the background is stored. Returns ------- S_h : array - the power spectral density from the background + the power spectral density of the background """ - return hc_bg**2 / (12 * np.pi**2 * freqs[:,np.newaxis]**3) + return background_hc**2 / (12 * np.pi**2 * freqs[:,np.newaxis]**3) def S_h(A, alpha, freqs): r""" diff --git a/hasasia/utils.py b/hasasia/utils.py index 8651e28..505173e 100644 --- a/hasasia/utils.py +++ b/hasasia/utils.py @@ -87,7 +87,7 @@ def pdf_F_signal(F, snr, Npsrs=None): N = int(4 * Npsrs) return ss.ncx2.pdf(2*F, N, snr**2) -def fdp(F0, snr, Npsrs=None, sky_ave=False): +def fdp(F0, snr, Npsrs=None, iota_psi_ave=False): ''' False detection probability of the F-statistic Use None for the Fe statistic and the number of pulsars for the Fp stat. @@ -96,21 +96,27 @@ def fdp(F0, snr, Npsrs=None, sky_ave=False): N = 4 elif isinstance(Npsrs,int): N = int(4 * Npsrs) - if sky_ave: + if iota_psi_ave: return ss.chi2.cdf(2*F0, df=N, loc=snr**2) else: return ss.ncx2.cdf(2*F0, df=N, nc=snr**2) +def detection_prob(F0, snr, Npsrs=None, iota_psi_ave=False): + ''' + Detection probability of the F-statistic + Use None for the Fe and the number of pulsars for the Fp stat. + ''' + return 1 - fdp(F0, snr, Npsrs, iota_psi_ave) def _solve_F_given_fap(fap0=0.003, Npsrs=None): return sopt.fsolve(lambda F :fap(F, Npsrs=Npsrs)-fap0, 10) -def _solve_F_given_fdp_snr(fdp0=0.05, snr=3, Npsrs=None, sky_ave=False): +def _solve_F_given_fdp_snr(fdp0=0.05, snr=3, Npsrs=None, iota_psi_ave=False): Npsrs = 1 if Npsrs is None else Npsrs F0 = (4*Npsrs+snr**2)/2 - return sopt.fsolve(lambda F :fdp(F, snr, Npsrs=Npsrs, sky_ave=sky_ave)-fdp0, F0) + return sopt.fsolve(lambda F :fdp(F, snr, Npsrs=Npsrs, iota_psi_ave=iota_psi_ave)-fdp0, F0) -def _solve_snr_given_fdp_F(fdp0=0.05, F=3, Npsrs=None, sky_ave=False): +def _solve_snr_given_fdp_F(fdp0=0.05, F=3, Npsrs=None, iota_psi_ave=False): Npsrs = 1 if Npsrs is None else Npsrs snr0 = np.sqrt(2*F-4*Npsrs) - return sopt.fsolve(lambda snr :fdp(F, snr, Npsrs=Npsrs, sky_ave=sky_ave)-fdp0, snr0) + return sopt.fsolve(lambda snr :fdp(F, snr, Npsrs=Npsrs, iota_psi_ave=iota_psi_ave)-fdp0, snr0) From 07cb940850b70b8a09941e4ffefa2d48976006e1 Mon Sep 17 00:00:00 2001 From: Jeremy Baier Date: Fri, 30 Jun 2023 17:48:56 -0700 Subject: [PATCH 04/32] forget what this was --- hasasia/utils.py | 18 +++++++++++++----- 1 file changed, 13 insertions(+), 5 deletions(-) diff --git a/hasasia/utils.py b/hasasia/utils.py index 505173e..df38019 100644 --- a/hasasia/utils.py +++ b/hasasia/utils.py @@ -87,9 +87,9 @@ def pdf_F_signal(F, snr, Npsrs=None): N = int(4 * Npsrs) return ss.ncx2.pdf(2*F, N, snr**2) -def fdp(F0, snr, Npsrs=None, iota_psi_ave=False): +def false_dismissal__prob(F0, snr, Npsrs=None, iota_psi_ave=False): ''' - False detection probability of the F-statistic + False dismissal probability of the F-statistic Use None for the Fe statistic and the number of pulsars for the Fp stat. ''' if Npsrs is None: @@ -106,7 +106,7 @@ def detection_prob(F0, snr, Npsrs=None, iota_psi_ave=False): Detection probability of the F-statistic Use None for the Fe and the number of pulsars for the Fp stat. ''' - return 1 - fdp(F0, snr, Npsrs, iota_psi_ave) + return 1 - false_dismissal__prob(F0, snr, Npsrs, iota_psi_ave) def _solve_F_given_fap(fap0=0.003, Npsrs=None): return sopt.fsolve(lambda F :fap(F, Npsrs=Npsrs)-fap0, 10) @@ -114,9 +114,17 @@ def _solve_F_given_fap(fap0=0.003, Npsrs=None): def _solve_F_given_fdp_snr(fdp0=0.05, snr=3, Npsrs=None, iota_psi_ave=False): Npsrs = 1 if Npsrs is None else Npsrs F0 = (4*Npsrs+snr**2)/2 - return sopt.fsolve(lambda F :fdp(F, snr, Npsrs=Npsrs, iota_psi_ave=iota_psi_ave)-fdp0, F0) + return sopt.fsolve(lambda F :false_dismissal__prob(F, snr, Npsrs=Npsrs, iota_psi_ave=iota_psi_ave)-fdp0, F0) def _solve_snr_given_fdp_F(fdp0=0.05, F=3, Npsrs=None, iota_psi_ave=False): Npsrs = 1 if Npsrs is None else Npsrs snr0 = np.sqrt(2*F-4*Npsrs) - return sopt.fsolve(lambda snr :fdp(F, snr, Npsrs=Npsrs, iota_psi_ave=iota_psi_ave)-fdp0, snr0) + return sopt.fsolve(lambda snr :false_dismissal__prob(F, snr, Npsrs=Npsrs, iota_psi_ave=iota_psi_ave)-fdp0, snr0) + +def _solve_F0_given_SNR(snr=3, Npsrs=None): + ''' + Returns the F0 (Fe stat threshold for a specified SNR) + Use None for the Fe and the number of pulsars for the Fp stat. + ''' + Npsrs = 1 if Npsrs is None else Npsrs + return 0.5*(4.*Npsrs+snr**2.) From 0e477f831df8b9fd0b1fd2093390d311b6c0b138 Mon Sep 17 00:00:00 2001 From: Jeremy Baier Date: Fri, 21 Jul 2023 13:18:21 -0700 Subject: [PATCH 05/32] delete underscore --- hasasia/utils.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/hasasia/utils.py b/hasasia/utils.py index df38019..4ff4072 100644 --- a/hasasia/utils.py +++ b/hasasia/utils.py @@ -87,7 +87,7 @@ def pdf_F_signal(F, snr, Npsrs=None): N = int(4 * Npsrs) return ss.ncx2.pdf(2*F, N, snr**2) -def false_dismissal__prob(F0, snr, Npsrs=None, iota_psi_ave=False): +def false_dismissal_prob(F0, snr, Npsrs=None, iota_psi_ave=False): ''' False dismissal probability of the F-statistic Use None for the Fe statistic and the number of pulsars for the Fp stat. From d2c26cbe8bbdca2490275e945249317c22f67034 Mon Sep 17 00:00:00 2001 From: Jeremy Baier Date: Fri, 28 Jul 2023 16:47:27 -0700 Subject: [PATCH 06/32] edit typo in utils --- hasasia/utils.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/hasasia/utils.py b/hasasia/utils.py index 4ff4072..c8171d5 100644 --- a/hasasia/utils.py +++ b/hasasia/utils.py @@ -106,7 +106,7 @@ def detection_prob(F0, snr, Npsrs=None, iota_psi_ave=False): Detection probability of the F-statistic Use None for the Fe and the number of pulsars for the Fp stat. ''' - return 1 - false_dismissal__prob(F0, snr, Npsrs, iota_psi_ave) + return 1 - false_dismissal_prob(F0, snr, Npsrs, iota_psi_ave) def _solve_F_given_fap(fap0=0.003, Npsrs=None): return sopt.fsolve(lambda F :fap(F, Npsrs=Npsrs)-fap0, 10) @@ -114,12 +114,12 @@ def _solve_F_given_fap(fap0=0.003, Npsrs=None): def _solve_F_given_fdp_snr(fdp0=0.05, snr=3, Npsrs=None, iota_psi_ave=False): Npsrs = 1 if Npsrs is None else Npsrs F0 = (4*Npsrs+snr**2)/2 - return sopt.fsolve(lambda F :false_dismissal__prob(F, snr, Npsrs=Npsrs, iota_psi_ave=iota_psi_ave)-fdp0, F0) + return sopt.fsolve(lambda F :false_dismissal_prob(F, snr, Npsrs=Npsrs, iota_psi_ave=iota_psi_ave)-fdp0, F0) def _solve_snr_given_fdp_F(fdp0=0.05, F=3, Npsrs=None, iota_psi_ave=False): Npsrs = 1 if Npsrs is None else Npsrs snr0 = np.sqrt(2*F-4*Npsrs) - return sopt.fsolve(lambda snr :false_dismissal__prob(F, snr, Npsrs=Npsrs, iota_psi_ave=iota_psi_ave)-fdp0, snr0) + return sopt.fsolve(lambda snr :false_dismissal_prob(F, snr, Npsrs=Npsrs, iota_psi_ave=iota_psi_ave)-fdp0, snr0) def _solve_F0_given_SNR(snr=3, Npsrs=None): ''' From ec377d227a8738774eb85fbbb2705f81e79275c5 Mon Sep 17 00:00:00 2001 From: Jeremy Baier Date: Thu, 1 Feb 2024 20:19:57 -0800 Subject: [PATCH 07/32] adding git ignore. not sure if i should --- .gitignore | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/.gitignore b/.gitignore index 84229f4..ef20f15 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,7 @@ +#MacBook stuff +.DS_Store +.gitignore + # Byte-compiled / optimized / DLL files __pycache__/ *.py[cod] From 072f8c050f8a0a441330923fed5a573bc5f5e64f Mon Sep 17 00:00:00 2001 From: Jeremy Baier Date: Thu, 1 Feb 2024 21:05:36 -0800 Subject: [PATCH 08/32] the commit before i add enterprise.pulsar stuff --- hasasia/sensitivity.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/hasasia/sensitivity.py b/hasasia/sensitivity.py index 76b37d3..b3fd730 100644 --- a/hasasia/sensitivity.py +++ b/hasasia/sensitivity.py @@ -315,12 +315,13 @@ class Pulsar(object): Earth-pulsar distance. Default units is kpc. """ - def __init__(self, toas, toaerrs, phi=None, theta=None, + def __init__(self, toas, toaerrs, phi=None, theta=None, name=None, designmatrix=None, N=None, pdist=1.0*u.kpc): self.toas = toas self.toaerrs = toaerrs self.phi = phi self.theta = theta + self.name = name self.pdist = make_quant(pdist,'kpc') if N is None: From 8c46ae83edd05caa0fbf14f011e7b80d9df16337 Mon Sep 17 00:00:00 2001 From: Jeremy Baier Date: Thu, 1 Feb 2024 21:15:59 -0800 Subject: [PATCH 09/32] added filter_data to hsen --- hasasia/sensitivity.py | 25 +++++++++++++++++++++++++ 1 file changed, 25 insertions(+) diff --git a/hasasia/sensitivity.py b/hasasia/sensitivity.py index b3fd730..8e28831 100644 --- a/hasasia/sensitivity.py +++ b/hasasia/sensitivity.py @@ -335,6 +335,31 @@ def __init__(self, toas, toaerrs, phi=None, theta=None, name=None, else: self.designmatrix = designmatrix + def filter_data(self, start_time=1, end_time=None): + """ + Parameters + ========== + start_time - float + MJD at which to begin data subset. + Defaults to 1 MJD. + end_time - float + MJD at which to end data subset. + + Filter data to create a time-slice of overall dataset. + Function adapted from enterprise.BasePulsar() class. + """ + if start_time is None and end_time is None: + mask = np.ones(self._toas.shape, dtype=bool) + else: + mask = np.logical_and(self._toas >= start_time * 86400, self._toas <= end_time * 86400) + + self._toas = self._toas[mask] + self._toaerrs = self._toaerrs[mask] + + self._designmatrix = self._designmatrix[mask, :] + dmx_mask = np.sum(self._designmatrix, axis=0) != 0.0 + self._designmatrix = self._designmatrix[:, dmx_mask] + @property def G(self): """Inverse Noise Weighted Transmission Function.""" From 0d1c5717b40c39a8ca11d848f74de331ed468298 Mon Sep 17 00:00:00 2001 From: Jeremy Baier Date: Fri, 2 Feb 2024 00:09:18 -0800 Subject: [PATCH 10/32] probably like some print statements --- hasasia/sensitivity.py | 20 +++++++++++--------- 1 file changed, 11 insertions(+), 9 deletions(-) diff --git a/hasasia/sensitivity.py b/hasasia/sensitivity.py index 8e28831..6e4a92e 100644 --- a/hasasia/sensitivity.py +++ b/hasasia/sensitivity.py @@ -335,13 +335,12 @@ def __init__(self, toas, toaerrs, phi=None, theta=None, name=None, else: self.designmatrix = designmatrix - def filter_data(self, start_time=1, end_time=None): + def filter_data(self, start_time=None, end_time=None): """ Parameters ========== start_time - float MJD at which to begin data subset. - Defaults to 1 MJD. end_time - float MJD at which to end data subset. @@ -349,16 +348,19 @@ def filter_data(self, start_time=1, end_time=None): Function adapted from enterprise.BasePulsar() class. """ if start_time is None and end_time is None: - mask = np.ones(self._toas.shape, dtype=bool) + mask = np.ones(self.toas.shape, dtype=bool) else: - mask = np.logical_and(self._toas >= start_time * 86400, self._toas <= end_time * 86400) + mask = np.logical_and(self.toas >= start_time * 86400, self.toas <= end_time * 86400) - self._toas = self._toas[mask] - self._toaerrs = self._toaerrs[mask] + self.toas = self.toas[mask] + self.toaerrs = self.toaerrs[mask] - self._designmatrix = self._designmatrix[mask, :] - dmx_mask = np.sum(self._designmatrix, axis=0) != 0.0 - self._designmatrix = self._designmatrix[:, dmx_mask] + self.designmatrix = self.designmatrix[mask, :] + #dmx_mask = np.sum(self.designmatrix, axis=0) != 0.0 + #self.designmatrix = self.designmatrix[:, dmx_mask] + self._G = G_matrix(designmatrix=self.designmatrix) + print("dim designmatrix = ", self.designmatrix.shape) + print("dim G matrxi = ", self.G.shape) @property def G(self): From 311da948902476fe8d4ecfa02945cc27946bc55c Mon Sep 17 00:00:00 2001 From: Jeremy Baier Date: Tue, 6 Feb 2024 20:39:15 -0800 Subject: [PATCH 11/32] changes for data slicing --- hasasia/sensitivity.py | 18 ++++++++++------ hasasia/sim.py | 10 +++++---- hasasia/utils.py | 49 ++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 67 insertions(+), 10 deletions(-) diff --git a/hasasia/sensitivity.py b/hasasia/sensitivity.py index 6e4a92e..43d3991 100644 --- a/hasasia/sensitivity.py +++ b/hasasia/sensitivity.py @@ -354,13 +354,15 @@ def filter_data(self, start_time=None, end_time=None): self.toas = self.toas[mask] self.toaerrs = self.toaerrs[mask] + self.N = self.N[mask, :][:, mask] - self.designmatrix = self.designmatrix[mask, :] + self.designmatrix = create_design_matrix(self.toas, RADEC=True, PROPER=True, PX=True) + #self.designmatrix = self.designmatrix[mask, :] #dmx_mask = np.sum(self.designmatrix, axis=0) != 0.0 #self.designmatrix = self.designmatrix[:, dmx_mask] self._G = G_matrix(designmatrix=self.designmatrix) - print("dim designmatrix = ", self.designmatrix.shape) - print("dim G matrxi = ", self.G.shape) + #print("dim designmatrix = ", self.designmatrix.shape) + #print("dim G matrxi = ", self.G.shape) @property def G(self): @@ -1050,9 +1052,13 @@ def get_Tspan(psrs): """ Returns the total timespan from a list or arry of Pulsar objects, psrs. """ - last = np.amax([p.toas.max() for p in psrs]) - first = np.amin([p.toas.min() for p in psrs]) - return last - first + try: + last = np.amax([p.toas.max() for p in psrs]) + first = np.amin([p.toas.min() for p in psrs]) + tspan = last-first + except ValueError: + tspan = 0 + return tspan def get_TspanIJ(psr1,psr2): """ diff --git a/hasasia/sim.py b/hasasia/sim.py index 2ad5d54..44ef18d 100644 --- a/hasasia/sim.py +++ b/hasasia/sim.py @@ -11,8 +11,8 @@ yr_sec = 365.25*24*3600 def sim_pta(timespan, cad, sigma, phi, theta, Npsrs=None, - A_rn=None, alpha=None, freqs=None, uneven=False, A_gwb=None, - fast=True, + A_rn=None, alpha=None, freqs=None, uneven=False, + A_gwb=None, gamma_gwb = 13/3, fast=True, psr_names = None, kwastro={'RADEC':True, 'PROPER':True, 'PX':True}): """ Make a simulated pulsar timing array. Using the available parameters, @@ -98,6 +98,8 @@ def sim_pta(timespan, cad, sigma, phi, theta, Npsrs=None, psrs = [] err_dim = pars['sigma'].ndim Timespan = np.amax(pars['timespan']) + if psr_names is None: + psr_names = [None for _ in range(Npsrs)] for ii in range(Npsrs): tspan = pars['timespan'][ii] Ntoas = int(np.floor(tspan*pars['cad'][ii])) @@ -121,13 +123,13 @@ def sim_pta(timespan, cad, sigma, phi, theta, Npsrs=None, if A_gwb is not None: gwb = red_noise_powerlaw(A=A_gwb, - alpha=-2/3., + alpha=gamma_gwb, freqs=freqs) N += corr_from_psd(freqs=freqs, psd=gwb, toas=toas, fast=fast) M = create_design_matrix(toas, **kwastro) - p = Pulsar(toas, toaerrs, phi=pars['phi'][ii], + p = Pulsar(toas, toaerrs, phi=pars['phi'][ii], name=psr_names[ii], theta=pars['theta'][ii], designmatrix=M, N=N) psrs.append(p) diff --git a/hasasia/utils.py b/hasasia/utils.py index c8171d5..1687f3c 100644 --- a/hasasia/utils.py +++ b/hasasia/utils.py @@ -5,6 +5,8 @@ import scipy.special as ssp import scipy.integrate as si import scipy.stats as ss +from astropy.coordinates import SkyCoord +import astropy.units as u __all__ = ['create_design_matrix', 'fap', @@ -128,3 +130,50 @@ def _solve_F0_given_SNR(snr=3, Npsrs=None): ''' Npsrs = 1 if Npsrs is None else Npsrs return 0.5*(4.*Npsrs+snr**2.) + +def theta_phi_to_SkyCoord(theta, phi): + """ + Parameters + ========== + phi : float, array of floats + The celestial longitude in solar system coordinates. + theta : float, array of floats + The celestial lattitude in solar system coordinates. + + Returns + ======= + skycoord - astropy.SkyCoord object + Can use this to convert to ra, dec, etc. + (e.g. SkyCoord.ra.deg) + + Converts an inputed longitude and lattitude into an `astropy.SkyCoord` object. + + """ + + return SkyCoord(phi*u.rad, ( theta - np.pi/2 )*u.rad) + +def skycoord_to_Jname(skycoord): + """ + Parameters + ========== + skycoord - astropy.SkyCoord object + Can use `theta_phi_to_SkyCoord()` to get this. + + Returns + ======= + Jname - string, array of strings + The traditional Jname of a pulsar with given coordinates + (eg. 'J1713+0747') + + """ + coord_pieces = [ + str(skycoord.ra.hms[0]).split('.')[0], + str(skycoord.ra.hms[1]).split('.')[0], + str(abs(skycoord.dec.hms[0])).split('.')[0], + str(abs(skycoord.dec.hms[1])).split('.')[0] + ] + sign = ['-','+'][int(skycoord.dec.hms[0]>0)] + for i, piece in enumerate(coord_pieces): + if len(str(abs(int(piece)))) < 2: + coord_pieces[i] = '0' + str(piece) + return 'J' + coord_pieces[0] + coord_pieces[1] + sign + coord_pieces[2] +coord_pieces[3] \ No newline at end of file From a52312fa1879f1867afec84b93ff5536f04d477f Mon Sep 17 00:00:00 2001 From: Jeremy Baier Date: Tue, 6 Feb 2024 22:27:42 -0800 Subject: [PATCH 12/32] adding function to altar cadence --- hasasia/sensitivity.py | 39 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 39 insertions(+) diff --git a/hasasia/sensitivity.py b/hasasia/sensitivity.py index 43d3991..4e15e74 100644 --- a/hasasia/sensitivity.py +++ b/hasasia/sensitivity.py @@ -354,6 +354,45 @@ def filter_data(self, start_time=None, end_time=None): self.toas = self.toas[mask] self.toaerrs = self.toaerrs[mask] + self.N = self.N[mask, :][:, mask] + + self.designmatrix = create_design_matrix(self.toas, RADEC=True, PROPER=True, PX=True) + #self.designmatrix = self.designmatrix[mask, :] + #dmx_mask = np.sum(self.designmatrix, axis=0) != 0.0 + #self.designmatrix = self.designmatrix[:, dmx_mask] + self._G = G_matrix(designmatrix=self.designmatrix) + #print("dim designmatrix = ", self.designmatrix.shape) + #print("dim G matrxi = ", self.G.shape) + + def change_cadence(self, start_time=None, end_time=None, + cadence=2): + """ + ***************** + WORK IN PROGRESSS + ***************** + Parameters + ========== + start_time - float + MJD at which to begin altared cadence. + end_time - float + MJD at which to end altared cadence. + cadence - float + factor by which to change the observing cadence. + + + Change observing cadence in a given time range. + Recalculate pulsar properties. + """ + if start_time is None and end_time is None: + mask = np.ones(self.toas.shape, dtype=bool) + else: + mask = np.logical_and(self.toas >= start_time * 86400, self.toas <= end_time * 86400) + # change cadence by changing number of toas + self.toas = self.toas[mask] + self.toaerrs = self.toaerrs[mask] + + + self.N = self.N[mask, :][:, mask] self.designmatrix = create_design_matrix(self.toas, RADEC=True, PROPER=True, PX=True) From 14f72db4bbc0ca472d5a4221beef1e2b714dd8ad Mon Sep 17 00:00:00 2001 From: Jeremy Baier Date: Sun, 11 Feb 2024 21:42:41 -0800 Subject: [PATCH 13/32] first working draft of change cadence function; additional changes to sim_pta required --- hasasia/sensitivity.py | 71 +++++++++++++++++++++++++++++++----------- hasasia/sim.py | 14 +++++++-- 2 files changed, 63 insertions(+), 22 deletions(-) diff --git a/hasasia/sensitivity.py b/hasasia/sensitivity.py index 4e15e74..8c53034 100644 --- a/hasasia/sensitivity.py +++ b/hasasia/sensitivity.py @@ -316,13 +316,16 @@ class Pulsar(object): """ def __init__(self, toas, toaerrs, phi=None, theta=None, name=None, - designmatrix=None, N=None, pdist=1.0*u.kpc): + designmatrix=None, N=None, pdist=1.0*u.kpc, A_rn=None, + alpha=None,): self.toas = toas self.toaerrs = toaerrs self.phi = phi self.theta = theta self.name = name self.pdist = make_quant(pdist,'kpc') + self.A_rn = A_rn + self.alpha = alpha if N is None: self.N = np.diag(toaerrs**2) #N ==> weights @@ -364,8 +367,10 @@ def filter_data(self, start_time=None, end_time=None): #print("dim designmatrix = ", self.designmatrix.shape) #print("dim G matrxi = ", self.G.shape) - def change_cadence(self, start_time=None, end_time=None, - cadence=2): + def change_cadence(self, start_time=0, end_time=1_000_000, + cadence=None, cadence_factor=2, uneven=False, + A_gwb=None, gamma_gwb=13/3., freqs=None, + fast=True,): """ ***************** WORK IN PROGRESSS @@ -383,25 +388,53 @@ def change_cadence(self, start_time=None, end_time=None, Change observing cadence in a given time range. Recalculate pulsar properties. """ - if start_time is None and end_time is None: - mask = np.ones(self.toas.shape, dtype=bool) + mask_before = self.toas <= start_time * 86400 + mask_after = self.toas >= end_time * 86400 + old_Ntoas = np.sum( + np.logical_and(self.toas >= start_time * 86400, + self.toas <= end_time * 86400) + ) + # store the old toas and errors + old_toas = self.toas + old_toaerrs = self.toaerrs + # calculate old cadence then modified cadence + if start_time < min(old_toas)/84600: + start_time = min(old_toas)/84600 + if end_time < min(old_toas)/84600: + print("trying to change non-existant campaign") + return 0 + duration = end_time - start_time # in MJD + old_cadence = old_Ntoas / duration * 365.25 # cad is Ntoas/year + if cadence is not None: + new_cadence = cadence else: - mask = np.logical_and(self.toas >= start_time * 86400, self.toas <= end_time * 86400) - # change cadence by changing number of toas - self.toas = self.toas[mask] - self.toaerrs = self.toaerrs[mask] - - - - self.N = self.N[mask, :][:, mask] - + new_cadence = old_cadence * cadence_factor + # create new toas and toa errors + campaign_Ntoas = int(np.floor( duration / 365.25 * new_cadence )) + campaign_toas = np.linspace(start_time, end_time, campaign_Ntoas) * 86400 + print(duration, old_Ntoas, new_cadence, campaign_Ntoas, len(campaign_toas)) + if uneven: + dt = duration / campaign_Ntoas / 4 * yr_sec + campaign_toas += np.random.uniform(-dt, dt, size=campaign_Ntoas) + self.toas = np.concatenate([old_toas[mask_before], campaign_toas, old_toas[mask_after]]) + campaign_toaerrs = np.mean(old_toaerrs)*np.ones(campaign_Ntoas) + self.toaerrs = np.concatenate([old_toaerrs[mask_before], campaign_toaerrs, old_toaerrs[mask_after]]) + # recalculate N, designmatrix, G, + N = np.diag(self.toaerrs**2) + if self.A_rn is not None: + plaw = red_noise_powerlaw(A=10**self.A_rn, + alpha=self.alpha, + freqs=freqs) + N += corr_from_psd(freqs=freqs, psd=plaw, toas=self.toas, fast=fast) + + if A_gwb is not None: + gwb = red_noise_powerlaw(A=A_gwb, + alpha=gamma_gwb, + freqs=freqs) + N += corr_from_psd(freqs=freqs, psd=gwb, toas=self.toas, fast=fast) self.designmatrix = create_design_matrix(self.toas, RADEC=True, PROPER=True, PX=True) - #self.designmatrix = self.designmatrix[mask, :] - #dmx_mask = np.sum(self.designmatrix, axis=0) != 0.0 - #self.designmatrix = self.designmatrix[:, dmx_mask] self._G = G_matrix(designmatrix=self.designmatrix) - #print("dim designmatrix = ", self.designmatrix.shape) - #print("dim G matrxi = ", self.G.shape) + self.N = N @property def G(self): diff --git a/hasasia/sim.py b/hasasia/sim.py index 44ef18d..eae59fe 100644 --- a/hasasia/sim.py +++ b/hasasia/sim.py @@ -128,9 +128,17 @@ def sim_pta(timespan, cad, sigma, phi, theta, Npsrs=None, N += corr_from_psd(freqs=freqs, psd=gwb, toas=toas, fast=fast) M = create_design_matrix(toas, **kwastro) - - p = Pulsar(toas, toaerrs, phi=pars['phi'][ii], name=psr_names[ii], - theta=pars['theta'][ii], designmatrix=M, N=N) + # FIXME to always loop over A_rn + if 'A_rn' in keys: + A_rn = pars['A_rn'][ii] + alpha = pars['alpha'][ii] + else: + A_rn = None + + p = Pulsar(toas, toaerrs, name=psr_names[ii], + phi=pars['phi'][ii], theta=pars['theta'][ii], + A_rn=A_rn, alpha=alpha, + designmatrix=M, N=N) psrs.append(p) return psrs From dbfe8b9c66dbf6e32794f9b7dab1261dbf27617c Mon Sep 17 00:00:00 2001 From: Jeremy Baier Date: Mon, 19 Feb 2024 12:39:38 -0800 Subject: [PATCH 14/32] pushing changes --- hasasia/sensitivity.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/hasasia/sensitivity.py b/hasasia/sensitivity.py index 8c53034..79e0b72 100644 --- a/hasasia/sensitivity.py +++ b/hasasia/sensitivity.py @@ -412,14 +412,15 @@ def change_cadence(self, start_time=0, end_time=1_000_000, # create new toas and toa errors campaign_Ntoas = int(np.floor( duration / 365.25 * new_cadence )) campaign_toas = np.linspace(start_time, end_time, campaign_Ntoas) * 86400 - print(duration, old_Ntoas, new_cadence, campaign_Ntoas, len(campaign_toas)) if uneven: - dt = duration / campaign_Ntoas / 4 * yr_sec + # FIXME check this with jeff to see what he was going for + # in sim_pta() + dt = duration / campaign_Ntoas / 8 * yr_sec campaign_toas += np.random.uniform(-dt, dt, size=campaign_Ntoas) self.toas = np.concatenate([old_toas[mask_before], campaign_toas, old_toas[mask_after]]) campaign_toaerrs = np.mean(old_toaerrs)*np.ones(campaign_Ntoas) self.toaerrs = np.concatenate([old_toaerrs[mask_before], campaign_toaerrs, old_toaerrs[mask_after]]) - # recalculate N, designmatrix, G, + # recalculate N, designmatrix, G with new toas N = np.diag(self.toaerrs**2) if self.A_rn is not None: plaw = red_noise_powerlaw(A=10**self.A_rn, From b450850efdd6ff5ff4a0b264745d9b69cd4fe813 Mon Sep 17 00:00:00 2001 From: Jeremy Baier Date: Wed, 20 Mar 2024 21:40:41 -0700 Subject: [PATCH 15/32] commiting changes --- hasasia/sensitivity.py | 98 ++++++++++++++++++++++++++++++++++++------ 1 file changed, 84 insertions(+), 14 deletions(-) diff --git a/hasasia/sensitivity.py b/hasasia/sensitivity.py index 79e0b72..e3067ac 100644 --- a/hasasia/sensitivity.py +++ b/hasasia/sensitivity.py @@ -364,17 +364,12 @@ def filter_data(self, start_time=None, end_time=None): #dmx_mask = np.sum(self.designmatrix, axis=0) != 0.0 #self.designmatrix = self.designmatrix[:, dmx_mask] self._G = G_matrix(designmatrix=self.designmatrix) - #print("dim designmatrix = ", self.designmatrix.shape) - #print("dim G matrxi = ", self.G.shape) def change_cadence(self, start_time=0, end_time=1_000_000, - cadence=None, cadence_factor=2, uneven=False, - A_gwb=None, gamma_gwb=13/3., freqs=None, + cadence=None, cadence_factor=4, uneven=False, + A_gwb=None, alpha_gwb=-2/3., freqs=None, fast=True,): """ - ***************** - WORK IN PROGRESSS - ***************** Parameters ========== start_time - float @@ -382,11 +377,23 @@ def change_cadence(self, start_time=0, end_time=1_000_000, end_time - float MJD at which to end altared cadence. cadence - float - factor by which to change the observing cadence. - + cadence for the modified campaign [toas/year] + cadence_facter - float + (instead of cadence) factor by which to modify the old cadence. + uneven - bool + whether or not to evenly space observation epochs + A_gwb - float + amplitude of injected gwb self-noise + alpha_gwb - float + spectral index of injected gwb self-noise. + note that this is residual space spectral index. + freqs - array + frequencies to construct the gwb noise and intrinsic noise + fast - bool + faster but slightly less accurate method to calculate noise injected in N. Change observing cadence in a given time range. - Recalculate pulsar properties. + Recalculate pulsar noise properties. """ mask_before = self.toas <= start_time * 86400 mask_after = self.toas >= end_time * 86400 @@ -418,19 +425,82 @@ def change_cadence(self, start_time=0, end_time=1_000_000, dt = duration / campaign_Ntoas / 8 * yr_sec campaign_toas += np.random.uniform(-dt, dt, size=campaign_Ntoas) self.toas = np.concatenate([old_toas[mask_before], campaign_toas, old_toas[mask_after]]) - campaign_toaerrs = np.mean(old_toaerrs)*np.ones(campaign_Ntoas) - self.toaerrs = np.concatenate([old_toaerrs[mask_before], campaign_toaerrs, old_toaerrs[mask_after]]) + campaign_toaerrs = np.median(old_toaerrs)*np.ones(campaign_Ntoas) + # FIXME can only use a fixed toaerr for the duration of the campaign + #self.toaerrs = np.concatenate([old_toaerrs[mask_before], campaign_toaerrs, old_toaerrs[mask_after]]) + self.toaerrs = np.ones(len(self.toas))*old_toaerrs[0] + print(f"old: {len(old_toaerrs)}, new: {len(self.toaerrs)}") # recalculate N, designmatrix, G with new toas N = np.diag(self.toaerrs**2) if self.A_rn is not None: - plaw = red_noise_powerlaw(A=10**self.A_rn, + plaw = red_noise_powerlaw(A=self.A_rn, alpha=self.alpha, freqs=freqs) N += corr_from_psd(freqs=freqs, psd=plaw, toas=self.toas, fast=fast) if A_gwb is not None: gwb = red_noise_powerlaw(A=A_gwb, - alpha=gamma_gwb, + alpha=alpha_gwb, + freqs=freqs) + N += corr_from_psd(freqs=freqs, psd=gwb, toas=self.toas, fast=fast) + self.designmatrix = create_design_matrix(self.toas, RADEC=True, PROPER=True, PX=True) + self._G = G_matrix(designmatrix=self.designmatrix) + self.N = N + + def change_sigma(self, start_time=0, end_time=1_000_000, + new_sigma=None, sigma_factor=4, uneven=False, + A_gwb=None, alpha_gwb=-2/3., freqs=None, + fast=True,): + """ + Parameters + ========== + start_time - float + MJD at which to begin altared cadence. + end_time - float + MJD at which to end altared cadence. + new_sigma - float + unvertainty of toas for the modified campaign [microseconds] + sigma_facter - float + (instead of sigmas) factor by which to modify the campaign sigmas. + uneven - bool + whether or not to evenly space observation epochs + A_gwb - float + amplitude of injected gwb self-noise + alpha_gwb - float + spectral index of injected gwb self-noise. + note that this is residual space spectral index (alpha). + freqs - array + frequencies to construct the gwb noise and intrinsic noise + fast - bool + faster but slightly less accurate method to calculate noise injected in N. + + Change observing cadence in a given time range. + Recalculate pulsar noise properties. + """ + mask_before = self.toas <= start_time * 86400 + mask_after = self.toas >= end_time * 86400 + campaign_mask = np.logical_and(self.toas >= start_time * 86400, + self.toas <= end_time * 86400) + campaign_ntoas = np.sum(campaign_mask) + # store the old toa errors + toaerrs_campaign = self.toaerrs[campaign_mask] + # modify the campaign toa errors + if sigma_factor is not None: + toaerrs_campaign = sigma_factor * toaerrs_campaign + elif sigma_factor is None and new_sigma is not None: + toaerrs_campaign = np.ones(campaign_ntoas)*new_sigma + self.toaerrs = np.concatenate([self.toaerrs[mask_before], toaerrs_campaign, self.toaerrs[mask_after]]) + # recalculate N, designmatrix, G with new toas + N = np.diag(self.toaerrs**2) + if self.A_rn is not None: + plaw = red_noise_powerlaw(A=self.A_rn, + alpha=self.alpha, + freqs=freqs) + N += corr_from_psd(freqs=freqs, psd=plaw, toas=self.toas, fast=fast) + + if A_gwb is not None: + gwb = red_noise_powerlaw(A=A_gwb, + alpha=alpha_gwb, freqs=freqs) N += corr_from_psd(freqs=freqs, psd=gwb, toas=self.toas, fast=fast) self.designmatrix = create_design_matrix(self.toas, RADEC=True, PROPER=True, PX=True) From 48e7358205e9c05c189ab1bfd234ceca7b35643b Mon Sep 17 00:00:00 2001 From: Jeremy Baier Date: Thu, 25 Jul 2024 22:42:02 -0700 Subject: [PATCH 16/32] enabling a polarization averaging with explicit inclination --- hasasia/skymap.py | 72 +++++++++++++++++++++++++++++------------------ 1 file changed, 44 insertions(+), 28 deletions(-) diff --git a/hasasia/skymap.py b/hasasia/skymap.py index 8b42e00..8585868 100644 --- a/hasasia/skymap.py +++ b/hasasia/skymap.py @@ -144,7 +144,7 @@ def h_thresh(self, SNR=1, iota=None, psi=None): An array representing the skymap of amplitudes needed to see the given signal with the SNR threshold specified. ''' - if iota is None or psi is None: + if iota is None: S_eff = self.S_eff else: S_eff = self.S_eff_full(iota, psi) @@ -174,11 +174,11 @@ def A_gwb(self, h_div_A, SNR=1): integrand = h_div_A[:,np.newaxis]**2 / self.S_eff return SNR / np.sqrt(np.trapz(integrand,x=self.freqs,axis=0 )) - def sky_response_full(self, iota, psi): + def sky_response_full(self, iota, psi=None): r''' Calculate the signal-to-noise ratio of a source given the strain amplitude. This is based on Equation (79) from Hazboun, et al., 2019 - `[1]`_. + `[1]`_, but modified so that you calculate it for a particular inclination angle. .. math:: \rho(\hat{n})=h_0\sqrt{\frac{T_{\rm obs}}{S_{\rm eff}(f_0 ,\hat{k})}} @@ -187,28 +187,36 @@ def sky_response_full(self, iota, psi): ''' Ap_sqr = (0.5 * (1 + np.cos(iota)**2))**2 Ac_sqr = (np.cos(iota))**2 - iota = iota if isinstance(iota, (int,float)) else np.array(iota) - psi = psi if isinstance(psi, (int,float)) else np.array(psi) - spsi = np.sin(2*np.array(psi)) - cpsi = np.cos(2*np.array(psi)) - if isinstance(Ap_sqr,np.ndarray) or isinstance(spsi,np.ndarray): - c1 = Ac_sqr[:,np.newaxis]*spsi**2 + Ap_sqr[:,np.newaxis]*cpsi**2 - c2 = (Ap_sqr[:,np.newaxis] - Ac_sqr[:,np.newaxis]) * cpsi * spsi - c3 = Ap_sqr[:,np.newaxis]*spsi**2 + Ac_sqr[:,np.newaxis]*cpsi**2 - term1 = self.Fplus[:,:,np.newaxis,np.newaxis]**2 * c1 - term2 = 2 * self.Fplus[:,:,np.newaxis,np.newaxis] * self.Fcross[:,:,np.newaxis,np.newaxis] * c2 - term3 = self.Fcross[:,:,np.newaxis,np.newaxis]**2 * c3 - elif isinstance(Ap_sqr,(int,float)) or isinstance(spsi,(int,float)): - c1 = Ac_sqr*spsi**2 + Ap_sqr*cpsi**2 - c2 = (Ap_sqr - Ac_sqr) * cpsi * spsi - c3 = Ap_sqr*spsi**2 + Ac_sqr*cpsi**2 - term1 = self.Fplus**2 * c1 - term2 = 2 * self.Fplus * self.Fcross * c2 - term3 = self.Fcross**2 * c3 - - return term1 + term2 + term3 - - def S_SkyI_full(self, iota, psi): + if psi is None: # case where we average over polarization but not inclination + # 0.5 is from averaging over polarization + # Fplus*Fcross term goes to zero + if isinstance(Ap_sqr,np.ndarray): + return (self.Fplus[:,:,np.newaxis]**2 + self.Fcross[:,:,np.newaxis]**2)* 0.5 * (Ap_sqr + Ac_sqr) + elif isinstance(Ap_sqr,(int,float)): + return (self.Fplus**2 + self.Fcross**2) * 0.5 * (Ac_sqr + Ap_sqr) + else: # case where we don't average over polarization or inclination + iota = iota if isinstance(iota, (int,float)) else np.array(iota) + psi = psi if isinstance(psi, (int,float)) else np.array(psi) + spsi = np.sin(2*np.array(psi)) + cpsi = np.cos(2*np.array(psi)) + if isinstance(Ap_sqr,np.ndarray) or isinstance(spsi,np.ndarray): + c1 = Ac_sqr[:,np.newaxis]*spsi**2 + Ap_sqr[:,np.newaxis]*cpsi**2 + c2 = (Ap_sqr[:,np.newaxis] - Ac_sqr[:,np.newaxis]) * cpsi * spsi + c3 = Ap_sqr[:,np.newaxis]*spsi**2 + Ac_sqr[:,np.newaxis]*cpsi**2 + term1 = self.Fplus[:,:,np.newaxis,np.newaxis]**2 * c1 + term2 = 2 * self.Fplus[:,:,np.newaxis,np.newaxis] * self.Fcross[:,:,np.newaxis,np.newaxis] * c2 + term3 = self.Fcross[:,:,np.newaxis,np.newaxis]**2 * c3 + elif isinstance(Ap_sqr,(int,float)) or isinstance(spsi,(int,float)): + c1 = Ac_sqr*spsi**2 + Ap_sqr*cpsi**2 + c2 = (Ap_sqr - Ac_sqr) * cpsi * spsi + c3 = Ap_sqr*spsi**2 + Ac_sqr*cpsi**2 + term1 = self.Fplus**2 * c1 + term2 = 2 * self.Fplus * self.Fcross * c2 + term3 = self.Fcross**2 * c3 + + return term1 + term2 + term3 + + def S_SkyI_full(self, iota, psi=None): """Per Pulsar Strain power sensitivity. """ t_I = self.T_I / self.Tspan RNcalInv = 3.0 * t_I[:,np.newaxis] / self.SnI @@ -221,14 +229,22 @@ def S_SkyI_full(self, iota, psi): if sky_resp.ndim == 2: self._S_SkyI_full = (RNcalInv.T[:,:,np.newaxis] * sky_resp[np.newaxis,:,:]) - if sky_resp.ndim == 4: + elif sky_resp.ndim == 3: + self._S_SkyI_full = (RNcalInv.T[:,:,np.newaxis,np.newaxis] + * sky_resp[np.newaxis,:,:,:]) + elif sky_resp.ndim == 4: self._S_SkyI_full = (RNcalInv.T[:,:,np.newaxis,np.newaxis,np.newaxis] * sky_resp[np.newaxis,:,:,:,:]) return self._S_SkyI_full - def S_eff_full(self, iota, psi): - """Strain power sensitivity. """ + def S_eff_full(self, iota, psi=None): + """ + Strain power sensitivity. + + Can calculate margninalized over polarization or not + with inclination explicit in both cases. + """ if self.pulsar_term == 'explicit': raise NotImplementedError('Currently cannot use pulsar term.') # self._S_eff_full = 1.0 / (4./5 * np.sum(self.S_SkyI, axis=1)) From 4a1ec8c8c927862c9229865f538fa35256de7a68 Mon Sep 17 00:00:00 2001 From: Jeremy Baier Date: Sun, 4 Aug 2024 21:50:56 -0700 Subject: [PATCH 17/32] adding detector volume calculate as an method of skymap and adding convenience function for D_L calculation --- hasasia/skymap.py | 35 +++++++++++++++++++++++++++++++++++ hasasia/utils.py | 23 +++++++++++++++++++++++ 2 files changed, 58 insertions(+) diff --git a/hasasia/skymap.py b/hasasia/skymap.py index 8585868..a7c732b 100644 --- a/hasasia/skymap.py +++ b/hasasia/skymap.py @@ -3,9 +3,11 @@ """Main module.""" import numpy as np import scipy.special as spec +import healpy as hp import astropy.units as u import astropy.constants as c from .sensitivity import DeterSensitivityCurve, resid_response, get_dt +from .utils import strain_and_chirp_mass_to_luminosity_distance __all__ = ['SkySensitivity', 'h_circ', @@ -307,6 +309,39 @@ def S_eff_mean(self): return self._S_eff_mean +def calculate_detection_volume(skymap, frequency, SNR_threshold, M_c): + """ + Calculates the detection volume of your PTA + at a given frequency or list of frequencies. + + Parameters + ========== + skymap - hasasia.skymap + the hasasia.skymap to use + frequency - float + the frequency [Hz] at which to calculate detection volume + SNR_threshold - float + the signal to noise to referene detection volume to + M_c - float + the chirp mass [Msun] at which to reference detection volume + Returns + ======= + volume - float + the detection volume in Mpc^3 + """ + NSIDE = hp.pixelfunc.npix2nside(skymap.S_eff.shape[1]) + dA = hp.pixelfunc.nside2pixarea(NSIDE, degrees=False) + if isinstance(frequency, (int,float)): + f_idx = np.array([np.argmin(abs(skymap.freqs - frequency))]) + elif isinstance(frequency, (np.ndarray, list)): + f_idx = np.array([np.argmin(abs(skymap.freqs - f)) for f in frequency]) + h0 = skymap.h_thresh(SNR=SNR_threshold) + volume = [dA*np.sum( + strain_and_chirp_mass_to_luminosity_distance(h0[fdx], M_c, skymap.freqs[fdx])**3, + axis=0).value for fdx in f_idx] + return volume[0] if len(volume)==1 else volume + + def h_circ(M_c, D_L, f0, Tspan, f): r""" Convenience function that returns the Fourier domain representation of a diff --git a/hasasia/utils.py b/hasasia/utils.py index 1687f3c..c23e908 100644 --- a/hasasia/utils.py +++ b/hasasia/utils.py @@ -7,6 +7,7 @@ import scipy.stats as ss from astropy.coordinates import SkyCoord import astropy.units as u +import astropy.constants as c __all__ = ['create_design_matrix', 'fap', @@ -131,6 +132,28 @@ def _solve_F0_given_SNR(snr=3, Npsrs=None): Npsrs = 1 if Npsrs is None else Npsrs return 0.5*(4.*Npsrs+snr**2.) +def strain_and_chirp_mass_to_luminosity_distance(h, M_c, f0): + r''' + Returns the luminosity distance to a source given the strain, chirp mass, and GW frequency. + + Parameters + ---------- + h : float + The strain of the source. + M_c : float + The chirp mass of the source [Msun]. + f0 : float + The GW frequency of the source [Hz]. + + Returns + ------- + d_L : float + The luminosity distance to the source [Mpc]. + ''' + return (4*c.c / (h * u.m/u.m) + * np.power(c.G * M_c * u.Msun/c.c**3, 5/3) + * np.power(np.pi * f0 * u.Hz, 2/3)).to('Mpc') + def theta_phi_to_SkyCoord(theta, phi): """ Parameters From 97cf01834b14c060158686bb2a7d0ec930a264d2 Mon Sep 17 00:00:00 2001 From: Jeremy Baier Date: Sun, 4 Aug 2024 21:51:26 -0700 Subject: [PATCH 18/32] updating the last commit --- hasasia/skymap.py | 44 +++++++++++++++++++++++++------------------- hasasia/utils.py | 2 +- 2 files changed, 26 insertions(+), 20 deletions(-) diff --git a/hasasia/skymap.py b/hasasia/skymap.py index a7c732b..63c82cc 100644 --- a/hasasia/skymap.py +++ b/hasasia/skymap.py @@ -273,7 +273,8 @@ def S_eff(self): @property def S_SkyI(self): - """Per Pulsar Strain power sensitivity. """ + """Per Pulsar Strain power sensitivity. + (Technically, 1 over this is the per pulsar strain power sensitivity.)""" if not hasattr(self, '_S_SkyI'): t_I = self.T_I / self.Tspan RNcalInv = t_I[:,np.newaxis] / self.SnI @@ -309,35 +310,40 @@ def S_eff_mean(self): return self._S_eff_mean -def calculate_detection_volume(skymap, frequency, SNR_threshold, M_c): - """ - Calculates the detection volume of your PTA +def calculate_detection_volume(self, f0, SNR_threshold=3.7145, M_c=1e9): + r""" + Calculates the detection volume of the PTA at a given frequency or list of frequencies. Parameters - ========== - skymap - hasasia.skymap - the hasasia.skymap to use - frequency - float + ---------- + + f0 : float the frequency [Hz] at which to calculate detection volume - SNR_threshold - float + + SNR_threshold : float the signal to noise to referene detection volume to - M_c - float + + M_c : float the chirp mass [Msun] at which to reference detection volume + Returns - ======= - volume - float + ------- + + volume : float the detection volume in Mpc^3 + """ - NSIDE = hp.pixelfunc.npix2nside(skymap.S_eff.shape[1]) + NSIDE = hp.pixelfunc.npix2nside(self.S_eff.shape[1]) dA = hp.pixelfunc.nside2pixarea(NSIDE, degrees=False) - if isinstance(frequency, (int,float)): - f_idx = np.array([np.argmin(abs(skymap.freqs - frequency))]) - elif isinstance(frequency, (np.ndarray, list)): - f_idx = np.array([np.argmin(abs(skymap.freqs - f)) for f in frequency]) - h0 = skymap.h_thresh(SNR=SNR_threshold) + if isinstance(f0, (int,float)): + f_idx = np.array([np.argmin(abs(self.freqs - f0))]) + elif isinstance(f0, (np.ndarray, list)): + f_idx = np.array([np.argmin(abs(self.freqs - f)) for f in f0]) + h0 = self.h_thresh(SNR=SNR_threshold) + # detection volume is is the sum of detection radius * pixel area over all pixels volume = [dA*np.sum( - strain_and_chirp_mass_to_luminosity_distance(h0[fdx], M_c, skymap.freqs[fdx])**3, + strain_and_chirp_mass_to_luminosity_distance(h0[fdx], M_c, self.freqs[fdx])**3, axis=0).value for fdx in f_idx] return volume[0] if len(volume)==1 else volume diff --git a/hasasia/utils.py b/hasasia/utils.py index c23e908..2e82942 100644 --- a/hasasia/utils.py +++ b/hasasia/utils.py @@ -147,7 +147,7 @@ def strain_and_chirp_mass_to_luminosity_distance(h, M_c, f0): Returns ------- - d_L : float + D_L : float The luminosity distance to the source [Mpc]. ''' return (4*c.c / (h * u.m/u.m) From c7c8531c9f3b6ebd0e21a44248242d95c0d1f09e Mon Sep 17 00:00:00 2001 From: Jeremy Baier Date: Sun, 4 Aug 2024 21:54:23 -0700 Subject: [PATCH 19/32] updating formatting --- hasasia/utils.py | 22 ++++++++++++---------- 1 file changed, 12 insertions(+), 10 deletions(-) diff --git a/hasasia/utils.py b/hasasia/utils.py index 2e82942..80c75d9 100644 --- a/hasasia/utils.py +++ b/hasasia/utils.py @@ -155,40 +155,42 @@ def strain_and_chirp_mass_to_luminosity_distance(h, M_c, f0): * np.power(np.pi * f0 * u.Hz, 2/3)).to('Mpc') def theta_phi_to_SkyCoord(theta, phi): - """ + r''' + Takes in a celestial longitude and lattitude and returns an `astropy.SkyCoord` object. + Parameters - ========== + ---------- phi : float, array of floats The celestial longitude in solar system coordinates. theta : float, array of floats The celestial lattitude in solar system coordinates. Returns - ======= + ------- skycoord - astropy.SkyCoord object Can use this to convert to ra, dec, etc. (e.g. SkyCoord.ra.deg) - Converts an inputed longitude and lattitude into an `astropy.SkyCoord` object. - - """ + ''' return SkyCoord(phi*u.rad, ( theta - np.pi/2 )*u.rad) def skycoord_to_Jname(skycoord): - """ + ''' + Takes in a SkyCoord object and returns the Jname of a pulsar with given coordinates. + Parameters - ========== + ---------- skycoord - astropy.SkyCoord object Can use `theta_phi_to_SkyCoord()` to get this. Returns - ======= + ------- Jname - string, array of strings The traditional Jname of a pulsar with given coordinates (eg. 'J1713+0747') - """ + ''' coord_pieces = [ str(skycoord.ra.hms[0]).split('.')[0], str(skycoord.ra.hms[1]).split('.')[0], From 69c6bbe4c1502b214e27b0a06fa1941eec589ed4 Mon Sep 17 00:00:00 2001 From: Jeremy Baier Date: Mon, 5 Aug 2024 16:02:50 -0700 Subject: [PATCH 20/32] adding not implemented error --- hasasia/skymap.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/hasasia/skymap.py b/hasasia/skymap.py index 63c82cc..78d6232 100644 --- a/hasasia/skymap.py +++ b/hasasia/skymap.py @@ -116,8 +116,10 @@ def SNR(self, h0, iota=None, psi=None): .. _[1]: https://arxiv.org/abs/1907.04341 ''' - if iota is None or psi is None: + if iota is None and psi is None: S_eff = self.S_eff + elif psi is not None and iota is None: + raise NotImplementedError('Currently cannot marginalize over phase but not inclination.') else: S_eff = self.S_eff_full(iota, psi) return h0 * np.sqrt(self.Tspan / S_eff) From f147c90143e011fd6516daa415b1bcf7f3535c77 Mon Sep 17 00:00:00 2001 From: Jeremy Baier Date: Tue, 6 Aug 2024 09:47:42 -0700 Subject: [PATCH 21/32] adding a factor of 2 which is wrong to see what happens --- hasasia/skymap.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/hasasia/skymap.py b/hasasia/skymap.py index 78d6232..f3b34af 100644 --- a/hasasia/skymap.py +++ b/hasasia/skymap.py @@ -189,6 +189,8 @@ def sky_response_full(self, iota, psi=None): .. _[1]: https://arxiv.org/abs/1907.04341 ''' + if iota is None and psi is not None: + raise NotImplementedError('Currently cannot marginalize over inclination but not phase.') Ap_sqr = (0.5 * (1 + np.cos(iota)**2))**2 Ac_sqr = (np.cos(iota))**2 if psi is None: # case where we average over polarization but not inclination @@ -197,7 +199,7 @@ def sky_response_full(self, iota, psi=None): if isinstance(Ap_sqr,np.ndarray): return (self.Fplus[:,:,np.newaxis]**2 + self.Fcross[:,:,np.newaxis]**2)* 0.5 * (Ap_sqr + Ac_sqr) elif isinstance(Ap_sqr,(int,float)): - return (self.Fplus**2 + self.Fcross**2) * 0.5 * (Ac_sqr + Ap_sqr) + return (self.Fplus**2 + self.Fcross**2) * 0.25 * (Ac_sqr + Ap_sqr) else: # case where we don't average over polarization or inclination iota = iota if isinstance(iota, (int,float)) else np.array(iota) psi = psi if isinstance(psi, (int,float)) else np.array(psi) From 2bca1b4ffeeb9a91bc6fd114263346794d5a03ce Mon Sep 17 00:00:00 2001 From: Jeremy Baier Date: Tue, 6 Aug 2024 10:07:39 -0700 Subject: [PATCH 22/32] rolling it back --- hasasia/skymap.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/hasasia/skymap.py b/hasasia/skymap.py index f3b34af..757fd67 100644 --- a/hasasia/skymap.py +++ b/hasasia/skymap.py @@ -199,7 +199,7 @@ def sky_response_full(self, iota, psi=None): if isinstance(Ap_sqr,np.ndarray): return (self.Fplus[:,:,np.newaxis]**2 + self.Fcross[:,:,np.newaxis]**2)* 0.5 * (Ap_sqr + Ac_sqr) elif isinstance(Ap_sqr,(int,float)): - return (self.Fplus**2 + self.Fcross**2) * 0.25 * (Ac_sqr + Ap_sqr) + return (self.Fplus**2 + self.Fcross**2) * 0.5 * (Ac_sqr + Ap_sqr) else: # case where we don't average over polarization or inclination iota = iota if isinstance(iota, (int,float)) else np.array(iota) psi = psi if isinstance(psi, (int,float)) else np.array(psi) From 9e8bf7026a8fd708718d340cc94c0f130f04457d Mon Sep 17 00:00:00 2001 From: Jeremy Baier Date: Tue, 6 Aug 2024 10:24:05 -0700 Subject: [PATCH 23/32] new factor --- hasasia/skymap.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/hasasia/skymap.py b/hasasia/skymap.py index 757fd67..4e96bce 100644 --- a/hasasia/skymap.py +++ b/hasasia/skymap.py @@ -197,9 +197,9 @@ def sky_response_full(self, iota, psi=None): # 0.5 is from averaging over polarization # Fplus*Fcross term goes to zero if isinstance(Ap_sqr,np.ndarray): - return (self.Fplus[:,:,np.newaxis]**2 + self.Fcross[:,:,np.newaxis]**2)* 0.5 * (Ap_sqr + Ac_sqr) + return (self.Fplus[:,:,np.newaxis]**2 + self.Fcross[:,:,np.newaxis]**2)* 4/5 * (Ap_sqr + Ac_sqr) elif isinstance(Ap_sqr,(int,float)): - return (self.Fplus**2 + self.Fcross**2) * 0.5 * (Ac_sqr + Ap_sqr) + return (self.Fplus**2 + self.Fcross**2) * 4/5 * (Ac_sqr + Ap_sqr) else: # case where we don't average over polarization or inclination iota = iota if isinstance(iota, (int,float)) else np.array(iota) psi = psi if isinstance(psi, (int,float)) else np.array(psi) From 52f20817634a62cb4436c14fc1ad21eec8bccf4b Mon Sep 17 00:00:00 2001 From: Jeremy Baier Date: Tue, 6 Aug 2024 10:30:17 -0700 Subject: [PATCH 24/32] removing the factor to see what it is --- hasasia/skymap.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/hasasia/skymap.py b/hasasia/skymap.py index 4e96bce..bc5d5d3 100644 --- a/hasasia/skymap.py +++ b/hasasia/skymap.py @@ -197,9 +197,9 @@ def sky_response_full(self, iota, psi=None): # 0.5 is from averaging over polarization # Fplus*Fcross term goes to zero if isinstance(Ap_sqr,np.ndarray): - return (self.Fplus[:,:,np.newaxis]**2 + self.Fcross[:,:,np.newaxis]**2)* 4/5 * (Ap_sqr + Ac_sqr) + return (self.Fplus[:,:,np.newaxis]**2 + self.Fcross[:,:,np.newaxis]**2) * (Ap_sqr + Ac_sqr) elif isinstance(Ap_sqr,(int,float)): - return (self.Fplus**2 + self.Fcross**2) * 4/5 * (Ac_sqr + Ap_sqr) + return (self.Fplus**2 + self.Fcross**2) * (Ac_sqr + Ap_sqr) else: # case where we don't average over polarization or inclination iota = iota if isinstance(iota, (int,float)) else np.array(iota) psi = psi if isinstance(psi, (int,float)) else np.array(psi) From 80e6178e44a8016b5801766526696c74fff91917 Mon Sep 17 00:00:00 2001 From: Jeremy Baier Date: Tue, 6 Aug 2024 10:34:28 -0700 Subject: [PATCH 25/32] trying a factor of 5/4 which i found emprically --- hasasia/skymap.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/hasasia/skymap.py b/hasasia/skymap.py index bc5d5d3..7f6d3ff 100644 --- a/hasasia/skymap.py +++ b/hasasia/skymap.py @@ -197,9 +197,9 @@ def sky_response_full(self, iota, psi=None): # 0.5 is from averaging over polarization # Fplus*Fcross term goes to zero if isinstance(Ap_sqr,np.ndarray): - return (self.Fplus[:,:,np.newaxis]**2 + self.Fcross[:,:,np.newaxis]**2) * (Ap_sqr + Ac_sqr) + return (self.Fplus[:,:,np.newaxis]**2 + self.Fcross[:,:,np.newaxis]**2) * 5/4 * (Ap_sqr + Ac_sqr) elif isinstance(Ap_sqr,(int,float)): - return (self.Fplus**2 + self.Fcross**2) * (Ac_sqr + Ap_sqr) + return (self.Fplus**2 + self.Fcross**2) * 5/4 * (Ac_sqr + Ap_sqr) else: # case where we don't average over polarization or inclination iota = iota if isinstance(iota, (int,float)) else np.array(iota) psi = psi if isinstance(psi, (int,float)) else np.array(psi) From 68b77f79fb36943e521e9833824d462b53562e4d Mon Sep 17 00:00:00 2001 From: Jeremy Baier Date: Tue, 6 Aug 2024 14:46:38 -0700 Subject: [PATCH 26/32] switching back to 1/2 --- hasasia/skymap.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/hasasia/skymap.py b/hasasia/skymap.py index 7f6d3ff..d8ab87f 100644 --- a/hasasia/skymap.py +++ b/hasasia/skymap.py @@ -197,9 +197,9 @@ def sky_response_full(self, iota, psi=None): # 0.5 is from averaging over polarization # Fplus*Fcross term goes to zero if isinstance(Ap_sqr,np.ndarray): - return (self.Fplus[:,:,np.newaxis]**2 + self.Fcross[:,:,np.newaxis]**2) * 5/4 * (Ap_sqr + Ac_sqr) + return (self.Fplus[:,:,np.newaxis]**2 + self.Fcross[:,:,np.newaxis]**2) * 0.5 * (Ap_sqr + Ac_sqr) elif isinstance(Ap_sqr,(int,float)): - return (self.Fplus**2 + self.Fcross**2) * 5/4 * (Ac_sqr + Ap_sqr) + return (self.Fplus**2 + self.Fcross**2) * 0.5 * (Ac_sqr + Ap_sqr) else: # case where we don't average over polarization or inclination iota = iota if isinstance(iota, (int,float)) else np.array(iota) psi = psi if isinstance(psi, (int,float)) else np.array(psi) @@ -222,7 +222,7 @@ def sky_response_full(self, iota, psi=None): return term1 + term2 + term3 - def S_SkyI_full(self, iota, psi=None): + def S_SkyI_full(self, iota, psi): """Per Pulsar Strain power sensitivity. """ t_I = self.T_I / self.Tspan RNcalInv = 3.0 * t_I[:,np.newaxis] / self.SnI @@ -244,7 +244,7 @@ def S_SkyI_full(self, iota, psi=None): return self._S_SkyI_full - def S_eff_full(self, iota, psi=None): + def S_eff_full(self, iota, psi): """ Strain power sensitivity. From 30d2453d0635ccdd4a04be31e45c12a1c547a6a6 Mon Sep 17 00:00:00 2001 From: Jeremy Baier Date: Wed, 7 Aug 2024 11:28:50 -0700 Subject: [PATCH 27/32] correcting gamma_gwb to alpha_gwb in sim_pta --- hasasia/sim.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/hasasia/sim.py b/hasasia/sim.py index eae59fe..a6bb561 100644 --- a/hasasia/sim.py +++ b/hasasia/sim.py @@ -12,7 +12,7 @@ def sim_pta(timespan, cad, sigma, phi, theta, Npsrs=None, A_rn=None, alpha=None, freqs=None, uneven=False, - A_gwb=None, gamma_gwb = 13/3, fast=True, psr_names = None, + A_gwb=None, alpha_gwb = -2/3, fast=True, psr_names = None, kwastro={'RADEC':True, 'PROPER':True, 'PX':True}): """ Make a simulated pulsar timing array. Using the available parameters, @@ -123,7 +123,7 @@ def sim_pta(timespan, cad, sigma, phi, theta, Npsrs=None, if A_gwb is not None: gwb = red_noise_powerlaw(A=A_gwb, - alpha=gamma_gwb, + alpha=alpha_gwb, freqs=freqs) N += corr_from_psd(freqs=freqs, psd=gwb, toas=toas, fast=fast) From 9dafb9b1daf0f220f1593e03575e69fb2e892f0d Mon Sep 17 00:00:00 2001 From: Jeremy Baier Date: Sat, 14 Sep 2024 15:27:52 -0700 Subject: [PATCH 28/32] bug fix: fixed missing in detector volume function --- hasasia/skymap.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/hasasia/skymap.py b/hasasia/skymap.py index d8ab87f..74fb39c 100644 --- a/hasasia/skymap.py +++ b/hasasia/skymap.py @@ -346,7 +346,7 @@ def calculate_detection_volume(self, f0, SNR_threshold=3.7145, M_c=1e9): f_idx = np.array([np.argmin(abs(self.freqs - f)) for f in f0]) h0 = self.h_thresh(SNR=SNR_threshold) # detection volume is is the sum of detection radius * pixel area over all pixels - volume = [dA*np.sum( + volume = [dA/3.*np.sum( strain_and_chirp_mass_to_luminosity_distance(h0[fdx], M_c, self.freqs[fdx])**3, axis=0).value for fdx in f_idx] return volume[0] if len(volume)==1 else volume From c389d16c015f2a743897b25e4503ee7cc634d1a8 Mon Sep 17 00:00:00 2001 From: Jeremy Baier Date: Sat, 14 Sep 2024 22:26:38 -0700 Subject: [PATCH 29/32] updates: moved the GWB to GWBSC rather than detSC --- hasasia/sensitivity.py | 66 ++++++++++++++++++++++++++++++++++++++++++ hasasia/skymap.py | 30 +++++++++++++++++++ 2 files changed, 96 insertions(+) diff --git a/hasasia/sensitivity.py b/hasasia/sensitivity.py index e3067ac..81d3b97 100644 --- a/hasasia/sensitivity.py +++ b/hasasia/sensitivity.py @@ -845,6 +845,72 @@ def S_effIJ(self): return self._S_effIJ + def dmu_dbeta_pl_gwb(self, A_gwb, alpha_gwb, components=14): + """Derivative of the GWB signal template (mu) w.r.t. the signal parameters (beta). + $$\mu_{\rm GWB} = \frac{A_{\rm GWB}}{12\pi^2} * \left(\frac{f}{f_{\rm ref}}\right)^{-2\alpha} \rm yr^3""" + # FIXME : how to deal with the frequency dependence of these derivatives ?? + fbins = np.linspace(1/self.Tspan, components/self.Tspan, components) + # dmu_dA = A_gwb/(12*np.pi**2)*(fbins/fyr)**(-2*alpha_gwb) * \ + # yr_sec**3 * np.log(10) + # dmu_dalpha = A_gwb/(12*np.pi**2)*(fbins/fyr)**(-2*alpha_gwb) * \ + # yr_sec**3 * np.log(fbins/fyr) + # As S_h=A^2*(f/fyr)^(2*alpha)/f + dmu_dA = A_gwb*(fbins/fyr)**(alpha_gwb) * np.log(10) + dmu_dalpha = A_gwb*(fbins/fyr)**(alpha_gwb) * np.log(fbins/fyr) + + return np.array([dmu_dA, dmu_dalpha]) + + def dmu_dbeta_freespec_gwb(self, log10_gw_rhos): + """Derivative of the GWB signal template (mu) w.r.t. the signal parameters (beta). + $$\mu_{\rm GWB} = \rho_{i}^2""" + dmu_drhos = [] + for i, _ in enumerate(log10_gw_rhos): + dmu_drho = np.zeros(len(log10_gw_rhos)) + # https://www.wolframalpha.com/input?i=derivative+of+10%5Ex%5E2+wrt+x + dmu_drho[i]+=10**log10_gw_rhos**2 * np.log(10) * 2*log10_gw_rhos + dmu_drhos.append(dmu_drho) + + return np.array(dmu_drhos) + + def GWBFisherMatrix(self, + Agwb, alpha_gwb, + log10_gw_rhos, + psd_type='powerlaw', + components=14): + """GWB Fisher Matrix""" + if psd_type == 'powerlaw': + assert(Agwb is not None) + assert(alpha_gwb is not None) + dmu_dbeta = self.dmu_dbeta_pl_gwb(Agwb, alpha_gwb, components=components) + elif psd_type == 'freespec': + assert(log10_gw_rhos is not None) + dmu_dbeta = self.dmu_dbeta_freespec_gwb(log10_gw_rhos) + else: + raise NotImplementedError(f"{psd_type} PSD not supported yet! \ + Choose \'powerlaw\' or \'freespec\'.") + fbins=np.linspace(1/self.Tspan, components/self.Tspan, components) + fidxs = np.array([np.argmin(abs(self.freqs - f)) for f in fbins]) + # FIXME: check the order of opperations here + mathcalF = np.matmul(dmu_dbeta.T, np.linalg.inv(self.S_eff[fidxs][None, :]), dmu_dbeta) + + return mathcalF + + def GWBFisherUncertainty(self, + Agwb=None, alpha_gwb=None, + log10_gw_rhos=None, + psd_type='powerlaw', + components=14): + """GWB Fisher Matrix Uncertainty""" + + mathcalF = self.GWBFisherMatrix(Agwb, alpha_gwb, + log10_gw_rhos, + psd_type=psd_type, + components=components) + print(mathcalF.shape) + mathcalFinv = np.linalg.inv(mathcalF) + + return np.sqrt(np.diag(mathcalFinv)) + class DeterSensitivityCurve(SensitivityCurve): ''' diff --git a/hasasia/skymap.py b/hasasia/skymap.py index 74fb39c..0e26444 100644 --- a/hasasia/skymap.py +++ b/hasasia/skymap.py @@ -313,6 +313,36 @@ def S_eff_mean(self): self._S_eff_mean = 1.0 / (12./5 * mean_sky) return self._S_eff_mean + @property + def dmu_dbeta_cw(self, Mc, D_l, f0, theta_cw, phi_cw, h0, + iota=None, psi=None, phi0=None): + """Derivative of the CGW signal template (mu) w.r.t. the signal parameters (beta). + $$\mu_{\rm CW} = h_+ F^+_I(\hat{k}) + h_x F^x_I(\hat{k})""" + hplus = jj + hcross = jj + h0 = h0_circ(Mc, D_l, f0) + dmu_dbeta = [] + + return np.array(dmu_dbeta) + + @property + def CWFisherMatrix(self, Mc, f0, theta_cw, phi_cw, h0, + iota=None, psi=None, phi0=None): + """Continuous GW Fisher Matrix""" + dmu_dbeta = self.dmu_dbeta_cw(Mc, f0, theta_cw, phi_cw, h0, iota=None, psi=None, phi0=None) + mathcalF = dmu_dbeta.T * self.NcalInv * dmu_dbeta + + return mathcalF + + @property + def CWFisherUncertainty(self, Mc, f0, theta_cw, phi_cw, h0, + iota=None, psi=None, phi0=None): + """GWB Fisher Matrix Uncertainty""" + mathcalF = self.CWFisherMatrix(Mc, f0, theta_cw, phi_cw, h0, + iota=None, psi=None, phi0=None) + mathcalFinv = np.linalg.inv(mathcalF) + + return np.sqrt(np.diag(mathcalFinv)) def calculate_detection_volume(self, f0, SNR_threshold=3.7145, M_c=1e9): r""" From 4ca390be2c0618e229d05dee9df5f8a72bcde8ad Mon Sep 17 00:00:00 2001 From: Jeremy Baier Date: Sun, 15 Sep 2024 14:42:48 -0700 Subject: [PATCH 30/32] add: add a few docstring comments --- hasasia/sensitivity.py | 32 +++++++++++++++++++++++++------- 1 file changed, 25 insertions(+), 7 deletions(-) diff --git a/hasasia/sensitivity.py b/hasasia/sensitivity.py index 81d3b97..088b875 100644 --- a/hasasia/sensitivity.py +++ b/hasasia/sensitivity.py @@ -265,6 +265,7 @@ def get_NcalInv(psr, nf=200, fmin=None, fmax=2e-7, freqs=None, NcalInv = np.linalg.inv(Ncal) #N_TOA-N_par x N_TOA-N_par TfN = np.matmul(np.conjugate(Gtilde),np.matmul(NcalInv,Gtilde.T)) / 2 + # N_freq x N_freq if return_Gtilde_Ncal: return np.real(TfN), Gtilde, Ncal elif full_matrix: @@ -848,7 +849,7 @@ def S_effIJ(self): def dmu_dbeta_pl_gwb(self, A_gwb, alpha_gwb, components=14): """Derivative of the GWB signal template (mu) w.r.t. the signal parameters (beta). $$\mu_{\rm GWB} = \frac{A_{\rm GWB}}{12\pi^2} * \left(\frac{f}{f_{\rm ref}}\right)^{-2\alpha} \rm yr^3""" - # FIXME : how to deal with the frequency dependence of these derivatives ?? + # FIXME : what units should this be in ?? fbins = np.linspace(1/self.Tspan, components/self.Tspan, components) # dmu_dA = A_gwb/(12*np.pi**2)*(fbins/fyr)**(-2*alpha_gwb) * \ # yr_sec**3 * np.log(10) @@ -876,8 +877,11 @@ def GWBFisherMatrix(self, Agwb, alpha_gwb, log10_gw_rhos, psd_type='powerlaw', - components=14): - """GWB Fisher Matrix""" + components=14,psr=None): + """ + GWB Fisher Matrix + Fisher matrix = signal derivative x noise covariance x signal derivative + """ if psd_type == 'powerlaw': assert(Agwb is not None) assert(alpha_gwb is not None) @@ -891,7 +895,16 @@ def GWBFisherMatrix(self, fbins=np.linspace(1/self.Tspan, components/self.Tspan, components) fidxs = np.array([np.argmin(abs(self.freqs - f)) for f in fbins]) # FIXME: check the order of opperations here - mathcalF = np.matmul(dmu_dbeta.T, np.linalg.inv(self.S_eff[fidxs][None, :]), dmu_dbeta) + #ncalinv = get_NcalInv(psr, freqs=self.freqs[fidxs]) + print("seff: ", np.diag(self.S_eff[fidxs]).shape) + print("dmu: ", dmu_dbeta.shape) + print("ncalinv: ", psr.NcalInv.shape) + #first = np.matmul(np.diag(psr.NcalInv[fidxs]), dmu_dbeta.T) + first = np.matmul(np.diag(self.S_eff[fidxs]), dmu_dbeta.T) + print(first.shape) + mathcalF = np.matmul(dmu_dbeta, first) + # + # return mathcalF @@ -899,15 +912,20 @@ def GWBFisherUncertainty(self, Agwb=None, alpha_gwb=None, log10_gw_rhos=None, psd_type='powerlaw', - components=14): - """GWB Fisher Matrix Uncertainty""" + components=14,psr=None): + """ + GWB Fisher Matrix Uncertainty + uncertainty of the i'th param is the square root of the the + (i,i) element of the inverse Fisher matrix + """ mathcalF = self.GWBFisherMatrix(Agwb, alpha_gwb, log10_gw_rhos, psd_type=psd_type, - components=components) + components=components,psr=psr) print(mathcalF.shape) mathcalFinv = np.linalg.inv(mathcalF) + #mathcalFinv = mathcalF return np.sqrt(np.diag(mathcalFinv)) From a08a32ab1ab80c95805a745e6ad12003e4c4ee04 Mon Sep 17 00:00:00 2001 From: Jeremy Baier Date: Sun, 15 Sep 2024 22:15:22 -0700 Subject: [PATCH 31/32] pushing something which is messy but i think it works --- hasasia/sensitivity.py | 55 +++++++++++++++++++++++++----------------- 1 file changed, 33 insertions(+), 22 deletions(-) diff --git a/hasasia/sensitivity.py b/hasasia/sensitivity.py index 088b875..5eedcfb 100644 --- a/hasasia/sensitivity.py +++ b/hasasia/sensitivity.py @@ -846,19 +846,21 @@ def S_effIJ(self): return self._S_effIJ - def dmu_dbeta_pl_gwb(self, A_gwb, alpha_gwb, components=14): + def dmu_dbeta_pl_gwb(self, A_gwb, gamma_gwb, components=20): """Derivative of the GWB signal template (mu) w.r.t. the signal parameters (beta). $$\mu_{\rm GWB} = \frac{A_{\rm GWB}}{12\pi^2} * \left(\frac{f}{f_{\rm ref}}\right)^{-2\alpha} \rm yr^3""" # FIXME : what units should this be in ?? fbins = np.linspace(1/self.Tspan, components/self.Tspan, components) - # dmu_dA = A_gwb/(12*np.pi**2)*(fbins/fyr)**(-2*alpha_gwb) * \ + #fbins = self.freqs + # dmu_dA = A_gwb/(12*np.pi**2)*(fbins/fyr)**(-2*gamma_gwb) * \ # yr_sec**3 * np.log(10) - # dmu_dalpha = A_gwb/(12*np.pi**2)*(fbins/fyr)**(-2*alpha_gwb) * \ + # dmu_dalpha = A_gwb/(12*np.pi**2)*(fbins/fyr)**(-2*gamma_gwb) * \ # yr_sec**3 * np.log(fbins/fyr) - # As S_h=A^2*(f/fyr)^(2*alpha)/f - dmu_dA = A_gwb*(fbins/fyr)**(alpha_gwb) * np.log(10) - dmu_dalpha = A_gwb*(fbins/fyr)**(alpha_gwb) * np.log(fbins/fyr) - + dmu_dA = A_gwb/(12*np.pi**2)*(fbins/fyr)**(-1*gamma_gwb) * \ + yr_sec**3 * np.log(10) + dmu_dalpha = A_gwb/(12*np.pi**2)*(fbins/fyr)**(-1*gamma_gwb) * \ + yr_sec**3 * np.log(fbins/fyr) + return np.array([dmu_dA, dmu_dalpha]) def dmu_dbeta_freespec_gwb(self, log10_gw_rhos): @@ -874,18 +876,18 @@ def dmu_dbeta_freespec_gwb(self, log10_gw_rhos): return np.array(dmu_drhos) def GWBFisherMatrix(self, - Agwb, alpha_gwb, + Agwb, gamma_gwb, log10_gw_rhos, psd_type='powerlaw', - components=14,psr=None): + components=20,psr=None): """ GWB Fisher Matrix Fisher matrix = signal derivative x noise covariance x signal derivative """ if psd_type == 'powerlaw': assert(Agwb is not None) - assert(alpha_gwb is not None) - dmu_dbeta = self.dmu_dbeta_pl_gwb(Agwb, alpha_gwb, components=components) + assert(gamma_gwb is not None) + dmu_dbeta = self.dmu_dbeta_pl_gwb(Agwb, gamma_gwb, components=components) elif psd_type == 'freespec': assert(log10_gw_rhos is not None) dmu_dbeta = self.dmu_dbeta_freespec_gwb(log10_gw_rhos) @@ -893,38 +895,47 @@ def GWBFisherMatrix(self, raise NotImplementedError(f"{psd_type} PSD not supported yet! \ Choose \'powerlaw\' or \'freespec\'.") fbins=np.linspace(1/self.Tspan, components/self.Tspan, components) + #fbins=self.freqs fidxs = np.array([np.argmin(abs(self.freqs - f)) for f in fbins]) # FIXME: check the order of opperations here #ncalinv = get_NcalInv(psr, freqs=self.freqs[fidxs]) print("seff: ", np.diag(self.S_eff[fidxs]).shape) print("dmu: ", dmu_dbeta.shape) - print("ncalinv: ", psr.NcalInv.shape) + #print("ncalinv: ", psr.NcalInv.shape) #first = np.matmul(np.diag(psr.NcalInv[fidxs]), dmu_dbeta.T) - first = np.matmul(np.diag(self.S_eff[fidxs]), dmu_dbeta.T) - print(first.shape) - mathcalF = np.matmul(dmu_dbeta, first) - # - # + #first = np.matmul(np.diag(self.S_eff[fidxs]), dmu_dbeta.T) + # sigmainv = np.diag(self.S_eff[fidxs]) + # first = np.matmul(sigmainv, dmu_dbeta.T) + # print("1st matmul", first.shape) + # mathcalF = np.matmul(dmu_dbeta, first) + # print("FishMat", mathcalF.shape) + fishers = np.zeros((2,2)) + for i in range(self.S_effIJ.shape[0]): + fisher = np.matmul(self.dmu_dbeta_pl_gwb(Agwb, + gamma_gwb), np.matmul(np.diag(self.S_effIJ[i][fidxs]), self.dmu_dbeta_pl_gwb(Agwb, + gamma_gwb).T)) + fishers = fisher + fishers - return mathcalF + return fishers def GWBFisherUncertainty(self, - Agwb=None, alpha_gwb=None, + Agwb=None, gamma_gwb=None, log10_gw_rhos=None, psd_type='powerlaw', - components=14,psr=None): + components=20,psr=None): """ GWB Fisher Matrix Uncertainty uncertainty of the i'th param is the square root of the the (i,i) element of the inverse Fisher matrix """ - mathcalF = self.GWBFisherMatrix(Agwb, alpha_gwb, + mathcalF = self.GWBFisherMatrix(Agwb, gamma_gwb, log10_gw_rhos, psd_type=psd_type, components=components,psr=psr) - print(mathcalF.shape) + mathcalFinv = np.linalg.inv(mathcalF) + print("FishMatInv", mathcalF.shape) #mathcalFinv = mathcalF return np.sqrt(np.diag(mathcalFinv)) From ceab77cd90b440c7946ff2881247e9155989d214 Mon Sep 17 00:00:00 2001 From: Jeremy Baier Date: Tue, 24 Dec 2024 14:52:30 -0800 Subject: [PATCH 32/32] docstring i think --- hasasia/sensitivity.py | 16 +++++++--------- 1 file changed, 7 insertions(+), 9 deletions(-) diff --git a/hasasia/sensitivity.py b/hasasia/sensitivity.py index 5eedcfb..230e093 100644 --- a/hasasia/sensitivity.py +++ b/hasasia/sensitivity.py @@ -879,7 +879,7 @@ def GWBFisherMatrix(self, Agwb, gamma_gwb, log10_gw_rhos, psd_type='powerlaw', - components=20,psr=None): + components=20): """ GWB Fisher Matrix Fisher matrix = signal derivative x noise covariance x signal derivative @@ -911,18 +911,16 @@ def GWBFisherMatrix(self, # print("FishMat", mathcalF.shape) fishers = np.zeros((2,2)) for i in range(self.S_effIJ.shape[0]): - fisher = np.matmul(self.dmu_dbeta_pl_gwb(Agwb, - gamma_gwb), np.matmul(np.diag(self.S_effIJ[i][fidxs]), self.dmu_dbeta_pl_gwb(Agwb, - gamma_gwb).T)) + fisher = np.matmul(dmu_dbeta, np.matmul(np.diag(self.S_effIJ[i][fidxs]), dmu_dbeta.T)) fishers = fisher + fishers return fishers def GWBFisherUncertainty(self, - Agwb=None, gamma_gwb=None, - log10_gw_rhos=None, - psd_type='powerlaw', - components=20,psr=None): + Agwb=None, gamma_gwb=None, + log10_gw_rhos=None, + psd_type='powerlaw', + components=20): """ GWB Fisher Matrix Uncertainty uncertainty of the i'th param is the square root of the the @@ -932,7 +930,7 @@ def GWBFisherUncertainty(self, mathcalF = self.GWBFisherMatrix(Agwb, gamma_gwb, log10_gw_rhos, psd_type=psd_type, - components=components,psr=psr) + components=components) mathcalFinv = np.linalg.inv(mathcalF) print("FishMatInv", mathcalF.shape)