From 7df55781ea80534f404be329c5c453346152baac Mon Sep 17 00:00:00 2001 From: Jeffrey Hazboun Date: Tue, 24 Oct 2023 13:19:25 -0700 Subject: [PATCH 1/6] allow alpha_gwb to be set by user --- hasasia/sim.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/hasasia/sim.py b/hasasia/sim.py index 2ad5d54..b1a2479 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, alpha_gwb=-2/3., fast=True, kwastro={'RADEC':True, 'PROPER':True, 'PX':True}): """ Make a simulated pulsar timing array. Using the available parameters, @@ -121,7 +121,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=-2/3., + alpha=alpha_gwb, freqs=freqs) N += corr_from_psd(freqs=freqs, psd=gwb, toas=toas, fast=fast) From d6deb3b7fc5e6a2a8fdba7431483094ffd554eeb Mon Sep 17 00:00:00 2001 From: GourlieK Date: Mon, 11 Mar 2024 18:11:16 -0700 Subject: [PATCH 2/6] Final changes for Jeremy --- hasasia/sensitivity.py | 30 ++++++++++++++++++------------ 1 file changed, 18 insertions(+), 12 deletions(-) diff --git a/hasasia/sensitivity.py b/hasasia/sensitivity.py index b60378f..16fb32e 100644 --- a/hasasia/sensitivity.py +++ b/hasasia/sensitivity.py @@ -7,6 +7,8 @@ import scipy.linalg as sl import os, pickle from astropy import units as u +import jax.numpy as jnp +import jax.scipy as jsc import hasasia from .utils import create_design_matrix @@ -62,14 +64,14 @@ def R_matrix(designmatrix, N): n,m = M.shape L = np.linalg.cholesky(N) Linv = np.linalg.inv(L) - U,s,_ = np.linalg.svd(np.matmul(Linv,M), full_matrices=True) + U,s,_ = np.linalg.svd(jnp.matmul(Linv,M), full_matrices=True) Id = np.eye(M.shape[0]) S = np.zeros_like(M) S[:m,:m] = np.diag(s) - inner = np.linalg.inv(np.matmul(S.T,S)) - outer = np.matmul(S,np.matmul(inner,S.T)) + inner = np.linalg.inv(jnp.matmul(S.T,S)) + outer = jnp.matmul(S,jnp.matmul(inner,S.T)) - return Id - np.matmul(L,np.matmul(np.matmul(U,outer),np.matmul(U.T,Linv))) + return Id - jnp.matmul(L,jnp.matmul(jnp.matmul(U,outer),jnp.matmul(U.T,Linv))) def G_matrix(designmatrix): """ @@ -169,7 +171,7 @@ def get_Tf(designmatrix, toas, N=None, nf=200, fmin=None, fmax=2e-7, m = G.shape[1] Gtilde = np.zeros((ff.size,G.shape[1]),dtype='complex128') Gtilde = np.dot(np.exp(1j*2*np.pi*ff[:,np.newaxis]*toas),G) - Tmat = np.matmul(np.conjugate(Gtilde),Gtilde.T)/N_TOA + Tmat = jnp.matmul(np.conjugate(Gtilde),Gtilde.T)/N_TOA if twofreqs: Tmat = np.real(Tmat) else: @@ -261,10 +263,14 @@ def get_NcalInv(psr, nf=200, fmin=None, fmax=2e-7, freqs=None, Gtilde = np.dot(np.exp(1j*2*np.pi*ff[:,np.newaxis]*toas),G) # N_freq x N_TOA-N_par - Ncal = np.matmul(G.T,np.matmul(psr.N,G)) #N_TOA-N_par x N_TOA-N_par - 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 + L = jsc.linalg.cholesky(psr.N) + A = jnp.matmul(L,G) + del L + Ncal = jnp.matmul(A.T,A) + del A + NcalInv = jnp.linalg.inv(Ncal) + + TfN = jnp.matmul(np.conjugate(Gtilde),jnp.matmul(NcalInv,Gtilde.T)) / 2 if return_Gtilde_Ncal: return np.real(TfN), Gtilde, Ncal elif full_matrix: @@ -810,7 +816,7 @@ def get_NcalInvIJ(psrs, A_GWB, freqs, full_matrix=False, # C_h = sl.block_diag(*[corr_from_psd(freqs=freqs, psd=psd, # toas=p.toas, fast=True) for p in psrs]) C = C_n + C_h - Ncal = np.matmul(G.T, np.matmul(C, G)) #N_TOA-N_par x N_TOA-N_par + Ncal = jnp.matmul(G.T, jnp.matmul(C, G)) #N_TOA-N_par x N_TOA-N_par NcalInv = np.linalg.inv(Ncal) #N_TOA-N_par x N_TOA-N_par TfN = NcalInv#np.matmul(G, np.matmul(NcalInv, G.T)) @@ -1066,7 +1072,7 @@ def corr_from_psd(freqs, psd, toas, fast=True): df = np.diff(freqs) df = np.append(df,df[-1]) tm = np.sqrt(psd*df)*np.exp(1j*2*np.pi*freqs*toas[:,np.newaxis]) - integrand = np.matmul(tm, np.conjugate(tm.T)) + integrand = jnp.matmul(tm, np.conjugate(tm.T)) return np.real(integrand) else: #Makes much larger arrays, but uses np.trapz t1, t2 = np.meshgrid(toas, toas, indexing='ij') @@ -1107,7 +1113,7 @@ def corr_from_psdIJ(freqs, psd, toasI, toasJ, fast=True): df = np.append(df,df[-1]) tmI = np.sqrt(psd*df)*np.exp(1j*2*np.pi*freqs*toasI[:,np.newaxis]) tmJ = np.sqrt(psd*df)*np.exp(1j*2*np.pi*freqs*toasJ[:,np.newaxis]) - integrand = np.matmul(tmI, np.conjugate(tmJ.T)) + integrand = jnp.matmul(tmI, np.conjugate(tmJ.T)) return np.real(integrand) else: #Makes much larger arrays, but uses np.trapz t1, t2 = np.meshgrid(toasI, toasJ, indexing='ij') From 1e2134e378bf290e8086b2867f9bb32762d48480 Mon Sep 17 00:00:00 2001 From: Jeff Hazboun Date: Tue, 3 Dec 2024 12:19:48 -0800 Subject: [PATCH 3/6] Angle-Averaged Detection Probability (#1) * allow alpha_gwb to be set by user * tweaks to get correct factors of 2 * remove functions to skymap * add fidx method * adding advanced detection prob funcs --- hasasia/sensitivity.py | 5 +- hasasia/skymap.py | 178 ++++++++++++++++++++++++++++++++++++----- hasasia/utils.py | 63 +-------------- 3 files changed, 161 insertions(+), 85 deletions(-) diff --git a/hasasia/sensitivity.py b/hasasia/sensitivity.py index e3067ac..514dbe8 100644 --- a/hasasia/sensitivity.py +++ b/hasasia/sensitivity.py @@ -96,7 +96,7 @@ def get_Tf(designmatrix, toas, N=None, nf=200, fmin=None, fmax=2e-7, freqs=None, exact_astro_freqs = False, from_G=True, twofreqs=False, Gmatrix=None): """ - Calculate the transmission function for a given pulsar design matrix, TOAs + the transmission function for a given pulsar design matrix, TOAs and TOA errors. Parameters @@ -751,6 +751,9 @@ def to_pickle(self, filepath): with open(filepath, "wb") as fout: pickle.dump(self, fout) + def fidx(self,f): + return np.argmin(np.abs(self.freqs-f)) + @property def S_eff(self): """Strain power sensitivity. """ diff --git a/hasasia/skymap.py b/hasasia/skymap.py index d8ab87f..f50d429 100644 --- a/hasasia/skymap.py +++ b/hasasia/skymap.py @@ -2,7 +2,10 @@ from __future__ import print_function """Main module.""" import numpy as np -import scipy.special as spec +import scipy.optimize as sopt +import scipy.special as ssp +import scipy.integrate as si +import scipy.stats as ss import healpy as hp import astropy.units as u import astropy.constants as c @@ -105,7 +108,7 @@ def __init__(self, spectra, theta_gw, phi_gw, self.sky_response = (0.5 * self.sky_response[np.newaxis,:,:] * self.pt_sqr) - def SNR(self, h0, iota=None, psi=None): + def SNR(self, h0, iota=None, psi=None, fidx=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 @@ -116,12 +119,17 @@ def SNR(self, h0, iota=None, psi=None): .. _[1]: https://arxiv.org/abs/1907.04341 ''' + 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.') + raise NotImplementedError('Currently cannot marginalize over polarization angle only.') else: S_eff = self.S_eff_full(iota, psi) + + if fidx is not None: + S_eff = S_eff[fidx,:] + return h0 * np.sqrt(self.Tspan / S_eff) @@ -191,8 +199,43 @@ def sky_response_full(self, iota, psi=None): ''' 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 + Ap_sqr = (1 + np.cos(iota)**2)**2 # (0.5 * (1 + np.cos(iota)**2))**2 + Ac_sqr = (2*np.cos(iota))**2 + 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 + c_pluplus, c_pluscross, c_crosscross = self._angle_coefficients(iota=iota, psi=psi) + if isinstance(Ap_sqr,np.ndarray) or isinstance(psi,np.ndarray): + term1 = self.Fplus[:,:,np.newaxis,np.newaxis]**2 * c_pluplus + term2 = 2 * self.Fplus[:,:,np.newaxis,np.newaxis] * self.Fcross[:,:,np.newaxis,np.newaxis] * c_pluscross + term3 = self.Fcross[:,:,np.newaxis,np.newaxis]**2 * c_crosscross + elif isinstance(Ap_sqr,(int,float)) or isinstance(psi,(int,float)): + term1 = self.Fplus**2 * c_pluplus + term2 = 2 * self.Fplus * self.Fcross * c_pluscross + term3 = self.Fcross**2 * c_crosscross + + return term1 + term2 + term3 + + def _angle_coefficients(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]`_, 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})}} + + .. _[1]: https://arxiv.org/abs/1907.04341 + ''' + if iota is None and psi is not None: + raise NotImplementedError('Currently can marginalize over inclination but not phase.') + Ap_sqr = (1 + np.cos(iota)**2)**2 # (0.5 * (1 + np.cos(iota)**2))**2 + Ac_sqr = (2*np.cos(iota))**2 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 @@ -206,21 +249,15 @@ def sky_response_full(self, iota, psi=None): 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 + c_pluplus = Ac_sqr[:,np.newaxis]*spsi**2 + Ap_sqr[:,np.newaxis]*cpsi**2 + c_pluscross = (Ap_sqr[:,np.newaxis] - Ac_sqr[:,np.newaxis]) * cpsi * spsi + c_crosscross = Ap_sqr[:,np.newaxis]*spsi**2 + Ac_sqr[:,np.newaxis]*cpsi**2 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 + c_pluplus = Ac_sqr*spsi**2 + Ap_sqr*cpsi**2 + c_pluscross = (Ap_sqr - Ac_sqr) * cpsi * spsi + c_crosscross = Ap_sqr*spsi**2 + Ac_sqr*cpsi**2 - return term1 + term2 + term3 + return c_pluplus, c_pluscross, c_crosscross def S_SkyI_full(self, iota, psi): """Per Pulsar Strain power sensitivity. """ @@ -262,17 +299,114 @@ def S_eff_full(self, iota, psi): self._S_eff_full = 1.0 / np.sum(self.S_SkyI_full(iota, psi), axis=1) return self._S_eff_full + + def fap(self, F, Npsrs=None): + ''' + False alarm probability of the F-statistic + Use None for the Fe statistic and the number of pulsars for the Fp stat. + ''' + if Npsrs is None: + N = [0,1] + elif isinstance(Npsrs,int): + N = np.arange((4*Npsrs)/2-1, dtype=float) + # else: + # raise ValueError('Npsrs must be an integer or None (for Fe)') + return np.exp(-F)*np.sum([(F**k)/np.math.factorial(k) for k in N]) + + def pdf_F_signal(self, F, snr, Npsrs=None): + if Npsrs is None: + N = 4 + elif isinstance(Npsrs,int): + N = int(4 * Npsrs) + return ss.ncx2.pdf(2*F, N, snr**2) + + def false_dismissal_prob(self, F_thresh, snr=None, Npsrs=None, ave=None, prob_kwargs={'h0':None,'fidx':None}): + ''' + False dismissal probability of the F-statistic + Use None for the Fe statistic and the number of pulsars for the Fp stat. + + Parameters + ---------- + F_thresh : float,array + F-statistic threshold + + snr : float, optional + Signal to noise ratio for a given single source. + + Npsrs : string, optional + Number of pulsar sin the array for using the Fp statistic. + Use None for the Fe statistic. + + Returns + ------- + fdp : float, array + False dismissal probability for a given F-statistic threshold. + + ''' + if Npsrs is None: + N = 4 + elif isinstance(Npsrs,int): + N = int(4 * Npsrs) + + if ave is None or ave=='none': + return ss.ncx2.cdf(2*F_thresh, df=N, nc=snr**2) + elif ave=='snr': + return ss.chi2.cdf(2*F_thresh, df=N, loc=snr**2) + elif ave=='prob': + if prob_kwargs['fidx'] is None: + raise ValueError('fidx must be set!') + return self._fdp_angle_averaged(F_thresh, prob_kwargs['h0'], prob_kwargs['fidx']) + + def detection_prob(self, F_thresh, snr=None, Npsrs=None, ave=None, prob_kwargs={'h0':None,'fidx':None}): + ''' + Detection probability of the F-statistic + Use None for the Fe and the number of pulsars for the Fp stat. + ''' + return 1 - self.false_dismissal_prob(F_thresh, snr, Npsrs, ave, prob_kwargs) + + def _fdp_angle_averaged(self, F_thresh, h0, fidx): + ''' + The angle-averaged false dismissal probablity. See arXiv.... + ''' + integrand = lambda psi, iota: np.sin(iota)/np.pi*ss.ncx2.cdf(2*F_thresh, df=4, + nc=self.SNR(h0,iota,psi,fidx).mean()**2) + return si.dblquad(integrand,0,np.pi,-np.pi/4,np.pi/4)[0] + + def _solve_F_given_fap(self, fap0=0.003, Npsrs=None): + return sopt.fsolve(lambda F :self.fap(F, Npsrs=Npsrs)-fap0, 10) + + def _solve_F_given_fdp_snr(self, 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 :self.false_dismissal_prob(F, snr, Npsrs=Npsrs, iota_psi_ave=iota_psi_ave)-fdp0, F0) + + def _solve_snr_given_fdp_F(self, 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 :self.false_dismissal_prob(F, snr, Npsrs=Npsrs, iota_psi_ave=iota_psi_ave)-fdp0, snr0) + + def _solve_F0_given_SNR(self, 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.) @property def S_eff(self): - """Strain power sensitivity. """ + """ + Strain power sensitivity. NOTE: The prefactors in these expressions are a factor of 4x larger than in + Hazboun, et al., 2019 `[1]` due to a redefinition of h0 to match the one in normal use in the PTA community. + .. _[1]: https://arxiv.org/abs/1907.04341 + """ if not hasattr(self, '_S_eff'): if self.pulsar_term == 'explicit': - self._S_eff = 1.0 / (4./5 * np.sum(self.S_SkyI, axis=1)) + self._S_eff = 1.0 / (16./5 * np.sum(self.S_SkyI, axis=1)) elif self.pulsar_term: - self._S_eff = 1.0 / (12./5 * np.sum(self.S_SkyI, axis=1)) + self._S_eff = 1.0 / (48./5 * np.sum(self.S_SkyI, axis=1)) else: - self._S_eff = 1.0 / (6./5 * np.sum(self.S_SkyI, axis=1)) + self._S_eff = 1.0 / (24./5 * np.sum(self.S_SkyI, axis=1)) return self._S_eff @property diff --git a/hasasia/utils.py b/hasasia/utils.py index 80c75d9..b9b7614 100644 --- a/hasasia/utils.py +++ b/hasasia/utils.py @@ -8,6 +8,7 @@ from astropy.coordinates import SkyCoord import astropy.units as u import astropy.constants as c +# from .skymap import SkySensitivity __all__ = ['create_design_matrix', 'fap', @@ -70,68 +71,6 @@ def create_design_matrix(toas, RADEC=False, PROPER=False, PX=False): return designmatrix -def fap(F, Npsrs=None): - ''' - False alarm probability of the F-statistic - Use None for the Fe statistic and the number of pulsars for the Fp stat. - ''' - if Npsrs is None: - N = [0,1] - elif isinstance(Npsrs,int): - N = np.arange((4*Npsrs)/2-1, dtype=float) - # else: - # raise ValueError('Npsrs must be an integer or None (for Fe)') - return np.exp(-F)*np.sum([(F**k)/np.math.factorial(k) for k in N]) - -def pdf_F_signal(F, snr, Npsrs=None): - if Npsrs is None: - N = 4 - elif isinstance(Npsrs,int): - 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): - ''' - 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: - N = 4 - elif isinstance(Npsrs,int): - N = int(4 * Npsrs) - 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 - 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) - -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) - -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) - -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.) - 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. From 932cfd37dc8201d9dbac598d456107f4436b67bb Mon Sep 17 00:00:00 2001 From: Jeremy Baier Date: Fri, 6 Dec 2024 12:51:40 -0800 Subject: [PATCH 4/6] add some utils and TDP function --- hasasia/sensitivity.py | 6 +++++- hasasia/skymap.py | 18 ++++++++++++++++ hasasia/utils.py | 49 +++++++++++++++++++++++++++++++++++++++++- 3 files changed, 71 insertions(+), 2 deletions(-) diff --git a/hasasia/sensitivity.py b/hasasia/sensitivity.py index 928e1e6..9e3efd5 100644 --- a/hasasia/sensitivity.py +++ b/hasasia/sensitivity.py @@ -758,7 +758,11 @@ def to_pickle(self, filepath): pickle.dump(self, fout) def fidx(self,f): - return np.argmin(np.abs(self.freqs-f)) + """Get the indices of a frequencies in the frequency array.""" + if isinstance(f, int) or isinstance(f, float): + f = np.array([f]) + f = np.asarray(f) + return np.array([np.argmin(abs(ff-self.freqs)) for ff in f]) @property def S_eff(self): diff --git a/hasasia/skymap.py b/hasasia/skymap.py index f50d429..e87370d 100644 --- a/hasasia/skymap.py +++ b/hasasia/skymap.py @@ -363,6 +363,24 @@ def detection_prob(self, F_thresh, snr=None, Npsrs=None, ave=None, prob_kwargs={ Use None for the Fe and the number of pulsars for the Fp stat. ''' return 1 - self.false_dismissal_prob(F_thresh, snr, Npsrs, ave, prob_kwargs) + + def total_detection_probability(self, F_thresh, snr=None, Npsrs=None, ave=None, prob_kwargs={'h0':None, 'fidx': None}): + ''' + See equation ## in Rosado. + Can be interpretted as the probability of detecting a single source + in *any* frequency bin. + + ''' + return 1. - np.prod( + [self.false_dismissal_prob( + F_thresh, + snr=snr, + Npsrs=Npsrs, + ave=ave, + prob_kwargs={'h0':h0,'fidx':fidx}) + for h0, fidx in zip(prob_kwargs['h0'],prob_kwargs['fidx']) + ], + axis=0) def _fdp_angle_averaged(self, F_thresh, h0, fidx): ''' diff --git a/hasasia/utils.py b/hasasia/utils.py index b9b7614..12c6b38 100644 --- a/hasasia/utils.py +++ b/hasasia/utils.py @@ -93,6 +93,28 @@ def strain_and_chirp_mass_to_luminosity_distance(h, M_c, f0): * 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 char_strain_to_strain_amp(hc, fc, df): + r''' + Calculate the strain amplitude of single sources given + their characteristic strains. + + Parameters + ---------- + hc : array_like + Characteristic strain of the single sources. + fc : array_like + Observed orbital frequency bin centers. + df : array_like + Observed orbital frequency bin widths. + + Returns + ------- + hs : + Strain amplitude of the single sources. + + ''' + return hc * np.sqrt(df/fc) + def theta_phi_to_SkyCoord(theta, phi): r''' Takes in a celestial longitude and lattitude and returns an `astropy.SkyCoord` object. @@ -140,4 +162,29 @@ def skycoord_to_Jname(skycoord): 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 + return 'J' + coord_pieces[0] + coord_pieces[1] + sign + coord_pieces[2] +coord_pieces[3] + +def distance_on_sphere(lat1, lon1, lat2, lon2): + '''' + Returns the distance between two points on the surface of the unit sphere. + + Parameters + ========== + lat1, lat2 - float + lattitude of first and second point respectively + lon1, lon2 - float + longitude of first and second point respectively + Returns + ======= + distance - float + distance between two points + ''' + dlat = lat2 - lat1 + dlon = lon2 - lon1 + + # Haversine formula for distance + a = np.sin(dlat / 2)**2 + np.cos(lat1) * np.cos(lat2) * np.sin(dlon / 2)**2 + c = 2 * np.arctan2(np.sqrt(a), np.sqrt(1 - a)) + distance = c + + return distance \ No newline at end of file From 1b428b81199935010f9ca98056c2870ec96d2368 Mon Sep 17 00:00:00 2001 From: Jeremy Baier Date: Fri, 6 Dec 2024 21:41:32 -0800 Subject: [PATCH 5/6] dropping a commit before i probably mess things --- hasasia/skymap.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/hasasia/skymap.py b/hasasia/skymap.py index e87370d..311faa2 100644 --- a/hasasia/skymap.py +++ b/hasasia/skymap.py @@ -386,6 +386,7 @@ def _fdp_angle_averaged(self, F_thresh, h0, fidx): ''' The angle-averaged false dismissal probablity. See arXiv.... ''' + snr = self.SNR(h0,0,0,fidx) integrand = lambda psi, iota: np.sin(iota)/np.pi*ss.ncx2.cdf(2*F_thresh, df=4, nc=self.SNR(h0,iota,psi,fidx).mean()**2) return si.dblquad(integrand,0,np.pi,-np.pi/4,np.pi/4)[0] @@ -395,7 +396,7 @@ def _solve_F_given_fap(self, fap0=0.003, Npsrs=None): def _solve_F_given_fdp_snr(self, 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 + F0 = (4*Npsrs+snr**2)/2 return sopt.fsolve(lambda F :self.false_dismissal_prob(F, snr, Npsrs=Npsrs, iota_psi_ave=iota_psi_ave)-fdp0, F0) def _solve_snr_given_fdp_F(self, fdp0=0.05, F=3, Npsrs=None, iota_psi_ave=False): From 6fafcccd9110c27ac43602997dbc8e0b42de616a Mon Sep 17 00:00:00 2001 From: Jeremy Baier Date: Sun, 8 Dec 2024 18:31:11 -0800 Subject: [PATCH 6/6] adding a trapz integration method --- hasasia/skymap.py | 56 +++++++++++++++++++++++++++++++++++++++++------ 1 file changed, 49 insertions(+), 7 deletions(-) diff --git a/hasasia/skymap.py b/hasasia/skymap.py index 311faa2..3864b4e 100644 --- a/hasasia/skymap.py +++ b/hasasia/skymap.py @@ -320,7 +320,7 @@ def pdf_F_signal(self, F, snr, Npsrs=None): N = int(4 * Npsrs) return ss.ncx2.pdf(2*F, N, snr**2) - def false_dismissal_prob(self, F_thresh, snr=None, Npsrs=None, ave=None, prob_kwargs={'h0':None,'fidx':None}): + def false_dismissal_prob(self, F_thresh, snr=None, Npsrs=None, ave=None, prob_kwargs={'int_method': 'dblquad','h0':None,'fidx':None}): ''' False dismissal probability of the F-statistic Use None for the Fe statistic and the number of pulsars for the Fp stat. @@ -355,7 +355,12 @@ def false_dismissal_prob(self, F_thresh, snr=None, Npsrs=None, ave=None, prob_kw elif ave=='prob': if prob_kwargs['fidx'] is None: raise ValueError('fidx must be set!') - return self._fdp_angle_averaged(F_thresh, prob_kwargs['h0'], prob_kwargs['fidx']) + if prob_kwargs['int_method'] == 'dblquad': + return self._fdp_angle_averaged_dblquad(F_thresh, prob_kwargs['h0'], prob_kwargs['fidx']) + elif prob_kwargs['int_method'] == 'trapz': + return self._fdp_angle_averaged_trapz(prob_kwargs['snr_grid'], F_thresh, prob_kwargs['h0'], prob_kwargs['fidx']) + else: + raise ValueError('int_method must be dblquad or trapz') def detection_prob(self, F_thresh, snr=None, Npsrs=None, ave=None, prob_kwargs={'h0':None,'fidx':None}): ''' @@ -371,25 +376,62 @@ def total_detection_probability(self, F_thresh, snr=None, Npsrs=None, ave=None, in *any* frequency bin. ''' + h0s = prob_kwargs['h0'] + fidxs = np.arange(len(prob_kwargs['fidx'])) return 1. - np.prod( [self.false_dismissal_prob( F_thresh, snr=snr, Npsrs=Npsrs, ave=ave, - prob_kwargs={'h0':h0,'fidx':fidx}) - for h0, fidx in zip(prob_kwargs['h0'],prob_kwargs['fidx']) + prob_kwargs={'h0':h0, + 'fidx':fidx, + 'int_method':prob_kwargs['int_method'], + 'snr_grid':prob_kwargs['snr_grid']}) + for h0, fidx in zip(h0s,fidxs) ], - axis=0) + axis=0) - def _fdp_angle_averaged(self, F_thresh, h0, fidx): + def _fdp_angle_averaged_dblquad(self, F_thresh, h0, fidx): ''' The angle-averaged false dismissal probablity. See arXiv.... ''' - snr = self.SNR(h0,0,0,fidx) integrand = lambda psi, iota: np.sin(iota)/np.pi*ss.ncx2.cdf(2*F_thresh, df=4, nc=self.SNR(h0,iota,psi,fidx).mean()**2) return si.dblquad(integrand,0,np.pi,-np.pi/4,np.pi/4)[0] + + def _fdp_angle_averaged_trapz(self, snrs, F_thresh, h0, fidx): + ''' + The angle-averaged false dismissal probablity. See arXiv.... + ''' + # define the points for integration + iotas = np.linspace(0, np.pi, snrs.shape[0]) + psis = np.linspace(-np.pi/4, np.pi/4, snrs.shape[1]) + # 2d array of integrand values + Z = np.sin(iotas)/np.pi*ss.ncx2.cdf(2*F_thresh, + df=4, + nc=(h0*snrs[:, :, fidx])**2) + # perform the 2D integration using the trapezoidal method twice + return np.trapz(np.trapz(Z, iotas, axis=1), psis, axis=0) + + def sky_ave_SNR_gridded(self, iota, psi, fidxs=None): + r''' + Calculate signal-to-noise ratio across a grid of iota and psi values. + **Note**: This isn't actually SNR but SNR divided by the strain amplitude. + This allows the values to be used for any signal value. + ''' + if fidxs is None: # trick so that all the frequencies get used + fidxs = np.arange(len(self.freqs)) + grid = [] + for i, iot in enumerate(iota): + column = [] + for j, ps in enumerate(psi): + column.append( + np.sqrt(self.Tspan/self.S_eff_full(iot, ps)[fidxs,:]).mean(axis=1) + ) + grid.append(column) + + return np.array(grid) def _solve_F_given_fap(self, fap0=0.003, Npsrs=None): return sopt.fsolve(lambda F :self.fap(F, Npsrs=Npsrs)-fap0, 10)